source: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/so2_to_h2so4.F90 @ 5403

Last change on this file since 5403 was 5101, checked in by abarral, 5 months ago

Handle DEBUG_IO in lmdz_cppkeys_wrapper.F90
Transform some files .F -> .[fF]90
[ne compile pas à cause de writefield_u non défini - en attente de réponse Laurent]
  • Property svn:keywords set to Id
File size: 8.1 KB
Line 
1
2! $Id: so2_to_h2so4.F90 5101 2024-07-23 06:22:55Z abarral $
3
4SUBROUTINE 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 lmdz_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)>0.0 .AND. SO2_lifetime(ilon,ilev)<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 >= 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 < (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)>0.0 .AND. H2SO4_lifetime(ilon,ilev)<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 < (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 
152END SUBROUTINE SO2_TO_H2SO4
Note: See TracBrowser for help on using the repository browser.