Context Navigation

new_cloud_venus.F @ 1543

Last change on this file since 1543 was 1442, checked in by slebonnois, 10 years ago
SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools
Property svn:executable set to ``*
File size: 57.9 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 :: RHODROPLET
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=20)
57
58	INTEGER :: FLAG
59	!----------------------------------------------------------------------------
60
61	!----------------------------------------------------------------------------
62	! External functions needed:
63	REAL :: IRFRMWV
64	!----------------------------------------------------------------------------
65
66
67	! >>> Program starts here:
68
69	!AM Venus
70	! These aerosols will then be given an equilibrium composition for the given size distribution
71
72	! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
73	! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
74	! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
75	!and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
76
77	!===========================================
78	! Debut boucle sur niveau et lat,lon
79	!===========================================
80	! Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0
81	NBRTOT(:,:,:)=0.0E+0
82	WH2SO4(:,:)=0.0E+0
83	rho_droplet(:,:)=0.0E+0
84
85	DO ilev=cloudmin, cloudmax
86	DO ilon=1, nblon
87
88	! Boucle sur les modes
89	RMODE=0.0E+0
90	K_SAV = 0.0
91
92	DO imode=1, nbr_mode
93	IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
94	! RMODE est le rayon modal de la distribution en volume du mode le plus
95	! representatif pour la Mtot
96	RMODE=R_MEDIAN(ilon,ilev,imode)*
97	& EXP(2.(DLOG(STDDEV(ilon,ilev,imode))*2.))
98	K_SAV=K_MASS(ilon,ilev,imode)
99	ENDIF
100	ENDDO ! FIN boucle imode
101
102	! Initialisation des bornes pour WV
103	WVMIN=1.E-90
104	WVMAX=mrt_wv(ilon,ilev)
105
106	! Accuracy de WVeq
107	WVACC=WVMAX*1.0E-3
108
109	! BRACWV borne la fonction f(WV) - WV = 0
110	! de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
111	! avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
112	! Elle fait appel à la fct/ssrtine ITERWV()
113
114	CALL BRACWV(WVMIN,WVMAX,NBRAC,RMODE,
115	& mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),TT(ilon,ilev),
116	& PP(ilon,ilev),FLAG,WSAFLAG,NBROOT)
117
118	SELECT CASE(FLAG)
119
120	CASE(1)
121	! Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
122	! IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
123	! Elle fait appel la fct/ssrtine ITERWV()
124
125	WH2SO4(ilon,ilev)=IRFRMWV(WVMIN,WVMAX,WVACC,MAXITE,RMODE,
126	& TT(ilon,ilev),PP(ilon,ilev),
127	& mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
128
129	rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
130	& TT(ilon,ilev))
131
132	! IF (rho_droplet(ilon,ilev).LT.1100.) THEN
133	! PRINT*,'PROBLEM RHO_DROPLET'
134	! PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
135	! PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
136	! PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev),
137	! & TT(ilon,ilev))
138	! PRINT*,'FLAG',FLAG,'NROOT',NBROOT
139	! STOP
140	! ENDIF
141
142	CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
143
144	NH2SO4=mrt_sa(ilon,ilev)*CONCM
145	NH2O=mrt_wv(ilon,ilev)*CONCM
146
147	CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
148	& rho_droplet(ilon,ilev),TT(ilon,ilev),
149	& H2SO4_liq,H2O_liq,MCONDTOT)
150
151	! Boucle sur les modes
152	DO imode=1, nbr_mode
153	IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
154	NBRTOT(ilon,ilev,imode)= 1.E-63./(4.RPI)*
155	& K_MASS(ilon,ilev,imode)MCONDTOT
156	& EXP(-4.5DLOG(STDDEV(ilon,ilev,imode))*2.)/
157	& (R_MEDIAN(ilon,ilev,imode)**3.)
158	ELSE
159	NBRTOT(ilon,ilev,imode)=0.0E+0
160	ENDIF
161	ENDDO
162
163	! Passage de #/m3 en VMR
164	H2O_liq=H2O_liq/CONCM
165	H2SO4_liq=H2SO4_liq/CONCM
166
167	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
168	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
169
170	! Problemes quand on a condense tout, on peut obtenir des -1e-24
171	! aprs la soustraction et conversion de ND VMR
172	IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
173	IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
174
175
176
177	CASE(2)
178	! Cas NROOT=0 mais proche de 0
179
180	WH2SO4(ilon,ilev)=WSAFLAG
181
182	rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
183	& TT(ilon,ilev))
184
185	! ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
186	! ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
187	! incohrentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
188	! Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
189	! donne par RHODROPLET (cf test externe en Python), sinon, la valeur est trop
190	! basse (de l'ordre de 1000 kg/m3) et correspond parfois la valeur avec
191	! WSA=0.1 (pas totalement sur)
192	! En tous cas, incoherent avec ce qui est attendue pour le WSA et T donnee
193	! La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
194	! la bonne valeur (tests externes Python confirment)
195
196	IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
197	& (WH2SO4(ilon,ilev).GT.0.1))THEN
198	PRINT*,'PROBLEM RHO_DROPLET'
199	PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
200	PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
201	PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev),
202	& TT(ilon,ilev))
203	PRINT*,'FLAG',FLAG,'NROOT',NBROOT
204	STOP
205	ENDIF
206
207
208	CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
209
210	NH2SO4=mrt_sa(ilon,ilev)*CONCM
211	NH2O=mrt_wv(ilon,ilev)*CONCM
212
213	CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
214	& rho_droplet(ilon,ilev),TT(ilon,ilev),
215	& H2SO4_liq,H2O_liq,MCONDTOT)
216
217	! Boucle sur les modes
218	DO imode=1, nbr_mode
219	IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
220	NBRTOT(ilon,ilev,imode)= 1.E-63./(4.RPI)*
221	& K_MASS(ilon,ilev,imode)MCONDTOT
222	& EXP(-4.5DLOG(STDDEV(ilon,ilev,imode))*2.)/
223	& (R_MEDIAN(ilon,ilev,imode)**3.)
224	ELSE
225	NBRTOT(ilon,ilev,imode)=0.0E+0
226	ENDIF
227	ENDDO
228
229	! Passage de #/m3 en VMR
230	H2O_liq=H2O_liq/CONCM
231	H2SO4_liq=H2SO4_liq/CONCM
232
233	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
234	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
235
236	! Problmes quand on a condense tout, on peut obtenir des -1e-24
237	! aprs la soustraction et conversion de ND VMR
238	IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
239	IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
240
241	CASE(3)
242	! Cas 0 NROOT
243	mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
244	mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
245	rho_droplet(ilon,ilev)=0.0E+0
246	WH2SO4(ilon,ilev)=0.0E+0
247	DO imode=1, nbr_mode
248	NBRTOT(ilon,ilev,imode)=0.0E+0
249	ENDDO
250
251	END SELECT
252	ENDDO !FIN boucle ilon
253	ENDDO !FIN boucle ilev
254
255	END SUBROUTINE new_cloud_venus
256
257
258	!*****************************************************************************
259	!* SUBROUTINE ITERWV()
260	SUBROUTINE ITERWV(WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,
261	+ TAIR,PAIR,RADIUS)
262	!*****************************************************************************
263	!* Cette routine est la solution par itration afin de trouver WSA pour un WV,
264	!* et donc LPPWV, donn. Ce qui nous donne egalement le WV correspondant au
265	!* WSA solution
266	!* For VenusGCM by A. Stolzenbach 07/2014
267	!* OUTPUT: WVEQ et WSAOUT
268
269	IMPLICIT NONE
270	REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS
271
272	REAl, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ
273
274	REAL :: WSAMIN, WSAMAX, WSAACC
275	PARAMETER(WSAACC=0.001)
276
277	REAL :: LPPWV
278
279	INTEGER :: MAXITSA, NBRACSA, NBROOT
280	PARAMETER(MAXITSA=20)
281	PARAMETER(NBRACSA=20)
282
283	LOGICAl :: FLAG1,FLAG2
284
285	! External Function
286	REAl :: IRFRMSA, WVCOND
287
288	IF (RADIUS.LT.1E-30) THEN
289	PRINT*,'RMODE == 0 FLAG 3'
290	STOP
291	ENDIF
292	! Initialisation WSA=[0.1,1.0]
293	WSAMIN=0.1
294	WSAMAX=1.0
295
296	LPPWV=DLOG(PAIR*WV)
297
298	! Appel Bracket de KEEQ
299	CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR,
300	& LPPWV,FLAG1,FLAG2,NBROOT)
301
302	IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN
303	! Appel Ridder's Method
304
305	WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA,
306	& RADIUS,TAIR,PAIR,LPPWV,NBROOT)
307	! IF (WSAOUT.EQ.1.0) WSAOUT=0.999999
308	! IF (WSAOUT.LT.0.1) WSAOUT=0.1
309
310	! Si BRACWSA ne trouve aucun ensemble solution KEEQ=0 on fixe WSA a 0.9999 ou 0.1
311	ELSE
312	IF (FLAG1.AND.(.NOT.FLAG2)) WSAOUT=0.999999
313	IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT=WSAMIN
314	IF (FLAG1.AND.FLAG2) THEN
315	PRINT*,'FLAGs BARCWSA tous TRUE'
316	STOP
317	ENDIF
318	ENDIF
319
320
321	! WVEQ output correspondant a WVliq lie a WSA calcule
322	WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT)
323	WVEQOUT=(WVLIQ+WV)/WVTOT-1.0
324
325	END SUBROUTINE ITERWV
326
327
328	!*****************************************************************************
329	!* SUBROUTINE BRACWV()
330	SUBROUTINE BRACWV(XA,XB,N,RADIUS,WVTOT,SATOT,TAIR,PAIR,
331	+ FLAGWV,WSAFLAG,NROOT)
332	!*****************************************************************************
333	!* Bracket de ITERWV
334	!* From Numerical Recipes
335	!* Adapted for VenusGCM A. Stolzenbach 07/2014
336	!* X est WVinput
337	!* OUTPUT: XA et XB
338
339	IMPLICIT NONE
340
341	REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR,PAIR
342	INTEGER, INTENT(IN) :: N
343
344	REAL, INTENT(INOUT) :: XA,XB
345	REAL, INTENT(OUT) :: WSAFLAG
346
347	INTEGER :: I,J
348
349	INTEGER, INTENT(OUT) :: NROOT
350
351	REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT
352	REAL :: XMAX,XMIN,WVEQACC
353
354	INTEGER, INTENT(OUT) :: FLAGWV
355
356	! WVEQACC est le seuil auquel on accorde un WSA correct meme
357	! si il ne fait pas partie d'une borne. Utile quand le modele
358	! s'approche de 0 mais ne l'atteint pas.
359	WVEQACC=1.0E-3
360
361	FLAGWV=1
362
363	NROOT=0
364
365	X=XA
366	XMAX=XB
367	XMIN=XA
368
369	! CAS 1 On borne la fonction (WVEQ=0)
370
371	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
372	FP=WVEQ
373
374	DO I=1,N-1
375	X=(1.-DLOG(REAL(N-I))/DLOG(REAL(N)))*XMAX
376	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
377	FC=WVEQ
378
379	IF ((FP*FC).LT.0.0) THEN
380	NROOT=NROOT+1
381	! Si NROOT>1 on place la borne sup output la borne min du calcul en i
382	IF (NROOT.GT.1) THEN
383	XB=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
384	ENDIF
385
386	IF (I.EQ.1) THEN
387	XA=XMIN
388	ELSE
389	XA=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
390	ENDIF
391	XB=X
392	ENDIF
393	FP=FC
394	ENDDO
395
396	! CAS 2 on refait la boucle pour tester si WVEQ est proche de 0
397	! avec le seuil WVEQACC
398	IF (NROOT.EQ.0) THEN
399	X=XMIN
400	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
401	+ TAIR,PAIR,RADIUS)
402	DO J=1,N-1
403	X=(1.-DLOG(REAL(N-J))/DLOG(REAL(N)))*XMAX
404	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
405	+ TAIR,PAIR,RADIUS)
406
407	IF (ABS(WVEQ).LE.WVEQACC) THEN
408	WSAFLAG=WSAOUT
409	FLAGWV=2
410	RETURN
411	ENDIF
412	ENDDO
413
414	! CAS 3 Pas de borne, WVEQ jamais proche de 0
415	FLAGWV=3
416	RETURN
417	ENDIF
418
419	END SUBROUTINE BRACWV
420
421	!*****************************************************************************
422	!* SUBROUTINE BRACWSA()
423	SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL,
424	+ NROOT)
425	!*****************************************************************************
426	!* Bracket de KEEQ
427	!* From Numerical Recipes
428	!* Adapted for VenusGCM A. Stolzenbach 07/2014
429
430	IMPLICIT NONE
431
432	!----------------------------------------------------------------------------
433	! External functions needed:
434	REAl KEEQ
435	!----------------------------------------------------------------------------
436
437	REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP
438	INTEGER, INTENT(IN) :: N
439
440	REAL, INTENT(INOUT) :: XA,XB
441
442	INTEGER, INTENT(OUT) :: NROOT
443
444	INTEGER :: I, J
445
446	REAL :: DX, FP, FC, X
447
448	LOGICAL, INTENT(OUT) :: FLAGH,FLAGL
449
450
451	FLAGL=.FALSE.
452	FLAGH=.FALSE.
453	NROOT=0
454	DX=(XB-XA)/N
455	X=XA
456	FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
457
458	DO I=1,N
459	X=X+DX
460	FC=KEEQ(RADIUS,TAIR,X,LPPWVINP)
461
462	IF ((FP*FC).LE.0.) THEN
463	NROOT=NROOT+1
464	XA=X-DX
465	XB=X
466	! RETURN
467	! IF (NROOT.GT.1) THEN
468	! PRINT*,'On a plus d1 intervalle KEEQ=0'
469	! PRINT*,'Probleme KEEQ=0 => 1 racine en theorie'
470	! X=X-(I*DX)
471	! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
472	! PRINT*,'KEEQ(WSA)',FP,X,TAIR
473	! DO J=1,N
474	! X=X+DX
475	! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
476	! PRINT*,'KEEQ(WSA)',FP,X
477	! ENDDO
478	! STOP
479	! ENDIF
480	ENDIF
481
482	FP=FC
483	ENDDO
484
485	IF (NROOT.EQ.0) THEN
486	! PRINT*,'On a 0 intervalle KEEQ=0'
487	! PRINT*,'Probleme KEEQ=0 => 1 racine en theorie'
488	! PRINT*,'XA',XA,'KEEQ',KEEQ(RADIUS,TAIR,XA,LPPWVINP)
489	! PRINT*,'XB',XB,'KEEQ',KEEQ(RADIUS,TAIR,XB,LPPWVINP)
490	! PRINT*,'TT',TAIR
491	! PRINT*,'RADIUS',RADIUS
492	! PRINT*,'NBRAC',N
493	! STOP
494
495	! X=XA
496	! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
497	! PRINT*,'KEEQ(WSA)',FP,X,TAIR
498	! DO I=1,N
499	! X=X+DX
500	! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
501	! PRINT*,'KEEQ(WSA)',FP,X,TAIR
502	! ENDDO
503
504
505	! Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]
506	IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
507	& ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE.
508	! On fixe flag low TRUE pour WSA = 0.1
509	IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
510	& ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).LT.0.0) FLAGL=.TRUE.
511	! STOP
512	ENDIF
513
514	END SUBROUTINE BRACWSA
515
516
517	!*****************************************************************************
518	!* REAL FUNCTION WVCOND()
519	REAL FUNCTION WVCOND(WSA,T,P,SAt)
520	!*****************************************************************************
521	!* Condensation de H2O selon WSA, T et P et H2SO4tot
522	!*
523	!* Adapted for VenusGCM A. Stolzenbach 07/2014
524	! INPUT:
525	! SAt : VMR of total H2SO4
526	! WSA: aerosol H2SO4 weight fraction (fraction)
527	! T: temperature (K)
528	! P: pressure (Pa)
529	! OUTPUT:
530	! WVCOND : VMR H2O condense
531
532	! USE chemparam_mod
533
534	IMPLICIT NONE
535
536	REAL, INTENT(IN) :: SAt, WSA
537	REAL, INTENT(IN) :: T, P
538
539	! working variables
540	REAL SA, WV
541	REAL DND2,pstand,lpar,acidps
542	REAL x1, satpacid
543	REAL , DIMENSION(2):: act
544	REAL CONCM
545	REAL NH2SO4
546	REAL H2OCOND, H2SO4COND
547
548
549	CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS!
550
551	NH2SO4=SAt*CONCM
552
553	pstand=1.01325E+5 !Pa 1 atm pressure
554
555	x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
556
557	CALL zeleznik(x1,T,act)
558
559	!pure acid satur vapor pressure
560	lpar= -11.695+DLOG(pstand) ! Zeleznik
561	acidps=1/360.15-1.0/T+0.38/545.
562	& *(1.0+DLOG(360.15/T)-360.15/T)
563	acidps = 10156.0*acidps +lpar
564	acidps = DEXP(acidps) !Pa
565
566	!acid sat.vap.PP over mixture (flat surface):
567	satpacid=act(2)*acidps ! Pa
568
569	! Conversion from Pa to N.D #/m3
570	DND2=satpacid/(1.3806488E-23*T)
571
572	! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
573	IF (NH2SO4.GT.DND2) THEN
574	H2SO4COND=NH2SO4-DND2
575	! calcul de H2O cond correspondant a H2SO4 cond
576	H2OCOND=H2SO4COND98.078(1.0-WSA)/(18.0153*WSA)
577
578	! Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30
579	ELSE
580	H2OCOND=1.0E-30*CONCM
581	END IF
582
583	!*****************************************************
584	! ATTENTION: Ici on ne prends pas en compte
585	! si H2O en defaut!
586	! On veut la situation theorique
587	! a l'equilibre
588	!*****************************************************
589	! Test si H2O en defaut H2Ocond>H2O dispo
590	! IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN
591
592	! On peut alors condenser tout le H2O dispo
593	! H2OCOND=NH2O
594	! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
595	! H2SO4COND=H2OCOND18.0153WSA/(98.078*(1.0-WSA))
596
597	! END IF
598
599	! Calcul de H2O condense VMR
600	WVCOND=H2OCOND/CONCM
601
602	END FUNCTION WVCOND
603
604	!*****************************************************************************
605	!* REAL FUNCTION IRFRMWV()
606	REAL FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,
607	+ WVTOT,SATOT,NROOT)
608	!*****************************************************************************
609	!* Iterative Root Finder Ridder's Method for Water Vapor calculus
610	!* From Numerical Recipes
611	!* Adapted for VenusGCM A. Stolzenbach 07/2014
612	!*
613	!* Les iterations sur [X1,X2] sont [WV1,WV2]
614	!* la variable X est WV
615	!* IRFRMWV sort en OUTPUT : WSALOC pour ITERWV=0 (ou WVEQ=0)
616
617	IMPLICIT NONE
618
619	REAL, INTENT(IN) :: X1, X2
620	REAL, INTENT(IN) :: XACC
621	INTEGER, INTENT(IN) :: MAXIT,NROOT
622
623	! LOCAL VARIABLES
624	REAL :: XL, XH, XM, XNEW, X
625	REAL :: WSALOC, WVEQ, WVLIQ
626	REAL :: FL, FH, FM, FNEW
627	REAL :: ANS, S, FSIGN
628	INTEGER i
629
630	! External variables needed:
631	REAL, INTENT(IN) :: TAIR,PAIR
632	REAL, INTENT(IN) :: WVTOT,SATOT
633	REAL, INTENT(IN) :: RADIUS
634
635
636	! Initialisation
637	X=X1
638	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
639	FL=WVEQ
640	X=X2
641	CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
642	FH=WVEQ
643
644	! Test Bracketed values
645	IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.
646	& ((FL.GT.0.).AND.(FH.LT.0.)))
647	& THEN
648	XL=X1
649	XH=X2
650	ANS=-9.99e99
651
652	DO i=1, MAXIT
653	XM=0.5*(XL+XH)
654	CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
655	& TAIR,PAIR,RADIUS)
656	FM=WVEQ
657	S=SQRT(FMFM-FLFH)
658
659	IF (S.EQ.0.0) THEN
660	IRFRMWV=WSALOC
661	RETURN
662	ENDIF
663
664	IF (FL.GT.FH) THEN
665	FSIGN=1.0
666	ELSE
667	FSIGN=-1.0
668	ENDIF
669
670	XNEW=XM+(XM-XL)(FSIGNFM/S)
671
672	IF (ABS(XNEW-ANS).LE.XACC) THEN
673	IRFRMWV=WSALOC
674	RETURN
675	ENDIF
676
677	ANS=XNEW
678	CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
679	& TAIR,PAIR,RADIUS)
680	FNEW=WVEQ
681
682	IF (FNEW.EQ.0.0) THEN
683	IRFRMWV=WSALOC
684	RETURN
685	ENDIF
686
687	IF (SIGN(FM, FNEW).NE.FM) THEN
688	XL=XM
689	FL=FM
690	XH=ANS
691	FH=FNEW
692	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
693	XH=ANS
694	FH=FNEW
695	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
696	XL=ANS
697	FL=FNEW
698	ELSE
699	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
700	PRINT*,'you shall not PAAAAAASS'
701	STOP
702	ENDIF
703	ENDDO
704	PRINT*,'Paaaaas bien MAXIT atteint'
705	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
706	PRINT*,'you shall not PAAAAAASS'
707	XL=X1
708	XH=X2
709	ANS=-9.99e99
710
711	DO i=1, MAXIT
712	XM=0.5*(XL+XH)
713	CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
714	& TAIR,PAIR,RADIUS)
715	FM=WVEQ
716	S=SQRT(FMFM-FLFH)
717	IF (FL.GT.FH) THEN
718	FSIGN=1.0
719	ELSE
720	FSIGN=-1.0
721	ENDIF
722
723	XNEW=XM+(XM-XL)(FSIGNFM/S)
724
725	ANS=XNEW
726	CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
727	& TAIR,PAIR,RADIUS)
728	FNEW=WVEQ
729	PRINT*,'WVliq',WVLIQ,'WVtot',WVTOT,'WVeq',WVEQ
730	PRINT*,'WSA',WSALOC,'SAtot',SATOT
731	PRINT*,'T',TAIR,'P',PAIR
732
733	IF (SIGN(FM, FNEW).NE.FM) THEN
734	XL=XM
735	FL=FM
736	XH=ANS
737	FH=FNEW
738	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
739	XH=ANS
740	FH=FNEW
741	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
742	XL=ANS
743	FL=FNEW
744	ELSE
745	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
746	PRINT*,'you shall not PAAAAAASS TWIIICE???'
747	STOP
748	ENDIF
749	ENDDO
750	STOP
751	ELSE
752	PRINT*,'IRFRMWV must be bracketed'
753	PRINT*,'NROOT de BRACWV', NROOT
754	IF (ABS(FL).LT.XACC) THEN
755	PRINT*,'IRFRMWV FL == 0',FL
756	PRINT*,'X1',X1,'X2',X2,'FH',FH
757	CALL ITERWV(X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
758	& TAIR,PAIR,RADIUS)
759	IRFRMWV=WSALOC
760	RETURN
761	ENDIF
762	IF (ABS(FH).LT.XACC) THEN
763	PRINT*,'IRFRMWV FH == 0',FH
764	PRINT*,'X1',X1,'X2',X2,'FL',FL
765	CALL ITERWV(X2,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
766	& TAIR,PAIR,RADIUS)
767	IRFRMWV=WSALOC
768	RETURN
769	ENDIF
770	IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN
771	PRINT*,'STOP dans IRFRMWV avec rien == 0'
772	PRINT*,'X1',X1,'X2',X2
773	PRINT*,'Fcalc',FL,FH
774	PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS
775	STOP
776	ENDIF
777	IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN
778	PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC'
779	PRINT*,FL,FH
780	STOP
781	ENDIF
782
783
784	ENDIF
785	! FIN Test Bracketed values
786
787	END FUNCTION IRFRMWV
788
789	!*****************************************************************************
790	!* REAL FUNCTION IRFRMSA()
791	REAL FUNCTION IRFRMSA(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,LPPWV,
792	+ NB)
793	!*****************************************************************************
794	!* Iterative Root Finder Ridder's Method for Sulfuric Acid calculus
795	!* From Numerical Recipes
796	!* Adapted for VenusGCM A. Stolzenbach 07/2014
797	!*
798	!* Les iterations sur [X1,X2] sont [WSA1,WSA2]
799	!* la variable X est WSA
800	!* IRFRMSA sort en OUTPUT : WSA pour KEEQ=0
801
802	IMPLICIT NONE
803
804	REAL, INTENT(IN) :: X1, X2
805	REAL, INTENT(IN) :: XACC
806	INTEGER, INTENT(IN) :: MAXIT, NB
807
808	! LOCAL VARIABLES
809	REAL XL, XH, XM, XNEW
810	REAL Fl, FH, FM, FNEW
811	REAL ANS, S, FSIGN
812	INTEGER i
813
814	! External variables needed:
815	REAL, INTENT(IN) :: TAIR,PAIR
816	REAL, INTENT(IN) :: LPPWV
817	REAL, INTENT(IN) :: RADIUS
818
819	! External functions needed:
820	REAL KEEQ
821
822
823
824	! Initialisation
825	FL=KEEQ(RADIUS,TAIR,X1,LPPWV)
826	FH=KEEQ(RADIUS,TAIR,X2,LPPWV)
827
828	! Test Bracketed values
829	IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.((FL.GT.0.).AND.(FH.LT.0.)))
830	& THEN
831	XL=X1
832	XH=X2
833	ANS=-9.99e99
834
835	DO i=1, MAXIT
836	XM=0.5*(XL+XH)
837	FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
838	S=SQRT(FMFM-FLFH)
839
840	IF (S.EQ.0.0) THEN
841	IRFRMSA=ANS
842	RETURN
843	ENDIF
844
845	IF (FL.GT.FH) THEN
846	FSIGN=1.0
847	ELSE
848	FSIGN=-1.0
849	ENDIF
850
851	XNEW=XM+(XM-XL)(FSIGNFM/S)
852
853	IF (ABS(XNEW-ANS).LE.XACC) THEN
854	IRFRMSA=ANS
855	RETURN
856	ENDIF
857
858	ANS=XNEW
859	FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
860
861	IF (FNEW.EQ.0.0) THEN
862	IRFRMSA=ANS
863	RETURN
864	ENDIF
865
866	IF (SIGN(FM, FNEW).NE.FM) THEN
867	XL=XM
868	FL=FM
869	XH=ANS
870	FH=FNEW
871	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
872	XH=ANS
873	FH=FNEW
874	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
875	XL=ANS
876	FL=FNEW
877	ELSE
878	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
879	PRINT*,'you shall not PAAAAAASS'
880	STOP
881	ENDIF
882	ENDDO
883	PRINT*,'Paaaaas bien MAXIT atteint'
884	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
885	PRINT*,'you shall not PAAAAAASS'
886	XL=X1
887	XH=X2
888	PRINT*,'Borne XL',XL,'XH',XH
889	ANS=-9.99e99
890
891	DO i=1, MAXIT
892	XM=0.5*(XL+XH)
893	FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
894	S=SQRT(FMFM-FLFH)
895
896	IF (FL.GT.FH) THEN
897	FSIGN=1.0
898	ELSE
899	FSIGN=-1.0
900	ENDIF
901
902	XNEW=XM+(XM-XL)(FSIGNFM/S)
903
904	ANS=XNEW
905	FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
906	PRINT*,'KEEQ result',FNEW,'T',TAIR,'R',RADIUS
907	IF (SIGN(FM, FNEW).NE.FM) THEN
908	XL=XM
909	FL=FM
910	XH=ANS
911	FH=FNEW
912	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
913	XH=ANS
914	FH=FNEW
915	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
916	XL=ANS
917	FL=FNEW
918	ELSE
919	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
920	PRINT*,'you shall not PAAAAAASS'
921	STOP
922	ENDIF
923	ENDDO
924	STOP
925	ELSE
926	PRINT*,'IRFRMSA must be bracketed'
927	IF (FL.EQ.0.0) THEN
928	PRINT*,'IRFRMSA FL == 0',Fl
929	IRFRMSA=X1
930	RETURN
931	ENDIF
932	IF (FH.EQ.0.0) THEN
933	PRINT*,'IRFRMSA FH == 0',FH
934	IRFRMSA=X2
935	RETURN
936	ENDIF
937	IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN
938	PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH
939	PRINT*,'X1',X1,'X2',X2
940	PRINT*,'Kind F', KIND(FL), KIND(FH)
941	PRINT*,'Kind X', KIND(X1), KIND(X2)
942	PRINT*,'Logical: ',(SIGN(FL,FH).NE.FL)
943	PRINT*,'Logical: ',(SIGN(FH,FL).NE.FH)
944	PRINT*,'nb root BRACWSA',NB
945	STOP
946	ENDIF
947
948	ENDIF
949	! FIN Test Bracketed values
950
951	END function IRFRMSA
952
953	!*****************************************************************************
954	!* REAL FUNCTION KEEQ()
955	REAL FUNCTION KEEQ(RADIUS,TAIR,WSA,LPPWV)
956	!*****************************************************************************
957	!* Kelvin Equation EQuality
958	!* ln(PPWV_eq) - (2Mh2o sigma)/(R T r rho) - ln(ph2osa) = 0
959	!*
960
961	IMPLICIT NONE
962
963	REAL, INTENT(IN) :: RADIUS,TAIR,WSA,LPPWV
964
965	! Physical constants:
966	REAL MH2O
967	REAL RGAS
968	PARAMETER(
969	! Molar weight of water (kg/mole)
970	+ MH2O=18.0153d-3,
971	! Universal gas constant (J/(mole K))
972	+ RGAS=8.314462175d0)
973	!
974	! External functions needed:
975	REAL PWVSAS_GV,SIGMADROPLET,RHODROPLET
976	! PWVSAS_GV: Natural logaritm of water vapor pressure over
977	! sulfuric acid solution
978	! SIGMADROPLET: Surface tension of sulfuric acid solution
979	! RHODROPLET: Density of sulfuric acid solution
980	!
981	! Auxiliary local variables:
982	REAL C1
983
984	PARAMETER(
985	+ C1=2.0d0*MH2O/RGAS)
986
987
988	KEEQ=LPPWV-C1*SIGMADROPLET(WSA,TAIR)/
989	& (TAIRRADIUSRHODROPLET(WSA,TAIR))-
990	& PWVSAS_GV(TAIR,WSA)
991
992	END FUNCTION KEEQ
993
994	*****************************************************************************
995	* REAL FUNCTION PWVSAS_GV(TAIR,WSA)
996	REAL FUNCTION PWVSAS_GV(TAIR,WSA)
997	*****************************************************************************
998	*
999	* Natural logaritm of saturated water vapor pressure over plane
1000	* sulfuric acid solution.
1001	*
1002	* Source: J.I.Gmitro & T.Vermeulen: A.I.Ch.E.J. 10,740,1964.
1003	* W.F.Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1004	*
1005	* The formula of Gmitro & Vermeulen for saturation pressure
1006	* is used:
1007	* ln(p) = A ln(298/T) + B/T + C + DT
1008	* with values of A,B,C and D given by Gmitro & Vermeulen,
1009	* and calculated from partial molal properties given by Giauque et al.
1010	*
1011	*
1012	*
1013	* Input: TAIR: Temperature (K)
1014	* WSA: Weight fraction of H2SO4 [0;1]
1015	* Output: Natural logaritm of water vapor pressure
1016	* over sulfuric acid solution ( ln(Pa) )
1017	*
1018	*
1019	* External functions needed for calculation of partial molal
1020	* properties of pure components at 25 ! as function of W.
1021	IMPLICIT NONE
1022
1023	REAL :: CPH2O,ALH2O,FFH2O,LH2O
1024	* CPH2O: Partial molal heat capacity of sulfuric acid solution.
1025	* ALH2O: Temparature derivative of CPH2O
1026	* FFH2O: Partial molal free energy of sulfuric acid solution.
1027	* LH2O: Partial molal enthalpy of sulfuric acid
1028	*
1029	!
1030	!
1031	REAL, INTENT(IN) :: TAIR,WSA
1032	REAL :: ADOT,BDOT,CDOT,DDOT
1033	REAL :: RGAS,MMHGPA
1034	REAL :: K1,K2
1035	REAL :: A,B,C,D,CP,L,F,ALFA
1036	! Physical constants given by Gmitro & Vermeulen:
1037	PARAMETER(
1038	+ ADOT=-3.67340,
1039	+ BDOT=-4143.5,
1040	+ CDOT=10.24353,
1041	+ DDOT=0.618943d-3)
1042	PARAMETER(
1043	! Gas constant (cal/(deg mole)):
1044	+ RGAS=1.98726,
1045	! Natural logarith of conversion factor between atm. and Pa:
1046	+ MMHGPA=11.52608845,
1047	+ K1=298.15,
1048	+ K2=K1*K1/2.0)
1049	!
1050	!
1051	CP=CPH2O(WSA)
1052	F=-FFH2O(WSA)
1053	L=-LH2O(WSA)
1054	ALFA=ALH2O(WSA)
1055	!
1056	A=ADOT+(CP-K1*ALFA)/RGAS
1057	B=BDOT+(L-K1CP+K2ALFA)/RGAS
1058	C=CDOT+(CP+(F-L)/K1)/RGAS
1059	D=DDOT-ALFA/(2.0d0*RGAS)
1060	!
1061	! WRITE(,) 'TAIR= ',TAIR,' WSA= ',WSA
1062	! WRITE(,) 'CPH2O(WSA)= ',CP
1063	! WRITE(,) 'ALFAH2O(WSA)= ',ALFA
1064	! WRITE(,) 'FFH2O(WSA)= ',F
1065	! WRITE(,) 'LH2O(WSA)= ',L
1066	!
1067	PWVSAS_GV=ADLOG(K1/TAIR)+B/TAIR+C+DTAIR+MMHGPA
1068
1069	END FUNCTION PWVSAS_GV
1070	*******************************************************************************
1071	* REAL FUNCTION CPH2O(W)
1072	REAL FUNCTION CPH2O(W)
1073	*******************************************************************************
1074	*
1075	* Relative partial molal heat capacity of water (cal/(deg mole) in
1076	* sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
1077	* calculated by cubic spline fitting.
1078	*
1079	* Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1080	*
1081	IMPLICIT NONE
1082
1083	INTEGER :: NPOINT,I
1084	PARAMETER(NPOINT=109)
1085	REAL, DIMENSION(NPOINT) :: WTAB(NPOINT),CPHTAB(NPOINT),
1086	+ Y2(NPOINT),YWORK(NPOINT)
1087	REAL, INTENT(IN):: W
1088	REAL :: CPH
1089	LOGICAL :: FIRST
1090	DATA (WTAB(I),I=1,NPOINT)/
1091	+0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
1092	+0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
1093	+0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
1094	+0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
1095	+0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
1096	+0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
1097	+0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
1098	+0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
1099	+0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
1100	+0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
1101	+0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
1102	+0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
1103	+0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
1104	+0.99908,0.99927,0.99945,0.99963,0.99982/
1105	DATA (CPHTAB(I),I=1,NPOINT)/
1106	+ 17.996, 17.896, 17.875, 17.858, 17.840, 17.820, 17.800, 17.791,
1107	+ 17.783, 17.777, 17.771, 17.769, 17.806, 17.891, 18.057, 18.248,
1108	+ 18.429, 18.567, 18.613, 18.640, 18.660, 18.660, 18.642, 18.592,
1109	+ 18.544, 18.468, 18.348, 18.187, 17.995, 17.782, 17.562, 17.352,
1110	+ 17.162, 16.993, 16.829, 16.657, 16.581, 16.497, 16.405, 16.302,
1111	+ 16.186, 16.053, 15.901, 15.730, 15.540, 15.329, 15.101, 14.853,
1112	+ 14.586, 14.296, 13.980, 13.638, 13.274, 12.896, 12.507, 12.111,
1113	+ 11.911, 11.711, 11.514, 11.320, 11.130, 10.940, 10.760, 10.570,
1114	+ 10.390, 10.200, 10.000, 9.8400, 9.7600, 9.7900, 9.9500, 10.310,
1115	+ 10.950, 11.960, 13.370, 15.060, 16.860, 18.550, 20.000, 21.170,
1116	+ 22.030, 22.570, 22.800, 22.750, 22.420, 21.850, 21.120, 20.280,
1117	+ 19.360, 18.350, 17.220, 15.940, 14.490, 12.840, 10.800, 9.8000,
1118	+ 7.8000, 3.8000,0.20000,-5.4000,-7.0000,-8.8000,-10.900,-13.500,
1119	+-17.000,-22.000,-29.000,-40.000,-59.000/
1120	DATA FIRST/.TRUE./
1121	SAVE FIRST,WTAB,CPHTAB,Y2
1122	!
1123	IF(FIRST) THEN
1124	FIRST=.FALSE.
1125	CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2)
1126	ENDIF
1127	CALL SPLINT(WTAB,CPHTAB,Y2,NPOINT,W,CPH)
1128	CPH2O=CPH
1129
1130	END FUNCTION CPH2O
1131	!
1132	*******************************************************************************
1133	REAL FUNCTION FFH2O(W)
1134	* REAL FUNCTION FFH2O(W)
1135	*******************************************************************************
1136	*
1137	* Relative partial molal free energy water (cal/mole) in
1138	* sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
1139	* calculated by cubic spline fitting.
1140	*
1141	* Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1142	*
1143	IMPLICIT NONE
1144
1145	INTEGER :: NPOINT,I
1146	PARAMETER(NPOINT=110)
1147	REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
1148	REAL, INTENT(IN) :: W
1149	REAL :: FF
1150	LOGICAL :: FIRST
1151	DATA (WTAB(I),I=1,NPOINT)/
1152	+0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
1153	+0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
1154	+0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
1155	+0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
1156	+0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
1157	+0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
1158	+0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
1159	+0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
1160	+0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
1161	+0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
1162	+0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
1163	+0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
1164	+0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
1165	+0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1166	DATA (FFTAB(I),I=1,NPOINT)/
1167	+0.00000, 22.840, 25.810, 29.250, 33.790, 39.970, 48.690, 54.560,
1168	+ 61.990, 71.790, 85.040, 103.70, 130.70, 145.20, 163.00, 184.50,
1169	+ 211.50, 245.60, 266.40, 290.10, 317.40, 349.00, 385.60, 428.40,
1170	+ 452.50, 478.80, 507.50, 538.80, 573.30, 611.60, 653.70, 700.50,
1171	+ 752.60, 810.60, 875.60, 948.60, 980.60, 1014.3, 1049.7, 1087.1,
1172	+ 1126.7, 1168.7, 1213.5, 1261.2, 1312.0, 1366.2, 1424.3, 1486.0,
1173	+ 1551.8, 1622.3, 1697.8, 1778.5, 1864.9, 1956.8, 2055.8, 2162.0,
1174	+ 2218.0, 2276.0, 2337.0, 2400.0, 2466.0, 2535.0, 2607.0, 2682.0,
1175	+ 2760.0, 2842.0, 2928.0, 3018.0, 3111.0, 3209.0, 3311.0, 3417.0,
1176	+ 3527.0, 3640.0, 3757.0, 3878.0, 4002.0, 4130.0, 4262.0, 4397.0,
1177	+ 4535.0, 4678.0, 4824.0, 4973.0, 5128.0, 5287.0, 5454.0, 5630.0,
1178	+ 5820.0, 6031.0, 6268.0, 6541.0, 6873.0, 7318.0, 8054.0, 8284.0,
1179	+ 8579.0, 8997.0, 9295.0, 9720.0, 9831.0, 9954.0, 10092., 10248.,
1180	+ 10423., 10618., 10838., 11099., 11460., 12014./
1181	DATA FIRST/.TRUE./
1182	SAVE FIRST,WTAB,FFTAB,Y2
1183	!
1184	IF(FIRST) THEN
1185	FIRST=.FALSE.
1186	CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2)
1187	ENDIF
1188	CALL SPLINT(WTAB,FFTAB,Y2,NPOINT,W,FF)
1189	FFH2O=FF
1190
1191	END FUNCTION FFH2O
1192	!
1193	*******************************************************************************
1194	REAL FUNCTION LH2O(W)
1195	* REAL FUNCTION LH2O(W)
1196	*******************************************************************************
1197	*
1198	* Relative partial molal heat content of water (cal/mole) in
1199	* sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
1200	* calculated by cubic spline fitting.
1201	*
1202	* Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1203	*
1204	IMPLICIT NONE
1205
1206	INTEGER :: NPOINT,I
1207	PARAMETER(NPOINT=110)
1208	REAL, DIMENSION(NPOINT) :: WTAB,LTAB,Y2,YWORK
1209	REAL, INTENT(IN) :: W
1210	REAL :: L
1211	LOGICAL :: FIRST
1212	DATA (WTAB(I),I=1,NPOINT)/
1213	+0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
1214	+0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
1215	+0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
1216	+0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
1217	+0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
1218	+0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
1219	+0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
1220	+0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
1221	+0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
1222	+0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
1223	+0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
1224	+0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
1225	+0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
1226	+0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1227	DATA (LTAB(I),I=1,NPOINT)/
1228	+0.00000, 5.2900, 6.1000, 7.1800, 8.7800, 11.210, 15.290, 18.680,
1229	+ 23.700, 31.180, 42.500, 59.900, 89.200, 106.70, 128.60, 156.00,
1230	+ 190.40, 233.80, 260.10, 290.00, 324.00, 362.50, 406.50, 456.10,
1231	+ 483.20, 512.40, 543.60, 577.40, 613.80, 653.50, 696.70, 744.50,
1232	+ 797.20, 855.80, 921.70, 995.70, 1028.1, 1062.3, 1098.3, 1136.4,
1233	+ 1176.7, 1219.3, 1264.7, 1313.0, 1364.3, 1418.9, 1477.3, 1539.9,
1234	+ 1607.2, 1679.7, 1757.9, 1842.7, 1934.8, 2035.4, 2145.5, 2267.0,
1235	+ 2332.0, 2401.0, 2473.0, 2550.0, 2631.0, 2716.0, 2807.0, 2904.0,
1236	+ 3007.0, 3118.0, 3238.0, 3367.0, 3507.0, 3657.0, 3821.0, 3997.0,
1237	+ 4186.0, 4387.0, 4599.0, 4819.0, 5039.0, 5258.0, 5476.0, 5694.0,
1238	+ 5906.0, 6103.0, 6275.0, 6434.0, 6592.0, 6743.0, 6880.0, 7008.0,
1239	+ 7133.0, 7255.0, 7376.0, 7497.0, 7618.0, 7739.0, 7855.0, 7876.0,
1240	+ 7905.0, 7985.0, 8110.0, 8415.0, 8515.0, 8655.0, 8835.0, 9125.0,
1241	+ 9575.0, 10325., 11575., 13500., 15200., 16125./
1242	DATA FIRST/.TRUE./
1243	SAVE FIRST,WTAB,LTAB,Y2
1244	!
1245	IF(FIRST) THEN
1246	FIRST=.FALSE.
1247	CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2)
1248	ENDIF
1249	CALL SPLINT(WTAB,LTAB,Y2,NPOINT,W,L)
1250	LH2O=L
1251
1252	END FUNCTION LH2O
1253	*******************************************************************************
1254	REAL FUNCTION ALH2O(W)
1255	* REAL FUNCTION ALH2O(W)
1256	*******************************************************************************
1257	*
1258	* Relative partial molal temperature derivative of heat capacity (water)
1259	* in sulfuric acid solution, (cal/deg**2), calculated by
1260	* cubic spline fitting.
1261	*
1262	* Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1263	*
1264	IMPLICIT NONE
1265
1266	INTEGER :: NPOINT,I
1267	PARAMETER(NPOINT=96)
1268	REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
1269	REAL, INTENT(IN) :: W
1270	REAL :: A
1271	LOGICAL :: FIRST
1272	DATA (WTAB(I),I=1,NPOINT)/
1273	+0.29517,0.31209,
1274	+0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
1275	+0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
1276	+0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
1277	+0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
1278	+0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
1279	+0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
1280	+0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
1281	+0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
1282	+0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
1283	+0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
1284	+0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
1285	+0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1286	DATA (ATAB(I),I=1,NPOINT)/
1287	+ 0.0190, 0.0182, 0.0180, 0.0177, 0.0174, 0.0169, 0.0167, 0.0164,
1288	+ 0.0172, 0.0212, 0.0239, 0.0264, 0.0276, 0.0273, 0.0259, 0.0238,
1289	+ 0.0213, 0.0190, 0.0170, 0.0155, 0.0143, 0.0133, 0.0129, 0.0124,
1290	+ 0.0120, 0.0114, 0.0106, 0.0097, 0.0084, 0.0067, 0.0047, 0.0024,
1291	+-0.0002,-0.0031,-0.0063,-0.0097,-0.0136,-0.0178,-0.0221,-0.0263,
1292	+-0.0303,-0.0340,-0.0352,-0.0360,-0.0362,-0.0356,-0.0343,-0.0321,
1293	+-0.0290,-0.0251,-0.0201,-0.0137,-0.0058, 0.0033, 0.0136, 0.0254,
1294	+ 0.0388, 0.0550, 0.0738, 0.0962, 0.1198, 0.1300, 0.1208, 0.0790,
1295	+ 0.0348, 0.0058,-0.0102,-0.0211,-0.0292,-0.0350,-0.0390,-0.0418,
1296	+-0.0432,-0.0436,-0.0429,-0.0411,-0.0384,-0.0346,-0.0292,-0.0220,
1297	+-0.0130,-0.0110,-0.0080,-0.0060,-0.0040,-0.0030,-0.0030,-0.0020,
1298	+-0.0020,-0.0020,-0.0020,-0.0010,-0.0010, 0.0000, 0.0000, 0.0000/
1299	DATA FIRST/.TRUE./
1300	SAVE FIRST,WTAB,ATAB,Y2
1301	!
1302	IF(FIRST) THEN
1303	FIRST=.FALSE.
1304	CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2)
1305	ENDIF
1306	CALL SPLINT(WTAB,ATAB,Y2,NPOINT,MAX(WTAB(1),W),A)
1307	ALH2O=A
1308
1309	END FUNCTION ALH2O
1310	!******************************************************************************
1311	SUBROUTINE SPLINE(X,Y,N,WORK,Y2)
1312	!******************************************************************************
1313	! Routine to calculate 2.nd derivatives of tabulated function
1314	! Y(i)=Y(Xi), to be used for cubic spline calculation.
1315	!
1316	IMPLICIT NONE
1317
1318	INTEGER, INTENT(IN) :: N
1319	INTEGER :: I
1320	REAL, DIMENSION(N), INTENT(IN) :: X,Y
1321	REAL, DIMENSION(N), INTENT(OUT) :: Y2,WORK
1322	REAL SIG,P,QN,UN,YP1,YPN
1323
1324	!AM Venus: Let's check the values
1325	! write(,) 'In spline, N ', N
1326
1327	YP1=(Y(2)-Y(1))/(X(2)-X(1))
1328	YPN=(Y(N)-Y(N-1))/(X(N)-X(N-1))
1329	IF(YP1.GT.99.0E+30) THEN
1330	Y2(1)=0.0
1331	WORK(1)=0.0
1332	ELSE
1333	Y2(1)=-0.5d0
1334	WORK(1)=(3.0d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
1335	ENDIF
1336	DO I=2,N-1
1337	! write(,) 'In spline, I ', I
1338	SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
1339	P=SIG*Y2(I-1)+2.0d0
1340	Y2(I)=(SIG-1.0d0)/P
1341	WORK(I)=(6.0d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
1342	+ /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*WORK(I-1))/P
1343	ENDDO
1344	IF(YPN.GT.99.0E+30) THEN
1345	QN=0.0
1346	UN=0.0
1347	ELSE
1348	QN=0.5d0
1349	UN=(3.0d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
1350	ENDIF
1351	Y2(N)=(UN-QNWORK(N-1))/(QNY2(N-1)+1.0d0)
1352	DO I=N-1,1,-1
1353	! write(,) 'In spline, I ', I
1354	Y2(I)=Y2(I)*Y2(I+1)+WORK(I)
1355	ENDDO
1356	!
1357	END SUBROUTINE SPLINE
1358
1359	!******************************************************************************
1360	SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
1361	!******************************************************************************
1362	! Cubic spline calculation
1363
1364	IMPLICIT NONE
1365
1366	INTEGER, INTENT(IN) :: N
1367	INTEGER :: KLO,KHI,K
1368	REAL, INTENT(IN), DIMENSION(N) :: XA,YA,Y2A
1369	REAL, INTENT(IN) :: X
1370	REAL, INTENT(OUT) :: Y
1371	REAL :: H,A,B
1372	!
1373	KLO=1
1374	KHI=N
1375	1 IF(KHI-KLO.GT.1) THEN
1376	K=(KHI+KLO)/2
1377	IF(XA(K).GT.X) THEN
1378	KHI=K
1379	ELSE
1380	KLO=K
1381	ENDIF
1382	GOTO 1
1383	ENDIF
1384	H=XA(KHI)-XA(KLO)
1385	A=(XA(KHI)-X)/H
1386	B=(X-XA(KLO))/H
1387	Y=AYA(KLO)+BYA(KHI)+
1388	+ ((A*3-A)Y2A(KLO)+(B*3-B)Y2A(KHI))(H*2)/6.0d0
1389	!
1390
1391	END SUBROUTINE SPLINT
1392	!******************************************************************
1393	SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
1394	+ T,H2SO4COND,H2OCOND,RMTOT)
1395
1396	! DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
1397	! FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
1398	! SIZE DISTRIBTUION
1399	! ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
1400	! ---------------------------------------------------------------
1401	! INPUT:
1402	! H2SO4: #/m3 of total H2SO4
1403	! H2O : #/m3 of total H2O
1404	! WSA: aerosol H2SO4 weight fraction (fraction)
1405	! DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume)
1406	! for total mass, almost same result with ro=1.67 gr/cm3
1407	! RSTDEV: standard deviation of aerosol distribution (no unit)
1408	! RADIUS: MEDIAN radius (m)
1409	! T: temperature (K)
1410	!
1411	! OUTPUT:
1412	! RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
1413	! mais rho_droplet et M_tot_distrib doivent tre de meme dimension
1414	! H2OCOND
1415	! H2SO4COND
1416
1417
1418
1419	IMPLICIT NONE
1420
1421	REAL, INTENT(IN) :: H2SO4, H2O, WSA
1422	REAL, INTENT(IN) :: DENSO4, T
1423	REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
1424	! working variables
1425	REAL :: RMH2S4
1426	REAL :: DND2,pstand,lpar,acidps
1427	REAL :: x1, satpacid
1428	REAL , DIMENSION(2):: act
1429	!
1430	! masse of an H2SO4 molecule (kg)
1431	RMH2S4=98.078/(6.02214129E+26)
1432
1433	pstand=1.01325E+5 !Pa 1 atm pressure
1434
1435	x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
1436
1437	call zeleznik(x1,t,act)
1438
1439	!pure acid satur vapor pressure
1440	lpar= -11.695+DLOG(pstand) ! Zeleznik
1441	acidps=1/360.15-1.0/t+0.38/545.
1442	+ *(1.0+DLOG(360.15/t)-360.15/t)
1443	acidps = 10156.0*acidps +lpar
1444	acidps = DEXP(acidps) !Pa
1445
1446	!acid sat.vap.PP over mixture (flat surface):
1447	satpacid=act(2)*acidps ! Pa
1448
1449	! Conversion from Pa to N.D #/m3
1450	DND2=satpacid/(1.3806488E-23*T)
1451	! Conversion from N.D #/m3 TO #/cm3
1452	! DND2=DND2*1.d-6
1453
1454	! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
1455	IF (H2SO4.GE.DND2) THEN
1456	H2SO4COND=H2SO4-DND2
1457	! calcul de H2O cond correspondant a H2SO4 cond
1458	H2OCOND=H2SO4COND98.078(1.0-WSA)/(18.0153*WSA)
1459
1460	! RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3
1461	! RMTOT=M_distrib/rho_droplet
1462
1463	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
1464
1465	! Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
1466	ELSE
1467	H2SO4COND=0.0E+0
1468	H2OCOND=0.0E+0
1469	RMTOT=0.0E+0
1470	END IF
1471
1472	! Test si H2O en defaut H2Ocond>H2O dispo
1473	IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
1474
1475	! Si H2O en dfaut, on as pas le bon WSA!
1476	! En effet, normalement, on a exactement le WSA correspondant a
1477	! WVg + WVl = WVtot
1478	! Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou
1479	! quand T est grand (Lower Haze), le modle reprsente mal le WSA
1480	! cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA
1481
1482	! PRINT*,'PROBLEM H2O EN DEFAUT'
1483	! PRINT*,'H2OCOND',H2OCOND,'H2O',H2O
1484	! PRINT*,'WSA',WSA,'RHO',DENSO4
1485	! STOP
1486
1487
1488	! On peut alors condenser tout le H2O dispo
1489	H2OCOND=H2O
1490	! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
1491	H2SO4COND=H2OCOND18.0153WSA/(98.078*(1.0-WSA))
1492
1493	! RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3
1494	! RMTOT=Volume of aerosol cm3 /cm3 of air
1495	! Volume of aerosol/cm3 air
1496
1497	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
1498
1499	END IF
1500
1501	END SUBROUTINE CALCM_SAT
1502
1503	SUBROUTINE Zeleznik(x,T,act)
1504
1505	!+++++++++++++++++++++++++++++++++++++++++++++++++++
1506	! Water and sulfuric acid activities in liquid
1507	! aqueous solutions.
1508	! Frank J. Zeleznik, Thermodynnamic properties
1509	! of the aqueous sulfuric acid system to 220K-350K,
1510	! mole fraction 0,...,1
1511	! J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991
1512	!+++++++++++++++++++++++++++++++++++++++++++++++++++
1513
1514	IMPLICIT NONE
1515
1516	REAL, INTENT(IN) :: x,T
1517	REAL :: activitya, activityw
1518	REAL, INTENT(OUT), DIMENSION(2):: act
1519	! REAL x,T, activitya, activityw
1520	! REAL, DIMENSION(2):: act
1521
1522
1523	! write(,) 'x, T ', x, T
1524
1525	act(2)=activitya(x,T)
1526	act(1)=activityw(x,T)
1527
1528	! write(,) 'act ', act
1529
1530	END SUBROUTINE Zeleznik
1531
1532	!start of functions related to zeleznik activities
1533
1534	REAL FUNCTION m111(T)
1535
1536	REAL, INTENT(IN) :: T
1537	m111=-23.524503387D0
1538	& +0.0406889449841D0*T
1539	& -0.151369362907D-4T*2+2961.44445015D0/T
1540	& +0.492476973663D0*dlog(T)
1541	END FUNCTION m111
1542
1543	REAL FUNCTION m121(T)
1544
1545	REAL, INTENT(IN) :: T
1546	m121=1114.58541077D0-1.1833078936D0*T
1547	& -0.00209946114412D0T*2-246749.842271D0/T
1548	& +34.1234558134D0*dlog(T)
1549	END FUNCTION m121
1550
1551	FUNCTION m221(T)
1552
1553	REAL, INTENT(IN) :: T
1554	m221=-80.1488100747D0-0.0116246143257D0*T
1555	& +0.606767928954D-5T*2+3092.72150882D0/T
1556	& +12.7601667471D0*dlog(T)
1557	END FUNCTION m221
1558
1559	REAL FUNCTION m122(T)
1560
1561	REAL, INTENT(IN) :: T
1562	m122=888.711613784D0-2.50531359687D0*T
1563	& +0.000605638824061D0T*2-196985.296431D0/T
1564	& +74.550064338D0*dlog(T)
1565	END FUNCTION m122
1566
1567	REAL FUNCTION e111(T)
1568
1569	REAL, INTENT(IN) :: T
1570	e111=2887.31663295D0-3.32602457749D0*T
1571	& -0.2820472833D-2T*2-528216.112353D0/T
1572	& +0.68699743564D0*dlog(T)
1573	END FUNCTION e111
1574
1575	REAL FUNCTION e121(T)
1576
1577	REAL, INTENT(IN) :: T
1578	e121=-370.944593249D0-0.690310834523D0*T
1579	& +0.56345508422D-3T*2-3822.52997064D0/T
1580	& +94.2682037574D0*dlog(T)
1581	END FUNCTION e121
1582
1583	REAL FUNCTION e211(T)
1584
1585	REAL, INTENT(IN) :: T
1586	e211=38.3025318809D0-0.0295997878789D0*T
1587	& +0.120999746782D-4T*2-3246.97498999D0/T
1588	& -3.83566039532D0*dlog(T)
1589	END FUNCTION e211
1590
1591	REAL FUNCTION e221(T)
1592
1593	REAL, INTENT(IN) :: T
1594	e221=2324.76399402D0-0.141626921317D0*T
1595	& -0.00626760562881D0T*2-450590.687961D0/T
1596	& -61.2339472744D0*dlog(T)
1597	END FUNCTION e221
1598
1599	REAL FUNCTION e122(T)
1600
1601	REAL, INTENT(IN) :: T
1602	e122=-1633.85547832D0-3.35344369968D0*T
1603	& +0.00710978119903D0T*2+198200.003569D0/T
1604	& +246.693619189D0*dlog(T)
1605	END FUNCTION e122
1606
1607	REAL FUNCTION e212(T)
1608
1609	REAL, INTENT(IN) :: T
1610	e212=1273.75159848D0+1.03333898148D0*T
1611	& +0.00341400487633D0T*2+195290.667051D0/T
1612	& -431.737442782D0*dlog(T)
1613	END FUNCTION e212
1614
1615	REAL FUNCTION lnAa(x1,T)
1616
1617	REAL, INTENT(IN) :: T,x1
1618	REAL ::
1619	& m111,m121,m221,m122
1620	& ,e111,e121,e211,e122,e212,e221
1621	lnAa=-(
1622	& (2m111(T)+e111(T)(2dlog(x1)+1))x1
1623	& +(2m121(T)+e211(T)dlog(1-x1)+e121(T)(dlog(x1)+1))(1-x1)
1624	& -(m111(T)+e111(T)(dlog(x1)+1))x1*x1
1625	& -(2m121(T)+e121(T)(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1)
1626	& -(2m122(T)+e122(T)dlog(x1)
1627	& +e212(T)dlog(1-x1))(1-x1))x1(1-x1)
1628	& -(m221(T)+e221(T)(dlog(1-x1)+1))(1-x1)**2
1629	& -x1(1-x1)(
1630	& (6m122(T)+e122(T)(3*dlog(x1)+1)
1631	& +e212(T)(3dlog(1-x1)+1)
1632	& )x1(1-x1)
1633	& -(2m122(T)+e122(T)(dlog(x1)+1)
1634	& +e212(T)*dlog(1-x1)
1635	& )*(1-x1))
1636	& )
1637	END FUNCTION lnAa
1638
1639	REAL FUNCTION lnAw(x1,T)
1640
1641	REAL, INTENT(IN) :: T,x1
1642	REAL ::
1643	& m111,m121,m221,m122
1644	& ,e111,e121,e211,e122,e212,e221
1645	lnAw=-(
1646	& (2m121(T)+e121(T)dlog(x1)+e211(T)(dlog(1-x1)+1))x1
1647	& +(2m221(T)+e221(T)(2dlog(1-x1)+1))(1-x1)
1648	& -(m111(T)+e111(T)(dlog(x1)+1))x1*x1
1649	& -(2m121(T)+e121(T)(dlog(x1)+1)
1650	& +e211(T)(dlog(1-x1)+1))x1*(1-x1)
1651	& -(m221(T)+e221(T)(dlog(1-x1)+1))(1-x1)**2
1652	& +x1(2m122(T)+e122(T)dlog(x1)+e212(T)dlog(1-x1))x1(1-x1)
1653	& +x1(1-x1)((2m122(T)+e122(T)dlog(x1)
1654	& +e212(T)(dlog(1-x1)+1))x1
1655	& -(6m122(T)+e122(T)(3*dlog(x1)+1)
1656	& +e212(T)(3dlog(1-x1)+1))(1-x1)x1)
1657	& )
1658	END FUNCTION lnAw
1659
1660	REAL FUNCTION activitya(xal,T)
1661
1662	REAL, INTENT(IN) :: T,xal
1663	REAL :: lnAa
1664	! & ,m111,m121,m221,m122 &
1665	! & ,e111,e121,e211,e122,e212,e221
1666
1667	! write(,) 'in activitya ', xal, T
1668	activitya=DEXP(lnAa(xal,T)-lnAa(1.D0-1.D-12,T))
1669	END FUNCTION activitya
1670
1671	FUNCTION activityw(xal,T)
1672
1673	REAL, INTENT(IN) :: T,xal
1674	REAL :: lnAw
1675
1676	activityw=DEXP(lnAw(xal,T)-lnAw(1.D-12,T))
1677	END FUNCTION activityw
1678
1679	! end of functions related to zeleznik activities
1680
1681
1682
1683
1684	FUNCTION SIGMADROPLET(xmass,t)
1685	! calculates the surface tension of the liquid in J/m^2
1686	! xmass=mass fraction of h2so4, t in kelvins
1687	! about 230-323 K , x=0,...,1
1688	!(valid down to the solid phase limit temp, which depends on molefraction)
1689	IMPLICIT NONE
1690	REAL :: SIGMADROPLET
1691	REAL, INTENT(IN):: xmass, t
1692	REAL :: a, b, t1, tc, xmole
1693	REAL, PARAMETER :: Msa=98.078
1694	REAL, PARAMETER :: Mwv=18.0153
1695
1696	IF (t .LT. 305.15) THEN
1697	!low temperature surface tension
1698	! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari
1699	! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
1700	! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
1701	!and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631
1702	a=0.11864+xmass(-0.11651+xmass(0.76852+xmass*
1703	& (-2.40909+xmass(2.95434-xmass1.25852))))
1704	b=-0.00015709+xmass(0.00040102+xmass(-0.00239950+xmass*
1705	& (0.007611235+xmass(-0.00937386+xmass0.00389722))))
1706	SIGMADROPLET=a+t*b
1707	ELSE
1708
1709	xmole = (xmass/Msa)*(1./((xmass/Msa)+(1.-xmass)/Mwv))
1710	! high temperature surface tension
1711	!H. Vehkam‰ki and M. Kulmala and K.E. J. lehtinen, 2003,
1712	!Modelling binary homogeneous nucleation of water-sulfuric acid vapours:
1713	! parameterisation for high temperature emissions, () Environ. Sci. Technol., 37, 3392-3398
1714
1715	tc= 647.15(1.0-xmole)(1.0-xmole)+900.0xmolexmole+
1716	& 3156.186xmole(1-xmole) !critical temperature
1717	t1=1.0-t/tc
1718	a= 0.2358+xmole(-0.529+xmole(4.073+xmole(-12.6707+xmole
1719	& (15.3552+xmole*(-6.3138)))))
1720	b=-0.14738+xmole(0.6253+xmole(-5.4808+xmole(17.2366+xmole
1721	& (-21.0487+xmole*(8.719)))))
1722	SIGMADROPLET=(a+bt1)t1**(1.256)
1723	END IF
1724
1725	RETURN
1726	END FUNCTION SIGMADROPLET
1727
1728	FUNCTION RHODROPLET(xmass,t)
1729	!
1730	! calculates the density of the liquid in kg/m^3
1731	! xmass=mass fraction of h2so4, t in kelvins
1732	! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari
1733	! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
1734	! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
1735	!and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631
1736
1737	! about 220-373 K , x=0,...,1
1738	!(valid down to the solid phase limit temp, which depends on molefraction)
1739
1740	IMPLICIT NONE
1741	REAL :: RHODROPLET
1742	REAL, INTENT(IN) :: xmass, t
1743	REAL :: a,b,c
1744
1745
1746	a=0.7681724+xmass(2.1847140+xmass(7.1630022+xmass*
1747	& (-44.31447+xmass*
1748	& (88.75606+xmass(-75.73729+xmass23.43228)))))
1749	b=1.808225e-3+xmass(-9.294656e-3+xmass(-0.03742148+
1750	& xmass(0.2565321+xmass(-0.5362872+xmass*
1751	& (0.4857736-xmass*0.1629592)))))
1752	c=-3.478524e-6+xmass(1.335867e-5+xmass
1753	& (5.195706e-5+xmass(-3.717636e-4+xmass
1754	& (7.990811e-4+xmass(-7.458060e-4+xmass2.58139e-4)))))
1755	RHODROPLET=a+t(b+ct) ! g/cm^3
1756	RHODROPLET= RHODROPLET*1.0e3 !kg/m^3
1757	RETURN
1758	END FUNCTION RHODROPLET
1759
1760
1761
1762

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

Download in other formats:

Original Format