source: LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 @ 5150

!
! $Id: so2_to_h2so4.F90 4601 2023-06-30 22:07:30Z 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, 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
  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

  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
           
           
           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)
        ENDIF
        ! IF (is_strato(ilon,ilev)) THEN
     ENDDO ! DO (klev)
  ENDDO ! DO (klon)
  
END SUBROUTINE SO2_TO_H2SO4