source: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/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
RevLine 
[1661]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!----------------------------------------------------------------------------
13SUBROUTINE  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
[1664]95!  write(*,*)'WSA ROSA NEW', RADIUS
[1661]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 prŽcisŽment 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 = C2*WSAS*RHOSA*RADIUS**3
127
128  CASE(2)
129     ! Cas NROOT=0 mais proche de 0
130     WSAS=WSAFLAG     
131     RHOSA=ROSAS(TAIR,WSAS)
132     MSAD=C2*WSAS*RHOSA*RADIUS**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     !     incohŽrentes 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 sžr)
140     !     En tous cas, incohŽrent 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
165END SUBROUTINE WSA_ROSA_NEW
166
167
168!*****************************************************************************                             
169SUBROUTINE ITERWV(TAIR,PAIR,WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,RADIUS)
170  !* Cette routine est la solution par itŽration 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
[1667]206  LPPWV = DLOG(PAIR*WV)
[1661]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
234END SUBROUTINE ITERWV
235
236
237!*****************************************************************************                             
238SUBROUTINE 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
[1664]265!  write(*,*)'BRACWV', RADIUS
[1661]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
[1667]295           XB=(1.-DLOG(DBLE(I+1))/DLOG(DBLE(N)))*XMAX
[1661]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
[1667]316        X=(1.-DLOG(DBLE(N-J))/DLOG(DBLE(N)))*XMAX
[1661]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
332END SUBROUTINE BRACWV
333
334
335!*****************************************************************************
336SUBROUTINE 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
[1667]367  FP=KEEQ(TAIR,RADIUS,X,LPPWVINP)
[1661]368
369  DO I=N,1,-1
370     X=X-DX
371
[1667]372     FC=KEEQ(TAIR,RADIUS,X,LPPWVINP)
[1661]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
393END SUBROUTINE BRACWSA
394
395
396!*****************************************************************************
397FUNCTION 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(FM*FM-FL*FH)
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)*(FSIGN*FM/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(FM*FM-FL*FH)
508        IF (FL.GT.FH) THEN
509           FSIGN=1.0D0
510        ELSE
511           FSIGN=-1.0D0
512        ENDIF
513
514        XNEW=XM+(XM-XL)*(FSIGN*FM/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
572END FUNCTION IRFRMWV
573
574!*****************************************************************************                             
575FUNCTION 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(FM*FM-FL*FH)
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)*(FSIGN*FM/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(FM*FM-FL*FH)
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)*(FSIGN*FM/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
739END function IRFRMSA
740
741!*****************************************************************************
742FUNCTION 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:
[1667]754  REAL :: KEEQ
[1661]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
[1667]766  C1=2.0D0*MWV/RGAS
[1661]767
768  KEEQ=LPPWV-C1*STSAS(TAIR,WX)/(TAIR*RADIUS*ROSAS(TAIR,WX))- &
769       &     PWVSAS_GV(TAIR,WX)
770
771END FUNCTION KEEQ
772
773
774!*****************************************************************************
775FUNCTION 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
[1667]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)
[1661]821  acidps = 10156.D0*acidps + lpar
[1667]822  acidps = DEXP(acidps)    !Pa
[1661]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=H2SO4COND*MSA*(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 thŽorique
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=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
851  !      END IF
852
853  !     Calcul de H2O condensŽe VMR         
854  WVCOND=H2OCOND/CONCM
855
856END FUNCTION WVCOND
857
858
859!*****************************************************************************
860FUNCTION 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
[1674]896      REAL :: RGAScal,MMHGPA
[1661]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(                   &
[1674]906!     Gas constant (cal/(deg mole)):
907           RGAScal=1.98726,        &
[1661]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!
[1674]949      A=ADOT+(CP-K1*ALFA)/RGAScal
950      B=BDOT+(L-K1*CP+K2*ALFA)/RGAScal
951      C=CDOT+(CP+(F-L)/K1)/RGAScal
952      Dd=DDOT-ALFA/(2.0d0*RGAScal)
[1661]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=A*DLOG(K1/T)+B/T+C+Dd*T+MMHGPA
961     
962      END FUNCTION PWVSAS_GV
963
964
965!*****************************************************************************
966REAL 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!*******************************************************************************
1038REAL 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!*******************************************************************************
1100REAL 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!*******************************************************************************
1162REAL 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!******************************************************************************
1229SUBROUTINE 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-QN*WORK(N-1))/(QN*Y2(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
1275END 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=A*YAlo+B*YAhi+((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.