Context Navigation

new_cloud_venus.F @ 3094

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

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

Download in other formats:

Original Format