var/da/da_transfer_model/da_transfer_wrftoxb.inc

   1 subroutine da_transfer_wrftoxb(xbx, grid, config_flags)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Transfers fields from WRF to first guess structure.
   5    !    Updated for Analysis on Arakawa-C grid
   6    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type (xbx_type), intent(inout)     :: xbx        ! Header & non-gridded vars.
  12
  13    type(domain), intent(inout)        :: grid
  14    type(grid_config_rec_type), intent(in) :: config_flags
  15
  16    integer :: i, j, k, ij
  17
  18    real    :: theta, tmpvar, alt
  19
  20    real, dimension(ims:ime,jms:jme) :: rgh_fac
  21
  22    character(len=19) :: current_date
  23
  24    real :: loc_psac_mean(1)
  25
  26    real, dimension(jds:jde) :: loc_latc_mean
  27
  28    integer :: size2d
  29
  30    real, dimension(kms:kme) :: DDT
  31
  32    real   :: qvf1, cvpm, cpovcv, p_surf, pfd, pfm, phm, pfu, ppb, ttb, albn, aln, height, temp
  33    real, allocatable :: arrayglobal(:,:)
  34
  35    logical:: no_ppb
  36    logical :: has_znt, has_regime
  37
  38    real, dimension(ims:ime,kms:kme) :: pf, pp
  39
  40 #ifdef A2C
  41    real   :: uu, vv
  42 #endif
  43
  44    real              :: zlclg(ids:ide,jds:jde)
  45    real              :: zlcl(ims:ime,jms:jme)
  46
  47
  48    call sfclayinit !for da_sfc_wtq
  49
  50    ! If grid%pb does not existed in FG (YRG, 08/26/2010):
  51      ppb = sum(grid%pb*grid%pb)
  52      no_ppb = ppb == 0.0
  53
  54    !---------------------------------------------------------------------------
  55    ! Set xb array range indices for processor subdomain.
  56    !---------------------------------------------------------------------------
  57
  58    if (trace_use) call da_trace_entry("da_transfer_wrftoxb")
  59
  60    grid%xb % map  = grid%map_proj
  61    grid%xb % ds   = grid%dx
  62
  63    grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
  64    grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
  65    grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
  66
  67    ! WHY?
  68    !rizvi   xbx%big_header%bhi(5,5) = grid%start_year
  69    !rizvi   xbx%big_header%bhi(6,5) = grid%start_month
  70    !rizvi   xbx%big_header%bhi(7,5) = grid%start_day
  71    !rizvi   xbx%big_header%bhi(8,5) = grid%start_hour
  72
  73    !---------------------------------------------------------------------------
  74    ! WRF-specific fields:
  75    !---------------------------------------------------------------------------
  76
  77    ptop = grid%p_top
  78
  79    grid%xb%sigmaf(kte+1) = grid%znw(kte+1)
  80
  81    grid%xb%znw(kte+1) = grid%znw(kte+1)
  82    grid%xb%znu(kte+1) = 0.0
  83
  84    do k=kts,kte
  85       grid%xb%sigmah(k) = grid%znu(k)
  86       grid%xb%sigmaf(k) = grid%znw(k)
  87
  88       grid%xb%znu(k) = grid%znu(k)
  89       grid%xb%znw(k) = grid%znw(k)
  90       grid%xb%dn(k)  = grid%dn(k)
  91       grid%xb%dnw(k) = grid%dnw(k)
  92    end do
  93
  94    grid%xb % ptop = ptop
  95
  96    !---------------------------------------------------------------------------
  97    ! Convert WRF fields to xb:
  98    !---------------------------------------------------------------------------
  99
 100    if (print_detail_xb) then
 101       write(unit=stdout, fmt='(3a, i8)') &
 102          'file:', __FILE__, ', line:', __LINE__
 103
 104       write(unit=stdout, fmt=*) 'its,ite=', its,ite
 105       write(unit=stdout, fmt=*) 'jts,jte=', jts,jte
 106       write(unit=stdout, fmt=*) 'kts,kte=', kts,kte
 107
 108       write(unit=stdout, fmt='(/5a/)') &
 109          'lvl         dnw                dn             rdnw       rdn'
 110
 111       do k=kts,kte+1
 112          write(unit=stdout, fmt='(i3,8f16.8)') k, &
 113             grid%dnw(k), grid%dn(k), grid%rdnw(k), grid%rdn(k)
 114       end do
 115
 116       write(unit=stdout, fmt='(/5a/)') &
 117          'lvl         znu                 znw           rdnw       rdn'
 118
 119       do k=kts,kte+1
 120          write(unit=stdout, fmt='(i3,8f16.8)') k, &
 121             grid%xb%sigmah(k), grid%xb%sigmaf(k), grid%rdnw(k), &
 122             grid%rdn(k)
 123       end do
 124
 125       write(unit=stdout, fmt='(/5a/)') &
 126          'lvl         phb                 ph_2'
 127
 128       do k=kts,kte
 129          write(unit=stdout, fmt='(i3,8e20.12)') k, &
 130                grid%phb(its,jts,k), grid%ph_2(its,jts,k)
 131       end do
 132
 133       write(unit=stdout, fmt=*) 'simple variables:'
 134
 135       if (jte == jde) then
 136          write(unit=stdout, fmt=*) ' '
 137
 138          do k=kts+5,kte,10
 139             do i=its,ite,10
 140                write(unit=stdout, fmt=*) &
 141                     '  grid%v_2(', i, ',', jde+1, ',', k, ')=', &
 142                        grid%v_2(i, jde+1,k)
 143             end do
 144             write(unit=stdout, fmt=*) ' '
 145          end do
 146       end if
 147
 148       if (ite == ide) then
 149          write(unit=stdout, fmt=*) ' '
 150
 151          do k=kts+5,kte,10
 152             do j=jts,jte,10
 153                write(unit=stdout, fmt=*) &
 154                   '  grid%u_2(', ide+1, ',', j, ',', k, ')=', &
 155                   grid%u_2(ide+1,j,k)
 156             end do
 157             write(unit=stdout, fmt=*) ' '
 158          end do
 159       end if
 160
 161       write(unit=stdout, fmt=*) 'simple variables:'
 162
 163       write(unit=stdout,fmt=*) &
 164          '  grid%u_1(its,jts,kts)=', grid%u_1(its,jts,kts)
 165       write(unit=stdout,fmt=*) &
 166          '  grid%v_1(its,jts,kts)=', grid%v_1(its,jts,kts)
 167       write(unit=stdout,fmt=*) &
 168          '  grid%w_1(its,jts,kts)=', grid%w_1(its,jts,kts)
 169       write(unit=stdout,fmt=*) &
 170          '  grid%t_1(its,jts,kts)=', grid%t_1(its,jts,kts)
 171       write(unit=stdout,fmt=*) &
 172          '  grid%ph_1(its,jts,kts)=',grid%ph_1(its,jts,kts)
 173
 174
 175       write(unit=stdout,fmt=*) &
 176          '  grid%u_2(its,jte,kts)=', grid%u_2(its,jte,kts)
 177       write(unit=stdout,fmt=*) &
 178          '  grid%v_2(ite,jts,kts)=', grid%v_2(ite,jts,kts)
 179       write(unit=stdout,fmt=*) &
 180          '  grid%w_2(its,jts,kts)=', grid%w_2(its,jts,kts)
 181       write(unit=stdout,fmt=*) &
 182          '  grid%t_2(its,jts,kts)=', grid%t_2(its,jts,kts)
 183       write(unit=stdout,fmt=*) &
 184          '  grid%ph_2(its,jts,kts)=',grid%ph_2(its,jts,kts)
 185       write(unit=stdout,fmt=*) &
 186          '  grid%phb(its,jts,kts)=', grid%phb(its,jts,kts)
 187
 188       write(unit=stdout, fmt=*) &
 189          '  grid%sm31,grid%em31,grid%sm32,grid%em32, grid%sm33,grid%em33=', &
 190          grid%sm31,grid%em31,grid%sm32,grid%em32,grid%sm33,grid%em33
 191
 192       write(unit=stdout, fmt=*) '  grid%p_top=', grid%p_top
 193       write(unit=stdout, fmt=*) '  grid%znu(kts)=', grid%znu(kts)
 194       write(unit=stdout, fmt=*) '  grid%mub(its,jts)=', grid%mub(its,jts)
 195       write(unit=stdout, fmt=*) '  grid%mu_2(its,jts)=', &
 196          grid%mu_2(its,jts)
 197
 198       write(unit=stdout, fmt=*) '  hbot(its,jts)=', grid%hbot(its,jts)
 199       write(unit=stdout, fmt=*) '  htop(its,jts)=', grid%htop(its,jts)
 200
 201       write(unit=stdout, fmt=*) '  grid%p_top=', grid%p_top
 202       write(unit=stdout, fmt=*) '  num_moist=', num_moist
 203       write(unit=stdout, fmt=*) '  P_QV=', P_QV
 204
 205       write(unit=stdout, fmt=*) '  moist(its,jts,kts,2)=', &
 206          grid%moist(its,jts,kts,2)
 207       write(unit=stdout, fmt=*) ' '
 208    end if
 209
 210    !---------------------------------------------------------------
 211    ! Need this to exchange values in the halo region.
 212    ! grid%xa%u and grid%xa%v are used as temporary arrays and so
 213    ! it is easy to use the existing exchange scheme.
 214    !
 215    ! Note, this is needed as u_2 and v_2 has no guarantee
 216    ! the most east column, and the most north row are
 217    ! properly initailized for each tile.
 218    !---------------------------------------------------------------
 219 #ifdef A2C
 220
 221    grid%xa%u(its:ite+1,jts:jte,kts:kte) = grid%u_2(its:ite+1,jts:jte,kts:kte)
 222    grid%xa%v(its:ite,jts:jte+1,kts:kte) = grid%v_2(its:ite,jts:jte+1,kts:kte)
 223    grid%xb%u(its:ite+1,jts:jte,kts:kte) = grid%u_2(its:ite+1,jts:jte,kts:kte)
 224    grid%xb%v(its:ite,jts:jte+1,kts:kte) = grid%v_2(its:ite,jts:jte+1,kts:kte)
 225
 226 !rizvi's update
 227   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 228        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 229      ipe = ipe + 1
 230      ide = ide + 1
 231   end if
 232
 233   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 234        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 235      jpe = jpe + 1
 236      jde = jde + 1
 237   end if
 238 !rizvi's update
 239 #ifdef DM_PARALLEL
 240 #include "HALO_XB_UV.inc"
 241 #endif
 242
 243 !rizvi's update
 244   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 245        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 246      ipe = ipe - 1
 247      ide = ide - 1
 248   end if
 249
 250   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 251        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 252      jpe = jpe - 1
 253      jde = jde - 1
 254   end if
 255
 256 #else
 257
 258    ! Fill the halo region for u and v.
 259
 260 #ifdef DM_PARALLEL
 261 #include "HALO_EM_C.inc"
 262 #endif
 263
 264 #endif
 265
 266    if (print_detail_xb) then
 267       write(unit=stdout, fmt=*) &
 268          ' ids,ide,jds,jde,kds,kde=', ids,ide,jds,jde,kds,kde
 269       write(unit=stdout, fmt=*) &
 270          ' its,ite,jts,jte,kts,kte=', its,ite,jts,jte,kts,kte
 271       write(unit=stdout, fmt=*) &
 272           ' ims,ime,jms,jme,kms,kme=', ims,ime,jms,jme,kms,kme
 273
 274       write(unit=stdout, fmt=*) &
 275          ' lbound(grid%xb%u)=',   lbound(grid%xb%u)
 276       write(unit=stdout, fmt=*) &
 277          ' lbound(grid%xb%v)=',   lbound(grid%xb%v)
 278       write(unit=stdout, fmt=*) &
 279          ' lbound(grid%u_2)=', lbound(grid%u_2)
 280       write(unit=stdout, fmt=*) &
 281          ' lbound(grid%v_2)=', lbound(grid%v_2)
 282       write(unit=stdout, fmt=*) &
 283          ' ubound(grid%xb%u)=',   ubound(grid%xb%u)
 284       write(unit=stdout, fmt=*) &
 285          ' ubound(grid%xb%v)=',   ubound(grid%xb%v)
 286       write(unit=stdout, fmt=*) &
 287          ' ubound(grid%u_2)=', ubound(grid%u_2)
 288       write(unit=stdout, fmt=*) &
 289          ' ubound(grid%v_2)=', ubound(grid%v_2)
 290    end if
 291
 292    !$OMP PARALLEL DO &
 293    !$OMP PRIVATE ( ij, i, j, k, cvpm, cpovcv, ppb, temp, ttb ) &
 294    !$OMP PRIVATE ( albn, qvf1, aln, theta ) &
 295    !$OMP PRIVATE ( p_surf, pfu, pfd, phm )
 296    do ij = 1 , grid%num_tiles
 297
 298    do j=grid%j_start(ij), grid%j_end(ij)
 299       k = kte+1
 300
 301       do i=its,ite
 302          grid%p(i,j,k) = 0.0
 303          grid%xb%map_factor(i,j) = grid%msft(i,j)
 304          grid%xb%cori(i,j) = grid%f(i,j)
 305          grid%xb%tgrn(i,j) = grid%sst(i,j)
 306          if (grid%xb%tgrn(i,j) < 100.0) &
 307             grid%xb%tgrn(i,j) = grid%tmn(i,j)
 308          grid%xb%lat(i,j) = grid%xlat(i,j)
 309          grid%xb%lon(i,j) = grid%xlong(i,j)
 310          grid%xb%terr(i,j) = grid%ht(i,j)
 311          grid%xb%snow(i,j) = grid%snowc(i,j)
 312          grid%xb%lanu(i,j) = grid%lu_index(i,j)
 313          grid%xb%landmask(i,j) = grid%landmask(i,j)
 314          grid%xb%xland(i,j) = grid%xland(i,j)
 315          ! Z. Liu below are variables used by RTTOV
 316          grid%xb%tsk(i,j) = grid%tsk(i,j)
 317          grid%xb%smois(i,j) = grid%smois(i,j,1)
 318          grid%xb%tslb(i,j) = grid%tslb(i,j,1)
 319          grid%xb%xice(i,j) = grid%xice(i,j)
 320          grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
 321          grid%xb%isltyp(i,j) = grid%isltyp(i,j)
 322          grid%xb%vegfra(i,j) = grid%vegfra(i,j)
 323          grid%xb%snowh(i,j) = grid%snowh(i,j)*1000.0 ! meter to mm
 324          if ( grid%xb%ivgtyp(i,j) == grid%iswater .and. &
 325               grid%xb%snow(i,j) > 0.0 )then
 326             grid%xb%snow(i,j) = 0.0
 327             grid%xb%snowh(i,j) = 0.0
 328          end if
 329       end do
 330    end do
 331
 332       ! WHY?
 333       ! Adapted the code from "real.init.code" by Y.-R. Guo 05/13/2004:
 334
 335       ! do i=its,ite
 336       !    k = kte
 337       !    qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV))
 338       !    qvf2 = 1.0/(1.0+qvf1)
 339       !    qvf1 = qvf1*qvf2
 340       !    grid%xb%p(i,j,k) = -0.5*(grid%mu_2(i,j)+qvf1* &
 341       !       grid%mub(i,j))/grid%rdnw(k)/qvf2
 342
 343       !    do k = kte-1,1,-1
 344       !       qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV))
 345       !       qvf2 = 1.0/(1.0+qvf1)
 346       !       qvf1 = qvf1*qvf2
 347       !       grid%p(i,j,k) = grid%p(i,j,k+1) - &
 348       !          (grid%mu_2(i,j)+qvf1*grid%mub(i,j))/qvf2/rdn(k+1)
 349       !    end do
 350       ! end do
 351
 352       ! Adapted the code from WRF module_big_step_utilitites_em.F ----
 353       !         subroutine calc_p_rho_phi      Y.-R. Guo (10/20/2004)
 354
 355       ! NOTE: as of V3.1, P (pressure perturbation) and PB (base state pressure)
 356       ! are included in the wrfinput file. However, P and PB are still
 357       ! re-calculated here.
 358
 359       ! NOTE: as of 7/01/2010, P and PB are directly from the first guess
 360       ! and no longer re-calculated here.
 361
 362       cvpm =  - (1.0 - gas_constant/cp)
 363       cpovcv = cp / (cp - gas_constant)
 364
 365       ! In case of var4d, pb etc. will be recalculated in start_em with realsize=8,
 366       ! However, the originals are computed with realsize=4.
 367       if ( no_ppb ) then
 368          do k=kts,kte
 369           do j=grid%j_start(ij), grid%j_end(ij)
 370             do i=its,ite
 371                ! The base specific volume (from real.init.code)
 372                ppb  = c3h(k) * grid%mub(i,j) + c4h(k) + ptop
 373                grid%pb(i,j,k) = ppb
 374                temp = MAX ( iso_temp, base_temp + base_lapse*log(ppb/base_pres) )
 375                if ( grid%pb(i,j,k) < base_pres_strat ) then
 376                   temp = iso_temp + base_lapse_strat*log(grid%pb(i,j,k)/base_pres_strat)
 377                end if
 378                ttb  = temp * (base_pres/ppb)**kappa
 379                ! ttb  = (base_temp + base_lapse*log(ppb/base_pres)) * &
 380                !   (base_pres/ppb)**kappa
 381                albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm
 382
 383                qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
 384                aln  = -1.0 / (c1h(k)*(grid%mub(i,j)+grid%mu_2(i,j))+c2h(k)) * &
 385                       (albn*(c1h(k)*grid%mu_2(i,j)) + grid%rdnw(k) * &
 386                       (grid%ph_2(i,j,k+1) - grid%ph_2(i,j,k)))
 387                ! total pressure:
 388                grid%xb%p(i,j,k) = base_pres * &
 389                                  ((gas_constant*(t0+grid%t_2(i,j,k))*qvf1) / &
 390                                  (base_pres*(aln+albn)))**cpovcv
 391                ! total density
 392                grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
 393                ! pressure perturbation:
 394                grid%p(i,j,k) = grid%xb%p(i,j,k) - ppb
 395             end do
 396            end do
 397          end do
 398       else
 399          do j=grid%j_start(ij), grid%j_end(ij)
 400           do i=its,ite
 401             p_surf = base_pres * EXP ( -base_temp/base_lapse + ( (base_temp/base_lapse)**2 - 2.*gravity*grid%ht(i,j)/base_lapse/gas_constant ) **0.5 )
 402             do k = kts, kte
 403                grid%pb(i,j,k) = c3h(k)*(p_surf - grid%p_top) + c4h(k) + grid%p_top
 404                temp = MAX ( iso_temp, base_temp + base_lapse*log(grid%pb(i,j,k)/base_pres) )
 405                if ( grid%pb(i,j,k) < base_pres_strat ) then
 406                   temp = iso_temp + base_lapse_strat*log(grid%pb(i,j,k)/base_pres_strat)
 407                end if
 408                grid%t_init(i,j,k) = temp*(base_pres/grid%pb(i,j,k))**(gas_constant/cp) - t0
 409                grid%alb(i,j,k) = (gas_constant/base_pres)*(grid%t_init(i,j,k)+t0)*(grid%pb(i,j,k)/base_pres)**cvpm
 410             end do
 411             grid%mub(i,j) = p_surf - grid%p_top
 412             grid%phb(i,j,1) = grid%ht(i,j) * gravity
 413             if (grid%hypsometric_opt == 1) then
 414                do k = kts,kte
 415                   grid%phb(i,j,k+1) = grid%phb(i,j,k) - grid%dnw(k)*(c1h(k)*grid%mub(i,j)+c2h(k))*grid%alb(i,j,k)
 416                end do
 417             else if (grid%hypsometric_opt == 2) then
 418                do k = kts,kte
 419                   pfu = grid%mub(i,j)*c3f(k+1) + c4f(k+1) + grid%p_top
 420                   pfd = grid%mub(i,j)*c3f(k  ) + c4f(k  ) + grid%p_top
 421                   phm = grid%mub(i,j)*c3h(k  ) + c4h(k  ) + grid%p_top
 422                   grid%phb(i,j,k+1) = grid%phb(i,j,k) + grid%alb(i,j,k)*phm*LOG(pfd/pfu)
 423                end do
 424             end if
 425           end do
 426          end do
 427
 428          if (grid%hypsometric_opt == 1) then
 429            do k=kts,kte
 430              do j=grid%j_start(ij), grid%j_end(ij)
 431                do i=its,ite
 432                 grid%al(i,j,k)=-1./(c1h(k)*(grid%mub(i,j)+grid%mu_2(i,j))+c2h(k))* &
 433                      (grid%alb(i,j,k)*(c1h(k)*grid%mu_2(i,j))  &
 434                      +grid%rdnw(k)*(grid%ph_2(i,j,k+1)-grid%ph_2(i,j,k)))
 435                 ! total density
 436                 grid%xb%rho(i,j,k)= 1.0 / (grid%alb(i,j,k)+grid%al(i,j,k))
 437               end do
 438             end do
 439            end do
 440          elseif (grid%hypsometric_opt == 2) then
 441            do k=kts,kte
 442              do j=grid%j_start(ij), grid%j_end(ij)
 443               do i=its,ite
 444                pfu = (grid%mub(i,j)+grid%mu_2(i,j))*c3f(k+1)+c4f(k+1)+grid%p_top
 445                pfd = (grid%mub(i,j)+grid%mu_2(i,j))*c3f(k)  +c4f(k)  +grid%p_top
 446                phm = (grid%mub(i,j)+grid%mu_2(i,j))*c3h(k)  +c4h(k)  +grid%p_top
 447                grid%al(i,j,k) = (grid%ph_2(i,j,k+1)-grid%ph_2(i,j,k)+grid%phb(i,j,k+1)-grid%phb(i,j,k)) &
 448                                  /phm/LOG(pfd/pfu)-grid%alb(i,j,k)
 449                ! total density
 450                grid%xb%rho(i,j,k)= 1.0 / (grid%alb(i,j,k)+grid%al(i,j,k))
 451             end do
 452            end do
 453           end do
 454          endif
 455          do k=kts,kte
 456            do j=grid%j_start(ij), grid%j_end(ij)
 457             do i=its,ite
 458               qvf1 = 1.+grid%moist(i,j,k,P_QV) / rd_over_rv
 459               grid%xb%p(i,j,k)=base_pres*( (gas_constant*(t0+grid%t_2(i,j,k))*qvf1)/       &
 460                          (base_pres*(grid%al(i,j,k)+grid%alb(i,j,k))) )**cpovcv
 461             end do
 462            end do
 463          end do
 464       endif
 465
 466       do k=kts,kte+1
 467         do j=grid%j_start(ij), grid%j_end(ij)
 468          do i=its,ite
 469             grid%xb%hf(i,j,k) = (grid%phb(i,j,k)+grid%ph_2(i,j,k))/gravity
 470             grid%xb%w (i,j,k) = grid%w_2(i,j,k)
 471          end do
 472         end do
 473       end do
 474
 475       do j=grid%j_start(ij), grid%j_end(ij)
 476
 477       if (grid%hypsometric_opt == 2) then
 478         ! Compute full-level pressure
 479         ! The full and half level dry hydrostatic pressure is easily computed (YRG, 12/20/2011):
 480         do i=its,ite
 481            do k=kts, kte+1
 482               pf(i,k) = (grid%mub(i,j)+grid%mu_2(i,j)) * c3f(k) + c4f(k) + ptop
 483               pp(i,k) = (grid%mub(i,j)+grid%mu_2(i,j)) * c3h(k) + c4h(k) + ptop
 484            enddo
 485         enddo
 486       endif
 487
 488       do k=kts,kte
 489          do i=its,ite
 490             grid%xb%u(i,j,k) = 0.5*(grid%u_2(i,j,k)+grid%u_2(i+1,j,k))
 491             grid%xb%v(i,j,k) = 0.5*(grid%v_2(i,j,k)+grid%v_2(i,j+1,k))
 492             grid%xb%wh(i,j,k)= 0.5*(grid%xb%w(i,j,k)+grid%xb%w(i,j,k+1))
 493             if (grid%hypsometric_opt == 1) then
 494                grid%xb%h(i,j,k) = 0.5*(grid%xb%hf(i,j,k)+grid%xb%hf(i,j,k+1))
 495             elseif (grid%hypsometric_opt == 2) then
 496             !  This is almost correct for pressure, but not for height.
 497             !  Arithmetic mean of full-level heights is always higher than the actual,
 498             !  leading to a biased model for height-based observations (e.g., GPS RO)
 499             !  and surface variables (2-meter TQ and 10-meter wind).
 500             !  wee 11/22/2011
 501
 502                grid%xb%h(i,j,k) = grid%xb%hf(i,j,k)+(grid%xb%hf(i,j,k+1)-grid%xb%hf(i,j,k)) &
 503                     *LOG(pf(i,k)/pp(i,k))/LOG(pf(i,k)/pf(i,k+1))
 504             endif
 505
 506             if ( num_pseudo == 0 ) then
 507                grid%moist(i,j,k,P_QV) = max(grid%moist(i,j,k,P_QV), qlimit)
 508             end if
 509             ! water vapor mixing ratio (kg h2o / kg dry air) from wrfinput file
 510             grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
 511
 512             theta = t0 + grid%t_2(i,j,k)
 513             grid%xb%t(i,j,k) = theta*(grid%xb%p(i,j,k)/base_pres)**kappa
 514
 515             ! Convert to specific humidity (kg h2o / kg moist air) from mixing ratio:
 516             grid%xb%q(i,j,k)=grid%xb%q(i,j,k)/(1.0+grid%xb%q(i,j,k))
 517
 518             ! Background qrn needed for radar radial velocity assimilation:
 519
 520             if (size(grid%moist,dim=4) >= 4) then
 521                grid%xb%qcw(i,j,k) = max(grid%moist(i,j,k,p_qc), 0.0)
 522                grid%xb%qrn(i,j,k) = max(grid%moist(i,j,k,p_qr), 0.0)
 523 !rizvi For doing single obs test with radiance one need to ensure non-zero hydrometeor
 524 !  For doing so uncomment following two cards below
 525 !               grid%xb%qcw(i,j,k) = 0.0001
 526 !               grid%xb%qrn(i,j,k) = 0.0001
 527
 528                grid%xb%qt (i,j,k) = grid%xb%q(i,j,k) + grid%xb%qcw(i,j,k) + &
 529                   grid%xb%qrn(i,j,k)
 530             end if
 531
 532             if (size(grid%moist,dim=4) >= 6) then
 533                grid%xb%qci(i,j,k) = max(grid%moist(i,j,k,p_qi), 0.0)
 534                grid%xb%qsn(i,j,k) = max(grid%moist(i,j,k,p_qs), 0.0)
 535 !rizvi For doing single obs test with radiance one need to ensure non-zero hydrometeor
 536 !  For doing so uncomment following two cards below
 537 !               grid%xb%qci(i,j,k) = 0.0001
 538 !               grid%xb%qsn(i,j,k) = 0.0001
 539             end if
 540
 541             if (size(grid%moist,dim=4) >= 7) then
 542                grid%xb%qgr(i,j,k) = max(grid%moist(i,j,k,p_qg), 0.0)
 543             end if
 544
 545             if ( grid%mp_physics == 3 ) then   ! WSM3-class scheme
 546                if ( grid%xb%t(i,j,k) <= t_kelvin ) then
 547                   grid%xb%qci(i,j,k) = grid%xb%qcw(i,j,k)
 548                   grid%xb%qcw(i,j,k) = 0.0
 549                   grid%xb%qsn(i,j,k) = grid%xb%qrn(i,j,k)
 550                   grid%xb%qrn(i,j,k) = 0.0
 551                end if
 552             end if
 553
 554          end do
 555       end do
 556    end do
 557
 558    do j=grid%j_start(ij), grid%j_end(ij)
 559       do i=its,ite
 560          grid%xb%psac(i,j) = grid%mub(i,j)+grid%mu_2(i,j)
 561 ! To make the Psfc consistent with WRF (YRG, 04/06/2010):
 562          grid%xb%psfc(i,j) = grid%psfc(i,j)
 563
 564          if (grid%xb%tgrn(i,j) < 100.0) &
 565             grid%xb%tgrn(i,j) = grid%xb%t(i,j,kts)+ &
 566             0.0065*(grid%xb%h(i,j,kts)-grid%xb%hf(i,j,kts))
 567       end do
 568    end do
 569
 570    end do
 571    !$OMP END PARALLEL DO
 572
 573 !   write (999,'("MABS=",e13.5)') sum(abs(grid%xb%psfc(:,:)-grid%psfc(:,:))) / &
 574 !                                       float((ite-its+1)*(jte-jts+1))
 575
 576    grid%xb%ztop = grid%xb%hf(its,jts,kte+1)
 577
 578    if (print_detail_xb) then
 579       write(unit=stdout, fmt=*) ' '
 580       if (print_detail_xb) then
 581          write(unit=stdout, fmt='(/5a/)') &
 582             'lvl         h                 p                t'
 583
 584          do k=kts,kte
 585             write(unit=stdout, fmt='(i3,8e20.12)') k, &
 586                grid%xb%h(its,jts,k), grid%xb%p(its,jts,k), grid%xb%t(its,jts,k)
 587          end do
 588       end if
 589
 590       write(unit=stdout,fmt=*) ' '
 591       write(unit=stdout,fmt=*) 'grid%xb%u(its,jte,kte)=', grid%xb%u(its,jte,kte)
 592       write(unit=stdout,fmt=*) 'grid%xb%v(ite,jts,kte)=', grid%xb%v(ite,jts,kte)
 593       write(unit=stdout,fmt=*) 'grid%xb%w(its,jts,kte)=', grid%xb%w(its,jts,kte)
 594       write(unit=stdout,fmt=*) 'grid%xb%t(its,jts,kte)=', grid%xb%t(its,jts,kte)
 595       write(unit=stdout,fmt=*) 'grid%xb%p(its,jts,kte)=', grid%xb%p(its,jts,kte)
 596       write(unit=stdout,fmt=*) 'grid%xb%q(its,jts,kte)=', grid%xb%q(its,jts,kte)
 597       write(unit=stdout,fmt=*) 'grid%xb%h(its,jts,kte)=', grid%xb%h(its,jts,kte)
 598       write(unit=stdout,fmt=*) &
 599          'grid%xb%hf(its,jts,kte)=', grid%xb%hf(its,jts,kte)
 600       write(unit=stdout,fmt=*) &
 601          'grid%xb%map_factor(its,jts)=', grid%xb%map_factor(its,jts)
 602       write(unit=stdout,fmt=*) 'grid%xb%cori(its,jts)=', grid%xb%cori(its,jts)
 603       write(unit=stdout,fmt=*) 'grid%xb%tgrn(its,jts)=', grid%xb%tgrn(its,jts)
 604       write(unit=stdout,fmt=*) 'grid%xb%lat(its,jts)=', grid%xb%lat(its,jts)
 605       write(unit=stdout,fmt=*) 'grid%xb%lon(its,jts)=', grid%xb%lon(its,jts)
 606       write(unit=stdout,fmt=*) 'grid%xb%terr(its,jts)=', grid%xb%terr(its,jts)
 607       write(unit=stdout,fmt=*) 'grid%xb%snow(its,jts)=', grid%xb%snow(its,jts)
 608       write(unit=stdout,fmt=*) 'grid%xb%lanu(its,jts)=', grid%xb%lanu(its,jts)
 609       write(unit=stdout,fmt=*) &
 610          'grid%xb%landmask(its,jts)=', grid%xb%landmask(its,jts)
 611       write(unit=stdout,fmt=*) '(ite,jte)=', ite,jte
 612       write(unit=stdout,fmt=*) 'grid%xb%lat(ite,jte)=', grid%xb%lat(ite,jte)
 613       write(unit=stdout,fmt=*) 'grid%xb%lon(ite,jte)=', grid%xb%lon(ite,jte)
 614       write(unit=stdout,fmt=*) ' '
 615    end if
 616
 617    !---------------------------------------------------------------------------
 618    ! [3.0] Calculate vertical inner product for use in vertical transform:
 619    !---------------------------------------------------------------------------
 620
 621    if (vertical_ip == vertical_ip_sqrt_delta_p) then
 622       ! Vertical inner product is sqrt(Delta p):
 623       do k=kts,kte
 624          grid%xb % vertical_inner_product(its:ite,jts:jte,k) = &
 625             sqrt(grid%xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k))
 626       end do
 627    else if (vertical_ip == vertical_ip_delta_p) then
 628
 629       ! Vertical inner product is Delta p:
 630       do k=1,grid%xb%mkz
 631          grid % xb % vertical_inner_product(its:ite,jts:jte,k) = &
 632          grid % xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k)
 633       end do
 634    end if
 635
 636    !---------------------------------------------------------------------------
 637    ! Roughness
 638    !---------------------------------------------------------------------------
 639
 640    current_date = 'yyyy-mm-dd_hh:mm:ss'
 641
 642    write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
 643       grid%start_year, '-', &
 644       grid%start_month, '-', &
 645       grid%start_day, '_', &
 646       grid%start_hour, ':', &
 647       grid%start_minute, ':', &
 648       grid%start_second
 649
 650    !xbx % mminlu = 'USGS'
 651    xbx % mminlu = trim(grid%mminlu)
 652
 653    has_regime = sum(grid%regime*grid%regime) > 0.0
 654    has_znt    = sum(grid%znt*grid%znt)       > 0.0
 655    if ( has_regime .neqv. has_znt ) then
 656       ! make sure they are consistent
 657       ! either both from model, or both calculated here
 658       has_regime = .false.
 659       has_znt    = .false.
 660    end if
 661    if ( has_znt ) then
 662       grid%xb%rough = grid%znt
 663    else
 664       ! calculate rough only when it is not available from fg
 665       call da_roughness_from_lanu(19, xbx % mminlu, current_date, &
 666          grid%xb % lanu, grid%xb % rough)
 667    end if
 668
 669    !---------------------------------------------------------------------------
 670    ! Calculate 1/grid box areas:
 671    !---------------------------------------------------------------------------
 672
 673    if (print_detail_xb) then
 674       write(unit=stdout, fmt='(/a, e24.16)') &
 675          'grid%xb % ds=', grid%xb % ds
 676
 677       write(unit=stdout, fmt='(a, e24.16/)') &
 678            'grid%xb % map_factor(its,jts)=', grid%xb % map_factor(its,jts)
 679    end if
 680
 681    !$OMP PARALLEL DO &
 682 #ifndef A2C
 683    !$OMP PRIVATE ( ij, i, j, tmpvar, height, message )
 684 #else
 685    !$OMP PRIVATE ( ij, i, j, tmpvar, height, message, uu, vv )
 686 #endif
 687    do ij = 1 , grid%num_tiles
 688
 689    do j=grid%j_start(ij),grid%j_end(ij)
 690       do i=its,ite
 691          if (grid%xb%ztop < grid%xb%hf(i,j,kte+1)) &
 692              grid%xb%ztop = grid%xb%hf(i,j,kte+1)
 693
 694          tmpvar = grid%xb%ds / grid%xb%map_factor(i,j)
 695
 696          grid%xb % grid_box_area(i,j) = tmpvar*tmpvar
 697
 698          ! Calculate surface variable(wind, moisture, temperature)
 699          ! sfc variables: 10-m wind, and 2-m T, Q, at cross points
 700
 701          height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
 702          if (height <= 0.0) then
 703             message(1) = "Negative height found"
 704             write (unit=message(2),FMT='(2I6,A,F10.2,A,F10.2)') &
 705                i,j,' ht = ',grid%xb%h(i,j,kts) ,' terr =  ',grid%xb%terr(i,j)
 706             call da_error(__FILE__,__LINE__, message(1:2))
 707          end if
 708
 709 #ifdef A2C
 710          uu = 0.5*(grid%xb%u(i,j,kts)+grid%xb%u(i+1,j,kts) )
 711          vv = 0.5*(grid%xb%v(i,j,kts)+grid%xb%v(i,j+1,kts) )
 712 #endif
 713          if ( has_regime ) then
 714             ! if fg contains valid regime info, use it
 715             grid%xb%regime(i,j) = grid%regime(i,j)
 716          else
 717             ! xb t2/q2/u10/v10 will be coming directly from fg.
 718             ! but still call da_sfc_wtq here in order to get regimes.
 719             ! regimes can not be zero for sfc_assi_options=2.
 720             ! regimes calculated here could be very different from WRF values
 721             ! due to the lack of some input sfc info
 722             call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
 723 #ifdef A2C
 724                grid%xb%p(i,j,kts), grid%xb%t(i,j,kts), grid%xb%q(i,j,kts),uu,vv, &
 725 #else
 726                grid%xb%p(i,j,kts), grid%xb%t(i,j,kts), grid%xb%q(i,j,kts), &
 727                grid%xb%u(i,j,kts), grid%xb%v(i,j,kts), &
 728 #endif
 729                height,  grid%xb%rough(i,j),grid%xb%xland(i,j), grid%xb%ds, &
 730                grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), &
 731                grid%xb%q2(i,j), grid%xb%regime(i,j))
 732          end if !if has_regime
 733
 734          ! use t2/q2/u10/v10 from fg
 735          grid%xb%u10(i,j) = grid%u10(i,j)
 736          grid%xb%v10(i,j) = grid%v10(i,j)
 737          grid%xb%t2(i,j) = grid%t2(i,j)
 738          grid%xb%q2(i,j) = grid%q2(i,j)
 739
 740       end do
 741    end do
 742
 743    end do
 744    !$OMP END PARALLEL DO
 745
 746 #ifdef VAR4D
 747    CALL set_physical_bc2d( grid%xb%grid_box_area, 't', config_flags, &
 748                            ids, ide, jds, jde,    &
 749                            ims, ime, jms, jme,    &
 750                            ips, ipe, jps, jpe,    &
 751                            its, ite, jts, jte    )
 752 #endif
 753
 754    !---------------------------------------------------------------------------
 755    ! Calculate saturation vapour pressure and relative humidity:
 756    !---------------------------------------------------------------------------
 757
 758    !$OMP PARALLEL DO &
 759    !$OMP PRIVATE ( ij, k, j, i )
 760    do ij = 1 , grid%num_tiles
 761       do k=kts,kte
 762          do j=grid%j_start(ij),grid%j_end(ij)
 763             do i=its,ite
 764                call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
 765                   grid%xb % q(i,j,k), grid%xb %es(i,j,k), grid%xb %qs(i,j,k), &
 766                   grid%xb %rh(i,j,k))
 767             end do
 768          end do
 769       end do
 770    end do
 771    !$OMP END PARALLEL DO
 772
 773    !---------------------------------------------------------------------------
 774    ! Calculate dew point temperature:
 775    !---------------------------------------------------------------------------
 776
 777    call da_trh_to_td (grid)
 778
 779    if (print_detail_xb) then
 780       i=its; j=jts; k=kts
 781
 782       write(unit=stdout, fmt=*) 'i,j,k=', i,j,k
 783       write(unit=stdout, fmt=*) 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
 784       write(unit=stdout, fmt=*) 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
 785       write(unit=stdout, fmt=*) 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
 786       write(unit=stdout, fmt=*) 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
 787       write(unit=stdout, fmt=*) ' '
 788    end if
 789
 790    !---------------------------------------------------------------------------
 791    ! Sea level pressure and total precipitable water
 792    !---------------------------------------------------------------------------
 793
 794    call da_wrf_tpq_2_slp (grid)
 795
 796    ! WHY?
 797    ! do j = jts,jte
 798    !    do i = its,ite
 799    !       call da_tpq_to_slp(grid%xb%t(i,j,:), grid%xb%q(i,j,:), &
 800    !          grid%xb%p(i,j,:), grid%xb%terr(i,j), &
 801    !          grid%xb%psfc(i,j), grid%xb%slp(i,j), grid%xp)
 802    !    end do
 803    ! end do
 804
 805    call da_integrat_dz(grid)
 806
 807    !---------------------------------------------------------------------------
 808    ! Surface wind speed
 809    !---------------------------------------------------------------------------
 810
 811    tmpvar = log(10.0/0.0001)
 812
 813    !$OMP PARALLEL DO &
 814 #ifndef A2C
 815    !$OMP PRIVATE (ij, i, j, height)
 816 #else
 817    !$OMP PRIVATE (ij, i, j, height, uu, vv)
 818 #endif
 819    do ij = 1, grid%num_tiles
 820
 821    do j=grid%j_start(ij), grid%j_end(ij)
 822       do i=its,ite
 823          height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
 824          rgh_fac(i,j) = 1.0/log(height/0.0001)
 825 #ifndef A2C
 826          grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,kts)*grid%xb%u(i,j,kts) &
 827                          + grid%xb%v(i,j,kts)*grid%xb%v(i,j,kts) + 1.0e-6) &
 828                     *tmpvar*rgh_fac(i,j)
 829 #else
 830          uu = 0.5*(grid%xb%u(i,j,kts)+grid%xb%u(i+1,j,kts) )
 831          vv = 0.5*(grid%xb%v(i,j,kts)+grid%xb%v(i,j+1,kts) )
 832          grid%xb%speed(i,j) = sqrt(uu*uu + vv*vv + 1.0e-6)*tmpvar*rgh_fac(i,j)
 833 #endif
 834       end do
 835    end do
 836
 837    end do
 838    !$OMP END PARALLEL DO
 839
 840    !---------------------------------------------------------------------------
 841    ! Brightness temperature SH Chen
 842    !---------------------------------------------------------------------------
 843
 844    if (use_ssmitbobs)   &
 845       call da_transform_xtotb(grid)
 846
 847    !---------------------------------------------------------------------------
 848    ! GPS Refractivity linked by Y.-R. Guo 05/28/2004
 849    !---------------------------------------------------------------------------
 850
 851    call da_transform_xtogpsref(grid)
 852
 853    !---------------------------------------------------------------------------
 854    ! Ground-based GPS ZTD must follow the GPS Refractivity calculation.
 855    !---------------------------------------------------------------------------
 856
 857    ! WHY? For certain computation method, not current one.
 858    if (use_gpsztdobs) then
 859       call da_transform_xtoztd(grid)
 860       if (print_detail_xb) then
 861         i=its; j=jts
 862         write(unit=stdout, fmt=*) 'grid%xb % tpw(i,j)=', grid%xb % tpw(i,j)
 863         write(unit=stdout, fmt=*) 'grid%xb % ztd(i,j)=', grid%xb % ztd(i,j)
 864         write(unit=stdout, fmt=*) ' '
 865       end if
 866    end if
 867
 868    !---------------------------------------------------------------------------
 869    ! Calculate means for later use in setting up background errors.
 870    !---------------------------------------------------------------------------
 871
 872    ! WHY?
 873    ! if (.not. associated(xbx % latc_mean)) then
 874    allocate (xbx % latc_mean(jds:jde))
 875    if (trace_use) call da_trace("da_transfer_wrftoxb",&
 876       message="allocated xbx%latc_mean")
 877    ! end if
 878
 879    size2d = (ide-ids+1)*(jde-jds+1)
 880
 881    tmpvar = 1.0/real(size2d)
 882
 883    ! Bitwitse-exact reduction preserves operation order of serial code for
 884    ! testing, at the cost of much-increased run-time.  Turn it off when not
 885    ! testing.  Thits will always be .false. for a serial or 1-process MPI run.
 886    if (test_dm_exact) then
 887       allocate(arrayglobal(ids:ide, jds:jde))
 888       call da_patch_to_global(grid,grid%xb%psac, arrayglobal)
 889       loc_psac_mean = tmpvar*sum(arrayglobal(ids:ide,jds:jde))
 890       deallocate(arrayglobal)
 891    else
 892       loc_latc_mean = 0.0
 893       loc_psac_mean = tmpvar*sum(grid%xb % psac(its:ite,jts:jte))
 894    end if
 895
 896    tmpvar = 1.0/real(ide-ids+1)
 897
 898    if (test_dm_exact) then
 899       allocate(arrayglobal(ids:ide, jds:jde))
 900       call da_patch_to_global(grid,grid%xb%lat, arrayglobal)
 901       do j=jds,jde
 902          loc_latc_mean(j) = tmpvar*sum(arrayglobal(ids:ide, j))
 903       end do
 904       deallocate(arrayglobal)
 905    else
 906       loc_latc_mean = 0.0
 907       do j=jts,jte
 908          loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(its:ite, j))
 909       end do
 910    end if
 911
 912    if (test_dm_exact) then
 913       ! Broadcast result from monitor to other tasks.
 914       call wrf_dm_bcast_real(loc_psac_mean, 1)
 915       xbx % psac_mean = loc_psac_mean(1)
 916       ! Broadcast result from monitor to other tasks.
 917       call wrf_dm_bcast_real(loc_latc_mean, (jde-jds+1))
 918       xbx % latc_mean = loc_latc_mean
 919    else
 920       xbx % psac_mean = wrf_dm_sum_real(loc_psac_mean(1))
 921       call wrf_dm_sum_reals(loc_latc_mean, xbx % latc_mean)
 922    end if
 923
 924    if (print_detail_xb) then
 925       ! write(unit=stdout, fmt=*) 'loc_psac_mean  =', loc_psac_mean
 926       write(unit=stdout, fmt=*) 'xbx % psac_mean=', xbx % psac_mean
 927
 928       ! write(unit=stdout, fmt=*) 'loc_latc_mean  =', loc_latc_mean(jts)
 929       write(unit=stdout, fmt=*) 'xbx % latc_mean=', xbx % latc_mean(jts)
 930    end if
 931
 932
 933    ! Fill the halo region for xb
 934
 935 #ifdef DM_PARALLEL
 936 #include "HALO_XB.inc"
 937 #endif
 938
 939    ! Calculate time step from one dimensional cloud model parameterization
 940
 941    if (dt_cloud_model) then
 942       do j = jts, jte
 943          do i = its, ite
 944             call da_cloud_model (grid%xb%t(I,J,:),  grid%xb%p(I,J,:), &
 945                grid%xb%q(I,J,:), grid%xb%qcw(I,J,:), grid%xb%qrn(I,J,:), &
 946                grid%xb%h(I,J,:), grid%xb%hf(I,J,:), ddt, kts, kte)
 947
 948             do k = kts, kte
 949                grid%xb%delt(i,j,k) = DDT(k)
 950             end do
 951          end do
 952       end do
 953    end if
 954
 955    deallocate (xbx % latc_mean)
 956
 957    ! calculate background/model LCL to be used by use_radar_rqv
 958    if ( use_radarobs .and. use_radar_rqv .and. cloudbase_calc_opt == 2 ) then
 959       do j = jts, jte
 960          do i = its, ite
 961             zlcl(i,j) = 125.0*(grid%xb%t(i,j,1)-grid%xb%td(i,j,1))
 962          end do
 963       end do
 964       ij = (ide-ids+1)*(jde-jds+1)
 965 #ifdef DM_PARALLEL
 966       call da_patch_to_global (grid, zlcl, zlclg)
 967       call wrf_dm_bcast_real( zlclg, ij )
 968 #else
 969       do j = jds, jde
 970          do i = ids, ide
 971             zlclg(i,j) = zlcl(i,j)
 972          end do
 973       end do
 974 #endif
 975       ! Calculate the domain grid mean LCL
 976       zlcl_mean = 0.0
 977       do j = jds, jde
 978          do i = ids, ide
 979             zlcl_mean = zlcl_mean + zlclg(i,j)
 980          end do
 981       end do
 982       zlcl_mean = zlcl_mean/float(ij)
 983    end if ! lcl for use_radar_rqv
 984
 985    if ( use_gpsephobs ) then
 986       call da_gpseph_init(grid)
 987    end if
 988
 989    if (trace_use) call da_trace_exit("da_transfer_wrftoxb")
 990
 991 end subroutine da_transfer_wrftoxb
 992