chem/module_cam_mam_cloudchem.F

   1 #define WRF_PORT
   2 !Future Modifications:
   3 !1. do_cam_sulfchem is hardwired in the current code. It should be obtained from
   4 !   mz_aerosol_intr subroutine
   5 !2. Cloud fraction compuations yeild only 0 or 1 currently, no fractional value for
   6 !   cloud fraction
   7
   8 module module_cam_mam_cloudchem
   9
  10   use shr_kind_mod,       only: r8 => shr_kind_r8
  11   use module_cam_support, only: pcnst =>pcnst_runtime, pcols, pver, fieldname_len, &
  12        gas_pcnst => gas_pcnst_modal_aero,iam
  13
  14   implicit none
  15   save
  16
  17   private
  18   public :: cam_mam_cloudchem_driver
  19   public :: cam_mam_cloudchem_inti
  20
  21   integer :: synoz_ndx, so4_ndx, h2o_ndx, o2_ndx, o_ndx, hno3_ndx, dst_ndx, cldice_ndx
  22   integer :: o3_ndx
  23   integer :: het1_ndx
  24   logical :: inv_o3, inv_oh, inv_no3, inv_ho2
  25   integer :: id_o3, id_oh, id_no3, id_ho2
  26   integer :: dgnum_idx       = 0
  27   integer :: dgnumwet_idx    = 0
  28   integer :: wetdens_ap_idx  = 0
  29
  30 contains
  31
  32   subroutine cam_mam_cloudchem_inti()
  33     use mo_setsox, only : sox_inti
  34     implicit none
  35
  36     !Call initialization for setsox
  37     call  sox_inti
  38
  39   end subroutine cam_mam_cloudchem_inti
  40
  41   subroutine cam_mam_cloudchem_driver(           &
  42        !Intent Outs
  43        dvmrdt_sv13d,dvmrcwdt_sv13d,              &
  44        !Intent in-outs
  45        chem,                                     &
  46        !Intent ins
  47        moist, scalar, p8w, prain3d, p_phy,       &
  48        t_phy, dtstepc, ktau,alt, f_ice_phy,      &
  49        f_rain_phy,cldfra, cldfra_mp_all,         &
  50        cldfrai, cldfral, is_CAMMGMP_used,        &
  51        ids,ide, jds,jde, kds,kde,                &
  52        ims,ime, jms,jme, kms,kme,                &
  53        its,ite, jts,jte, kts,kte                 )
  54     !!-----------------------------------------------------------------------
  55     !!     ... Chem_solver advances the volumetric mixing ratio
  56     !!         forward one time step via a combination of explicit,
  57     !!         ebi, hov, fully implicit, and/or rodas algorithms.
  58     !!-----------------------------------------------------------------------
  59     use module_configure,          only: grid_config_rec_type
  60     use module_state_description,  only: num_moist, num_chem, p_qv, p_qc,    &
  61          p_qi, p_qs, p_qnc, p_qni, param_first_scalar, num_scalar, f_qc,     &
  62          f_qi,  f_qv, f_qs
  63     use module_cam_support,        only: pcnst =>pcnst_runtime, pcols, pver, &
  64          pcnst_non_chem => pcnst_non_chem_modal_aero, nfs,                   &
  65          gas_pcnst => gas_pcnst_modal_aero,                                  &
  66          gas_pcnst_pos => gas_pcnst_modal_aero_pos
  67     use constituents,              only: cnst_get_ind
  68     use module_data_cam_mam_asect, only: lptr_chem_to_q, lptr_chem_to_qqcw,  &
  69          factconv_chem_to_q, mw_q_array
  70     use physconst,                 only: mwdry, avogad
  71
  72     use mo_setsox,                 only: setsox, has_sox
  73     use modal_aero_data,           only: ntot_amode
  74     use infnan,                    only: nan
  75
  76     !
  77     implicit none
  78
  79     logical, intent(in) :: is_CAMMGMP_used
  80     !Scalar Intent-ins
  81     integer, intent(in) :: ktau       !Time step number
  82     integer, intent(in) ::          &
  83          ids,ide, jds,jde, kds,kde, &
  84          ims,ime, jms,jme, kms,kme, &
  85          its,ite, jts,jte, kts,kte
  86
  87     real, intent(in) :: dtstepc       !Chemistry time step in seconds(s)
  88
  89     !3D Intent-ins
  90     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: p8w     !Hydrostatic Pressure at level interface (Pa)
  91     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: prain3d !Rate of conversion of condensate to precipitation (kg/kg/s)
  92     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: p_phy   !Hydrostatic pressure(Pa)
  93     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: t_phy   !Temperature (K)
  94     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: alt
  95     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: F_ICE_PHY   !Fraction of ice
  96     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: F_RAIN_PHY  !Fraction of rain
  97     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfra   !cloud fraction
  98     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfra_mp_all   !cloud fraction from MGMP micrpphysics
  99     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfrai
 100     real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfral
 101     !4D Intent ins
 102     real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_moist )  :: moist  !Mixing ratios (kg/kg for mass species )
 103     real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_scalar ) :: scalar !Mixing ratios (#/kg for number species)
 104     !4D Intent-inouts
 105     real, intent(inout), dimension( ims:ime, kms:kme, jms:jme, 1:num_chem )   :: chem !Chem array
 106     !4D Intent-outs
 107     real, intent(out), dimension( ims:ime, kms:kme, jms:jme, 1:gas_pcnst_pos )   :: dvmrdt_sv13d,dvmrcwdt_sv13d
 108
 109
 110     !!-----------------------------------------------------------------------
 111     !!        ... Dummy arguments
 112     !!-----------------------------------------------------------------------
 113     !Arguments which were Intent-in in the original CAM's interface
 114     integer            :: lchnk                         ! chunk index
 115     integer            :: ncol                          ! number columns in chunk
 116     integer            :: imozart                       ! gas phase start index in q
 117     real(r8)           :: delt                          ! timestep (s)
 118     real(r8)           :: tfld(pcols,kte)               ! midpoint temperature (K)
 119     real(r8)           :: pmid(pcols,kte)               ! midpoint pressures (Pa)
 120     real(r8)           :: pdel(pcols,kte)               ! pressure delta about midpoints (Pa)
 121     real(r8)           :: cldw(pcols,kte)               ! cloud water (kg/kg)
 122     real(r8)           :: ncldwtr(pcols,kte)            ! droplet number concentration (#/kg)
 123
 124     !!-----------------------------------------------------------------------
 125     !!          ... Local variables
 126     !!-----------------------------------------------------------------------
 127     integer      ::  i, k, m, n
 128     real(r8)     ::  invariants(pcols,kte,nfs)
 129     real(r8)     ::  vmr(pcols,kte,gas_pcnst)              ! xported species (vmr)
 130     real(r8)     ::  vmrcw(pcols,kte,gas_pcnst)            ! cloud-borne aerosol (vmr)
 131     real(r8), dimension(pcols,kte) :: &
 132          mbar                                           ! mean wet atmospheric mass ( amu )
 133     real(r8), dimension(pcols,kte) :: &
 134          cwat, &                                           ! cloud water mass mixing ratio (kg/kg)
 135          cldnum, &                                         ! droplet number concentration (#/kg)
 136          cldfr, &                                          ! cloud fraction
 137          prain
 138
 139     real(r8) :: qh2o(pcols,kte)               ! specific humidity (kg/kg)
 140     real(r8) :: dvmrdt_sv1(pcols,pver,gas_pcnst_pos)
 141     real(r8) :: dvmrcwdt_sv1(pcols,pver,gas_pcnst_pos)
 142
 143     !Variables needed for porting CAM parameterization
 144
 145     logical, parameter :: do_cam_sulfchem = .FALSE. !Forced it to FALSE so that setsox can execute,in CAM it is obtained from mz_aerosols_intr.F90
 146
 147     integer      ::  imozart_m1, icol, itsm1, itile_len
 148     integer      ::  iw, jw, kw, ktep1, kflip, l, l2, l3
 149     integer      ::  l_mo_h2so4, l_mo_soag, p1st, ichem
 150     real(r8)     ::  dp, multFrc, fconv
 151     real(r8)     ::  xhnm(pcols,kte)
 152     real(r8)     ::  state_q(pcols,kte,pcnst)
 153     real(r8)     ::  qqcw(pcols,kte,pcnst)      !cloud-borne aerosol
 154
 155
 156     !Time step for chemistry (following module_cam_mam_aerchem_driver.F)
 157     delt = dtstepc
 158     pver = kte
 159
 160
 161     !Following imozart computation was directly taken from module_cam_mam_aerchem_driver.F
 162     imozart_m1 = pcnst_non_chem
 163     imozart = imozart_m1 + 1
 164
 165     !Following h2so4 and soag computations were directly taken from module_cam_mam_aerchem_driver.F
 166     call cnst_get_ind( 'h2so4', l_mo_h2so4, .false. )
 167     l_mo_h2so4 = l_mo_h2so4 - imozart_m1
 168     if ((l_mo_h2so4 < 1) .or. (l_mo_h2so4 > gas_pcnst)) &
 169          call wrf_error_fatal('cam_mam_cloudchem error -- no h2so4 species' )
 170     write(*,*) 'l_mo_h2so4 = ', l_mo_h2so4
 171
 172     call cnst_get_ind( 'soag', l_mo_soag, .false. )
 173     l_mo_soag = l_mo_soag - imozart_m1
 174     if ((l_mo_soag < 1) .or. (l_mo_soag > gas_pcnst)) &
 175          call  wrf_error_fatal( 'cam_mam_cloudchem error -- no soag species' )
 176     write(*,*) 'l_mo_soag = ', l_mo_soag
 177
 178
 179     !Required assignments
 180     p1st = param_first_scalar ! Obtain CHEM array's first element's index
 181
 182
 183     ncol = pcols
 184     icol = ncol !Used in some CAM variables
 185
 186     !This subroutine requires that ncol == 1
 187     if(ncol .NE. 1) then
 188        call wrf_error_fatal('Number of CAM Columns (NCOL) in CAM_MAM_CLOUDCHEM scheme must be 1')
 189     endif
 190
 191     !Divide domain in chuncks and map WRF variables into CAM
 192     !Loop counters are named iw,jw,kw to represent that they index WRF sided variables
 193
 194     itsm1     = its - 1
 195     itile_len = ite - itsm1
 196     do jw     = jts , jte
 197        do iw  = its , ite
 198
 199           lchnk   = (jw - jts) * itile_len + (iw - itsm1)             !1-D index location from a 2-D tile
 200           ktep1   = kte + 1
 201
 202           !Flip vertically quantities computed at the mid points
 203           do kw  = kts, kte
 204              kflip          = ktep1 - kw
 205              pmid(1,kflip)  = p_phy(iw,kw,jw)                   !Pressure at the mid-points (Pa) [state%pmid in CAM]
 206              dp             = p8w(iw,kw,jw) - p8w(iw,kw+1,jw)   !Change in pressure (Pa)
 207              pdel(1,kflip)  = dp
 208              tfld(1,kflip)  = t_phy(iw,kw,jw)                   !Temprature at the mid points (K) [state%t in CAM]
 209
 210              !Following three formulas are obtained from ported CAM's ZM cumulus scheme
 211              !Values of 0 cause a crash in entropy
 212              multFrc              = 1._r8/(1._r8 + moist(iw,kw,jw,P_QV))
 213              state_q(1,kflip,1)   = max( moist(iw,kw,jw,P_QV)*multFrc, 1.0e-30_r8 ) !Specific humidity                       [state%q(:,:,1) in CAM]
 214              state_q(1,kflip,2)   = moist(iw,kw,jw,P_QC)*multFrc                    !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
 215              state_q(1,kflip,3)   = moist(iw,kw,jw,P_QI)*multFrc                    !cloud ice                               [state%q(:,:,3) in CAM]
 216              state_q(1,kflip,4)   = scalar(iw,kw,jw,P_QNC)*multFrc                  !Liquid cloud number
 217              state_q(1,kflip,5)   = scalar(iw,kw,jw,P_QNI)*multFrc                  !Ice cloud number
 218
 219              !Followiing formulas are obtained from Chemistry.F90 of CAM
 220              qh2o(1,kflip)        = state_q(1,kflip,1)
 221              cwat(1,kflip)        = state_q(1,kflip,2) + state_q(1,kflip,3)
 222              cldnum(1,kflip)      = state_q(1,kflip,4)
 223
 224              !populate state_q and qqcw arrays
 225              !Following Do-Loop is obtained from chem/module_cam_mam_aerchem_driver.F
 226              do l = p1st, num_chem
 227                 l2 = lptr_chem_to_q(l)
 228                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 229                    state_q(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)
 230                 end if
 231                 l2 = lptr_chem_to_qqcw(l)
 232                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 233                    qqcw(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)     !Cloud borne aerosols
 234                 end if
 235              end do ! l
 236
 237              prain(1,kflip)        = prain3d(iw,kw,jw)                      !Rate of conversion of condensate to precipitation (kg/kg/s)
 238              if(is_CAMMGMP_used) then
 239                 cldfr(1,kflip)        = cldfral(iw,kw,jw)                       !Cloud fraction from CAMMGMP
 240              else
 241                 cldfr(1,kflip)        = cldfra(iw,kw,jw)                        !Cloud fraction
 242              endif
 243              cldfr(1,kflip)        = min(max(cldfr(1,kflip),0._r8),1._r8)
 244              !invariants array is NEVER used in the computations when state_q is defined. Therefore it is set to nan
 245              invariants(1,kflip,:) = nan
 246              mbar(1,kflip)         = mwdry
 247
 248              !xhnm is air density in molecules/cm3 (Formulated by RCE, codded by balwinder.singh@pnnl.gov)
 249              xhnm(1,kflip) = (avogad*1.0e-6_r8)/(mwdry*alt(iw,kw,jw)) !Note that avogad is in kmoles/moles
 250              !initialize dvmrdt_sv13d and dvmrcwdt_sv13d to zero
 251              dvmrdt_sv13d(iw,kw,jw,:)   = 0.0
 252              dvmrcwdt_sv13d(iw,kw,jw,:) = 0.0
 253
 254              !initialize vmr and vmrc
 255              vmr(icol,kflip,:)          = 0.0_r8
 256              vmrcw(icol,kflip,:)        = 0.0_r8
 257
 258              !Following vmr and vmrcw computation was directly taken from module_cam_mam_aerchem_driver.F
 259              !BSINGH - I have changed the loop counters and structure to avoid accessing qqcw array for
 260              !indices where qqcw is NOT defined (qqcw is defined only for _c1, _c2 and _c3 aerosols, rest
 261              !of the array has junk values)
 262
 263              !do l2 = pcnst_non_chem + 1, pcnst !Original loop counters
 264              do ichem = p1st , num_chem
 265                 l2 = lptr_chem_to_q(ichem)
 266                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 267                    l3                   = l2 - pcnst_non_chem
 268                    fconv                = mwdry/mw_q_array(l2)
 269                    vmr(icol,kflip,l3)   = state_q(icol,kflip,l2)*fconv
 270                 endif
 271
 272                 if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_h2so4) then
 273                    write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3)
 274                 endif
 275                 if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_soag) then
 276                    write(*,'(a,1p,2e11.3)') '1,1,1 soag  q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3)
 277                 endif
 278
 279                 l2 = lptr_chem_to_qqcw(ichem)
 280                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 281                    l3                   = l2 - pcnst_non_chem
 282                    fconv                = mwdry/mw_q_array(l2)
 283                    vmrcw(icol,kflip,l3) = qqcw(icol,kflip,l2)*fconv
 284                 endif
 285              end do
 286
 287           enddo !enddo for kw=kts,kte loop
 288
 289           dvmrdt_sv1 = vmr
 290           dvmrcwdt_sv1 = vmrcw
 291
 292           !Note: Only vmrcw and vm are the outputs from the setsox call
 293           if( has_sox .and. (.not. do_cam_sulfchem) ) then
 294              call setsox( ncol,   &
 295                   pmid,   &
 296                   delt,   &
 297                   tfld,   &
 298                   qh2o,   &
 299                   cwat,   &
 300                   lchnk,  &
 301                   pdel,   &
 302                   mbar,   &
 303                   prain,  &
 304                   cldfr,  &
 305                   cldnum, &
 306                   vmrcw,    &
 307                   imozart-1,&
 308                   xhnm, & !In original call, invariants(1,1,indexm), is being passed but it is replaced here with xhnm
 309                   vmr, &
 310                   invariants )
 311
 312           endif
 313           dvmrdt_sv1 = (vmr - dvmrdt_sv1)/delt
 314           dvmrcwdt_sv1 = (vmrcw - dvmrcwdt_sv1)/delt
 315
 316           !Post processing of the output from CAM's parameterization
 317           do l2 = pcnst_non_chem+1, pcnst
 318              l3                 = l2 - pcnst_non_chem
 319              fconv              = mw_q_array(l2)/mwdry
 320              state_q(icol,:,l2) = vmr(icol,:,l3)*fconv
 321              qqcw(icol,:,l2)    = vmrcw(icol,:,l3)*fconv
 322              if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_h2so4) &
 323                   write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3)
 324              if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_soag) &
 325                   write(*,'(a,1p,2e11.3)') '1,1,1 soag  q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3)
 326           end do
 327
 328
 329           do kw = kts , kte
 330              kflip = kte-kw+1
 331
 332              do l = p1st, num_chem
 333                 l2 = lptr_chem_to_q(l)
 334                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 335                    chem(iw,kw,jw,l)         = state_q(1,kflip,l2)/factconv_chem_to_q(l)
 336                 end if
 337                 l2 = lptr_chem_to_qqcw(l)
 338                 if ((l2 >= 1) .and. (l2 <= pcnst)) then
 339                    chem(iw,kw,jw,l) = qqcw(1,kflip,l2)/factconv_chem_to_q(l)
 340                 end if
 341              end do ! l
 342              do l = 1 , gas_pcnst
 343                 dvmrdt_sv13d(iw,kw,jw,l)   = dvmrdt_sv1(1,kflip,l)
 344                 dvmrcwdt_sv13d(iw,kw,jw,l) = dvmrcwdt_sv1(1,kflip,l)
 345              enddo
 346
 347           end do !kw post processing do -loop
 348
 349
 350        enddo !iw do-loop
 351     enddo    !jw do-loop
 352
 353
 354
 355   end subroutine cam_mam_cloudchem_driver
 356
 357 end module module_cam_mam_cloudchem