#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: EBISCHEME ! ! !DESCRIPTION: Eulerian Backward Iteration (EBI) is a chemistry solver for ! the modified CB05 scheme. !\\ !\\ ! !INTERFACE: ! module EBISCHEME ! ! !USES: ! use GO, only : gol, goPr, goErr use dims, only : nregions use tm5_distgrid, only : dgrid, Get_DistGrid use chem_param #ifdef with_budgets use budget_global, only : nbudg, nbud_vg use reaction_data, only : nreac, nreacw use photolysis_data, only : nj #endif implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: EBI ! ! !PUBLIC DATA MEMBERS: ! #ifdef with_budgets real, dimension(nregions), public :: sum_wetchem real, dimension(nregions), public :: sum_chemistry real, dimension(nregions), public :: sum_deposition type, public :: buddep_data real, dimension(:,:,:), pointer :: dry ! i,j,ntrace end type buddep_data type(buddep_data), dimension(nregions), public :: buddep_dat real, dimension(nbudg, nbud_vg, nreac ), public :: budrrg real, dimension(nbudg, nbud_vg, nj ), public :: budrjg real, dimension(nbudg, nbud_vg, nreacw), public :: budrwg real, dimension(nbudg, nbud_vg, nmark ), public :: budmarkg !allow saving of production of aqueous, gas-phase sulfate, elvoc, svoc type,public :: production_data real, dimension(:,:,:,:), pointer :: prod ! i,j,lev,nprod end type production_data integer,parameter,public :: nprod_AC_o3=2 integer,parameter,public :: n_loss=2 !loss of ch4 and co integer,parameter,public :: nprod=12 type(production_data), dimension(nregions),public:: diag_prod !The chemical production container for AerChemMIP output, needs own container, since the variable is zeroed when ! writing out and could possibly conflict with the general output routine. integer,parameter,public :: nprod_AC=4 ! N=4: GASSO4, AQSO4, SOA(3d monthly),SOA(2d hourly) integer,parameter,public :: iprod_gasso4=1, iprod_aqso4=2,iprod_soa3dmon=3,iprod_soa2dhour=4 integer,parameter,public :: iloss_ch4=1, iloss_co=2 type(production_data), dimension(nregions),public:: AC_diag_prod type(production_data), dimension(nregions),public:: AC_O3_lp type(production_data), dimension(nregions),public:: AC_loss real,dimension(:,:), allocatable :: temp_prod_so4 real,dimension(:,:), allocatable :: temp_loss real,dimension(:,:), allocatable :: temp_prod_vocs #endif logical, public::isoprene_on character(len=*), parameter :: mname = 'ebischeme' ! ! !REVISION HISTORY: ! ! Feb 2007 - Elisabetta Vignati - changed for inclusion of cloud chemistry on modes ! 22 Mar 2012 - Philippe Le Sager - adapted for lon-lat MPI domain decomposition ! 21 mar 2014 - Jason Williams - port the modified CB05 chemistry scheme ! ! !REMARKS: ! (1) initializations are made in chemistry/chemistry_init ! ! !TODO : FIXME check sum_chemistry and sum_deposition ! !EOP !------------------------------------------------------------------------ contains !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EBI ! ! !DESCRIPTION: perform Eulerian backwards chemistry at one model layer ! level in one region !\\ !\\ ! !INTERFACE: ! subroutine EBI( region, level, is, js, rj, rr, y, ye, status) ! ! !USES: ! use dims, only : im, jm, nchem, tref, okdebug use Dims, only : CheckShape use binas, only : avog use chem_param, only : xmso4 use global_data, only : region_dat #ifndef without_dry_deposition use dry_deposition, only : vd #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region, level, is, js real, intent(in) :: rr(is:,js:,:) ! array of reaction rates real, intent(in) :: rj(is:,js:,:) ! array of photolysis rates ! ! !INPUT/OUTPUT PARAMETERS: ! real, intent(inout) :: y(is:,js:,:) ! array of concentration real, intent(inout) :: ye(is:,js:,:) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/ebi' #ifdef with_budgets real, dimension(:,:,:), allocatable :: cr2, cr3, cr4 ! reaction budgets #endif real, dimension(:,:,:), allocatable :: y0 real :: dtime, dt, dt2, fnoy integer :: iterebi, i, j, ib, maxit integer :: io, sfstart, sfend !SOA parameters real,parameter :: total_soa_yield_isoprene=0.01 real,parameter :: total_soa_yield_terpenes=0.15 real,parameter :: y_oh_isop_elvoc=0.0003 real,parameter :: y_o3_isop_elvoc=0.0001 real,parameter :: y_oh_terp_elvoc=0.01 real,parameter :: y_o3_terp_elvoc=0.05 real :: y_oh_isop_svoc ! =0.01-0.0003 real :: y_o3_isop_svoc ! =0.01-0.0001 real :: y_oh_terp_svoc ! =0.15-0.01 real :: y_o3_terp_svoc ! =0.15-0.05 ! For vectorization/blocking .... ! npts can be varied to optimize cache memory management. integer, parameter :: npts=15 integer, dimension(npts) :: ipts, jpts real, dimension(npts, nreac ) :: rrv real, dimension(npts, nj ) :: rjv real, dimension(npts, ntrace) :: vdv ! deposition velocities real, dimension(npts) :: emino real, dimension(npts, maxtrace) :: yv, y0v integer :: iv, itrace, ivt, n integer :: i1, j1, i2, j2 ! --- BEGIN -------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! check arguments ... call CheckShape( (/i2-i1+1, j2-j1+1, nreac /), shape(rr), status ) IF_NOTOK_RETURN(status=1) call CheckShape( (/i2-i1+1, j2-j1+1, nj /), shape(rj), status ) IF_NOTOK_RETURN(status=1) call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status ) IF_NOTOK_RETURN(status=1) call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status ) IF_NOTOK_RETURN(status=1) dtime=nchem/(2*tref(region)) !CMK iterebi=max(1,nint(dtime/2400)) !needed if nchem <2400 iterebi=max(1,nint(dtime/1350)) !needed if nchem <2400 dt=dtime/iterebi ! define svoc yield as remainder total soa yield after elvoc yields ! total yields at present (to check see variable defs): ! 1% for isoprene ! 15% for terpenes y_oh_isop_svoc=total_soa_yield_isoprene - y_oh_isop_elvoc y_o3_isop_svoc=total_soa_yield_isoprene - y_o3_isop_elvoc y_oh_terp_svoc=total_soa_yield_terpenes - y_oh_terp_elvoc y_o3_terp_svoc=total_soa_yield_terpenes - y_o3_terp_elvoc allocate( y0 (i1:i2, j1:j2, maxtrace)) #ifdef with_budgets allocate( cr2(i1:i2, j1:j2, nj )) allocate( cr3(i1:i2, j1:j2, nreac )) allocate( cr4(i1:i2, j1:j2, nreacw )) allocate( temp_prod_so4(npts,2)) allocate( temp_loss(npts,2)) allocate( temp_prod_vocs(npts,8)) #endif !*** SCALING OF NOx, which has changed due to transport/deposition ! This does not yet work. TODO: Make this working. (FIXME : ask VH what's this is about???) do j = j1, j2 do i = i1, i2 y(i,j,ino) =max(1e-3,y(i,j,ino)) y(i,j,ino2) =max(1e-3,y(i,j,ino2)) y(i,j,ino3) =max(1e-3,y(i,j,ino3)) y(i,j,in2o5)=max(1e-3,y(i,j,in2o5)) y(i,j,ihno4)=max(1e-3,y(i,j,ihno4)) fnoy=y(i,j,ino)+y(i,j,ino2)+y(i,j,ino3)+2.*y(i,j,in2o5)+y(i,j,ihno4) fnoy=y(i,j,inox)/fnoy y(i,j,ino) =fnoy*y(i,j,ino) y(i,j,ino2) =fnoy*y(i,j,ino2) y(i,j,ino3) =fnoy*y(i,j,ino3) y(i,j,in2o5)=fnoy*y(i,j,in2o5) y(i,j,ihno4)=fnoy*y(i,j,ihno4) end do end do ! ! set budget accumulators to zero ! #ifdef with_budgets cr2 = 0.0 cr3 = 0.0 cr4 = 0.0 temp_loss = 0.0 temp_prod_so4 = 0.0 temp_prod_vocs = 0.0 #endif !=========================================================== ! ** Start iterating over CHEMISTRY !=========================================================== do ib=1,iterebi maxit=8 !CMKTEMP if(level<=3) maxit = maxit*2 ! lowest layers more iterations y0 = y !------------------------------- ! wet sulphur/ammonia chemistry !------------------------------ #ifdef with_budgets call wetS(region, level, i1, j1, y0, dt, y, ye, cr4, status) #else call wetS(region, level, i1, j1, y0, dt, y, ye, status) #endif !------------------------------------- ! gasphase chemistry using EBI solver !------------------------------------- y0 = y ! reset initial concentrations ! ______do EBI solver_______ dt2 = dt*dt ! block the input for EBI for efficiency: ! copy values with faster running index in inside loop iv = 0 do j = j1, j2 do i = i1, i2 iv = iv+1 ipts(iv) = i jpts(iv) = j if(iv==npts) then ! copy reaction rates... do itrace=1,nreac do ivt=1,npts rrv(ivt,itrace) = rr(ipts(ivt),jpts(ivt),itrace) end do end do ! copy photolysis rates.... do itrace=1,nj do ivt=1,npts rjv(ivt,itrace) = rj(ipts(ivt),jpts(ivt),itrace) end do end do ! copy tracer concentrations .... do itrace=1,maxtrace do ivt=1,npts yv(ivt,itrace) = y(ipts(ivt),jpts(ivt),itrace) y0v(ivt,itrace) = y0(ipts(ivt),jpts(ivt),itrace) end do end do ! deposition .... #ifndef without_dry_deposition if(level == 1) then do itrace=1,ntrace do ivt=1,npts vdv(ivt,itrace) = & vd(region,itrace)%surf(ipts(ivt),jpts(ivt)) & / ye(ipts(ivt),jpts(ivt),idz) !1/s end do end do else vdv(:,:) = 0.0 end if #endif ! copy nox emissions.... do ivt=1,npts emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno) end do ! do ebi solver.... call DO_EBI(iv, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v) do itrace=1,maxtrace do ivt=1,npts y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace) end do end do #ifdef with_budgets #ifdef with_m7 do ivt=1,npts ! add change in gas-phase so4 production to diagnostic and change to mole->mass kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4!kg AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,iprod_gasso4)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,iprod_gasso4)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)+temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)+temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)+temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)+temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)+temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)+temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)+temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc!kg diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)+temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc!kg end do #endif #endif iv=0 end if end do end do ! do the 'remaining' points... if(iv > 0) then do itrace=1,nreac do ivt=1,iv rrv(ivt,itrace) = rr(ipts(ivt),jpts(ivt),itrace) end do end do do itrace=1,nj do ivt=1,iv rjv(ivt,itrace) = rj(ipts(ivt),jpts(ivt),itrace) end do end do do itrace=1,maxtrace do ivt=1,iv yv(ivt,itrace) = y(ipts(ivt),jpts(ivt),itrace) y0v(ivt,itrace) = y0(ipts(ivt),jpts(ivt),itrace) end do end do ! deposition .... #ifndef without_dry_deposition if(level == 1) then do itrace=1,ntrace do ivt=1,iv vdv(ivt,itrace) = & vd(region,itrace)%surf(ipts(ivt),jpts(ivt)) & / ye(ipts(ivt),jpts(ivt),idz) !1/s end do end do else vdv(:,:) = 0.0 end if #endif do ivt=1,iv emino(ivt) = ye(ipts(ivt),jpts(ivt),ieno) end do !the actual EBI solver on remaining cells call DO_EBI(iv, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v) do itrace=1,maxtrace do ivt=1,iv y(ipts(ivt),jpts(ivt),itrace)=yv(ivt,itrace) end do end do #ifdef with_budgets #ifdef with_m7 do ivt=1,iv ! add change in so4 production to diagnostic and change to molec->mass kg ! diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)= diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,1)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4!kg AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,iprod_gasso4)=AC_diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,iprod_gasso4)+temp_prod_so4(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmso4!kg ! SAVE SOA twice, one for hourly output and once fro monthly. Two variables needed because zeroing happens at different time steps diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,5)+temp_prod_vocs(ivt,1)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,6)+temp_prod_vocs(ivt,2)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,7)+temp_prod_vocs(ivt,3)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,8)+temp_prod_vocs(ivt,4)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmelvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,9)+temp_prod_vocs(ivt,5)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,10)+temp_prod_vocs(ivt,6)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,11)+temp_prod_vocs(ivt,7)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)=diag_prod(region)%prod(ipts(ivt),jpts(ivt),level,12)+temp_prod_vocs(ivt,8)/y(ipts(ivt),jpts(ivt),iair)*ye(ipts(ivt),jpts(ivt),iairm)/xmair*xmsvoc end do #endif #endif end if call NOYmass !------------------------------------- ! marked tracers ! apply after correction of nitrogen components !---------------------------------------------- call MARK_TRAC(region, level, is, js, y, rr, rj, dt, ye) !------------------------------------------------------------ ! increase budget accumulators cr2 and cr3 (cr4 is done in wetS) !------------------------------------------------------------ #ifdef with_budgets call incc2c3 #endif !=========================================================== ! ** END iterating over CHEMISTRY !=========================================================== end do !iterebi #ifdef with_budgets call REACBUD #endif deallocate(y0) #ifdef with_budgets deallocate(cr2) deallocate(cr3) deallocate(cr4) deallocate(temp_prod_so4) deallocate(temp_prod_vocs) deallocate(temp_loss) #endif ! ok status = 0 contains subroutine DO_EBI(lvec, npts, maxit, dt, dt2, rrv, rjv, vdv, emino, yv, y0v) ! INPUT/OUTPUT integer,intent(in) :: lvec, npts, maxit real, intent(in) :: dt, dt2, rrv(npts,nreac), rjv(npts,nj), & vdv(npts,ntrace), emino(npts), & y0v(npts,maxtrace) real, intent(out) :: yv(npts,maxtrace) ! Local integer :: ivec,iter real :: r57, r89, p1, r12, r21, xl1, p2, xl2, p3, p32, & xl3, x1, x2, x3, c1, c2, c3, y2, xjt,r21t, & r12t, acub, bcub,ccub,cubdet,dno2,r56, r65, r75,p5, & xl5, r66, x5, p6, xl6, x6, c6, x17, y1, c7, & r98, p8, xl8, x4, c5, xl9, r101,r102,xl7, & c8 , x101, r1920, r1919, p19, xl19,r2019, xl20, & xlhno3, ph2o2, xlh2o2, pch2o, pco, phno3, xlch2o, & pch3o2, xlch3o2, pch3o2h, xlch3o2h, pald2, & xlald2, pmgly, xlmgly, pole, xleth, xlole, & xlisop, prxpar, xlrxpar, ppar, xlpar, pror, & xlror, pxo2,xlxo2,pxo2n, xlxo2n, prooh,xlterp, & xlethooh,xlethp,phcooh,pmcooh,xlc3h6, xlrooh, & porgntr, xlorgntr, xlco, qdms, pso2, qso2, qso2d, & qnh3, qnh2o2, ppnh3, dnh3, pnh2, qnh2, qdms1, qdms2, & pmsa, pnh3, pispd,xlispd,xlacet,paco2,xlaco2, & pch3oh, pc3h7o2,phypro2,xlc3h7o2,xlhypro2,pacet, & pelvoc,psvoc!RM do iter=1,maxit do ivec=1,lvec ! --- Short living compounds & groups ! --- First group: NO NO2 O3 P1=rjv(ivec,jbno3)*yv(ivec,ino3)+rjv(ivec,jn2o5b)*yv(ivec,in2o5)& +rjv(ivec,jhono)*yv(ivec,ihono)+rrv(ivec,knh2o2)*yv(ivec,inh2)+emino(ivec) R12=0. R21=rrv(ivec,kho2no)*yv(ivec,iho2)+rrv(ivec,kmo2no)*yv(ivec,ich3o2)& +rrv(ivec,kc79)*yv(ivec,ixo2)+rrv(ivec,kc46)*yv(ivec,ic2o3) & +rrv(ivec,kaco2no)*yv(ivec,iaco2)+rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)& +rrv(ivec,khypo2no)*yv(ivec,ihypro2)+rrv(ivec,knh2o2no)*yv(ivec,inh2o2) XL1=rrv(ivec,knono3)*yv(ivec,ino3)+rrv(ivec,kc81)*yv(ivec,ixo2n)& +rrv(ivec,knh2no)*yv(ivec,inh2)+rrv(ivec,khono)*yv(ivec,ioh) XL1 = XL1 + vdv(ivec,ino) P2=rjv(ivec,jhno3)*yv(ivec,ihno3)+rjv(ivec,jn2o5a)*yv(ivec,in2o5)& +rrv(ivec,kn2o5)*yv(ivec,in2o5)+rjv(ivec,jano3)*yv(ivec,ino3)& +yv(ivec,ihno4)*(rjv(ivec,jhno4)+rrv(ivec,khno4m)+rrv(ivec,khno4oh)*yv(ivec,ioh))& +rrv(ivec,kohhono)*yv(ivec,ioh)*yv(ivec,ihono)& +2.*rrv(ivec,knono3)*yv(ivec,ino3)*yv(ivec,ino)& +rrv(ivec,kc48)*yv(ivec,ipan)+rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +rjv(ivec,jorgn)*yv(ivec,iorgntr)& +rjv(ivec,jpana)*yv(ivec,ipan)& +0.2*rrv(ivec,kc78) *yv(ivec,iisop)*yv(ivec,ino3)& +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)*yv(ivec,ino3)& +rrv(ivec,kno3c2o3)*yv(ivec,ic2o3)*yv(ivec,ino3)& +rrv(ivec,kno3xo2)*yv(ivec,ixo2)*yv(ivec,ino3)& +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +rjv(ivec,jmo2no2a)*yv(ivec,ich3o2no2)& +rrv(ivec,kch3o2no2m)*yv(ivec,ich3o2no2) XL2=rrv(ivec,kno2oh)*yv(ivec,ioh)+rrv(ivec,kno2no3)*yv(ivec,ino3)& +rrv(ivec,kno2ho2)*yv(ivec,iho2)+rrv(ivec,kno2o3)*yv(ivec,io3)& +rrv(ivec,kc47)*yv(ivec,ic2o3)+rrv(ivec,knh2no2)*yv(ivec,inh2)& +rrv(ivec,kmo2no2)*yv(ivec,ich3o2) XL2 = XL2 + vdv(ivec,ino2) P3=rjv(ivec,jano3)*yv(ivec,ino3)+rjv(ivec,jo2) ! JEW : O3P + O2 = O3 XL3=rrv(ivec,ko3ho2)*yv(ivec,iho2)+rrv(ivec,ko3oh)*yv(ivec,ioh)& +rrv(ivec,kno2o3)*yv(ivec,ino2)+rjv(ivec,jo3d)+rrv(ivec,ko3po3)& +rrv(ivec,kc58)*yv(ivec,iole)& +rrv(ivec,kc62)*yv(ivec,ieth)& +rrv(ivec,kc77)*yv(ivec,iisop)& +rrv(ivec,ko3c3h6)*yv(ivec,ic3h6)& +rrv(ivec,ko3terp)*yv(ivec,iterp)& +rrv(ivec,ko3ispd)*yv(ivec,iispd)& +rrv(ivec,knh2o3)*yv(ivec,inh2)& +rrv(ivec,knh2o2o3)*yv(ivec,inh2o2) XL3 = XL3 + vdv(ivec,io3) X1=y0v(ivec,ino)+P1*DT X2=y0v(ivec,ino2)+P2*DT X3=y0v(ivec,io3)+P3*DT C1=1.+XL1*DT C2=1.+XL2*DT C3=1.+XL3*DT Y1=rrv(ivec,knoo3)*DT R21T=R21*DT R12T=R12*DT XJT=rjv(ivec,jno2)*DT P32=0.4*rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2) ! --- Solve unknown: x ACUB=-2.*Y1*(C2+R12T+C2*R21T/C1) BCUB=2.*C1*C2*C3+2.*C1*C3*(R12T+XJT)+2.*C2*C3*R21T+& 2.*Y1*(R12T*(X1-X2)+2.*C2*R21T*X1/C1+C2*(X1+X3)) CCUB=2.*C1*C3*X2*(R12T+XJT)-2.*C2*C3*X1*R21T+2.*Y1*X1*& (X2*R12T-C2*X3-C2*R21T*X1/C1) CUBDET=BCUB*BCUB-4.*ACUB*CCUB DNO2=(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB) dno2=min(x1,dno2) yv(ivec,ino2)=(X2+DNO2)/C2 yv(ivec,ino)=(X1-DNO2)/C1 yv(ivec,io3)=(X3+(P32*dt)+XJT*yv(ivec,ino2))/(C3+Y1*yv(ivec,ino)) ! --- Second group: yv(ivec,iho2) yv(ivec,ioh) yv(ivec,ihno4) yv(ivec,ihono) R57=rjv(ivec,jhno4)+rrv(ivec,khno4m) R56=rrv(ivec,kcooh)*yv(ivec,ico)+rrv(ivec,ko3oh)*yv(ivec,io3)+rrv(ivec,khpoh)& *yv(ivec,ih2o2)+rrv(ivec,kfrmoh)*yv(ivec,ich2o)+rrv(ivec,kh2oh)& +rrv(ivec,kso2oh)*yv(ivec,iso2) p5=2.*rjv(ivec,jbch2o)*yv(ivec,ich2o)& +rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)& +rjv(ivec,jmepe)*yv(ivec,ich3o2h)& +2.0*rjv(ivec,jald2)*yv(ivec,iald2)& +rjv(ivec,jmgly)*yv(ivec,imgly)+0.11*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)& +0.94*rrv(ivec,kc53)*yv(ivec,iror)+rrv(ivec,kc54)*yv(ivec,iror)& +1.57*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +0.76*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)& +0.56*rrv(ivec,kc59)*yv(ivec,ino3)*yv(ivec,iole)& +rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)& +0.22*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)& +0.066*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)& +0.9*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +0.9*rjv(ivec,jorgn)*yv(ivec,iorgntr)& +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)& +0.74*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2)& +0.9*rjv(ivec,jrooh)*yv(ivec,irooh)& +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)*yv(ivec,ino3)& +0.19*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +0.28*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.75*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +0.154*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.925*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +1.033*rjv(ivec,jispd)*yv(ivec,iispd)& +rrv(ivec,kohch3oh)*yv(ivec,ich3oh)& +rrv(ivec,kohhcooh)*yv(ivec,ihcooh)& +rrv(ivec,kohethoh)*yv(ivec,iethoh)& +1.22*rrv(ivec,kohterp)*yv(ivec,iterp)& +rrv(ivec,kohc2h6)*yv(ivec,ic2h6)& +0.5*rrv(ivec,kc76)*yv(ivec,ioh)*yv(ivec,iisop)& ! reduced according to archibald et al, AE, 2011 +0.503*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.5*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)& +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)& +rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)*yv(ivec,ino)& +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2) XL5=rrv(ivec,kho2no) *yv(ivec,ino)& +rrv(ivec,kno2ho2) *yv(ivec,ino2)& +rrv(ivec,ko3ho2) *yv(ivec,io3)& +rrv(ivec,kmo2ho2a) *yv(ivec,ich3o2)& +rrv(ivec,kmo2ho2b) *yv(ivec,ich3o2)& +rrv(ivec,kc50) *yv(ivec,ic2o3)& +rrv(ivec,kho2oh) *yv(ivec,ioh)& +rrv(ivec,kc82) *yv(ivec,ixo2)& +rrv(ivec,kc85) *yv(ivec,ixo2n)& +rrv(ivec,kno3ho2) *yv(ivec,ino3)& +rrv(ivec,kaco2ho2)*yv(ivec,iaco2)& +rrv(ivec,kc3h7o2ho2)*yv(ivec,ic3h7o2)& +rrv(ivec,khypo2ho2)*yv(ivec,ihypro2)& +rrv(ivec,knh2ho2)*yv(ivec,inh2)& +rrv(ivec,knh2o2ho2)*yv(ivec,inh2o2)& +rrv(ivec,kho2_aer) & +rrv(ivec,kho2_l) R66=2.*rrv(ivec,kho2ho2) X5=y0v(ivec,iho2)+P5*DT R65=rrv(ivec,kho2no)*yv(ivec,ino)+rrv(ivec,ko3ho2)*yv(ivec,io3) P6=rjv(ivec,jhno3)*yv(ivec,ihno3)& +2.*rjv(ivec,jo3d)*yv(ivec,io3)& +2.*rjv(ivec,jh2o2)*yv(ivec,ih2o2)& +rjv(ivec,jmepe)*yv(ivec,ich3o2h)& +0.1*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)& +0.266*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +rjv(ivec,jrooh)*yv(ivec,irooh)& +0.12*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)& +0.33*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +0.57*rrv(ivec,ko3terp)*yv(ivec,iterp)*yv(ivec,io3)& +0.268*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd) XL6=rrv(ivec,khno4oh)*yv(ivec,ihno4)& +rrv(ivec,kohhono)*yv(ivec,ihono)& +rrv(ivec,kho2oh)*yv(ivec,iho2)& +rrv(ivec,kno2oh)*yv(ivec,ino2)& +rrv(ivec,khono)*yv(ivec,ino)& +rrv(ivec,kohhno3)*yv(ivec,ihno3)& +rrv(ivec,kcooh)*yv(ivec,ico)& +rrv(ivec,ko3oh)*yv(ivec,io3)& +rrv(ivec,khpoh)*yv(ivec,ih2o2)& +rrv(ivec,kfrmoh)*yv(ivec,ich2o)& +rrv(ivec,kch4oh)*yv(ivec,ich4)& +0.7*rrv(ivec,kohmper)*yv(ivec,ich3o2h)& +rrv(ivec,kc43)*yv(ivec,iald2)& +rrv(ivec,kc73)*yv(ivec,imgly)& +rrv(ivec,kc52)*yv(ivec,ipar)& +rrv(ivec,kc57)*yv(ivec,iole)& +rrv(ivec,kc61)*yv(ivec,ieth)& +rrv(ivec,kc76)*yv(ivec,iisop)& +0.77*rrv(ivec,kohrooh)*yv(ivec,irooh)& ! note: change from '0.7' to '0.77' +rrv(ivec,kc84)*yv(ivec,iorgntr)& +rrv(ivec,kh2oh)& +rrv(ivec,kso2oh)*yv(ivec,iso2) & ! bug found by Jason 01/2008 +(rrv(ivec,kdmsoha)+rrv(ivec,kdmsohb)) *yv(ivec,idms) & +rrv(ivec,knh3oh)*yv(ivec,inh3)& +rrv(ivec,knh2oh)*yv(ivec,inh2)& +rrv(ivec,kohch3oh)*yv(ivec,ich3oh)& +rrv(ivec,kohhcooh)*yv(ivec,ihcooh)& +rrv(ivec,kohethoh)*yv(ivec,iethoh)& +rrv(ivec,kohterp)*yv(ivec,iterp)& +rrv(ivec,kohispd)*yv(ivec,iispd)& +rrv(ivec,kohmcooh)*yv(ivec,imcooh)& +rrv(ivec,kohc2h6)*yv(ivec,ic2h6)& +rrv(ivec,kohc3h8)*yv(ivec,ic3h8)& +rrv(ivec,kohc3h6)*yv(ivec,ic3h6)& +rrv(ivec,kohacet)*yv(ivec,iacet) R101=rjv(ivec,jhono) X6=y0v(ivec,ioh)+(P6*DT)+(R101*y0v(ivec,ihono)*DT) C6=1.+XL6*DT R75=rrv(ivec,kno2ho2)*yv(ivec,ino2) R102=rrv(ivec,khono)*yv(ivec,ino) X101=rjv(ivec,jhono)+rrv(ivec,kohhono)*yv(ivec,ioh) X101= X101 + vdv(ivec,ihono) C8=1.+X101*DT XL7=rjv(ivec,jhno4)+rrv(ivec,khno4oh)*yv(ivec,ioh)+rrv(ivec,khno4m) XL7 = XL7 + vdv(ivec,ihno4) C7=1.+XL7*DT Y1=R57/C7 Y2=R56/C6 ACUB=R66*DT BCUB=1.+XL5*DT-DT2*(Y1*R75+Y2*R65) CCUB=-1.*X5-DT*(Y1*y0v(ivec,ihno4)+Y2*X6) CUBDET=BCUB*BCUB-4.*ACUB*CCUB CUBDET=MAX(CUBDET,1.E-20) yv(ivec,iho2)=max(0.1,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB)) yv(ivec,ioh)=(X6+R65*yv(ivec,iho2)*DT)/C6 yv(ivec,ihono)=(y0v(ivec,ihono)+R102*DT*yv(ivec,ioh))/C8 yv(ivec,ihno4)=(y0v(ivec,ihno4)+R75*DT*yv(ivec,iho2))/C7 ! ! --- Third group: NO3 N2O5 ! R89=rjv(ivec,jn2o5a)+rjv(ivec,jn2o5b)+rrv(ivec,kn2o5) P8=rrv(ivec,kohhno3)*yv(ivec,ihno3)*yv(ivec,ioh)+rrv(ivec,kno2o3)*yv(ivec,ino2)*yv(ivec,io3)& +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2)+rjv(ivec,jpanb)*yv(ivec,ipan) XL8=rjv(ivec,jbno3)+rjv(ivec,jano3)& +rrv(ivec,knono3)*yv(ivec,ino)& +rrv(ivec,kno2no3)*yv(ivec,ino2)& +rrv(ivec,kc44)*yv(ivec,iald2)& +rrv(ivec,kc59)*yv(ivec,iole)& +rrv(ivec,kc78)*yv(ivec,iisop)& +rrv(ivec,kc41)*yv(ivec,ich2o)& +rrv(ivec,kdmsno3)*yv(ivec,idms)& +rrv(ivec,kno3_aer)& +rrv(ivec,kno3ho2)*yv(ivec,iho2)& +rrv(ivec,kno3mo2)*yv(ivec,ich3o2)& +rrv(ivec,kno3c2o3)*yv(ivec,ic2o3)& +rrv(ivec,kno3xo2)*yv(ivec,ixo2)& +rrv(ivec,kno3c3h6)*yv(ivec,ic3h6)& +rrv(ivec,kno3terp)*yv(ivec,iterp)& +rrv(ivec,kno3ispd)*yv(ivec,iispd) XL8 = XL8 + vdv(ivec,ino3) X4=y0v(ivec,ino3)+P8*DT C5=1.+XL8*DT R98=rrv(ivec,kno2no3)*yv(ivec,ino2) XL9=rjv(ivec,jn2o5a)+rjv(ivec,jn2o5b)+rrv(ivec,kn2o5)+rrv(ivec,kn2o5_aer)+rrv(ivec,kn2o5l) !cmk rates now idependent from y XL9 = XL9 + vdv(ivec,in2o5) C6=1.+XL9*DT C7=(C5*C6-R89*R98*DT2) yv(ivec,in2o5)=(C5*y0v(ivec,in2o5)+R98*DT*X4)/C7 yv(ivec,ino3)=(C6*X4+R89*DT*y0v(ivec,in2o5))/C7 ! ! --- Fourth group: C2O3 PAN ! R1920=rrv(ivec,kc48)+rjv(ivec,jpana) R1919=rrv(ivec,kc49) P19=rrv(ivec,kc43)*yv(ivec,iald2)*yv(ivec,ioh)& +rrv(ivec,kc44)*yv(ivec,iald2)*yv(ivec,ino3)& +rrv(ivec,kc73)*yv(ivec,imgly)*yv(ivec,ioh) & +rjv(ivec,jmgly)*yv(ivec,imgly)& +0.2*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +0.74*rjv(ivec,jrooh)*yv(ivec,irooh)& +0.74*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +0.74*rjv(ivec,jorgn)*yv(ivec,iorgntr)& +0.39*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.498*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.114*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.075*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.967*rjv(ivec,jispd)*yv(ivec,iispd)& +0.3*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)& +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)& +rjv(ivec,jb_acet)*yv(ivec,iacet)& +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2) XL19=rrv(ivec,kc46)*yv(ivec,ino)+rrv(ivec,kc50)*yv(ivec,iho2)& +rrv(ivec,kc47)*yv(ivec,ino2)& +rrv(ivec,kno3c2o3)*yv(ivec,ino3) XL19 = XL19 + vdv(ivec,ic2o3) R2019=rrv(ivec,kc47)*yv(ivec,ino2) XL20=rrv(ivec,kc48)+rjv(ivec,jpana)+rjv(ivec,jpanb) XL20 = XL20 + vdv(ivec,ipan) ACUB=2*R1919*DT*(1+XL20*DT) BCUB=(1+XL20*DT)*(1+XL19*DT)-R1920*DT*R2019*DT CCUB=(1+XL20*DT)*(y0v(ivec,ic2o3)+P19*DT)+R1920*DT*y0v(ivec,ipan) CUBDET=BCUB*BCUB+4.*ACUB*CCUB yv(ivec,ic2o3)=max(1e-8,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB)) !cmk put max here.... yv(ivec,ipan)=(y0v(ivec,ipan)+R2019*yv(ivec,ic2o3)*DT)/(1.+XL20*DT) ! ! ---- Fifth group : CH3O2 CH3O2NO2 R1920=rrv(ivec,kch3o2no2m)+rjv(ivec,jmo2no2a)+rjv(ivec,jmo2no2b) R1919=rrv(ivec,kmo2mo2) P19=rrv(ivec,kch4oh)*yv(ivec,ich4)*yv(ivec,ioh) & +0.7*rrv(ivec,kohmper)*yv(ivec,ioh)*yv(ivec,ich3o2h)& +rrv(ivec,kno3c2o3)*yv(ivec,ino3)*yv(ivec,ic2o3)& +rrv(ivec,kc46)*yv(ivec,ino)*yv(ivec,ic2o3)& +2*rrv(ivec,kc49)*yv(ivec,ic2o3)*yv(ivec,ic2o3)& +rjv(ivec,jpanb)*yv(ivec,ipan)& +rjv(ivec,jald2)*yv(ivec,iald2)& +0.74*rjv(ivec,jrooh)*yv(ivec,irooh)& +0.74*rrv(ivec,kc84)*yv(ivec,iorgntr)& +0.74*rjv(ivec,jorgn)*yv(ivec,iorgntr)& +rrv(ivec,kohmcooh)*yv(ivec,ioh)*yv(ivec,imcooh)& +0.31*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +0.39*rrv(ivec,ko3terp)*yv(ivec,iterp)*yv(ivec,io3)& +2.0*rjv(ivec,ja_acet)*yv(ivec,iacet)+rjv(ivec,jb_acet)*yv(ivec,iacet) XL19=rrv(ivec,kmo2no)*yv(ivec,ino)+rrv(ivec,kmo2ho2a)*yv(ivec,iho2)& +rrv(ivec,kmo2ho2b)*yv(ivec,iho2)+rrv(ivec,kno3mo2)*yv(ivec,ino3)& +rrv(ivec,kaco2mo2)*yv(ivec,iaco2)+rrv(ivec,kmo2no2)*yv(ivec,ino2) R2019=rrv(ivec,kmo2no2)*yv(ivec,ino2) XL20=rrv(ivec,kch3o2no2m)+rjv(ivec,jmo2no2a)+rjv(ivec,jmo2no2b) XL20 = XL20 + vdv(ivec,ich3o2no2) ACUB=2*R1919*DT*(1+XL20*DT) BCUB=(1+XL20*DT)*(1+XL19*DT)-R1920*DT*R2019*DT CCUB=(1+XL20*DT)*(y0v(ivec,ich3o2)+P19*DT)+R1920*DT*y0v(ivec,ich3o2no2) CUBDET=BCUB*BCUB+4.*ACUB*CCUB yv(ivec,ich3o2)=max(1e-8,(-1.*BCUB+SQRT(CUBDET))/(2.*ACUB)) !cmk put max here.... yv(ivec,ich3o2no2)=(y0v(ivec,ich3o2no2)+R2019*yv(ivec,ich3o2)*DT)/(1.+XL20*DT) ! ! -------- ISPD chemistry ! pispd=0.912*rrv(ivec,kc76)*yv(ivec,ioh)*yv(ivec,iisop)& +0.65*rrv(ivec,kc77)*yv(ivec,io3)*yv(ivec,iisop)& +0.2*rrv(ivec,kc78)*yv(ivec,ino3)*yv(ivec,iisop) xlispd=rrv(ivec,kohispd)*yv(ivec,ioh)+rrv(ivec,ko3ispd)*yv(ivec,io3)& +rrv(ivec,kno3ispd)*yv(ivec,ino3)+rjv(ivec,jispd)+vdv(ivec,iispd) yv(ivec,iispd)=(y0v(ivec,iispd)+pispd*DT)/(1.+xlispd*DT) ! ! -------- ACO2 chemistry ! paco2=rrv(ivec,kohacet)*yv(ivec,iacet)*yv(ivec,ioh) xlaco2=rrv(ivec,kaco2ho2)*yv(ivec,iho2)+rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)& +rrv(ivec,kaco2no)*yv(ivec,ino)+rrv(ivec,kaco2xo2)*yv(ivec,ixo2) yv(ivec,iaco2)=(y0v(ivec,iaco2)+paco2*DT)/(1.+xlaco2*DT) ! ! -------- C3H7O2 chemistry ! pc3h7o2=rrv(ivec,kohc3h8)*yv(ivec,ic3h8)*yv(ivec,ioh) xlc3h7o2=rrv(ivec,kc3h7o2no)*yv(ivec,ino)+rrv(ivec,kc3h7o2ho2)*yv(ivec,iho2) yv(ivec,ic3h7o2)=(y0v(ivec,ic3h7o2)+pc3h7o2*DT)/(1.+xlc3h7o2*DT) ! ! -------- HYPRO2 chemistry ! phypro2=rrv(ivec,kohc3h6)*yv(ivec,ic3h6)*yv(ivec,ioh) xlhypro2=rrv(ivec,khypo2no)*yv(ivec,ino)+rrv(ivec,khypo2ho2)*yv(ivec,iho2) yv(ivec,ihypro2)=(y0v(ivec,ihypro2)+phypro2*DT)/(1.+xlhypro2*DT) ! ! ! --- CBM4 chem.(short living compounds & operators) ! PRXPAR=0.11*rrv(ivec,kc52)*yv(ivec,ioh)*yv(ivec,ipar)& +2.1*rrv(ivec,kc53)*yv(ivec,iror)& +0.7*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)& +rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)& +1.98*rjv(ivec,jrooh)*yv(ivec,irooh)& +1.98*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +1.98*rjv(ivec,jorgn)*yv(ivec,iorgntr) XLRXPAR=rrv(ivec,kc83)*yv(ivec,ipar) yv(ivec,irxpar)=(y0v(ivec,irxpar)+PRXPAR*DT)/(1.+XLRXPAR*DT) XLISOP=rrv(ivec,kc76)*yv(ivec,ioh)+rrv(ivec,kc77)*yv(ivec,io3)+rrv(ivec,kc78)*yv(ivec,ino3) yv(ivec,iisop)=y0v(ivec,iisop)/(1.+XLISOP*DT) PROR=0.76*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)+0.02*rrv(ivec,kc53)*yv(ivec,iror) XLROR=rrv(ivec,kc53)+rrv(ivec,kc54) yv(ivec,iror)=(y0v(ivec,iror)+PROR*DT)/(1.+XLROR*DT) xlterp=rrv(ivec,kohterp)*yv(ivec,ioh)+rrv(ivec,ko3terp)*yv(ivec,io3)& +rrv(ivec,kno3terp)*yv(ivec,ino3) yv(ivec,iterp)=y0v(ivec,iterp)/(1.+xlterp*dt) #ifdef with_m7 !ELVOC ! contribution from monoterpene pelvoc=y_oh_terp_elvoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)+y_o3_terp_elvoc*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp) #ifdef with_budgets temp_prod_vocs(ivec,1)=y_oh_terp_elvoc * rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)*DT temp_prod_vocs(ivec,2)=y_o3_terp_elvoc * rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)*DT #endif !Contribution from isoprene if (isoprene_on)then pelvoc=pelvoc+ y_oh_isop_elvoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)+y_o3_isop_elvoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop) #ifdef with_budgets temp_prod_vocs(ivec,3) = y_oh_isop_elvoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)*DT temp_prod_vocs(ivec,4) = y_o3_isop_elvoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)*DT #endif endif yv(ivec,ielvoc)=y0v(ivec,ielvoc)+pelvoc*dt !SVOC ! contribution from monoterpene psvoc=y_oh_terp_svoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)+y_o3_terp_svoc*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp) #ifdef with_budgets temp_prod_vocs(ivec,5) = y_oh_terp_svoc*rrv(ivec,kohterp)*yv(ivec,ioh)*y0v(ivec,iterp)*DT temp_prod_vocs(ivec,6) = y_o3_terp_svoc*rrv(ivec,ko3terp)*yv(ivec,io3)*y0v(ivec,iterp)*DT #endif !Contribution from isoprene if (isoprene_on)then psvoc=psvoc + y_oh_isop_svoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)+y_o3_isop_svoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop) #ifdef with_budgets temp_prod_vocs(ivec,7)= y_oh_isop_svoc*rrv(ivec,kc76)*yv(ivec,ioh)*y0v(ivec,iisop)*DT temp_prod_vocs(ivec,8)= y_o3_isop_svoc*rrv(ivec,kc77)*yv(ivec,io3)*y0v(ivec,iisop)*DT #endif end if yv(ivec,isvoc)=y0v(ivec,isvoc)+psvoc*dt #endif pxo2=rrv(ivec,kc73)*yv(ivec,imgly)*yv(ivec,ioh)& +0.87*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)& +0.96*rrv(ivec,kc53)*yv(ivec,iror)& +0.8*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +0.22*rrv(ivec,kc58)*yv(ivec,io3)*yv(ivec,iole)& +0.91*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)& +0.991*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)& +rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)& +0.77*rrv(ivec,kohrooh)*yv(ivec,irooh)*yv(ivec,ioh)& +0.5*rjv(ivec,jrooh)*yv(ivec,irooh)& +0.51*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +0.51*rjv(ivec,jorgn)*yv(ivec,iorgntr)& +0.2*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +0.1*rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)& +0.991*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)& +1.25*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +0.76*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +1.03*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +0.713*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.064*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.075*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.7*rjv(ivec,jispd)*yv(ivec,iispd) xlxo2=rrv(ivec,kc79)*yv(ivec,ino)+2.*rrv(ivec,kc80)*yv(ivec,ixo2)& +rrv(ivec,kc82)*yv(ivec,iho2)+rrv(ivec,kno3xo2)*yv(ivec,ino3)& +rrv(ivec,kaco2xo2)*yv(ivec,iaco2)+rrv(ivec,kxo2xo2n)*yv(ivec,ixo2n) yv(ivec,ixo2)=(y0v(ivec,ixo2)+pxo2*dt)/(1.+xlxo2*dt) pxo2n=0.13*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)& +0.04*rrv(ivec,kc53)*yv(ivec,iror)& +0.09*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +0.009*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)& +0.088*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)& +0.25*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +0.18*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.25*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp) xlxo2n=rrv(ivec,kc81)*yv(ivec,ino)+rrv(ivec,kc85)*yv(ivec,iho2)& +rrv(ivec,kxo2xo2n)*yv(ivec,ixo2n)+rrv(ivec,kxo2n)*yv(ivec,ixo2n) yv(ivec,ixo2n)=(y0v(ivec,ixo2n)+pxo2n*dt)/(1.+xlxo2n*dt) end do !ivec if ( mod(iter,2) == 0 ) then do ivec=1,lvec ! --- Species with intermediate lifetimes ! --- Inorganic compounds (HNO3 H2O2) ! PHNO3=rrv(ivec,kno2oh)*yv(ivec,ino2)*yv(ivec,ioh)& +2.*(rrv(ivec,kn2o5_aer)+rrv(ivec,kn2o5l))*yv(ivec,in2o5)& +rrv(ivec,kc44)*yv(ivec,iald2)*yv(ivec,ino3)& +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)& +rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +rrv(ivec,kno3ho2)*yv(ivec,ino3)*yv(ivec,iho2)& +0.15*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +rrv(ivec,kdmsno3)*yv(ivec,ino3)*yv(ivec,idms)& +rrv(ivec,kno3_aer)*yv(ivec,ino3) XLHNO3=rjv(ivec,jhno3)+rrv(ivec,kohhno3)*yv(ivec,ioh) XLHNO3=XLHNO3 + vdv(ivec,ihno3) yv(ivec,ihno3)=(y0v(ivec,ihno3)+PHNO3*DT)/(1.+XLHNO3*DT) PH2O2=rrv(ivec,kho2ho2)*yv(ivec,iho2)*yv(ivec,iho2)& +0.5*rrv(ivec,kho2_aer)*yv(ivec,iho2) XLH2O2=rjv(ivec,jh2o2)+rrv(ivec,khpoh)*yv(ivec,ioh) XLH2O2=XLH2O2 + vdv(ivec,ih2o2) yv(ivec,ih2o2)=(y0v(ivec,ih2o2)+PH2O2*DT)/(1.+XLH2O2*DT) ! --- CH4-chemistry (methyl peroxide formaldehyde) PCH3O2H=rrv(ivec,kmo2ho2a)*yv(ivec,ich3o2)*yv(ivec,iho2) XLCH3O2H=rrv(ivec,kohmper)*yv(ivec,ioh)+rjv(ivec,jmepe) XLCH3O2H=XLCH3O2H + vdv(ivec,ich3o2h) yv(ivec,ich3o2h)=(y0v(ivec,ich3o2h)+PCH3O2H*DT)/(1.+XLCH3O2H*DT) pch2o=0.3*rrv(ivec,kohmper)*yv(ivec,ich3o2h)*yv(ivec,ioh)& +rrv(ivec,kmo2no)*yv(ivec,ich3o2)*yv(ivec,ino)& +rrv(ivec,kmo2ho2b)*yv(ivec,ich3o2)*yv(ivec,iho2)& +1.37*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2)& +rjv(ivec,jmepe)*yv(ivec,ich3o2h)& +0.8*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +0.74*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)& +rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +1.56*rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)& +rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)& +0.629*rrv(ivec,kc76)*yv(ivec,iisop)*yv(ivec,ioh)& +0.6*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +rrv(ivec,kohch3oh)*yv(ivec,ioh)*yv(ivec,ich3oh)& +rrv(ivec,kno3mo2)*yv(ivec,ino3)*yv(ivec,ich3o2)& +0.1*rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)& +0.54*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +1.22*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +1.8*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.167*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.15*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.282*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.9*rjv(ivec,jispd)*yv(ivec,iispd)& +rrv(ivec,kaco2no)*yv(ivec,iaco2)*yv(ivec,ino)& +rrv(ivec,khypo2no)*yv(ivec,ino)*yv(ivec,ihypro2)& +rjv(ivec,jmo2no2b)*yv(ivec,ich3o2no2) XLCH2O=rjv(ivec,jach2o)+rjv(ivec,jbch2o)+yv(ivec,ioh)*rrv(ivec,kfrmoh)& +rrv(ivec,kc41)*yv(ivec,ino3) XLCH2O=XLCH2O + vdv(ivec,ich2o) yv(ivec,ich2o)=(y0v(ivec,ich2o)+PCH2O*DT)/(1.+XLCH2O*DT) ! --- CBIV-elements for higher HC-chemistry: ALD2 MGLY ! --- ETH OLE ISOP ROOH ORGNTR PALD2=0.11*rrv(ivec,kc52)*yv(ivec,ipar)*yv(ivec,ioh)& +1.1*rrv(ivec,kc53)*yv(ivec,iror)& +0.95*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +0.5*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)& +0.91*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +0.22*rrv(ivec,kc61)*yv(ivec,ieth)*yv(ivec,ioh)& +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)& +0.04*rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)& +0.991*rrv(ivec,kohc2h6)*yv(ivec,ioh)*yv(ivec,ic2h6)& +0.3*rjv(ivec,jrooh)*yv(ivec,irooh)& +0.3*rrv(ivec,kc84)*yv(ivec,ioh)*yv(ivec,iorgntr)& +0.3*rjv(ivec,jorgn)*yv(ivec,iorgntr)& +rrv(ivec,kohethoh)*yv(ivec,ioh)*yv(ivec,iethoh)& +0.5*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +0.47*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +0.21*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +0.273*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.02*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.357*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.067*rjv(ivec,jispd)*yv(ivec,iispd)& +0.7*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2)& +0.27*rrv(ivec,kc3h7o2no)*yv(ivec,ic3h7o2)*yv(ivec,ino)& +rrv(ivec,khypo2no)*yv(ivec,ihypro2)*yv(ivec,ino) XLALD2=rrv(ivec,kc43)*yv(ivec,ioh)+rrv(ivec,kc44)*yv(ivec,ino3)& +rjv(ivec,jald2) XLALD2=XLALD2 + vdv(ivec,iald2) yv(ivec,iald2)=(y0v(ivec,iald2)+PALD2*DT)/(1.+XLALD2*DT) PMGLY=0.19*rrv(ivec,kohrooh)*yv(ivec,ioh)*yv(ivec,irooh)& +0.168*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.85*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.5*rrv(ivec,kaco2mo2)*yv(ivec,iaco2)*yv(ivec,ich3o2) XLMGLY=rrv(ivec,kc73)*yv(ivec,ioh)+rjv(ivec,jmgly) yv(ivec,imgly)=(y0v(ivec,imgly)+PMGLY*DT)/(1.+XLMGLY*DT) XLETH=rrv(ivec,kc61)*yv(ivec,ioh)+rrv(ivec,kc62)*yv(ivec,io3) yv(ivec,ieth)=y0v(ivec,ieth)/(1.+XLETH*DT) POLE=0. XLOLE=rrv(ivec,kc57)*yv(ivec,ioh)+rrv(ivec,kc58)*yv(ivec,io3)+rrv(ivec,kc59)*yv(ivec,ino3) yv(ivec,iole)=(y0v(ivec,iole)+POLE*DT)/(1.+XLOLE*DT) PROOH=rrv(ivec,kc82)*yv(ivec,ixo2)*yv(ivec,iho2)& +rrv(ivec,kc85)*yv(ivec,iho2)*yv(ivec,ixo2n)& +rrv(ivec,kaco2ho2)*yv(ivec,iaco2)*yv(ivec,iho2)& +rrv(ivec,kc3h7o2ho2)*yv(ivec,ic3h7o2)*yv(ivec,iho2)& +rrv(ivec,khypo2ho2)*yv(ivec,ihypro2)*yv(ivec,iho2) XLROOH=rjv(ivec,jrooh)+rrv(ivec,kohrooh)*yv(ivec,ioh) XLROOH = XLROOH + vdv(ivec,irooh) yv(ivec,irooh)=(y0v(ivec,irooh)+PROOH*DT)/(1.+XLROOH*DT) PORGNTR=rrv(ivec,kc81)*yv(ivec,ino)*yv(ivec,ixo2n)& +0.8*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)& +rrv(ivec,kno3c3h6)*yv(ivec,ic3h6)*yv(ivec,ino3)& +0.53*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +0.85*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd) XLORGNTR=rrv(ivec,kc84)*yv(ivec,ioh)+rjv(ivec,jorgn) XLORGNTR=XLORGNTR+vdv(ivec,iorgntr) yv(ivec,iorgntr)=(y0v(ivec,iorgntr)+PORGNTR*DT)/(1.+XLORGNTR*DT) PACET=0.82*rrv(ivec,kc3h7o2no)*yv(ivec,ino)*yv(ivec,ic3h7o2) XLACET=rjv(ivec,ja_acet)+rjv(ivec,jb_acet)+rrv(ivec,kohacet)*yv(ivec,ioh) XLACET=XLACET+vdv(ivec,iacet) yv(ivec,iacet)=(y0v(ivec,iacet)+PACET*DT)/(1.+XLACET*DT) ! gas phase sulfur & ammonia qdms1=rrv(ivec,kdmsoha)*yv(ivec,ioh)+rrv(ivec,kdmsno3)*yv(ivec,ino3) qdms2=rrv(ivec,kdmsohb)*yv(ivec,ioh) qdms=qdms1+qdms2 yv(ivec,idms)=y0v(ivec,idms)/(1.+qdms*DT) pso2=yv(ivec,idms)*(qdms1+0.75*qdms2) pmsa=yv(ivec,idms)*0.25*qdms2 qso2=rrv(ivec,kso2oh)*yv(ivec,ioh) qso2d=qso2 + vdv(ivec,iso2) yv(ivec,iso2)=(y0v(ivec,iso2)+pso2*DT) /(1.+qso2d*DT) !qso2d includes deposition #ifdef with_m7 ! Dry deposition of gas-phase H2SO4 is not important and therefore neglected: yv(ivec,iso4)=(y0v(ivec,iso4)+qso2*yv(ivec,iso2)*DT) ! Do not apply dry deposition to bulk aerosols (NO3_a, NH4 and MSA). ! When M7 is active, this is done in the sedimentation routine. #ifdef with_budgets ! leave out dt to get s-1 values temp_prod_so4(ivec,1)=qso2*yv(ivec,iso2)*DT #endif #else !VH: Do apply dry deposition to SO4_A : This deposition velocity does represent aerosol deposition !VH: Use the same aerosol deposition velocity for NO3_A deposition. yv(ivec,iso4)=(y0v(ivec,iso4)+qso2*yv(ivec,iso2)*DT) /(1. + vdv(ivec,iso4)*DT) !corrected CMK qso2/qso2d ! Without M7, apply dry deposition to other bulk aerosols: yv(ivec,ino3_a)=y0v(ivec,ino3_a) /(1.+vdv(ivec,ino3_a)*DT) yv(ivec,inh4)=y0v(ivec,inh4)/(1.+vdv(ivec,inh4)*DT) yv(ivec,imsa)=(y0v(ivec,imsa)+pmsa*DT) /(1.+vdv(ivec,imsa)*DT) #endif pnh2=yv(ivec,inh2o2)*yv(ivec,ino)*rrv(ivec,knh2o2no)& +yv(ivec,inh2o2)*yv(ivec,io3)*rrv(ivec,knh2o2o3)& +yv(ivec,inh2o2)*yv(ivec,iho2)*rrv(ivec,knh2o2ho2) pnh3=rrv(ivec,Knh2ho2)*yv(ivec,inh2)*yv(ivec,iho2) dnh3=yv(ivec,ioh)*rrv(ivec,knh3oh) + vdv(ivec,inh3) yv(ivec,inh3)=(y0v(ivec,inh3)+pnh3*DT)/(1.+dnh3*DT) ppnh3=yv(ivec,ioh)*yv(ivec,inh3)*rrv(ivec,knh3oh) qnh2= rrv(ivec,knh2oh)*yv(ivec,ioh)+rrv(ivec,knh2no)*yv(ivec,ino)& +rrv(ivec,knh2no2)*yv(ivec,ino2)+rrv(ivec,knh2ho2)*yv(ivec,iho2)& +rrv(ivec,knh2o3)*yv(ivec,io3)& +rrv(ivec,knh2o2) yv(ivec,inh2)=(y0v(ivec,inh2)+ppnh3*dt+pnh2*dt)/(1.+qnh2*dt) ! ! Now nh2o2 radical ! qnh2o2=rrv(ivec,knh2o2no)*yv(ivec,ino)+rrv(ivec,knh2o2o3)*yv(ivec,io3)+rrv(ivec,knh2o2ho2)*yv(ivec,iho2) yv(ivec,inh2o2)=(y0v(ivec,inh2o2)+(rrv(ivec,knh2o3)*yv(ivec,io3)*yv(ivec,inh2)*DT))/(1.+qnh2o2*DT) end do !ivec end if if ( mod(iter,maxit) == 0 ) then ! --- Long living compounds do ivec=1,lvec yv(ivec,ich4)=y0v(ivec,ich4)/(1.+rrv(ivec,kch4oh)*yv(ivec,ioh)*DT) #ifdef with_budgets !ch4loss=(1.0-1.0/((1.+rrv(ivec,kch4oh)*yv(ivec,ioh)*DT)))*y0v(ivec,ich4) temp_loss(ivec,1)=(1.0-1.0/((1.+rrv(ivec,kch4oh)*yv(ivec,ioh)*DT)))*y0v(ivec,ich4) #endif !methane loss? PCO=yv(ivec,ich2o)*(rjv(ivec,jach2o)+rjv(ivec,jbch2o)& +yv(ivec,ioh)*rrv(ivec,kfrmoh))& +rjv(ivec,jald2)*yv(ivec,iald2)& +rjv(ivec,jmgly)*yv(ivec,imgly)& +0.62*rrv(ivec,kc57)*yv(ivec,iole)*yv(ivec,ioh)& +0.65*rrv(ivec,kc58)*yv(ivec,iole)*yv(ivec,io3)& +0.56*rrv(ivec,kc59)*yv(ivec,iole)*yv(ivec,ino3)& +0.24*rrv(ivec,kc62)*yv(ivec,ieth)*yv(ivec,io3)& +0.066*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +rrv(ivec,kc41)*yv(ivec,ich2o)*yv(ivec,ino3)& +0.56*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)& +0.47*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +0.211*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +0.47*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +0.334*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.225*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +0.643*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.333*rjv(ivec,jispd)*yv(ivec,iispd)& +rjv(ivec,ja_acet)*yv(ivec,iacet) XLCO = rrv(ivec,kcooh)*yv(ivec,ioh) XLCO = XLCO + vdv(ivec,ico) yv(ivec,ico)=(y0v(ivec,ico)+PCO*DT)/(1.+XLCO*DT) #ifdef with_budgets !carbon monoxide loss? change=(t-1)-t temp_loss(ivec,2)=y0v(ivec,ico)-(y0v(ivec,ico)+PCO*DT)/(1.+XLCO*DT) #endif pch3oh=(0.63*rrv(ivec,kmo2mo2)*yv(ivec,ich3o2)*yv(ivec,ich3o2))+& (0.5*rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)*yv(ivec,iaco2)) yv(ivec,ich3oh)=(y0v(ivec,ich3oh)+pch3oh*dt)/& (1.+(vdv(ivec,ich3oh)+rrv(ivec,kohch3oh)*yv(ivec,ioh))*dt) phcooh=(0.52*rrv(ivec,kc62)*yv(ivec,io3)*yv(ivec,ieth)) +& (0.25*rrv(ivec,ko3c3h6)*yv(ivec,io3)*yv(ivec,ic3h6)) yv(ivec,ihcooh)=(y0v(ivec,ihcooh)+(phcooh*dt))/& (1.+rrv(ivec,kohhcooh)*yv(ivec,ioh)*dt) pmcooh=0.4*rrv(ivec,kc50)*yv(ivec,ic2o3)*yv(ivec,iho2) yv(ivec,imcooh)=(y0v(ivec,imcooh)+(pmcooh*dt))/& (1.+(vdv(ivec,imcooh)+rrv(ivec,kohmcooh)*yv(ivec,ioh))*dt) yv(ivec,ic2h6)=y0v(ivec,ic2h6)/(1.+rrv(ivec,kohc2h6)*yv(ivec,ioh)*dt) yv(ivec,iethoh)=(y0v(ivec,iethoh)/(1.+(vdv(ivec,iethoh)+rrv(ivec,kohethoh)*yv(ivec,ioh))*dt)) yv(ivec,ic3h8)=(y0v(ivec,ic3h8)/(1.+rrv(ivec,kohc3h8)*yv(ivec,ioh)*dt)) xlc3h6=rrv(ivec,kohc3h6)*yv(ivec,ioh)& +rrv(ivec,ko3c3h6)*yv(ivec,io3)& +rrv(ivec,kno3c3h6)*yv(ivec,ino3) yv(ivec,ic3h6)=(y0v(ivec,ic3h6)/(1.+xlc3h6*dt)) ppar=0.35*rrv(ivec,kc77)*yv(ivec,iisop)*yv(ivec,io3)& +2.4*rrv(ivec,kc78)*yv(ivec,iisop)*yv(ivec,ino3)& +5.0*rrv(ivec,kohterp)*yv(ivec,ioh)*yv(ivec,iterp)& +6.0*rrv(ivec,ko3terp)*yv(ivec,io3)*yv(ivec,iterp)& +6.0*rrv(ivec,kno3terp)*yv(ivec,ino3)*yv(ivec,iterp)& +1.565*rrv(ivec,kohispd)*yv(ivec,ioh)*yv(ivec,iispd)& +0.36*rrv(ivec,ko3ispd)*yv(ivec,io3)*yv(ivec,iispd)& +1.282*rrv(ivec,kno3ispd)*yv(ivec,ino3)*yv(ivec,iispd)& +0.832*rjv(ivec,jispd)*yv(ivec,iispd)& +0.6*rrv(ivec,kaco2mo2)*yv(ivec,ich3o2)*yv(ivec,iaco2) xlpar=rrv(ivec,kc52)*yv(ivec,ioh)+rrv(ivec,kc83)*yv(ivec,irxpar) yv(ivec,ipar)=(y0v(ivec,ipar)+ppar*dt)/(1.+xlpar*dt) !cmk ____added rn222 chemistry in EBI language yv(ivec,irn222) = y0v(ivec,irn222)/(1.+rrv(ivec,krn222)*dt) yv(ivec,ipb210) = y0v(ivec,ipb210)+y0v(ivec,irn222)-yv(ivec,irn222) end do !ivec end if end do !ITER end subroutine DO_EBI subroutine NOYmass integer i,j,imax real :: ncormax,ncorav,totn,totn0,fnoy,fnoy1 real :: ncorr,ncorr1,ncorr2,ncorr3,ncorr4,ncorr5, totdep logical :: nerr ncormax=0. ncorav=0. nerr=.false. imax = 0 do j=j1,j2 do i=i1,i2 imax = imax + 1 ! !** Guarantee exact mass conservation of NOY ! (this may matter a few percent) ! fnoy=y(i,j,ino)+y(i,j,ino2)+y(i,j,ino3)+2.*y(i,j,in2o5)+y(i,j,ihno4) if (level == 1) then #ifndef without_dry_deposition totdep = (y(i,j,ino) *vd(region,ino )%surf(i,j) + & y(i,j,ino2)*vd(region,ino2)%surf(i,j) + & y(i,j,ino3)*vd(region,ino3)%surf(i,j) + & y(i,j,ihno3)*vd(region,ihno3)%surf(i,j) + & y(i,j,ipan)*vd(region,ipan)%surf(i,j) + & y(i,j,iorgntr)*vd(region,iorgntr)%surf(i,j) + & y(i,j,ich3o2no2)*vd(region,ich3o2no2)%surf(i,j) + & 2*y(i,j,in2o5)*vd(region,in2o5)%surf(i,j) + & y(i,j,ihno4)*vd(region,ihno4)%surf(i,j) + & y(i,j,ihono)*vd(region,ihono)%surf(i,j))*dt/ye(i,j,idz) #else totdep = 0.0 #endif else totdep = 0.0 end if totn0=y0(i,j,inox)+y0(i,j,ihno3)+y0(i,j,ipan)+y0(i,j,ihono)+ & y0(i,j,iorgntr) + y0(i,j,ich3o2no2) + ye(i,j,ieno)*dt - totdep ! note that emino is added here and the total deposition is subtracted ! ! totn0 contains all nitrogen at beginning of timestep + nox emissions ! ! ! totn contains all nitrogen at end of timestep ! totn=fnoy+y(i,j,ihno3)+y(i,j,ipan)+y(i,j,iorgntr)+y(i,j,ihono)+ y(i,j,ich3o2no2) ! correction factor for all nitrogen compounds ncorr=totn-totn0 if ( totn < tiny(totn) ) cycle if ( (abs(ncorr)/totn) > 0.05 ) then !CMK changed from 0.1 to 0.05 nerr=.true. !AJS>>> ! print *,'NOYmass: N-error....',region,level,i,j,totn0,totn ! print *,'NOYmass: emino ',ye(i,j,ieno)*dt/y0(i,j,iair)*1e9 ! print *,'NOYmass: NO(0) ', & ! y0(i,j,ino)/y0(i,j,iair)*1e9,y(i,j,ino)/y(i,j,iair)*1e9 ! print *,'NOYmass: NO2(0) ', & ! y0(i,j,ino2)/y0(i,j,iair)*1e9,y(i,j,ino2)/y(i,j,iair)*1e9 ! print *,'NOYmass: O3(0) ', & ! y0(i,j,io3)/y0(i,j,iair)*1e9,y(i,j,io3)/y(i,j,iair)*1e9 ! print *,'NOYmass: NO3(0) ', & ! y0(i,j,ino3)/y0(i,j,iair)*1e9,y(i,j,ino3)/y(i,j,iair)*1e9 ! print *,'NOYmass: N2O5(0) ', & ! y0(i,j,in2o5)/y0(i,j,iair)*1e9,y(i,j,in2o5)/y(i,j,iair)*1e9 ! print *,'NOYmass: HNO4(0) ', & ! y0(i,j,ihno4)/y0(i,j,iair)*1e9,y(i,j,ihno4)/y(i,j,iair)*1e9 ! print *,'NOYmass: HNO3(0) ', & ! y0(i,j,ihno3)/y0(i,j,iair)*1e9,y(i,j,ihno3)/y(i,j,iair)*1e9 ! print *,'NOYmass: PAN(0) ', & ! y0(i,j,ipan)/y0(i,j,iair)*1e9,y(i,j,ipan)/y(i,j,iair)*1e9 ! print *,'NOYmass: ORGNT(0) ', & ! y0(i,j,iorgntr)/y0(i,j,iair)*1e9,y(i,j,iorgntr)/y(i,j,iair)*1e9 ! print *,'NOYmass: NOx(0) ', & ! y0(i,j,inox)/y0(i,j,iair)*1e9,y(i,j,inox)/y(i,j,iair)*1e9 ! print *,'NOYmass: ',rj(i,j,jhno3),rr(i,j,kohhno3)*y(i,j,ioh), & ! y(i,j,ioh)/y(i,j,iair)*1e9 !<<5% ")'); call goPr write (gol,'(" Maximum correction : ",f8.2)') ncormax; call goPr write (gol,'(" Average correction in this loop (imax) : ",f8.2," (",i6,")")') ncorav/imax, imax; call goPr end if end subroutine NOYmass #ifdef with_budgets !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: ! ! !DESCRIPTION: increase reaction budgets for each reaction ! arrays nrr and nrj determine which species are ! involved in a reaction !\\ !\\ ! !INTERFACE: ! subroutine INCC2C3 ! #ifdef with_tendencies use TRACER_DATA, only : PLC_AddValue, plc_ipr_lddep, plc_kg_from_tm #endif ! integer, intent(out) :: status integer :: i01,n1,n2,jl,i,j ! nrj and nrr used for reaction budget calculations integer,dimension(nj),parameter :: nrj=(/io3,ino2,ih2o2,ihno3,ihno4,in2o5,in2o5,ich2o,ich2o, & ich3o2h,ino3,ino3,ipan,ipan,iorgntr,iald2,imgly,irooh,io2,iispd,iacet,iacet, & ihono,ich3o2no2,ich3o2no2/) integer,dimension(nreac,2),parameter :: nrr = reshape((/ & ino ,iho2 ,ich3o2,ino2 ,ioh ,ino2 ,ino ,ino2 ,in2o5, ihno4 , & ino2 ,ihno4 ,iair ,ih2o ,io3 ,ico ,io3 ,ih2o2 ,ich2o, ich4 , & ioh ,ioh ,ich3o2,ich3o2 ,ich3o2 ,iho2 ,iho2 ,in2o5 ,ioh , ich2o , & iald2 ,iald2 ,ic2o3 ,ic2o3 ,ipan ,ic2o3 ,ic2o3 ,ipar ,iror , iror , & ioh ,io3 ,ino3 ,ioh ,io3 ,ioh ,ioh ,io3 ,ino3 , ixo2 , & ixo2 ,ixo2n ,ixo2 ,irxpar ,iorgntr,ixo2n ,idms ,idms ,idms , iso2 , & inh3 ,inh3 ,inh2 ,inh2 ,inh2 ,inh2 ,inh2 ,inh2 ,ioh ,ioh , ino3 , & ino3 ,ino3 ,ino3 ,ioh ,ioh ,ioh ,ioh ,ioh ,io3 , ino3 , & ioh ,io3 ,ino3 ,ioh ,io3 ,ino3 ,irn222,io3 ,iair , iacet , & iaco2 ,iaco2 ,iaco2 ,iaco2 ,ixo2 ,ixo2n ,ino ,iho2 ,ino , iho2 , & in2o5 ,ino3 ,iho2 ,iho2 ,ino ,io3 ,iho2, ioh ,ioh ,ich3o2,ich3o2no2, & !second reaction partner (if monmolecular = 0) io3 ,ino ,ino ,ioh ,ihno3 ,io3 ,ino3 ,ino3 ,0 ,ioh , & iho2 ,0 ,0 ,0 ,iho2 ,ioh ,ioh ,ioh ,ioh ,ioh , & ich3o2h,irooh ,iho2 ,iho2 ,ich3o2 ,ioh ,iho2 ,0 ,0 ,ino3 , & ioh ,ino3 ,ino ,ino2 ,0 ,ic2o3 ,iho2 ,ioh ,0 ,0 , & iole ,iole ,iole ,ieth ,ieth ,imgly ,iisop ,iisop ,iisop ,ino , & ixo2 ,ino ,iho2 ,ipar ,ioh ,iho2 ,ioh ,ioh ,ino3 ,ioh , & iacid ,ioh ,ioh ,ino ,ino2 ,iho2 ,0 ,io3 ,ich3oh ,ihcooh ,iho2 , & ich3o2 ,ic2o3 ,ixo2 ,imcooh ,ic2h6 ,iethoh,ic3h8 ,ic3h6 ,ic3h6 ,ic3h6 , & iterp ,iterp ,iterp ,iispd ,iispd ,iispd ,0 ,iair ,0 ,ioh , & iho2 ,ich3o2,ino ,ixo2 ,ixo2n ,ixo2n ,ic3h7o2,ic3h7o2,ihypro2,ihypro2, & 0 ,0 ,0 ,0 ,inh2o2 ,inh2o2,inh2o2, ino ,ihono ,ino2 ,0/),(/nreac,2/)) real :: c1,xdep c1=dt*1000./xmair !conversion to moles... ! Reaction budgets do i01=1,nj !photolysis rates n1=nrj(i01) do j=j1,j2 do i=i1,i2 if(n1 > 0) cr2(i,j,i01)=cr2(i,j,i01)+rj(i,j,i01)*y(i,j,n1) end do end do end do!i01=1,nj ! do i01=1,nreac !reactions n1=nrr(i01,1) !make sure n1 > 0 n2=nrr(i01,2) if (n2 > 0.) then do j=j1,j2 do i=i1,i2 cr3(i,j,i01)= cr3(i,j,i01)+y(i,j,n1)*y(i,j,n2)*rr(i,j,i01) end do end do else do j=j1,j2 do i=i1,i2 cr3(i,j,i01)= cr3(i,j,i01)+y(i,j,n1)*rr(i,j,i01) end do end do end if end do !i01=1,nreac if ( level == 1 ) then ! deposition budget do i01=1,ntrace if (fscale(i01) > 0) then do j=j1,j2 do i=i1,i2 #ifndef without_dry_deposition if (i01 .ne. inox) then xdep = y(i,j,i01)*vd(region,i01)%surf(i,j)/ye(i,j,idz)* & c1*ye(i,j,iairm)/y(i,j,iair) !from updated concentrations else ! compute nox deposition from other contributors! xdep = (y(i,j,ino) *vd(region,ino)%surf(i,j) + & y(i,j,ino2 )*vd(region,ino2)%surf(i,j)+ & y(i,j,ino3) *vd(region,ino3)%surf(i,j)+ & y(i,j,ihno3)*vd(region,ihno3)%surf(i,j)+ & 2*y(i,j,in2o5)*vd(region,in2o5)%surf(i,j) + & y(i,j,ihono)*vd(region,ihono)%surf(i,j)) & /ye(i,j,idz)* & c1*ye(i,j,iairm)/y(i,j,iair) !from updated concentrations endif #else xdep = 0.0 #endif buddep_dat(region)%dry(i,j,i01) = & buddep_dat(region)%dry(i,j,i01) + xdep if ( i01 == 1 ) then !seperate deposition from 'other' chemistry #ifndef without_dry_deposition sum_deposition(region) = sum_deposition(region) - & xdep*ra(1)*1e-3 ! in kg #endif sum_chemistry(region) = sum_chemistry(region) + & (y(i,j,1)-y0(i,j,1))/y(i,j,iair)* & ye(i,j,iairm)/xmair*ra(1) + xdep*ra(1)*1e-3 end if ! FIXME TENDENCIES #ifdef with_tendencies ! Add deposition budget in kg to tendencies; ! (mole tm tracer) * (g/mole) * (kg/g) = kg tm tracer call PLC_AddValue( region, plc_ipr_lddep, i, j, 1, i01, & xdep * ra(i01) * 1e-3 * plc_kg_from_tm(i01), & ! kg plc tracer status ) ! IF_NOTOK_RETURN(status=1) #endif end do end do endif end do !i01 else ! other layers do j=j1,j2 do i=i1,i2 sum_chemistry(region) = sum_chemistry(region) + & (y(i,j,1)-y0(i,j,1))/y(i,j,iair)*ye(i,j,iairm)/xmair*ra(1) end do end do end if !level ==1 end subroutine INCC2C3 !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: REACBUD ! ! !DESCRIPTION: accumulate budgets, o3 P/L !\\ !\\ ! !INTERFACE: ! SUBROUTINE REACBUD ! ! !USES: ! USE BUDGET_GLOBAL, ONLY : budg_dat, nzon_vg ! !EOP !------------------------------------------------------------------------ !BOC integer :: i01, i, j, nzone, nzone_v real :: c1 c1=dt*1000./xmair !conversion to moles do j=j1,j2 do i=i1,i2 nzone=budg_dat(region)%nzong(i,j) !global budget nzone_v=nzon_vg(level) !level is passed to ebi... do i01=1,nj budrjg(nzone,nzone_v,i01)=budrjg(nzone,nzone_v,i01)+ & cr2(i,j,i01)*c1*ye(i,j,iairm)/y(i,j,iair) !units mole end do !nj ! ozone production on a 3D grid: ! defined as: HO2 + NO, CH3O2 + NO, XO2 + NO, C2O3 + NO o3p(region)%d3(i,j,level) = o3p(region)%d3(i,j,level) + & (cr3(i,j,2) + cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50)) & *c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox !AERCHEMMIP output #ifdef with_m7 AC_O3_lp(region)%prod(i,j,level,1)=AC_O3_lp(region)%prod(i,j,level,1)+& (cr3(i,j,2)+ cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50)) & *c1*ye(i,j,iairm)/y(i,j,iair) #endif o3l(region)%d3(i,j,level) = o3l(region)%d3(i,j,level) + & (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + & cr3(i,j,15) + cr3(i,j,17) + cr3(i,j,42) - cr3(i,j,43) + & cr3(i,j,45) + cr3(i,j,48) - 0.1*cr3(i,j,49) - cr3(i,j,55) + & cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) & *c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox !AERCHEMMIP output #ifdef with_m7 AC_O3_lp(region)%prod(i,j,level,2)=AC_O3_lp(region)%prod(i,j,level,2)+& (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + & cr3(i,j,15) + cr3(i,j,17) + cr3(i,j,42) - cr3(i,j,43) + & cr3(i,j,45) + cr3(i,j,48) - 0.1*cr3(i,j,49) - cr3(i,j,55) + & cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) & *c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox #endif o3l(region)%d3(i,j,level) = o3l(region)%d3(i,j,level) - & cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair) !O3 + SO2 (note negative) #ifdef with_m7 AC_O3_lp(region)%prod(i,j,level,2)=AC_O3_lp(region)%prod(i,j,level,2)-& cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair) !O3 + SO2 (note negative) ! ch4+oh loss AC_loss(region)%prod(i,j,level,iloss_ch4)= AC_loss(region)%prod(i,j,level,iloss_ch4) +& cr3(i,j,20)*c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox !co+oh loss AC_loss(region)%prod(i,j,level,iloss_co)= AC_loss(region)%prod(i,j,level,iloss_co) +& cr3(i,j,16)*c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox #endif !write(4321,*) AC_O3_lp(region)%prod(i,j,level,2) !write(4322,*) (cr3(i,j,2) + cr3(i,j,3) + cr3(i,j,33) + cr3(i,j,50)) & ! *c1*ye(i,j,iairm)/y(i,j,iair) !write(4323,*) (cr3(i,j,4)-2*cr3(i,j,7) + 2*cr3(i,j,6) + cr3(i,j,8) - cr3(i,j,9) + & ! cr3(i,j,15) + cr3(i,j,17) + cr3(i,j,42) - cr3(i,j,43) + & ! cr3(i,j,45) + cr3(i,j,48) - 0.1*cr3(i,j,49) - cr3(i,j,55) + & ! cr2(i,j,1) - cr2(i,j,4) - cr2(i,j,6) - cr2(i,j,13) - 2.*cr2(i,j,10)) & ! *c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox !write(4324,*) cr4(i,j,1)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair) !O3 + SO2 (note negative) !PLS ! ch4+oh !PLS ch4oh(region)%d3(i,j,level) = ch4oh(region)%d3(i,j,level) + & !PLS cr3(i,j,20)*c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox !PLS ! S gas phase production (OH + SO2---> SO4, OH + DMS = 0.25 MSA) !PLS so4pg(region)%d3(i,j,level) = so4pg(region)%d3(i,j,level) + & !PLS (0.25*cr3(i,j,58) + cr3(i,j,60))*c1*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox OH + SO2 !PLS ! SO4 wet production !PLS so4pa(region)%d3(i,j,level) = so4pa(region)%d3(i,j,level) - & !PLS (cr4(i,j,1)+cr4(i,j,2))*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair) ! acc. moles/gridbox note neg. do i01=1,nreac budrrg(nzone,nzone_v,i01)=budrrg(nzone,nzone_v,i01)+ & cr3(i,j,i01)*c1*ye(i,j,iairm)/y(i,j,iair) !units mole end do do i01=1,nreacw budrwg(nzone,nzone_v,i01)=budrwg(nzone,nzone_v,i01)- & cr4(i,j,i01)*(1000./xmair)*ye(i,j,iairm)/y(i,j,iair) ! mole !note: changed sign to get 'positive' budget, just a ! matter of definition, !CMK end do end do end do !sum up global budgets (done in chemistry/chemistry_done) end subroutine REACBUD !EOC #endif end subroutine EBI !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: WETS ! ! !DESCRIPTION: aqueous phase chemistry of sulfur (and other): oxidation of ! SO2 and uptake of other gases in the aqueous phase ! Method : implicit solution of oxidation of SO2 !\\ !\\ ! !INTERFACE: ! #ifdef with_budgets subroutine wetS(region, level, is, js, y0, dt, y, ye, c4, status) #else subroutine wetS(region, level, is, js, y0, dt, y, ye, status) #endif ! ! !USES: ! use Dims, only: CheckShape, idatei use global_data, only: region_dat use reaction_data, only: nreacw, ntlow, kso2hp, kso2o3 use chem_param use dims, only: im, jm use Binas, only: Avog use boundary, only: LCMIP6_CO2, co2_glob ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region !region region under operation (provides im,jm,lm via chemistry.mod) integer,intent(in) :: level, is, js ! vertical level, tile start indices real, intent(in) :: dt ! chemistry timestep real, intent(in) :: y0(is:,js:,:) ! initial concentration ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: y(is:,js:,:) ! concentrations at time is t integer,intent(out) :: status ! ! !INPUT/OUTPUT PARAMETERS: ! real, intent(inout) :: ye(is:,js:,:) ! extra fields (temp, cc, pH) #ifdef with_budgets real, intent(inout) :: c4(is:,js:,:) ! budget accumulatior #endif ! ! !REVISION HISTORY: ! ???? - Ad Jeuken (KNMI) and Frank Dentener (IMAU) - v0 ! Jan 2002 - Maarten Krol (IMAU) - adapted for TM5 ! Feb 2007 - Elisabetta Vignati (JRC) - adapted for TM5/M7 ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Wets' integer n,i,j,l,itemp,iter real x1,x2,x3,b1,b2,so2x,dh2o2,dso2,disc,dnh3,dn2o5,xso2o3a,xso2o3b,co2 real,parameter :: rg=0.08314 real,dimension(:,:),allocatable :: hkso2 ! Henry's constant for sulfur dioxide real,dimension(:,:),allocatable :: hkh2o2 ! Henry's constant for hydroperoxide real,dimension(:,:),allocatable :: hko3 ! Henry's constant for ozone real,dimension(:,:),allocatable :: dkso2 ! Dissociation constant for SO2 real,dimension(:,:),allocatable :: dkhso3 ! Dissociation constant for HSO3- real,dimension(:,:),allocatable :: dkh2o ! dissociation constant water real,dimension(:,:),allocatable :: dknh3 ! dissociation constant ammonia real,dimension(:,:),allocatable :: hknh3 ! Henry's constant ammonia real,dimension(:,:),allocatable :: hkco2 ! Henry's constant CO2 real,dimension(:,:),allocatable :: dkco2 ! Dissociation constant CO2 real phs4 ! effective dissolvation of S(IV) real phso2 ! effective dissolvation of SO2 real phh2o2 ! effective dissolvation of H2O2 real phozone ! effective dissolvation of O3 real,dimension(:,:),allocatable :: hplus !concentration h+ real,dimension(:,:),allocatable :: sulph !accumul+coarse mode sulphate real a1,a2,a,b,c,z,ft ! help variables real xcov,xliq,xl,temp,rt,ztr,h2o,air,press ! meteo real,dimension(:,:,:),allocatable :: rw ! reaction rates logical,dimension(:,:),allocatable :: cloudy integer :: i1, i2, j1, j2 real :: l_sum_wet ! --- begin -------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! check arguments ... call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y0), status ) IF_NOTOK_RETURN(status=1) call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status ) IF_NOTOK_RETURN(status=1) call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status ) IF_NOTOK_RETURN(status=1) #ifdef with_budgets call CheckShape( (/i2-i1+1, j2-j1+1, nreacw /), shape(c4), status ) IF_NOTOK_RETURN(status=1) #endif allocate(hkso2 (i1:i2, j1:j2)) allocate(hkh2o2 (i1:i2, j1:j2)) allocate(hko3 (i1:i2, j1:j2)) allocate(dkso2 (i1:i2, j1:j2)) allocate(dkhso3 (i1:i2, j1:j2)) allocate(dkh2o (i1:i2, j1:j2)) allocate(dknh3 (i1:i2, j1:j2)) allocate(hknh3 (i1:i2, j1:j2)) allocate(hkco2 (i1:i2, j1:j2)) allocate(dkco2 (i1:i2, j1:j2)) allocate(hplus (i1:i2, j1:j2)) allocate(rw (i1:i2, j1:j2, nreacw)) allocate(cloudy (i1:i2, j1:j2)) allocate(sulph (i1:i2, j1:j2)) !----------------------------- ! wet phase reactions !----------------------------- rw =0.0 hplus=0.0 if (LCMIP6_CO2) then co2=co2_glob else ! ! JEW: now scaled to 2000 to account for annual growth of ~2ppbv/yr-1 ! co2=3.69e-4 ! was parameter co2=3.75e-4, endif #if defined (with_budgets) l_sum_wet = 0.0 #endif do j = j1, j2 do i = i1, i2 cloudy(i,j)=.false. ! lwc is dimensionless if ((ye(i,j,ilwc) > 1e-10).and.(ye(i,j,icc) > 0.02)) then cloudy(i,j)=.true. TEMP=ye(i,j,i_temp) ZTR=(1./TEMP-1./298) RT=TEMP*rg ITEMP=nint(TEMP-float(ntlow)) ! !CEV sulph is the initial total sulphate content in accumulation+ !coarse mode, the incloud production is calculated on bulk !characteristics, and redistributed on the modes depending on their !particle numbers #ifdef with_m7 !Stelios: small contributions from nucleation and Aitken mode ! as well as gas-phase should be added for pH calculation sulph(i,j)=y0(i,j,iso4acs)+y0(i,j,iso4cos)+& y0(i,j,iso4nus)+y0(i,j,iso4ais)+y0(i,j,iso4) #else sulph(i,j)=y0(i,j,iso4) #endif ! ! Henry and dissociation equilibria ! dkh2o(i,j) =1.01e-14*exp(-6706.0 *ztr) !h2o<=>hplus+so3-- !bug hkco2(i,j) =3.4e-2*(2420.*ztr) ! is already dimensionless hkco2(i,j) =3.4e-2*exp(2420.*ztr) ! is already dimensionless dkco2(i,j) =4.5E-7*exp(-1000.*ztr) !co2aq<=>hco3- + hplus hkso2(i,j) =henry(iso2,itemp)*rt !dimensionless dknh3(i,j) =1.8e-5*exp(-450.*ztr) !nh3<=>nh4+ + OH- hknh3(i,j) =henry(inh3,itemp)*rt !dimensionless hkh2o2(i,j)=henry(ih2o2,itemp)*rt !dimensionless hko3(i,j) =henry(io3,itemp)*rt !dimensionless dkso2(i,j) =1.7e-2*exp(2090.*ztr) !so2<=>hso3m+hplus dkhso3(i,j)=6.6e-8*exp(1510.*ztr) !hso3m<=>so3-- + hplus ! ! calculate H+ from initial sulfate, ammonium, hno3, and nh3 ! if solution is strongly acidic no further calculations are performed ! xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc) !x1 is initial strong acidity from SO4 and NO3 ! !acidity from strong acids alone ! !CMK hplus(i,j)=(2.*y0(i,j,iso4)+y0(i,j,imsa)-y0(i,j,inh4)+ & !CMK y0(i,j,ihno3)+y0(i,j,ino3_a))/xl hplus(i,j)=(2.*sulph(i,j) + & y0(i,j,imsa)-y0(i,j,inh4)+ & y0(i,j,ihno3)+y0(i,j,ino3_a))/xl end if end do end do do iter=1,10 do j=j1, j2 do i=i1, i2 ! only if solution pH>4.5 if ( cloudy(i,j) .and. hplus(i,j) < 3e-5 ) then xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc) !CEV y0(i,j,iso4)---> sulph(i,j) x1=(2.*sulph(i,j)+y0(i,j,imsa)+y0(i,j,ihno3)+ & y0(i,j,ino3_a))/xl !x2 is initial total NHx x2=(y0(i,j,inh3)+y0(i,j,inh4))/xl !x3 is combined dissolution and solubility const for CO2 x3=dkco2(i,j)*hkco2(i,j)*co2 a1=dkh2o(i,j)/dknh3(i,j)*(1.+1./hknh3(i,j)) ! integration constant a2=y0(i,j,iso2)/xl !initial SO2 ! trap division by zero ... if ( hplus(i,j) == 0.0 ) then z = 0.0 else z = a2/( hplus(i,j)/dkso2(i,j)*(1.0+1.0/hkso2(i,j)) + dkhso3(i,j)/hplus(i,j) + 1.0 ) end if ! solve quadratic equation for new H+ concentration: a=1.+x2/(a1+hplus(i,j)) b=-x1-z c=-x3-2.*dkhso3(i,j)*z z=max(0.,(b*b-4.*a*c)) hplus(i,j)=max(1.e-10,(-b+sqrt(z))/(2.*a)) end if end do ! end do ! i,j loop end do !iter do j=j1,j2 do i=i1,i2 if (cloudy(i,j)) then temp=ye(i,j,i_temp) ZTR=(1./TEMP-1./298) xliq=ye(i,j,ilwc)/ye(i,j,icc) xl=ye(i,j,ilwc)*Avog*1e-3/ye(i,j,icc) ye(i,j,iph)=-log10(hplus(i,j)) ! pH for diagnostics ! phase factor ratio of aqueous phase to gas phase concentration phs4 =hkso2(i,j) *(1.+dkso2(i,j)/hplus(i,j)+ & dkso2(i,j)*dkhso3(i,j)/hplus(i,j)/hplus(i,j))*xliq phso2 =hkso2(i,j) *xliq phh2o2 =hkh2o2(i,j)*xliq phozone=hko3(i,j) *xliq ! the original rate equations could be partly in look-up table rw(i,j,KSO2HP) =8e4*exp(-3560.*ztr)/(0.1+hplus(i,j)) XSO2O3A=4.39e11*exp(-4131./temp)+2.56e3*exp(-966./temp) !S(IV) XSO2O3B=2.56e3*exp(-966./temp)/hplus(i,j) !divide by [H+]!S(IV) ! make rate constants dimensionless by multiplying ! by (1./xliq/avo=6e20) ! multiply with the fractions of concentrations residing ! in the aqueous phase rw(i,j,KSO2HP)=rw(i,j,KSO2HP)/xl*phso2/(1.+phs4)*phh2o2/(1.+phh2o2) rw(i,j,KSO2O3)=(XSO2O3A+XSO2O3B)/xl*phs4/(1.+phs4)*phozone/ & (1.+phozone) end if !cloudy end do ! end do ! I,J, LOOP !------------- Start main loop do j=j1,j2 do i=i1,i2 ! ! only cloud chemistry if substantial amount of clouds are present ! if (cloudy(i,j)) then ! ! oxidation of S(IV) by O3 ! so2x=y0(i,j,iso2) xcov=ye(i,j,icc) x1=min(100.,rw(i,j,kso2o3)*y0(i,j,io3)*dt) dso2=y0(i,j,iso2)*xcov*(exp(-x1)-1.) !only applied to xcov part of cloud !CMK print *, i,j, xcov, x1, y0(i,j,iso2), dso2 dso2=max(-y0(i,j,io3)*xcov,dso2)! limit to o3 availability !CEV y(i,j,iso2)=y0(i,j,iso2)+dso2 !NOTE CMK: parallel MPI should take care here! #ifdef with_m7 ft = y0(i,j,iacs_n) + y0(i,j,icos_n) if (ft > tiny(ft)) then y(i,j,iso4acs)=y0(i,j,iso4acs)-dso2*(y0(i,j,iacs_n)/ft) y(i,j,iso4cos)=y0(i,j,iso4cos)-dso2*(y0(i,j,icos_n)/ft) y(i,j,iso2)=y0(i,j,iso2)+dso2 y(i,j,io3)=y0(i,j,io3)+dso2 !AERHEMMIP ! Production of liquid phase so4 ! D_prod_liq(i,j)=-dso2 #ifdef with_budgets c4(i,j,1)=c4(i,j,1)+dso2 !conversion 1e-3 g->kg 1e6 cm-3 ->1m-3 diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4 AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)=AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4!kg #endif else #ifdef with_budgets diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)+0.0 AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)=AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)+0.0 #endif !CEV y(i,j,iso4)=y0(i,j,iso4)-dso2 y(i,j,iso4acs)=y0(i,j,iso4acs) y(i,j,iso4cos)=y0(i,j,iso4cos) y(i,j,iso2)=y0(i,j,iso2) y(i,j,io3)=y0(i,j,io3) endif !CEV y(i,j,io3)=y0(i,j,io3)+dso2 #else ! gas phase chemistry: ft = 0. y(i,j,iso4)=y0(i,j,iso4)-dso2 y(i,j,iso2)=y0(i,j,iso2)+dso2 y(i,j,io3)=y0(i,j,io3)+dso2 #ifdef with_budgets c4(i,j,1)=c4(i,j,1)+dso2 #endif #endif #ifdef with_budgets if ( io3 == 1 ) l_sum_wet = l_sum_wet- & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) if ( iso2 == 1 ) l_sum_wet = l_sum_wet+ & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) if ( iso4 == 1 ) l_sum_wet = l_sum_wet- & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) !CEV c4(i,j,1)=c4(i,j,1)+dso2 #endif xliq=ye(i,j,ilwc)/ye(i,j,icc) ! ! oxidation of S(IV) by H2O2 ! !*** here we explicitly solve the dv: ! y'=P-Q*y-R*y*y (P and Q are 0=>b3=0.) ! so2x=y(i,j,iso2) if ( so2x > tiny(so2x) ) then b1=rw(i,j,kso2hp) b2=b1*(y0(i,j,ih2o2)-so2x) disc=min(100.,sqrt(b2*b2)) ! disc is b2 for b3=0.0 x1=(b2-disc)/(-2.*b1) ! in this case x1 =0. x2=(b2+disc)/(-2.*b1) x3=(so2x-x1)/(so2x-x2)*exp(-disc*dt) so2x=(x1-x2*x3)/(1.-x3) dso2=(so2x-y(i,j,iso2))*xcov dso2=max(dso2,-y0(i,j,ih2o2)*xcov) !CEV y(i,j,iso2) =y(i,j,iso2)+dso2 ! dso2 is loss of SO2 and H2O2 ! divide produced SO4 over CO/ACC #ifdef with_m7 ft = y0(i,j,iacs_n) + y0(i,j,icos_n) if (ft > tiny(ft)) then y(i,j,iso4acs)=y(i,j,iso4acs)-dso2*(y0(i,j,iacs_n)/ft) y(i,j,iso4cos)=y(i,j,iso4cos)-dso2*(y0(i,j,icos_n)/ft) y(i,j,iso2)=y(i,j,iso2)+dso2 y(i,j,ih2o2)=y0(i,j,ih2o2)+dso2 #ifdef with_budgets c4(i,j,2)=c4(i,j,2)+dso2 ! add amount liquid so4 production to diagnostic and change to molec->mass kg ! diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4!kg AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)=AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)-dso2/y(i,j,iair)*ye(i,j,iairm)/xmair*xmso4!kg #endif else #ifdef with_budgets diag_prod(region)%prod(i,j,level,2)=diag_prod(region)%prod(i,j,level,2)+0.0 AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)=AC_diag_prod(region)%prod(i,j,level,iprod_aqso4)+0.0 #endif y(i,j,ih2o2) =y0(i,j,ih2o2) endif #else ! gas - phase chemistry y(i,j,iso4)=y(i,j,iso4)-dso2 y(i,j,iso2)=y(i,j,iso2)+dso2 y(i,j,ih2o2)=y0(i,j,ih2o2)+dso2 #ifdef with_budgets c4(i,j,2)=c4(i,j,2)+dso2 #endif #endif #ifdef with_budgets if ( ih2o2 == 1 ) l_sum_wet = l_sum_wet- & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) if ( iso2 == 1 ) l_sum_wet = l_sum_wet- & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) if ( iso4 == 1 ) l_sum_wet = l_sum_wet+ & dso2 *ye(i,j,iairm)/ (fscale(1)*y(i,j,iair)) !CEV c4(i,j,2)=c4(i,j,2)+dso2 #endif end if ! ! NH3 uptake in cloud droplets is limited by SO4 availability ! no HNO3 is considered at this point ! assume instantaneous uptake of NH3 incloud only in cloudy part ! ! EQSAM gives hell to any NH3/NH4 interchange. It first drops both in one heap ! and then redistributes. So this cloud chemistry reaction does not matter. end if !cloudy end do ! i,j,loop end do ! #ifdef with_budgets sum_wetchem(region) = sum_wetchem(region) + l_sum_wet #endif !free memory deallocate(hkso2 ) deallocate(hkh2o2 ) deallocate(hko3 ) deallocate(dkso2 ) deallocate(dkhso3 ) deallocate(dkh2o ) deallocate(dknh3 ) deallocate(hknh3 ) deallocate(hkco2 ) deallocate(dkco2 ) deallocate(hplus ) deallocate(rw ) deallocate(cloudy ) deallocate(sulph ) status = 0 end subroutine WETS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: MARK_TRAC ! ! !DESCRIPTION: calculate nox/pan/orgn/hno3 analogous to ebi scheme ! ozone production from marked nox ! simple nhx chemistry, scaled to real nhx !\\ !\\ ! !INTERFACE: ! subroutine MARK_TRAC( region, level, is, js, y, rr, rj, dt, ye) ! ! !USES: ! use chem_param use Dims, only : CheckShape use global_data, only : region_dat use dims, only : at, bt, im, jm #ifdef with_budgets use budget_global, only : budg_dat, nzon_vg #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region, level, is, js real :: dt ! time step ! ! !INPUT/OUTPUT PARAMETERS: ! real, intent(inout):: y (is:,js:,:) ! concentrations real, intent(in) :: rr(is:,js:,:) ! reaction rates real, intent(in) :: rj(is:,js:,:) ! photolysis rates real, intent(in) :: ye(is:,js:,:) ! help fields ( air masses ) ! ! !REVISION HISTORY: ! Jan 1999 - fjd ! Jan 2002 - MK - adapted for TM5 ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC integer :: status, i1, i2, j1, j2 ! --- begin -------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! check arguments ... call CheckShape( (/i2-i1+1, j2-j1+1, maxtrace/), shape(y ), status ) call CheckShape( (/i2-i1+1, j2-j1+1, nreac /), shape(rr), status ) call CheckShape( (/i2-i1+1, j2-j1+1, nj /), shape(rj), status ) call CheckShape( (/i2-i1+1, j2-j1+1, n_extra /), shape(ye), status ) call MARK_O3S ! ! more marked tracers possible here ! contains subroutine mark_o3s !--------------------------------------------------- ! marked tracer O3S stratospheric ozone !--------------------------------------------------- #ifndef without_dry_deposition use dry_deposition, only: vd #endif integer :: i, j, nzone, nzone_v real :: p3, xl3, o3old do j = j1, j2 do i = i1, i2 if (at(level+1)+bt(level+1)*1e5<= 14000) then ! ! well, you want to count all layers below 140 hPa ! (given surface pressure of 1e5 Pa) ! in the current model setup (25 layers) this means ! 12077 + 1e5*0.00181 = 12258 Pa and above... ! p3: production of o3 in stratosphere P3 = rj(i,j,jano3)*y(i,j,ino3)+ & rj(i,j,jno2)*y(i,j,ino2) XL3= rr(i,j,ko3ho2)*y(i,j,iho2)+& rr(i,j,ko3oh)*y(i,j,ioh)+ & rr(i,j,kno2o3)*y(i,j,ino2)+& rj(i,j,jo3d)+& rr(i,j,knoo3)*y(i,j,ino)+& rr(i,j,kc62)*y(i,j,ieth)+& rr(i,j,kc58)*y(i,j,iole)+& rr(i,j,kc77)*y(i,j,iisop)+& rr(i,j,ko3c3h6)*y(i,j,ic3h6)+& rr(i,j,ko3terp)*y(i,j,iterp)+& rr(i,j,ko3ispd)*y(i,j,iispd) else ! ! these are only the net destruction reactions ! P3 = 0. XL3= rr(i,j,ko3ho2)*y(i,j,iho2)+& rr(i,j,ko3oh)*y(i,j,ioh)+& rj(i,j,jo3d)+& rr(i,j,kc62)*y(i,j,ieth)+& rr(i,j,kc58)*y(i,j,iole)+& rr(i,j,kc77)*y(i,j,iisop)+& rr(i,j,ko3c3h6)*y(i,j,ic3h6)+& rr(i,j,ko3terp)*y(i,j,iterp)+& rr(i,j,ko3ispd)*y(i,j,iispd)+& rr(i,j,knh2o3)*y(i,j,inh2)+& rr(i,j,knh2o2o3)*y(i,j,inh2o2) ! add up deposition.... #ifndef without_dry_deposition if ( level == 1 ) & XL3 = XL3 + vd(region,io3)%surf(i,j)/ye(i,j,idz) #endif end if o3old=y(i,j,io3s) y(i,j,io3s)=(o3old+p3*dt)/(1.+xl3*dt) #ifdef with_budgets nzone=budg_dat(region)%nzong(i,j) ! global budget nzone_v=nzon_vg(level) budmarkg(nzone,nzone_v,1)=budmarkg(nzone,nzone_v,1)+ & (y(i,j,io3s)-o3old)*ye(i,j,iairm)*1000./xmair/y(i,j,iair) #endif end do end do !i,j end subroutine MARK_O3S end subroutine MARK_TRAC !EOC end module EBISCHEME