#define WRF_PORT !Future Modifications: !1. do_cam_sulfchem is hardwired in the current code. It should be obtained from ! mz_aerosol_intr subroutine !2. Cloud fraction compuations yeild only 0 or 1 currently, no fractional value for ! cloud fraction module module_cam_mam_cloudchem use shr_kind_mod, only: r8 => shr_kind_r8 use module_cam_support, only: pcnst =>pcnst_runtime, pcols, pver, fieldname_len, & gas_pcnst => gas_pcnst_modal_aero,iam implicit none save private public :: cam_mam_cloudchem_driver public :: cam_mam_cloudchem_inti integer :: synoz_ndx, so4_ndx, h2o_ndx, o2_ndx, o_ndx, hno3_ndx, dst_ndx, cldice_ndx integer :: o3_ndx integer :: het1_ndx logical :: inv_o3, inv_oh, inv_no3, inv_ho2 integer :: id_o3, id_oh, id_no3, id_ho2 integer :: dgnum_idx = 0 integer :: dgnumwet_idx = 0 integer :: wetdens_ap_idx = 0 contains subroutine cam_mam_cloudchem_inti() use mo_setsox, only : sox_inti implicit none !Call initialization for setsox call sox_inti end subroutine cam_mam_cloudchem_inti subroutine cam_mam_cloudchem_driver( & !Intent Outs dvmrdt_sv13d,dvmrcwdt_sv13d, & !Intent in-outs chem, & !Intent ins moist, scalar, p8w, prain3d, p_phy, & t_phy, dtstepc, ktau,alt, f_ice_phy, & f_rain_phy,cldfra, cldfra_mp_all, & cldfrai, cldfral, is_CAMMGMP_used, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !!----------------------------------------------------------------------- !! ... Chem_solver advances the volumetric mixing ratio !! forward one time step via a combination of explicit, !! ebi, hov, fully implicit, and/or rodas algorithms. !!----------------------------------------------------------------------- use module_configure, only: grid_config_rec_type use module_state_description, only: num_moist, num_chem, p_qv, p_qc, & p_qi, p_qs, p_qnc, p_qni, param_first_scalar, num_scalar, f_qc, & f_qi, f_qv, f_qs use module_cam_support, only: pcnst =>pcnst_runtime, pcols, pver, & pcnst_non_chem => pcnst_non_chem_modal_aero, nfs, & gas_pcnst => gas_pcnst_modal_aero, & gas_pcnst_pos => gas_pcnst_modal_aero_pos use constituents, only: cnst_get_ind use module_data_cam_mam_asect, only: lptr_chem_to_q, lptr_chem_to_qqcw, & factconv_chem_to_q, mw_q_array use physconst, only: mwdry, avogad use mo_setsox, only: setsox, has_sox use modal_aero_data, only: ntot_amode use infnan, only: nan ! implicit none logical, intent(in) :: is_CAMMGMP_used !Scalar Intent-ins integer, intent(in) :: ktau !Time step number integer, intent(in) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real, intent(in) :: dtstepc !Chemistry time step in seconds(s) !3D Intent-ins real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: p8w !Hydrostatic Pressure at level interface (Pa) real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: prain3d !Rate of conversion of condensate to precipitation (kg/kg/s) real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: p_phy !Hydrostatic pressure(Pa) real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: t_phy !Temperature (K) real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: alt real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: F_ICE_PHY !Fraction of ice real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: F_RAIN_PHY !Fraction of rain real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfra !cloud fraction real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfra_mp_all !cloud fraction from MGMP micrpphysics real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfrai real, intent(in), dimension( ims:ime, kms:kme, jms:jme ) :: cldfral !4D Intent ins real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: moist !Mixing ratios (kg/kg for mass species ) real, intent(in), dimension( ims:ime, kms:kme, jms:jme, 1:num_scalar ) :: scalar !Mixing ratios (#/kg for number species) !4D Intent-inouts real, intent(inout), dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: chem !Chem array !4D Intent-outs real, intent(out), dimension( ims:ime, kms:kme, jms:jme, 1:gas_pcnst_pos ) :: dvmrdt_sv13d,dvmrcwdt_sv13d !!----------------------------------------------------------------------- !! ... Dummy arguments !!----------------------------------------------------------------------- !Arguments which were Intent-in in the original CAM's interface integer :: lchnk ! chunk index integer :: ncol ! number columns in chunk integer :: imozart ! gas phase start index in q real(r8) :: delt ! timestep (s) real(r8) :: tfld(pcols,kte) ! midpoint temperature (K) real(r8) :: pmid(pcols,kte) ! midpoint pressures (Pa) real(r8) :: pdel(pcols,kte) ! pressure delta about midpoints (Pa) real(r8) :: cldw(pcols,kte) ! cloud water (kg/kg) real(r8) :: ncldwtr(pcols,kte) ! droplet number concentration (#/kg) !!----------------------------------------------------------------------- !! ... Local variables !!----------------------------------------------------------------------- integer :: i, k, m, n real(r8) :: invariants(pcols,kte,nfs) real(r8) :: vmr(pcols,kte,gas_pcnst) ! xported species (vmr) real(r8) :: vmrcw(pcols,kte,gas_pcnst) ! cloud-borne aerosol (vmr) real(r8), dimension(pcols,kte) :: & mbar ! mean wet atmospheric mass ( amu ) real(r8), dimension(pcols,kte) :: & cwat, & ! cloud water mass mixing ratio (kg/kg) cldnum, & ! droplet number concentration (#/kg) cldfr, & ! cloud fraction prain real(r8) :: qh2o(pcols,kte) ! specific humidity (kg/kg) real(r8) :: dvmrdt_sv1(pcols,pver,gas_pcnst_pos) real(r8) :: dvmrcwdt_sv1(pcols,pver,gas_pcnst_pos) !Variables needed for porting CAM parameterization 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 integer :: imozart_m1, icol, itsm1, itile_len integer :: iw, jw, kw, ktep1, kflip, l, l2, l3 integer :: l_mo_h2so4, l_mo_soag, p1st, ichem real(r8) :: dp, multFrc, fconv real(r8) :: xhnm(pcols,kte) real(r8) :: state_q(pcols,kte,pcnst) real(r8) :: qqcw(pcols,kte,pcnst) !cloud-borne aerosol !Time step for chemistry (following module_cam_mam_aerchem_driver.F) delt = dtstepc pver = kte !Following imozart computation was directly taken from module_cam_mam_aerchem_driver.F imozart_m1 = pcnst_non_chem imozart = imozart_m1 + 1 !Following h2so4 and soag computations were directly taken from module_cam_mam_aerchem_driver.F call cnst_get_ind( 'h2so4', l_mo_h2so4, .false. ) l_mo_h2so4 = l_mo_h2so4 - imozart_m1 if ((l_mo_h2so4 < 1) .or. (l_mo_h2so4 > gas_pcnst)) & call wrf_error_fatal('cam_mam_cloudchem error -- no h2so4 species' ) write(*,*) 'l_mo_h2so4 = ', l_mo_h2so4 call cnst_get_ind( 'soag', l_mo_soag, .false. ) l_mo_soag = l_mo_soag - imozart_m1 if ((l_mo_soag < 1) .or. (l_mo_soag > gas_pcnst)) & call wrf_error_fatal( 'cam_mam_cloudchem error -- no soag species' ) write(*,*) 'l_mo_soag = ', l_mo_soag !Required assignments p1st = param_first_scalar ! Obtain CHEM array's first element's index ncol = pcols icol = ncol !Used in some CAM variables !This subroutine requires that ncol == 1 if(ncol .NE. 1) then call wrf_error_fatal('Number of CAM Columns (NCOL) in CAM_MAM_CLOUDCHEM scheme must be 1') endif !Divide domain in chuncks and map WRF variables into CAM !Loop counters are named iw,jw,kw to represent that they index WRF sided variables itsm1 = its - 1 itile_len = ite - itsm1 do jw = jts , jte do iw = its , ite lchnk = (jw - jts) * itile_len + (iw - itsm1) !1-D index location from a 2-D tile ktep1 = kte + 1 !Flip vertically quantities computed at the mid points do kw = kts, kte kflip = ktep1 - kw pmid(1,kflip) = p_phy(iw,kw,jw) !Pressure at the mid-points (Pa) [state%pmid in CAM] dp = p8w(iw,kw,jw) - p8w(iw,kw+1,jw) !Change in pressure (Pa) pdel(1,kflip) = dp tfld(1,kflip) = t_phy(iw,kw,jw) !Temprature at the mid points (K) [state%t in CAM] !Following three formulas are obtained from ported CAM's ZM cumulus scheme !Values of 0 cause a crash in entropy multFrc = 1._r8/(1._r8 + moist(iw,kw,jw,P_QV)) state_q(1,kflip,1) = max( moist(iw,kw,jw,P_QV)*multFrc, 1.0e-30_r8 ) !Specific humidity [state%q(:,:,1) in CAM] state_q(1,kflip,2) = moist(iw,kw,jw,P_QC)*multFrc !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM] state_q(1,kflip,3) = moist(iw,kw,jw,P_QI)*multFrc !cloud ice [state%q(:,:,3) in CAM] state_q(1,kflip,4) = scalar(iw,kw,jw,P_QNC)*multFrc !Liquid cloud number state_q(1,kflip,5) = scalar(iw,kw,jw,P_QNI)*multFrc !Ice cloud number !Followiing formulas are obtained from Chemistry.F90 of CAM qh2o(1,kflip) = state_q(1,kflip,1) cwat(1,kflip) = state_q(1,kflip,2) + state_q(1,kflip,3) cldnum(1,kflip) = state_q(1,kflip,4) !populate state_q and qqcw arrays !Following Do-Loop is obtained from chem/module_cam_mam_aerchem_driver.F do l = p1st, num_chem l2 = lptr_chem_to_q(l) if ((l2 >= 1) .and. (l2 <= pcnst)) then state_q(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l) end if l2 = lptr_chem_to_qqcw(l) if ((l2 >= 1) .and. (l2 <= pcnst)) then qqcw(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l) !Cloud borne aerosols end if end do ! l prain(1,kflip) = prain3d(iw,kw,jw) !Rate of conversion of condensate to precipitation (kg/kg/s) if(is_CAMMGMP_used) then cldfr(1,kflip) = cldfral(iw,kw,jw) !Cloud fraction from CAMMGMP else cldfr(1,kflip) = cldfra(iw,kw,jw) !Cloud fraction endif cldfr(1,kflip) = min(max(cldfr(1,kflip),0._r8),1._r8) !invariants array is NEVER used in the computations when state_q is defined. Therefore it is set to nan invariants(1,kflip,:) = nan mbar(1,kflip) = mwdry !xhnm is air density in molecules/cm3 (Formulated by RCE, codded by balwinder.singh@pnnl.gov) xhnm(1,kflip) = (avogad*1.0e-6_r8)/(mwdry*alt(iw,kw,jw)) !Note that avogad is in kmoles/moles !initialize dvmrdt_sv13d and dvmrcwdt_sv13d to zero dvmrdt_sv13d(iw,kw,jw,:) = 0.0 dvmrcwdt_sv13d(iw,kw,jw,:) = 0.0 !initialize vmr and vmrc vmr(icol,kflip,:) = 0.0_r8 vmrcw(icol,kflip,:) = 0.0_r8 !Following vmr and vmrcw computation was directly taken from module_cam_mam_aerchem_driver.F !BSINGH - I have changed the loop counters and structure to avoid accessing qqcw array for !indices where qqcw is NOT defined (qqcw is defined only for _c1, _c2 and _c3 aerosols, rest !of the array has junk values) !do l2 = pcnst_non_chem + 1, pcnst !Original loop counters do ichem = p1st , num_chem l2 = lptr_chem_to_q(ichem) if ((l2 >= 1) .and. (l2 <= pcnst)) then l3 = l2 - pcnst_non_chem fconv = mwdry/mw_q_array(l2) vmr(icol,kflip,l3) = state_q(icol,kflip,l2)*fconv endif if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_h2so4) then write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3) endif if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_soag) then write(*,'(a,1p,2e11.3)') '1,1,1 soag q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3) endif l2 = lptr_chem_to_qqcw(ichem) if ((l2 >= 1) .and. (l2 <= pcnst)) then l3 = l2 - pcnst_non_chem fconv = mwdry/mw_q_array(l2) vmrcw(icol,kflip,l3) = qqcw(icol,kflip,l2)*fconv endif end do enddo !enddo for kw=kts,kte loop dvmrdt_sv1 = vmr dvmrcwdt_sv1 = vmrcw !Note: Only vmrcw and vm are the outputs from the setsox call if( has_sox .and. (.not. do_cam_sulfchem) ) then call setsox( ncol, & pmid, & delt, & tfld, & qh2o, & cwat, & lchnk, & pdel, & mbar, & prain, & cldfr, & cldnum, & vmrcw, & imozart-1,& xhnm, & !In original call, invariants(1,1,indexm), is being passed but it is replaced here with xhnm vmr, & invariants ) endif dvmrdt_sv1 = (vmr - dvmrdt_sv1)/delt dvmrcwdt_sv1 = (vmrcw - dvmrcwdt_sv1)/delt !Post processing of the output from CAM's parameterization do l2 = pcnst_non_chem+1, pcnst l3 = l2 - pcnst_non_chem fconv = mw_q_array(l2)/mwdry state_q(icol,:,l2) = vmr(icol,:,l3)*fconv qqcw(icol,:,l2) = vmrcw(icol,:,l3)*fconv if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_h2so4) & write(*,'(a,1p,2e11.3)') '1,1,1 h2so4 q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3) if (iw*jw == 1 .and. kw == kts .and. l3 == l_mo_soag) & write(*,'(a,1p,2e11.3)') '1,1,1 soag q8 & vmr8', state_q(icol,pver,l2), vmr(icol,pver,l3) end do do kw = kts , kte kflip = kte-kw+1 do l = p1st, num_chem l2 = lptr_chem_to_q(l) if ((l2 >= 1) .and. (l2 <= pcnst)) then chem(iw,kw,jw,l) = state_q(1,kflip,l2)/factconv_chem_to_q(l) end if l2 = lptr_chem_to_qqcw(l) if ((l2 >= 1) .and. (l2 <= pcnst)) then chem(iw,kw,jw,l) = qqcw(1,kflip,l2)/factconv_chem_to_q(l) end if end do ! l do l = 1 , gas_pcnst dvmrdt_sv13d(iw,kw,jw,l) = dvmrdt_sv1(1,kflip,l) dvmrcwdt_sv13d(iw,kw,jw,l) = dvmrcwdt_sv1(1,kflip,l) enddo end do !kw post processing do -loop enddo !iw do-loop enddo !jw do-loop end subroutine cam_mam_cloudchem_driver end module module_cam_mam_cloudchem