subroutine da_transfer_xatowrf(grid) !--------------------------------------------------------------------------- ! Purpose: Convert analysis increments into WRF increments ! Updated for Analysis on Arakawa-C grid ! Author: Syed RH Rizvi, MMM/ESSL/NCAR, Date: 10/22/2008 ! ! The following WRF fields are modified: ! grid%u_2 ! grid%v_2 ! grid%w_2 ! grid%mu_2 ! grid%ph_2 ! grid%t_2 ! grid%moist ! grid%p ! grid%psfc ! grid%t2, grid%q2, grid%u10, grid%v10, grid%th2 ! !--------------------------------------------------------------------------- implicit none type(domain), intent(inout) :: grid integer :: i, j, k real :: sdmd, s1md ! arrays to hold wrf increments on the c-grid #ifdef A2C real, dimension(ims:ime,jms:jme, kms:kme) :: q_cgrid, ph_cgrid #else real, dimension(ims:ime,jms:jme, kms:kme) :: & u_cgrid, v_cgrid, q_cgrid, ph_cgrid #endif real, dimension(ims:ime,jms:jme) :: mu_cgrid real :: t_full, p_full, rho_dry, q_full, ph_full, ph_xb_hd, & qvf1, qvf2, qvf1_b, qvf2_b real :: uu, vv, ps1, ps2, ts1, ts2, qv1, qv2, height real :: zs1, zs2, pf2, theta, thetam, mu_full, pfu, pfd, phm, cvpm, p2m real, dimension(kms:kme) :: ald, ph if (trace_use) call da_trace_entry("da_transfer_xatowrf") ! To keep the background PH perturbation: do j=jts,jte do i=its,ite do k=kts, kte+1 ph_cgrid(i,j,k) = grid%ph_2(i,j,k) end do end do end do !--------------------------------------------------------------------------- ! [1.0] Get the mixing ratio of moisture first, as it its easy. !--------------------------------------------------------------------------- do k=kts,kte do j=jts,jte do i=its,ite if ((grid%xb%q(i,j,k)+grid%xa%q(i,j,k)) < 0.0) then q_cgrid(i,j,k) =-grid%xb%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2 else q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2 end if end do end do end do !--------------------------------------------------------------------------- ! [2.0] compute increments of dry-column air mass per unit area !--------------------------------------------------------------------------- do j=jts,jte do i=its,ite sdmd=0.0 s1md=0.0 do k=kts,kte sdmd=sdmd+q_cgrid(i,j,k)*grid%dnw(k) s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%dnw(k) end do mu_cgrid(i,j)=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md end do end do !--------------------------------------------------------------------------- ! [3.0] compute pressure increments !--------------------------------------------------------------------------- if ( .not. var4d ) then ! for 4dvar, it is usefulless ! Tangent linear code for grid%xa%p (based on WRF "real.init.code") ! developed by Y.-R. Guo 05/13/2004: do j=jts,jte do i=its,ite k = kte qvf1 = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k)) qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV)) qvf2 = - qvf1 / ((1.0+qvf1_b)*(1.0+qvf1_b)) qvf2_b = 1.0/(1.0+qvf1_b) qvf1 = qvf1*qvf2_b + qvf1_b*qvf2 qvf1_b = qvf1_b*qvf2_b grid%xa%p(i,j,k) = (-0.5/grid%rdnw(k)) * & ((mu_cgrid(i,j)+qvf1*grid%mub(i,j)) / qvf2_b & -(grid%mu_2(i,j)+qvf1_b*grid%mub(i,j))*qvf2/(qvf2_b*qvf2_b)) do k = kte-1,1,-1 qvf1 = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k+1)) qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV)) qvf2 = - qvf1 / ((1.0+qvf1_b)*(1.0+qvf1_b)) qvf2_b = 1.0/(1.0+qvf1_b) qvf1 = qvf1*qvf2_b + qvf1_b*qvf2 qvf1_b = qvf1_b*qvf2_b grid%xa%p(i,j,k) = grid%xa%p(i,j,k+1) & - (1.0/grid%rdn(k+1)) * & ((mu_cgrid(i,j)+qvf1*grid%mub(i,j)) / qvf2_b & -(grid%mu_2(i,j)+qvf1_b*grid%mub(i,j))*qvf2/(qvf2_b*qvf2_b)) end do end do end do else do k=kts, kte do j=jts,jte do i=its,ite grid%xa%p(i,j,k) = 0. end do end do end do endif ! update perturbation pressure do k=kts, kte do j=jts,jte do i=its,ite grid%p(i,j,k) = grid%p(i,j,k) + grid%xa%p(i,j,k) end do end do end do ! Adjust grid%xa%q to make grid%xb%q + grid%xa%q > 0.0 if (check_rh == check_rh_tpw) then ! Shu-Hua~s TPW conservation: call da_check_rh(grid) else if (check_rh == check_rh_simple) then ! Simple resetting to max/min values: call da_check_rh_simple(grid) end if do k=kts,kte do j=jts,jte do i=its,ite q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2 end do end do end do !--------------------------------------------------------------------------- ! [4.0] Convert temperature increments into theta increments ! Evaluate also the increments of (1/rho) and geopotential !--------------------------------------------------------------------------- if (print_detail_xa) then write(unit=stdout, fmt='(a, e24.12)') & 'sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte)))=', & sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte))), & 'sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte)))=', & sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte))), & 'sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte)))=', & sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte))), & 'sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte)))=', & sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte))), & 'sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte)))=', & sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte))) write(unit=stdout, fmt='(2(2x, a, e20.12))') & 'maxval(grid%xa%u(its:ite,jts:jte,kts:kte))=', & maxval(grid%xa%u(its:ite,jts:jte,kts:kte)), & 'minval(grid%xa%u(its:ite,jts:jte,kts:kte))=', & minval(grid%xa%u(its:ite,jts:jte,kts:kte)), & 'maxval(grid%xa%v(its:ite,jts:jte,kts:kte))=', & maxval(grid%xa%v(its:ite,jts:jte,kts:kte)), & 'minval(grid%xa%v(its:ite,jts:jte,kts:kte))=', & minval(grid%xa%v(its:ite,jts:jte,kts:kte)), & 'maxval(grid%xa%t(its:ite,jts:jte,kts:kte))=', & maxval(grid%xa%t(its:ite,jts:jte,kts:kte)), & 'minval(grid%xa%t(its:ite,jts:jte,kts:kte))=', & minval(grid%xa%t(its:ite,jts:jte,kts:kte)), & 'maxval(grid%xa%q(its:ite,jts:jte,kts:kte))=', & maxval(grid%xa%q(its:ite,jts:jte,kts:kte)), & 'minval(grid%xa%q(its:ite,jts:jte,kts:kte))=', & minval(grid%xa%q(its:ite,jts:jte,kts:kte)), & 'maxval(grid%xa%p(its:ite,jts:jte,kts:kte))=', & maxval(grid%xa%p(its:ite,jts:jte,kts:kte)), & 'minval(grid%xa%p(its:ite,jts:jte,kts:kte))=', & minval(grid%xa%p(its:ite,jts:jte,kts:kte)), & 'maxval(grid%xa%psfc(its:ite,jts:jte)) =', & maxval(grid%xa%psfc(its:ite,jts:jte)), & 'minval(grid%xa%psfc(its:ite,jts:jte)) =', & minval(grid%xa%psfc(its:ite,jts:jte)) end if IF (grid%hypsometric_opt == 1) THEN do j=jts,jte do i=its,ite ph_full = grid%ht(i,j) * gravity ph_xb_hd = grid%ht(i,j) * gravity do k = kts, kte ! To obtain all of the full fitelds: t, p, q(mixing ratio), rho t_full = grid%xa%t(i,j,k) + grid%xb%t(i,j,k) p_full = grid%xa%p(i,j,k) + grid%xb%p(i,j,k) q_full = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k) ! Note: According to WRF, this is the dry air density used to ! compute the geopotential height: rho_dry = p_full / (gas_constant*t_full*(1.0+q_full/rd_over_rv)) ! To compute the theta increment with the full fields: grid%t_2(i,j,k) = t_full*(base_pres/p_full)**kappa - t0 ! The full fiteld of analysis ph: ph_full = ph_full & - grid%xb%dnw(k) * (grid%xb%psac(i,j)+mu_cgrid(i,j)) / rho_dry ! background hydrostatic phi: ph_xb_hd = ph_xb_hd & - grid%xb%dnw(k) * grid%xb%psac(i,j) / grid%xb%rho(i,j,k) ! The analysis perturbation = Hydro_ph - base_ph + nonhydro_xb_ph: grid%ph_2(i,j,k+1) = ph_full - grid%phb(i,j,k+1) & + (grid%xb%hf(i,j,k+1)*gravity - ph_xb_hd) end do end do end do ELSE IF (grid%hypsometric_opt == 2) THEN ! Geopotential increment reflecting hypsometric_opt: wee 11/29/2011 cvpm = - (1.0 - gas_constant/cp) DO j=jts,jte DO i=its,ite ! dry air density mu_full = grid%mub(i,j)+grid%mu_2(i,j)+mu_cgrid(i,j) ph = 0.0 ! Compute geopotential (using dry inverse density and dry pressure) ph(kts) = grid%ht(i,j) * gravity DO k = kts, kte p_full = grid%xb%p(i,j,k) + grid%xa%p(i,j,k) t_full = grid%xb%t(i,j,k) + grid%xa%t(i,j,k) q_full = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k) theta = t_full*(base_pres/p_full)**kappa ! Update potential temperature grid%t_2(i,j,k) = theta - t0 ! Compute dry inverse density using the equation of state thetam = theta*(1.0 + q_full/rd_over_rv) ald(k) = (gas_constant/base_pres)*thetam*(p_full/base_pres)**cvpm ! Dry mass is purely hydrostatic: Native approach in WRF pfu = mu_full*grid%znw(k+1) + grid%p_top pfd = mu_full*grid%znw(k) + grid%p_top phm = mu_full*grid%znu(k) + grid%p_top ph(k+1) = ph(k) + ald(k)*phm*LOG(pfd/pfu) grid%ph_2(i,j,k+1) = ph(k+1) - grid%phb(i,j,k+1) END DO ! Update geopotential perturbation ! grid%ph_2(i,j,:) = 0.0 ! DO k= kts, kte+1 ! grid%ph_2(i,j,k) = ph(k) - grid%phb(i,j,k) ! END DO END DO END DO ENDIF ! hypsometric_opt ! To compute the geopotential height increment: do k=kts, kte+1 do j=jts,jte do i=its,ite ph_cgrid(i,j,k) = grid%ph_2(i,j,k) - ph_cgrid(i,j,k) end do end do end do ! ======================== ! Write out the increment: ! ======================== if (write_increments) then write(unit=stdout,fmt='(/"Write out increment for plotting......")') call da_write_increments (grid, q_cgrid, mu_cgrid, ph_cgrid) end if #ifdef A2C if ((fg_format==fg_format_wrf_arw_regional .or. & fg_format==fg_format_wrf_arw_global ) .and. ide == ipe ) then ipe = ipe + 1 ide = ide + 1 end if if ((fg_format==fg_format_wrf_arw_regional .or. & fg_format==fg_format_wrf_arw_global ) .and. jde == jpe ) then jpe = jpe + 1 jde = jde + 1 end if #else ! CONVERT FROM A-GRID TO C-GRID #endif #ifdef DM_PARALLEL #include "HALO_XA_A.inc" #endif #ifdef A2C if ((fg_format==fg_format_wrf_arw_regional .or. & fg_format==fg_format_wrf_arw_global ) .and. ide == ipe ) then ipe = ipe - 1 ide = ide - 1 end if if ((fg_format==fg_format_wrf_arw_regional .or. & fg_format==fg_format_wrf_arw_global ) .and. jde == jpe ) then jpe = jpe - 1 jde = jde - 1 end if #else ! Fill the boundary ! The southern boundary if (jts == jds) then grid%xa%v(its:ite,jts-1,kts:kte)=2.0*grid%xa%v(its:ite,jts ,kts:kte) & - grid%xa%v(its:ite,jts+1,kts:kte) end if ! The northern boundary if (jte == jde) then grid%xa%v(its:ite,jte+1,kts:kte)=2.0*grid%xa%v(its:ite,jte ,kts:kte) & - grid%xa%v(its:ite,jte-1,kts:kte) end if ! The western boundary if (its == ids) then grid%xa%u(its-1,jts:jte,kts:kte)=2.0*grid%xa%u(its ,jts:jte,kts:kte) & - grid%xa%u(its+1,jts:jte,kts:kte) end if ! The eastern boundary if (ite == ide) then grid%xa%u(ite+1,jts:jte,kts:kte)=2.0*grid%xa%u(ite ,jts:jte,kts:kte) & - grid%xa%u(ite-1,jts:jte,kts:kte) end if do k=kts,kte do j=jts,jte+1 do i=its,ite+1 u_cgrid(i,j,k)=0.5*(grid%xa%u(i-1,j ,k)+grid%xa%u(i,j,k)) v_cgrid(i,j,k)=0.5*(grid%xa%v(i ,j-1,k)+grid%xa%v(i,j,k)) end do end do end do !------------------------------------------------------------------------ ! For later plot and comparation Purpose only, zero out the unused var. !------------------------------------------------------------------------ ! The northern boundary if (jte == jde) then u_cgrid(its:ite+1,jte+1,kts:kte)=0.0 end if ! The eastern boundary if (ite == ide) then v_cgrid(ite+1,jts:jte+1,kts:kte)=0.0 end if #endif !--------------------------------------------------------------------------- ! [5.0] add increment to the original guess and update "grid" !--------------------------------------------------------------------------- do j=jts,jte do i=its,ite grid%mu_2(i,j) = grid%mu_2(i,j) + mu_cgrid(i,j) grid%w_2(i,j,kte+1)= grid%w_2(i,j,kte+1) + grid%xa%w(i,j,kte+1) grid%psfc(i,j) = grid%psfc(i,j) + grid%xa%psfc(i,j) end do do k=kts,kte do i=its,ite #ifndef A2C grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k) grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k) #endif grid%w_2(i,j,k) = grid%w_2(i,j,k) + grid%xa%w(i,j,k) ! (xb%q+xa%q in specific humidity) >= 0.0 ! does not guarantee that (qv+q_cgrid in mixing ratio) >= 0.0 ! for example, when xa%q is negative, q_cgrid is a lager negative value. ! impose a minimum value to prevent negative final qv if ( num_pseudo == 0 ) then grid%moist(i,j,k,P_QV) = max(grid%moist(i,j,k,P_QV)+q_cgrid(i,j,k), qlimit) else grid%moist(i,j,k,P_QV) = grid%moist(i,j,k,P_QV)+q_cgrid(i,j,k) end if if (size(grid%moist,dim=4) >= 4) then grid%moist(i,j,k,p_qc) = max(grid%moist(i,j,k,p_qc) + grid%xa%qcw(i,j,k), 0.0) grid%moist(i,j,k,p_qr) = max(grid%moist(i,j,k,p_qr) + grid%xa%qrn(i,j,k), 0.0) end if if (size(grid%moist,dim=4) >= 6) then grid%moist(i,j,k,p_qi) = max(grid%moist(i,j,k,p_qi) + grid%xa%qci(i,j,k), 0.0) grid%moist(i,j,k,p_qs) = max(grid%moist(i,j,k,p_qs) + grid%xa%qsn(i,j,k), 0.0) end if if (size(grid%moist,dim=4) >= 7) then grid%moist(i,j,k,p_qg) = max(grid%moist(i,j,k,p_qg) + grid%xa%qgr(i,j,k), 0.0) end if end do end do end do #ifdef A2C do j=jts,jte do k=kts,kte do i=its,ite+1 grid%u_2(i,j,k) = grid%u_2(i,j,k) + grid%xa%u(i,j,k) end do end do end do do j=jts,jte+1 do k=kts,kte do i=its,ite grid%v_2(i,j,k) = grid%v_2(i,j,k) + grid%xa%v(i,j,k) end do end do end do #else ! The northern boundary if (jte == jde) then j=jte+1 do k=kts,kte do i=its,ite grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k) end do end do end if ! The eastern boundary if (ite == ide) then i=ite+1 do k=kts,kte do j=jts,jte grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k) end do end do end if #endif #ifdef DM_PARALLEL #include "HALO_EM_C.inc" #endif ! re-calculate T2, Q2, U10, V10, TH2 using updated fields do j=jts,jte do i=its,ite uu = 0.5*(grid%u_2(i,j,kts)+grid%u_2(i+1,j,kts) ) vv = 0.5*(grid%v_2(i,j,kts)+grid%v_2(i,j+1,kts) ) ps1 = grid%p(i,j,kts) + grid%pb(i,j,kts) ps2 = grid%p(i,j,kts+1) + grid%pb(i,j,kts+1) ts1 = (t0+grid%t_2(i,j,kts))*(ps1/base_pres)**kappa ts2 = (t0+grid%t_2(i,j,kts+1))*(ps2/base_pres)**kappa qv1 = grid%moist(i,j,kts, p_qv)/(1.0+grid%moist(i,j,kts,p_qv)) qv2 = grid%moist(i,j,kts+1, p_qv)/(1.0+grid%moist(i,j,kts+1,p_qv)) if (grid%hypsometric_opt == 1) then height = 0.5*(grid%phb(i,j,kts)+grid%ph_2(i,j,kts)+ & grid%phb(i,j,kts+1)+grid%ph_2(i,j,kts+1))/gravity height = height - grid%ht(i,j) elseif (grid%hypsometric_opt == 2) then ! Height is in proportion to log pressure: wee 11/22/2011 zs1 = (grid%phb(i,j,kts)+grid%ph_2(i,j,kts))/gravity zs2 = (grid%phb(i,j,kts+1)+grid%ph_2(i,j,kts+1))/gravity mu_full = grid%mub(i,j)+grid%mu_2(i,j) pfu = mu_full*grid%znw(kts+1) + grid%p_top pfd = mu_full*grid%znw(kts) + grid%p_top phm = mu_full*grid%znu(kts) + grid%p_top height = (zs2-zs1)*LOG(pfd/phm)/LOG(pfd/pfu) endif if (height <= 0.0) then message(1) = "Negative height found" write (unit=message(2),FMT='(2I6,A,F10.2,A,F10.2)') & i,j,' ht = ',height ,' terr = ',grid%ht(i,j) call da_error(__FILE__,__LINE__, message(1:2)) end if call da_sfc_wtq(grid%psfc(i,j), grid%tsk(i,j), & ps1, ts1, qv1, uu, vv, & ps2, ts2, qv2, & height, grid%xb%rough(i,j),grid%xb%xland(i,j), grid%xb%ds, & grid%u10(i,j), grid%v10(i,j), grid%t2(i,j), & grid%q2(i,j), grid%xb%regime(i,j)) ! This is obviously incorrect: wee 11/22/2011 ! grid%th2(i,j) = grid%t2(i,j)*(base_pres/ps1)**kappa ! 2-m pressure: Ground level and first half level are used p2m = grid%psfc(i,j)*EXP(-2.0/height*LOG(grid%psfc(i,j)/ps1)) grid%th2(i,j) = grid%t2(i,j)*(base_pres/p2m)**kappa end do end do if (print_detail_xa) then write(unit=stdout, fmt=*) 'simple variables:' if (ite == ide) then write (unit=stdout,fmt=*) ' ' do k=kts+5,kte,10 do j=jts,jte,10 write(unit=stdout, fmt=*) & ' grid%u_2(', ide+1, ',', j, ',', k, ')=', & grid%u_2(ide+1,j,k) end do write(unit=stdout, fmt=*) ' ' end do end if if (jte == jde) then write(unit=stdout, fmt=*) ' ' do k=kts+5,kte,10 do i=its,ite,10 write(unit=stdout, fmt=*) & ' grid%v_2(', i, ',', jde+1, ',', k, ')=', & grid%v_2(i, jde+1,k) end do write(unit=stdout, fmt=*) ' ' end do end if #ifdef A2C write(unit=stdout, fmt='(2(2x, a, e20.12))') & 'maxval(mu_cgrid(its:ite,jts:jte)) =', & maxval(mu_cgrid(its:ite,jts:jte)), & 'minval(mu_cgrid(its:ite,jts:jte)) =', & minval(mu_cgrid(its:ite,jts:jte)), & 'maxval(q_cgrid(its:ite,jts:jte,kts:kte)) =', & maxval(q_cgrid(its:ite,jts:jte,kts:kte)), & 'minval(q_cgrid(its:ite,jts:jte,kts:kte)) =', & minval(q_cgrid(its:ite,jts:jte,kts:kte)) #else write(unit=stdout, fmt='(2(2x, a, e20.12))') & 'maxval(mu_cgrid(its:ite,jts:jte)) =', & maxval(mu_cgrid(its:ite,jts:jte)), & 'minval(mu_cgrid(its:ite,jts:jte)) =', & minval(mu_cgrid(its:ite,jts:jte)), & 'maxval(u_cgrid(its:ite,jts:jte,kts:kte)) =', & maxval(u_cgrid(its:ite,jts:jte,kts:kte)), & 'minval(u_cgrid(its:ite,jts:jte,kts:kte)) =', & minval(u_cgrid(its:ite,jts:jte,kts:kte)), & 'maxval(v_cgrid(its:ite,jts:jte,kts:kte)) =', & maxval(v_cgrid(its:ite,jts:jte,kts:kte)), & 'minval(v_cgrid(its:ite,jts:jte,kts:kte)) =', & minval(v_cgrid(its:ite,jts:jte,kts:kte)), & 'maxval(q_cgrid(its:ite,jts:jte,kts:kte)) =', & maxval(q_cgrid(its:ite,jts:jte,kts:kte)), & 'minval(q_cgrid(its:ite,jts:jte,kts:kte)) =', & minval(q_cgrid(its:ite,jts:jte,kts:kte)) #endif do k=kts,kte write(unit=stdout, fmt='(a, i3)') 'k=', k #ifdef A2C write(unit=stdout, fmt='(2(2x, a, e20.12))') & 'maxval(q_cgrid(its:ite,jts:jte,k)) =', maxval(q_cgrid(its:ite,jts:jte,k)), & 'minval(q_cgrid(its:ite,jts:jte,k)) =', minval(q_cgrid(its:ite,jts:jte,k)) #else write(unit=stdout, fmt='(2(2x, a, e20.12))') & 'maxval(u_cgrid(its:ite,jts:jte,k)) =', maxval(u_cgrid(its:ite,jts:jte,k)), & 'minval(u_cgrid(its:ite,jts:jte,k)) =', minval(u_cgrid(its:ite,jts:jte,k)), & 'maxval(v_cgrid(its:ite,jts:jte,k)) =', maxval(v_cgrid(its:ite,jts:jte,k)), & 'minval(v_cgrid(its:ite,jts:jte,k)) =', minval(v_cgrid(its:ite,jts:jte,k)), & 'maxval(q_cgrid(its:ite,jts:jte,k)) =', maxval(q_cgrid(its:ite,jts:jte,k)), & 'minval(q_cgrid(its:ite,jts:jte,k)) =', minval(q_cgrid(its:ite,jts:jte,k)) #endif end do end if #ifdef VAR4D IF ( var4d ) THEN ! We need to transfer LBC perturbation from model_grid to grid for output grid%u_bxs = grid%u_bxs + model_grid%g_u_bxs grid%u_bxe = grid%u_bxe + model_grid%g_u_bxe grid%u_bys = grid%u_bys + model_grid%g_u_bys grid%u_bye = grid%u_bye + model_grid%g_u_bye grid%v_bxs = grid%v_bxs + model_grid%g_v_bxs grid%v_bxe = grid%v_bxe + model_grid%g_v_bxe grid%v_bys = grid%v_bys + model_grid%g_v_bys grid%v_bye = grid%v_bye + model_grid%g_v_bye ! grid%w_bxs = grid%w_bxs + model_grid%g_w_bxs ! grid%w_bxe = grid%w_bxe + model_grid%g_w_bxe ! grid%w_bys = grid%w_bys + model_grid%g_w_bys ! grid%w_bye = grid%w_bye + model_grid%g_w_bye grid%t_bxs = grid%t_bxs + model_grid%g_t_bxs grid%t_bxe = grid%t_bxe + model_grid%g_t_bxe grid%t_bys = grid%t_bys + model_grid%g_t_bys grid%t_bye = grid%t_bye + model_grid%g_t_bye grid%mu_bxs = grid%mu_bxs + model_grid%g_mu_bxs grid%mu_bxe = grid%mu_bxe + model_grid%g_mu_bxe grid%mu_bys = grid%mu_bys + model_grid%g_mu_bys grid%mu_bye = grid%mu_bye + model_grid%g_mu_bye grid%ph_bxs = grid%ph_bxs + model_grid%g_ph_bxs grid%ph_bxe = grid%ph_bxe + model_grid%g_ph_bxe grid%ph_bys = grid%ph_bys + model_grid%g_ph_bys grid%ph_bye = grid%ph_bye + model_grid%g_ph_bye grid%moist_bxs = grid%moist_bxs + model_grid%g_moist_bxs grid%moist_bxe = grid%moist_bxe + model_grid%g_moist_bxe grid%moist_bys = grid%moist_bys + model_grid%g_moist_bys grid%moist_bye = grid%moist_bye + model_grid%g_moist_bye grid%u_btxs = grid%u_btxs + model_grid%g_u_btxs grid%u_btxe = grid%u_btxe + model_grid%g_u_btxe grid%u_btys = grid%u_btys + model_grid%g_u_btys grid%u_btye = grid%u_btye + model_grid%g_u_btye grid%v_btxs = grid%v_btxs + model_grid%g_v_btxs grid%v_btxe = grid%v_btxe + model_grid%g_v_btxe grid%v_btys = grid%v_btys + model_grid%g_v_btys grid%v_btye = grid%v_btye + model_grid%g_v_btye ! grid%w_btxs = grid%w_btxs + model_grid%g_w_btxs ! grid%w_btxe = grid%w_btxe + model_grid%g_w_btxe ! grid%w_btys = grid%w_btys + model_grid%g_w_btys ! grid%w_btye = grid%w_btye + model_grid%g_w_btye grid%t_btxs = grid%t_btxs + model_grid%g_t_btxs grid%t_btxe = grid%t_btxe + model_grid%g_t_btxe grid%t_btys = grid%t_btys + model_grid%g_t_btys grid%t_btye = grid%t_btye + model_grid%g_t_btye grid%mu_btxs = grid%mu_btxs + model_grid%g_mu_btxs grid%mu_btxe = grid%mu_btxe + model_grid%g_mu_btxe grid%mu_btys = grid%mu_btys + model_grid%g_mu_btys grid%mu_btye = grid%mu_btye + model_grid%g_mu_btye grid%ph_btxs = grid%ph_btxs + model_grid%g_ph_btxs grid%ph_btxe = grid%ph_btxe + model_grid%g_ph_btxe grid%ph_btys = grid%ph_btys + model_grid%g_ph_btys grid%ph_btye = grid%ph_btye + model_grid%g_ph_btye grid%moist_btxs = grid%moist_btxs + model_grid%g_moist_btxs grid%moist_btxe = grid%moist_btxe + model_grid%g_moist_btxe grid%moist_btys = grid%moist_btys + model_grid%g_moist_btys grid%moist_btye = grid%moist_btye + model_grid%g_moist_btye ENDIF #endif if (trace_use) call da_trace_exit("da_transfer_xatowrf") end subroutine da_transfer_xatowrf