var/da/da_transfer_model/da_transfer_wrftoxb_lite.inc

   1 subroutine da_transfer_wrftoxb_lite(xbx, grid, config_flags)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Transfers fields from WRF to first guess structure.
   5    !    Author: Xin Zhang,  MMM/NESL/NCAR,  Date: 10/10/2011
   6    !---------------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type (xbx_type), intent(inout)     :: xbx        ! Header & non-gridded vars.
  11
  12    type(domain), intent(inout)        :: grid
  13    type(grid_config_rec_type), intent(in) :: config_flags
  14
  15    integer :: i, j, k, ij
  16
  17    real    :: theta, tmpvar
  18
  19    real, dimension(ims:ime,jms:jme) :: rgh_fac
  20
  21    character(len=19) :: current_date
  22
  23    real :: loc_psac_mean
  24
  25    real, dimension(jds:jde) :: loc_latc_mean
  26
  27    integer :: size2d
  28
  29    real, dimension(kms:kme) :: DDT
  30
  31    real   :: qvf1, cvpm, cpovcv, ppb, ttb, albn, aln, height, temp
  32    real, allocatable :: arrayglobal(:,:)
  33
  34    logical:: no_ppb
  35
  36
  37    ! If grid%pb does not existed in FG (YRG, 08/26/2010):
  38      ppb = sum(grid%pb*grid%pb)
  39      no_ppb = ppb == 0.0
  40
  41    !---------------------------------------------------------------------------
  42    ! Set xb array range indices for processor subdomain.
  43    !---------------------------------------------------------------------------
  44
  45    if (trace_use) call da_trace_entry("da_transfer_wrftoxb_lite")
  46
  47    grid%xb % map  = grid%map_proj
  48    grid%xb % ds   = grid%dx
  49
  50    grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
  51    grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
  52    grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
  53
  54    !---------------------------------------------------------------------------
  55    ! WRF-specific fitelds:
  56    !---------------------------------------------------------------------------
  57
  58    ptop = grid%p_top
  59
  60    grid%xb%sigmaf(kte+1) = grid%znw(kte+1)
  61
  62    grid%xb%znw(kte+1) = grid%znw(kte+1)
  63    grid%xb%znu(kte+1) = 0.0
  64
  65    do k=kts,kte
  66       grid%xb%sigmah(k) = grid%znu(k)
  67       grid%xb%sigmaf(k) = grid%znw(k)
  68
  69       grid%xb%znu(k) = grid%znu(k)
  70       grid%xb%znw(k) = grid%znw(k)
  71       grid%xb%dn(k)  = grid%dn(k)
  72       grid%xb%dnw(k) = grid%dnw(k)
  73    end do
  74
  75    grid%xb % ptop = ptop
  76
  77    !---------------------------------------------------------------------------
  78    ! Convert WRF fitelds to xb:
  79    !---------------------------------------------------------------------------
  80
  81    !---------------------------------------------------------------
  82    ! Need this to exchange values in the halo region.
  83    ! grid%xa%u and grid%xa%v are used as temporary arrays and so
  84    ! it is easy to use the existing exchange scheme.
  85    !
  86    ! Note, this is needed as u_2 and v_2 has no guarantee
  87    ! the most east column, and the most north row are
  88    ! properly initailized for each tile.
  89    !---------------------------------------------------------------
  90
  91    ! Fill the halo region for u and v.
  92
  93 #ifdef DM_PARALLEL
  94 #include "HALO_EM_C.inc"
  95 #endif
  96
  97    !$OMP PARALLEL DO &
  98    !$OMP PRIVATE ( ij, i, j, k, cvpm, cpovcv, ppb, temp, ttb ) &
  99    !$OMP PRIVATE ( albn, qvf1, aln, theta )
 100    do ij = 1 , grid%num_tiles
 101
 102    do j=grid%j_start(ij), grid%j_end(ij)
 103       k = kte+1
 104
 105       do i=its,ite
 106          grid%p(i,j,k) = 0.0
 107          grid%xb%map_factor(i,j) = grid%msft(i,j)
 108          grid%xb%cori(i,j) = grid%f(i,j)
 109          grid%xb%tgrn(i,j) = grid%sst(i,j)
 110          if (grid%xb%tgrn(i,j) < 100.0) &
 111             grid%xb%tgrn(i,j) = grid%tmn(i,j)
 112          grid%xb%lat(i,j) = grid%xlat(i,j)
 113          grid%xb%lon(i,j) = grid%xlong(i,j)
 114          grid%xb%terr(i,j) = grid%ht(i,j)
 115          grid%xb%snow(i,j) = grid%snowc(i,j)
 116          grid%xb%lanu(i,j) = grid%lu_index(i,j)
 117          grid%xb%landmask(i,j) = grid%landmask(i,j)
 118          grid%xb%xland(i,j) = grid%xland(i,j)
 119          ! Z. Liu below are variables used by RTTOV
 120          grid%xb%tsk(i,j) = grid%tsk(i,j)
 121          grid%xb%smois(i,j) = grid%smois(i,j,1)
 122          grid%xb%tslb(i,j) = grid%tslb(i,j,1)
 123          grid%xb%xice(i,j) = grid%xice(i,j)
 124          grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
 125          grid%xb%isltyp(i,j) = grid%isltyp(i,j)
 126          grid%xb%vegfra(i,j) = grid%vegfra(i,j)
 127          grid%xb%snowh(i,j) = grid%snowh(i,j)*1000.0 ! meter to mm
 128       end do
 129
 130       cvpm =  - (1.0 - gas_constant/cp)
 131       cpovcv = cp / (cp - gas_constant)
 132
 133       ! In case of var4d, pb etc. will be recalculated in start_em with realsize=8,
 134       ! Howvwer, the originals are computed with realsize=4.
 135       if ( no_ppb ) then
 136          do k=kts,kte
 137             do i=its,ite
 138                ! The base specific volume (from real.init.code)
 139                ppb  = grid%znu(k) * grid%mub(i,j) + ptop
 140                grid%pb(i,j,k) = ppb
 141                temp = MAX ( iso_temp, base_temp + base_lapse*log(ppb/base_pres) )
 142                ttb  = temp * (base_pres/ppb)**kappa
 143                ! ttb  = (base_temp + base_lapse*log(ppb/base_pres)) * &
 144                !   (base_pres/ppb)**kappa
 145                albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm
 146
 147                qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
 148                aln  = -1.0 / (grid%mub(i,j)+grid%mu_2(i,j)) * &
 149                       (albn*grid%mu_2(i,j) + grid%rdnw(k) * &
 150                       (grid%ph_2(i,j,k+1) - grid%ph_2(i,j,k)))
 151                ! total pressure:
 152                grid%xb%p(i,j,k) = base_pres * &
 153                                  ((gas_constant*(t0+grid%t_2(i,j,k))*qvf1) / &
 154                                  (base_pres*(aln+albn)))**cpovcv
 155                ! total density
 156                grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
 157                ! pressure purtubation:
 158                grid%p(i,j,k) = grid%xb%p(i,j,k) - ppb
 159             end do
 160          end do
 161       else
 162          do k=kts,kte
 163             do i=its,ite
 164                ppb = grid%pb(i,j,k)
 165                temp = MAX ( iso_temp, base_temp + base_lapse*log(ppb/base_pres) )
 166                ttb  = temp * (base_pres/ppb)**kappa
 167                ! ttb  = (base_temp + base_lapse*log(ppb/base_pres)) * &
 168                !   (base_pres/ppb)**kappa
 169                albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm
 170
 171                qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
 172                aln  = -1.0 / (grid%mub(i,j)+grid%mu_2(i,j)) * &
 173                       (albn*grid%mu_2(i,j) + grid%rdnw(k) * &
 174                       (grid%ph_2(i,j,k+1) - grid%ph_2(i,j,k)))
 175                grid%xb%p(i,j,k) = grid%pb(i,j,k) + grid%p(i,j,k)
 176                ! total density
 177                grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
 178             end do
 179          end do
 180       endif
 181
 182       do k=kts,kte+1
 183          do i=its,ite
 184             grid%xb%hf(i,j,k) = (grid%phb(i,j,k)+grid%ph_2(i,j,k))/gravity
 185             grid%xb%w (i,j,k) = grid%w_2(i,j,k)
 186          end do
 187       end do
 188
 189       do k=kts,kte
 190          do i=its,ite
 191             grid%xb%u(i,j,k) = 0.5*(grid%u_2(i,j,k)+grid%u_2(i+1,j,k))
 192             grid%xb%v(i,j,k) = 0.5*(grid%v_2(i,j,k)+grid%v_2(i,j+1,k))
 193             grid%xb%wh(i,j,k)= 0.5*(grid%xb%w(i,j,k)+grid%xb%w(i,j,k+1))
 194             grid%xb%h(i,j,k) = 0.5*(grid%xb%hf(i,j,k)+grid%xb%hf(i,j,k+1))
 195
 196             if ( num_pseudo == 0 ) then
 197                grid%moist(i,j,k,P_QV) = max(grid%moist(i,j,k,P_QV), qlimit)
 198             end if
 199             grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
 200
 201             theta = t0 + grid%t_2(i,j,k)
 202             grid%xb%t(i,j,k) = theta*(grid%xb%p(i,j,k)/base_pres)**kappa
 203
 204             ! Convert to specific humidity from mixing ratio of water vapor:
 205             grid%xb%q(i,j,k)=grid%xb%q(i,j,k)/(1.0+grid%xb%q(i,j,k))
 206
 207             ! Background qrn needed for radar radial velocity assmiilation:
 208
 209             if (size(grid%moist,dim=4) >= 4) then
 210                grid%xb%qcw(i,j,k) = max(grid%moist(i,j,k,p_qc), 0.0)
 211                grid%xb%qrn(i,j,k) = max(grid%moist(i,j,k,p_qr), 0.0)
 212                grid%xb%qt (i,j,k) = grid%xb%q(i,j,k) + grid%xb%qcw(i,j,k) + &
 213                   grid%xb%qrn(i,j,k)
 214             end if
 215
 216             if (size(grid%moist,dim=4) >= 6) then
 217                grid%xb%qci(i,j,k) = max(grid%moist(i,j,k,p_qi), 0.0)
 218                grid%xb%qsn(i,j,k) = max(grid%moist(i,j,k,p_qs), 0.0)
 219             end if
 220
 221             if (size(grid%moist,dim=4) >= 7) then
 222                grid%xb%qgr(i,j,k) = max(grid%moist(i,j,k,p_qg), 0.0)
 223             end if
 224
 225             if ( config_flags%mp_physics == 3 ) then   ! WSM3-class scheme
 226                if ( grid%xb%t(i,j,k) <= t_kelvin ) then
 227                   grid%xb%qci(i,j,k) = grid%xb%qcw(i,j,k)
 228                   grid%xb%qcw(i,j,k) = 0.0
 229                   grid%xb%qsn(i,j,k) = grid%xb%qrn(i,j,k)
 230                   grid%xb%qrn(i,j,k) = 0.0
 231                end if
 232             end if
 233
 234          end do
 235       end do
 236
 237       do i=its,ite
 238          grid%xb%psac(i,j) = grid%mub(i,j)+grid%mu_2(i,j)
 239 ! To make the Psfc consistent with WRF (YRG, 04/06/2010):
 240          grid%xb%psfc(i,j) = grid%psfc(i,j)
 241
 242       end do
 243    end do
 244
 245    end do
 246    !$OMP END PARALLEL DO
 247
 248    !---------------------------------------------------------------------------
 249    ! [3.0] Calculate vertical inner product for use in vertical transform:
 250    !---------------------------------------------------------------------------
 251
 252    if (vertical_ip == vertical_ip_sqrt_delta_p) then
 253       ! Vertical inner product is sqrt(Delta p):
 254       do k=kts,kte
 255          grid%xb % vertical_inner_product(its:ite,jts:jte,k) = &
 256             sqrt(grid%xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k))
 257       end do
 258    else if (vertical_ip == vertical_ip_delta_p) then
 259
 260       ! Vertical inner product is Delta p:
 261       do k=1,grid%xb%mkz
 262          grid % xb % vertical_inner_product(its:ite,jts:jte,k) = &
 263          grid % xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k)
 264       end do
 265    end if
 266
 267    !---------------------------------------------------------------------------
 268    ! Calculate saturation vapour pressure and relative humidity:
 269    !---------------------------------------------------------------------------
 270
 271    !$OMP PARALLEL DO &
 272    !$OMP PRIVATE ( ij, k, j, i )
 273    do ij = 1 , grid%num_tiles
 274       do k=kts,kte
 275          do j=grid%j_start(ij),grid%j_end(ij)
 276             do i=its,ite
 277                call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
 278                   grid%xb % q(i,j,k), grid%xb %es(i,j,k), grid%xb %qs(i,j,k), &
 279                   grid%xb %rh(i,j,k))
 280             end do
 281          end do
 282       end do
 283    end do
 284    !$OMP END PARALLEL DO
 285
 286    ! Fill the halo region for xb
 287
 288 #ifdef DM_PARALLEL
 289 #include "HALO_XB.inc"
 290 #endif
 291
 292    ! Calculate time step from one dimensional cloud model parameterization
 293
 294    if (dt_cloud_model) then
 295       do j = jts, jte
 296          do i = its, ite
 297             call da_cloud_model (grid%xb%t(I,J,:),  grid%xb%p(I,J,:), &
 298                grid%xb%q(I,J,:), grid%xb%qcw(I,J,:), grid%xb%qrn(I,J,:), &
 299                grid%xb%h(I,J,:), grid%xb%hf(I,J,:), ddt, kts, kte)
 300
 301             do k = kts, kte
 302                grid%xb%delt(i,j,k) = DDT(k)
 303             end do
 304          end do
 305       end do
 306    end if
 307
 308    if (trace_use) call da_trace_exit("da_transfer_wrftoxb_lite")
 309
 310 end subroutine da_transfer_wrftoxb_lite
 311