Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

wsa_new.F90

Last change on this file was 1674, checked in by slebonnois, 8 years ago
SL: correction of a bug in microphysics
File size: 39.1 KB

Line
1
2	! SUBROUTINE WSA_ROSA_NEW
3	! SUBROUTINE ITERWV WSA pour un WV
4	! SUBROUTINE BRACWV Bracket de ITERWV
5	! SUBROUTINE BRACWSA Bracket de KEEQ
6	! FUNCTION IRFRMWV Iterative Root Finder Ridder's Method for WV
7	! FUNCTION IRFRMSA Iterative Root Finder Ridder's Method for SA
8	! FUNCTION KEEQ Kelvin Equation EQuality
9	! FUNCTION WVCOND H2O Condensation with WSA, T, P and H2SO4tot
10
11
12	!----------------------------------------------------------------------------
13	SUBROUTINE WSA_ROSA_NEW(TAIR,PAIR,RADIUS,WSAS,MSAD)
14
15	!* This subroutine calculates the acid mass fraction, density, and
16	!* mass of sulfuric acid in a single aerosol droplet of a specified
17	!* radius in equilibrium with ambient water vapor partial pressure
18	!* and temperature.
19	!*
20	!* The calculation is performed by iteration of
21	!* ln(PPWV) - [(2Mh2o sigma)/(R T r rho) - ln(ph2osa)] = 0
22	!* using the secant method. Vapor pressures by Gmitro and Vermeulen
23	!* (PWVSAS_GV) are used.
24	!* Zeleznik valid only up to 350 K
25	!*
26	!* Input/output variables:
27	!* REAL(KIND=4) RADIUS,TAIR,PPWV,WSAS,RHOSA,MSA
28	!*
29	!* Input:
30	!* RADIUS: m Radius of aerosol droplet
31	!* TAIR: K Temperature of ambient air
32	!* PPWV: Pa Partial pressure of ambient water vapor
33	!*
34	!* Output:
35	!* WSAS: mass fraction of sulfuric acid. [0.1;1]
36	!* RHOSA: kg/m**3 Density of sulfuric acid solution droplet
37	!* MSAD: kg Mass of sulfuric acid in droplet
38	!* modified from
39	!* PROGRAM PSC_MODEL_E
40	!* by A. Määttänen & Slimane Bekki
41	!* subroutine for LMDZ+photochemistry VENUS
42	!* by A. Stolzenbach
43	!*
44	!* Modified by S.Guilbon for microphysical module to Venus GCM
45	!*
46
47	USE donnees
48	USE free_param
49
50	IMPLICIT NONE
51
52	! Inputs
53	REAL, intent(in) :: RADIUS, TAIR, PAIR
54	! Outputs
55	REAL, intent(out) :: WSAS, MSAD
56	! Auxilary variables:
57	REAL :: mrt_wv, mrt_sa
58	REAL :: N_H2SO4, N_H2O
59	REAL :: H2SO4_liq, H2O_liq
60	REAL :: CONCM
61	REAL :: MCONDTOT
62	REAL :: RMODE
63	REAL :: WSAFLAG
64	REAL :: power
65	! Ridder's Method variables:
66	REAL :: WVMIN, WVMAX, WVACC
67
68	INTEGER :: NBROOT
69
70	INTEGER :: MAXITE
71	PARAMETER(MAXITE=20)
72
73	INTEGER :: NBRAC
74	PARAMETER(NBRAC=5)
75
76	INTEGER :: FLAG
77
78	! External functions needed:
79	REAL :: IRFRMWV
80
81	! Physical constants:
82	REAL :: MH2O
83
84	! External functions needed:
85	REAL :: PWVSAS_GV,STSAS,ROSAS
86	! PWVSAS_GV: Natural logaritm of water vapor pressure over
87	! sulfuric acid solution
88	! STSAS: Surface tension of sulfuric acid solution
89	! ROSAS: Density of sulfuric acid solution
90
91	! Auxiliary local variables:
92	REAL :: DELW,DELLP,C1,C2,W0,W1,W2,F0,F1,WGUESS,LPPWV,RO
93	REAL :: psatwv,watact
94	INTEGER :: ITERAT
95	! write(,)'WSA ROSA NEW', RADIUS
96	MH2O=MWV
97
98	C1=2.0D0*MH2O/RGAS
99	C2=4.0D0*PI/3.0D0
100
101	mrt_sa=ppsa/pair
102	mrt_wv=ppwv/pair
103
104	! Initialisation des bornes pour WV
105	WVMIN=1.D-35
106	WVMAX=mrt_wv
107
108	! Accuracy de WVeq
109	WVACC=WVMAX*1.0D-3
110
111	! BRACWV borne la fonction f(WV) - WV = 0
112	! de WV=0 WV=WVtot on cherche l'intervalle o f(WV) - WV = 0
113	! avec prcisment f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
114	! Elle fait appel la fct/ssrtine ITERWV()
115	CALL BRACWV(TAIR,PAIR,WVMIN,WVMAX,NBRAC,RADIUS,mrt_wv,mrt_sa,FLAG,WSAFLAG,NBROOT)
116
117	SELECT CASE(FLAG)
118
119	CASE(1)
120	! Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
121	! IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
122	! Elle fait appel la fct/ssrtine ITERWV()
123
124	WSAS=IRFRMWV(TAIR,PAIR,WVMIN,WVMAX,WVACC,MAXITE,RADIUS,mrt_wv,mrt_sa,NBROOT)
125	RHOSA = ROSAS(TAIR,WSAS)
126	MSAD = C2WSASRHOSARADIUS*3
127
128	CASE(2)
129	! Cas NROOT=0 mais proche de 0
130	WSAS=WSAFLAG
131	RHOSA=ROSAS(TAIR,WSAS)
132	MSAD=C2WSASRHOSARADIUS*3
133	! ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
134	! ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
135	! incohrentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
136	! Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
137	! donne par ROSAS (cf test externe en Python), sinon, la valeur est trop
138	! basse (de l'ordre de 1000 kg/m3) et correspond parfois la valeur avec
139	! WSA=0.1 (pas totalement sr)
140	! En tous cas, incohrent avec ce qui est attendue pour le WSA et T donn
141	! La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
142	! la bonne valeur (tests externes Python confirment)
143
144	IF ((RHOSA.LT.1100.0D0).AND. (WSAS.GT.0.1D0))THEN
145	PRINT*,'PROBLEM RHO_DROPLET'
146	PRINT*,'rho_droplet',RHOSA
147	PRINT*,'T',TAIR,'WSA',WSAS
148	PRINT*,'ROSAS',ROSAS(TAIR, WSAS)
149	PRINT*,'FLAG',FLAG,'NROOT',NBROOT
150	STOP
151	ENDIF
152
153	CASE(3)
154	write(,)'Case 0 NROOT'
155	RHOSA=0.0D+0
156	WSAS=0.0D+0
157	MSAD=0.0D+0
158
159	END SELECT
160
161
162
163	RETURN
164
165	END SUBROUTINE WSA_ROSA_NEW
166
167
168	!*****************************************************************************
169	SUBROUTINE ITERWV(TAIR,PAIR,WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,RADIUS)
170	!* Cette routine est la solution par itration afin de trouver WSA pour un WV,
171	!* et donc LPPWV, donn. Ce qui nous donne egalement le WV correspondant au
172	!* WSA solution
173	!* For VenusGCM by A. Stolzenbach 07/2014
174	!* OUTPUT: WVEQ et WSAOUT
175
176	USE donnees
177	USE free_param
178	IMPLICIT NONE
179
180	REAL, INTENT(IN) :: TAIR, PAIR
181	REAL, INTENT(IN) :: WV, WVTOT, SATOT, RADIUS
182	REAL, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ
183
184	REAL :: LPPWV
185
186	REAL :: WSAMIN, WSAMAX, WSAACC
187	PARAMETER(WSAACC=0.01D0)
188
189	INTEGER :: MAXITSA, NBRACSA, NBROOT
190	PARAMETER(MAXITSA=20)
191	PARAMETER(NBRACSA=5)
192
193	LOGICAl :: FLAG1,FLAG2
194
195	! External Function
196	REAL :: IRFRMSA, WVCOND
197
198	IF (RADIUS.LT.1D-30) THEN
199	PRINT*,'RMODE == 0 FLAG 3', RADIUS
200	STOP
201	ENDIF
202
203	! Initialisation WSA=[0.1,1.0]
204	WSAMIN = 0.1D0
205	WSAMAX = 1.0D0
206	LPPWV = DLOG(PAIR*WV)
207
208	! Appel Bracket de KEEQ
209	CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR,LPPWV,FLAG1,FLAG2,NBROOT)
210
211	IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN
212
213	WSAOUT=IRFRMSA(TAIR,PAIR,WSAMIN,WSAMAX,WSAACC,MAXITSA,RADIUS,LPPWV,NBROOT)
214
215	!!$ ! AM uncommented the two following lines to avoid problems with nucleation
216	!!$ IF (WSAOUT.GT.1.0) WSAOUT=0.999999
217	!!$ IF (WSAOUT.LT.0.1) WSAOUT=0.1
218	!!$ write(,) 'in 1 wsaout 2', WSAOUT
219
220	! Si BRACWSA ne trouve aucun ensemble solution KEEQ=0 on fixe WSA a 0.9999 ou 0.1
221	ELSE
222	IF (FLAG1.AND.(.NOT.FLAG2)) WSAOUT = 0.999999D0
223	IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT = WSAMIN
224	IF (FLAG1.AND.FLAG2) THEN
225	PRINT*,'FLAGs BARCWSA tous TRUE'
226	STOP
227	ENDIF
228	ENDIF
229
230	! WVEQ output correspondant a WVliq lie a WSA calcule
231	WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT)
232	WVEQOUT=(WVLIQ+WV)/WVTOT-1.0D0
233
234	END SUBROUTINE ITERWV
235
236
237	!*****************************************************************************
238	SUBROUTINE BRACWV(TAIR,PAIR,XA,XB,N,RADIUS,WVTOT,SATOT,FLAGWV,WSAFLAG,NROOT)
239
240	!* Bracket de ITERWV
241	!* From Numerical Recipes
242	!* Adapted for VenusGCM A. Stolzenbach 07/2014
243	!* X est WVinput
244	!* OUTPUT: XA et XB
245
246	USE donnees
247	USE free_param
248
249	IMPLICIT NONE
250
251	REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR, PAIR
252	INTEGER, INTENT(IN) :: N
253
254	REAL, INTENT(INOUT) :: XA,XB
255	REAL, INTENT(OUT) :: WSAFLAG
256
257	INTEGER :: I,J
258
259	INTEGER, INTENT(OUT) :: NROOT
260
261	REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT
262	REAL :: XMAX,XMIN,WVEQACC
263
264	INTEGER, INTENT(OUT) :: FLAGWV
265	! write(,)'BRACWV', RADIUS
266	! WVEQACC est le seuil auquel on accorde un WSA correct meme
267	! si il ne fait pas partie d'une borne. Utile quand le modele
268	! s'approche de 0 mais ne l'atteint pas.
269	WVEQACC = 1.0D-3
270	FLAGWV = 1
271	NROOT = 0
272
273	! 25/11/2016
274	! On change ordre on va du max au min
275	X = XB
276	XMAX = XB
277	XMIN = XA
278
279	! CAS 1 On borne la fonction (WVEQ=0)
280	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,RADIUS)
281
282	FP=WVEQ
283
284	DO I=N-1,1,-1
285	X=(1.-DLOG(DBLE(N-I))/DLOG(DBLE(N)))*XMAX
286
287	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,RADIUS)
288
289	FC=WVEQ
290
291	IF ((FP*FC).LT.0.D0) THEN
292	NROOT=NROOT+1
293	! Si NROOT>1 on place la borne sup output la borne min du calcul en i
294	IF (NROOT.GT.1) THEN
295	XB=(1.-DLOG(DBLE(I+1))/DLOG(DBLE(N)))*XMAX
296	ENDIF
297
298	IF (I.EQ.1) THEN
299	XA=XMIN
300	ELSE
301	XA=X
302	ENDIF
303	RETURN
304	ENDIF
305	FP=FC
306	ENDDO
307
308	! CAS 2 on refait la boucle pour tester si WVEQ est proche de 0
309	! avec le seuil WVEQACC
310	IF (NROOT.EQ.0) THEN
311	X=XMAX
312
313	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,RADIUS)
314
315	DO J=N-1,1,-1
316	X=(1.-DLOG(DBLE(N-J))/DLOG(DBLE(N)))*XMAX
317	! write(,) 'BRACWV, bf 4th ITERWV (cas 2) '
318	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,RADIUS)
319
320	IF (ABS(WVEQ).LE.WVEQACC) THEN
321	WSAFLAG=WSAOUT
322	FLAGWV=2
323	RETURN
324	ENDIF
325	ENDDO
326
327	! CAS 3 Pas de borne, WVEQ jamais proche de 0
328	FLAGWV=3
329	RETURN
330	ENDIF
331
332	END SUBROUTINE BRACWV
333
334
335	!*****************************************************************************
336	SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL,NROOT)
337
338	!* Bracket de KEEQ
339	!* From Numerical Recipes
340	!* Adapted for Venus GCM A. Stolzenbach 07/2014
341
342	USE donnees
343	USE free_param
344	IMPLICIT NONE
345
346	! External functions needed:
347	REAL :: KEEQ
348
349	REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP
350	INTEGER, INTENT(IN) :: N
351
352	REAL, INTENT(INOUT) :: XA,XB
353
354	INTEGER, INTENT(OUT) :: NROOT
355
356	INTEGER :: I, J
357
358	REAL :: DX, FP, FC, X
359
360	LOGICAL, INTENT(OUT) :: FLAGH,FLAGL
361
362	FLAGL=.FALSE.
363	FLAGH=.FALSE.
364	NROOT=0
365	DX=(XB-XA)/N
366	X=XB
367	FP=KEEQ(TAIR,RADIUS,X,LPPWVINP)
368
369	DO I=N,1,-1
370	X=X-DX
371
372	FC=KEEQ(TAIR,RADIUS,X,LPPWVINP)
373
374	IF ((FP*FC).LE.0.) THEN
375	NROOT=NROOT+1
376	XA=X
377	XB=X+DX
378	RETURN
379	ENDIF
380
381	FP=FC
382	ENDDO
383
384	IF (NROOT.EQ.0) THEN
385	! Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]
386	IF ((ABS(KEEQ(TAIR,RADIUS,XA,LPPWVINP))- &
387	& ABS(KEEQ(TAIR,RADIUS,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE.
388	! On fixe flag low TRUE pour WSA = 0.1
389	IF ((ABS(KEEQ(TAIR,RADIUS,XA,LPPWVINP))- &
390	& ABS(KEEQ(TAIR,RADIUS,XB,LPPWVINP))).LT.0.0) FLAGL=.TRUE.
391	ENDIF
392
393	END SUBROUTINE BRACWSA
394
395
396	!*****************************************************************************
397	FUNCTION IRFRMWV(TAIR,PAIR,X1,X2,XACC,MAXIT,RADIUS,WVTOT,SATOT,NROOT)
398
399	!* Iterative Root Finder Ridder's Method for Water Vapor calculus
400	!* From Numerical Recipes
401	!* Adapted for VenusGCM A. Stolzenbach 07/2014
402
403	!* Les iterations sur [X1,X2] sont [WV1,WV2]
404	!* la variable X est WV
405	!* IRFRMWV sort en OUTPUT : WSALOC pour ITERWV=0 (ou WVEQ=0)
406
407	USE donnees
408	USE free_param
409	IMPLICIT NONE
410
411	REAL, INTENT(IN) :: TAIR, PAIR
412	REAL, INTENT(IN) :: X1, X2
413	REAL, INTENT(IN) :: XACC
414	REAL :: IRFRMWV
415	INTEGER, INTENT(IN) :: MAXIT,NROOT
416
417	! LOCAL VARIABLES
418	REAL :: XL, XH, XM, XNEW, X
419	REAL :: WSALOC, WVEQ, WVLIQ
420	REAL :: FL, FH, FM, FNEW
421	REAL :: ANS, S, FSIGN
422	INTEGER i
423	LOGICAL :: FLAGH,FLAGL
424
425	! External variables needed:
426	REAL, INTENT(IN) :: WVTOT,SATOT
427	REAL, INTENT(IN) :: RADIUS
428
429	! Initialisation
430	X=X1
431	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
432	FL=WVEQ
433	X=X2
434	CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
435	FH=WVEQ
436
437	! Test Bracketed values
438	IF (((FL.LT.0.).AND.(FH.GT.0.)).OR. &
439	& ((FL.GT.0.).AND.(FH.LT.0.))) &
440	& THEN
441	XL=X1
442	XH=X2
443	ANS=-1.D38
444
445	DO i=1, MAXIT
446	XM=0.5D0*(XL+XH)
447	CALL ITERWV(TAIR,PAIR,XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
448	FM=WVEQ
449	S=SQRT(FMFM-FLFH)
450
451	IF (S.EQ.0.0) THEN
452	IRFRMWV=WSALOC
453	RETURN
454	ENDIF
455
456	IF (FL.GT.FH) THEN
457	FSIGN=1.0D0
458	ELSE
459	FSIGN=-1.0D0
460	ENDIF
461
462	XNEW=XM+(XM-XL)(FSIGNFM/S)
463
464	IF (ABS(XNEW-ANS).LE.XACC) THEN
465	IRFRMWV=WSALOC
466	RETURN
467	ENDIF
468
469	ANS=XNEW
470	CALL ITERWV(TAIR,PAIR,ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
471	FNEW=WVEQ
472
473	IF (FNEW.EQ.0.D0) THEN
474	IRFRMWV=WSALOC
475	RETURN
476	ENDIF
477
478	IF (SIGN(FM, FNEW).NE.FM) THEN
479	XL=XM
480	FL=FM
481	XH=ANS
482	FH=FNEW
483	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
484	XH=ANS
485	FH=FNEW
486	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
487	XL=ANS
488	FL=FNEW
489	ELSE
490	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
491	PRINT*,'you shall not PAAAAAASS'
492	STOP
493	ENDIF
494	ENDDO
495	PRINT*,'Paaaaas bien MAXIT atteint'
496	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
497	PRINT*,'you shall not PAAAAAASS'
498	XL=X1
499	XH=X2
500	! ANS=-9.99e99
501	ANS=-1.D38
502
503	DO i=1, MAXIT
504	XM=0.5D0*(XL+XH)
505	CALL ITERWV(TAIR,PAIR,XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
506	FM=WVEQ
507	S=SQRT(FMFM-FLFH)
508	IF (FL.GT.FH) THEN
509	FSIGN=1.0D0
510	ELSE
511	FSIGN=-1.0D0
512	ENDIF
513
514	XNEW=XM+(XM-XL)(FSIGNFM/S)
515
516	ANS=XNEW
517	CALL ITERWV(TAIR,PAIR,ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
518	FNEW=WVEQ
519	PRINT*,'WVliq',WVLIQ,'WVtot',WVTOT,'WVeq',WVEQ
520	PRINT*,'WSA',WSALOC,'SAtot',SATOT
521	PRINT*,'T',TAIR,'P',PAIR
522
523	IF (SIGN(FM, FNEW).NE.FM) THEN
524	XL=XM
525	FL=FM
526	XH=ANS
527	FH=FNEW
528	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
529	XH=ANS
530	FH=FNEW
531	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
532	XL=ANS
533	FL=FNEW
534	ELSE
535	PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
536	PRINT*,'you shall not PAAAAAASS TWIIICE???'
537	STOP
538	ENDIF
539	ENDDO
540	STOP
541	ELSE
542	PRINT*,'IRFRMWV must be bracketed'
543	PRINT*,'NROOT de BRACWV', NROOT
544	IF (ABS(FL).LT.XACC) THEN
545	PRINT*,'IRFRMWV FL == 0',FL
546	PRINT*,'X1',X1,'X2',X2,'FH',FH
547	CALL ITERWV(TAIR,PAIR,X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
548	IRFRMWV=WSALOC
549	RETURN
550	ENDIF
551	IF (ABS(FH).LT.XACC) THEN
552	PRINT*,'IRFRMWV FH == 0',FH
553	PRINT*,'X1',X1,'X2',X2,'FL',FL
554	CALL ITERWV(TAIR,PAIR,X2,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,RADIUS)
555	IRFRMWV=WSALOC
556	RETURN
557	ENDIF
558	IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN
559	PRINT*,'STOP dans IRFRMWV avec rien == 0'
560	PRINT*,'X1',X1,'X2',X2
561	PRINT*,'Fcalc',FL,FH
562	PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS
563	STOP
564	ENDIF
565	IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN
566	PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC'
567	PRINT*,FL,FH
568	STOP
569	ENDIF
570
571	end IF
572	END FUNCTION IRFRMWV
573
574	!*****************************************************************************
575	FUNCTION IRFRMSA(TAIR,PAIR,X1,X2,XACC,MAXIT,RADIUS,LPPWV,NB)
576
577	!* Iterative Root Finder Ridder's Method for Sulfuric Acid calculus
578	!* From Numerical Recipes
579	!* Adapted for VenusGCM A. Stolzenbach 07/2014
580	!*
581	!* Les iterations sur [X1,X2] sont [WSA1,WSA2]
582	!* la variable X est WSA
583	!* IRFRMSA sort en OUTPUT : WSA pour KEEQ=0
584
585	use donnees
586	use free_param
587	IMPLICIT NONE
588
589	REAL, INTENT(IN) :: TAIR, PAIR
590	REAL, INTENT(IN) :: X1, X2
591	REAL, INTENT(IN) :: XACC
592	INTEGER, INTENT(IN) :: MAXIT, NB
593
594	! LOCAL VARIABLES
595	REAL :: IRFRMSA
596	REAL :: XL, XH, XM, XNEW
597	REAL :: Fl, FH, FM, FNEW
598	REAL :: ANS, S, FSIGN
599	INTEGER i
600
601	! External variables needed:
602	REAL, INTENT(IN) :: LPPWV
603	REAL, INTENT(IN) :: RADIUS
604
605	! External functions needed:
606	REAL :: KEEQ
607
608	! Initialisation
609	FL=KEEQ(TAIR,RADIUS,X1,LPPWV)
610	FH=KEEQ(TAIR,RADIUS,X2,LPPWV)
611
612	! Test Bracketed values
613	IF (((FL.LT.0.D0).AND.(FH.GT.0.D0)).OR.((FL.GT.0.D0).AND.(FH.LT.0.D0))) THEN
614
615	XL=X1
616	XH=X2
617
618	ANS=-1.D38
619
620	DO i=1, MAXIT
621	XM=0.5D0*(XL+XH)
622	FM=KEEQ(TAIR,RADIUS,XM,LPPWV)
623
624	S=SQRT(FMFM-FLFH)
625
626	IF (S.EQ.0.D0) THEN
627	IRFRMSA=ANS
628	RETURN
629	ENDIF
630
631	IF (FL.GT.FH) THEN
632	FSIGN=1.0D0
633	ELSE
634	FSIGN=-1.0D0
635	ENDIF
636
637	XNEW=XM+(XM-XL)(FSIGNFM/S)
638
639	IF (ABS(XNEW-ANS).LE.XACC) THEN
640	IRFRMSA=ANS
641
642	RETURN
643	ENDIF
644
645	ANS=XNEW
646	FNEW=KEEQ(TAIR,RADIUS,ANS,LPPWV)
647
648	IF (FNEW.EQ.0.D0) THEN
649	IRFRMSA=ANS
650	RETURN
651	ENDIF
652
653	IF (SIGN(FM, FNEW).NE.FM) THEN
654	XL=XM
655	FL=FM
656	XH=ANS
657	FH=FNEW
658	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
659	XH=ANS
660	FH=FNEW
661	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
662	XL=ANS
663	FL=FNEW
664	ELSE
665	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
666	PRINT*,'you shall not PAAAAAASS'
667	STOP
668	ENDIF
669
670	ENDDO
671	PRINT*,'Paaaaas bien MAXIT atteint'
672	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
673	PRINT*,'you shall not PAAAAAASS'
674	XL=X1
675	XH=X2
676	PRINT*,'Borne XL',XL,'XH',XH
677
678	ANS=-1.D38
679
680	DO i=1, MAXIT
681	XM=0.5D0*(XL+XH)
682	FM=KEEQ(TAIR,RADIUS,XM,LPPWV)
683	S=SQRT(FMFM-FLFH)
684
685	IF (FL.GT.FH) THEN
686	FSIGN=1.0D0
687	ELSE
688	FSIGN=-1.0D0
689	ENDIF
690
691	XNEW=XM+(XM-XL)(FSIGNFM/S)
692	ANS=XNEW
693	FNEW=KEEQ(TAIR,RADIUS,ANS,LPPWV)
694	PRINT*,'KEEQ result',FNEW,'T',TAIR,'R',RADIUS
695	IF (SIGN(FM, FNEW).NE.FM) THEN
696	XL=XM
697	FL=FM
698	XH=ANS
699	FH=FNEW
700	ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
701	XH=ANS
702	FH=FNEW
703	ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
704	XL=ANS
705	FL=FNEW
706	ELSE
707	PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
708	PRINT*,'you shall not PAAAAAASS'
709	STOP
710	ENDIF
711	ENDDO
712	STOP
713
714	ELSE
715	PRINT*,'IRFRMSA must be bracketed'
716	IF (FL.EQ.0.D0) THEN
717	PRINT*,'IRFRMSA FL == 0',Fl
718	IRFRMSA=X1
719	RETURN
720	ENDIF
721	IF (FH.EQ.0.D0) THEN
722	PRINT*,'IRFRMSA FH == 0',FH
723	IRFRMSA=X2
724	RETURN
725	ENDIF
726	IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN
727	PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH
728	PRINT*,'X1',X1,'X2',X2
729	PRINT*,'Kind F', KIND(FL), KIND(FH)
730	PRINT*,'Kind X', KIND(X1), KIND(X2)
731	PRINT*,'Logical: ',(SIGN(FL,FH).NE.FL)
732	PRINT*,'Logical: ',(SIGN(FH,FL).NE.FH)
733	PRINT*,'nb root BRACWSA',NB
734	STOP
735	ENDIF
736
737	ENDIF
738
739	END function IRFRMSA
740
741	!*****************************************************************************
742	FUNCTION KEEQ(TAIR,RADIUS,WX,LPPWV)
743
744	!* Kelvin Equation EQuality
745	!* ln(PPWV_eq) - (2Mh2o sigma)/(R T r rho) - ln(ph2osa) = 0
746
747	use donnees
748	use free_param
749	IMPLICIT NONE
750
751	REAL, INTENT(IN) :: RADIUS,WX,LPPWV,TAIR
752
753	! Physical constants:
754	REAL :: KEEQ
755
756	! External functions needed:
757	REAL :: PWVSAS_GV, STSAS, ROSAS
758	! PWVSAS_GV: Natural logaritm of water vapor pressure over
759	! sulfuric acid solution
760	! STSAS: Surface tension of sulfuric acid solution
761	! ROSAS: Density of sulfuric acid solution
762	!
763	! Auxiliary local variables:
764	REAL :: C1
765
766	C1=2.0D0*MWV/RGAS
767
768	KEEQ=LPPWV-C1STSAS(TAIR,WX)/(TAIRRADIUS*ROSAS(TAIR,WX))- &
769	& PWVSAS_GV(TAIR,WX)
770
771	END FUNCTION KEEQ
772
773
774	!*****************************************************************************
775	FUNCTION WVCOND(WX,T,P,SAt)
776
777	!* Condensation de H2O selon WSA, T et P et H2SO4tot
778
779	!* Adapted for VenusGCM A. Stolzenbach 07/2014
780	! INPUT:
781	! SAt : VMR of total H2SO4
782	! WX: aerosol H2SO4 weight fraction (fraction)
783	! T: temperature (K)
784	! P: pressure (Pa)
785	! OUTPUT:
786	! WVCOND : VMR H2O condense
787
788	! USE chemparam_mod
789
790	use donnees
791	use free_param
792
793	IMPLICIT NONE
794
795	REAL, INTENT(IN) :: SAt, WX
796	REAL, INTENT(IN) :: T, P
797
798	! working variables
799	REAL :: WVCOND
800	REAL :: SA, WV, KBOLTZ
801	REAL :: DND2,pstand,lpar,acidps
802	REAL :: x1, satpacid
803	REAL, DIMENSION(2):: act
804	REAL :: CONCM
805	REAL :: NH2SO4
806	REAL :: H2OCOND, H2SO4COND
807
808	KBOLTZ=KBZ
809	CONCM= (P)/(KBOLTZ*T) !air number density, molec/m3? CHECK UNITS!
810
811	NH2SO4=SAt*CONCM
812	pstand=1.01325D+5 !Pa 1 atm pressure
813
814	x1=(WX/MSA)/(WX/MSA + ((1.-WX)/MWV))
815
816	CALL zeleznik(x1,T,act)
817
818	!pure acid satur vapor pressure
819	lpar= -11.695D0 + DLOG(pstand) ! Zeleznik
820	acidps = 1/360.15D0 - 1.0/T + 0.38D0/545.D0*(1.0+DLOG(360.15D0/T)-360.15D0/T)
821	acidps = 10156.D0*acidps + lpar
822	acidps = DEXP(acidps) !Pa
823
824	!acid sat.vap.PP over mixture (flat surface):
825	satpacid=act(2)*acidps ! Pa
826
827	! Conversion from Pa to N.D #/m3
828	DND2=satpacid/(KBOLTZ*T)
829	! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
830	IF (NH2SO4.GT.DND2) THEN
831	H2SO4COND=NH2SO4-DND2
832	! calcul de H2O cond correspondant a H2SO4 cond
833	H2OCOND=H2SO4CONDMSA(1.0D0-WX)/(MWV*WX)
834	! Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30
835	ELSE
836	H2OCOND=1.0D-30*CONCM
837	END IF
838
839	!*****************************************************
840	! ATTENTION: Ici on ne prends pas en compte
841	! si H2O en defaut!
842	! On veut la situation thorique
843	! l'equilibre
844	!*****************************************************
845	! Test si H2O en defaut H2Ocond>H2O dispo
846	! IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN
847	! On peut alors condenser tout le H2O dispo
848	! H2OCOND=NH2O
849	! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
850	! H2SO4COND=H2OCOND18.0153WSA/(98.078*(1.0-WSA))
851	! END IF
852
853	! Calcul de H2O condense VMR
854	WVCOND=H2OCOND/CONCM
855
856	END FUNCTION WVCOND
857
858
859	!*****************************************************************************
860	FUNCTION PWVSAS_GV(T,W)
861
862	!* Natural logaritm of saturated water vapor pressure over plane
863	!* sulfuric acid solution.
864	!*
865	!* Source: J.I.Gmitro & T.Vermeulen: A.I.Ch.E.J. 10,740,1964.
866	!* W.F.Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
867	!*
868	!* The formula of Gmitro & Vermeulen for saturation pressure
869	!* is used:
870	!* ln(p) = A ln(298/T) + B/T + C + DT
871	!* with values of A,B,C and D given by Gmitro & Vermeulen,
872	!* and calculated from partial molal properties given by Giauque et al.
873	!*
874	!* Input: T: Temperature (K)
875	!* W: Weight fraction of H2SO4 [0;1]
876	!* Output: Natural logaritm of water vapor pressure
877	!* over sulfuric acid solution ( ln(Pa) )
878	!*
879	!* External functions needed for calculation of partial molal
880	!* properties of pure components at 25 C as function of W.
881
882	use donnees
883	IMPLICIT NONE
884
885	REAL :: CPH2O,ALH2O,FFH2O,LH2O
886	! CPH2O: Partial molal heat capacity of sulfuric acid solution.
887	! ALH2O: Temparature derivative of CPH2O
888	! FFH2O: Partial molal free energy of sulfuric acid solution.
889	! LH2O: Partial molal enthalpy of sulfuric acid
890	!
891	!
892	!
893	REAL, INTENT(IN) :: T,W
894	REAL :: PWVSAS_GV
895	REAL :: ADOT,BDOT,CDOT,DDOT
896	REAL :: RGAScal,MMHGPA
897	REAL :: K1,K2
898	REAL :: A,B,C,Dd,CP,L,F,ALFA
899	! Physical constants given by Gmitro & Vermeulen:
900	PARAMETER( &
901	ADOT=-3.67340, &
902	BDOT=-4143.5, &
903	CDOT=10.24353, &
904	DDOT=0.618943d-3)
905	PARAMETER( &
906	! Gas constant (cal/(deg mole)):
907	RGAScal=1.98726, &
908	! Natural logarith of conversion factor between atm. and Pa:
909	MMHGPA=11.52608845, &
910	K1=298.15, &
911	K2=K1*K1/2.0)
912	!
913	INTEGER :: KLO,KHI,K,I,NPOINT
914	PARAMETER(NPOINT=110)
915	REAL, DIMENSION(NPOINT) :: WTAB(NPOINT)
916	DATA (WTAB(I),I=1,NPOINT)/ &
917	0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, &
918	0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, &
919	0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, &
920	0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, &
921	0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, &
922	0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, &
923	0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, &
924	0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, &
925	0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, &
926	0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, &
927	0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, &
928	0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, &
929	0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, &
930	0.99908,0.99927,0.99945,0.99963,0.99982,1.0000/
931	!
932	KLO=1
933	KHI=NPOINT
934	1 IF(KHI-KLO.GT.1) THEN
935	K=(KHI+KLO)/2
936	IF(WTAB(K).GT.MAX(WTAB(1),W)) THEN
937	KHI=K
938	ELSE
939	KLO=K
940	ENDIF
941	GOTO 1
942	ENDIF
943	!
944	CP=CPH2O(W,KHI,KLO)
945	F=-FFH2O(W,KHI,KLO)
946	L=-LH2O(W,KHI,KLO)
947	ALFA=ALH2O(W,KHI,KLO)
948	!
949	A=ADOT+(CP-K1*ALFA)/RGAScal
950	B=BDOT+(L-K1CP+K2ALFA)/RGAScal
951	C=CDOT+(CP+(F-L)/K1)/RGAScal
952	Dd=DDOT-ALFA/(2.0d0*RGAScal)
953	!
954	! WRITE(,) 'TAIR= ',T,' WSA= ',W
955	! WRITE(,) 'CPH2O(W)= ',CP
956	! WRITE(,) 'ALFAH2O(W)= ',ALFA
957	! WRITE(,) 'FFH2O(W)= ',F
958	! WRITE(,) 'LH2O(W)= ',L
959	!
960	PWVSAS_GV=ADLOG(K1/T)+B/T+C+DdT+MMHGPA
961
962	END FUNCTION PWVSAS_GV
963
964
965	!*****************************************************************************
966	REAL FUNCTION CPH2O(W,khi_in,klo_in)
967
968	! Relative partial molal heat capacity of water (cal/(deg mole) in
969	! sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
970	! calculated by cubic spline fitting.
971	!
972	! Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
973
974	IMPLICIT NONE
975
976	INTEGER :: NPOINT,I
977	PARAMETER(NPOINT=109)
978	REAL, DIMENSION(NPOINT) :: WTAB(NPOINT),CPHTAB(NPOINT), &
979	Y2(NPOINT),YWORK(NPOINT)
980	REAL, INTENT(IN):: W
981	INTEGER, INTENT(IN):: khi_in,klo_in
982	INTEGER :: khi,klo
983	REAL :: CPH
984	LOGICAL :: FIRST
985	DATA (WTAB(I),I=1,NPOINT)/ &
986	0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, &
987	0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, &
988	0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, &
989	0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, &
990	0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, &
991	0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, &
992	0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, &
993	0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, &
994	0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, &
995	0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, &
996	0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, &
997	0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, &
998	0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, &
999	0.99908,0.99927,0.99945,0.99963,0.99982/
1000	DATA (CPHTAB(I),I=1,NPOINT)/ &
1001	17.996, 17.896, 17.875, 17.858, 17.840, 17.820, 17.800, 17.791, &
1002	17.783, 17.777, 17.771, 17.769, 17.806, 17.891, 18.057, 18.248, &
1003	18.429, 18.567, 18.613, 18.640, 18.660, 18.660, 18.642, 18.592, &
1004	18.544, 18.468, 18.348, 18.187, 17.995, 17.782, 17.562, 17.352, &
1005	17.162, 16.993, 16.829, 16.657, 16.581, 16.497, 16.405, 16.302, &
1006	16.186, 16.053, 15.901, 15.730, 15.540, 15.329, 15.101, 14.853, &
1007	14.586, 14.296, 13.980, 13.638, 13.274, 12.896, 12.507, 12.111, &
1008	11.911, 11.711, 11.514, 11.320, 11.130, 10.940, 10.760, 10.570, &
1009	10.390, 10.200, 10.000, 9.8400, 9.7600, 9.7900, 9.9500, 10.310, &
1010	10.950, 11.960, 13.370, 15.060, 16.860, 18.550, 20.000, 21.170, &
1011	22.030, 22.570, 22.800, 22.750, 22.420, 21.850, 21.120, 20.280, &
1012	19.360, 18.350, 17.220, 15.940, 14.490, 12.840, 10.800, 9.8000, &
1013	7.8000, 3.8000,0.20000,-5.4000,-7.0000,-8.8000,-10.900,-13.500, &
1014	-17.000,-22.000,-29.000,-40.000,-59.000/
1015	DATA FIRST/.TRUE./
1016	SAVE FIRST,WTAB,CPHTAB,Y2
1017	!
1018	IF(FIRST) THEN
1019	FIRST=.FALSE.
1020	CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2)
1021	ENDIF
1022
1023	if(khi_in.GT.NPOINT) then
1024	khi=NPOINT
1025	klo=NPOINT-1
1026	else
1027	khi=khi_in
1028	klo=klo_in
1029	endif
1030
1031	CALL SPLINT(WTAB(khi),WTAB(klo),CPHTAB(khi),CPHTAB(klo),Y2(khi),Y2(klo),W,CPH)
1032	CPH2O=CPH
1033
1034	END FUNCTION CPH2O
1035
1036
1037	!*******************************************************************************
1038	REAL FUNCTION FFH2O(W,khi,klo)
1039
1040	! Relative partial molal free energy water (cal/mole) in
1041	! sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
1042	! calculated by cubic spline fitting.
1043
1044	! Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1045
1046	IMPLICIT NONE
1047
1048	INTEGER :: NPOINT,I
1049	PARAMETER(NPOINT=110)
1050	REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
1051	REAL, INTENT(IN) :: W
1052	INTEGER, INTENT(IN):: khi,klo
1053	REAL :: FF
1054	LOGICAL :: FIRST
1055	DATA (WTAB(I),I=1,NPOINT)/ &
1056	0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, &
1057	0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, &
1058	0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, &
1059	0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, &
1060	0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, &
1061	0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, &
1062	0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, &
1063	0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, &
1064	0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, &
1065	0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, &
1066	0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, &
1067	0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, &
1068	0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, &
1069	0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1070	DATA (FFTAB(I),I=1,NPOINT)/ &
1071	0.00000, 22.840, 25.810, 29.250, 33.790, 39.970, 48.690, 54.560, &
1072	61.990, 71.790, 85.040, 103.70, 130.70, 145.20, 163.00, 184.50, &
1073	211.50, 245.60, 266.40, 290.10, 317.40, 349.00, 385.60, 428.40, &
1074	452.50, 478.80, 507.50, 538.80, 573.30, 611.60, 653.70, 700.50, &
1075	752.60, 810.60, 875.60, 948.60, 980.60, 1014.3, 1049.7, 1087.1, &
1076	1126.7, 1168.7, 1213.5, 1261.2, 1312.0, 1366.2, 1424.3, 1486.0, &
1077	1551.8, 1622.3, 1697.8, 1778.5, 1864.9, 1956.8, 2055.8, 2162.0, &
1078	2218.0, 2276.0, 2337.0, 2400.0, 2466.0, 2535.0, 2607.0, 2682.0, &
1079	2760.0, 2842.0, 2928.0, 3018.0, 3111.0, 3209.0, 3311.0, 3417.0, &
1080	3527.0, 3640.0, 3757.0, 3878.0, 4002.0, 4130.0, 4262.0, 4397.0, &
1081	4535.0, 4678.0, 4824.0, 4973.0, 5128.0, 5287.0, 5454.0, 5630.0, &
1082	5820.0, 6031.0, 6268.0, 6541.0, 6873.0, 7318.0, 8054.0, 8284.0, &
1083	8579.0, 8997.0, 9295.0, 9720.0, 9831.0, 9954.0, 10092., 10248., &
1084	10423., 10618., 10838., 11099., 11460., 12014./
1085	DATA FIRST/.TRUE./
1086	SAVE FIRST,WTAB,FFTAB,Y2
1087	!
1088	IF(FIRST) THEN
1089	FIRST=.FALSE.
1090	CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2)
1091	ENDIF
1092
1093	CALL SPLINT(WTAB(khi),WTAB(klo),FFTAB(khi),FFTAB(klo),Y2(khi),Y2(klo),W,FF)
1094	FFH2O=FF
1095
1096	END FUNCTION FFH2O
1097
1098
1099	!*******************************************************************************
1100	REAL FUNCTION LH2O(W,khi,klo)
1101
1102	! Relative partial molal heat content of water (cal/mole) in
1103	! sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
1104	! calculated by cubic spline fitting.
1105
1106	! Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1107
1108	IMPLICIT NONE
1109
1110	INTEGER :: NPOINT,I
1111	PARAMETER(NPOINT=110)
1112	REAL, DIMENSION(NPOINT) :: WTAB,LTAB,Y2,YWORK
1113	REAL, INTENT(IN) :: W
1114	INTEGER, INTENT(IN):: khi,klo
1115	REAL :: L
1116	LOGICAL :: FIRST
1117	DATA (WTAB(I),I=1,NPOINT)/ &
1118	0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, &
1119	0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, &
1120	0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, &
1121	0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, &
1122	0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, &
1123	0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, &
1124	0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, &
1125	0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, &
1126	0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, &
1127	0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, &
1128	0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, &
1129	0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, &
1130	0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, &
1131	0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1132	DATA (LTAB(I),I=1,NPOINT)/ &
1133	0.00000, 5.2900, 6.1000, 7.1800, 8.7800, 11.210, 15.290, 18.680, &
1134	23.700, 31.180, 42.500, 59.900, 89.200, 106.70, 128.60, 156.00, &
1135	190.40, 233.80, 260.10, 290.00, 324.00, 362.50, 406.50, 456.10, &
1136	483.20, 512.40, 543.60, 577.40, 613.80, 653.50, 696.70, 744.50, &
1137	797.20, 855.80, 921.70, 995.70, 1028.1, 1062.3, 1098.3, 1136.4, &
1138	1176.7, 1219.3, 1264.7, 1313.0, 1364.3, 1418.9, 1477.3, 1539.9, &
1139	1607.2, 1679.7, 1757.9, 1842.7, 1934.8, 2035.4, 2145.5, 2267.0, &
1140	2332.0, 2401.0, 2473.0, 2550.0, 2631.0, 2716.0, 2807.0, 2904.0, &
1141	3007.0, 3118.0, 3238.0, 3367.0, 3507.0, 3657.0, 3821.0, 3997.0, &
1142	4186.0, 4387.0, 4599.0, 4819.0, 5039.0, 5258.0, 5476.0, 5694.0, &
1143	5906.0, 6103.0, 6275.0, 6434.0, 6592.0, 6743.0, 6880.0, 7008.0, &
1144	7133.0, 7255.0, 7376.0, 7497.0, 7618.0, 7739.0, 7855.0, 7876.0, &
1145	7905.0, 7985.0, 8110.0, 8415.0, 8515.0, 8655.0, 8835.0, 9125.0, &
1146	9575.0, 10325., 11575., 13500., 15200., 16125./
1147	DATA FIRST/.TRUE./
1148	SAVE FIRST,WTAB,LTAB,Y2
1149	!
1150	IF(FIRST) THEN
1151	FIRST=.FALSE.
1152	CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2)
1153	ENDIF
1154
1155	CALL SPLINT(WTAB(khi),WTAB(klo),LTAB(khi),LTAB(klo),Y2(khi),Y2(klo),W,L)
1156	LH2O=L
1157
1158	END FUNCTION LH2O
1159
1160
1161	!*******************************************************************************
1162	REAL FUNCTION ALH2O(W,khi_in,klo_in)
1163
1164	! Relative partial molal temperature derivative of heat capacity (water)
1165	! in sulfuric acid solution, (cal/deg**2), calculated by
1166	! cubic spline fitting.
1167
1168	! Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
1169
1170	IMPLICIT NONE
1171
1172	INTEGER :: NPOINT,I
1173	PARAMETER(NPOINT=96)
1174	REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
1175	REAL, INTENT(IN) :: W
1176	INTEGER, INTENT(IN):: khi_in,klo_in
1177	INTEGER :: khi,klo
1178	REAL :: A
1179	LOGICAL :: FIRST
1180	DATA (WTAB(I),I=1,NPOINT)/ &
1181	0.29517,0.31209, &
1182	0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, &
1183	0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, &
1184	0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, &
1185	0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, &
1186	0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, &
1187	0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, &
1188	0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, &
1189	0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, &
1190	0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, &
1191	0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, &
1192	0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, &
1193	0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
1194	DATA (ATAB(I),I=1,NPOINT)/ &
1195	0.0190, 0.0182, 0.0180, 0.0177, 0.0174, 0.0169, 0.0167, 0.0164, &
1196	0.0172, 0.0212, 0.0239, 0.0264, 0.0276, 0.0273, 0.0259, 0.0238, &
1197	0.0213, 0.0190, 0.0170, 0.0155, 0.0143, 0.0133, 0.0129, 0.0124, &
1198	0.0120, 0.0114, 0.0106, 0.0097, 0.0084, 0.0067, 0.0047, 0.0024, &
1199	-0.0002,-0.0031,-0.0063,-0.0097,-0.0136,-0.0178,-0.0221,-0.0263, &
1200	-0.0303,-0.0340,-0.0352,-0.0360,-0.0362,-0.0356,-0.0343,-0.0321, &
1201	-0.0290,-0.0251,-0.0201,-0.0137,-0.0058, 0.0033, 0.0136, 0.0254, &
1202	0.0388, 0.0550, 0.0738, 0.0962, 0.1198, 0.1300, 0.1208, 0.0790, &
1203	0.0348, 0.0058,-0.0102,-0.0211,-0.0292,-0.0350,-0.0390,-0.0418, &
1204	-0.0432,-0.0436,-0.0429,-0.0411,-0.0384,-0.0346,-0.0292,-0.0220, &
1205	-0.0130,-0.0110,-0.0080,-0.0060,-0.0040,-0.0030,-0.0030,-0.0020, &
1206	-0.0020,-0.0020,-0.0020,-0.0010,-0.0010, 0.0000, 0.0000, 0.0000/
1207	DATA FIRST/.TRUE./
1208	SAVE FIRST,WTAB,ATAB,Y2
1209	!
1210	IF(FIRST) THEN
1211	FIRST=.FALSE.
1212	CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2)
1213	ENDIF
1214
1215	if(klo_in.LT.15) then
1216	khi=2
1217	klo=1
1218	else
1219	khi=khi_in-14
1220	klo=klo_in-14
1221	endif
1222
1223	CALL SPLINT(WTAB(khi),WTAB(klo),ATAB(khi),ATAB(klo),Y2(khi),Y2(klo),W,A)
1224	ALH2O=A
1225
1226	END FUNCTION ALH2O
1227
1228	!******************************************************************************
1229	SUBROUTINE SPLINE(X,Y,N,WORK,Y2)
1230
1231	! Routine to calculate 2.nd derivatives of tabulated function
1232	! Y(i)=Y(Xi), to be used for cubic spline calculation.
1233
1234	IMPLICIT NONE
1235
1236	INTEGER N,I
1237	REAL,intent(in) :: X(N),Y(N)
1238	REAL,intent(out) :: WORK(N),Y2(N)
1239	REAL :: SIG,P,QN,UN,YP1,YPN
1240
1241	YP1=(Y(2)-Y(1))/(X(2)-X(1))
1242	YPN=(Y(N)-Y(N-1))/(X(N)-X(N-1))
1243
1244	IF(YP1.GT.99.0D+30) THEN
1245	Y2(1)=0.0
1246	WORK(1)=0.0
1247	ELSE
1248	Y2(1)=-0.5D0
1249	WORK(1)=(3.0D0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
1250	ENDIF
1251
1252	DO I=2,N-1
1253	SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
1254	P=SIG*Y2(I-1)+2.0D0
1255	Y2(I)=(SIG-1.0D0)/P
1256	WORK(I)=(6.0D0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) &
1257	& /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*WORK(I-1))/P
1258	ENDDO
1259
1260	IF(YPN.GT.99.0D+30) THEN
1261	QN=0.0
1262	UN=0.0
1263	ELSE
1264	QN=0.5D0
1265	UN=(3.0D0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
1266	ENDIF
1267
1268	Y2(N)=(UN-QNWORK(N-1))/(QNY2(N-1)+1.0D0)
1269
1270	DO I=N-1,1,-1
1271	Y2(I)=Y2(I)*Y2(I+1)+WORK(I)
1272	ENDDO
1273
1274	RETURN
1275	END SUBROUTINE SPLINE
1276
1277
1278	!******************************************************************************
1279	SUBROUTINE SPLINT(XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo,X,Y)
1280
1281	! Cubic spline calculation
1282
1283	IMPLICIT NONE
1284
1285	REAL, INTENT(IN) :: XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo
1286	REAL, INTENT(IN) :: X
1287	REAL, INTENT(OUT) :: Y
1288	REAL :: H,A,B
1289	!
1290	H=XAhi-XAlo
1291	A=(XAhi-X)/H
1292	B=(X-XAlo)/H
1293	Y=AYAlo+BYAhi+((A*3-A)Y2Alo+(B*3-B)Y2Ahi)(H*2)/6.0d0
1294	!
1295
1296	END SUBROUTINE SPLINT

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

Download in other formats:

Original Format