1 | ! |
---|
2 | ! $Id: so2_to_h2so4.F90 4743 2023-10-30 16:37:00Z dcugnet $ |
---|
3 | ! |
---|
4 | SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
5 | |
---|
6 | USE dimphy, ONLY : klon,klev |
---|
7 | USE aerophys |
---|
8 | USE infotrac_phy |
---|
9 | USE YOMCST, ONLY : RG, RD |
---|
10 | ! lifetime (sec) et O3_clim (VMR) |
---|
11 | USE phys_local_var_mod, ONLY : SO2_lifetime, H2SO4_lifetime, O3_clim, budg_3D_so2_to_h2so4, budg_so2_to_h2so4 |
---|
12 | USE strataer_local_var_mod, ONLY : flag_OH_reduced, flag_H2SO4_photolysis, flag_min_rreduce |
---|
13 | |
---|
14 | IMPLICIT NONE |
---|
15 | |
---|
16 | !-------------------------------------------------------- |
---|
17 | |
---|
18 | ! transfer variables when calling this routine |
---|
19 | REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
---|
20 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
---|
21 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
---|
22 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
---|
23 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
---|
24 | LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato ! stratospheric flag |
---|
25 | |
---|
26 | ! local variables in coagulation routine |
---|
27 | INTEGER :: i,j,k,nb,ilon,ilev |
---|
28 | REAL :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh |
---|
29 | REAL :: rmairdens,rrak0,rrak1,rrate,rreduce |
---|
30 | |
---|
31 | !--convert SO2 to H2SO4 |
---|
32 | budg_3D_so2_to_h2so4(:,:)=0.0 |
---|
33 | budg_so2_to_h2so4(:)=0.0 |
---|
34 | |
---|
35 | DO ilon=1, klon |
---|
36 | DO ilev=1, klev |
---|
37 | !only in the stratosphere |
---|
38 | IF (is_strato(ilon,ilev)) THEN |
---|
39 | IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN |
---|
40 | IF (flag_OH_reduced) THEN |
---|
41 | !--convert SO2 to H2SO4 (slimane) |
---|
42 | ! Reduce OH because of SO2 |
---|
43 | ! d[OH]/dt = k1[HO2][O] +k2[HO2][O3] +k3[HO2][NO] |
---|
44 | ! − (k4[O] +k5[O3] +k6[CO] +ks[SO2])[OH] |
---|
45 | ! In steady-state: d[OH]/dt = 0, so the ratio [OH]/[HO2] = |
---|
46 | ! (k1[O] +k2[O3] +k3[NO]) / (k4[O] +k5[O3] +k6[CO] +ks[SO2]) |
---|
47 | ! In the lower and middle stratosphere, |
---|
48 | ! [OH]/[HO2] ~(k2[O3] +k3[NO]) / (k5[O3] +k6[CO] +ks[SO2]) |
---|
49 | ! CO is only important near the tropopause. NO should not dominate over O3. |
---|
50 | ! So when consideringthe big terms: [OH]/[HO2] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) |
---|
51 | ! For c: control run c (no effect of SO2), [OHc]/[HO2c] ~ ~ k2[O3] / k5[O3] = 1/ Rc |
---|
52 | ! For p: perturbed run (high SO2), [OHp]/[HO2p] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) =1/Rp |
---|
53 | ! If HOx = OH + HO2 = constante, |
---|
54 | ! OHc + HO2c = OHp + HO2p |
---|
55 | ! OHc (1+Rc) = OHp (1+Rp) |
---|
56 | ! OHp = OHc (1+Rc) / (1+Rp) |
---|
57 | ! 1+Rc = (k2[O3] +k5[O3] ) / k2[O3] |
---|
58 | ! 1+Rp = (k2[O3] +k5[O3] +ks[SO2]) / k2[O3] |
---|
59 | ! (1+Rc) / (1+Rp)= (k2[O3] +k5[O3] )/ (k2[O3] +k5[O3] +ks[SO2]) |
---|
60 | ! OHp = OHc * (k(HO2+OH)*[O3] +k(OH+O3)[O3] ) / |
---|
61 | ! (k(HO2+OH)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) |
---|
62 | ! Final: OHp = OHc * (k(HO2+O3)*[O3] +k(OH+O3)[O3] ) / |
---|
63 | ! (k(HO2+O3)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) |
---|
64 | |
---|
65 | ! latest jpl 2020 |
---|
66 | rkho2o3 = 1.0e-14*exp( -490./t_seri(ilon,ilev) ) |
---|
67 | rkoho3 = 1.7e-12*exp( -940./t_seri(ilon,ilev) ) |
---|
68 | |
---|
69 | ! rkso2oh= so2 + oh + m ->hso3 + m-> h2so4 + m |
---|
70 | ! air density M(molecule/cm3) = Pa/(K*1.38e-19) |
---|
71 | rmairdens = pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19) |
---|
72 | ! JPL 2020 |
---|
73 | rrak0 = 3.3e-31*( (300./t_seri(ilon,ilev))**4.3 ) |
---|
74 | rrak1 = 1.6e-12 |
---|
75 | rrate = (rrak0*rmairdens) / ( 1. + rrak0*rmairdens/rrak1 ) |
---|
76 | rreduce = 1./( 1. + log10((rrak0*rmairdens)/rrak1)**2 ) |
---|
77 | rkso2oh = rrate*(0.6**rreduce) |
---|
78 | |
---|
79 | ! O3 (molec/cm3): O3_clim in vmr |
---|
80 | rrak0 = O3_clim(ilon,ilev)*rmairdens |
---|
81 | ! SO2 (molec/cm3): convert from kg/kgA |
---|
82 | rrak1 = tr_seri(ilon,ilev,id_SO2_strat) & |
---|
83 | & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mSO2mol |
---|
84 | |
---|
85 | IF (rrak1 .GE. 0.0) THEN |
---|
86 | rreduce =( (rkho2o3+rkoho3)*rrak0 + rkso2oh*rrak1 ) / & |
---|
87 | ( (rkho2o3+rkoho3)*rrak0 ) |
---|
88 | ! reduce OH, so extend SO2 lifetime |
---|
89 | rreduce = rreduce*SO2_lifetime(ilon,ilev) |
---|
90 | ELSE |
---|
91 | rreduce = SO2_lifetime(ilon,ilev) |
---|
92 | ENDIF |
---|
93 | ! end of IF (rrak1) THEN |
---|
94 | ELSE |
---|
95 | rreduce = SO2_lifetime(ilon,ilev) |
---|
96 | ENDIF |
---|
97 | ! end of IF (flag_OH_reduced) THEN |
---|
98 | |
---|
99 | ! Check lifetime rreduce < timestep*1.5 (such as SO2 loss > 0.5*SO2) with exp(-1/1.5)=0.52 |
---|
100 | ! Check lifetime rreduce < timestep*3 (such as SO2 loss > 0.28*SO2) with exp(-1/3)=0.72 |
---|
101 | IF(flag_min_rreduce) THEN |
---|
102 | IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys |
---|
103 | ENDIF |
---|
104 | budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce)) |
---|
105 | tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev) |
---|
106 | tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev) |
---|
107 | ENDIF |
---|
108 | ! IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN |
---|
109 | |
---|
110 | |
---|
111 | IF (flag_H2SO4_photolysis) THEN |
---|
112 | !--convert H2SO4 to SO2 using HCl photolysis (Rinsland et al, 1995) |
---|
113 | ! 4.6e-8 sec-1: J strato taken from Feierabend, et al., Chem. Phys. Lett., 2006. |
---|
114 | ! error before: ...exp(-pdtphys/4.6e-8) so negligible photolysis, |
---|
115 | ! should have been exp(-pdtphys*4.6e-8) |
---|
116 | ! Also: Lane and Kjaergaard: ,J. Phys. Chem. A, 2008 |
---|
117 | ! We find that the dominant photodissociation mechanism of H2SO4 |
---|
118 | ! below 70 km is absorption in the visible region by OH stretching overtone |
---|
119 | ! transitions, whereas above 70 km, absorption of Lyman-R ... |
---|
120 | ! In 2D model, 0.3 Rinsland et al, 1995 but Burkholder, Mills and McKeen say 0.016 |
---|
121 | ! JH2SO4=JHCL*0.3 (*0.016 would be too slow) |
---|
122 | ! test: quick sequential integration for SO2 to H2SO4 and reverse |
---|
123 | |
---|
124 | IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN |
---|
125 | |
---|
126 | rreduce = H2SO4_lifetime(ilon,ilev) |
---|
127 | dummyso4toso2 = 0.0 |
---|
128 | |
---|
129 | ! Check lifetime rreduce < timestep*1.5 (such as H2SO4 loss > 0.5*H2SO4) with exp(-1/1.5)=0.52 |
---|
130 | ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72 |
---|
131 | IF(flag_min_rreduce) THEN |
---|
132 | IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys |
---|
133 | ENDIF |
---|
134 | dummyso4toso2 = (mSO2mol/mH2SO4mol)*tr_seri(ilon,ilev,id_H2SO4_strat)*(1.0-exp(-pdtphys/rreduce)) |
---|
135 | budg_3D_so2_to_h2so4(ilon,ilev) = budg_3D_so2_to_h2so4(ilon,ilev) + dummyso4toso2 |
---|
136 | tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + dummyso4toso2 |
---|
137 | tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) - (mH2SO4mol/mSO2mol)*dummyso4toso2 |
---|
138 | ENDIF |
---|
139 | ! IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN |
---|
140 | ENDIF |
---|
141 | ! end of IF (flag_H2SO4_photolysis) THEN |
---|
142 | |
---|
143 | !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic |
---|
144 | budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol* & |
---|
145 | (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys |
---|
146 | budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev) |
---|
147 | ENDIF |
---|
148 | ! IF (is_strato(ilon,ilev)) THEN |
---|
149 | ENDDO ! DO (klev) |
---|
150 | ENDDO ! DO (klon) |
---|
151 | |
---|
152 | END SUBROUTINE SO2_TO_H2SO4 |
---|