[3526] | 1 | ! |
---|
| 2 | ! $Id: so2_to_h2so4.F90 4601 2023-06-30 22:07:30Z jyg $ |
---|
| 3 | ! |
---|
[2752] | 4 | SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
[2690] | 5 | |
---|
| 6 | USE dimphy, ONLY : klon,klev |
---|
| 7 | USE aerophys |
---|
[3677] | 8 | USE infotrac_phy |
---|
[4601] | 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 | |
---|
[2690] | 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) |
---|
[2695] | 24 | LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato ! stratospheric flag |
---|
[2690] | 25 | |
---|
| 26 | ! local variables in coagulation routine |
---|
| 27 | INTEGER :: i,j,k,nb,ilon,ilev |
---|
[4601] | 28 | REAL :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh |
---|
| 29 | REAL :: rmairdens,rrak0,rrak1,rrate,rreduce |
---|
| 30 | |
---|
[2690] | 31 | !--convert SO2 to H2SO4 |
---|
[2752] | 32 | budg_3D_so2_to_h2so4(:,:)=0.0 |
---|
| 33 | budg_so2_to_h2so4(:)=0.0 |
---|
| 34 | |
---|
[2690] | 35 | DO ilon=1, klon |
---|
[4601] | 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 | |
---|
[2690] | 152 | END SUBROUTINE SO2_TO_H2SO4 |
---|