var/da/da_transfer_model/da_transfer_xatowrf.inc

   1 subroutine da_transfer_xatowrf(grid, config_flags)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Convert analysis increments into WRF increments
   5    !    Updated for Analysis on Arakawa-C grid
   6    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   7    !
   8    !  The following WRF fields are modified:
   9    !    grid%u_2
  10    !    grid%v_2
  11    !    grid%w_2
  12    !    grid%mu_2
  13    !    grid%ph_2
  14    !    grid%t_2
  15    !    grid%moist
  16    !    grid%p
  17    !    grid%psfc
  18    !    grid%t2, grid%q2, grid%u10, grid%v10, grid%th2
  19    !
  20    !---------------------------------------------------------------------------
  21
  22 #if (WRF_CHEM == 1)
  23    USE module_state_description, ONLY : num_chem
  24 #endif
  25
  26    implicit none
  27
  28    type(domain), intent(inout)            :: grid
  29    type(grid_config_rec_type), intent(in) :: config_flags
  30
  31    integer :: i, j, k
  32 #if (WRF_CHEM == 1)
  33    integer :: ic
  34 #endif
  35
  36    real    :: sdmd, s1md
  37
  38    ! arrays to hold wrf increments on the c-grid
  39
  40 #ifdef A2C
  41    real, dimension(ims:ime,jms:jme, kms:kme) :: q_cgrid, ph_cgrid
  42 #else
  43    real, dimension(ims:ime,jms:jme, kms:kme) :: &
  44       u_cgrid, v_cgrid, q_cgrid, ph_cgrid
  45 #endif
  46
  47    real, dimension(ims:ime,jms:jme) :: mu_cgrid
  48
  49    real :: t_full, p_full, rho_dry, q_full, ph_full, ph_xb_hd, &
  50            qvf1, qvf2, qvf1_b, qvf2_b
  51
  52    real :: uu, vv, ps1, ts1, qv1, height
  53    real :: zs1, zs2, pf2, theta, thetam, mu_full, pfu, pfd, phm, cvpm, p2m
  54    real, dimension(kms:kme) :: ald, ph
  55    logical :: has_lsm_info
  56
  57    if (trace_use) call da_trace_entry("da_transfer_xatowrf")
  58
  59    has_lsm_info = .false.
  60    if ( grid%sf_surface_physics == 2 ) then
  61        if ( sum(grid%hfx*grid%hfx)   > 0.0 .and. &
  62             sum(grid%qfx*grid%qfx)   > 0.0 ) then
  63           has_lsm_info = .true.
  64        end if
  65    end if
  66
  67    ! To keep the background PH perturbation:
  68
  69    do j=jts,jte
  70       do i=its,ite
  71          do k=kts, kte+1
  72             ph_cgrid(i,j,k) = grid%ph_2(i,j,k)
  73          end do
  74       end do
  75    end do
  76
  77    !---------------------------------------------------------------------------
  78    ! [1.0] Get the mixing ratio of moisture first, as it its easy.
  79    !---------------------------------------------------------------------------
  80
  81    do k=kts,kte
  82       do j=jts,jte
  83          do i=its,ite
  84             if ((grid%xb%q(i,j,k)+grid%xa%q(i,j,k)) < 0.0) then
  85                q_cgrid(i,j,k) =-grid%xb%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
  86             else
  87                q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
  88             end if
  89          end do
  90       end do
  91    end do
  92
  93    !---------------------------------------------------------------------------
  94    ! [2.0] compute increments of dry-column air mass per unit area
  95    !---------------------------------------------------------------------------
  96
  97    do j=jts,jte
  98       do i=its,ite
  99          sdmd=0.0
 100          s1md=0.0
 101          do k=kts,kte
 102             sdmd=sdmd+q_cgrid(i,j,k)*grid%dnw(k)
 103             s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%dnw(k)
 104          end do
 105
 106          mu_cgrid(i,j)=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md
 107       end do
 108    end do
 109
 110    !---------------------------------------------------------------------------
 111    ! [3.0] compute pressure increments
 112    !---------------------------------------------------------------------------
 113
 114    if ( .not. var4d ) then ! for 4dvar, it is usefulless
 115
 116    ! Tangent linear code for grid%xa%p (based on WRF "real.init.code")
 117    ! developed by Y.-R. Guo 05/13/2004:
 118
 119    do j=jts,jte
 120       do i=its,ite
 121
 122          k = kte
 123          qvf1   = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k))
 124          qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV))
 125          qvf2   = - qvf1 / ((1.0+qvf1_b)*(1.0+qvf1_b))
 126          qvf2_b = 1.0/(1.0+qvf1_b)
 127          qvf1   = qvf1*qvf2_b + qvf1_b*qvf2
 128          qvf1_b = qvf1_b*qvf2_b
 129          grid%xa%p(i,j,k) = (-0.5/grid%rdnw(k)) * &
 130                     ((c1f(k)*mu_cgrid(i,j)+qvf1*c1f(k)*grid%mub(i,j)) / qvf2_b &
 131                      -(c1f(k)*grid%mu_2(i,j)+qvf1_b*c1f(k)*grid%mub(i,j))*qvf2/(qvf2_b*qvf2_b))
 132
 133          do k = kte-1,1,-1
 134             qvf1   = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k+1))
 135             qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV))
 136             qvf2   = - qvf1 / ((1.0+qvf1_b)*(1.0+qvf1_b))
 137             qvf2_b = 1.0/(1.0+qvf1_b)
 138             qvf1   = qvf1*qvf2_b + qvf1_b*qvf2
 139             qvf1_b = qvf1_b*qvf2_b
 140             grid%xa%p(i,j,k) = grid%xa%p(i,j,k+1)  &
 141                        - (1.0/grid%rdn(k+1)) * &
 142                        ((c1f(k)*mu_cgrid(i,j)+qvf1*c1f(k)*grid%mub(i,j)) / qvf2_b &
 143                         -(c1f(k)*grid%mu_2(i,j)+qvf1_b*c1f(k)*grid%mub(i,j))*qvf2/(qvf2_b*qvf2_b))
 144          end do
 145
 146       end do
 147    end do
 148
 149    else
 150
 151       do k=kts, kte
 152         do j=jts,jte
 153            do i=its,ite
 154               grid%xa%p(i,j,k) = 0.
 155            end do
 156         end do
 157       end do
 158
 159    endif
 160
 161    ! update perturbation pressure
 162
 163    do k=kts, kte
 164      do j=jts,jte
 165         do i=its,ite
 166            grid%p(i,j,k) = grid%p(i,j,k) + grid%xa%p(i,j,k)
 167         end do
 168      end do
 169    end do
 170
 171    ! Adjust grid%xa%q to make grid%xb%q + grid%xa%q > 0.0
 172
 173    if (check_rh == check_rh_tpw) then
 174       ! Shu-Hua~s TPW conservation:
 175       call da_check_rh(grid)
 176    else if (check_rh == check_rh_simple) then
 177       ! Simple resetting to max/min values:
 178       call da_check_rh_simple(grid)
 179    end if
 180
 181    do k=kts,kte
 182       do j=jts,jte
 183          do i=its,ite
 184             q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
 185          end do
 186       end do
 187    end do
 188
 189    !---------------------------------------------------------------------------
 190    ! [4.0] Convert temperature increments into theta increments
 191    !       Evaluate also the increments of (1/rho) and geopotential
 192    !---------------------------------------------------------------------------
 193
 194    if (print_detail_xa) then
 195       write(unit=stdout, fmt='(a, e24.12)') &
 196          'sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte)))=', &
 197          sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte))), &
 198          'sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte)))=', &
 199          sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte))), &
 200          'sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte)))=', &
 201          sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte))), &
 202          'sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte)))=', &
 203          sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte))), &
 204          'sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte)))=', &
 205          sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte)))
 206
 207        write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 208           'maxval(grid%xa%u(its:ite,jts:jte,kts:kte))=', &
 209           maxval(grid%xa%u(its:ite,jts:jte,kts:kte)), &
 210           'minval(grid%xa%u(its:ite,jts:jte,kts:kte))=', &
 211           minval(grid%xa%u(its:ite,jts:jte,kts:kte)), &
 212           'maxval(grid%xa%v(its:ite,jts:jte,kts:kte))=', &
 213           maxval(grid%xa%v(its:ite,jts:jte,kts:kte)), &
 214           'minval(grid%xa%v(its:ite,jts:jte,kts:kte))=', &
 215           minval(grid%xa%v(its:ite,jts:jte,kts:kte)), &
 216           'maxval(grid%xa%t(its:ite,jts:jte,kts:kte))=', &
 217           maxval(grid%xa%t(its:ite,jts:jte,kts:kte)), &
 218           'minval(grid%xa%t(its:ite,jts:jte,kts:kte))=', &
 219           minval(grid%xa%t(its:ite,jts:jte,kts:kte)), &
 220           'maxval(grid%xa%q(its:ite,jts:jte,kts:kte))=', &
 221           maxval(grid%xa%q(its:ite,jts:jte,kts:kte)), &
 222           'minval(grid%xa%q(its:ite,jts:jte,kts:kte))=', &
 223           minval(grid%xa%q(its:ite,jts:jte,kts:kte)), &
 224           'maxval(grid%xa%p(its:ite,jts:jte,kts:kte))=', &
 225           maxval(grid%xa%p(its:ite,jts:jte,kts:kte)), &
 226           'minval(grid%xa%p(its:ite,jts:jte,kts:kte))=', &
 227           minval(grid%xa%p(its:ite,jts:jte,kts:kte)), &
 228           'maxval(grid%xa%psfc(its:ite,jts:jte))   =', &
 229           maxval(grid%xa%psfc(its:ite,jts:jte)), &
 230           'minval(grid%xa%psfc(its:ite,jts:jte))   =', &
 231           minval(grid%xa%psfc(its:ite,jts:jte))
 232    end if
 233
 234    IF (grid%hypsometric_opt == 1) THEN
 235
 236    do j=jts,jte
 237       do i=its,ite
 238
 239          ph_full  = grid%ht(i,j) * gravity
 240          ph_xb_hd = grid%ht(i,j) * gravity
 241          do k = kts, kte
 242             ! To obtain all of the full fitelds: t, p, q(mixing ratio), rho
 243             t_full   = grid%xa%t(i,j,k) + grid%xb%t(i,j,k)
 244             p_full   = grid%xa%p(i,j,k) + grid%xb%p(i,j,k)
 245             q_full   = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k)
 246
 247             ! Note: According to WRF, this is the dry air density used to
 248             !       compute the geopotential height:
 249             rho_dry = p_full / (gas_constant*t_full*(1.0+q_full/rd_over_rv))
 250
 251             ! To compute the theta increment with the full fields:
 252             grid%t_2(i,j,k) = t_full*(base_pres/p_full)**kappa - t0
 253
 254             ! The full field of analysis ph:
 255             ph_full  = ph_full  &
 256                        - grid%xb%dnw(k) * &
 257                        ((c1h(k)*grid%xb%psac(i,j)+c2h(k))+c1h(k)*mu_cgrid(i,j)) &
 258                        / rho_dry
 259
 260             ! background hydrostatic phi:
 261             ph_xb_hd  = ph_xb_hd  &
 262                        - grid%xb%dnw(k) * &
 263                        (c1h(k)*grid%xb%psac(i,j)+c2h(k)) / grid%xb%rho(i,j,k)
 264
 265             ! The analysis perturbation = Hydro_ph - base_ph + nonhydro_xb_ph:
 266             grid%ph_2(i,j,k+1) = ph_full - grid%phb(i,j,k+1) &
 267                             + (grid%xb%hf(i,j,k+1)*gravity - ph_xb_hd)
 268          end do
 269       end do
 270    end do
 271
 272    ELSE IF  (grid%hypsometric_opt == 2) THEN
 273
 274    ! Geopotential increment reflecting hypsometric_opt: wee 11/29/2011
 275
 276    cvpm =  - (1.0 - gas_constant/cp)
 277
 278    DO j=jts,jte
 279    DO i=its,ite
 280
 281       mu_full = grid%mub(i,j)+grid%mu_2(i,j)+mu_cgrid(i,j)
 282
 283       ! Compute geopotential (using dry inverse density and dry pressure)
 284       ph      = 0.0
 285       ph(kts) = grid%ht(i,j) * gravity
 286       DO k = kts, kte
 287          p_full = grid%xb%p(i,j,k) + grid%xa%p(i,j,k)
 288          t_full = grid%xb%t(i,j,k) + grid%xa%t(i,j,k)
 289          q_full = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k)
 290          theta  = t_full*(base_pres/p_full)**kappa
 291
 292          ! Update potential temperature
 293          grid%t_2(i,j,k) = theta - t0
 294
 295          ! Compute dry inverse density using the equation of state
 296          thetam = theta*(1.0 + q_full/rd_over_rv)
 297          ald(k) = (gas_constant/base_pres)*thetam*(p_full/base_pres)**cvpm
 298
 299       ! Dry mass is purely hydrostatic: Native approach in WRF
 300
 301           pfu = mu_full*c3f(k+1) + c4f(k+1) + grid%p_top
 302           pfd = mu_full*c3f(k)   + c4f(k)   + grid%p_top
 303           phm = mu_full*c3h(k)   + c4h(k)   + grid%p_top
 304           ph(k+1) = ph(k) + ald(k)*phm*LOG(pfd/pfu)
 305           grid%ph_2(i,j,k+1) = ph(k+1) - grid%phb(i,j,k+1)
 306        END DO
 307
 308       ! Update geopotential perturbation
 309
 310 !      grid%ph_2(i,j,:) = 0.0
 311 !      DO k= kts, kte+1
 312 !         grid%ph_2(i,j,k) = ph(k) - grid%phb(i,j,k)
 313 !      END DO
 314
 315    END DO
 316    END DO
 317
 318    ENDIF ! hypsometric_opt
 319
 320    ! To compute the geopotential height increment:
 321
 322    do k=kts, kte+1
 323      do j=jts,jte
 324         do i=its,ite
 325            ph_cgrid(i,j,k) = grid%ph_2(i,j,k) - ph_cgrid(i,j,k)
 326         end do
 327      end do
 328    end do
 329
 330    ! ========================
 331    ! Write out the increment:
 332    ! ========================
 333
 334    if (write_increments) then
 335       write(unit=stdout,fmt='(/"Write out increment for plotting......")')
 336       call da_write_increments (grid, q_cgrid, mu_cgrid, ph_cgrid)
 337    end if
 338
 339 #ifdef A2C
 340
 341   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 342        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 343      ipe = ipe + 1
 344      ide = ide + 1
 345   end if
 346
 347   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 348        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 349      jpe = jpe + 1
 350      jde = jde + 1
 351   end if
 352 #else
 353    ! CONVERT FROM A-GRID TO C-GRID
 354
 355 #endif
 356
 357 #ifdef DM_PARALLEL
 358 #include "HALO_XA_A.inc"
 359 #endif
 360
 361 #ifdef A2C
 362   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 363        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 364      ipe = ipe - 1
 365      ide = ide - 1
 366   end if
 367
 368   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 369        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 370      jpe = jpe - 1
 371      jde = jde - 1
 372   end if
 373 #else
 374    ! Fill the boundary
 375
 376    ! The southern boundary
 377    if (jts == jds) then
 378       grid%xa%v(its:ite,jts-1,kts:kte)=2.0*grid%xa%v(its:ite,jts  ,kts:kte) &
 379                             -    grid%xa%v(its:ite,jts+1,kts:kte)
 380    end if
 381
 382    ! The northern boundary
 383    if (jte == jde) then
 384       grid%xa%v(its:ite,jte+1,kts:kte)=2.0*grid%xa%v(its:ite,jte  ,kts:kte) &
 385                             -    grid%xa%v(its:ite,jte-1,kts:kte)
 386    end if
 387
 388    ! The western boundary
 389    if (its == ids) then
 390       grid%xa%u(its-1,jts:jte,kts:kte)=2.0*grid%xa%u(its  ,jts:jte,kts:kte) &
 391                             -    grid%xa%u(its+1,jts:jte,kts:kte)
 392    end if
 393
 394    ! The eastern boundary
 395    if (ite == ide) then
 396       grid%xa%u(ite+1,jts:jte,kts:kte)=2.0*grid%xa%u(ite  ,jts:jte,kts:kte) &
 397                             -    grid%xa%u(ite-1,jts:jte,kts:kte)
 398    end if
 399
 400    do k=kts,kte
 401       do j=jts,jte+1
 402          do i=its,ite+1
 403             u_cgrid(i,j,k)=0.5*(grid%xa%u(i-1,j  ,k)+grid%xa%u(i,j,k))
 404             v_cgrid(i,j,k)=0.5*(grid%xa%v(i  ,j-1,k)+grid%xa%v(i,j,k))
 405          end do
 406       end do
 407    end do
 408
 409    !------------------------------------------------------------------------
 410    ! For later plot and comparation Purpose only, zero out the unused var.
 411    !------------------------------------------------------------------------
 412
 413    ! The northern boundary
 414    if (jte == jde) then
 415       u_cgrid(its:ite+1,jte+1,kts:kte)=0.0
 416    end if
 417
 418    ! The eastern boundary
 419    if (ite == ide) then
 420       v_cgrid(ite+1,jts:jte+1,kts:kte)=0.0
 421    end if
 422 #endif
 423    !---------------------------------------------------------------------------
 424    ! [5.0] add increment to the original guess and update "grid"
 425    !---------------------------------------------------------------------------
 426
 427    do j=jts,jte
 428       do i=its,ite
 429          grid%mu_2(i,j) = grid%mu_2(i,j) + mu_cgrid(i,j)
 430          grid%w_2(i,j,kte+1)=  grid%w_2(i,j,kte+1) + grid%xa%w(i,j,kte+1)
 431          grid%psfc(i,j) = grid%psfc(i,j) + grid%xa%psfc(i,j)
 432       end do
 433
 434       do k=kts,kte
 435          do i=its,ite
 436 #ifndef A2C
 437             grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k)
 438             grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k)
 439 #endif
 440             grid%w_2(i,j,k) = grid%w_2(i,j,k) + grid%xa%w(i,j,k)
 441
 442             ! (xb%q+xa%q in specific humidity) >= 0.0
 443             ! does not guarantee that (qv+q_cgrid in mixing ratio) >= 0.0
 444             ! for example, when xa%q is negative, q_cgrid is a lager negative value.
 445             ! impose a minimum value to prevent negative final qv
 446             if ( num_pseudo == 0 ) then
 447                grid%moist(i,j,k,P_QV) = max(grid%moist(i,j,k,P_QV)+q_cgrid(i,j,k), qlimit)
 448             else
 449                grid%moist(i,j,k,P_QV) = grid%moist(i,j,k,P_QV)+q_cgrid(i,j,k)
 450             end if
 451
 452             if (size(grid%moist,dim=4) >= 4 .and. cloud_cv_options >= 1) then
 453                grid%moist(i,j,k,p_qc) = max(grid%moist(i,j,k,p_qc) + grid%xa%qcw(i,j,k), 0.0)
 454                grid%moist(i,j,k,p_qr) = max(grid%moist(i,j,k,p_qr) + grid%xa%qrn(i,j,k), 0.0)
 455             end if
 456
 457             if (size(grid%moist,dim=4) >= 6 .and. cloud_cv_options >= 2) then
 458                grid%moist(i,j,k,p_qi) = max(grid%moist(i,j,k,p_qi) + grid%xa%qci(i,j,k), 0.0)
 459                grid%moist(i,j,k,p_qs) = max(grid%moist(i,j,k,p_qs) + grid%xa%qsn(i,j,k), 0.0)
 460             end if
 461
 462             if (size(grid%moist,dim=4) >= 7 .and. cloud_cv_options >= 2) then
 463                grid%moist(i,j,k,p_qg) = max(grid%moist(i,j,k,p_qg) + grid%xa%qgr(i,j,k), 0.0)
 464             end if
 465          end do
 466       end do
 467    end do
 468
 469 #ifdef A2C
 470    do j=jts,jte
 471       do k=kts,kte
 472          do i=its,ite+1
 473           grid%u_2(i,j,k) = grid%u_2(i,j,k) + grid%xa%u(i,j,k)
 474          end do
 475       end do
 476    end do
 477
 478    do j=jts,jte+1
 479       do k=kts,kte
 480          do i=its,ite
 481           grid%v_2(i,j,k) = grid%v_2(i,j,k) + grid%xa%v(i,j,k)
 482          end do
 483       end do
 484    end do
 485 #else
 486    ! The northern boundary
 487    if (jte == jde) then
 488       j=jte+1
 489       do k=kts,kte
 490          do i=its,ite
 491             grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k)
 492          end do
 493       end do
 494    end if
 495
 496    ! The eastern boundary
 497    if (ite == ide) then
 498       i=ite+1
 499       do k=kts,kte
 500          do j=jts,jte
 501             grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k)
 502          end do
 503       end do
 504    end if
 505 #endif
 506
 507 #ifdef DM_PARALLEL
 508 #include "HALO_EM_C.inc"
 509 #endif
 510 ! re-calculate T2, Q2, U10, V10, TH2 using updated fields
 511
 512    do j=jts,jte
 513       do i=its,ite
 514          uu = 0.5*(grid%u_2(i,j,kts)+grid%u_2(i+1,j,kts) )
 515          vv = 0.5*(grid%v_2(i,j,kts)+grid%v_2(i,j+1,kts) )
 516          ps1 = grid%p(i,j,kts)   + grid%pb(i,j,kts)
 517          ts1 = (t0+grid%t_2(i,j,kts))*(ps1/base_pres)**kappa
 518          qv1 = grid%moist(i,j,kts, p_qv) !qv1, input to da_sfc_wtq, is mixing ratio
 519          !hcl-07/2015 comment out below
 520          !if (grid%hypsometric_opt == 1) then
 521          !   height = 0.5*(grid%phb(i,j,kts)+grid%ph_2(i,j,kts)+ &
 522          !                 grid%phb(i,j,kts+1)+grid%ph_2(i,j,kts+1))/gravity
 523          !   height = height - grid%ht(i,j)
 524          !elseif (grid%hypsometric_opt == 2) then
 525          ! ! Height is in proportion to log pressure: wee 11/22/2011
 526          !   zs1 = (grid%phb(i,j,kts)+grid%ph_2(i,j,kts))/gravity
 527          !   zs2 = (grid%phb(i,j,kts+1)+grid%ph_2(i,j,kts+1))/gravity
 528          !
 529          !   mu_full = grid%mub(i,j)+grid%mu_2(i,j)
 530          !   pfu = mu_full*grid%znw(kts+1) + grid%p_top
 531          !   pfd = mu_full*grid%znw(kts)   + grid%p_top
 532          !   phm = mu_full*grid%znu(kts)   + grid%p_top
 533          !   height = (zs2-zs1)*LOG(pfd/phm)/LOG(pfd/pfu)
 534          !endif
 535          !hcl-07/2015 to be consistent with the height calculation done in wrf
 536          height = 0.5*(grid%phb(i,j,kts)+grid%ph_2(i,j,kts)+ &
 537                        grid%phb(i,j,kts+1)+grid%ph_2(i,j,kts+1))/gravity
 538          height = height - grid%ht(i,j)
 539          if (height <= 0.0) then
 540             message(1) = "Negative height found"
 541             write (unit=message(2),FMT='(2I6,A,F10.2,A,F10.2)') &
 542                i,j,' ht = ',height ,' terr =  ',grid%ht(i,j)
 543             call da_error(__FILE__,__LINE__, message(1:2))
 544          end if
 545          if ( update_sfcdiags ) then
 546             if ( use_wrf_sfcinfo ) then
 547                call da_sfc_wtq(grid%psfc(i,j), grid%tsk(i,j),                 &
 548                   ps1, ts1, qv1, uu, vv,                                      &
 549                   height,  grid%xb%rough(i,j),grid%xb%xland(i,j), grid%xb%ds, &
 550                   grid%u10(i,j), grid%v10(i,j), grid%t2(i,j),                 &
 551                   grid%q2(i,j), grid%xb%regime(i,j),                          &
 552                   has_lsm_info, regime_wrf=grid%regime(i,j),                  &
 553                   qsfc_wrf=grid%qsfc(i,j), znt_wrf=grid%znt(i,j),             &
 554                   ust_wrf=grid%ust(i,j), mol_wrf=grid%mol(i,j),               &
 555                   hfx=grid%hfx(i,j), qfx=grid%qfx(i,j), pblh=grid%pblh(i,j) )
 556             else
 557                call da_sfc_wtq(grid%psfc(i,j), grid%tsk(i,j),                 &
 558                   ps1, ts1, qv1, uu, vv,                                      &
 559                   height,  grid%xb%rough(i,j),grid%xb%xland(i,j), grid%xb%ds, &
 560                   grid%u10(i,j), grid%v10(i,j), grid%t2(i,j),                 &
 561                   grid%q2(i,j), grid%xb%regime(i,j))
 562             end if
 563
 564             ! 2-m pressure: Ground level and first half level are used
 565             !hcl p2m = grid%psfc(i,j)*EXP(-2.0/height*LOG(grid%psfc(i,j)/ps1))
 566             !hcl grid%th2(i,j) = grid%t2(i,j)*(base_pres/p2m)**kappa
 567
 568             !hcl-07/2015 to be consistent with the th2 calculation done in wrf
 569             grid%th2(i,j) = grid%t2(i,j)*(base_pres/grid%psfc(i,j))**kappa
 570          end if ! if update_sfcdiags
 571       end do
 572    end do
 573
 574    if (print_detail_xa) then
 575       write(unit=stdout, fmt=*) 'simple variables:'
 576
 577       if (ite == ide) then
 578          write (unit=stdout,fmt=*)  ' '
 579
 580          do k=kts+5,kte,10
 581             do j=jts,jte,10
 582                write(unit=stdout, fmt=*) &
 583                     '  grid%u_2(', ide+1, ',', j, ',', k, ')=', &
 584                        grid%u_2(ide+1,j,k)
 585             end do
 586             write(unit=stdout, fmt=*) ' '
 587          end do
 588       end if
 589
 590       if (jte == jde) then
 591          write(unit=stdout, fmt=*) ' '
 592
 593          do k=kts+5,kte,10
 594             do i=its,ite,10
 595                write(unit=stdout, fmt=*) &
 596                     '  grid%v_2(', i, ',', jde+1, ',', k, ')=', &
 597                        grid%v_2(i, jde+1,k)
 598             end do
 599             write(unit=stdout, fmt=*) ' '
 600          end do
 601       end if
 602
 603 #ifdef A2C
 604       write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 605          'maxval(mu_cgrid(its:ite,jts:jte))       =', &
 606          maxval(mu_cgrid(its:ite,jts:jte)), &
 607          'minval(mu_cgrid(its:ite,jts:jte))       =', &
 608          minval(mu_cgrid(its:ite,jts:jte)), &
 609          'maxval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 610          maxval(q_cgrid(its:ite,jts:jte,kts:kte)), &
 611          'minval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 612          minval(q_cgrid(its:ite,jts:jte,kts:kte))
 613 #else
 614       write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 615          'maxval(mu_cgrid(its:ite,jts:jte))       =', &
 616          maxval(mu_cgrid(its:ite,jts:jte)), &
 617          'minval(mu_cgrid(its:ite,jts:jte))       =', &
 618          minval(mu_cgrid(its:ite,jts:jte)), &
 619         'maxval(u_cgrid(its:ite,jts:jte,kts:kte))  =', &
 620          maxval(u_cgrid(its:ite,jts:jte,kts:kte)), &
 621          'minval(u_cgrid(its:ite,jts:jte,kts:kte))  =', &
 622          minval(u_cgrid(its:ite,jts:jte,kts:kte)), &
 623          'maxval(v_cgrid(its:ite,jts:jte,kts:kte))  =', &
 624          maxval(v_cgrid(its:ite,jts:jte,kts:kte)), &
 625          'minval(v_cgrid(its:ite,jts:jte,kts:kte))  =', &
 626          minval(v_cgrid(its:ite,jts:jte,kts:kte)), &
 627          'maxval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 628          maxval(q_cgrid(its:ite,jts:jte,kts:kte)), &
 629          'minval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 630          minval(q_cgrid(its:ite,jts:jte,kts:kte))
 631 #endif
 632
 633       do k=kts,kte
 634          write(unit=stdout, fmt='(a, i3)') 'k=', k
 635
 636 #ifdef A2C
 637          write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 638             'maxval(q_cgrid(its:ite,jts:jte,k))  =', maxval(q_cgrid(its:ite,jts:jte,k)), &
 639             'minval(q_cgrid(its:ite,jts:jte,k))  =', minval(q_cgrid(its:ite,jts:jte,k))
 640 #else
 641          write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 642             'maxval(u_cgrid(its:ite,jts:jte,k))  =', maxval(u_cgrid(its:ite,jts:jte,k)), &
 643             'minval(u_cgrid(its:ite,jts:jte,k))  =', minval(u_cgrid(its:ite,jts:jte,k)), &
 644             'maxval(v_cgrid(its:ite,jts:jte,k))  =', maxval(v_cgrid(its:ite,jts:jte,k)), &
 645             'minval(v_cgrid(its:ite,jts:jte,k))  =', minval(v_cgrid(its:ite,jts:jte,k)), &
 646             'maxval(q_cgrid(its:ite,jts:jte,k))  =', maxval(q_cgrid(its:ite,jts:jte,k)), &
 647             'minval(q_cgrid(its:ite,jts:jte,k))  =', minval(q_cgrid(its:ite,jts:jte,k))
 648 #endif
 649       end do
 650    end if
 651 #ifdef VAR4D
 652    IF ( var4d ) THEN
 653
 654       ! We need to transfer LBC perturbation from model_grid to grid for output
 655
 656       grid%u_bxs = grid%u_bxs + model_grid%g_u_bxs
 657       grid%u_bxe = grid%u_bxe + model_grid%g_u_bxe
 658       grid%u_bys = grid%u_bys + model_grid%g_u_bys
 659       grid%u_bye = grid%u_bye + model_grid%g_u_bye
 660       grid%v_bxs = grid%v_bxs + model_grid%g_v_bxs
 661       grid%v_bxe = grid%v_bxe + model_grid%g_v_bxe
 662       grid%v_bys = grid%v_bys + model_grid%g_v_bys
 663       grid%v_bye = grid%v_bye + model_grid%g_v_bye
 664 !     grid%w_bxs = grid%w_bxs + model_grid%g_w_bxs
 665 !     grid%w_bxe = grid%w_bxe + model_grid%g_w_bxe
 666 !     grid%w_bys = grid%w_bys + model_grid%g_w_bys
 667 !     grid%w_bye = grid%w_bye + model_grid%g_w_bye
 668       grid%t_bxs = grid%t_bxs + model_grid%g_t_bxs
 669       grid%t_bxe = grid%t_bxe + model_grid%g_t_bxe
 670       grid%t_bys = grid%t_bys + model_grid%g_t_bys
 671       grid%t_bye = grid%t_bye + model_grid%g_t_bye
 672       grid%mu_bxs = grid%mu_bxs + model_grid%g_mu_bxs
 673       grid%mu_bxe = grid%mu_bxe + model_grid%g_mu_bxe
 674       grid%mu_bys = grid%mu_bys + model_grid%g_mu_bys
 675       grid%mu_bye = grid%mu_bye + model_grid%g_mu_bye
 676       grid%ph_bxs = grid%ph_bxs + model_grid%g_ph_bxs
 677       grid%ph_bxe = grid%ph_bxe + model_grid%g_ph_bxe
 678       grid%ph_bys = grid%ph_bys + model_grid%g_ph_bys
 679       grid%ph_bye = grid%ph_bye + model_grid%g_ph_bye
 680       grid%moist_bxs = grid%moist_bxs + model_grid%g_moist_bxs
 681       grid%moist_bxe = grid%moist_bxe + model_grid%g_moist_bxe
 682       grid%moist_bys = grid%moist_bys + model_grid%g_moist_bys
 683       grid%moist_bye = grid%moist_bye + model_grid%g_moist_bye
 684
 685       grid%u_btxs = grid%u_btxs + model_grid%g_u_btxs
 686       grid%u_btxe = grid%u_btxe + model_grid%g_u_btxe
 687       grid%u_btys = grid%u_btys + model_grid%g_u_btys
 688       grid%u_btye = grid%u_btye + model_grid%g_u_btye
 689       grid%v_btxs = grid%v_btxs + model_grid%g_v_btxs
 690       grid%v_btxe = grid%v_btxe + model_grid%g_v_btxe
 691       grid%v_btys = grid%v_btys + model_grid%g_v_btys
 692       grid%v_btye = grid%v_btye + model_grid%g_v_btye
 693 !     grid%w_btxs = grid%w_btxs + model_grid%g_w_btxs
 694 !     grid%w_btxe = grid%w_btxe + model_grid%g_w_btxe
 695 !     grid%w_btys = grid%w_btys + model_grid%g_w_btys
 696 !     grid%w_btye = grid%w_btye + model_grid%g_w_btye
 697       grid%t_btxs = grid%t_btxs + model_grid%g_t_btxs
 698       grid%t_btxe = grid%t_btxe + model_grid%g_t_btxe
 699       grid%t_btys = grid%t_btys + model_grid%g_t_btys
 700       grid%t_btye = grid%t_btye + model_grid%g_t_btye
 701       grid%mu_btxs = grid%mu_btxs + model_grid%g_mu_btxs
 702       grid%mu_btxe = grid%mu_btxe + model_grid%g_mu_btxe
 703       grid%mu_btys = grid%mu_btys + model_grid%g_mu_btys
 704       grid%mu_btye = grid%mu_btye + model_grid%g_mu_btye
 705       grid%ph_btxs = grid%ph_btxs + model_grid%g_ph_btxs
 706       grid%ph_btxe = grid%ph_btxe + model_grid%g_ph_btxe
 707       grid%ph_btys = grid%ph_btys + model_grid%g_ph_btys
 708       grid%ph_btye = grid%ph_btye + model_grid%g_ph_btye
 709       grid%moist_btxs = grid%moist_btxs + model_grid%g_moist_btxs
 710       grid%moist_btxe = grid%moist_btxe + model_grid%g_moist_btxe
 711       grid%moist_btys = grid%moist_btys + model_grid%g_moist_btys
 712       grid%moist_btye = grid%moist_btye + model_grid%g_moist_btye
 713
 714    ENDIF
 715 #else !! not var4d !!
 716
 717 #if (WRF_CHEM == 1)
 718       do ic = PARAM_FIRST_SCALAR, num_chem
 719          do j=jts,jte
 720             do i=its,ite
 721                do k=kts, kte
 722                 if (grid%xbchem%chem_ic(i,j,k,ic) + grid%xachem%chem_ic(i,j,k,ic).lt.0.0) then
 723                   grid%chem(i,j,k,ic) = grid%xbchem%chem_ic(i,j,k,ic)
 724                 else
 725                   grid%chem(i,j,k,ic) = grid%xbchem%chem_ic(i,j,k,ic) + grid%xachem%chem_ic(i,j,k,ic)
 726                 end if
 727                end do
 728             end do
 729          end do
 730       end do
 731
 732 #endif
 733
 734 #endif
 735
 736    if (trace_use) call da_trace_exit("da_transfer_xatowrf")
 737
 738 end subroutine da_transfer_xatowrf
 739
 740