chem/module_gocart_chem.F

   1 MODULE MODULE_GOCART_CHEM
   2
   3 CONTAINS
   4
   5   subroutine gocart_chem_driver(curr_secs,dt,config_flags,            &
   6          gmt,julday,t_phy,moist,                                      &
   7          chem,rho_phy,dz8w,p8w,backg_oh,backg_h2o2,backg_no3,         &
   8          gd_cldf,dx,g,xlat,xlong,ttday,tcosz, &
   9          ids,ide, jds,jde, kds,kde,                                        &
  10          ims,ime, jms,jme, kms,kme,                                        &
  11          its,ite, jts,jte, kts,kte                                         )
  12   USE module_configure
  13   USE module_state_description
  14   USE module_phot_mad, only : calc_zenith
  15   IMPLICIT NONE
  16    TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  17
  18    INTEGER,      INTENT(IN   ) :: julday,                                  &
  19                                   ids,ide, jds,jde, kds,kde,               &
  20                                   ims,ime, jms,jme, kms,kme,               &
  21                                   its,ite, jts,jte, kts,kte
  22     REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
  23          INTENT(IN ) ::                                   moist
  24    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
  25          INTENT(INOUT ) ::                                   chem
  26    REAL,  DIMENSION( ims:ime , jms:jme ),                        &
  27           INTENT(IN   ) ::                                                 &
  28               xlat,xlong,ttday,tcosz
  29    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  30           OPTIONAL,                                                        &
  31           INTENT(IN   ) ::                     gd_cldf
  32    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  33           INTENT(IN   ) ::                     t_phy,                      &
  34                               backg_oh,backg_h2o2,backg_no3,dz8w,p8w,      &
  35                                               rho_phy
  36   REAL(KIND=8), INTENT(IN) :: curr_secs
  37   REAL, INTENT(IN   ) :: dt,dx,g,gmt
  38   integer :: nmx,i,j,k,imx,jmx,lmx
  39   real*8, DIMENSION (1,1,1) :: tmp,airden,airmas,oh,xno3,h2o2,chldms_oh,    &
  40                                chldms_no3,chldms_x,chpso2,chpmsa,chpso4,    &
  41                                chlso2_oh,chlso2_aq,cldf
  42   real*8, DIMENSION (1,1,4) :: tdry
  43   real*8, DIMENSION (1,1) :: cossza
  44   real, DIMENSION (1,1) :: sza,cosszax
  45   real*8, DIMENSION (1,1,1,4) :: tc,bems
  46   real*8, dimension (1) :: dxy
  47   real(kind=8) :: xtime,xhour
  48   real:: rlat,xlonn
  49   real :: zenith,zenita,azimuth,xmin,xtimin,gmtp
  50   integer(kind=8) :: ixhour
  51
  52   imx=1
  53   jmx=1
  54   lmx=1
  55   nmx=4
  56   tdry=0.d0
  57
  58   xtime=curr_secs/60._8
  59   ixhour=int(gmt+.01,8)+int(xtime/60._8,8)
  60   xhour=real(ixhour,8)
  61   xmin=60.*gmt+real(xtime-xhour*60._8,8)
  62   gmtp=mod(xhour,24._8)
  63   gmtp=gmtp+xmin/60.
  64
  65   dxy(1)=dx*dx
  66 !
  67 ! following arrays for busget stuff only
  68 !
  69 !
  70 !
  71        chem_select: SELECT CASE(config_flags%chem_opt)
  72           CASE (GOCART_SIMPLE)
  73            CALL wrf_debug(15,'calling gocart chemistry ')
  74        do j=jts,jte
  75        do i=its,ite
  76         zenith=0.
  77         zenita=0.
  78         azimuth=0.
  79         rlat=xlat(i,j)*3.1415926535590/180.
  80         xlonn=xlong(i,j)
  81         CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat)
  82         cossza(1,1)=cosszax(1,1)
  83 !
  84        do k=kts,kte-1
  85        chldms_oh=0.
  86        chldms_no3=0.
  87        chldms_x=0.
  88        chpso2=0.
  89        chpmsa=0.
  90        chpso4=0.
  91        chlso2_oh=0.
  92        chlso2_aq=0.
  93        if (present(gd_cldf) ) then
  94           cldf(1,1,1)=gd_cldf(i,k,j)
  95        else
  96           cldf(1,1,1)=0.
  97        endif
  98           if(p_qc.gt.1 .and. p_qi.gt.1)then
  99           if(moist(i,k,j,p_qc).gt.0.or.moist(i,k,j,p_qi).gt.0.)cldf(1,1,1)=1.
 100           elseif(p_qc.gt.1 .and. p_qi.le.1)then
 101           if(moist(i,k,j,p_qc).gt.0.)cldf(1,1,1)=1.
 102           endif
 103           tc(1,1,1,1)=chem(i,k,j,p_dms)*1.d-6
 104           tc(1,1,1,2)=chem(i,k,j,p_so2)*1.d-6
 105           tc(1,1,1,3)=chem(i,k,j,p_sulf)*1.d-6
 106           tc(1,1,1,4)=chem(i,k,j,p_msa)*1.d-6
 107           airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*dx*dx/g
 108           airden(1,1,1)=rho_phy(i,k,j)
 109           tmp(1,1,1)=t_phy(i,k,j)
 110           oh(1,1,1)=86400./dt*cossza(1,1)*backg_oh(i,k,j)/tcosz(i,j)
 111           h2o2(1,1,1)=backg_h2o2(i,k,j)
 112            IF (COSSZA(1,1) > 0.0) THEN
 113               XNO3(1,1,1) = 0.0
 114            ELSE
 115               ! -- Fraction of night
 116               ! fnight       = 1.0 - TTDAY(i,j)/86400.0
 117               ! The original xno3 values have been averaged over daytime
 118               ! as well => divide by fnight to get the appropriate night-time
 119               ! fraction from the monthly average
 120               ! fnight/=0.0 (for fnight=0: all cosszax (including current
 121               ! cossza) > 0.0)
 122               xno3(1,1,1) = backg_no3(i,k,j) / (1.0 - TTDAY(i,j)/86400.)
 123            END IF
 124 !         if(i.eq.19.and.j.eq.19.and.k.eq.kts)then
 125 !          write(0,*)backg_oh(i,k,j),backg_no3(i,k,j),ttday(i,j),tcosz(i,j)
 126 !         endif
 127
 128           call chmdrv_su( imx,jmx,lmx,&
 129                nmx, dt, tmp, airden, airmas, &
 130                oh, xno3, h2o2, cldf, tc, tdry,cossza,  &
 131                chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, &
 132                chlso2_oh, chlso2_aq)
 133           chem(i,k,j,p_dms)=tc(1,1,1,1)*1.e6
 134           chem(i,k,j,p_so2)=tc(1,1,1,2)*1.e6
 135           chem(i,k,j,p_sulf)=tc(1,1,1,3)*1.e6
 136           chem(i,k,j,p_msa)=tc(1,1,1,4)*1.e6
 137        enddo
 138        enddo
 139        enddo
 140      CASE (GOCARTRACM_KPP,GOCARTRADM2)
 141        CALL wrf_debug(15,'calling gocart chemistry in addition to racm_kpp')
 142        do j=jts,jte
 143        do i=its,ite
 144         zenith=0.
 145         zenita=0.
 146         azimuth=0.
 147         rlat=xlat(i,j)*3.1415926535590/180.
 148         xlonn=xlong(i,j)
 149         CALL szangle(1, 1, julday, gmtp, sza, cosszax,xlonn,rlat)
 150         cossza(1,1)=cosszax(1,1)
 151        do k=kts,kte-1
 152        chldms_oh=0.
 153        chldms_no3=0.
 154        chldms_x=0.
 155        chpso2=0.
 156        chpmsa=0.
 157        chpso4=0.
 158        chlso2_oh=0.
 159        chlso2_aq=0.
 160           if( present(gd_cldf) ) then
 161             cldf(1,1,1)=gd_cldf(i,k,j)
 162           else
 163             cldf(1,1,1)=0.
 164           endif
 165           if(p_qi.gt.1)then
 166           if(moist(i,k,j,p_qc).gt.0.or.moist(i,k,j,p_qi).gt.0.)cldf(1,1,1)=1.
 167           elseif(p_qc.gt.1)then
 168           if(moist(i,k,j,p_qc).gt.0.)cldf(1,1,1)=1.
 169           endif
 170           tc(1,1,1,1)=chem(i,k,j,p_dms)*1.d-6
 171           tc(1,1,1,2)=chem(i,k,j,p_so2)*1.d-6
 172           tc(1,1,1,3)=chem(i,k,j,p_sulf)*1.d-6
 173           tc(1,1,1,4)=chem(i,k,j,p_msa)*1.d-6
 174           airmas(1,1,1)=-(p8w(i,k+1,j)-p8w(i,k,j))*dx*dx/g
 175           airden(1,1,1)=rho_phy(i,k,j)
 176           tmp(1,1,1)=t_phy(i,k,j)
 177           oh(1,1,1)=chem(i,k,j,p_ho)*1.d-6
 178           h2o2(1,1,1)=chem(i,k,j,p_h2o2)*1.d-6
 179           xno3(1,1,1) = chem(i,k,j,p_no3)*1.d-6
 180           IF (COSSZA(1,1) > 0.0)xno3(1,1,1) = 0.
 181 !         if(i.eq.19.and.j.eq.19.and.k.eq.kts)then
 182 !          write(0,*)backg_oh(i,k,j),backg_no3(i,k,j),ttday(i,j),tcosz(i,j)
 183 !         endif
 184
 185           call chmdrv_su( imx,jmx,lmx,&
 186                nmx, dt, tmp, airden, airmas, &
 187                oh, xno3, h2o2, cldf, tc, tdry,cossza,  &
 188                chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, &
 189                chlso2_oh, chlso2_aq)
 190           chem(i,k,j,p_dms)=tc(1,1,1,1)*1.e6
 191           chem(i,k,j,p_so2)=tc(1,1,1,2)*1.e6
 192           chem(i,k,j,p_sulf)=tc(1,1,1,3)*1.e6
 193           chem(i,k,j,p_msa)=tc(1,1,1,4)*1.e6
 194           chem(i,k,j,p_h2o2)=h2o2(1,1,1)*1.e6
 195        enddo
 196        enddo
 197        enddo
 198    END SELECT chem_select
 199 end subroutine gocart_chem_driver
 200
 201 !SUBROUTINE chmdrv_su( &
 202 !     imx, jmx, lmx, nmx, ndt1, tmp, drydf, airden, airmas, &
 203 !     oh, xno3, h2o2, cldf, tc, tdry, depso2, depso4, depmsa, &
 204 !     chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, &
 205 !     chlso2_oh, chlso2_aq)
 206
 207 !We don't apply losses due to dry deposition here, this is done in vertical mixing
 208 SUBROUTINE chmdrv_su( imx,jmx,lmx,&
 209      nmx, dt1, tmp, airden, airmas, &
 210      oh, xno3, h2o2, cldf, tc, tdry,cossza,  &
 211      chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa, chpso4, &
 212      chlso2_oh, chlso2_aq)
 213
 214 ! ****************************************************************************
 215 ! **                                                                        **
 216 ! **  Chemistry subroutine.  For tracers with dry deposition, the loss      **
 217 ! **  rate of dry dep is combined in chem loss term.                        **
 218 ! **                                                                        **
 219 ! ****************************************************************************
 220
 221 ! USE module_data_gocart
 222
 223   IMPLICIT NONE
 224
 225   INTEGER, INTENT(IN) :: nmx,imx,jmx,lmx
 226   integer :: ndt1
 227   real, intent(in) :: dt1
 228   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas
 229   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: oh, xno3, cldf
 230   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2
 231 !  REAL*8, INTENT(IN)    :: drydf(imx,jmx,nmx)
 232   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 233   REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx)
 234   real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza
 235 !  REAL*8, DIMENSION(imx,jmx),     INTENT(INOUT) :: depso2, depso4, depmsa
 236   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_oh, chldms_no3
 237   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_x, chpso2, chpmsa
 238   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chpso4, chlso2_oh, chlso2_aq
 239
 240   REAL*8, DIMENSION(imx,jmx,lmx) :: pso2_dms, pmsa_dms, pso4_so2
 241
 242   ! executable statements
 243   ndt1=ifix(dt1)
 244   if(ndt1.le.0)stop
 245
 246      CALL chem_dms(imx,jmx,lmx,nmx, ndt1, tmp, airden, airmas, oh, xno3, &
 247           tc, chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa,cossza, &
 248           pso2_dms, pmsa_dms)
 249 !     WRITE(*,*) 'after CHEM_DMS'
 250      CALL chem_so2(imx,jmx,lmx,nmx, ndt1, tmp, airden, airmas, &
 251           cldf, oh, h2o2, tc, tdry, cossza,&
 252           chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2)
 253 !          depso2, chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2)
 254 !     WRITE(*,*) 'after CHEM_SO2'
 255      CALL chem_so4(imx,jmx,lmx,nmx, ndt1, airmas, tc, tdry,cossza, &
 256           pso4_so2)
 257 !          depso4, pso4_so2)
 258 !     WRITE(*,*) 'after CHEM_SO4'
 259      CALL chem_msa(imx,jmx,lmx,nmx, ndt1, airmas, tc, tdry, cossza,&
 260           pmsa_dms)
 261 !          depmsa, pmsa_dms)
 262 !     WRITE(*,*) 'after CHEM_MSA'
 263
 264 END SUBROUTINE chmdrv_su
 265
 266 !=============================================================================
 267 SUBROUTINE chem_dms( imx,jmx,lmx,&
 268      nmx, ndt1, tmp, airden, airmas, oh, xno3, &
 269      tc, chldms_oh, chldms_no3, chldms_x, chpso2, chpmsa,cossza, &
 270      pso2_dms, pmsa_dms)
 271
 272 ! ****************************************************************************
 273 ! *                                                                          *
 274 ! *  This is DMS chemistry subroutine.                                       *
 275 ! *                                                                          *
 276 ! *  R1:    DMS + OH  -> a*SO2 + b*MSA                OH addition channel    *
 277 ! *         k1 = { 1.7e-42*exp(7810/T)*[O2] / (1+5.5e-31*exp(7460/T)*[O2] }  *
 278 ! *         a = 0.75, b = 0.25                                               *
 279 ! *                                                                          *
 280 ! *  R2:    DMS + OH  ->   SO2 + ...                  OH abstraction channel *
 281 ! *         k2 = 1.2e-11*exp(-260/T)                                         *
 282 ! *                                                                          *
 283 ! *     DMS_OH = DMS0 * exp(-(r1+r2)*NDT1)                                   *
 284 ! *         where DMS0 is the DMS concentration at the beginning,            *
 285 ! *         r1 = k1*[OH], r2 = k2*[OH].                                      *
 286 ! *                                                                          *
 287 ! *  R3:    DMS + NO3 ->   SO2 + ...                                         *
 288 ! *         k3 = 1.9e-13*exp(500/T)                                          *
 289 ! *                                                                          *
 290 ! *     DMS = DMS_OH * exp(-r3*NDT1)                                         *
 291 ! *         where r3 = k3*[NO3].                                             *
 292 ! *                                                                          *
 293 ! *  R4:    DMS + X   ->   SO2 + ...                                         *
 294 ! *         assume to be at the rate of DMS+OH and DMS+NO3 combined.         *
 295 ! *                                                                          *
 296 ! *  The production of SO2 and MSA here, PSO2_DMS and PMSA_DMS, are saved    *
 297 ! *  for use in CHEM_SO2 and CHEM_MSA subroutines as a source term.  They    *
 298 ! *  are in unit of MixingRatio/timestep.                                    *
 299 ! *                                                                          *
 300 ! ****************************************************************************
 301
 302   USE module_data_gocartchem
 303
 304   IMPLICIT NONE
 305
 306   INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx
 307   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas
 308   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: oh, xno3
 309   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 310   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_oh, chldms_no3
 311   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chldms_x, chpso2, chpmsa
 312   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(OUT)   :: pso2_dms, pmsa_dms
 313   real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza
 314
 315   REAL*8, PARAMETER :: fx = 1.0
 316   REAL*8, PARAMETER :: a = 0.75
 317   REAL*8, PARAMETER :: b = 0.25
 318
 319   ! From D4: only 0.8 efficiency, also some goes to DMSO and lost.
 320   ! So we assume 0.75 efficiency for DMS addtion channel to form
 321   ! products.
 322
 323   REAL*8, PARAMETER :: eff = 1.0
 324   ! -- Factor to convert AIRDEN from kgair/m3 to molecules/cm3:
 325   REAL*8, PARAMETER :: f = 1000.0 / airmw * 6.022D23 * 1.0D-6
 326   INTEGER :: i, j, l
 327   REAL(KIND=8) :: tk, o2, dms0, rk1, rk2, rk3, dms_oh, dms, xoh, xn3, xx
 328
 329   ! executable statements
 330
 331   DO l = 1,lmx
 332 !CMIC$ doall autoscope
 333      DO j = 1,jmx
 334         DO i = 1,imx
 335
 336            tk = tmp(i,j,l)
 337            o2 = airden(i,j,l) * f * 0.21
 338            dms0 = tc(i,j,l,NDMS)
 339
 340 ! ****************************************************************************
 341 ! *  (1) DMS + OH:  RK1 - addition channel;  RK2 - abstraction channel.      *
 342 ! ****************************************************************************
 343
 344            rk1 = 0.0d0
 345            rk2 = 0.0d0
 346            rk3 = 0.0d0
 347
 348            IF (oh(i,j,l) > 0.0) THEN
 349 !             IF (TRIM(oh_units) == 'mol/mol') THEN
 350                  ! mozech: oh is in mol/mol
 351                  ! convert to molecules/cm3
 352                  rk1 = (1.7D-42 * EXP(7810.0/tk) * o2) / &
 353                       (1.0 + 5.5D-31 * EXP(7460.0/tk) * o2 ) * oh(i,j,l) * &
 354                       airden(i,j,l)*f
 355                  rk2 = 1.2D-11*EXP(-260.0/tk) * oh(i,j,l)*airden(i,j,l)*f
 356 !             ELSE
 357 !                rk1 = (1.7D-42 * EXP(7810.0/tk) * o2) / &
 358 !                     (1.0 + 5.5D-31 * EXP(7460.0/tk) * o2 ) * oh(i,j,l)
 359 !                rk2 = 1.2D-11*EXP(-260.0/tk) * oh(i,j,l)
 360 !             END IF
 361            END IF
 362
 363 ! ****************************************************************************
 364 ! *  (2) DMS + NO3 (only happens at night):                                  *
 365 ! ****************************************************************************
 366
 367            IF (cossza(i,j) <= 0.0) THEN
 368
 369 !             IF (TRIM(no3_units) == 'cm-3') THEN
 370 !                ! IMAGES: XNO3 is in molecules/cm3.
 371 !                rk3 = 1.9D-13 * EXP(500.0/tk) * xno3(i,j,l)
 372
 373 !             ELSE
 374                  ! GEOSCHEM (mergechem) and mozech: XNO3 is in mol/mol (v/v)
 375                  ! convert xno3 from volume mixing ratio to molecules/cm3
 376                  rk3 = 1.9D-13 * EXP(500.0/tk) * xno3(i,j,l) * &
 377                       airden(i,j,l) * f
 378 !             END IF
 379
 380            END IF
 381
 382 ! ****************************************************************************
 383 ! *  Update DMS concentrations after reaction with OH and NO3, and also      *
 384 ! *  account for DMS + X assuming at a rate as (DMS+OH)*Fx in the day and    *
 385 ! *  (DMS+NO3)*Fx at night:                                                  *
 386 ! *       DMS_OH       :  DMS concentration after reaction with OH           *
 387 ! *       DMS          :  DMS concentration after reaction with NO3          *
 388 ! *                           (min(DMS) = 1.0E-32)                           *
 389 ! ****************************************************************************
 390
 391            dms_oh = dms0   * EXP( -(rk1 + rk2) * fx * REAL(ndt1) )
 392            dms    = dms_oh * EXP( -(rk3) * fx * REAL(ndt1) )
 393            dms    = MAX(dms, 1.0D-32)
 394
 395            tc(i,j,l,NDMS) = dms
 396
 397 ! ****************************************************************************
 398 ! *  Save SO2 and MSA production from DMS oxidation                          *
 399 ! *  (in MixingRatio/timestep):                                              *
 400 ! *                                                                          *
 401 ! *  SO2 is formed in DMS + OH addition (0.85) and abstraction (1.0)         *
 402 ! *      channels as well as DMS + NO3 reaction.  We also assume that        *
 403 ! *      SO2 yield from DMS + X is 1.0.                                      *
 404 ! *  MSA is formed in DMS + OH addition (0.15) channel.                      *
 405 ! ****************************************************************************
 406
 407            IF ((rk1 + rk2) == 0.0) THEN
 408               pmsa_dms(i,j,l) = 0.0D0
 409            ELSE
 410 !       pmsa_dms(i,j,l) = (dms0 - dms_oh) * b*rk1/((rk1+rk2)*fx)
 411               pmsa_dms(i,j,l) = max(0.0D0,(dms0 - dms_oh) * b*rk1/((rk1+rk2) * fx) * eff)
 412            END IF
 413        pso2_dms(i,j,l) =  max(0.0D0,dms0 - dms - pmsa_dms(i,j,l))
 414 !      pso2_dms(i,j,l) =  (dms0 - dms - pmsa_dms(i,j,l)/eff) * eff
 415
 416            !    ------------------------------------------------------------
 417            !    DIAGNOSTICS:      DMS loss       (kgS/timstep)
 418            !                      SO2 production (kgS/timestep)
 419            !                      MSA production (kgS/timestep)
 420            !    ------------------------------------------------------------
 421            xoh  = (dms0   - dms_oh) / fx  * airmas(i,j,l)/airmw*smw
 422            xn3  = (dms_oh - dms)    / fx  * airmas(i,j,l)/airmw*smw
 423            xx   = (dms0 - dms) * airmas(i,j,l)/airmw*smw - xoh - xn3
 424
 425            chldms_oh (i,j,l) = chldms_oh (i,j,l) + xoh
 426            chldms_no3(i,j,l) = chldms_no3(i,j,l) + xn3
 427            chldms_x  (i,j,l) = chldms_x  (i,j,l) + xx
 428
 429            chpso2(i,j,l) = chpso2(i,j,l) + pso2_dms(i,j,l) &
 430                 * airmas(i,j,l) / airmw * smw
 431            chpmsa(i,j,l) = chpmsa(i,j,l) + pmsa_dms(i,j,l) &
 432                 * airmas(i,j,l) / airmw * smw
 433
 434         END DO
 435      END DO
 436   END DO
 437
 438 END SUBROUTINE chem_dms
 439
 440 !=============================================================================
 441
 442 SUBROUTINE chem_so2( imx,jmx,lmx,&
 443      nmx, ndt1, tmp, airden, airmas, &
 444      cldf, oh, h2o2, tc, tdry, cossza,&
 445      chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2)
 446 !     depso2, chpso4, chlso2_oh, chlso2_aq, pso2_dms, pso4_so2)
 447
 448 ! ****************************************************************************
 449 ! *                                                                          *
 450 ! *  This is SO2 chemistry subroutine.                                       *
 451 ! *                                                                          *
 452 ! *  SO2 production:                                                         *
 453 ! *    DMS + OH, DMS + NO3 (saved in CHEM_DMS)                               *
 454 ! *                                                                          *
 455 ! *  SO2 loss:                                                               *
 456 ! *    SO2 + OH  -> SO4                                                      *
 457 ! *    SO2       -> drydep (NOT USED IN WRF/CHEM                             *
 458 ! *    SO2 + H2O2 or O3 (aq) -> SO4                                          *
 459 ! *                                                                          *
 460 ! *  SO2 = SO2_0 * exp(-bt)                                                  *
 461 ! *      + PSO2_DMS/bt * [1-exp(-bt)]                                        *
 462 ! *    where b is the sum of the reaction rate of SO2 + OH and the dry       *
 463 ! *    deposition rate of SO2, PSO2_DMS is SO2 production from DMS in        *
 464 ! *    MixingRatio/timestep.                                                 *
 465 ! *                                                                          *
 466 ! *  If there is cloud in the gridbox (fraction = fc), then the aqueous      *
 467 ! *  phase chemistry also takes place in cloud. The amount of SO2 oxidized   *
 468 ! *  by H2O2 in cloud is limited by the available H2O2; the rest may be      *
 469 ! *  oxidized due to additional chemistry, e.g, reaction with O3 or O2       *
 470 ! *  (catalyzed by trace metal).                                             *
 471 ! *                                                                          *
 472 ! ****************************************************************************
 473   USE module_data_gocartchem
 474
 475   IMPLICIT NONE
 476
 477   INTEGER, INTENT(IN) ::  nmx, ndt1,imx,jmx,lmx
 478   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: tmp, airden, airmas
 479   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: cldf, oh
 480   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: h2o2
 481   real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza
 482 !  REAL*8, INTENT(IN)    :: drydf(imx,jmx,nmx)
 483   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 484   REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx)
 485
 486 !  REAL*8, DIMENSION(imx,jmx),     INTENT(INOUT) :: depso2
 487   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(INOUT) :: chpso4, chlso2_oh, chlso2_aq
 488   REAL*8, INTENT(IN)  :: pso2_dms(imx,jmx,lmx)
 489   REAL*8, INTENT(OUT) :: pso4_so2(imx,jmx,lmx)
 490
 491   REAL*8 ::  k0, kk, m, l1, l2, ld
 492   ! Factor to convert AIRDEN from kgair/m3 to molecules/cm3:
 493   REAL*8, PARAMETER :: f  = 1000. / airmw * 6.022D23 * 1.0D-6
 494   REAL*8, PARAMETER :: ki = 1.5D-12
 495   INTEGER :: i, j, l
 496   REAL*8 :: so20, tk, f1, rk1, rk2, rk, rkt, so2_cd, fc, so2
 497
 498   ! executable statements
 499
 500   DO l = 1,lmx
 501      DO j = 1,jmx
 502         DO i = 1,imx
 503
 504            so20 = tc(i,j,l,NSO2)
 505
 506            ! RK1: SO2 + OH(g), in s-1
 507            tk = tmp(i,j,l)
 508            k0 = 3.0D-31 * (300.0/tk)**3.3
 509            m  = airden(i,j,l) * f
 510            kk = k0 * m / ki
 511            f1 = ( 1.0+ ( LOG10(kk) )**2 )**(-1)
 512 !          IF (TRIM(oh_units) == 'mol/mol') THEN
 513               ! mozech: oh is in mol/mol
 514               ! convert to molecules/cm3
 515               rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * &
 516                    oh(i,j,l)*airden(i,j,l)*f
 517 !          ELSE
 518 !             rk1 = ( k0 * m / (1.0 + kk) ) * 0.6**f1 * oh(i,j,l)
 519 !          END IF
 520
 521            ! RK2: SO2 drydep frequency, s-1
 522 !           IF (l == 1) THEN ! at the surface
 523 !              rk2 = drydf(i,j,NSO2)
 524 !           ELSE
 525               rk2 = 0.0
 526 !           END IF
 527
 528            rk  = (rk1 + rk2)
 529            rkt =  rk * REAL(ndt1)
 530
 531 ! ****************************************************************************
 532 ! *  Update SO2 concentration after gas phase chemistry and deposition.      *
 533 ! ****************************************************************************
 534
 535            IF (rk > 0.0) THEN
 536               so2_cd = so20 * EXP(-rkt) &
 537                    + pso2_dms(i,j,l) * (1.0 - EXP(-rkt)) / rkt
 538               l1     = (so20 - so2_cd + pso2_dms(i,j,l)) * rk1/rk
 539               IF (l == 1) THEN
 540                  ld    = (so20 - so2_cd + pso2_dms(i,j,l)) * rk2/rk
 541               ELSE
 542                  ld    = 0.0
 543               END IF
 544            ELSE
 545               so2_cd = so20
 546               l1 = 0.0
 547            END IF
 548
 549 ! ****************************************************************************
 550 ! *  Update SO2 concentration after cloud chemistry.                         *
 551 ! *  SO2 chemical loss rate  = SO4 production rate (MixingRatio/timestep).   *
 552 ! ****************************************************************************
 553
 554            ! Cloud chemistry (above 258K):
 555            fc = cldf(i,j,l)
 556            IF (fc > 0.0 .AND. so2_cd > 0.0 .AND. tk > 258.0) THEN
 557
 558               IF (so2_cd > h2o2(i,j,l)) THEN
 559                  fc = fc * (h2o2(i,j,l)/so2_cd)
 560                  h2o2(i,j,l) = h2o2(i,j,l) * (1.0 - cldf(i,j,l))
 561               ELSE
 562                  h2o2(i,j,l) = h2o2(i,j,l) * &
 563                       (1.0 - cldf(i,j,l)*so2_cd/h2o2(i,j,l))
 564               END IF
 565               so2 = so2_cd * (1.0 - fc)
 566               ! Aqueous phase SO2 loss rate (MixingRatio/timestep):
 567               l2  = so2_cd * fc
 568            ELSE
 569               so2 = so2_cd
 570               l2 = 0.0
 571            END IF
 572
 573            so2    = MAX(so2, 1.0D-32)
 574            tc(i,j,l,NSO2) = so2
 575
 576 ! ****************************************************************************
 577 ! *  SO2 chemical loss rate  = SO4 production rate (MixingRatio/timestep).   *
 578 ! ****************************************************************************
 579
 580            pso4_so2(i,j,l) = max(0.0D0,l1 + l2)
 581
 582            !    ---------------------------------------------------------------
 583            !    DIAGNOSTICS:      SO2 gas-phase loss       (kgS/timestep)
 584            !                      SO2 aqueous-phase loss   (kgS/timestep)
 585            !                      SO2 dry deposition loss  (kgS/timestep)
 586            !                      SO4 production           (kgS/timestep)
 587            !    ---------------------------------------------------------------
 588            chlso2_oh(i,j,l) = chlso2_oh(i,j,l) &
 589                 + l1 * airmas(i,j,l) / airmw * smw
 590            chlso2_aq(i,j,l) = chlso2_aq(i,j,l) &
 591                 + l2 * airmas(i,j,l) / airmw * smw
 592            IF (l == 1) &
 593 !                depso2(i,j) = depso2(i,j) + ld * airmas(i,j,l) / airmw * smw
 594
 595            chpso4(i,j,l) = chpso4(i,j,l) + pso4_so2(i,j,l) &
 596                 * airmas(i,j,l) / airmw * smw
 597
 598         END DO
 599      END DO
 600   END DO
 601
 602 !  tdry(:,:,NSO2) = depso2(:,:)*tcmw(NSO2)/smw ! kg of SO2
 603
 604 END SUBROUTINE chem_so2
 605
 606 !=============================================================================
 607
 608 SUBROUTINE chem_so4( imx,jmx,lmx,&
 609      nmx, ndt1, airmas, tc, tdry, cossza,&
 610      pso4_so2)
 611 !     depso4, pso4_so2)
 612
 613 ! ****************************************************************************
 614 ! *                                                                          *
 615 ! *  This is SO4 chemistry subroutine.                                       *
 616 ! *                                                                          *
 617 ! *  The Only production is from SO2 oxidation (save in CHEM_SO2), and the   *
 618 ! *  only loss is dry depsition here.  Wet deposition will be treated in     *
 619 ! *  WETDEP subroutine.                                                      *
 620 ! *                                                                          *
 621 ! *  SO4 = SO4_0 * exp(-kt) + PSO4_SO2/kt * (1.-exp(-kt))                    *
 622 ! *    where k = dry deposition.                                             *
 623 ! *                                                                          *
 624 ! ****************************************************************************
 625   USE module_data_gocartchem
 626
 627   IMPLICIT NONE
 628
 629   INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx
 630   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: airmas
 631 !  REAL*8, INTENT(IN)    :: drydf(imx,jmx,nmx)
 632   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 633   REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx)
 634
 635 !  REAL*8, DIMENSION(imx,jmx), INTENT(INOUT) :: depso4
 636   REAL*8, INTENT(IN) :: pso4_so2(imx,jmx,lmx)
 637   real*8, DIMENSION (imx,jmx),INTENT(IN) :: cossza
 638
 639   INTEGER :: i, j, l
 640   REAL*8 :: so40, rk, rkt, so4
 641
 642   ! executable statements
 643
 644   DO l = 1,lmx
 645      DO j = 1,jmx
 646         DO i = 1,imx
 647
 648            so40 = tc(i,j,l,NSO4)
 649
 650            ! RK: SO4 drydep frequency, s-1
 651 !           IF (l == 1) THEN
 652 !              rk  = drydf(i,j,NSO4)
 653 !              rkt = rk * REAL(ndt1)
 654 !
 655 !              so4 = so40 * EXP(-rkt) + pso4_so2(i,j,l)/rkt * (1.0 - EXP(-rkt))
 656 !           ELSE
 657               so4 = so40 + pso4_so2(i,j,l)
 658 !           END IF
 659            if(pso4_so2(i,j,l).lt.0.)then
 660              write(0,*)'so4 routine, pso4 = ',pso4_so2(i,j,l),so4,so40
 661            endif
 662
 663            so4    = MAX(so4, 1.0D-32)
 664            tc(i,j,l,NSO4) = so4
 665
 666            !  --------------------------------------------------------------
 667            !  DIAGNOSTICS:      SO4 dry deposition  (kgS/timestep)
 668            !  --------------------------------------------------------------
 669 !           IF (l == 1) &
 670 !                depso4(i,j) = depso4(i,j) + (so40 - so4 + pso4_so2(i,j,l)) &
 671 !                * airmas(i,j,l) / airmw * smw
 672
 673         END DO
 674      END DO
 675   END DO
 676
 677  ! tdry(:,:,NSO4) = depso4(:,:)*tcmw(NSO4)/smw ! kg of SO4
 678
 679 END SUBROUTINE chem_so4
 680
 681 !=============================================================================
 682
 683 SUBROUTINE chem_msa( imx,jmx,lmx,&
 684      nmx, ndt1, airmas, tc, tdry, cossza,&
 685      pmsa_dms)
 686 !     depmsa, pmsa_dms)
 687
 688 ! ****************************************************************************
 689 ! *                                                                          *
 690 ! *  This is MSA chemistry subroutine.                                       *
 691 ! *                                                                          *
 692 ! *  The Only production is from DMS oxidation (save in CHEM_DMS), and the   *
 693 ! *  only loss is dry depsition here.  Wet deposition will be treated in     *
 694 ! *  WETDEP subroutine.                                                      *
 695 ! *                                                                          *
 696 ! *  MSA = MSA_0 * exp(-dt) + PMSA_DMS/kt * (1.-exp(-kt))                    *
 697 ! *    where k = dry deposition.                                             *
 698 ! *                                                                          *
 699 ! ****************************************************************************
 700   USE module_data_gocartchem
 701
 702   IMPLICIT NONE
 703
 704   INTEGER, INTENT(IN) :: nmx, ndt1,imx,jmx,lmx
 705   REAL*8, DIMENSION(imx,jmx,lmx), INTENT(IN) :: airmas
 706   REAL*8, DIMENSION(imx,jmx), INTENT(IN) :: cossza
 707 !  REAL, INTENT(IN)    :: drydf(imx,jmx,nmx)
 708   REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
 709   REAL*8, INTENT(INOUT) :: tdry(imx,jmx,nmx)
 710 !  REAL, DIMENSION(imx,jmx), INTENT(INOUT) :: depmsa
 711   REAL*8, INTENT(IN) :: pmsa_dms(imx,jmx,lmx)
 712
 713   REAL*8 :: msa0, msa, rk, rkt
 714   INTEGER :: i, j, l
 715
 716   ! executable statements
 717
 718   DO l = 1,lmx
 719      DO j = 1,jmx
 720         DO i = 1,imx
 721
 722            msa0 = tc(i,j,l,NMSA)
 723
 724            ! RK: MSA drydep frequency, s-1
 725 !           IF (l == 1) THEN
 726 !              rk  = drydf(i,j,NMSA)
 727 !              rkt = rk * REAL(ndt1)
 728 !
 729 !              msa = msa0 * EXP(-rkt) &
 730 !                   + pmsa_dms(i,j,l)/rkt * (1.0 - EXP(-rkt))
 731 !
 732 !           ELSE
 733               msa = msa0 + pmsa_dms(i,j,l)
 734 !           END IF
 735
 736            msa    = MAX(msa, 1.0D-32)
 737            tc(i,j,l,NMSA) = msa
 738
 739            !  --------------------------------------------------------------
 740            !  DIAGNOSTICS:      MSA dry deposition  (kgS/timestep)
 741            !  --------------------------------------------------------------
 742 !           IF (l == 1) &
 743 !                depmsa(i,j) = depmsa(i,j) + (msa0 - msa + pmsa_dms(i,j,l)) &
 744 !                * airmas(i,j,l) / airmw * smw
 745
 746         END DO
 747      END DO
 748   END DO
 749
 750 !  tdry(:,:,NMSA) = depmsa(:,:)*tcmw(NMSA)/smw ! kg of MSA
 751
 752 END SUBROUTINE chem_msa
 753 SUBROUTINE szangle(imx, jmx, doy, xhour, sza, cossza,xlon,rlat)
 754
 755 !
 756 ! ****************************************************************************
 757 ! **                                                                        **
 758 ! **  This subroutine computes solar zenith angle (SZA):                    **
 759 ! **                                                                        **
 760 ! **      cos(SZA) = sin(LAT)*sin(DEC) + cos(LAT)*cos(DEC)*cos(AHR)         **
 761 ! **                                                                        **
 762 ! **  where LAT is the latitude angle, DEC is the solar declination angle,  **
 763 ! **  and AHR is the hour angle, all in radius.                             **
 764 ! **                                                                        **
 765 ! **  DOY = day-of-year, XHOUR = UT time (hrs).                             **
 766 ! **  XLON = longitude in degree, RLAT = latitude in radian.                **
 767 ! ****************************************************************************
 768 !
 769
 770   IMPLICIT NONE
 771
 772   INTEGER, INTENT(IN)    :: imx, jmx
 773   INTEGER, INTENT(IN)    :: doy
 774   REAL,    INTENT(IN)    :: xhour
 775   REAL,    INTENT(OUT)   :: sza(imx,jmx), cossza(imx,jmx)
 776
 777   REAL    :: a0, a1, a2, a3, b1, b2, b3, r, dec, timloc, ahr,xlon,rlat
 778   real, parameter :: pi=3.14
 779   INTEGER :: i, j
 780
 781   ! executable statements
 782
 783   ! ***************************************************************************
 784   ! *  Solar declination angle:                                               *
 785   ! ***************************************************************************
 786   a0 = 0.006918
 787   a1 = 0.399912
 788   a2 = 0.006758
 789   a3 = 0.002697
 790   b1 = 0.070257
 791   b2 = 0.000907
 792   b3 = 0.000148
 793   r  = 2.0* pi * REAL(doy-1)/365.0
 794   !
 795   dec = a0 - a1*COS(  r)   + b1*SIN(  r)   &
 796            - a2*COS(2.0*r) + b2*SIN(2.0*r) &
 797            - a3*COS(3.0*r) + b3*SIN(3.0*r)
 798   !
 799   DO i = 1,imx
 800      ! ************************************************************************
 801      ! *  Hour angle (AHR) is a function of longitude.  AHR is zero at        *
 802      ! *  solar noon, and increases by 15 deg for every hour before or        *
 803      ! *  after solar noon.                                                   *
 804      ! ************************************************************************
 805      ! -- Local time in hours
 806      timloc  = xhour + xlon/15.0
 807      !      IF (timloc < 0.0) timloc = 24.0 + timloc
 808      IF (timloc > 24.0) timloc = timloc - 24.0
 809      !
 810      ! -- Hour angle
 811      ahr = ABS(timloc - 12.0) * 15.0 * pi/180.0
 812      !
 813      DO j = 1,jmx
 814         ! -- Solar zenith angle
 815         cossza(i,j) = SIN(rlat) * SIN(dec) + &
 816                       COS(rlat) * COS(dec) * COS(ahr)
 817         sza(i,j)    = ACOS(cossza(i,j)) * 180.0/pi
 818         IF (cossza(i,j) < 0.0)   cossza(i,j) = 0.0
 819         !
 820      END do
 821   END DO
 822
 823 END subroutine szangle
 824
 825 END MODULE MODULE_GOCART_CHEM