chem/module_cam_mam_gas_wetdep_driver.F

   1 #define WRF_PORT
   2 !Future Modifications:
   3
   4 module module_cam_mam_gas_wetdep_driver
   5
   6   use shr_kind_mod,       only: r8 => shr_kind_r8
   7   use module_cam_support, only: pcnst =>pcnst_runtime, pcols, pver
   8
   9   implicit none
  10   save
  11
  12   private
  13   public :: cam_mam_gas_wetdep_driver
  14   public :: cam_mam_gas_wetdep_inti
  15
  16 contains
  17
  18   subroutine cam_mam_gas_wetdep_inti()
  19     use mo_sethet, only : sethet_inti
  20     implicit none
  21
  22     !Call initialization for sethet
  23     call  sethet_inti
  24
  25   end subroutine cam_mam_gas_wetdep_inti
  26
  27   subroutine cam_mam_gas_wetdep_driver(          &
  28        !Intent in-outs
  29        chem,                                     &
  30        !Intent ins
  31        dtstepc, config_flags, ht, XLAT, nevapr3d,&
  32        rprdsh3d, rprddp3d, prain3d, z_sea_level, &
  33        p_phy, t_phy, alt, moist, scalar,         &
  34        ids,ide, jds,jde, kds,kde,                &
  35        ims,ime, jms,jme, kms,kme,                &
  36        its,ite, jts,jte, kts,kte                 )
  37
  38
  39     use module_state_description,  only: num_moist, num_chem, p_qv, p_qc,    &
  40          p_qi, p_qnc, p_qni, param_first_scalar, num_scalar, CAMUWSHCUSCHEME,&
  41          CAMZMSCHEME
  42
  43     use module_cam_support,        only: pcnst =>pcnst_runtime, pcols, pver, &
  44          pcnst_non_chem => pcnst_non_chem_modal_aero, nfs,                   &
  45          gas_pcnst => gas_pcnst_modal_aero, hetcnt => gas_pcnst_modal_aero !Hetcnt is same as gas_pcnst in CAM5.1
  46
  47     use module_data_cam_mam_asect, only: lptr_chem_to_q, lptr_chem_to_qqcw,  &
  48          factconv_chem_to_q, mw_q_array
  49     use physconst,                 only: mwdry, avogad, gravit, rga
  50
  51     use mo_sethet,                 only: sethet
  52     use infnan,                    only: nan
  53     use module_configure,          only: grid_config_rec_type
  54
  55     !
  56     implicit none
  57
  58     TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  59
  60     !Scalar Intent-ins
  61     integer, intent(in) ::          &
  62          ids,ide, jds,jde, kds,kde, &
  63          ims,ime, jms,jme, kms,kme, &
  64          its,ite, jts,jte, kts,kte
  65
  66     real, intent(in) :: dtstepc       !Chemistry time step in seconds(s)
  67
  68     !2D Intent-ins
  69     real, dimension( ims:ime, jms:jme ), intent(in) :: ht     !Terrain height (m)
  70     real, dimension( ims:ime, jms:jme ), intent(in) :: XLAT   !Latitude
  71
  72     !3D Intent-ins
  73     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: prain3d     !Rate of conversion of condensate to precipitation (kg/kg/s)
  74     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: p_phy       !Hydrostatic pressure(Pa)
  75     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: t_phy       !Temperature (K)
  76     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: alt         !inverse density(m3/kg)
  77     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: z_sea_level !Height above sea level at mid-level (m)
  78     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: nevapr3d    !Evaporation rate of rain + snow (kg/kg/s)
  79     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: rprdsh3d    !dq/dt due to deep and shallow convective rainout(kg/kg/s)
  80     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: rprddp3d    !dq/dt due to deep convective rainout (kg/kg/s)
  81
  82     !4D Intent ins
  83     real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_moist )  :: moist  !Mixing ratios (kg/kg for mass species )
  84     real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_scalar ) :: scalar !Mixing ratios (#/kg for number species)
  85
  86     !4D Intent-inouts
  87     real, intent(inout), dimension( ims:ime, kms:kme, jms:jme, 1:num_chem )   :: chem !Chem array
  88
  89
  90     !!-----------------------------------------------------------------------
  91     !!        ... Dummy arguments
  92     !!-----------------------------------------------------------------------
  93     !Arguments which were Intent-in in the original CAM's interface
  94     integer            :: lchnk                         ! chunk index
  95     integer            :: ncol                          ! number columns in chunk
  96     real(r8)           :: delt                          ! timestep (s)
  97     real(r8)           :: tfld(pcols,kte)               ! midpoint temperature (K)
  98     real(r8)           :: pmid(pcols,kte)               ! midpoint pressures (Pa)
  99
 100     !!-----------------------------------------------------------------------
 101     !!          ... Local variables
 102     !!-----------------------------------------------------------------------
 103     real(r8), parameter :: m2km  = 1.e-3_r8
 104
 105     integer  :: i, k, m, n
 106     real(r8) :: zsurf(pcols)                          ! surface height (m)
 107     real(r8) :: phis(pcols)
 108     real(r8) :: cmfdqr(pcols, kte)
 109     real(r8) :: nevapr(pcols, kte)
 110     real(r8) :: rprdsh(pcols, kte)
 111     real(r8) :: rprddp(pcols, kte)
 112     real(r8) :: state_zm(pcols,kte)
 113     real(r8) :: zmid(pcols,kte)
 114     real(r8) :: invariants(pcols,kte,nfs)
 115     real(r8) :: vmr(pcols,kte,gas_pcnst)              ! xported species (vmr)
 116     real(r8) :: het_rates(pcols,kte,hetcnt)           !washout rate (1/s)
 117     real(r8) :: prain(pcols,kte)
 118
 119     !Variables needed for porting CAM parameterization in WRF
 120     integer  ::  icol, itsm1, itile_len
 121     integer  ::  iw, jw, kw, ktep1, kflip, l, l2, l3
 122     integer  ::  p1st, ichem
 123     real(r8) ::  dp, multFrc, fconv
 124     real(r8) ::  rlat(pcols)
 125     real(r8) ::  xhnm(pcols,kte)
 126     real(r8) ::  state_q(pcols,kte,pcnst)
 127
 128
 129     !Time step for chemistry (following module_cam_mam_aerchem_driver.F)
 130     delt = dtstepc ! converts to r8
 131     pver = kte
 132
 133     !Required assignments
 134     p1st = param_first_scalar ! Obtain CHEM array's first element's index
 135
 136     ncol = pcols
 137     icol = ncol !Used in some CAM variables
 138
 139     !This subroutine requires that ncol == 1
 140     if(ncol .NE. 1) then
 141        call wrf_error_fatal('Number of CAM Columns (NCOL) in CAM_MAM_CLOUDCHEM scheme must be 1')
 142     endif
 143
 144     !Divide domain in chuncks and map WRF variables into CAM
 145     !Loop counters are named iw,jw,kw to represent that they index WRF sided variables
 146
 147     itsm1     = its - 1
 148     itile_len = ite - itsm1
 149     do jw     = jts , jte
 150        do iw  = its , ite
 151
 152           lchnk   = (jw - jts) * itile_len + (iw - itsm1)             !1-D index location from a 2-D tile
 153           ktep1   = kte + 1
 154
 155           phis(1) = ht(iw,jw)  * gravit
 156
 157           !Flip vertically quantities computed at the mid points
 158           do kw  = kts, kte
 159              kflip             = ktep1 - kw
 160              pmid(1,kflip)     = p_phy(iw,kw,jw)                   !Pressure at the mid-points (Pa) [state%pmid in CAM]
 161              tfld(1,kflip)     = t_phy(iw,kw,jw)                   !Temprature at the mid points (K) [state%t in CAM]
 162              state_zm(1,kflip) = z_sea_level(iw,kw,jw) - ht(iw,jw) !Height above the ground at the midpoints (m) [state%zm in CAM]
 163              rlat(1)           = xlat(iw,jw)
 164
 165              !Following three formulas are obtained from ported CAM's ZM cumulus scheme
 166              !Values of 0 cause a crash in entropy
 167              multFrc              = 1._r8/(1._r8 + moist(iw,kw,jw,P_QV))
 168              state_q(1,kflip,1)   = max( moist(iw,kw,jw,P_QV)*multFrc, 1.0e-30_r8 ) !Specific humidity                       [state%q(:,:,1) in CAM]
 169              state_q(1,kflip,2)   = moist(iw,kw,jw,P_QC)*multFrc                    !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
 170              state_q(1,kflip,3)   = moist(iw,kw,jw,P_QI)*multFrc                    !cloud ice                               [state%q(:,:,3) in CAM]
 171              state_q(1,kflip,4)   = scalar(iw,kw,jw,P_QNC)*multFrc                  !Liquid cloud number
 172              state_q(1,kflip,5)   = scalar(iw,kw,jw,P_QNI)*multFrc                  !Ice cloud number
 173
 174              !populate state_q
 175              !Following Do-Loop is obtained from chem/module_cam_mam_aerchem_driver.F
 176              do l = p1st, num_chem
 177                 l2 = lptr_chem_to_q(l)
 178                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 179                    state_q(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)
 180                 end if
 181              end do ! l
 182              prain(1,kflip)        = prain3d(iw,kw,jw)
 183              nevapr(1,kflip)       = nevapr3d(iw,kw,jw)                           !Evaporation rate of rain + snow (kg/kg/s)
 184
 185              rprdsh(1,kflip)            = 0.0_r8
 186
 187              if(config_flags%shcu_physics==CAMUWSHCUSCHEME) then
 188                 !inputs from shallow convection
 189                 rprdsh(1,kflip)         = rprdsh3d(iw,kw,jw)                     !dq/dt due to deep and shallow convective rainout(kg/kg/s)
 190              endif
 191
 192              rprddp(1,kflip)            = 0.0_r8
 193              if(config_flags%cu_physics==CAMZMSCHEME)then
 194                 !inputs from deep convection
 195                 rprddp(1,kflip)         = rprddp3d(iw,kw,jw)                     !dq/dt due to deep convective rainout (kg/kg/s)
 196              endif
 197
 198              !xhnm is air density in molecules/cm3 (Formulated by RCE, codded by balwinder.singh@pnnl.gov)
 199              xhnm(1,kflip) = (avogad*1.0e-6_r8)/(mwdry*alt(iw,kw,jw)) !Note that avogad is in kmoles/moles
 200
 201              !invariants array is NEVER used in the computations when state_q is defined. Therefore it is set to nan
 202              invariants(1,kflip,:) = nan
 203
 204              do ichem = p1st , num_chem
 205                 l2 = lptr_chem_to_q(ichem)
 206                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 207                    l3                   = l2 - pcnst_non_chem
 208                    fconv                = mwdry/mw_q_array(l2)
 209                    vmr(icol,kflip,l3)   = state_q(icol,kflip,l2)*fconv
 210                 endif
 211              end do
 212
 213           enddo !enddo for kw=kts, kte loop
 214
 215           zsurf(:ncol) = rga * phis(:ncol)
 216           do k = 1,pver
 217              zmid(:ncol,k) = m2km * (state_zm(:ncol,k) + zsurf(:ncol))
 218           end do
 219
 220
 221           cmfdqr(:ncol,:) = rprddp(:ncol,:)  + rprdsh(:ncol,:)
 222           !Note: Only het_rates is the output from the sethet call
 223           call sethet( het_rates, pmid, zmid, phis, tfld, &
 224                    cmfdqr, prain, nevapr, delt, xhnm, & !In original call, invariants(1,1,indexm), is being passed but it is replaced here with xhnm
 225                    vmr, ncol, lchnk,rlat)
 226
 227           do l2 = pcnst_non_chem+1, pcnst
 228              l3                 = l2 - pcnst_non_chem
 229              het_rates(icol,:,l3) = min(max(0._r8, het_rates(icol,:,l3)),(1._r8-1.e-5_r8)/delt) !BSINGH - PMA advised to introduce limit on het_rates to keep state_q positive
 230              state_q(icol,:,l2) = state_q(icol,:,l2) - het_rates(icol,:,l3)*state_q(icol,:,l2)*delt
 231           end do
 232
 233           !Post processing of the output from CAM's parameterization
 234           do kw = kts , kte
 235              kflip = kte-kw+1
 236
 237              do l = p1st, num_chem
 238                 l2 = lptr_chem_to_q(l)
 239                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 240                    chem(iw,kw,jw,l) = state_q(1,kflip,l2)/factconv_chem_to_q(l)
 241                 end if
 242              end do ! l
 243           end do !kw post processing do -loop
 244
 245
 246        enddo !iw do-loop
 247     enddo    !jw do-loop
 248
 249
 250
 251   end subroutine cam_mam_gas_wetdep_driver
 252
 253 end module module_cam_mam_gas_wetdep_driver