subroutine da_transfer_kmatoxb(xbx, grid) 2,15
!---------------------------------------------------------------------------
! Purpose: Transfers fields from KMA to first guess (xb) structure.
!---------------------------------------------------------------------------
implicit none
type (xbx_type), intent(inout) :: xbx ! Header & non-gridded vars.
type(domain), intent(inout) :: grid
integer :: i, j, k
integer :: is, ie, js, je, ks, ke
real :: tmpvar, earth_radius_sq,conv
real :: tmpp,tmp_p,tmpps,tmp_ps
! real :: rgh_fac(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme)
character(len=19) :: current_date
real :: TV(grid%xp%kms:grid%xp%kme)
real :: ALPHA(grid%xp%kms:grid%xp%kme)
real :: PHALF(grid%xp%kms:grid%xp%kme)
real :: PHALFL(grid%xp%kms:grid%xp%kme)
real :: pu, pd, coef, delp, hydro, rgasg, shgt
real, dimension(jds:jde) :: loc_latc_mean
real, allocatable :: arrayglobal(:,:)
if (trace_use) call da_trace_entry
("da_transfer_kmatoxb")
if (use_ssmiretrievalobs) then
call da_error(__FILE__,__LINE__, &
(/"Cannot use ssmi obs with kma global runs"/))
end if
conv = pi/180.0
earth_radius_sq = earth_radius*1000.0 * earth_radius*1000.0 * &
conv*(360.0/(coarse_ix-1))*(180.0/(coarse_jy-2))*conv
COEF=0.6080
RGASG = gas_constant/gravity
!---------------------------------------------------------------------------
! Set xb array range indices for processor subdomain.
!---------------------------------------------------------------------------
is = grid%xp % its
ie = grid%xp % ite
js = grid%xp % jts
je = grid%xp % jte
ks = grid%xp % kts
ke = grid%xp % kte
grid%xb % ds = grid%dx
grid%xb % kma_a = grid%kma_a
grid%xb % kma_b = grid%kma_b
if (print_detail_xb) then
write(unit=stdout, fmt='(/a/)') &
'lvl kma_a kma_b'
do k=ks,ke
write(unit=stdout, fmt='(i3,8e20.12)') k, grid%xb%kma_a(k), grid%xb%kma_b(k)
end do
end if
grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
!---------------------------------------------------------------------------
! KMA-specific fields:
!---------------------------------------------------------------------------
! Fix ptop as 0.4 hPa.
ptop = 40.0
grid%xb % ptop = ptop
!---------------------------------------------------------------------------
! Convert KMA fields to xb:
!---------------------------------------------------------------------------
grid%xb%lat(is:ie,js:je) = grid%xlat(is:ie,js:je)
grid%xb%lon(is:ie,js:je) = grid%xlong(is:ie,js:je)
#ifdef DM_PARALLEL
#include "HALO_PSICHI_UV_ADJ.inc"
#endif
tmpvar = 1.0/real(ide-ids+1)
!---------------------------------------------------------------
! transfer u,v,t & q(specific humidity in g/g)
!---------------------------------------------------------------
do k=ks,ke
do j=js,je
do i=is,ie
grid%xb%u(i,j,k) = grid%u_2(i,j,k)
grid%xb%v(i,j,k) = grid%v_2(i,j,k)
grid%xb%t(i,j,k) = grid%t_2(i,j,k)
grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
end do
end do
end do
!---------------------------------------------------------------------------
! Fix wind at Poles
!---------------------------------------------------------------------------
call da_get_vpoles(grid%xb%u,grid%xb%v, &
ids,ide,jds,jde, &
grid%xp%ims,grid%xp%ime,grid%xp%jms,grid%xp%jme,grid%xp%kms,grid%xp%kme, &
grid%xp%its,grid%xp%ite,grid%xp%jts,grid%xp%jte,grid%xp%kts,grid%xp%kte )
if (print_detail_xb) then
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%its, grid%xp%ite=', grid%xp%its, grid%xp%ite
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jts, grid%xp%jte=', grid%xp%jts, grid%xp%jte
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kts, grid%xp%kte=', grid%xp%kts, grid%xp%kte
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%ims, grid%xp%ime=', grid%xp%ims, grid%xp%ime
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jms, grid%xp%jme=', grid%xp%jms, grid%xp%jme
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kms, grid%xp%kme=', grid%xp%kms, grid%xp%kme
write(unit=stdout,fmt='(A,2I5)') 'ids, ide=', ids, ide
write(unit=stdout,fmt='(A,2I5)') 'jds, jde=', jds, jde
write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kds, grid%xp%kde=', grid%xp%kds, grid%xp%kde
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=1)=', size(grid%xb%u, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=2)=', size(grid%xb%u, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=3)=', size(grid%xb%u, dim=3)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=1)=', size(grid%xb%v, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=2)=', size(grid%xb%v, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=3)=', size(grid%xb%v, dim=3)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=1)=', size(grid%xa%u, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=2)=', size(grid%xa%u, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=3)=', size(grid%xa%u, dim=3)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=1)=', size(grid%xa%v, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=2)=', size(grid%xa%v, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=3)=', size(grid%xa%v, dim=3)
write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=1)=', size(grid%u_2, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=2)=', size(grid%u_2, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=3)=', size(grid%u_2, dim=3)
write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=1)=', size(grid%v_2, dim=1)
write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=2)=', size(grid%v_2, dim=2)
write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=3)=', size(grid%v_2, dim=3)
end if
ALPHA(ke) = LOG(2.0)
do j=js,je
do i=is,ie
grid%xb%cori(i,j) = grid%f(i,j)
grid%xb%psfc(i,j) = grid%psfc(i,j)
grid%xb%terr(i,j) = grid%ht(i,j)
! Zhiquan Liu add some RTTOV variables
!--------------------------------------
grid%xb%t2(i,j) = grid%t2(i,j)
grid%xb%q2(i,j) = grid%q2(i,j)
grid%xb%u10(i,j) = grid%u10(i,j)
grid%xb%v10(i,j) = grid%v10(i,j)
grid%xb%tsk(i,j) = grid%tsk(i,j)
! currently KMA guess have no SST, replaced by TSK
! grid%xb%tgrn(i,j) = grid%sst(i,j)
grid%xb%tgrn(i,j) = grid%tsk(i,j)
grid%xb%landmask(i,j) = grid%landmask(i,j)
grid%xb%xland(i,j) = grid%xland(i,j)
grid%xb%smois(i,j) = grid%smois(i,j,1)
grid%xb%tslb(i,j) = grid%tslb(i,j,1)
grid%xb%xice(i,j) = grid%xice(i,j)
grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
grid%xb%isltyp(i,j) = grid%isltyp(i,j)
grid%xb%vegfra(i,j) = grid%vegfra(i,j)
! convert snow water equivalent kg/m2 to snow depth(mm)
grid%xb%snowh(i,j) = grid%snow(i,j)*10.0
!---------------------------------------------------------------------
! Syed RH Rizvi
!
! The following code is to be activated after
! getting sst, land use, land mask and snow information for KMA grid
! This infor needs to be packed in "KMA2NETCDF" convertor
!
!---------------------------------------------------------------------
! grid%xb%tgrn(i,j) = grid%sst(i,j)
! if (grid%xb%tgrn(i,j) < 100.0) grid%xb%tgrn(i,j) = tmn(i,j)
! grid%xb%lanu(i,j) = grid%lu_index(i,j)
! grid%xb%landmask(i,j) = grid%landmask(i,j)
! grid%xb%xland(i,j) = grid%xland(i,j)
! grid%xb%snow(i,j) = grid%snowc(i,j)
! Since data is on full levels get full level pr & ht.
! Note p = A + B * Psfc formula needs Psfc to be in hPa
do K = ks,ke-1
PU = grid%xb%KMA_A(K+1) + grid%xb%KMA_B(K+1)*grid%xb%psfc(i,j)/100.0
PD = grid%xb%KMA_A(K ) + grid%xb%KMA_B(K )*grid%xb%psfc(i,j)/100.0
if (PU == PD)then
message(1)='PU, PD equal and so denominator will be zero'
write(unit=message(2),fmt='(A,3I5)') ' i, j, k = ',i,j,k
write(unit=message(3),fmt='(A,2F10.2)') &
' KMA_A ',grid%KMA_A(k),grid%KMA_A(K+1)
write(unit=message(4),fmt='(A,2F10.2)') &
' KMA_B ',grid%KMA_B(k),grid%KMA_B(K+1)
write(unit=message(5),fmt='(A,2F10.2)') &
' psfc(Pa) = ',grid%xb%psfc(i,j)
call da_error(__FILE__,__LINE__,message(1:5))
end if
grid%xb%p(i,j,k)= 100.0 * exp ((PD*LOG(PD)-PU*LOG(PU))/(PD-PU) -1.0)
end do
grid%xb%p(i,j,ke) = 100.0*(grid%xb%KMA_A(ke)+grid%xb%KMA_B(ke)*grid%xb%psfc(i,j)/100.0)/2.0
do k=ks,ke
if (grid%xb%t(i,j,k) <= 0.0) then
write(unit=message(1),fmt='(A,3I5,F10.2)') &
'Problem in temp = ',i,j,k,grid%xb%t(i,j,k)
call da_error(__FILE__,__LINE__,message(1:1))
end if
grid%xb%rho(i,j,k)=grid%xb%p(i,j,k)/(gas_constant* &
(1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K))
end do
! compute full level height
do K = ks,ke
PHALF(K) = grid%xb%KMA_A(K) + grid%xb%KMA_B(K)*grid%xb%psfc(I,J)/100.0
TV (K) = (1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K)*rgasg
end do
do K = ks,ke-1
DELP = PHALF(K) - PHALF(K+1)
PHALFL(K)= LOG(PHALF(K)/PHALF(K+1))
ALPHA(K) = 1.0-PHALF(K+1)*PHALFL(K)/DELP
end do
SHGT = grid%xb%terr(i,j)
do K = ks, ke
grid%xb%h(I,J,K) = SHGT + ALPHA(K)*TV(K)
end do
HYDRO = 0.0
do K = ks+1, ke
HYDRO = HYDRO + TV(K-1)*PHALFL(K-1)
grid%xb%h(I,J,K) = grid%xb%h(I,J,K) + HYDRO
end do
end do
end do
!---------------------------------------------------------------------------
! Sigma values are needed for interpolating
! background error statistics in vertical
! Fix sigmah values based on global average surface pressure
! and level wise pressure
!---------------------------------------------------------------------------
tmp_ps = sum(grid%xb%psfc(is:ie,js:je)) /real((ide-ids+1)*(jde-jds+1))
tmpps = wrf_dm_sum_real(tmp_ps)
do k=ks,ke
tmp_p = sum(grid%xb%p(is:ie,js:je,k)) /real((ide-ids+1)*(jde-jds+1))
tmpp = wrf_dm_sum_real(tmp_p)
! 0.01 is added to see that sigmah should not become negative
grid%xb%sigmah(ke+ks-k) = (tmpp - ptop+0.01) / (tmpps- ptop+0.01)
if (grid%xb%sigmah(ke+ks-k) < 0.0) then
write(unit=message(1),fmt='(A,F10.2)') &
' average surface pressure = ',tmpps
write(unit=message(2),fmt='(A,I5,A,F10.2)') &
' average pressure at this level= ',k,' = ',tmpp
write(unit=message(3),fmt='(A,I5,A,F10.2)') &
' negative sigmah(',ke+ks-k,') = ',grid%xb%sigmah(ke+ks-k)
call da_error(__FILE__,__LINE__,message(1:3))
end if
end do
! Fix ztop
grid%xb%ztop = grid%xb%h(is,js,ke)
if (print_detail_xb) then
write(unit=stdout, fmt='(/5a/)') &
'lvl h p t'
do k=ks,ke
write(unit=stdout, fmt='(i3,8e20.12)') k, &
grid%xb%h(is,js,k), grid%xb%p(is,js,k), grid%xb%t(is,js,k)
end do
write (unit=stdout,fmt=*) ' '
write (unit=stdout,fmt=*) 'grid%xb%u(is,je,ke)=', grid%xb%u(is,je,ke)
write (unit=stdout,fmt=*) 'grid%xb%v(ie,js,ke)=', grid%xb%v(ie,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%w(is,js,ke)=', grid%xb%w(is,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%t(is,js,ke)=', grid%xb%t(is,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%p(is,js,ke)=', grid%xb%p(is,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%q(is,js,ke)=', grid%xb%q(is,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%hf(is,js,ke)=', grid%xb%hf(is,js,ke)
write (unit=stdout,fmt=*) 'grid%xb%map_factor(is,js)=', grid%xb%map_factor(is,js)
write (unit=stdout,fmt=*) 'grid%xb%cori(is,js)=', grid%xb%cori(is,js)
write (unit=stdout,fmt=*) 'grid%xb%tgrn(is,js)=', grid%xb%tgrn(is,js)
write (unit=stdout,fmt=*) 'grid%xb%lat(is,js)=', grid%xb%lat(is,js)
write (unit=stdout,fmt=*) 'grid%xb%lon(is,js)=', grid%xb%lon(is,js)
write (unit=stdout,fmt=*) 'grid%xb%terr(is,js)=', grid%xb%terr(is,js)
write (unit=stdout,fmt=*) 'grid%xb%snow(is,js)=', grid%xb%snow(is,js)
write (unit=stdout,fmt=*) 'grid%xb%lanu(is,js)=', grid%xb%lanu(is,js)
write (unit=stdout,fmt=*) 'grid%xb%landmask(is,js)=', grid%xb%landmask(is,js)
write (unit=stdout,fmt=*) ' '
end if
write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
grid%start_year, '-', &
grid%start_month, '-', &
grid%start_day, '_', &
grid%start_hour, ':', &
grid%start_minute,':', &
grid%start_second
write(unit=stdout,fmt='(2A)') 'Current date is ',current_date
!---------------------------------------------------------------------------
! Syed RH Rizvi
!
! Following code for calculating roughness length needs to be activated
! after getting land use data for KMA grid
!---------------------------------------------------------------------------
!
! Set proper value for "mminlu" if using KMA info
!xbx % mminlu = 'USGS'
!call da_roughness_from_lanu(19, xbx % mminlu, current_date, grid%xp, &
! grid%xb % lanu(grid%xp%ims,grid%xp%jms), grid%xb % rough(grid%xp%ims,grid%xp%jms))
do j=js,je
do i=is,ie
if (grid%xb%ztop < grid%xb%hf(i,j,ke+1)) grid%xb%ztop = grid%xb%hf(i,j,ke+1)
! Calculate grid_box area and vertical inner product for use in
! vertical transform:
grid%xb % grid_box_area(i,j) = earth_radius_sq * cos(grid%xlat(i,j)*conv)
if (vertical_ip == vertical_ip_sqrt_delta_p) then
! Vertical inner product is sqrt(Delta p):
do k=1,grid%xb%mkz
grid%xb % vertical_inner_product(i,j,k) = &
sqrt(grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1))
end do
else if (vertical_ip == vertical_ip_delta_p) then
! Vertical inner product is Delta p:
do k=1,grid%xb%mkz
grid%xb % vertical_inner_product(i,j,k) = grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1)
end do
end if
!------------------------------------------------------------------------------
! Syed RH Rizvi
!
! Following code for calculating roughness length needs to be activated
! to calculate surface variables (10-m wind, and 2-m T, Q) ,
! After testing KMA-WRFVAR for SFC_ASSIM_OPTIONS = 2
!
!------------------------------------------------------------------------------
! call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
! grid%xb%p(i,j,ks), grid%xb%t(i,j,ks), grid%xb%q(i,j,ks), &
! grid%xb%u(i,j,ks), grid%xb%v(i,j,ks), &
! grid%xb%p(i,j,ks+1), grid%xb%t(i,j,ks+1), grid%xb%q(i,j,ks+1), &
! grid%xb%h(i,j,ks), grid%xb%rough(i,j),grid%xb%landmask(i,j), &
! grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), grid%xb%q2(i,j), &
! grid%xb%regime(i,j))
! write (unit=stdout,fmt='(7I5)') &
! i,j,grid%xb%psfc(i,j),grid%xb%t2(i,j),grid%xb%q2(i,j), &
! grid%xb%u10(i,j),grid%xb%v10(i,j)
! ------------------------------------------------------------------------------
end do
end do
!---------------------------------------------------------------------------
! Calculate saturation vapour pressure and relative humidity:
!---------------------------------------------------------------------------
do j=js,je
do k=ks,ke
do i=is,ie
call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
grid%xb % q(i,j,k), &
grid%xb %es(i,j,k), grid%xb %qs(i,j,k), grid%xb %rh(i,j,k))
end do
end do
end do
!---------------------------------------------------------------------------
! Calculate dew point temperature:
!---------------------------------------------------------------------------
call da_trh_to_td (grid)
if (print_detail_xb) then
i=is; j=js; k=ks
write (unit=stdout,fmt='(A,3I5)') 'i,j,k=', i,j,k
write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
write (unit=stdout,fmt='(A)') ' '
end if
!---------------------------------------------------------------------------
! Sea level pressure and total precipitable water
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
! Following code for calculating roughness length needs to be activated
! for calculating roughness length if sea level pressure is desired
!---------------------------------------------------------------------------
! call da_wrf_tpq_2_slp (grid)
call da_integrat_dz(grid)
!---------------------------------------------------------------------------
! Following code for calculating roughness length needs to be activated
! for surface wind speed for SFC_ASSIM_OPTIONS = 2
!
! It is not working because KMA terrain height output from wave to grid
! transform is negative at many grid points
!---------------------------------------------------------------------------
! tmpvar = log(10.0/0.0001)
! do j=js,je
! do i=is,ie
! if (grid%xb%hf(i,j,ks) <= 0) then
! write(unit=message(1), fmt=*) ' zero hf at i/j ',i,j,' ks= ',ks
! call da_error(__FILE__,__LINE__,message(1:1))
! end if
! rgh_fac(i,j) = 1.0/log(grid%xb%hf(i,j,ks)/0.0001)
! grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,ks)*grid%xb%u(i,j,ks) &
! + grid%xb%v(i,j,ks)*grid%xb%v(i,j,ks) + 1.0e-6) &
! *tmpvar*rgh_fac(i,j)
! end do
! end do
!---------------------------------------------------------------------------
! Following code for calculating roughness length needs to be activated
! if SSMI brightness temperature are used
!---------------------------------------------------------------------------
! call da_transform_xtotb(xa)
!---------------------------------------------------------------------------
! Calculate means for later use in setting up background errors.
!---------------------------------------------------------------------------
allocate (xbx % latc_mean(jds:jde))
tmpvar = 1.0/real(ide-ids+1)
loc_latc_mean(:) = 0.0
! Bitwise-exact reduction preserves operation order of serial code for
! testing, at the cost of much-increased run-time. Turn it off when not
! testing. This will always be .false. for a serial or 1-process MPI run.
if (test_dm_exact) then
allocate(arrayglobal(ids:ide, jds:jde))
call da_patch_to_global(grid, grid%xb%lat, arrayglobal)
do j=jds,jde
loc_latc_mean(j) = tmpvar*sum(arrayglobal(ids:ide, j))
end do
deallocate(arrayglobal)
! Broadcast result from monitor to other tasks.
call wrf_dm_bcast_real(loc_latc_mean, (jde-jds+1))
xbx % latc_mean = loc_latc_mean
else
do j=js,je
loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(is:ie, j))
end do
call wrf_dm_sum_reals (loc_latc_mean, xbx % latc_mean)
end if
! WHY?
! do j=js,je
! do i=is,ie
! write(unit=stdout,fmt='(2i4,3f12.2,e12.4,2f12.2)') &
! j,i,grid%xb%psfc(i,j),grid%xb%tsk(i,j),grid%xb%t2(i,j), &
! grid%xb%q2(i,j),grid%xb%u10(i,j),grid%xb%v10(i,j)
! end do
! end do
if (trace_use) call da_trace_exit("da_transfer_kmatoxb")
end subroutine da_transfer_kmatoxb