! ! $Id: so2_to_h2so4.f90 5605 2025-04-07 14:05:50Z lebasn $ ! SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) USE dimphy, ONLY : klon,klev USE aerophys USE infotrac_phy USE yomcst_mod_h, ONLY : RG, RD ! lifetime (sec) et O3_clim (VMR) USE phys_local_var_mod, ONLY : SO2_lifetime,H2SO4_lifetime,O3_clim,budg_3D_so2_to_h2so4,budg_so2_to_h2so4,SO2_chlm USE strataer_local_var_mod, ONLY : flag_OH_reduced, flag_H2SO4_photolysis, flag_min_rreduce IMPLICIT NONE !-------------------------------------------------------- ! transfer variables when calling this routine REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato ! stratospheric flag ! local variables in coagulation routine INTEGER :: i,j,k,nb,ilon,ilev REAL :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh REAL :: rmairdens,rrak0,rrak1,rrate,rreduce !--convert SO2 to H2SO4 budg_3D_so2_to_h2so4(:,:)=0.0 budg_so2_to_h2so4(:)=0.0 SO2_chlm(:,:)=0.0 DO ilon=1, klon DO ilev=1, klev !only in the stratosphere IF (is_strato(ilon,ilev)) THEN IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN IF (flag_OH_reduced) THEN !--convert SO2 to H2SO4 (slimane) ! Reduce OH because of SO2 ! d[OH]/dt = k1[HO2][O] +k2[HO2][O3] +k3[HO2][NO] ! − (k4[O] +k5[O3] +k6[CO] +ks[SO2])[OH] ! In steady-state: d[OH]/dt = 0, so the ratio [OH]/[HO2] = ! (k1[O] +k2[O3] +k3[NO]) / (k4[O] +k5[O3] +k6[CO] +ks[SO2]) ! In the lower and middle stratosphere, ! [OH]/[HO2] ~(k2[O3] +k3[NO]) / (k5[O3] +k6[CO] +ks[SO2]) ! CO is only important near the tropopause. NO should not dominate over O3. ! So when consideringthe big terms: [OH]/[HO2] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) ! For c: control run c (no effect of SO2), [OHc]/[HO2c] ~ ~ k2[O3] / k5[O3] = 1/ Rc ! For p: perturbed run (high SO2), [OHp]/[HO2p] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) =1/Rp ! If HOx = OH + HO2 = constante, ! OHc + HO2c = OHp + HO2p ! OHc (1+Rc) = OHp (1+Rp) ! OHp = OHc (1+Rc) / (1+Rp) ! 1+Rc = (k2[O3] +k5[O3] ) / k2[O3] ! 1+Rp = (k2[O3] +k5[O3] +ks[SO2]) / k2[O3] ! (1+Rc) / (1+Rp)= (k2[O3] +k5[O3] )/ (k2[O3] +k5[O3] +ks[SO2]) ! OHp = OHc * (k(HO2+OH)*[O3] +k(OH+O3)[O3] ) / ! (k(HO2+OH)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) ! Final: OHp = OHc * (k(HO2+O3)*[O3] +k(OH+O3)[O3] ) / ! (k(HO2+O3)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) ! latest jpl 2020 rkho2o3 = 1.0e-14*exp( -490./t_seri(ilon,ilev) ) rkoho3 = 1.7e-12*exp( -940./t_seri(ilon,ilev) ) ! rkso2oh= so2 + oh + m ->hso3 + m-> h2so4 + m ! air density M(molecule/cm3) = Pa/(K*1.38e-19) rmairdens = pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19) ! JPL 2020 rrak0 = 3.3e-31*( (300./t_seri(ilon,ilev))**4.3 ) rrak1 = 1.6e-12 rrate = (rrak0*rmairdens) / ( 1. + rrak0*rmairdens/rrak1 ) rreduce = 1./( 1. + log10((rrak0*rmairdens)/rrak1)**2 ) rkso2oh = rrate*(0.6**rreduce) ! O3 (molec/cm3): O3_clim in vmr rrak0 = O3_clim(ilon,ilev)*rmairdens ! SO2 (molec/cm3): convert from kg/kgA rrak1 = tr_seri(ilon,ilev,id_SO2_strat) & & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mSO2mol IF (rrak1 .GE. 0.0) THEN rreduce =( (rkho2o3+rkoho3)*rrak0 + rkso2oh*rrak1 ) / & ( (rkho2o3+rkoho3)*rrak0 ) ! reduce OH, so extend SO2 lifetime rreduce = rreduce*SO2_lifetime(ilon,ilev) ELSE rreduce = SO2_lifetime(ilon,ilev) ENDIF ! end of IF (rrak1) THEN ELSE rreduce = SO2_lifetime(ilon,ilev) ENDIF ! end of IF (flag_OH_reduced) THEN ! Check lifetime rreduce < timestep*1.5 (such as SO2 loss > 0.5*SO2) with exp(-1/1.5)=0.52 ! Check lifetime rreduce < timestep*3 (such as SO2 loss > 0.28*SO2) with exp(-1/3)=0.72 IF(flag_min_rreduce) THEN IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys ENDIF budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce)) tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev) tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev) ENDIF ! IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN SO2_chlm(ilon,ilev) = tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce)) * & pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19) / pdtphys !SO2 loss rate [mole.cm-3.s-1] IF (flag_H2SO4_photolysis) THEN !--convert H2SO4 to SO2 using HCl photolysis (Rinsland et al, 1995) ! 4.6e-8 sec-1: J strato taken from Feierabend, et al., Chem. Phys. Lett., 2006. ! error before: ...exp(-pdtphys/4.6e-8) so negligible photolysis, ! should have been exp(-pdtphys*4.6e-8) ! Also: Lane and Kjaergaard: ,J. Phys. Chem. A, 2008 ! We find that the dominant photodissociation mechanism of H2SO4 ! below 70 km is absorption in the visible region by OH stretching overtone ! transitions, whereas above 70 km, absorption of Lyman-R ... ! In 2D model, 0.3 Rinsland et al, 1995 but Burkholder, Mills and McKeen say 0.016 ! JH2SO4=JHCL*0.3 (*0.016 would be too slow) ! test: quick sequential integration for SO2 to H2SO4 and reverse IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN rreduce = H2SO4_lifetime(ilon,ilev) dummyso4toso2 = 0.0 ! Check lifetime rreduce < timestep*1.5 (such as H2SO4 loss > 0.5*H2SO4) with exp(-1/1.5)=0.52 ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72 IF(flag_min_rreduce) THEN IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys ENDIF dummyso4toso2 = (mSO2mol/mH2SO4mol)*tr_seri(ilon,ilev,id_H2SO4_strat)*(1.0-exp(-pdtphys/rreduce)) budg_3D_so2_to_h2so4(ilon,ilev) = budg_3D_so2_to_h2so4(ilon,ilev) + dummyso4toso2 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + dummyso4toso2 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) - (mH2SO4mol/mSO2mol)*dummyso4toso2 ENDIF ! IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN ENDIF ! end of IF (flag_H2SO4_photolysis) THEN !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol* & (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev) ELSE ! troposphere ! SO2 tropospheric lifetime (in sec) set to 5 days rreduce = 5.0*24.0*60.0*60.0 rrate =tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce)) tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - rrate SO2_chlm(ilon,ilev) = rrate * & pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19) / pdtphys !SO2 loss rate [moleccm-3s-1] ENDIF ! IF (is_strato(ilon,ilev)) THEN ENDDO ! DO (klev) ENDDO ! DO (klon) END SUBROUTINE SO2_TO_H2SO4