source: LMDZ5/trunk/libf/phylmd/StratAer/micphy_tstep.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by oboucher, 8 years ago

Adding a module for stratospheric aerosols with a bin scheme.
The module gets activated with -strataer true compiling option.
May not quite work yet, more testing needed, but should not affect
the rest of LMDz as everything is under a CPP_StratAer key.
File size: 6.6 KB
Line 
1SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
2
3  USE dimphy, ONLY : klon,klev
4  USE aerophys
5  USE infotrac
6  USE phys_local_var_mod, ONLY: mdw, sulf_nucl, sulf_cond_evap, R2SO4, DENSO4, f_r_wet
7  USE nucleation_tstep_mod
8  USE cond_evap_tstep_mod
9  USE sulfate_aer_mod, ONLY : STRAACT
10  USE YOMCST
11
12  IMPLICIT NONE
13
14  !--------------------------------------------------------
15
16  ! transfer variables when calling this routine
17  REAL,INTENT(IN)                               :: pdtphys ! Pas d'integration pour la physique (seconde)
18  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
19  REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
20  REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression pour le mileu de chaque couche (en Pa)
21  REAL,DIMENSION(klon,klev+1),INTENT(IN)        :: paprs   ! pression pour chaque inter-couche (en Pa)
22  REAL,DIMENSION(klon,klev),INTENT(IN)          :: rh      ! humidite relative
23  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
24
25  ! local variables in coagulation routine
26  INTEGER, PARAMETER        :: nbtstep=4  ! Max number of time steps in microphysics per time step in physics
27  INTEGER                   :: it,ilon,ilev,IK,count_tstep
28  REAL                      :: rhoa !H2SO4 number density [molecules/cm3]
29  REAL                      :: ntot !total number of molecules in the critical cluster (ntot>4)
30  REAL                      :: x    ! molefraction of H2SO4 in the critical cluster     
31  REAL Vbin(nbtr_bin)
32  REAL a_xm, b_xm, c_xm
33  REAL PDT, dt
34  REAL H2SO4_init
35  REAL ACTSO4(klon,klev)
36  REAL RRSI(nbtr_bin)
37  REAL nucl_rate
38  REAL cond_evap_rate
39  REAL evap_rate
40  REAL FL(nbtr_bin)
41  REAL ASO4(nbtr_bin)
42  REAL DNDR(nbtr_bin)
43  REAL H2SO4_sat(nbtr_bin)
44
45  DO IK=1,nbtr_bin
46    Vbin(IK)=4.0*RPI*((mdw(IK)/2.)**3)/3.0
47  ENDDO
48
49  !coefficients for H2SO4 density parametrization used for nucleation if ntot<4
50  a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + &
51       & 1.*(88.75606 + 1.*(-75.73729 + 1.*23.43228)))))
52  b_xm = 1.808225e-3 + 1.*(-9.294656e-3 + 1.*(-0.03742148 + 1.*(0.2565321 + &
53       & 1.*(-0.5362872 + 1.*(0.4857736 - 1.*0.1629592)))))
54  c_xm = -3.478524e-6 + 1.*(1.335867e-5 + 1.*(5.195706e-5 + 1.*(-3.717636e-4 + &
55       & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 )))))
56
57  ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap
58  CALL STRAACT(ACTSO4)
59
60  ! compute particle radius in cm RRSI from diameter in m
61  DO it=1,nbtr_bin
62    RRSI(it)=mdw(it)/2.*100.
63  ENDDO
64
65  DO ilon=1, klon
66  DO ilev=1, klev
67  ! only in the stratosphere
68  IF (is_strato(ilon,ilev)) THEN
69    ! initialize sulfur fluxes
70    sulf_nucl(ilon,ilev)=0.0
71    sulf_cond_evap(ilon,ilev)=0.0
72    H2SO4_init=tr_seri(ilon,ilev,id_H2SO4_strat)
73    ! adaptive timestep for nucleation and condensation
74    PDT=pdtphys
75    count_tstep=0
76    DO while(PDT>0.0)
77      count_tstep=count_tstep+1
78      if( count_tstep .GT. nbtstep )  exit
79      ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3)
80      rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) &
81          & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
82      ! compute nucleation rate in kg(H2SO4)/kgA/s
83      CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), &
84             & a_xm,b_xm,c_xm,nucl_rate,ntot,x)
85      ! compute cond/evap rate in kg(H2SO4)/kgA/s
86      CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
87             & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
88             & RRSI,Vbin,FL,ASO4,DNDR)
89      ! consider only condensation (positive FL)
90      DO IK=1,nbtr_bin
91        FL(IK)=MAX(FL(IK),0.)
92      ENDDO
93      ! compute total H2SO4 cond flux for all particles
94      cond_evap_rate=0.0
95      DO IK=1, nbtr_bin
96        cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol
97      ENDDO
98      ! determine appropriate time step
99      dt=(H2SO4_init-H2SO4_sat(nbtr_bin))/float(nbtstep)/MAX(1.e-30, nucl_rate+cond_evap_rate) !cond_evap_rate pos. for cond. and neg. for evap.
100      IF (dt.LT.0.0) THEN
101        dt=PDT
102      ENDIF
103      dt=MIN(dt,PDT)
104      ! update H2SO4 concentration
105      tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt)
106      ! apply cond to bins
107      CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
108      ! apply nucl. to bins
109      CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:))
110      ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
111      sulf_cond_evap(ilon,ilev)=sulf_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
112               & *cond_evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys
113      sulf_nucl(ilon,ilev)=sulf_nucl(ilon,ilev)+mSatom/mH2SO4mol &
114               & *nucl_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys
115      ! update time step
116      PDT=PDT-dt
117    ENDDO
118    ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3)
119    rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) &
120        & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol
121    ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys)
122    CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), &
123           & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), &
124           & RRSI,Vbin,FL,ASO4,DNDR)
125    ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet
126    DO IK=1,nbtr_bin
127      FL(IK)=MAX(FL(IK)*pdtphys,0.-ASO4(IK))/pdtphys
128      ! consider only evap (negative FL)
129      FL(IK)=MIN(FL(IK),0.)
130    ENDDO
131    ! compute total H2SO4 evap flux for all particles
132    evap_rate=0.0
133    DO IK=1, nbtr_bin
134      evap_rate=evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol
135    ENDDO
136    ! update H2SO4 concentration after evap
137    tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys)
138    ! apply evap to bins
139    CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))
140    ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond)
141    sulf_cond_evap(ilon,ilev)=sulf_cond_evap(ilon,ilev)+mSatom/mH2SO4mol &
142             & *evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
143  ENDIF
144  ENDDO
145  ENDDO
146
147  IF (MINVAL(tr_seri(:,:,:)).LT.0.0) THEN
148    DO ilon=1, klon
149    DO ilev=1, klev   
150    DO IK=1, nbtr
151      IF (tr_seri(ilon,ilev,IK).LT.0.0) THEN
152        PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,IK), ilon, ilev, IK
153      ENDIF
154    ENDDO
155    ENDDO
156    ENDDO
157  ENDIF
158
159END SUBROUTINE micphy_tstep
Note: See TracBrowser for help on using the repository browser.