source: LMDZ6/branches/blowing_snow/libf/phylmd/StratAer/micphy_tstep.F90

Last change on this file was 4293, checked in by dcugnet, 2 years ago

Commit for Nicolas: fixes for StratAer?.
  • Property svn:keywords set to Id
File size: 7.7 KB
Line 
1!
2! $Id: micphy_tstep.F90 4293 2022-10-04 19:42:41Z aborella $
3!
4SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
5
6  USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain
7  USE dimphy, ONLY : klon,klev
8  USE aerophys
9  USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat
10  USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet
11  USE nucleation_tstep_mod
12  USE cond_evap_tstep_mod
13  USE sulfate_aer_mod, ONLY : STRAACT
14  USE YOMCST, ONLY : RPI, RD, RG
15  USE print_control_mod, ONLY: lunout
16  USE strataer_mod
17 
18  IMPLICIT NONE
19
20  !--------------------------------------------------------
21
22  ! transfer variables when calling this routine
23  REAL,INTENT(IN)                               :: pdtphys ! Pas d'integration pour la physique (seconde)
24  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
25  REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
26  REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression pour le mileu de chaque couche (en Pa)
27  REAL,DIMENSION(klon,klev+1),INTENT(IN)        :: paprs   ! pression pour chaque inter-couche (en Pa)
28  REAL,DIMENSION(klon,klev),INTENT(IN)          :: rh      ! humidite relative
29  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
30
31  ! local variables in coagulation routine
32  INTEGER, PARAMETER        :: nbtstep=4  ! Max number of time steps in microphysics per time step in physics
33  INTEGER                   :: it,ilon,ilev,count_tstep
34  REAL                      :: rhoa !H2SO4 number density [molecules/cm3]
35  REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
36  REAL                      :: x    ! molefraction of H2SO4 in the critical cluster     
37  REAL Vbin(nbtr_bin)
38  REAL a_xm, b_xm, c_xm
39  REAL PDT, dt
40  REAL H2SO4_init
41  REAL ACTSO4(klon,klev)
42  REAL RRSI(nbtr_bin)
43  REAL nucl_rate
44  REAL cond_evap_rate
45  REAL evap_rate
46  REAL FL(nbtr_bin)
47  REAL ASO4(nbtr_bin)
48  REAL DNDR(nbtr_bin)
49  REAL H2SO4_sat
50
51  DO it=1,nbtr_bin
52    Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0
53  ENDDO
54
55  !coefficients for H2SO4 density parametrization used for nucleation if ntot<4
56  a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + &
57       & 1.*(88.75606 + 1.*(-75.73729 + 1.*23.43228)))))
58  b_xm = 1.808225e-3 + 1.*(-9.294656e-3 + 1.*(-0.03742148 + 1.*(0.2565321 + &
59       & 1.*(-0.5362872 + 1.*(0.4857736 - 1.*0.1629592)))))
60  c_xm = -3.478524e-6 + 1.*(1.335867e-5 + 1.*(5.195706e-5 + 1.*(-3.717636e-4 + &
61       & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
62
63  ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
64  CALL STRAACT(ACTSO4)
65
66  ! compute particle radius in cm RRSI from diameter in m
67  DO it=1,nbtr_bin
68    RRSI(it)=mdw(it)/2.*100.
69  ENDDO
70
71  DO ilon=1, klon
72!
73!--initialisation of diagnostic
74  budg_h2so4_to_part(ilon)=0.0
75!
76  DO ilev=1, klev
77!
78!--initialisation of diagnostic
79  budg_3D_nucl(ilon,ilev)=0.0
80  budg_3D_cond_evap(ilon,ilev)=0.0
81!
82  ! only in the stratosphere
83  IF (is_strato(ilon,ilev)) THEN
84    ! initialize sulfur fluxes
85    H2SO4_init=tr_seri(ilon,ilev,id_H2SO4_strat)
86    ! adaptive timestep for nucleation and condensation
87    PDT=pdtphys
88    count_tstep=0
89    DO WHILE (PDT>0.0)
90      count_tstep=count_tstep+1
91      IF (count_tstep .GT. nbtstep)  EXIT
92      ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3)
93      rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) &
94          & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
95      ! compute nucleation rate in kg(H2SO4)/kgA/s
96      CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), &
97           & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
98      !NL - add nucleation box (if flag on)
99      IF (flag_nuc_rate_box) THEN
100         IF (latitude_deg(ilon).LE.nuclat_min .OR. latitude_deg(ilon).GE.nuclat_max &
101              .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev).LE.nucpres_min) THEN
102            nucl_rate=0.0
103         ENDIF
104      ENDIF
105      ! compute cond/evap rate in kg(H2SO4)/kgA/s
106      CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
107             & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
108             & RRSI,Vbin,FL,ASO4,DNDR)
109      ! Compute H2SO4 saturate vapor for big particules
110      H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol)
111      ! consider only condensation (positive FL)
112      DO it=1,nbtr_bin
113        FL(it)=MAX(FL(it),0.)
114      ENDDO
115      ! compute total H2SO4 cond flux for all particles
116      cond_evap_rate=0.0
117      DO it=1, nbtr_bin
118        cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol
119      ENDDO
120      ! determine appropriate time step
121      dt=(H2SO4_init-H2SO4_sat)/float(nbtstep)/MAX(1.e-30, nucl_rate+cond_evap_rate) !cond_evap_rate pos. for cond. and neg. for evap.
122      IF (dt.LT.0.0) THEN
123        dt=PDT
124      ENDIF
125      dt=MIN(dt,PDT)
126      ! update H2SO4 concentration
127      tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt)
128      ! apply cond to bins
129      CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
130      ! apply nucl. to bins
131      CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:))
132      ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
133      budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
134               & *cond_evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys
135      budg_3D_nucl(ilon,ilev)=budg_3D_nucl(ilon,ilev)+mSatom/mH2SO4mol &
136               & *nucl_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys
137      ! update time step
138      PDT=PDT-dt
139    ENDDO
140    ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3)
141    rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) &
142        & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
143    ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys)
144    CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
145           & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
146           & RRSI,Vbin,FL,ASO4,DNDR)
147    ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
148    DO it=1,nbtr_bin
149      FL(it)=MAX(FL(it)*pdtphys,0.-ASO4(it))/pdtphys
150      ! consider only evap (negative FL)
151      FL(it)=MIN(FL(it),0.)
152    ENDDO
153    ! compute total H2SO4 evap flux for all particles
154    evap_rate=0.0
155    DO it=1, nbtr_bin
156      evap_rate=evap_rate+tr_seri(ilon,ilev,it+nbtr_sulgas)*FL(it)*mH2SO4mol
157    ENDDO
158    ! update H2SO4 concentration after evap
159    tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys)
160    ! apply evap to bins
161    CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
162    ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
163    budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
164             & *evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
165    ! compute vertically integrated flux due to the net effect of nucleation and condensation/evaporation
166    budg_h2so4_to_part(ilon)=budg_h2so4_to_part(ilon)+(H2SO4_init-tr_seri(ilon,ilev,id_H2SO4_strat)) &
167             & *mSatom/mH2SO4mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
168  ENDIF
169  ENDDO
170  ENDDO
171
172  IF (MINVAL(tr_seri).LT.0.0) THEN
173    DO ilon=1, klon
174    DO ilev=1, klev   
175    DO it=1, nbtr
176      IF (tr_seri(ilon,ilev,it).LT.0.0) THEN
177         WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it
178      ENDIF
179    ENDDO
180    ENDDO
181    ENDDO
182  ENDIF
183
184END SUBROUTINE micphy_tstep
Note: See TracBrowser for help on using the repository browser.