Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

new_cloud_venus.F

Last change on this file was 3323, checked in by flefevre, 8 months ago
1) revised cloud microphysical parameters (this changes the results) 2) the number of modes is passed as an argument to cloud routines (this does not change the results) 3) some cosmetics to chemparam_mod.F90 (this does not change the results)
Property svn:executable set to ``*
File size: 13.7 KB

Line
1	!******************************************************************************
2	!* venus_SAS_composition SUBROUTINE
3	!* modified from
4	!* PROGRAM PSC_MODEL_E
5	!* by A. Määttänen
6	!* subroutine for LMDZ+photochemistry VENUS
7	!* by A. Stolzenbach
8	!*
9	!* Input/Output files:
10	!* -------------------
11	!*
12	!----------------------------------------------------------------------------
13	SUBROUTINE new_cloud_venus(nbr_mode,
14	+ nblev, nblon,
15	+ TT,PP,
16	+ mrt_wv,mrt_sa,
17	+ mr_wv,mr_sa)
18
19	USE chemparam_mod
20	IMPLICIT NONE
21
22	#include "YOMCST.h"
23
24	integer :: nbr_mode
25	INTEGER, INTENT(IN) :: nblon ! nombre de points horizontaux
26	INTEGER, INTENT(IN) :: nblev ! nombre de couches verticales
27
28	!----------------------------------------------------------------------------
29	! Ambient air state variables:
30	REAL, INTENT(IN), DIMENSION(nblon,nblev) :: mrt_wv,mrt_sa,
31	+ TT,PP
32	REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: mr_wv,mr_sa
33	!----------------------------------------------------------------------------
34	INTEGER :: ilon, ilev, imode
35	!----------------------------------------------------------------------------
36	! Thermodynamic functions:
37	REAL :: ROSAS
38	!----------------------------------------------------------------------------
39	! Auxilary variables:
40	REAL :: NH2SO4,NH2O
41	REAL :: H2SO4_liq,H2O_liq
42	REAL :: CONCM
43	REAL :: MCONDTOT
44	REAL :: RMODE
45	REAL :: WSAFLAG
46	REAL :: K_SAV
47	!----------------------------------------------------------------------------
48	! Ridder's Method variables:
49	REAL :: WVMIN, WVMAX, WVACC
50
51	INTEGER :: NBROOT
52
53	INTEGER :: MAXITE
54	PARAMETER(MAXITE=20)
55
56	INTEGER :: NBRAC
57	PARAMETER(NBRAC=5)
58
59	INTEGER :: FLAG
60
61	REAL :: qmass
62	c masse volumique du coeur (kg.m-3)
63	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_sedim !
64	REAL, PARAMETER :: rho_core = 2500.0
65	!----------------------------------------------------------------------------
66	! External functions needed:
67	REAL :: IRFRMWV
68	!----------------------------------------------------------------------------
69
70
71	! >>> Program starts here:
72
73	!AM Venus
74	! These aerosols will then be given an equilibrium composition for the given size distribution
75
76	! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
77	! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
78	! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
79	!and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
80
81	!===========================================
82	! Debut boucle sur niveau et lat,lon
83	!===========================================
84	! Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0
85	NBRTOT(:,:,:)=0.0E+0
86	WH2SO4(:,:)=0.0E+0
87	rho_droplet(:,:)=0.0E+0
88
89	DO ilev=cloudmin, cloudmax
90	DO ilon=1, nblon
91
92	! Boucle sur les modes
93	RMODE=0.0E+0
94	K_SAV = 0.0
95
96	DO imode=1, nbr_mode
97	IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
98	! RMODE est le rayon modal de la distribution en volume du mode le plus
99	! representatif pour la Mtot
100	RMODE=R_MEDIAN(ilon,ilev,imode)*
101	& EXP(2.(DLOG(STDDEV(ilon,ilev,imode))*2.))
102	K_SAV=K_MASS(ilon,ilev,imode)
103	ENDIF
104	ENDDO ! FIN boucle imode
105
106	! Initialisation des bornes pour WV
107	WVMIN=1.E-90
108	WVMAX=mrt_wv(ilon,ilev)
109
110	! Accuracy de WVeq
111	WVACC=WVMAX*1.0E-3
112
113	! BRACWV borne la fonction f(WV) - WV = 0
114	! de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
115	! avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
116	! Elle fait appel à la fct/ssrtine ITERWV()
117
118	CALL BRACWV(TT(ilon,ilev),PP(ilon,ilev),WVMIN,WVMAX,NBRAC,
119	& RMODE,mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),FLAG,WSAFLAG,NBROOT)
120
121	SELECT CASE(FLAG)
122
123	CASE(1)
124	! Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
125	! IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
126	! Elle fait appel la fct/ssrtine ITERWV()
127
128	WH2SO4(ilon,ilev)=IRFRMWV(TT(ilon,ilev),PP(ilon,ilev),
129	& WVMIN,WVMAX,WVACC,MAXITE,RMODE,
130	& mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
131
132	rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
133
134	! IF (rho_droplet(ilon,ilev).LT.1100.) THEN
135	! PRINT*,'PROBLEM RHO_DROPLET'
136	! PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
137	! PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
138	! PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
139	! PRINT*,'FLAG',FLAG,'NROOT',NBROOT
140	! STOP
141	! ENDIF
142
143	CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
144
145	NH2SO4=mrt_sa(ilon,ilev)*CONCM
146	NH2O=mrt_wv(ilon,ilev)*CONCM
147
148	CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
149	& rho_droplet(ilon,ilev),TT(ilon,ilev),
150	& H2SO4_liq,H2O_liq,MCONDTOT)
151
152	! Boucle sur les modes
153	! *********************************************
154	! AVEC MODE 3 type J. Cimino 1982, Icarus
155	! *********************************************
156	! Mode 1 et 2 avec les Kmass modifies
157	DO imode=1, nbr_mode-1
158	! calcul qmass
159	qmass = (rho_coreqrad*3)/
160	& (rho_coreqrad3 + rho_droplet(ilon,ilev)(1.-qrad**3))
161
162	IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
163	NBRTOT(ilon,ilev,imode)= 1.E-63./(4.RPI)*
164	& K_MASS(ilon,ilev,imode)/(1.-qmassK_MASS(ilon,ilev,3))
165	& MCONDTOTEXP(-4.5DLOG(STDDEV(ilon,ilev,imode))**2.)/
166	& (R_MEDIAN(ilon,ilev,imode)**3.)
167	ELSE
168	NBRTOT(ilon,ilev,imode)=0.0E+0
169	ENDIF
170	ENDDO
171	! Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
172	IF (K_MASS(ilon,ilev,3).GT.0.) THEN
173	NBRTOT(ilon,ilev,3)= 1.E-63./(4.RPI)*
174	& K_MASS(ilon,ilev,3)MCONDTOT
175	& EXP(-4.5DLOG(STDDEV(ilon,ilev,3))*2.)/
176	& (R_MEDIAN(ilon,ilev,3)**3.)
177	ELSE
178	NBRTOT(ilon,ilev,3)=0.0E+0
179	ENDIF
180
181	! Passage de #/m3 en VMR
182	H2O_liq=H2O_liq/CONCM
183	H2SO4_liq=H2SO4_liq/CONCM
184
185	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
186	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
187
188	! Problemes quand on a condense tout, on peut obtenir des -1e-24
189	! aprs la soustraction et conversion de ND VMR
190	IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
191	IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
192
193
194
195	CASE(2)
196	! Cas NROOT=0 mais proche de 0
197
198	WH2SO4(ilon,ilev)=WSAFLAG
199
200	rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
201
202	! ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
203	! ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
204	! incoherentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
205	! Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
206	! donne par ROSAS (cf test externe en Python), sinon, la valeur est trop
207	! basse (de l'ordre de 1000 kg/m3) et correspond parfois a la valeur avec
208	! WSA=0.1 (pas totalement sur)
209	! En tous cas, incoherent avec ce qui est attendue pour le WSA et T donnee
210	! La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
211	! la bonne valeur (tests externes Python confirment)
212
213	IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
214	& (WH2SO4(ilon,ilev).GT.0.1))THEN
215	PRINT*,'PROBLEM RHO_DROPLET'
216	PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
217	PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
218	PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
219	PRINT*,'FLAG',FLAG,'NROOT',NBROOT
220	STOP
221	ENDIF
222
223
224	CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
225
226	NH2SO4=mrt_sa(ilon,ilev)*CONCM
227	NH2O=mrt_wv(ilon,ilev)*CONCM
228
229	CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
230	& rho_droplet(ilon,ilev),TT(ilon,ilev),
231	& H2SO4_liq,H2O_liq,MCONDTOT)
232
233	! Boucle sur les modes
234	! *********************************************
235	! AVEC MODE 3 type J. Cimino 1982, Icarus
236	! *********************************************
237	! Mode 1 et 2 avec alcul coeff*Kmass
238	DO imode=1, nbr_mode-1
239	! calcul qmass
240	qmass = (rho_coreqrad*3)/
241	& (rho_coreqrad3 + rho_droplet(ilon,ilev)(1.-qrad**3))
242
243	IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
244	NBRTOT(ilon,ilev,imode)= 1.E-63./(4.RPI)*
245	& K_MASS(ilon,ilev,imode)/(1.-qmassK_MASS(ilon,ilev,3))
246	& MCONDTOTEXP(-4.5DLOG(STDDEV(ilon,ilev,imode))**2.)/
247	& (R_MEDIAN(ilon,ilev,imode)**3.)
248	ELSE
249	NBRTOT(ilon,ilev,imode)=0.0E+0
250	ENDIF
251	ENDDO
252	! Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
253	IF (K_MASS(ilon,ilev,3).GT.0.) THEN
254	NBRTOT(ilon,ilev,3)= 1.E-63./(4.RPI)*
255	& K_MASS(ilon,ilev,3)MCONDTOT
256	& EXP(-4.5DLOG(STDDEV(ilon,ilev,3))*2.)/
257	& (R_MEDIAN(ilon,ilev,3)**3.)
258	ELSE
259	NBRTOT(ilon,ilev,3)=0.0E+0
260	ENDIF
261
262
263	! Passage de #/m3 en VMR
264	H2O_liq=H2O_liq/CONCM
265	H2SO4_liq=H2SO4_liq/CONCM
266
267	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
268	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
269
270	! Problmes quand on a condense tout, on peut obtenir des -1e-24
271	! aprs la soustraction et conversion de ND VMR
272	IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
273	IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
274
275	CASE(3)
276	! Cas 0 NROOT
277	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
278	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
279	rho_droplet(ilon,ilev)=0.0E+0
280	WH2SO4(ilon,ilev)=0.0E+0
281	DO imode=1, nbr_mode
282	NBRTOT(ilon,ilev,imode)=0.0E+0
283	ENDDO
284
285	END SELECT
286	ENDDO !FIN boucle ilon
287	ENDDO !FIN boucle ilev
288
289	END SUBROUTINE new_cloud_venus
290
291
292	!******************************************************************
293	SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
294	+ T,H2SO4COND,H2OCOND,RMTOT)
295
296	! DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
297	! FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
298	! SIZE DISTRIBTUION
299	! ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
300	! ---------------------------------------------------------------
301	! INPUT:
302	! H2SO4: #/m3 of total H2SO4
303	! H2O : #/m3 of total H2O
304	! WSA: aerosol H2SO4 weight fraction (fraction)
305	! DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume)
306	! for total mass, almost same result with ro=1.67 gr/cm3
307	! RSTDEV: standard deviation of aerosol distribution (no unit)
308	! RADIUS: MEDIAN radius (m)
309	! T: temperature (K)
310	!
311	! OUTPUT:
312	! RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
313	! mais rho_droplet et M_tot_distrib doivent etre de meme dimension
314	! H2OCOND
315	! H2SO4COND
316
317
318
319	IMPLICIT NONE
320
321	REAL, INTENT(IN) :: H2SO4, H2O, WSA
322	REAL, INTENT(IN) :: DENSO4, T
323	REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
324	! working variables
325	REAL :: RMH2S4
326	REAL :: DND2,pstand,lpar,acidps
327	REAL :: x1, satpacid
328	REAL , DIMENSION(2):: act
329	!
330	! masse of an H2SO4 molecule (kg)
331	RMH2S4=98.078/(6.02214129E+26)
332
333	pstand=1.01325E+5 !Pa 1 atm pressure
334
335	x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
336
337	call zeleznik(x1,t,act)
338
339	!pure acid satur vapor pressure
340	lpar= -11.695+DLOG(pstand) ! Zeleznik
341	acidps=1/360.15-1.0/t+0.38/545.
342	+ *(1.0+DLOG(360.15/t)-360.15/t)
343	acidps = 10156.0*acidps +lpar
344	acidps = DEXP(acidps) !Pa
345
346	!acid sat.vap.PP over mixture (flat surface):
347	satpacid=act(2)*acidps ! Pa
348
349	! Conversion from Pa to N.D #/m3
350	DND2=satpacid/(1.3806488E-23*T)
351	! Conversion from N.D #/m3 TO #/cm3
352	! DND2=DND2*1.d-6
353
354	! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
355	IF (H2SO4.GE.DND2) THEN
356	H2SO4COND=H2SO4-DND2
357	! calcul de H2O cond correspondant a H2SO4 cond
358	H2OCOND=H2SO4COND98.078(1.0-WSA)/(18.0153*WSA)
359
360	! RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
361	! RMTOT = Mtot/rho_droplet
362	! RMTOT=M_distrib/rho_droplet
363
364	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
365
366	! Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
367	ELSE
368	H2SO4COND=0.0E+0
369	H2OCOND=0.0E+0
370	RMTOT=0.0E+0
371	END IF
372
373	! Test si H2O en defaut H2Ocond>H2O dispo
374	IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
375
376	! Si H2O en dfaut, on as pas le bon WSA!
377	! En effet, normalement, on a exactement le WSA correspondant a
378	! WVg + WVl = WVtot
379	! Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou
380	! quand T est grand (Lower Haze), le modle reprsente mal le WSA
381	! cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA
382
383	! PRINT*,'PROBLEM H2O EN DEFAUT'
384	! PRINT*,'H2OCOND',H2OCOND,'H2O',H2O
385	! PRINT*,'WSA',WSA,'RHO',DENSO4
386	! STOP
387
388
389	! On peut alors condenser tout le H2O dispo
390	H2OCOND=H2O
391	! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
392	H2SO4COND=H2OCOND18.0153WSA/(98.078*(1.0-WSA))
393
394	! RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
395	! RMTOT=Volume of aerosol m3 /m3 of air
396	! Volume of aerosol/m3 air
397
398	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
399
400	END IF
401
402	END SUBROUTINE CALCM_SAT
403

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

Download in other formats:

Original Format