Context Navigation

new_cloud_venus.F @ 1453

Last change on this file since 1453 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

Rev	Line
[1305]	1	!******************************************************************************
	2	!* venus_SAS_composition SUBROUTINE
	3	!* modified from
	4	!* PROGRAM PSC_MODEL_E
[1442]	5	!* by A. Määttänen
[1305]	6	!* subroutine for LMDZ+photochemistry VENUS
	7	!* by A. Stolzenbach
	8	!*
	9	!* Input/Output files:
	10	!* -------------------
	11	!*
	12	!----------------------------------------------------------------------------
[1442]	13	SUBROUTINE new_cloud_venus(
	14	+ nblev, nblon,
	15	+ TT,PP,
[1305]	16	+ mrt_wv,mrt_sa,
[1442]	17	+ mr_wv,mr_sa)
[1305]	18
	19	USE chemparam_mod
	20	IMPLICIT NONE
[1442]	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
[1305]	27	!----------------------------------------------------------------------------
	28	! Ambient air state variables:
[1442]	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
[1305]	32	!----------------------------------------------------------------------------
[1442]	33	INTEGER :: ilon, ilev, imode
	34	!----------------------------------------------------------------------------
[1305]	35	! Thermodynamic functions:
[1442]	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
[1305]	49
[1442]	50	INTEGER :: NBROOT
	51
	52	INTEGER :: MAXITE
	53	PARAMETER(MAXITE=20)
	54
	55	INTEGER :: NBRAC
	56	PARAMETER(NBRAC=20)
	57
	58	INTEGER :: FLAG
[1305]	59	!----------------------------------------------------------------------------
[1442]	60
[1305]	61	!----------------------------------------------------------------------------
[1442]	62	! External functions needed:
	63	REAL :: IRFRMWV
	64	!----------------------------------------------------------------------------
[1305]	65
[1442]	66
[1305]	67	! >>> Program starts here:
	68
	69	!AM Venus
[1442]	70	! These aerosols will then be given an equilibrium composition for the given size distribution
[1305]	71
[1442]	72	! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
[1305]	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
[1442]	75	!and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
[1305]	76
[1442]	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
[1305]	101
[1442]	102	! Initialisation des bornes pour WV
	103	WVMIN=1.E-90
	104	WVMAX=mrt_wv(ilon,ilev)
[1305]	105
[1442]	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)
[1305]	128
[1442]	129	rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
	130	& TT(ilon,ilev))
[1305]	131
[1442]	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
[1305]	143
[1442]	144	NH2SO4=mrt_sa(ilon,ilev)*CONCM
	145	NH2O=mrt_wv(ilon,ilev)*CONCM
[1305]	146
[1442]	147	CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
	148	& rho_droplet(ilon,ilev),TT(ilon,ilev),
	149	& H2SO4_liq,H2O_liq,MCONDTOT)
[1305]	150
[1442]	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
[1305]	162
[1442]	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)
[1305]	195
[1442]	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
[1305]	207
[1442]	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
[1305]	240
[1442]	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
[1305]	254
[1442]	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
[1305]	268
[1442]	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
[1305]	284
[1442]	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
[1305]	339	IMPLICIT NONE
	340
[1442]	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
[1305]	355
[1442]	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
[1305]	362
[1442]	363	NROOT=0
[1305]	364
[1442]	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
[1305]	432	!----------------------------------------------------------------------------
[1442]	433	! External functions needed:
	434	REAl KEEQ
	435	!----------------------------------------------------------------------------
[1305]	436
[1442]	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
[1305]	512	ENDIF
[1442]	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
[1305]	531
[1442]	532	! USE chemparam_mod
	533
	534	IMPLICIT NONE
[1305]	535
[1442]	536	REAL, INTENT(IN) :: SAt, WSA
	537	REAL, INTENT(IN) :: T, P
[1305]	538
[1442]	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
[1305]	548
[1442]	549	CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS!
[1305]	550
[1442]	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
[1305]	604	!*****************************************************************************
[1442]	605	!* REAL FUNCTION IRFRMWV()
	606	REAL FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,
	607	+ WVTOT,SATOT,NROOT)
[1305]	608	!*****************************************************************************
[1442]	609	!* Iterative Root Finder Ridder's Method for Water Vapor calculus
	610	!* From Numerical Recipes
	611	!* Adapted for VenusGCM A. Stolzenbach 07/2014
[1305]	612	!*
[1442]	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
[1305]	797	!*
[1442]	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
[1305]	802	IMPLICIT NONE
[1442]	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
[1305]	821
	822
	823
[1442]	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	!*
[1305]	960
	961	IMPLICIT NONE
	962
[1442]	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)
[1305]	973	!
[1442]	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
[1305]	983
[1442]	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
[1305]	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
[1442]	1023	REAL :: CPH2O,ALH2O,FFH2O,LH2O
[1305]	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	!
[1442]	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
[1305]	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
[1442]	1068
	1069	END FUNCTION PWVSAS_GV
[1305]	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
[1442]	1083	INTEGER :: NPOINT,I
[1305]	1084	PARAMETER(NPOINT=109)
[1442]	1085	REAL, DIMENSION(NPOINT) :: WTAB(NPOINT),CPHTAB(NPOINT),
	1086	+ Y2(NPOINT),YWORK(NPOINT)
	1087	REAL, INTENT(IN):: W
	1088	REAL :: CPH
	1089	LOGICAL :: FIRST
[1305]	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
[1442]	1129
	1130	END FUNCTION CPH2O
[1305]	1131	!
	1132	*******************************************************************************
[1442]	1133	REAL FUNCTION FFH2O(W)
[1305]	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
[1442]	1145	INTEGER :: NPOINT,I
[1305]	1146	PARAMETER(NPOINT=110)
[1442]	1147	REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
	1148	REAL, INTENT(IN) :: W
	1149	REAL :: FF
	1150	LOGICAL :: FIRST
[1305]	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
[1442]	1190
	1191	END FUNCTION FFH2O
[1305]	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
[1442]	1206	INTEGER :: NPOINT,I
[1305]	1207	PARAMETER(NPOINT=110)
[1442]	1208	REAL, DIMENSION(NPOINT) :: WTAB,LTAB,Y2,YWORK
	1209	REAL, INTENT(IN) :: W
	1210	REAL :: L
	1211	LOGICAL :: FIRST
[1305]	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
[1442]	1251
	1252	END FUNCTION LH2O
[1305]	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
[1442]	1266	INTEGER :: NPOINT,I
[1305]	1267	PARAMETER(NPOINT=96)
[1442]	1268	REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
	1269	REAL, INTENT(IN) :: W
	1270	REAL :: A
	1271	LOGICAL :: FIRST
[1305]	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
[1442]	1308
	1309	END FUNCTION ALH2O
[1305]	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
[1442]	1318	INTEGER, INTENT(IN) :: N
	1319	INTEGER :: I
	1320	REAL, DIMENSION(N), INTENT(IN) :: X,Y
	1321	REAL, DIMENSION(N), INTENT(OUT) :: Y2,WORK
[1305]	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	!
[1442]	1357	END SUBROUTINE SPLINE
[1305]	1358
	1359	!******************************************************************************
	1360	SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
	1361	!******************************************************************************
	1362	! Cubic spline calculation
[1442]	1363
[1305]	1364	IMPLICIT NONE
	1365
[1442]	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
[1305]	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	!
[1442]	1390
	1391	END SUBROUTINE SPLINT
[1305]	1392	!******************************************************************
[1442]	1393	SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
	1394	+ T,H2SO4COND,H2OCOND,RMTOT)
[1305]	1395
	1396	! DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
	1397	! FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
	1398	! SIZE DISTRIBTUION
[1442]	1399	! ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
[1305]	1400	! ---------------------------------------------------------------
	1401	! INPUT:
[1442]	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
[1305]	1407	! RSTDEV: standard deviation of aerosol distribution (no unit)
[1442]	1408	! RADIUS: MEDIAN radius (m)
[1305]	1409	! T: temperature (K)
	1410	!
	1411	! OUTPUT:
[1442]	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
[1305]	1414	! H2OCOND
	1415	! H2SO4COND
	1416
[1442]	1417
	1418
[1305]	1419	IMPLICIT NONE
	1420
[1442]	1421	REAL, INTENT(IN) :: H2SO4, H2O, WSA
	1422	REAL, INTENT(IN) :: DENSO4, T
	1423	REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
[1305]	1424	! working variables
[1442]	1425	REAL :: RMH2S4
	1426	REAL :: DND2,pstand,lpar,acidps
	1427	REAL :: x1, satpacid
[1305]	1428	REAL , DIMENSION(2):: act
	1429	!
[1442]	1430	! masse of an H2SO4 molecule (kg)
	1431	RMH2S4=98.078/(6.02214129E+26)
[1305]	1432
[1442]	1433	pstand=1.01325E+5 !Pa 1 atm pressure
[1305]	1434
[1442]	1435	x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
[1305]	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
[1442]	1446	!acid sat.vap.PP over mixture (flat surface):
[1305]	1447	satpacid=act(2)*acidps ! Pa
	1448
	1449	! Conversion from Pa to N.D #/m3
[1442]	1450	DND2=satpacid/(1.3806488E-23*T)
[1305]	1451	! Conversion from N.D #/m3 TO #/cm3
[1442]	1452	! DND2=DND2*1.d-6
[1305]	1453
[1442]	1454	! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
[1305]	1455	IF (H2SO4.GE.DND2) THEN
	1456	H2SO4COND=H2SO4-DND2
	1457	! calcul de H2O cond correspondant a H2SO4 cond
[1442]	1458	H2OCOND=H2SO4COND98.078(1.0-WSA)/(18.0153*WSA)
[1305]	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
[1442]	1463	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
[1305]	1464
	1465	! Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
	1466	ELSE
[1442]	1467	H2SO4COND=0.0E+0
	1468	H2OCOND=0.0E+0
	1469	RMTOT=0.0E+0
[1305]	1470	END IF
	1471
	1472	! Test si H2O en defaut H2Ocond>H2O dispo
[1442]	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
[1305]	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
[1442]	1491	H2SO4COND=H2OCOND18.0153WSA/(98.078*(1.0-WSA))
[1305]	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
[1442]	1497	RMTOT=H2SO4CONDRMH2S4/(DENSO4WSA)
	1498
[1305]	1499	END IF
[1442]	1500
	1501	END SUBROUTINE CALCM_SAT
[1305]	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
[1442]	1511	! J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991
[1305]	1512	!+++++++++++++++++++++++++++++++++++++++++++++++++++
	1513
	1514	IMPLICIT NONE
	1515
[1442]	1516	REAL, INTENT(IN) :: x,T
	1517	REAL :: activitya, activityw
	1518	REAL, INTENT(OUT), DIMENSION(2):: act
[1305]	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
[1442]	1534	REAL FUNCTION m111(T)
[1305]	1535
[1442]	1536	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1543	REAL FUNCTION m121(T)
[1305]	1544
[1442]	1545	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1553	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1559	REAL FUNCTION m122(T)
[1305]	1560
[1442]	1561	REAL, INTENT(IN) :: T
[1305]	1562	m122=888.711613784D0-2.50531359687D0*T
	1563	& +0.000605638824061D0T*2-196985.296431D0/T
	1564	& +74.550064338D0*dlog(T)
	1565	END FUNCTION m122
	1566
[1442]	1567	REAL FUNCTION e111(T)
[1305]	1568
[1442]	1569	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1575	REAL FUNCTION e121(T)
[1305]	1576
[1442]	1577	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1583	REAL FUNCTION e211(T)
[1305]	1584
[1442]	1585	REAL, INTENT(IN) :: T
[1305]	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
[1442]	1591	REAL FUNCTION e221(T)
[1305]	1592
[1442]	1593	REAL, INTENT(IN) :: T
[1305]	1594	e221=2324.76399402D0-0.141626921317D0*T
	1595	& -0.00626760562881D0T*2-450590.687961D0/T
	1596	& -61.2339472744D0*dlog(T)
	1597	END FUNCTION e221
	1598
[1442]	1599	REAL FUNCTION e122(T)
[1305]	1600
[1442]	1601	REAL, INTENT(IN) :: T
[1305]	1602	e122=-1633.85547832D0-3.35344369968D0*T
	1603	& +0.00710978119903D0T*2+198200.003569D0/T
	1604	& +246.693619189D0*dlog(T)
	1605	END FUNCTION e122
	1606
[1442]	1607	REAL FUNCTION e212(T)
[1305]	1608
[1442]	1609	REAL, INTENT(IN) :: T
[1305]	1610	e212=1273.75159848D0+1.03333898148D0*T
	1611	& +0.00341400487633D0T*2+195290.667051D0/T
	1612	& -431.737442782D0*dlog(T)
	1613	END FUNCTION e212
	1614
[1442]	1615	REAL FUNCTION lnAa(x1,T)
[1305]	1616
[1442]	1617	REAL, INTENT(IN) :: T,x1
	1618	REAL ::
	1619	& m111,m121,m221,m122
[1305]	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
[1442]	1639	REAL FUNCTION lnAw(x1,T)
[1305]	1640
[1442]	1641	REAL, INTENT(IN) :: T,x1
	1642	REAL ::
	1643	& m111,m121,m221,m122
[1305]	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
[1442]	1660	REAL FUNCTION activitya(xal,T)
[1305]	1661
[1442]	1662	REAL, INTENT(IN) :: T,xal
	1663	REAL :: lnAa
[1305]	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
[1442]	1673	REAL, INTENT(IN) :: T,xal
	1674	REAL :: lnAw
	1675
[1305]	1676	activityw=DEXP(lnAw(xal,T)-lnAw(1.D-12,T))
	1677	END FUNCTION activityw
	1678
	1679	! end of functions related to zeleznik activities
[1442]	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