!****************************************************************************** !* venus_SAS_composition SUBROUTINE !* modified from !* PROGRAM PSC_MODEL_E !* by A. Määttänen !* subroutine for LMDZ+photochemistry VENUS !* by A. Stolzenbach !* !* Input/Output files: !* ------------------- !* !---------------------------------------------------------------------------- SUBROUTINE new_cloud_venus( + nblev, nblon, + TT,PP, + mrt_wv,mrt_sa, + mr_wv,mr_sa) USE chemparam_mod IMPLICIT NONE #include "YOMCST.h" INTEGER, INTENT(IN) :: nblon ! nombre de points horizontaux INTEGER, INTENT(IN) :: nblev ! nombre de couches verticales !---------------------------------------------------------------------------- ! Ambient air state variables: REAL, INTENT(IN), DIMENSION(nblon,nblev) :: mrt_wv,mrt_sa, + TT,PP REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: mr_wv,mr_sa !---------------------------------------------------------------------------- INTEGER :: ilon, ilev, imode !---------------------------------------------------------------------------- ! Thermodynamic functions: REAL :: RHODROPLET !---------------------------------------------------------------------------- ! Auxilary variables: REAL :: NH2SO4,NH2O REAL :: H2SO4_liq,H2O_liq REAL :: CONCM REAL :: MCONDTOT REAL :: RMODE REAL :: WSAFLAG REAL :: K_SAV !---------------------------------------------------------------------------- ! Ridder's Method variables: REAL :: WVMIN, WVMAX, WVACC INTEGER :: NBROOT INTEGER :: MAXITE PARAMETER(MAXITE=20) INTEGER :: NBRAC PARAMETER(NBRAC=20) INTEGER :: FLAG !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- ! External functions needed: REAL :: IRFRMWV !---------------------------------------------------------------------------- ! >>> Program starts here: !AM Venus ! These aerosols will then be given an equilibrium composition for the given size distribution ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002, ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631 !=========================================== ! Debut boucle sur niveau et lat,lon !=========================================== ! Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0 NBRTOT(:,:,:)=0.0E+0 WH2SO4(:,:)=0.0E+0 rho_droplet(:,:)=0.0E+0 DO ilev=cloudmin, cloudmax DO ilon=1, nblon ! Boucle sur les modes RMODE=0.0E+0 K_SAV = 0.0 DO imode=1, nbr_mode IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN ! RMODE est le rayon modal de la distribution en volume du mode le plus ! representatif pour la Mtot RMODE=R_MEDIAN(ilon,ilev,imode)* & EXP(2.*(DLOG(STDDEV(ilon,ilev,imode))**2.)) K_SAV=K_MASS(ilon,ilev,imode) ENDIF ENDDO ! FIN boucle imode ! Initialisation des bornes pour WV WVMIN=1.E-90 WVMAX=mrt_wv(ilon,ilev) ! Accuracy de WVeq WVACC=WVMAX*1.0E-3 ! BRACWV borne la fonction f(WV) - WV = 0 ! de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0 ! avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0 ! Elle fait appel à la fct/ssrtine ITERWV() CALL BRACWV(WVMIN,WVMAX,NBRAC,RMODE, & mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),TT(ilon,ilev), & PP(ilon,ilev),FLAG,WSAFLAG,NBROOT) SELECT CASE(FLAG) CASE(1) ! Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant) ! IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0 ! Elle fait appel ˆ la fct/ssrtine ITERWV() WH2SO4(ilon,ilev)=IRFRMWV(WVMIN,WVMAX,WVACC,MAXITE,RMODE, & TT(ilon,ilev),PP(ilon,ilev), & mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT) rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev), & TT(ilon,ilev)) ! IF (rho_droplet(ilon,ilev).LT.1100.) THEN ! PRINT*,'PROBLEM RHO_DROPLET' ! PRINT*,'rho_droplet',rho_droplet(ilon,ilev) ! PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev) ! PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev), ! & TT(ilon,ilev)) ! PRINT*,'FLAG',FLAG,'NROOT',NBROOT ! STOP ! ENDIF CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3 NH2SO4=mrt_sa(ilon,ilev)*CONCM NH2O=mrt_wv(ilon,ilev)*CONCM CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev), & rho_droplet(ilon,ilev),TT(ilon,ilev), & H2SO4_liq,H2O_liq,MCONDTOT) ! Boucle sur les modes DO imode=1, nbr_mode IF (K_MASS(ilon,ilev,imode).GT.0.) THEN NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)* & K_MASS(ilon,ilev,imode)*MCONDTOT* & EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/ & (R_MEDIAN(ilon,ilev,imode)**3.) ELSE NBRTOT(ilon,ilev,imode)=0.0E+0 ENDIF ENDDO ! Passage de #/m3 en VMR H2O_liq=H2O_liq/CONCM H2SO4_liq=H2SO4_liq/CONCM mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq ! Problemes quand on a condense tout, on peut obtenir des -1e-24 ! aprs la soustraction et conversion de ND ˆ VMR IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30 IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30 CASE(2) ! Cas NROOT=0 mais proche de 0 WH2SO4(ilon,ilev)=WSAFLAG rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev), & TT(ilon,ilev)) ! ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation ! ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont ! incohŽrentes avec TT et WH2SO4 (a priori lorsque NTOT=0) ! Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur ! donne par RHODROPLET (cf test externe en Python), sinon, la valeur est trop ! basse (de l'ordre de 1000 kg/m3) et correspond parfois ˆ la valeur avec ! WSA=0.1 (pas totalement sur) ! En tous cas, incoherent avec ce qui est attendue pour le WSA et T donneeŽ ! La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a ! la bonne valeur (tests externes Python confirment) IF ((rho_droplet(ilon,ilev).LT.1100.).AND. & (WH2SO4(ilon,ilev).GT.0.1))THEN PRINT*,'PROBLEM RHO_DROPLET' PRINT*,'rho_droplet',rho_droplet(ilon,ilev) PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev) PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev), & TT(ilon,ilev)) PRINT*,'FLAG',FLAG,'NROOT',NBROOT STOP ENDIF CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3 NH2SO4=mrt_sa(ilon,ilev)*CONCM NH2O=mrt_wv(ilon,ilev)*CONCM CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev), & rho_droplet(ilon,ilev),TT(ilon,ilev), & H2SO4_liq,H2O_liq,MCONDTOT) ! Boucle sur les modes DO imode=1, nbr_mode IF (K_MASS(ilon,ilev,imode).GT.0.) THEN NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)* & K_MASS(ilon,ilev,imode)*MCONDTOT* & EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/ & (R_MEDIAN(ilon,ilev,imode)**3.) ELSE NBRTOT(ilon,ilev,imode)=0.0E+0 ENDIF ENDDO ! Passage de #/m3 en VMR H2O_liq=H2O_liq/CONCM H2SO4_liq=H2SO4_liq/CONCM mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq ! Problmes quand on a condense tout, on peut obtenir des -1e-24 ! aprs la soustraction et conversion de ND ˆ VMR IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30 IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30 CASE(3) ! Cas 0 NROOT mr_wv(ilon,ilev)=mrt_wv(ilon,ilev) mr_sa(ilon,ilev)=mrt_sa(ilon,ilev) rho_droplet(ilon,ilev)=0.0E+0 WH2SO4(ilon,ilev)=0.0E+0 DO imode=1, nbr_mode NBRTOT(ilon,ilev,imode)=0.0E+0 ENDDO END SELECT ENDDO !FIN boucle ilon ENDDO !FIN boucle ilev END SUBROUTINE new_cloud_venus !***************************************************************************** !* SUBROUTINE ITERWV() SUBROUTINE ITERWV(WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT, + TAIR,PAIR,RADIUS) !***************************************************************************** !* Cette routine est la solution par itŽration afin de trouver WSA pour un WV, !* et donc LPPWV, donnŽ. Ce qui nous donne egalement le WV correspondant au !* WSA solution !* For VenusGCM by A. Stolzenbach 07/2014 !* OUTPUT: WVEQ et WSAOUT IMPLICIT NONE REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS REAl, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ REAL :: WSAMIN, WSAMAX, WSAACC PARAMETER(WSAACC=0.001) REAL :: LPPWV INTEGER :: MAXITSA, NBRACSA, NBROOT PARAMETER(MAXITSA=20) PARAMETER(NBRACSA=20) LOGICAl :: FLAG1,FLAG2 ! External Function REAl :: IRFRMSA, WVCOND IF (RADIUS.LT.1E-30) THEN PRINT*,'RMODE == 0 FLAG 3' STOP ENDIF ! Initialisation WSA=[0.1,1.0] WSAMIN=0.1 WSAMAX=1.0 LPPWV=DLOG(PAIR*WV) ! Appel Bracket de KEEQ CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR, & LPPWV,FLAG1,FLAG2,NBROOT) IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN ! Appel Ridder's Method WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA, & RADIUS,TAIR,PAIR,LPPWV,NBROOT) ! IF (WSAOUT.EQ.1.0) WSAOUT=0.999999 ! IF (WSAOUT.LT.0.1) WSAOUT=0.1 ! Si BRACWSA ne trouve aucun ensemble solution KEEQ=0 on fixe WSA a 0.9999 ou 0.1 ELSE IF (FLAG1.AND.(.NOT.FLAG2)) WSAOUT=0.999999 IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT=WSAMIN IF (FLAG1.AND.FLAG2) THEN PRINT*,'FLAGs BARCWSA tous TRUE' STOP ENDIF ENDIF ! WVEQ output correspondant a WVliq lie a WSA calcule WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT) WVEQOUT=(WVLIQ+WV)/WVTOT-1.0 END SUBROUTINE ITERWV !***************************************************************************** !* SUBROUTINE BRACWV() SUBROUTINE BRACWV(XA,XB,N,RADIUS,WVTOT,SATOT,TAIR,PAIR, + FLAGWV,WSAFLAG,NROOT) !***************************************************************************** !* Bracket de ITERWV !* From Numerical Recipes !* Adapted for VenusGCM A. Stolzenbach 07/2014 !* X est WVinput !* OUTPUT: XA et XB IMPLICIT NONE REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR,PAIR INTEGER, INTENT(IN) :: N REAL, INTENT(INOUT) :: XA,XB REAL, INTENT(OUT) :: WSAFLAG INTEGER :: I,J INTEGER, INTENT(OUT) :: NROOT REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT REAL :: XMAX,XMIN,WVEQACC INTEGER, INTENT(OUT) :: FLAGWV ! WVEQACC est le seuil auquel on accorde un WSA correct meme ! si il ne fait pas partie d'une borne. Utile quand le modele ! s'approche de 0 mais ne l'atteint pas. WVEQACC=1.0E-3 FLAGWV=1 NROOT=0 X=XA XMAX=XB XMIN=XA ! CAS 1 On borne la fonction (WVEQ=0) CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS) FP=WVEQ DO I=1,N-1 X=(1.-DLOG(REAL(N-I))/DLOG(REAL(N)))*XMAX CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS) FC=WVEQ IF ((FP*FC).LT.0.0) THEN NROOT=NROOT+1 ! Si NROOT>1 on place la borne sup output ˆ la borne min du calcul en i IF (NROOT.GT.1) THEN XB=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX ENDIF IF (I.EQ.1) THEN XA=XMIN ELSE XA=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX ENDIF XB=X ENDIF FP=FC ENDDO ! CAS 2 on refait la boucle pour tester si WVEQ est proche de 0 ! avec le seuil WVEQACC IF (NROOT.EQ.0) THEN X=XMIN CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT, + TAIR,PAIR,RADIUS) DO J=1,N-1 X=(1.-DLOG(REAL(N-J))/DLOG(REAL(N)))*XMAX CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT, + TAIR,PAIR,RADIUS) IF (ABS(WVEQ).LE.WVEQACC) THEN WSAFLAG=WSAOUT FLAGWV=2 RETURN ENDIF ENDDO ! CAS 3 Pas de borne, WVEQ jamais proche de 0 FLAGWV=3 RETURN ENDIF END SUBROUTINE BRACWV !***************************************************************************** !* SUBROUTINE BRACWSA() SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL, + NROOT) !***************************************************************************** !* Bracket de KEEQ !* From Numerical Recipes !* Adapted for VenusGCM A. Stolzenbach 07/2014 IMPLICIT NONE !---------------------------------------------------------------------------- ! External functions needed: REAl KEEQ !---------------------------------------------------------------------------- REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP INTEGER, INTENT(IN) :: N REAL, INTENT(INOUT) :: XA,XB INTEGER, INTENT(OUT) :: NROOT INTEGER :: I, J REAL :: DX, FP, FC, X LOGICAL, INTENT(OUT) :: FLAGH,FLAGL FLAGL=.FALSE. FLAGH=.FALSE. NROOT=0 DX=(XB-XA)/N X=XA FP=KEEQ(RADIUS,TAIR,X,LPPWVINP) DO I=1,N X=X+DX FC=KEEQ(RADIUS,TAIR,X,LPPWVINP) IF ((FP*FC).LE.0.) THEN NROOT=NROOT+1 XA=X-DX XB=X ! RETURN ! IF (NROOT.GT.1) THEN ! PRINT*,'On a plus d1 intervalle KEEQ=0' ! PRINT*,'Probleme KEEQ=0 => 1 racine en theorie' ! X=X-(I*DX) ! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP) ! PRINT*,'KEEQ(WSA)',FP,X,TAIR ! DO J=1,N ! X=X+DX ! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP) ! PRINT*,'KEEQ(WSA)',FP,X ! ENDDO ! STOP ! ENDIF ENDIF FP=FC ENDDO IF (NROOT.EQ.0) THEN ! PRINT*,'On a 0 intervalle KEEQ=0' ! PRINT*,'Probleme KEEQ=0 => 1 racine en theorie' ! PRINT*,'XA',XA,'KEEQ',KEEQ(RADIUS,TAIR,XA,LPPWVINP) ! PRINT*,'XB',XB,'KEEQ',KEEQ(RADIUS,TAIR,XB,LPPWVINP) ! PRINT*,'TT',TAIR ! PRINT*,'RADIUS',RADIUS ! PRINT*,'NBRAC',N ! STOP ! X=XA ! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP) ! PRINT*,'KEEQ(WSA)',FP,X,TAIR ! DO I=1,N ! X=X+DX ! FP=KEEQ(RADIUS,TAIR,X,LPPWVINP) ! PRINT*,'KEEQ(WSA)',FP,X,TAIR ! ENDDO ! Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX] IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))- & ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE. ! On fixe flag low TRUE pour WSA = 0.1 IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))- & ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).LT.0.0) FLAGL=.TRUE. ! STOP ENDIF END SUBROUTINE BRACWSA !***************************************************************************** !* REAL FUNCTION WVCOND() REAL FUNCTION WVCOND(WSA,T,P,SAt) !***************************************************************************** !* Condensation de H2O selon WSA, T et P et H2SO4tot !* !* Adapted for VenusGCM A. Stolzenbach 07/2014 ! INPUT: ! SAt : VMR of total H2SO4 ! WSA: aerosol H2SO4 weight fraction (fraction) ! T: temperature (K) ! P: pressure (Pa) ! OUTPUT: ! WVCOND : VMR H2O condense ! USE chemparam_mod IMPLICIT NONE REAL, INTENT(IN) :: SAt, WSA REAL, INTENT(IN) :: T, P ! working variables REAL SA, WV REAL DND2,pstand,lpar,acidps REAL x1, satpacid REAL , DIMENSION(2):: act REAL CONCM REAL NH2SO4 REAL H2OCOND, H2SO4COND CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS! NH2SO4=SAt*CONCM pstand=1.01325E+5 !Pa 1 atm pressure x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153)) CALL zeleznik(x1,T,act) !pure acid satur vapor pressure lpar= -11.695+DLOG(pstand) ! Zeleznik acidps=1/360.15-1.0/T+0.38/545. & *(1.0+DLOG(360.15/T)-360.15/T) acidps = 10156.0*acidps +lpar acidps = DEXP(acidps) !Pa !acid sat.vap.PP over mixture (flat surface): satpacid=act(2)*acidps ! Pa ! Conversion from Pa to N.D #/m3 DND2=satpacid/(1.3806488E-23*T) ! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat IF (NH2SO4.GT.DND2) THEN H2SO4COND=NH2SO4-DND2 ! calcul de H2O cond correspondant a H2SO4 cond H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA) ! Si on a H2SO4H2O dispo ! IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN ! On peut alors condenser tout le H2O dispo ! H2OCOND=NH2O ! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond ! H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA)) ! END IF ! Calcul de H2O condensŽe VMR WVCOND=H2OCOND/CONCM END FUNCTION WVCOND !***************************************************************************** !* REAL FUNCTION IRFRMWV() REAL FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR, + WVTOT,SATOT,NROOT) !***************************************************************************** !* Iterative Root Finder Ridder's Method for Water Vapor calculus !* From Numerical Recipes !* Adapted for VenusGCM A. Stolzenbach 07/2014 !* !* Les iterations sur [X1,X2] sont [WV1,WV2] !* la variable X est WV !* IRFRMWV sort en OUTPUT : WSALOC pour ITERWV=0 (ou WVEQ=0) IMPLICIT NONE REAL, INTENT(IN) :: X1, X2 REAL, INTENT(IN) :: XACC INTEGER, INTENT(IN) :: MAXIT,NROOT ! LOCAL VARIABLES REAL :: XL, XH, XM, XNEW, X REAL :: WSALOC, WVEQ, WVLIQ REAL :: FL, FH, FM, FNEW REAL :: ANS, S, FSIGN INTEGER i ! External variables needed: REAL, INTENT(IN) :: TAIR,PAIR REAL, INTENT(IN) :: WVTOT,SATOT REAL, INTENT(IN) :: RADIUS ! Initialisation X=X1 CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS) FL=WVEQ X=X2 CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS) FH=WVEQ ! Test Bracketed values IF (((FL.LT.0.).AND.(FH.GT.0.)).OR. & ((FL.GT.0.).AND.(FH.LT.0.))) & THEN XL=X1 XH=X2 ANS=-9.99e99 DO i=1, MAXIT XM=0.5*(XL+XH) CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) FM=WVEQ S=SQRT(FM*FM-FL*FH) IF (S.EQ.0.0) THEN IRFRMWV=WSALOC RETURN ENDIF IF (FL.GT.FH) THEN FSIGN=1.0 ELSE FSIGN=-1.0 ENDIF XNEW=XM+(XM-XL)*(FSIGN*FM/S) IF (ABS(XNEW-ANS).LE.XACC) THEN IRFRMWV=WSALOC RETURN ENDIF ANS=XNEW CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) FNEW=WVEQ IF (FNEW.EQ.0.0) THEN IRFRMWV=WSALOC RETURN ENDIF IF (SIGN(FM, FNEW).NE.FM) THEN XL=XM FL=FM XH=ANS FH=FNEW ELSEIF (SIGN(FL, FNEW).NE.FL) THEN XH=ANS FH=FNEW ELSEIF (SIGN(FH, FNEW).NE.FH) THEN XL=ANS FL=FNEW ELSE PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS' STOP ENDIF ENDDO PRINT*,'Paaaaas bien MAXIT atteint' PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS' XL=X1 XH=X2 ANS=-9.99e99 DO i=1, MAXIT XM=0.5*(XL+XH) CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) FM=WVEQ S=SQRT(FM*FM-FL*FH) IF (FL.GT.FH) THEN FSIGN=1.0 ELSE FSIGN=-1.0 ENDIF XNEW=XM+(XM-XL)*(FSIGN*FM/S) ANS=XNEW CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) FNEW=WVEQ PRINT*,'WVliq',WVLIQ,'WVtot',WVTOT,'WVeq',WVEQ PRINT*,'WSA',WSALOC,'SAtot',SATOT PRINT*,'T',TAIR,'P',PAIR IF (SIGN(FM, FNEW).NE.FM) THEN XL=XM FL=FM XH=ANS FH=FNEW ELSEIF (SIGN(FL, FNEW).NE.FL) THEN XH=ANS FH=FNEW ELSEIF (SIGN(FH, FNEW).NE.FH) THEN XL=ANS FL=FNEW ELSE PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS TWIIICE???' STOP ENDIF ENDDO STOP ELSE PRINT*,'IRFRMWV must be bracketed' PRINT*,'NROOT de BRACWV', NROOT IF (ABS(FL).LT.XACC) THEN PRINT*,'IRFRMWV FL == 0',FL PRINT*,'X1',X1,'X2',X2,'FH',FH CALL ITERWV(X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) IRFRMWV=WSALOC RETURN ENDIF IF (ABS(FH).LT.XACC) THEN PRINT*,'IRFRMWV FH == 0',FH PRINT*,'X1',X1,'X2',X2,'FL',FL CALL ITERWV(X2,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT, & TAIR,PAIR,RADIUS) IRFRMWV=WSALOC RETURN ENDIF IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN PRINT*,'STOP dans IRFRMWV avec rien == 0' PRINT*,'X1',X1,'X2',X2 PRINT*,'Fcalc',FL,FH PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS STOP ENDIF IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC' PRINT*,FL,FH STOP ENDIF ENDIF ! FIN Test Bracketed values END FUNCTION IRFRMWV !***************************************************************************** !* REAL FUNCTION IRFRMSA() REAL FUNCTION IRFRMSA(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,LPPWV, + NB) !***************************************************************************** !* Iterative Root Finder Ridder's Method for Sulfuric Acid calculus !* From Numerical Recipes !* Adapted for VenusGCM A. Stolzenbach 07/2014 !* !* Les iterations sur [X1,X2] sont [WSA1,WSA2] !* la variable X est WSA !* IRFRMSA sort en OUTPUT : WSA pour KEEQ=0 IMPLICIT NONE REAL, INTENT(IN) :: X1, X2 REAL, INTENT(IN) :: XACC INTEGER, INTENT(IN) :: MAXIT, NB ! LOCAL VARIABLES REAL XL, XH, XM, XNEW REAL Fl, FH, FM, FNEW REAL ANS, S, FSIGN INTEGER i ! External variables needed: REAL, INTENT(IN) :: TAIR,PAIR REAL, INTENT(IN) :: LPPWV REAL, INTENT(IN) :: RADIUS ! External functions needed: REAL KEEQ ! Initialisation FL=KEEQ(RADIUS,TAIR,X1,LPPWV) FH=KEEQ(RADIUS,TAIR,X2,LPPWV) ! Test Bracketed values IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.((FL.GT.0.).AND.(FH.LT.0.))) & THEN XL=X1 XH=X2 ANS=-9.99e99 DO i=1, MAXIT XM=0.5*(XL+XH) FM=KEEQ(RADIUS,TAIR,XM,LPPWV) S=SQRT(FM*FM-FL*FH) IF (S.EQ.0.0) THEN IRFRMSA=ANS RETURN ENDIF IF (FL.GT.FH) THEN FSIGN=1.0 ELSE FSIGN=-1.0 ENDIF XNEW=XM+(XM-XL)*(FSIGN*FM/S) IF (ABS(XNEW-ANS).LE.XACC) THEN IRFRMSA=ANS RETURN ENDIF ANS=XNEW FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV) IF (FNEW.EQ.0.0) THEN IRFRMSA=ANS RETURN ENDIF IF (SIGN(FM, FNEW).NE.FM) THEN XL=XM FL=FM XH=ANS FH=FNEW ELSEIF (SIGN(FL, FNEW).NE.FL) THEN XH=ANS FH=FNEW ELSEIF (SIGN(FH, FNEW).NE.FH) THEN XL=ANS FL=FNEW ELSE PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS' STOP ENDIF ENDDO PRINT*,'Paaaaas bien MAXIT atteint' PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS' XL=X1 XH=X2 PRINT*,'Borne XL',XL,'XH',XH ANS=-9.99e99 DO i=1, MAXIT XM=0.5*(XL+XH) FM=KEEQ(RADIUS,TAIR,XM,LPPWV) S=SQRT(FM*FM-FL*FH) IF (FL.GT.FH) THEN FSIGN=1.0 ELSE FSIGN=-1.0 ENDIF XNEW=XM+(XM-XL)*(FSIGN*FM/S) ANS=XNEW FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV) PRINT*,'KEEQ result',FNEW,'T',TAIR,'R',RADIUS IF (SIGN(FM, FNEW).NE.FM) THEN XL=XM FL=FM XH=ANS FH=FNEW ELSEIF (SIGN(FL, FNEW).NE.FL) THEN XH=ANS FH=FNEW ELSEIF (SIGN(FH, FNEW).NE.FH) THEN XL=ANS FL=FNEW ELSE PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus' PRINT*,'you shall not PAAAAAASS' STOP ENDIF ENDDO STOP ELSE PRINT*,'IRFRMSA must be bracketed' IF (FL.EQ.0.0) THEN PRINT*,'IRFRMSA FL == 0',Fl IRFRMSA=X1 RETURN ENDIF IF (FH.EQ.0.0) THEN PRINT*,'IRFRMSA FH == 0',FH IRFRMSA=X2 RETURN ENDIF IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH PRINT*,'X1',X1,'X2',X2 PRINT*,'Kind F', KIND(FL), KIND(FH) PRINT*,'Kind X', KIND(X1), KIND(X2) PRINT*,'Logical: ',(SIGN(FL,FH).NE.FL) PRINT*,'Logical: ',(SIGN(FH,FL).NE.FH) PRINT*,'nb root BRACWSA',NB STOP ENDIF ENDIF ! FIN Test Bracketed values END function IRFRMSA !***************************************************************************** !* REAL FUNCTION KEEQ() REAL FUNCTION KEEQ(RADIUS,TAIR,WSA,LPPWV) !***************************************************************************** !* Kelvin Equation EQuality !* ln(PPWV_eq) - (2Mh2o sigma)/(R T r rho) - ln(ph2osa) = 0 !* IMPLICIT NONE REAL, INTENT(IN) :: RADIUS,TAIR,WSA,LPPWV ! Physical constants: REAL MH2O REAL RGAS PARAMETER( ! Molar weight of water (kg/mole) + MH2O=18.0153d-3, ! Universal gas constant (J/(mole K)) + RGAS=8.314462175d0) ! ! External functions needed: REAL PWVSAS_GV,SIGMADROPLET,RHODROPLET ! PWVSAS_GV: Natural logaritm of water vapor pressure over ! sulfuric acid solution ! SIGMADROPLET: Surface tension of sulfuric acid solution ! RHODROPLET: Density of sulfuric acid solution ! ! Auxiliary local variables: REAL C1 PARAMETER( + C1=2.0d0*MH2O/RGAS) KEEQ=LPPWV-C1*SIGMADROPLET(WSA,TAIR)/ & (TAIR*RADIUS*RHODROPLET(WSA,TAIR))- & PWVSAS_GV(TAIR,WSA) END FUNCTION KEEQ ***************************************************************************** * REAL FUNCTION PWVSAS_GV(TAIR,WSA) REAL FUNCTION PWVSAS_GV(TAIR,WSA) ***************************************************************************** * * Natural logaritm of saturated water vapor pressure over plane * sulfuric acid solution. * * Source: J.I.Gmitro & T.Vermeulen: A.I.Ch.E.J. 10,740,1964. * W.F.Giauque et al.: J. Amer. Chem. Soc. 82,62,1960. * * The formula of Gmitro & Vermeulen for saturation pressure * is used: * ln(p) = A ln(298/T) + B/T + C + DT * with values of A,B,C and D given by Gmitro & Vermeulen, * and calculated from partial molal properties given by Giauque et al. * * * * Input: TAIR: Temperature (K) * WSA: Weight fraction of H2SO4 [0;1] * Output: Natural logaritm of water vapor pressure * over sulfuric acid solution ( ln(Pa) ) * * * External functions needed for calculation of partial molal * properties of pure components at 25 ! as function of W. IMPLICIT NONE REAL :: CPH2O,ALH2O,FFH2O,LH2O * CPH2O: Partial molal heat capacity of sulfuric acid solution. * ALH2O: Temparature derivative of CPH2O * FFH2O: Partial molal free energy of sulfuric acid solution. * LH2O: Partial molal enthalpy of sulfuric acid * ! ! REAL, INTENT(IN) :: TAIR,WSA REAL :: ADOT,BDOT,CDOT,DDOT REAL :: RGAS,MMHGPA REAL :: K1,K2 REAL :: A,B,C,D,CP,L,F,ALFA ! Physical constants given by Gmitro & Vermeulen: PARAMETER( + ADOT=-3.67340, + BDOT=-4143.5, + CDOT=10.24353, + DDOT=0.618943d-3) PARAMETER( ! Gas constant (cal/(deg mole)): + RGAS=1.98726, ! Natural logarith of conversion factor between atm. and Pa: + MMHGPA=11.52608845, + K1=298.15, + K2=K1*K1/2.0) ! ! CP=CPH2O(WSA) F=-FFH2O(WSA) L=-LH2O(WSA) ALFA=ALH2O(WSA) ! A=ADOT+(CP-K1*ALFA)/RGAS B=BDOT+(L-K1*CP+K2*ALFA)/RGAS C=CDOT+(CP+(F-L)/K1)/RGAS D=DDOT-ALFA/(2.0d0*RGAS) ! ! WRITE(*,*) 'TAIR= ',TAIR,' WSA= ',WSA ! WRITE(*,*) 'CPH2O(WSA)= ',CP ! WRITE(*,*) 'ALFAH2O(WSA)= ',ALFA ! WRITE(*,*) 'FFH2O(WSA)= ',F ! WRITE(*,*) 'LH2O(WSA)= ',L ! PWVSAS_GV=A*DLOG(K1/TAIR)+B/TAIR+C+D*TAIR+MMHGPA END FUNCTION PWVSAS_GV ******************************************************************************* * REAL FUNCTION CPH2O(W) REAL FUNCTION CPH2O(W) ******************************************************************************* * * Relative partial molal heat capacity of water (cal/(deg mole) in * sulfuric acid solution, as a function of H2SO4 weight fraction [0;1], * calculated by cubic spline fitting. * * Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960. * IMPLICIT NONE INTEGER :: NPOINT,I PARAMETER(NPOINT=109) REAL, DIMENSION(NPOINT) :: WTAB(NPOINT),CPHTAB(NPOINT), + Y2(NPOINT),YWORK(NPOINT) REAL, INTENT(IN):: W REAL :: CPH LOGICAL :: FIRST DATA (WTAB(I),I=1,NPOINT)/ +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, +0.99908,0.99927,0.99945,0.99963,0.99982/ DATA (CPHTAB(I),I=1,NPOINT)/ + 17.996, 17.896, 17.875, 17.858, 17.840, 17.820, 17.800, 17.791, + 17.783, 17.777, 17.771, 17.769, 17.806, 17.891, 18.057, 18.248, + 18.429, 18.567, 18.613, 18.640, 18.660, 18.660, 18.642, 18.592, + 18.544, 18.468, 18.348, 18.187, 17.995, 17.782, 17.562, 17.352, + 17.162, 16.993, 16.829, 16.657, 16.581, 16.497, 16.405, 16.302, + 16.186, 16.053, 15.901, 15.730, 15.540, 15.329, 15.101, 14.853, + 14.586, 14.296, 13.980, 13.638, 13.274, 12.896, 12.507, 12.111, + 11.911, 11.711, 11.514, 11.320, 11.130, 10.940, 10.760, 10.570, + 10.390, 10.200, 10.000, 9.8400, 9.7600, 9.7900, 9.9500, 10.310, + 10.950, 11.960, 13.370, 15.060, 16.860, 18.550, 20.000, 21.170, + 22.030, 22.570, 22.800, 22.750, 22.420, 21.850, 21.120, 20.280, + 19.360, 18.350, 17.220, 15.940, 14.490, 12.840, 10.800, 9.8000, + 7.8000, 3.8000,0.20000,-5.4000,-7.0000,-8.8000,-10.900,-13.500, +-17.000,-22.000,-29.000,-40.000,-59.000/ DATA FIRST/.TRUE./ SAVE FIRST,WTAB,CPHTAB,Y2 ! IF(FIRST) THEN FIRST=.FALSE. CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2) ENDIF CALL SPLINT(WTAB,CPHTAB,Y2,NPOINT,W,CPH) CPH2O=CPH END FUNCTION CPH2O ! ******************************************************************************* REAL FUNCTION FFH2O(W) * REAL FUNCTION FFH2O(W) ******************************************************************************* * * Relative partial molal free energy water (cal/mole) in * sulfuric acid solution, as a function of H2SO4 weight fraction [0;1], * calculated by cubic spline fitting. * * Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960. * IMPLICIT NONE INTEGER :: NPOINT,I PARAMETER(NPOINT=110) REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK REAL, INTENT(IN) :: W REAL :: FF LOGICAL :: FIRST DATA (WTAB(I),I=1,NPOINT)/ +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/ DATA (FFTAB(I),I=1,NPOINT)/ +0.00000, 22.840, 25.810, 29.250, 33.790, 39.970, 48.690, 54.560, + 61.990, 71.790, 85.040, 103.70, 130.70, 145.20, 163.00, 184.50, + 211.50, 245.60, 266.40, 290.10, 317.40, 349.00, 385.60, 428.40, + 452.50, 478.80, 507.50, 538.80, 573.30, 611.60, 653.70, 700.50, + 752.60, 810.60, 875.60, 948.60, 980.60, 1014.3, 1049.7, 1087.1, + 1126.7, 1168.7, 1213.5, 1261.2, 1312.0, 1366.2, 1424.3, 1486.0, + 1551.8, 1622.3, 1697.8, 1778.5, 1864.9, 1956.8, 2055.8, 2162.0, + 2218.0, 2276.0, 2337.0, 2400.0, 2466.0, 2535.0, 2607.0, 2682.0, + 2760.0, 2842.0, 2928.0, 3018.0, 3111.0, 3209.0, 3311.0, 3417.0, + 3527.0, 3640.0, 3757.0, 3878.0, 4002.0, 4130.0, 4262.0, 4397.0, + 4535.0, 4678.0, 4824.0, 4973.0, 5128.0, 5287.0, 5454.0, 5630.0, + 5820.0, 6031.0, 6268.0, 6541.0, 6873.0, 7318.0, 8054.0, 8284.0, + 8579.0, 8997.0, 9295.0, 9720.0, 9831.0, 9954.0, 10092., 10248., + 10423., 10618., 10838., 11099., 11460., 12014./ DATA FIRST/.TRUE./ SAVE FIRST,WTAB,FFTAB,Y2 ! IF(FIRST) THEN FIRST=.FALSE. CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2) ENDIF CALL SPLINT(WTAB,FFTAB,Y2,NPOINT,W,FF) FFH2O=FF END FUNCTION FFH2O ! ******************************************************************************* REAL FUNCTION LH2O(W) * REAL FUNCTION LH2O(W) ******************************************************************************* * * Relative partial molal heat content of water (cal/mole) in * sulfuric acid solution, as a function of H2SO4 weight fraction [0;1], * calculated by cubic spline fitting. * * Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960. * IMPLICIT NONE INTEGER :: NPOINT,I PARAMETER(NPOINT=110) REAL, DIMENSION(NPOINT) :: WTAB,LTAB,Y2,YWORK REAL, INTENT(IN) :: W REAL :: L LOGICAL :: FIRST DATA (WTAB(I),I=1,NPOINT)/ +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525, +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209, +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/ DATA (LTAB(I),I=1,NPOINT)/ +0.00000, 5.2900, 6.1000, 7.1800, 8.7800, 11.210, 15.290, 18.680, + 23.700, 31.180, 42.500, 59.900, 89.200, 106.70, 128.60, 156.00, + 190.40, 233.80, 260.10, 290.00, 324.00, 362.50, 406.50, 456.10, + 483.20, 512.40, 543.60, 577.40, 613.80, 653.50, 696.70, 744.50, + 797.20, 855.80, 921.70, 995.70, 1028.1, 1062.3, 1098.3, 1136.4, + 1176.7, 1219.3, 1264.7, 1313.0, 1364.3, 1418.9, 1477.3, 1539.9, + 1607.2, 1679.7, 1757.9, 1842.7, 1934.8, 2035.4, 2145.5, 2267.0, + 2332.0, 2401.0, 2473.0, 2550.0, 2631.0, 2716.0, 2807.0, 2904.0, + 3007.0, 3118.0, 3238.0, 3367.0, 3507.0, 3657.0, 3821.0, 3997.0, + 4186.0, 4387.0, 4599.0, 4819.0, 5039.0, 5258.0, 5476.0, 5694.0, + 5906.0, 6103.0, 6275.0, 6434.0, 6592.0, 6743.0, 6880.0, 7008.0, + 7133.0, 7255.0, 7376.0, 7497.0, 7618.0, 7739.0, 7855.0, 7876.0, + 7905.0, 7985.0, 8110.0, 8415.0, 8515.0, 8655.0, 8835.0, 9125.0, + 9575.0, 10325., 11575., 13500., 15200., 16125./ DATA FIRST/.TRUE./ SAVE FIRST,WTAB,LTAB,Y2 ! IF(FIRST) THEN FIRST=.FALSE. CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2) ENDIF CALL SPLINT(WTAB,LTAB,Y2,NPOINT,W,L) LH2O=L END FUNCTION LH2O ******************************************************************************* REAL FUNCTION ALH2O(W) * REAL FUNCTION ALH2O(W) ******************************************************************************* * * Relative partial molal temperature derivative of heat capacity (water) * in sulfuric acid solution, (cal/deg**2), calculated by * cubic spline fitting. * * Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960. * IMPLICIT NONE INTEGER :: NPOINT,I PARAMETER(NPOINT=96) REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK REAL, INTENT(IN) :: W REAL :: A LOGICAL :: FIRST DATA (WTAB(I),I=1,NPOINT)/ +0.29517,0.31209, +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749, +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126, +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195, +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037, +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133, +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286, +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939, +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188, +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156, +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270, +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890, +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/ DATA (ATAB(I),I=1,NPOINT)/ + 0.0190, 0.0182, 0.0180, 0.0177, 0.0174, 0.0169, 0.0167, 0.0164, + 0.0172, 0.0212, 0.0239, 0.0264, 0.0276, 0.0273, 0.0259, 0.0238, + 0.0213, 0.0190, 0.0170, 0.0155, 0.0143, 0.0133, 0.0129, 0.0124, + 0.0120, 0.0114, 0.0106, 0.0097, 0.0084, 0.0067, 0.0047, 0.0024, +-0.0002,-0.0031,-0.0063,-0.0097,-0.0136,-0.0178,-0.0221,-0.0263, +-0.0303,-0.0340,-0.0352,-0.0360,-0.0362,-0.0356,-0.0343,-0.0321, +-0.0290,-0.0251,-0.0201,-0.0137,-0.0058, 0.0033, 0.0136, 0.0254, + 0.0388, 0.0550, 0.0738, 0.0962, 0.1198, 0.1300, 0.1208, 0.0790, + 0.0348, 0.0058,-0.0102,-0.0211,-0.0292,-0.0350,-0.0390,-0.0418, +-0.0432,-0.0436,-0.0429,-0.0411,-0.0384,-0.0346,-0.0292,-0.0220, +-0.0130,-0.0110,-0.0080,-0.0060,-0.0040,-0.0030,-0.0030,-0.0020, +-0.0020,-0.0020,-0.0020,-0.0010,-0.0010, 0.0000, 0.0000, 0.0000/ DATA FIRST/.TRUE./ SAVE FIRST,WTAB,ATAB,Y2 ! IF(FIRST) THEN FIRST=.FALSE. CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2) ENDIF CALL SPLINT(WTAB,ATAB,Y2,NPOINT,MAX(WTAB(1),W),A) ALH2O=A END FUNCTION ALH2O !****************************************************************************** SUBROUTINE SPLINE(X,Y,N,WORK,Y2) !****************************************************************************** ! Routine to calculate 2.nd derivatives of tabulated function ! Y(i)=Y(Xi), to be used for cubic spline calculation. ! IMPLICIT NONE INTEGER, INTENT(IN) :: N INTEGER :: I REAL, DIMENSION(N), INTENT(IN) :: X,Y REAL, DIMENSION(N), INTENT(OUT) :: Y2,WORK REAL SIG,P,QN,UN,YP1,YPN !AM Venus: Let's check the values ! write(*,*) 'In spline, N ', N YP1=(Y(2)-Y(1))/(X(2)-X(1)) YPN=(Y(N)-Y(N-1))/(X(N)-X(N-1)) IF(YP1.GT.99.0E+30) THEN Y2(1)=0.0 WORK(1)=0.0 ELSE Y2(1)=-0.5d0 WORK(1)=(3.0d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO I=2,N-1 ! write(*,*) 'In spline, I ', I SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2.0d0 Y2(I)=(SIG-1.0d0)/P WORK(I)=(6.0d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) + /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*WORK(I-1))/P ENDDO IF(YPN.GT.99.0E+30) THEN QN=0.0 UN=0.0 ELSE QN=0.5d0 UN=(3.0d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*WORK(N-1))/(QN*Y2(N-1)+1.0d0) DO I=N-1,1,-1 ! write(*,*) 'In spline, I ', I Y2(I)=Y2(I)*Y2(I+1)+WORK(I) ENDDO ! END SUBROUTINE SPLINE !****************************************************************************** SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) !****************************************************************************** ! Cubic spline calculation IMPLICIT NONE INTEGER, INTENT(IN) :: N INTEGER :: KLO,KHI,K REAL, INTENT(IN), DIMENSION(N) :: XA,YA,Y2A REAL, INTENT(IN) :: X REAL, INTENT(OUT) :: Y REAL :: H,A,B ! KLO=1 KHI=N 1 IF(KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X) THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ + ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0d0 ! END SUBROUTINE SPLINT !****************************************************************** SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4, + T,H2SO4COND,H2OCOND,RMTOT) ! DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION) ! FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL ! SIZE DISTRIBTUION ! ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED ! --------------------------------------------------------------- ! INPUT: ! H2SO4: #/m3 of total H2SO4 ! H2O : #/m3 of total H2O ! WSA: aerosol H2SO4 weight fraction (fraction) ! DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume) ! for total mass, almost same result with ro=1.67 gr/cm3 ! RSTDEV: standard deviation of aerosol distribution (no unit) ! RADIUS: MEDIAN radius (m) ! T: temperature (K) ! ! OUTPUT: ! RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension ! mais rho_droplet et M_tot_distrib doivent tre de meme dimension ! H2OCOND ! H2SO4COND IMPLICIT NONE REAL, INTENT(IN) :: H2SO4, H2O, WSA REAL, INTENT(IN) :: DENSO4, T REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT ! working variables REAL :: RMH2S4 REAL :: DND2,pstand,lpar,acidps REAL :: x1, satpacid REAL , DIMENSION(2):: act ! ! masse of an H2SO4 molecule (kg) RMH2S4=98.078/(6.02214129E+26) pstand=1.01325E+5 !Pa 1 atm pressure x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153)) call zeleznik(x1,t,act) !pure acid satur vapor pressure lpar= -11.695+DLOG(pstand) ! Zeleznik acidps=1/360.15-1.0/t+0.38/545. + *(1.0+DLOG(360.15/t)-360.15/t) acidps = 10156.0*acidps +lpar acidps = DEXP(acidps) !Pa !acid sat.vap.PP over mixture (flat surface): satpacid=act(2)*acidps ! Pa ! Conversion from Pa to N.D #/m3 DND2=satpacid/(1.3806488E-23*T) ! Conversion from N.D #/m3 TO #/cm3 ! DND2=DND2*1.d-6 ! H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat IF (H2SO4.GE.DND2) THEN H2SO4COND=H2SO4-DND2 ! calcul de H2O cond correspondant a H2SO4 cond H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA) ! RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3 ! RMTOT=M_distrib/rho_droplet RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA) ! Si on a H2SO4H2O dispo IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN ! Si H2O en dŽfaut, on as pas le bon WSA! ! En effet, normalement, on a exactement le WSA correspondant a ! WVg + WVl = WVtot ! Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou ! quand T est grand (Lower Haze), le modle reprŽsente mal le WSA ! cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA ! PRINT*,'PROBLEM H2O EN DEFAUT' ! PRINT*,'H2OCOND',H2OCOND,'H2O',H2O ! PRINT*,'WSA',WSA,'RHO',DENSO4 ! STOP ! On peut alors condenser tout le H2O dispo H2OCOND=H2O ! On met alors egalement a jour le H2SO4 cond correspondant au H2O cond H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA)) ! RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3 ! RMTOT=Volume of aerosol cm3 /cm3 of air ! Volume of aerosol/cm3 air RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA) END IF END SUBROUTINE CALCM_SAT SUBROUTINE Zeleznik(x,T,act) !+++++++++++++++++++++++++++++++++++++++++++++++++++ ! Water and sulfuric acid activities in liquid ! aqueous solutions. ! Frank J. Zeleznik, Thermodynnamic properties ! of the aqueous sulfuric acid system to 220K-350K, ! mole fraction 0,...,1 ! J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991 !+++++++++++++++++++++++++++++++++++++++++++++++++++ IMPLICIT NONE REAL, INTENT(IN) :: x,T REAL :: activitya, activityw REAL, INTENT(OUT), DIMENSION(2):: act ! REAL x,T, activitya, activityw ! REAL, DIMENSION(2):: act ! write(*,*) 'x, T ', x, T act(2)=activitya(x,T) act(1)=activityw(x,T) ! write(*,*) 'act ', act END SUBROUTINE Zeleznik !start of functions related to zeleznik activities REAL FUNCTION m111(T) REAL, INTENT(IN) :: T m111=-23.524503387D0 & +0.0406889449841D0*T & -0.151369362907D-4*T**2+2961.44445015D0/T & +0.492476973663D0*dlog(T) END FUNCTION m111 REAL FUNCTION m121(T) REAL, INTENT(IN) :: T m121=1114.58541077D0-1.1833078936D0*T & -0.00209946114412D0*T**2-246749.842271D0/T & +34.1234558134D0*dlog(T) END FUNCTION m121 FUNCTION m221(T) REAL, INTENT(IN) :: T m221=-80.1488100747D0-0.0116246143257D0*T & +0.606767928954D-5*T**2+3092.72150882D0/T & +12.7601667471D0*dlog(T) END FUNCTION m221 REAL FUNCTION m122(T) REAL, INTENT(IN) :: T m122=888.711613784D0-2.50531359687D0*T & +0.000605638824061D0*T**2-196985.296431D0/T & +74.550064338D0*dlog(T) END FUNCTION m122 REAL FUNCTION e111(T) REAL, INTENT(IN) :: T e111=2887.31663295D0-3.32602457749D0*T & -0.2820472833D-2*T**2-528216.112353D0/T & +0.68699743564D0*dlog(T) END FUNCTION e111 REAL FUNCTION e121(T) REAL, INTENT(IN) :: T e121=-370.944593249D0-0.690310834523D0*T & +0.56345508422D-3*T**2-3822.52997064D0/T & +94.2682037574D0*dlog(T) END FUNCTION e121 REAL FUNCTION e211(T) REAL, INTENT(IN) :: T e211=38.3025318809D0-0.0295997878789D0*T & +0.120999746782D-4*T**2-3246.97498999D0/T & -3.83566039532D0*dlog(T) END FUNCTION e211 REAL FUNCTION e221(T) REAL, INTENT(IN) :: T e221=2324.76399402D0-0.141626921317D0*T & -0.00626760562881D0*T**2-450590.687961D0/T & -61.2339472744D0*dlog(T) END FUNCTION e221 REAL FUNCTION e122(T) REAL, INTENT(IN) :: T e122=-1633.85547832D0-3.35344369968D0*T & +0.00710978119903D0*T**2+198200.003569D0/T & +246.693619189D0*dlog(T) END FUNCTION e122 REAL FUNCTION e212(T) REAL, INTENT(IN) :: T e212=1273.75159848D0+1.03333898148D0*T & +0.00341400487633D0*T**2+195290.667051D0/T & -431.737442782D0*dlog(T) END FUNCTION e212 REAL FUNCTION lnAa(x1,T) REAL, INTENT(IN) :: T,x1 REAL :: & m111,m121,m221,m122 & ,e111,e121,e211,e122,e212,e221 lnAa=-( & (2*m111(T)+e111(T)*(2*dlog(x1)+1))*x1 & +(2*m121(T)+e211(T)*dlog(1-x1)+e121(T)*(dlog(x1)+1))*(1-x1) & -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 & -(2*m121(T)+e121(T)*(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1) & -(2*m122(T)+e122(T)*dlog(x1) & +e212(T)*dlog(1-x1))*(1-x1))*x1*(1-x1) & -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 & -x1*(1-x1)*( & (6*m122(T)+e122(T)*(3*dlog(x1)+1) & +e212(T)*(3*dlog(1-x1)+1) & )*x1*(1-x1) & -(2*m122(T)+e122(T)*(dlog(x1)+1) & +e212(T)*dlog(1-x1) & )*(1-x1)) & ) END FUNCTION lnAa REAL FUNCTION lnAw(x1,T) REAL, INTENT(IN) :: T,x1 REAL :: & m111,m121,m221,m122 & ,e111,e121,e211,e122,e212,e221 lnAw=-( & (2*m121(T)+e121(T)*dlog(x1)+e211(T)*(dlog(1-x1)+1))*x1 & +(2*m221(T)+e221(T)*(2*dlog(1-x1)+1))*(1-x1) & -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 & -(2*m121(T)+e121(T)*(dlog(x1)+1) & +e211(T)*(dlog(1-x1)+1))*x1*(1-x1) & -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 & +x1*(2*m122(T)+e122(T)*dlog(x1)+e212(T)*dlog(1-x1))*x1*(1-x1) & +x1*(1-x1)*((2*m122(T)+e122(T)*dlog(x1) & +e212(T)*(dlog(1-x1)+1))*x1 & -(6*m122(T)+e122(T)*(3*dlog(x1)+1) & +e212(T)*(3*dlog(1-x1)+1))*(1-x1)*x1) & ) END FUNCTION lnAw REAL FUNCTION activitya(xal,T) REAL, INTENT(IN) :: T,xal REAL :: lnAa ! & ,m111,m121,m221,m122 & ! & ,e111,e121,e211,e122,e212,e221 ! write(*,*) 'in activitya ', xal, T activitya=DEXP(lnAa(xal,T)-lnAa(1.D0-1.D-12,T)) END FUNCTION activitya FUNCTION activityw(xal,T) REAL, INTENT(IN) :: T,xal REAL :: lnAw activityw=DEXP(lnAw(xal,T)-lnAw(1.D-12,T)) END FUNCTION activityw ! end of functions related to zeleznik activities FUNCTION SIGMADROPLET(xmass,t) ! calculates the surface tension of the liquid in J/m^2 ! xmass=mass fraction of h2so4, t in kelvins ! about 230-323 K , x=0,...,1 !(valid down to the solid phase limit temp, which depends on molefraction) IMPLICIT NONE REAL :: SIGMADROPLET REAL, INTENT(IN):: xmass, t REAL :: a, b, t1, tc, xmole REAL, PARAMETER :: Msa=98.078 REAL, PARAMETER :: Mwv=18.0153 IF (t .LT. 305.15) THEN !low temperature surface tension ! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002, ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric !and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631 a=0.11864+xmass*(-0.11651+xmass*(0.76852+xmass* & (-2.40909+xmass*(2.95434-xmass*1.25852)))) b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass* & (0.007611235+xmass*(-0.00937386+xmass*0.00389722)))) SIGMADROPLET=a+t*b ELSE xmole = (xmass/Msa)*(1./((xmass/Msa)+(1.-xmass)/Mwv)) ! high temperature surface tension !H. Vehkam‰ki and M. Kulmala and K.E. J. lehtinen, 2003, !Modelling binary homogeneous nucleation of water-sulfuric acid vapours: ! parameterisation for high temperature emissions, () Environ. Sci. Technol., 37, 3392-3398 tc= 647.15*(1.0-xmole)*(1.0-xmole)+900.0*xmole*xmole+ & 3156.186*xmole*(1-xmole) !critical temperature t1=1.0-t/tc a= 0.2358+xmole*(-0.529+xmole*(4.073+xmole*(-12.6707+xmole* & (15.3552+xmole*(-6.3138))))) b=-0.14738+xmole*(0.6253+xmole*(-5.4808+xmole*(17.2366+xmole* & (-21.0487+xmole*(8.719))))) SIGMADROPLET=(a+b*t1)*t1**(1.256) END IF RETURN END FUNCTION SIGMADROPLET FUNCTION RHODROPLET(xmass,t) ! ! calculates the density of the liquid in kg/m^3 ! xmass=mass fraction of h2so4, t in kelvins ! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002, ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric !and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631 ! about 220-373 K , x=0,...,1 !(valid down to the solid phase limit temp, which depends on molefraction) IMPLICIT NONE REAL :: RHODROPLET REAL, INTENT(IN) :: xmass, t REAL :: a,b,c a=0.7681724+xmass*(2.1847140+xmass*(7.1630022+xmass* & (-44.31447+xmass* & (88.75606+xmass*(-75.73729+xmass*23.43228))))) b=1.808225e-3+xmass*(-9.294656e-3+xmass*(-0.03742148+ & xmass*(0.2565321+xmass*(-0.5362872+xmass* & (0.4857736-xmass*0.1629592))))) c=-3.478524e-6+xmass*(1.335867e-5+xmass* & (5.195706e-5+xmass*(-3.717636e-4+xmass* & (7.990811e-4+xmass*(-7.458060e-4+xmass*2.58139e-4))))) RHODROPLET=a+t*(b+c*t) ! g/cm^3 RHODROPLET= RHODROPLET*1.0e3 !kg/m^3 RETURN END FUNCTION RHODROPLET