MODULE MODULE_GOCART_CHEM CONTAINS subroutine gocart_chem_driver(curr_secs,dt,config_flags, & gmt,julday,t_phy,moist, & chem,rho_phy,dz8w,p8w,backg_oh,backg_h2o2,backg_no3, & gd_cldf,dx,g,xlat,xlong,ttday,tcosz, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_configure USE module_state_description USE module_phot_mad, only : calc_zenith IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: julday, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & INTENT(IN ) :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: & xlat,xlong,ttday,tcosz REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & OPTIONAL, & INTENT(IN ) :: gd_cldf REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: t_phy, & backg_oh,backg_h2o2,backg_no3,dz8w,p8w, & rho_phy REAL(KIND=8), INTENT(IN) :: curr_secs REAL, INTENT(IN ) :: dt,dx,g,gmt integer :: nmx,i,j,k,imx,jmx,lmx real*8, DIMENSION (1,1,1) :: tmp,airden,airmas,oh,xno3,h2o2,chldms_oh, & chldms_no3,chldms_x,chpso2,chpmsa,chpso4, & chlso2_oh,chlso2_aq,cldf real*8, DIMENSION (1,1,4) :: tdry real*8, DIMENSION (1,1) :: cossza real, DIMENSION (1,1) :: sza,cosszax real*8, DIMENSION (1,1,1,4) :: tc,bems real*8, dimension (1) :: dxy real(kind=8) :: xtime,xhour real:: rlat,xlonn real :: zenith,zenita,azimuth,xmin,xtimin,gmtp integer(kind=8) :: ixhour imx=1 jmx=1 lmx=1 nmx=4 tdry=0.d0 xtime=curr_secs/60._8 ixhour=int(gmt+.01,8)+int(xtime/60._8,8) xhour=real(ixhour,8) xmin=60.*gmt+real(xtime-xhour*60._8,8) gmtp=mod(xhour,24._8) gmtp=gmtp+xmin/60. dxy(1)=dx*dx ! ! following arrays for busget stuff only ! ! ! chem_select: SELECT CASE(config_flags%chem_opt) CASE (GOCART_SIMPLE) CALL wrf_debug(15,'calling gocart chemistry ') do j=jts,jte do i=its,ite zenith=0. zenita=0. azimuth=0. rlat=xlat(i,j)*3.1415926535590/180. xlonn=xlong(i,j) CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat) cossza(1,1)=cosszax(1,1) ! do k=kts,kte-1 chldms_oh=0. chldms_no3=0. chldms_x=0. chpso2=0. chpmsa=0. chpso4=0. chlso2_oh=0. chlso2_aq=0. if (present(gd_cldf) ) then cldf(1,1,1)=gd_cldf(i,k,j) else cldf(1,1,1)=0. endif if(p_qc.gt.1 .and. p_qi.gt.1)then if(moist(i,k,j,p_qc).gt.0.or.moist(i,k,j,p_qi).gt.0.)cldf(1,1,1)=1. elseif(p_qc.gt.1 .and. p_qi.le.1)then if(moist(i,k,j,p_qc).gt.0.)cldf(1,1,1)=1. endif tc(1,1,1,1)=chem(i,k,j,p_dms)*1.d-6 tc(1,1,1,2)=chem(i,k,j,p_so2)*1.d-6 tc(1,1,1,3)=chem(i,k,j,p_sulf)*1.d-6 tc(1,1,1,4)=chem(i,k,j,p_msa)*1.d-6 airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*dx*dx/g airden(1,1,1)=rho_phy(i,k,j) tmp(1,1,1)=t_phy(i,k,j) oh(1,1,1)=86400./dt*cossza(1,1)*backg_oh(i,k,j)/tcosz(i,j) h2o2(1,1,1)=backg_h2o2(i,k,j) IF (COSSZA(1,1) > 0.0) THEN XNO3(1,1,1) = 0.0 ELSE ! -- Fraction of night ! fnight = 1.0 - TTDAY(i,j)/86400.0 ! The original xno3 values have been averaged over daytime ! as well => divide by fnight to get the appropriate night-time ! fraction from the monthly average ! fnight/=0.0 (for fnight=0: all cosszax (including current ! cossza) > 0.0) xno3(1,1,1) = backg_no3(i,k,j) / (1.0 - TTDAY(i,j)/86400.) END IF ! if(i.eq.19.and.j.eq.19.and.k.eq.kts)then ! write(0,*)backg_oh(i,k,j),backg_no3(i,k,j),ttday(i,j),tcosz(i,j) ! endif call chmdrv_su( imx,jmx,lmx,& nmx, dt, tmp, airden, airmas, & oh, xno3, h2o2, cldf, tc, tdry,cossza, & chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, & chlso2_oh, chlso2_aq) chem(i,k,j,p_dms)=tc(1,1,1,1)*1.e6 chem(i,k,j,p_so2)=tc(1,1,1,2)*1.e6 chem(i,k,j,p_sulf)=tc(1,1,1,3)*1.e6 chem(i,k,j,p_msa)=tc(1,1,1,4)*1.e6 enddo enddo enddo CASE (GOCARTRACM_KPP,GOCARTRADM2) CALL wrf_debug(15,'calling gocart chemistry in addition to racm_kpp') do j=jts,jte do i=its,ite zenith=0. zenita=0. azimuth=0. rlat=xlat(i,j)*3.1415926535590/180. xlonn=xlong(i,j) CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat) cossza(1,1)=cosszax(1,1) do k=kts,kte-1 chldms_oh=0. chldms_no3=0. chldms_x=0. chpso2=0. chpmsa=0. chpso4=0. chlso2_oh=0. chlso2_aq=0. if( present(gd_cldf) ) then cldf(1,1,1)=gd_cldf(i,k,j) else cldf(1,1,1)=0. endif if(p_qi.gt.1)then if(moist(i,k,j,p_qc).gt.0.or.moist(i,k,j,p_qi).gt.0.)cldf(1,1,1)=1. elseif(p_qc.gt.1)then if(moist(i,k,j,p_qc).gt.0.)cldf(1,1,1)=1. endif tc(1,1,1,1)=chem(i,k,j,p_dms)*1.d-6 tc(1,1,1,2)=chem(i,k,j,p_so2)*1.d-6 tc(1,1,1,3)=chem(i,k,j,p_sulf)*1.d-6 tc(1,1,1,4)=chem(i,k,j,p_msa)*1.d-6 airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*dx*dx/g airden(1,1,1)=rho_phy(i,k,j) tmp(1,1,1)=t_phy(i,k,j) oh(1,1,1)=chem(i,k,j,p_ho)*1.d-6 h2o2(1,1,1)=chem(i,k,j,p_h2o2)*1.d-6 xno3(1,1,1) = chem(i,k,j,p_no3)*1.d-6 IF (COSSZA(1,1) > 0.0)xno3(1,1,1) = 0. ! if(i.eq.19.and.j.eq.19.and.k.eq.kts)then ! write(0,*)backg_oh(i,k,j),backg_no3(i,k,j),ttday(i,j),tcosz(i,j) ! endif call chmdrv_su( imx,jmx,lmx,& nmx, dt, tmp, airden, airmas, & oh, xno3, h2o2, cldf, tc, tdry,cossza, & chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, & chlso2_oh, chlso2_aq) chem(i,k,j,p_dms)=tc(1,1,1,1)*1.e6 chem(i,k,j,p_so2)=tc(1,1,1,2)*1.e6 chem(i,k,j,p_sulf)=tc(1,1,1,3)*1.e6 chem(i,k,j,p_msa)=tc(1,1,1,4)*1.e6 chem(i,k,j,p_h2o2)=h2o2(1,1,1)*1.e6 enddo enddo enddo END SELECT chem_select end subroutine gocart_chem_driver !SUBROUTINE chmdrv_su( & ! imx, jmx, lmx, nmx, ndt1, tmp, drydf, airden, airmas, & ! oh, xno3, h2o2, cldf, tc, tdry, depso2, depso4, depmsa, & ! chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, & ! chlso2_oh, chlso2_aq) !We don't apply losses due to dry deposition here, this is done in vertical mixing SUBROUTINE chmdrv_su( imx,jmx,lmx,& nmx, dt1, tmp, airden, airmas, & oh, xno3, h2o2, cldf, tc, tdry,cossza, & chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, & chlso2_oh, chlso2_aq) ! **************************************************************************** ! ** ** ! ** Chemistry subroutine. For tracers with dry deposition, the loss ** ! ** rate of dry dep is combined in chem loss term. ** ! ** ** ! **************************************************************************** ! USE module_data_gocart IMPLICIT NONE INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx integer :: ndt1 real, intent(in) :: dt1 REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: oh, xno3, cldf REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2 ! REAL*8, INTENT(IN) :: drydf(imx,jmx,nmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx) real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza ! REAL*8, DIMENSION(imx,jmx), INTENT(INOUT) :: depso2, depso4, depmsa REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_oh, chldms_no3 REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_x, chpso2, chpmsa REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chpso4, chlso2_oh, chlso2_aq REAL*8, DIMENSION(imx,jmx,lmx) :: pso2_dms, pmsa_dms, pso4_so2 ! executable statements ndt1=ifix(dt1) if(ndt1.le.0)stop CALL chem_dms(imx,jmx,lmx,nmx, ndt1, tmp, airden, airmas, oh, xno3, & tc, chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa,cossza, & pso2_dms, pmsa_dms) ! WRITE(*,*) 'after CHEM_DMS' CALL chem_so2(imx,jmx,lmx,nmx, ndt1, tmp, airden, airmas, & cldf, oh, h2o2, tc, tdry, cossza,& chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2) ! depso2, chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2) ! WRITE(*,*) 'after CHEM_SO2' CALL chem_so4(imx,jmx,lmx,nmx, ndt1, airmas, tc, tdry,cossza, & pso4_so2) ! depso4, pso4_so2) ! WRITE(*,*) 'after CHEM_SO4' CALL chem_msa(imx,jmx,lmx,nmx, ndt1, airmas, tc, tdry, cossza,& pmsa_dms) ! depmsa, pmsa_dms) ! WRITE(*,*) 'after CHEM_MSA' END SUBROUTINE chmdrv_su !============================================================================= SUBROUTINE chem_dms( imx,jmx,lmx,& nmx, ndt1, tmp, airden, airmas, oh, xno3, & tc, chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa,cossza, & pso2_dms, pmsa_dms) ! **************************************************************************** ! * * ! * This is DMS chemistry subroutine. * ! * * ! * R1: DMS + OH -> a*SO2 + b*MSA OH addition channel * ! * k1 = { 1.7e-42*exp(7810/T)*[O2] / (1+5.5e-31*exp(7460/T)*[O2] } * ! * a = 0.75, b = 0.25 * ! * * ! * R2: DMS + OH -> SO2 + ... OH abstraction channel * ! * k2 = 1.2e-11*exp(-260/T) * ! * * ! * DMS_OH = DMS0 * exp(-(r1+r2)*NDT1) * ! * where DMS0 is the DMS concentration at the beginning, * ! * r1 = k1*[OH], r2 = k2*[OH]. * ! * * ! * R3: DMS + NO3 -> SO2 + ... * ! * k3 = 1.9e-13*exp(500/T) * ! * * ! * DMS = DMS_OH * exp(-r3*NDT1) * ! * where r3 = k3*[NO3]. * ! * * ! * R4: DMS + X -> SO2 + ... * ! * assume to be at the rate of DMS+OH and DMS+NO3 combined. * ! * * ! * The production of SO2 and MSA here, PSO2_DMS and PMSA_DMS, are saved * ! * for use in CHEM_SO2 and CHEM_MSA subroutines as a source term. They * ! * are in unit of MixingRatio/timestep. * ! * * ! **************************************************************************** USE module_data_gocartchem IMPLICIT NONE INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: oh, xno3 REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_oh, chldms_no3 REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_x, chpso2, chpmsa REAL*8, DIMENSION(imx,jmx,lmx), INTENT(OUT) :: pso2_dms, pmsa_dms real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza REAL*8, PARAMETER :: fx = 1.0 REAL*8, PARAMETER :: a = 0.75 REAL*8, PARAMETER :: b = 0.25 ! From D4: only 0.8 efficiency, also some goes to DMSO and lost. ! So we assume 0.75 efficiency for DMS addtion channel to form ! products. REAL*8, PARAMETER :: eff = 1.0 ! -- Factor to convert AIRDEN from kgair/m3 to molecules/cm3: REAL*8, PARAMETER :: f = 1000.0 / airmw * 6.022D23 * 1.0D-6 INTEGER :: i, j, l REAL(KIND=8) :: tk, o2, dms0, rk1, rk2, rk3, dms_oh, dms, xoh, xn3, xx ! executable statements DO l = 1,lmx !CMIC$ doall autoscope DO j = 1,jmx DO i = 1,imx tk = tmp(i,j,l) o2 = airden(i,j,l) * f * 0.21 dms0 = tc(i,j,l,NDMS) ! **************************************************************************** ! * (1) DMS + OH: RK1 - addition channel; RK2 - abstraction channel. * ! **************************************************************************** rk1 = 0.0d0 rk2 = 0.0d0 rk3 = 0.0d0 IF (oh(i,j,l) > 0.0) THEN ! IF (TRIM(oh_units) == 'mol/mol') THEN ! mozech: oh is in mol/mol ! convert to molecules/cm3 rk1 = (1.7D-42 * EXP(7810.0/tk) * o2) / & (1.0 + 5.5D-31 * EXP(7460.0/tk) * o2 ) * oh(i,j,l) * & airden(i,j,l)*f rk2 = 1.2D-11*EXP(-260.0/tk) * oh(i,j,l)*airden(i,j,l)*f ! ELSE ! rk1 = (1.7D-42 * EXP(7810.0/tk) * o2) / & ! (1.0 + 5.5D-31 * EXP(7460.0/tk) * o2 ) * oh(i,j,l) ! rk2 = 1.2D-11*EXP(-260.0/tk) * oh(i,j,l) ! END IF END IF ! **************************************************************************** ! * (2) DMS + NO3 (only happens at night): * ! **************************************************************************** IF (cossza(i,j) <= 0.0) THEN ! IF (TRIM(no3_units) == 'cm-3') THEN ! ! IMAGES: XNO3 is in molecules/cm3. ! rk3 = 1.9D-13 * EXP(500.0/tk) * xno3(i,j,l) ! ELSE ! GEOSCHEM (mergechem) and mozech: XNO3 is in mol/mol (v/v) ! convert xno3 from volume mixing ratio to molecules/cm3 rk3 = 1.9D-13 * EXP(500.0/tk) * xno3(i,j,l) * & airden(i,j,l) * f ! END IF END IF ! **************************************************************************** ! * Update DMS concentrations after reaction with OH and NO3, and also * ! * account for DMS + X assuming at a rate as (DMS+OH)*Fx in the day and * ! * (DMS+NO3)*Fx at night: * ! * DMS_OH : DMS concentration after reaction with OH * ! * DMS : DMS concentration after reaction with NO3 * ! * (min(DMS) = 1.0E-32) * ! **************************************************************************** dms_oh = dms0 * EXP( -(rk1 + rk2) * fx * REAL(ndt1) ) dms = dms_oh * EXP( -(rk3) * fx * REAL(ndt1) ) dms = MAX(dms, 1.0D-32) tc(i,j,l,NDMS) = dms ! **************************************************************************** ! * Save SO2 and MSA production from DMS oxidation * ! * (in MixingRatio/timestep): * ! * * ! * SO2 is formed in DMS + OH addition (0.85) and abstraction (1.0) * ! * channels as well as DMS + NO3 reaction. We also assume that * ! * SO2 yield from DMS + X is 1.0. * ! * MSA is formed in DMS + OH addition (0.15) channel. * ! **************************************************************************** IF ((rk1 + rk2) == 0.0) THEN pmsa_dms(i,j,l) = 0.0D0 ELSE ! pmsa_dms(i,j,l) = (dms0 - dms_oh) * b*rk1/((rk1+rk2)*fx) pmsa_dms(i,j,l) = max(0.0D0,(dms0 - dms_oh) * b*rk1/((rk1+rk2) * fx) * eff) END IF pso2_dms(i,j,l) = max(0.0D0,dms0 - dms - pmsa_dms(i,j,l)) ! pso2_dms(i,j,l) = (dms0 - dms - pmsa_dms(i,j,l)/eff) * eff ! ------------------------------------------------------------ ! DIAGNOSTICS: DMS loss (kgS/timstep) ! SO2 production (kgS/timestep) ! MSA production (kgS/timestep) ! ------------------------------------------------------------ xoh = (dms0 - dms_oh) / fx * airmas(i,j,l)/airmw*smw xn3 = (dms_oh - dms) / fx * airmas(i,j,l)/airmw*smw xx = (dms0 - dms) * airmas(i,j,l)/airmw*smw - xoh - xn3 chldms_oh (i,j,l) = chldms_oh (i,j,l) + xoh chldms_no3(i,j,l) = chldms_no3(i,j,l) + xn3 chldms_x (i,j,l) = chldms_x (i,j,l) + xx chpso2(i,j,l) = chpso2(i,j,l) + pso2_dms(i,j,l) & * airmas(i,j,l) / airmw * smw chpmsa(i,j,l) = chpmsa(i,j,l) + pmsa_dms(i,j,l) & * airmas(i,j,l) / airmw * smw END DO END DO END DO END SUBROUTINE chem_dms !============================================================================= SUBROUTINE chem_so2( imx,jmx,lmx,& nmx, ndt1, tmp, airden, airmas, & cldf, oh, h2o2, tc, tdry, cossza,& chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2) ! depso2, chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2) ! **************************************************************************** ! * * ! * This is SO2 chemistry subroutine. * ! * * ! * SO2 production: * ! * DMS + OH, DMS + NO3 (saved in CHEM_DMS) * ! * * ! * SO2 loss: * ! * SO2 + OH -> SO4 * ! * SO2 -> drydep (NOT USED IN WRF/CHEM * ! * SO2 + H2O2 or O3 (aq) -> SO4 * ! * * ! * SO2 = SO2_0 * exp(-bt) * ! * + PSO2_DMS/bt * [1-exp(-bt)] * ! * where b is the sum of the reaction rate of SO2 + OH and the dry * ! * deposition rate of SO2, PSO2_DMS is SO2 production from DMS in * ! * MixingRatio/timestep. * ! * * ! * If there is cloud in the gridbox (fraction = fc), then the aqueous * ! * phase chemistry also takes place in cloud. The amount of SO2 oxidized * ! * by H2O2 in cloud is limited by the available H2O2; the rest may be * ! * oxidized due to additional chemistry, e.g, reaction with O3 or O2 * ! * (catalyzed by trace metal). * ! * * ! **************************************************************************** USE module_data_gocartchem IMPLICIT NONE INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: cldf, oh REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2 real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza ! REAL*8, INTENT(IN) :: drydf(imx,jmx,nmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx) ! REAL*8, DIMENSION(imx,jmx), INTENT(INOUT) :: depso2 REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chpso4, chlso2_oh, chlso2_aq REAL*8, INTENT(IN) :: pso2_dms(imx,jmx,lmx) REAL*8, INTENT(OUT) :: pso4_so2(imx,jmx,lmx) REAL*8 :: k0, kk, m, l1, l2, ld ! Factor to convert AIRDEN from kgair/m3 to molecules/cm3: REAL*8, PARAMETER :: f = 1000. / airmw * 6.022D23 * 1.0D-6 REAL*8, PARAMETER :: ki = 1.5D-12 INTEGER :: i, j, l REAL*8 :: so20, tk, f1, rk1, rk2, rk, rkt, so2_cd, fc, so2 ! executable statements DO l = 1,lmx DO j = 1,jmx DO i = 1,imx so20 = tc(i,j,l,NSO2) ! RK1: SO2 + OH(g), in s-1 tk = tmp(i,j,l) k0 = 3.0D-31 * (300.0/tk)**3.3 m = airden(i,j,l) * f kk = k0 * m / ki f1 = ( 1.0+ ( LOG10(kk) )**2 )**(-1) ! IF (TRIM(oh_units) == 'mol/mol') THEN ! mozech: oh is in mol/mol ! convert to molecules/cm3 rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * & oh(i,j,l)*airden(i,j,l)*f ! ELSE ! rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * oh(i,j,l) ! END IF ! RK2: SO2 drydep frequency, s-1 ! IF (l == 1) THEN ! at the surface ! rk2 = drydf(i,j,NSO2) ! ELSE rk2 = 0.0 ! END IF rk = (rk1 + rk2) rkt = rk * REAL(ndt1) ! **************************************************************************** ! * Update SO2 concentration after gas phase chemistry and deposition. * ! **************************************************************************** IF (rk > 0.0) THEN so2_cd = so20 * EXP(-rkt) & + pso2_dms(i,j,l) * (1.0 - EXP(-rkt)) / rkt l1 = (so20 - so2_cd + pso2_dms(i,j,l)) * rk1/rk IF (l == 1) THEN ld = (so20 - so2_cd + pso2_dms(i,j,l)) * rk2/rk ELSE ld = 0.0 END IF ELSE so2_cd = so20 l1 = 0.0 END IF ! **************************************************************************** ! * Update SO2 concentration after cloud chemistry. * ! * SO2 chemical loss rate = SO4 production rate (MixingRatio/timestep). * ! **************************************************************************** ! Cloud chemistry (above 258K): fc = cldf(i,j,l) IF (fc > 0.0 .AND. so2_cd > 0.0 .AND. tk > 258.0) THEN IF (so2_cd > h2o2(i,j,l)) THEN fc = fc * (h2o2(i,j,l)/so2_cd) h2o2(i,j,l) = h2o2(i,j,l) * (1.0 - cldf(i,j,l)) ELSE h2o2(i,j,l) = h2o2(i,j,l) * & (1.0 - cldf(i,j,l)*so2_cd/h2o2(i,j,l)) END IF so2 = so2_cd * (1.0 - fc) ! Aqueous phase SO2 loss rate (MixingRatio/timestep): l2 = so2_cd * fc ELSE so2 = so2_cd l2 = 0.0 END IF so2 = MAX(so2, 1.0D-32) tc(i,j,l,NSO2) = so2 ! **************************************************************************** ! * SO2 chemical loss rate = SO4 production rate (MixingRatio/timestep). * ! **************************************************************************** pso4_so2(i,j,l) = max(0.0D0,l1 + l2) ! --------------------------------------------------------------- ! DIAGNOSTICS: SO2 gas-phase loss (kgS/timestep) ! SO2 aqueous-phase loss (kgS/timestep) ! SO2 dry deposition loss (kgS/timestep) ! SO4 production (kgS/timestep) ! --------------------------------------------------------------- chlso2_oh(i,j,l) = chlso2_oh(i,j,l) & + l1 * airmas(i,j,l) / airmw * smw chlso2_aq(i,j,l) = chlso2_aq(i,j,l) & + l2 * airmas(i,j,l) / airmw * smw IF (l == 1) & ! depso2(i,j) = depso2(i,j) + ld * airmas(i,j,l) / airmw * smw chpso4(i,j,l) = chpso4(i,j,l) + pso4_so2(i,j,l) & * airmas(i,j,l) / airmw * smw END DO END DO END DO ! tdry(:,:,NSO2) = depso2(:,:)*tcmw(NSO2)/smw ! kg of SO2 END SUBROUTINE chem_so2 !============================================================================= SUBROUTINE chem_so4( imx,jmx,lmx,& nmx, ndt1, airmas, tc, tdry, cossza,& pso4_so2) ! depso4, pso4_so2) ! **************************************************************************** ! * * ! * This is SO4 chemistry subroutine. * ! * * ! * The Only production is from SO2 oxidation (save in CHEM_SO2), and the * ! * only loss is dry depsition here. Wet deposition will be treated in * ! * WETDEP subroutine. * ! * * ! * SO4 = SO4_0 * exp(-kt) + PSO4_SO2/kt * (1.-exp(-kt)) * ! * where k = dry deposition. * ! * * ! **************************************************************************** USE module_data_gocartchem IMPLICIT NONE INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: airmas ! REAL*8, INTENT(IN) :: drydf(imx,jmx,nmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx) ! REAL*8, DIMENSION(imx,jmx), INTENT(INOUT) :: depso4 REAL*8, INTENT(IN) :: pso4_so2(imx,jmx,lmx) real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza INTEGER :: i, j, l REAL*8 :: so40, rk, rkt, so4 ! executable statements DO l = 1,lmx DO j = 1,jmx DO i = 1,imx so40 = tc(i,j,l,NSO4) ! RK: SO4 drydep frequency, s-1 ! IF (l == 1) THEN ! rk = drydf(i,j,NSO4) ! rkt = rk * REAL(ndt1) ! ! so4 = so40 * EXP(-rkt) + pso4_so2(i,j,l)/rkt * (1.0 - EXP(-rkt)) ! ELSE so4 = so40 + pso4_so2(i,j,l) ! END IF if(pso4_so2(i,j,l).lt.0.)then write(0,*)'so4 routine, pso4 = ',pso4_so2(i,j,l),so4,so40 endif so4 = MAX(so4, 1.0D-32) tc(i,j,l,NSO4) = so4 ! -------------------------------------------------------------- ! DIAGNOSTICS: SO4 dry deposition (kgS/timestep) ! -------------------------------------------------------------- ! IF (l == 1) & ! depso4(i,j) = depso4(i,j) + (so40 - so4 + pso4_so2(i,j,l)) & ! * airmas(i,j,l) / airmw * smw END DO END DO END DO ! tdry(:,:,NSO4) = depso4(:,:)*tcmw(NSO4)/smw ! kg of SO4 END SUBROUTINE chem_so4 !============================================================================= SUBROUTINE chem_msa( imx,jmx,lmx,& nmx, ndt1, airmas, tc, tdry, cossza,& pmsa_dms) ! depmsa, pmsa_dms) ! **************************************************************************** ! * * ! * This is MSA chemistry subroutine. * ! * * ! * The Only production is from DMS oxidation (save in CHEM_DMS), and the * ! * only loss is dry depsition here. Wet deposition will be treated in * ! * WETDEP subroutine. * ! * * ! * MSA = MSA_0 * exp(-dt) + PMSA_DMS/kt * (1.-exp(-kt)) * ! * where k = dry deposition. * ! * * ! **************************************************************************** USE module_data_gocartchem IMPLICIT NONE INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: airmas REAL*8, DIMENSION(imx,jmx), INTENT(IN) :: cossza ! REAL, INTENT(IN) :: drydf(imx,jmx,nmx) REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx) REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx) ! REAL, DIMENSION(imx,jmx), INTENT(INOUT) :: depmsa REAL*8, INTENT(IN) :: pmsa_dms(imx,jmx,lmx) REAL*8 :: msa0, msa, rk, rkt INTEGER :: i, j, l ! executable statements DO l = 1,lmx DO j = 1,jmx DO i = 1,imx msa0 = tc(i,j,l,NMSA) ! RK: MSA drydep frequency, s-1 ! IF (l == 1) THEN ! rk = drydf(i,j,NMSA) ! rkt = rk * REAL(ndt1) ! ! msa = msa0 * EXP(-rkt) & ! + pmsa_dms(i,j,l)/rkt * (1.0 - EXP(-rkt)) ! ! ELSE msa = msa0 + pmsa_dms(i,j,l) ! END IF msa = MAX(msa, 1.0D-32) tc(i,j,l,NMSA) = msa ! -------------------------------------------------------------- ! DIAGNOSTICS: MSA dry deposition (kgS/timestep) ! -------------------------------------------------------------- ! IF (l == 1) & ! depmsa(i,j) = depmsa(i,j) + (msa0 - msa + pmsa_dms(i,j,l)) & ! * airmas(i,j,l) / airmw * smw END DO END DO END DO ! tdry(:,:,NMSA) = depmsa(:,:)*tcmw(NMSA)/smw ! kg of MSA END SUBROUTINE chem_msa SUBROUTINE szangle(imx, jmx, doy, xhour, sza, cossza,xlon,rlat) ! ! **************************************************************************** ! ** ** ! ** This subroutine computes solar zenith angle (SZA): ** ! ** ** ! ** cos(SZA) = sin(LAT)*sin(DEC) + cos(LAT)*cos(DEC)*cos(AHR) ** ! ** ** ! ** where LAT is the latitude angle, DEC is the solar declination angle, ** ! ** and AHR is the hour angle, all in radius. ** ! ** ** ! ** DOY = day-of-year, XHOUR = UT time (hrs). ** ! ** XLON = longitude in degree, RLAT = latitude in radian. ** ! **************************************************************************** ! IMPLICIT NONE INTEGER, INTENT(IN) :: imx, jmx INTEGER, INTENT(IN) :: doy REAL, INTENT(IN) :: xhour REAL, INTENT(OUT) :: sza(imx,jmx), cossza(imx,jmx) REAL :: a0, a1, a2, a3, b1, b2, b3, r, dec, timloc, ahr,xlon,rlat real, parameter :: pi=3.14 INTEGER :: i, j ! executable statements ! *************************************************************************** ! * Solar declination angle: * ! *************************************************************************** a0 = 0.006918 a1 = 0.399912 a2 = 0.006758 a3 = 0.002697 b1 = 0.070257 b2 = 0.000907 b3 = 0.000148 r = 2.0* pi * REAL(doy-1)/365.0 ! dec = a0 - a1*COS( r) + b1*SIN( r) & - a2*COS(2.0*r) + b2*SIN(2.0*r) & - a3*COS(3.0*r) + b3*SIN(3.0*r) ! DO i = 1,imx ! ************************************************************************ ! * Hour angle (AHR) is a function of longitude. AHR is zero at * ! * solar noon, and increases by 15 deg for every hour before or * ! * after solar noon. * ! ************************************************************************ ! -- Local time in hours timloc = xhour + xlon/15.0 ! IF (timloc < 0.0) timloc = 24.0 + timloc IF (timloc > 24.0) timloc = timloc - 24.0 ! ! -- Hour angle ahr = ABS(timloc - 12.0) * 15.0 * pi/180.0 ! DO j = 1,jmx ! -- Solar zenith angle cossza(i,j) = SIN(rlat) * SIN(dec) + & COS(rlat) * COS(dec) * COS(ahr) sza(i,j) = ACOS(cossza(i,j)) * 180.0/pi IF (cossza(i,j) < 0.0) cossza(i,j) = 0.0 ! END do END DO END subroutine szangle END MODULE MODULE_GOCART_CHEM