Context Navigation

micphy_tstep.F90 @ 5425

Last change on this file since 5425 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 jyg $
3	!
4	SUBROUTINE 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.0RPI((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
184	END SUBROUTINE micphy_tstep

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format