wrfv2_fire/var/da/da_transfer_model/da_transfer_kmatoxb.inc

   1 subroutine da_transfer_kmatoxb(xbx, grid)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Transfers fields from KMA to first guess (xb) structure.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (xbx_type), intent(inout)     :: xbx          ! Header & non-gridded vars.
  10
  11    type(domain), intent(inout)        :: grid
  12
  13    integer :: i, j, k
  14    integer :: is, ie, js, je, ks, ke
  15    real    :: tmpvar, earth_radius_sq,conv
  16    real    :: tmpp,tmp_p,tmpps,tmp_ps
  17    ! real    :: rgh_fac(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme)
  18    character(len=19) :: current_date
  19
  20    real :: TV(grid%xp%kms:grid%xp%kme)
  21    real :: ALPHA(grid%xp%kms:grid%xp%kme)
  22    real :: PHALF(grid%xp%kms:grid%xp%kme)
  23    real :: PHALFL(grid%xp%kms:grid%xp%kme)
  24
  25    real :: pu, pd, coef, delp, hydro, rgasg, shgt
  26
  27    real, dimension(jds:jde) :: loc_latc_mean
  28
  29    real, allocatable :: arrayglobal(:,:)
  30
  31    if (trace_use) call da_trace_entry("da_transfer_kmatoxb")
  32
  33    if (use_ssmiretrievalobs) then
  34       call da_error(__FILE__,__LINE__, &
  35          (/"Cannot use ssmi obs with kma global runs"/))
  36    end if
  37
  38    conv = pi/180.0
  39    earth_radius_sq = earth_radius*1000.0 * earth_radius*1000.0 * &
  40                      conv*(360.0/(coarse_ix-1))*(180.0/(coarse_jy-2))*conv
  41    COEF=0.6080
  42    RGASG = gas_constant/gravity
  43
  44    !---------------------------------------------------------------------------
  45    ! Set xb array range indices for processor subdomain.
  46    !---------------------------------------------------------------------------
  47
  48    is = grid%xp % its
  49    ie = grid%xp % ite
  50    js = grid%xp % jts
  51    je = grid%xp % jte
  52    ks = grid%xp % kts
  53    ke = grid%xp % kte
  54
  55    grid%xb % ds    = grid%dx
  56    grid%xb % kma_a = grid%kma_a
  57    grid%xb % kma_b = grid%kma_b
  58
  59    if (print_detail_xb) then
  60        write(unit=stdout, fmt='(/a/)') &
  61           'lvl         kma_a                 kma_b'
  62
  63        do k=ks,ke
  64           write(unit=stdout, fmt='(i3,8e20.12)') k, grid%xb%kma_a(k), grid%xb%kma_b(k)
  65       end do
  66    end if
  67
  68    grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
  69    grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
  70    grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
  71
  72    !---------------------------------------------------------------------------
  73    ! KMA-specific fields:
  74    !---------------------------------------------------------------------------
  75
  76    ! Fix ptop as 0.4 hPa.
  77    ptop = 40.0
  78
  79    grid%xb % ptop = ptop
  80
  81    !---------------------------------------------------------------------------
  82    ! Convert KMA fields to xb:
  83    !---------------------------------------------------------------------------
  84
  85    grid%xb%lat(is:ie,js:je) =  grid%xlat(is:ie,js:je)
  86    grid%xb%lon(is:ie,js:je) = grid%xlong(is:ie,js:je)
  87
  88 #ifdef DM_PARALLEL
  89 #include "HALO_PSICHI_UV_ADJ.inc"
  90 #endif
  91
  92    tmpvar  = 1.0/real(ide-ids+1)
  93
  94    !---------------------------------------------------------------
  95    ! transfer u,v,t & q(specific humidity in g/g)
  96    !---------------------------------------------------------------
  97
  98    do k=ks,ke
  99       do j=js,je
 100          do i=is,ie
 101             grid%xb%u(i,j,k) = grid%u_2(i,j,k)
 102             grid%xb%v(i,j,k) = grid%v_2(i,j,k)
 103             grid%xb%t(i,j,k) = grid%t_2(i,j,k)
 104             grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
 105          end do
 106       end do
 107    end do
 108
 109    !---------------------------------------------------------------------------
 110    ! Fix wind at Poles
 111    !---------------------------------------------------------------------------
 112
 113     call da_get_vpoles(grid%xb%u,grid%xb%v,          &
 114        ids,ide,jds,jde, &
 115        grid%xp%ims,grid%xp%ime,grid%xp%jms,grid%xp%jme,grid%xp%kms,grid%xp%kme, &
 116        grid%xp%its,grid%xp%ite,grid%xp%jts,grid%xp%jte,grid%xp%kts,grid%xp%kte )
 117
 118    if (print_detail_xb) then
 119       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%its, grid%xp%ite=', grid%xp%its, grid%xp%ite
 120       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jts, grid%xp%jte=', grid%xp%jts, grid%xp%jte
 121       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kts, grid%xp%kte=', grid%xp%kts, grid%xp%kte
 122       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%ims, grid%xp%ime=', grid%xp%ims, grid%xp%ime
 123       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jms, grid%xp%jme=', grid%xp%jms, grid%xp%jme
 124       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kms, grid%xp%kme=', grid%xp%kms, grid%xp%kme
 125       write(unit=stdout,fmt='(A,2I5)') 'ids, ide=', ids, ide
 126       write(unit=stdout,fmt='(A,2I5)') 'jds, jde=', jds, jde
 127       write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kds, grid%xp%kde=', grid%xp%kds, grid%xp%kde
 128
 129       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=1)=', size(grid%xb%u, dim=1)
 130       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=2)=', size(grid%xb%u, dim=2)
 131       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=3)=', size(grid%xb%u, dim=3)
 132       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=1)=', size(grid%xb%v, dim=1)
 133       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=2)=', size(grid%xb%v, dim=2)
 134       write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=3)=', size(grid%xb%v, dim=3)
 135
 136       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=1)=', size(grid%xa%u, dim=1)
 137       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=2)=', size(grid%xa%u, dim=2)
 138       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=3)=', size(grid%xa%u, dim=3)
 139       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=1)=', size(grid%xa%v, dim=1)
 140       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=2)=', size(grid%xa%v, dim=2)
 141       write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=3)=', size(grid%xa%v, dim=3)
 142
 143       write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=1)=', size(grid%u_2, dim=1)
 144       write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=2)=', size(grid%u_2, dim=2)
 145       write(unit=stdout,fmt='(A,I5)') 'size(grid%u_2, dim=3)=', size(grid%u_2, dim=3)
 146       write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=1)=', size(grid%v_2, dim=1)
 147       write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=2)=', size(grid%v_2, dim=2)
 148       write(unit=stdout,fmt='(A,I5)') 'size(grid%v_2, dim=3)=', size(grid%v_2, dim=3)
 149    end if
 150
 151    ALPHA(ke) = LOG(2.0)
 152
 153    do j=js,je
 154       do i=is,ie
 155          grid%xb%cori(i,j) = grid%f(i,j)
 156          grid%xb%psfc(i,j) = grid%psfc(i,j)
 157          grid%xb%terr(i,j) = grid%ht(i,j)
 158
 159          ! Zhiquan Liu add some RTTOV variables
 160          !--------------------------------------
 161          grid%xb%t2(i,j) = grid%t2(i,j)
 162          grid%xb%q2(i,j) = grid%q2(i,j)
 163          grid%xb%u10(i,j) = grid%u10(i,j)
 164          grid%xb%v10(i,j) = grid%v10(i,j)
 165          grid%xb%tsk(i,j) = grid%tsk(i,j)
 166          ! currently KMA guess have no SST, replaced by TSK
 167          ! grid%xb%tgrn(i,j) = grid%sst(i,j)
 168          grid%xb%tgrn(i,j) = grid%tsk(i,j)
 169          grid%xb%landmask(i,j) = grid%landmask(i,j)
 170          grid%xb%xland(i,j) = grid%xland(i,j)
 171
 172          grid%xb%smois(i,j) = grid%smois(i,j,1)
 173          grid%xb%tslb(i,j) = grid%tslb(i,j,1)
 174          grid%xb%xice(i,j) = grid%xice(i,j)
 175          grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
 176          grid%xb%isltyp(i,j) = grid%isltyp(i,j)
 177          grid%xb%vegfra(i,j) = grid%vegfra(i,j)
 178          ! convert snow water equivalent kg/m2 to snow depth(mm)
 179          grid%xb%snowh(i,j) = grid%snow(i,j)*10.0
 180
 181          !---------------------------------------------------------------------
 182          !  Syed RH Rizvi
 183          !
 184          ! The following code is to be activated after
 185          ! getting sst, land use, land mask and snow information for KMA grid
 186          ! This infor needs to be packed in "KMA2NETCDF" convertor
 187          !
 188          !---------------------------------------------------------------------
 189
 190          ! grid%xb%tgrn(i,j) = grid%sst(i,j)
 191          ! if (grid%xb%tgrn(i,j) < 100.0) grid%xb%tgrn(i,j) = tmn(i,j)
 192          ! grid%xb%lanu(i,j)     = grid%lu_index(i,j)
 193          ! grid%xb%landmask(i,j) = grid%landmask(i,j)
 194          ! grid%xb%xland(i,j) = grid%xland(i,j)
 195          ! grid%xb%snow(i,j)     = grid%snowc(i,j)
 196
 197          ! Since data is on full levels get full level pr & ht.
 198          ! Note p = A + B * Psfc formula needs Psfc to be in hPa
 199
 200          do K = ks,ke-1
 201             PU  = grid%xb%KMA_A(K+1) + grid%xb%KMA_B(K+1)*grid%xb%psfc(i,j)/100.0
 202             PD  = grid%xb%KMA_A(K ) + grid%xb%KMA_B(K )*grid%xb%psfc(i,j)/100.0
 203             if (PU == PD)then
 204                message(1)='PU, PD equal and so denominator will be zero'
 205                write(unit=message(2),fmt='(A,3I5)') ' i, j, k = ',i,j,k
 206                write(unit=message(3),fmt='(A,2F10.2)') &
 207                   ' KMA_A ',grid%KMA_A(k),grid%KMA_A(K+1)
 208                write(unit=message(4),fmt='(A,2F10.2)') &
 209                   ' KMA_B ',grid%KMA_B(k),grid%KMA_B(K+1)
 210                write(unit=message(5),fmt='(A,2F10.2)') &
 211                   ' psfc(Pa) = ',grid%xb%psfc(i,j)
 212                call da_error(__FILE__,__LINE__,message(1:5))
 213             end if
 214             grid%xb%p(i,j,k)= 100.0 * exp ((PD*LOG(PD)-PU*LOG(PU))/(PD-PU) -1.0)
 215          end do
 216
 217          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
 218
 219          do k=ks,ke
 220             if (grid%xb%t(i,j,k) <= 0.0) then
 221                write(unit=message(1),fmt='(A,3I5,F10.2)') &
 222                   'Problem in  temp = ',i,j,k,grid%xb%t(i,j,k)
 223                call da_error(__FILE__,__LINE__,message(1:1))
 224             end if
 225
 226             grid%xb%rho(i,j,k)=grid%xb%p(i,j,k)/(gas_constant*  &
 227                (1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K))
 228          end do
 229
 230          ! compute full level height
 231
 232          do K = ks,ke
 233             PHALF(K) = grid%xb%KMA_A(K) + grid%xb%KMA_B(K)*grid%xb%psfc(I,J)/100.0
 234             TV   (K) = (1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K)*rgasg
 235          end do
 236
 237          do K = ks,ke-1
 238             DELP     = PHALF(K) - PHALF(K+1)
 239             PHALFL(K)= LOG(PHALF(K)/PHALF(K+1))
 240             ALPHA(K) = 1.0-PHALF(K+1)*PHALFL(K)/DELP
 241          end do
 242
 243          SHGT = grid%xb%terr(i,j)
 244          do K = ks, ke
 245             grid%xb%h(I,J,K) = SHGT + ALPHA(K)*TV(K)
 246          end do
 247
 248          HYDRO = 0.0
 249          do K = ks+1, ke
 250             HYDRO = HYDRO + TV(K-1)*PHALFL(K-1)
 251             grid%xb%h(I,J,K) = grid%xb%h(I,J,K) + HYDRO
 252          end do
 253       end do
 254    end do
 255
 256    !---------------------------------------------------------------------------
 257    ! Sigma values are needed for interpolating
 258    ! background error statistics in vertical
 259    ! Fix sigmah values based on global average surface pressure
 260    !                    and level wise pressure
 261    !---------------------------------------------------------------------------
 262
 263    tmp_ps = sum(grid%xb%psfc(is:ie,js:je)) /real((ide-ids+1)*(jde-jds+1))
 264
 265    tmpps = wrf_dm_sum_real(tmp_ps)
 266
 267    do k=ks,ke
 268       tmp_p = sum(grid%xb%p(is:ie,js:je,k)) /real((ide-ids+1)*(jde-jds+1))
 269       tmpp = wrf_dm_sum_real(tmp_p)
 270
 271       ! 0.01 is added to see that sigmah should not become negative
 272       grid%xb%sigmah(ke+ks-k) = (tmpp - ptop+0.01) / (tmpps- ptop+0.01)
 273       if (grid%xb%sigmah(ke+ks-k) < 0.0) then
 274          write(unit=message(1),fmt='(A,F10.2)') &
 275             ' average surface pressure = ',tmpps
 276          write(unit=message(2),fmt='(A,I5,A,F10.2)') &
 277             ' average pressure at this level= ',k,' = ',tmpp
 278          write(unit=message(3),fmt='(A,I5,A,F10.2)') &
 279             ' negative sigmah(',ke+ks-k,') = ',grid%xb%sigmah(ke+ks-k)
 280          call da_error(__FILE__,__LINE__,message(1:3))
 281       end if
 282    end do
 283
 284    ! Fix ztop
 285
 286    grid%xb%ztop = grid%xb%h(is,js,ke)
 287
 288    if (print_detail_xb) then
 289       write(unit=stdout, fmt='(/5a/)') &
 290          'lvl         h                 p                t'
 291
 292       do k=ks,ke
 293          write(unit=stdout, fmt='(i3,8e20.12)') k, &
 294             grid%xb%h(is,js,k), grid%xb%p(is,js,k), grid%xb%t(is,js,k)
 295       end do
 296
 297       write (unit=stdout,fmt=*) ' '
 298       write (unit=stdout,fmt=*) 'grid%xb%u(is,je,ke)=', grid%xb%u(is,je,ke)
 299       write (unit=stdout,fmt=*) 'grid%xb%v(ie,js,ke)=', grid%xb%v(ie,js,ke)
 300       write (unit=stdout,fmt=*) 'grid%xb%w(is,js,ke)=', grid%xb%w(is,js,ke)
 301       write (unit=stdout,fmt=*) 'grid%xb%t(is,js,ke)=', grid%xb%t(is,js,ke)
 302       write (unit=stdout,fmt=*) 'grid%xb%p(is,js,ke)=', grid%xb%p(is,js,ke)
 303       write (unit=stdout,fmt=*) 'grid%xb%q(is,js,ke)=', grid%xb%q(is,js,ke)
 304       write (unit=stdout,fmt=*) 'grid%xb%hf(is,js,ke)=', grid%xb%hf(is,js,ke)
 305       write (unit=stdout,fmt=*) 'grid%xb%map_factor(is,js)=', grid%xb%map_factor(is,js)
 306       write (unit=stdout,fmt=*) 'grid%xb%cori(is,js)=', grid%xb%cori(is,js)
 307       write (unit=stdout,fmt=*) 'grid%xb%tgrn(is,js)=', grid%xb%tgrn(is,js)
 308       write (unit=stdout,fmt=*) 'grid%xb%lat(is,js)=', grid%xb%lat(is,js)
 309       write (unit=stdout,fmt=*) 'grid%xb%lon(is,js)=', grid%xb%lon(is,js)
 310       write (unit=stdout,fmt=*) 'grid%xb%terr(is,js)=', grid%xb%terr(is,js)
 311       write (unit=stdout,fmt=*) 'grid%xb%snow(is,js)=', grid%xb%snow(is,js)
 312       write (unit=stdout,fmt=*) 'grid%xb%lanu(is,js)=', grid%xb%lanu(is,js)
 313       write (unit=stdout,fmt=*) 'grid%xb%landmask(is,js)=', grid%xb%landmask(is,js)
 314       write (unit=stdout,fmt=*) ' '
 315    end if
 316
 317    write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
 318       grid%start_year, '-', &
 319       grid%start_month, '-', &
 320       grid%start_day,   '_', &
 321       grid%start_hour,  ':', &
 322       grid%start_minute,':', &
 323       grid%start_second
 324
 325    write(unit=stdout,fmt='(2A)') 'Current date is ',current_date
 326
 327    !---------------------------------------------------------------------------
 328    ! Syed RH Rizvi
 329    !
 330    ! Following code for calculating roughness length needs to be activated
 331    ! after getting land use data for KMA grid
 332    !---------------------------------------------------------------------------
 333
 334    !
 335    ! Set proper value for "mminlu" if using KMA info
 336
 337    !xbx % mminlu = 'USGS'
 338    !call da_roughness_from_lanu(19, xbx % mminlu, current_date, grid%xp, &
 339    !   grid%xb % lanu(grid%xp%ims,grid%xp%jms), grid%xb % rough(grid%xp%ims,grid%xp%jms))
 340
 341    do j=js,je
 342       do i=is,ie
 343          if (grid%xb%ztop < grid%xb%hf(i,j,ke+1)) grid%xb%ztop = grid%xb%hf(i,j,ke+1)
 344          ! Calculate grid_box area and vertical inner product for use in
 345          ! vertical transform:
 346          grid%xb % grid_box_area(i,j) = earth_radius_sq * cos(grid%xlat(i,j)*conv)
 347
 348          if (vertical_ip == vertical_ip_sqrt_delta_p) then
 349             ! Vertical inner product is sqrt(Delta p):
 350             do k=1,grid%xb%mkz
 351                grid%xb % vertical_inner_product(i,j,k) = &
 352                   sqrt(grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1))
 353             end do
 354          else if (vertical_ip == vertical_ip_delta_p) then
 355             ! Vertical inner product is Delta p:
 356             do k=1,grid%xb%mkz
 357             grid%xb % vertical_inner_product(i,j,k) = grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1)
 358             end do
 359          end if
 360
 361          !------------------------------------------------------------------------------
 362          !  Syed RH Rizvi
 363          !
 364          ! Following code for calculating roughness length needs to be activated
 365          ! to calculate surface variables (10-m wind, and 2-m T, Q) ,
 366          ! After testing KMA-WRFVAR for SFC_ASSIM_OPTIONS = 2
 367          !
 368          !------------------------------------------------------------------------------
 369
 370          ! call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
 371          !    grid%xb%p(i,j,ks), grid%xb%t(i,j,ks), grid%xb%q(i,j,ks), &
 372          !    grid%xb%u(i,j,ks), grid%xb%v(i,j,ks), &
 373          !    grid%xb%p(i,j,ks+1), grid%xb%t(i,j,ks+1), grid%xb%q(i,j,ks+1), &
 374          !    grid%xb%h(i,j,ks), grid%xb%rough(i,j),grid%xb%landmask(i,j), &
 375          !    grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), grid%xb%q2(i,j), &
 376          !    grid%xb%regime(i,j))
 377
 378          ! write (unit=stdout,fmt='(7I5)') &
 379          !    i,j,grid%xb%psfc(i,j),grid%xb%t2(i,j),grid%xb%q2(i,j), &
 380          !    grid%xb%u10(i,j),grid%xb%v10(i,j)
 381          ! ------------------------------------------------------------------------------
 382
 383       end do
 384    end do
 385
 386    !---------------------------------------------------------------------------
 387    ! Calculate saturation vapour pressure and relative humidity:
 388    !---------------------------------------------------------------------------
 389
 390    do j=js,je
 391       do k=ks,ke
 392          do i=is,ie
 393             call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
 394                grid%xb % q(i,j,k), &
 395                grid%xb %es(i,j,k), grid%xb %qs(i,j,k), grid%xb %rh(i,j,k))
 396          end do
 397       end do
 398    end do
 399
 400    !---------------------------------------------------------------------------
 401    ! Calculate dew point temperature:
 402    !---------------------------------------------------------------------------
 403
 404    call da_trh_to_td (grid)
 405
 406    if (print_detail_xb) then
 407       i=is; j=js; k=ks
 408
 409       write (unit=stdout,fmt='(A,3I5)') 'i,j,k=', i,j,k
 410       write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
 411       write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
 412       write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
 413       write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
 414       write (unit=stdout,fmt='(A)') ' '
 415    end if
 416
 417    !---------------------------------------------------------------------------
 418    ! Sea level pressure and total precipitable water
 419    !---------------------------------------------------------------------------
 420
 421    !---------------------------------------------------------------------------
 422    ! Following code for calculating roughness length needs to be activated
 423    ! for calculating roughness length if sea level pressure is desired
 424    !---------------------------------------------------------------------------
 425
 426    ! call da_wrf_tpq_2_slp (grid)
 427
 428    call da_integrat_dz(grid)
 429    !---------------------------------------------------------------------------
 430    ! Following code for calculating roughness length needs to be activated
 431    ! for surface wind speed for SFC_ASSIM_OPTIONS = 2
 432    !
 433    ! It is not working because KMA terrain height output from wave to grid
 434    ! transform is negative at many grid points
 435    !---------------------------------------------------------------------------
 436
 437    ! tmpvar = log(10.0/0.0001)
 438    ! do j=js,je
 439    !    do i=is,ie
 440    !       if (grid%xb%hf(i,j,ks) <= 0) then
 441    !          write(unit=message(1), fmt=*) ' zero hf at i/j ',i,j,' ks= ',ks
 442    !          call da_error(__FILE__,__LINE__,message(1:1))
 443    !       end if
 444    !       rgh_fac(i,j) = 1.0/log(grid%xb%hf(i,j,ks)/0.0001)
 445
 446    !       grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,ks)*grid%xb%u(i,j,ks) &
 447    !          + grid%xb%v(i,j,ks)*grid%xb%v(i,j,ks) + 1.0e-6) &
 448    !          *tmpvar*rgh_fac(i,j)
 449    !    end do
 450    ! end do
 451
 452    !---------------------------------------------------------------------------
 453    ! Following code for calculating roughness length needs to be activated
 454    ! if SSMI brightness temperature are used
 455    !---------------------------------------------------------------------------
 456
 457    ! call da_transform_xtotb(xa)
 458
 459    !---------------------------------------------------------------------------
 460    ! Calculate means for later use in setting up background errors.
 461    !---------------------------------------------------------------------------
 462
 463    allocate (xbx % latc_mean(jds:jde))
 464
 465    tmpvar = 1.0/real(ide-ids+1)
 466    loc_latc_mean(:) = 0.0
 467
 468    ! Bitwise-exact reduction preserves operation order of serial code for
 469    ! testing, at the cost of much-increased run-time.  Turn it off when not
 470    ! testing.  This will always be .false. for a serial or 1-process MPI run.
 471    if (test_dm_exact) then
 472       allocate(arrayglobal(ids:ide, jds:jde))
 473       call da_patch_to_global(grid, grid%xb%lat, arrayglobal)
 474       do j=jds,jde
 475          loc_latc_mean(j) = tmpvar*sum(arrayglobal(ids:ide, j))
 476       end do
 477       deallocate(arrayglobal)
 478       ! Broadcast result from monitor to other tasks.
 479       call wrf_dm_bcast_real(loc_latc_mean, (jde-jds+1))
 480       xbx % latc_mean = loc_latc_mean
 481    else
 482       do j=js,je
 483          loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(is:ie, j))
 484       end do
 485       call wrf_dm_sum_reals (loc_latc_mean, xbx % latc_mean)
 486    end if
 487
 488    ! WHY?
 489    !  do j=js,je
 490    !    do i=is,ie
 491    !       write(unit=stdout,fmt='(2i4,3f12.2,e12.4,2f12.2)') &
 492    !          j,i,grid%xb%psfc(i,j),grid%xb%tsk(i,j),grid%xb%t2(i,j), &
 493    !          grid%xb%q2(i,j),grid%xb%u10(i,j),grid%xb%v10(i,j)
 494    !    end do
 495    ! end do
 496
 497    if (trace_use) call da_trace_exit("da_transfer_kmatoxb")
 498
 499 end subroutine da_transfer_kmatoxb
 500
 501