Changeset 1661


Ignore:
Timestamp:
Feb 8, 2017, 4:40:26 PM (8 years ago)
Author:
slebonnois
Message:

SL: Cloud model for Venus. Not validated yet.

Location:
trunk
Files:
14 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/makelmdz

    r1650 r1661  
    376376   OPTIMC="${coptim}"
    377377   INCLUDEC='-I$(LIBF)/grid -I.'
     378fi
     379#########
     380
     381######### CAS PARTICULIER NUAGES VENUS
     382if [[ "$physique" == "venus" ]]
     383then
     384   src_dirs="$src_dirs phy${physique}/cloudvenus"
    378385fi
    379386#########

  • trunk/LMDZ.COMMON/makelmdz_fcm

    r1650 r1661  
    5959COSP_PATH=$LMDGCM/.void_dir
    6060CHEM_PATH=$LMDGCM/.void_dir
     61CLOUD_PATH=$LMDGCM/.void_dir
    6162AERONO_PATH=$LMDGCM/.void_dir
    6263# Path to fcm utility:
    … …  
    473474fi
    474475
     476# for Venus (but could be used by others as well), there is also "phyvenus/cloudvenus"
     477if [[ -d ${LIBFGCM}/phy${physique}/cloud${physique} ]]
     478then
     479   CLOUD_PATH="${LIBFGCM}/phy${physique}/cloud${physique}"
     480   INCLUDE="$INCLUDE -I${LIBFGCM}/phy${physique}/cloud${physique}"
     481fi
     482
    475483# for Mars (but could be used by others as well), there is also "aeronomars"
    476484if [[ -d ${LIBFGCM}/aerono${physique} ]]
    … …  
    715723echo "%COSP          $COSP_PATH"     >> $config_fcm
    716724echo "%CHEM          $CHEM_PATH"     >> $config_fcm
     725echo "%CLOUD         $CLOUD_PATH"    >> $config_fcm
    717726echo "%AERONO        $AERONO_PATH"   >> $config_fcm
    718727echo "%CPP_KEY       $CPP_KEY"       >> $config_fcm

  • trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90

    r1621 r1661  
    11MODULE chemparam_mod
    22
    3 !MODULE qui définit les indices des traceurs et leurs masses molaires.
     3!MODULE qui definit les indices des traceurs et leurs masses molaires.
    44! utilise aussi pour variables communes nuages/photochimie
    55
    … …  
    1313                 i_cocl2, i_s, i_so, i_so2, i_so3,       &
    1414                 i_s2o2, i_ocs, i_hso3, i_h2so4, i_s2,   &
    15                  i_clso2, i_oscl, i_h2oliq, i_h2so4liq,  &
    16                  i_n2
     15                 i_clso2, i_oscl, i_n2
    1716                 
     17INTEGER, SAVE :: i_h2oliq, i_h2so4liq
     18
     19INTEGER, SAVE :: i_m0_aer, i_m3_aer,                       &
     20                 i_m0_mode1drop, i_m0_mode1ccn,            &
     21                 i_m3_mode1sa, i_m3_mode1w, i_m3_mode1ccn, &
     22                 i_m0_mode2drop, i_m0_mode2ccn,            &
     23                 i_m3_mode2sa, i_m3_mode2w, i_m3_mode2ccn
     24
     25INTEGER, SAVE :: nmicro  ! number of microphysical tracers
     26
    1827REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr
    1928 
    2029
    2130!----------------------------------------------------------------------------
     31! DEF FOR CL_SCHEME = 1 (AURELIEN)
     32
    2233!     number of clouds mode modelized
    2334      INTEGER, PARAMETER :: nbr_mode = 3
    … …  
    3445      REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: WH2SO4
    3546      REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: rho_droplet
     47!----------------------------------------------------------------------------
     48! DEF FOR CL_SCHEME = 2 (FULL MICROPHYS)
     49
    3650!----------------------------------------------------------------------------
    3751
    … …  
    670684  END SUBROUTINE cloud_ini
    671685 
     686! ===========================================================
     687
    672688  SUBROUTINE chemparam_ini
    673689  USE infotrac_phy, ONLY: nqtot, tname
    … …  
    811827                PRINT*,'oscl',i_oscl
    812828                M_tr(i_oscl)=83.517
     829                CASE('n2')
     830                i_n2=i
     831                M_tr(i_n2)=28.013
     832! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
    813833                CASE('h2oliq')
    814834                i_h2oliq=i
    … …  
    819839                PRINT*,'h2so4liq',i_h2so4liq
    820840                M_tr(i_h2so4liq)=98.078
    821                 CASE('n2')
    822                 i_n2=i
    823                 M_tr(i_n2)=28.013
     841! MICROPHYSICAL TRACERS FOR CL_SCHEME=2
     842                CASE('M0_aer')
     843                i_m0_aer=i
     844                PRINT*,'M0_aer',i_m0_aer
     845                CASE('M3_aer')
     846                i_m3_aer=i
     847                PRINT*,'M3_aer',i_m3_aer
     848                CASE('M0_m1drop')
     849                i_m0_mode1drop=i
     850                PRINT*,'M0_m1drop',i_m0_mode1drop
     851                CASE('M0_m1ccn')
     852                i_m0_mode1ccn=i
     853                PRINT*,'M0_m1ccn',i_m0_mode1ccn
     854                CASE('M3_m1sa')
     855                i_m3_mode1sa=i
     856                PRINT*,'M3_m1sa',i_m3_mode1sa
     857                CASE('M3_m1w')
     858                i_m3_mode1w=i
     859                PRINT*,'M3_m1w',i_m3_mode1w
     860                CASE('M3_m1ccn')
     861                i_m3_mode1ccn=i
     862                PRINT*,'M3_m1ccn',i_m3_mode1ccn
     863                CASE('M0_m2drop')
     864                i_m0_mode2drop=i
     865                PRINT*,'M0_m2drop',i_m0_mode2drop
     866                CASE('M0_m2ccn')
     867                i_m0_mode2ccn=i
     868                PRINT*,'M0_m2ccn',i_m0_mode2ccn
     869                CASE('M3_m2sa')
     870                i_m3_mode2sa=i
     871                PRINT*,'M3_m2sa',i_m3_mode2sa
     872                CASE('M3_m2w')
     873                i_m3_mode2w=i
     874                PRINT*,'M3_m2w',i_m3_mode2w
     875                CASE('M3_m2ccn')
     876                i_m3_mode2ccn=i
     877                PRINT*,'M3_m2ccn',i_m3_mode2ccn
    824878        END SELECT
    825879       
    … …  
    828882 
    829883  END SUBROUTINE chemparam_ini
    830  
     884
     885! ===========================================================
     886
     887  SUBROUTINE vapors4muphy_ini(nlon,nlev,trac)
     888  USE infotrac_phy, ONLY: nqtot, tname
     889  IMPLICIT NONE
     890
     891  integer :: nlon, nlev
     892  real    :: trac(nlon,nlev,nqtot) ! traceur ( en vmr)
     893
     894  integer :: i
     895  real    :: trac1d(nlev,2) ! traceur lu ( en vmr)
     896 
     897! lecture d'un fichier texte contenant les profils de trac1d(:1) = H2O et trac1d(:,2) = H2SO4
     898
     899  DO i=1,nlon
     900     trac(i,:,i_h2o) = trac1d(:,1)
     901     trac(i,:,i_h2so4) = trac1d(:,2)
     902  ENDDO
     903
     904
     905  END SUBROUTINE vapors4muphy_ini
     906
    831907END MODULE chemparam_mod
    832908

  • trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h

    r1642 r1661  
    1515       INTEGER nbapp_rad, nbapp_chim, iflag_con, iflag_ajs
    1616       INTEGER lev_histins, lev_histday, lev_histmth
    17        INTEGER tr_scheme
     17       INTEGER tr_scheme, cl_scheme
    1818       INTEGER nircorr, nltemodel, solvarmod
    1919       REAL    ecriphy
    … …  
    3232     &     iflag_con, iflag_ajs,                                        &
    3333     &     lev_histins, lev_histday, lev_histmth, tr_scheme,            &
    34      &     nircorr, nltemodel, solvarmod, nb_mode
     34     &     cl_scheme, nircorr, nltemodel, solvarmod, nb_mode
    3535
    3636       COMMON/clesphys_r/ ecriphy, solaire, z0, lmixmin,                &

  • trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90

    r1642 r1661  
    334334  ok_cloud = .FALSE.
    335335  call getin('ok_cloud',ok_cloud)
     336
     337!
     338!Config Key  = cl_scheme
     339!Config Desc =
     340!Config Def  = 2
     341!Config Help =
     342!
     343! 1   = Simple microphysics (Aurelien Stolzenbach's PhD)
     344! 2   = Full microphysics (momentum scheme, Sabrina Guilbon's PhD)
     345
     346  cl_scheme = 2
     347  call getin('cl_scheme',cl_scheme)
    336348
    337349!

  • trunk/LMDZ.VENUS/libf/phyvenus/ini_histins.h

    r1530 r1661  
    141141     .                "ins(X)", zsto,zout)
    142142c
    143           if (ok_cloud) THEN
     143          if (ok_cloud.and.(cl_scheme.eq.1)) THEN
    144144
    145145          if (nb_mode.GE.1) THEN
    … …  
    235235                ENDIF
    236236               
    237           if (ok_sedim) THEN
     237          if (ok_sedim.and.(cl_scheme.eq.1)) THEN
    238238c
    239239         CALL histdef(nid_ins, "d_tr_sed_H2SO4", "var mmr from sedim",
    … …  
    257257c
    258258                ENDIF
     259
    259260      ENDIF !lev_histins.GE.2
    260261c

  • trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F

    r1639 r1661  
    3434!----------------------------------------------------------------------------
    3535!     Thermodynamic functions:
    36       REAL :: RHODROPLET
     36      REAL :: ROSAS
    3737!----------------------------------------------------------------------------
    3838!     Auxilary variables:
    … …  
    134134     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
    135135
    136           rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
    137      &    TT(ilon,ilev))
     136          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
    138137
    139138!          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
    … …  
    141140!            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
    142141!            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
    143 !            PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev),
    144 !     &      TT(ilon,ilev))
     142!            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
    145143!            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
    146144!            STOP
    … …  
    204202          WH2SO4(ilon,ilev)=WSAFLAG
    205203
    206           rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
    207      &    TT(ilon,ilev))
     204          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
    208205
    209206!     ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
    … …  
    211208!     incoherentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
    212209!     Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
    213 !     donne par RHODROPLET (cf test externe en Python), sinon, la valeur est trop
     210!     donne par ROSAS (cf test externe en Python), sinon, la valeur est trop
    214211!     basse (de l'ordre de 1000 kg/m3) et correspond parfois a la valeur avec
    215212!     WSA=0.1 (pas totalement sur)
    … …  
    223220            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
    224221            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
    225             PRINT*,'RHODROPLET',RHODROPLET(WH2SO4(ilon,ilev),
    226      &      TT(ilon,ilev))
     222            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
    227223            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
    228224            STOP
    … …  
    298294
    299295
    300 !*****************************************************************************
    301 !*    SUBROUTINE ITERWV()
    302       SUBROUTINE ITERWV(WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,
    303      + TAIR,PAIR,RADIUS)
    304 !*****************************************************************************
    305 !* Cette routine est la solution par iteration afin de trouver WSA pour un WV,
    306 !* et donc LPPWV, donne. Ce qui nous donne egalement le WV correspondant au
    307 !* WSA solution
    308 !* For VenusGCM by A. Stolzenbach 07/2014
    309 !* OUTPUT: WVEQ et WSAOUT
    310 
    311       IMPLICIT NONE
    312       REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS
    313 
    314       REAl, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ
    315 
    316       REAL :: WSAMIN, WSAMAX, WSAACC
    317       PARAMETER(WSAACC=0.01)
    318 
    319       REAL :: LPPWV
    320 
    321       INTEGER :: MAXITSA, NBRACSA, NBROOT
    322       PARAMETER(MAXITSA=20)
    323       PARAMETER(NBRACSA=5)
    324 
    325       LOGICAl :: FLAG1,FLAG2
    326 
    327 !     External Function
    328       REAl :: IRFRMSA, WVCOND
    329 
    330       IF (RADIUS.LT.1E-30) THEN
    331         PRINT*,'RMODE == 0 FLAG 3'
    332         STOP
    333       ENDIF
    334 !     Initialisation WSA=[0.1,1.0]
    335       WSAMIN=0.1
    336       WSAMAX=1.0
    337 
    338       LPPWV=DLOG(PAIR*WV)
    339 
    340 !     Appel Bracket de KEEQ
    341       CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR,
    342      &     LPPWV,FLAG1,FLAG2,NBROOT)
    343 
    344       IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN
    345 !          Appel Ridder's Method
    346 
    347         WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA,
    348      &  RADIUS,TAIR,PAIR,LPPWV,NBROOT)
    349 !              IF (WSAOUT.EQ.1.0) WSAOUT=0.999999
    350 !              IF (WSAOUT.LT.0.1) WSAOUT=0.1
    351 
    352 !       Si BRACWSA ne trouve aucun ensemble solution KEEQ=0 on fixe WSA a 0.9999 ou 0.1
    353       ELSE
    354         IF (FLAG1.AND.(.NOT.FLAG2)) WSAOUT=0.999999
    355         IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT=WSAMIN
    356         IF (FLAG1.AND.FLAG2) THEN
    357           PRINT*,'FLAGs BRACWSA tous TRUE'
    358           STOP
    359         ENDIF
    360       ENDIF
    361 
    362 
    363 !     WVEQ output correspondant a WVliq lie a WSA calcule
    364       WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT)
    365       WVEQOUT=(WVLIQ+WV)/WVTOT-1.0
    366 
    367       END SUBROUTINE ITERWV
    368 
    369 
    370 !*****************************************************************************
    371 !*    SUBROUTINE BRACWV()
    372       SUBROUTINE BRACWV(XA,XB,N,RADIUS,WVTOT,SATOT,TAIR,PAIR,
    373      +           FLAGWV,WSAFLAG,NROOT)
    374 !*****************************************************************************
    375 !* Bracket de ITERWV
    376 !* From Numerical Recipes
    377 !* Adapted for VenusGCM A. Stolzenbach 07/2014
    378 !* X est WVinput
    379 !* OUTPUT: XA et XB
    380 
    381       IMPLICIT NONE
    382 
    383       REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR,PAIR
    384       INTEGER, INTENT(IN) :: N
    385 
    386       REAL, INTENT(INOUT) :: XA,XB
    387       REAL, INTENT(OUT) :: WSAFLAG
    388 
    389       INTEGER :: I,J
    390 
    391       INTEGER, INTENT(OUT) :: NROOT
    392 
    393       REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT
    394       REAL :: XMAX,XMIN,WVEQACC
    395 
    396       INTEGER, INTENT(OUT) :: FLAGWV
    397 
    398 !     WVEQACC est le seuil auquel on accorde un WSA correct meme
    399 !     si il ne fait pas partie d'une borne. Utile quand le modele
    400 !     s'approche de 0 mais ne l'atteint pas.
    401       WVEQACC=1.0E-3
    402 
    403       FLAGWV=1
    404 
    405       NROOT=0
    406 
    407 !     25/11/2016
    408 !     On change ordre on va du max au min
    409       X=XB
    410       XMAX=XB
    411       XMIN=XA
    412 
    413 !     CAS 1 On borne la fonction (WVEQ=0)
    414 
    415       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
    416       FP=WVEQ
    417 
    418       DO I=N-1,1,-1
    419          X=(1.-DLOG(REAL(N-I))/DLOG(REAL(N)))*XMAX
    420          CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
    421          FC=WVEQ
    422 
    423           IF ((FP*FC).LT.0.0) THEN
    424             NROOT=NROOT+1
    425 !           Si NROOT>1 on place la borne sup output ˆ la borne min du calcul en i
    426             IF (NROOT.GT.1) THEN
    427                XB=(1.-DLOG(REAL(I+1))/DLOG(REAL(N)))*XMAX
    428             ENDIF
    429 
    430             IF (I.EQ.1) THEN
    431               XA=XMIN
    432             ELSE
    433               XA=X
    434             ENDIF
    435             RETURN
    436          ENDIF
    437          FP=FC
    438       ENDDO
    439 
    440 !     CAS 2 on refait la boucle pour tester si WVEQ est proche de 0
    441 !     avec le seuil WVEQACC
    442       IF (NROOT.EQ.0) THEN
    443          X=XMAX
    444          CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
    445      +   TAIR,PAIR,RADIUS)
    446           DO J=N-1,1,-1
    447              X=(1.-DLOG(REAL(N-J))/DLOG(REAL(N)))*XMAX
    448              CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
    449      +       TAIR,PAIR,RADIUS)
    450 
    451              IF (ABS(WVEQ).LE.WVEQACC) THEN
    452                 WSAFLAG=WSAOUT
    453                 FLAGWV=2
    454                 RETURN
    455              ENDIF
    456           ENDDO
    457 
    458 !     CAS 3 Pas de borne, WVEQ jamais proche de 0
    459           FLAGWV=3
    460           RETURN
    461        ENDIF
    462 
    463       END SUBROUTINE BRACWV
    464 
    465 !*****************************************************************************
    466 !*    SUBROUTINE BRACWSA()
    467       SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL,
    468      +           NROOT)
    469 !*****************************************************************************
    470 !* Bracket de KEEQ
    471 !* From Numerical Recipes
    472 !* Adapted for VenusGCM A. Stolzenbach 07/2014
    473 
    474       IMPLICIT NONE
    475 
    476 !----------------------------------------------------------------------------
    477 !     External functions needed:
    478       REAl KEEQ
    479 !----------------------------------------------------------------------------
    480 
    481       REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP
    482       INTEGER, INTENT(IN) :: N
    483 
    484       REAL, INTENT(INOUT) :: XA,XB
    485 
    486       INTEGER, INTENT(OUT) ::  NROOT
    487 
    488       INTEGER :: I, J
    489 
    490       REAL :: DX, FP, FC, X
    491 
    492       LOGICAL, INTENT(OUT) :: FLAGH,FLAGL
    493 
    494 
    495       FLAGL=.FALSE.
    496       FLAGH=.FALSE.
    497       NROOT=0
    498       DX=(XB-XA)/N
    499       X=XB
    500       FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    501 
    502       DO I=N,1,-1
    503          X=X-DX
    504          FC=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    505 
    506          IF ((FP*FC).LE.0.) THEN
    507             NROOT=NROOT+1
    508             XA=X
    509             XB=X+DX
    510             RETURN
    511 !            IF (NROOT.GT.1) THEN
    512 !               PRINT*,'On a plus d1 intervalle KEEQ=0'
    513 !               PRINT*,'Probleme KEEQ=0 => 1 racine en theorie'
    514 !               X=X-(I*DX)
    515 !               FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    516 !               PRINT*,'KEEQ(WSA)',FP,X,TAIR
    517 !               DO J=1,N
    518 !                 X=X+DX
    519 !                FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    520 !                PRINT*,'KEEQ(WSA)',FP,X
    521 !               ENDDO
    522 !               STOP
    523 !            ENDIF
    524          ENDIF
    525 
    526          FP=FC
    527       ENDDO
    528 
    529       IF (NROOT.EQ.0) THEN
    530 !         PRINT*,'On a 0 intervalle KEEQ=0'
    531 !         PRINT*,'Probleme KEEQ=0 => 1 racine en theorie'
    532 !         PRINT*,'XA',XA,'KEEQ',KEEQ(RADIUS,TAIR,XA,LPPWVINP)
    533 !         PRINT*,'XB',XB,'KEEQ',KEEQ(RADIUS,TAIR,XB,LPPWVINP)
    534 !         PRINT*,'TT',TAIR
    535 !         PRINT*,'RADIUS',RADIUS
    536 !         PRINT*,'NBRAC',N
    537 !         STOP
    538 
    539 !         X=XA
    540 !         FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    541 !         PRINT*,'KEEQ(WSA)',FP,X,TAIR
    542 !         DO I=1,N
    543 !           X=X+DX
    544 !           FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    545 !           PRINT*,'KEEQ(WSA)',FP,X,TAIR
    546 !         ENDDO
    547 
    548 
    549 !        Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]
    550          IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
    551      &    ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE.
    552 !        On fixe flag low TRUE pour WSA = 0.1
    553          IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
    554      &    ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).LT.0.0) FLAGL=.TRUE.
    555 !        STOP
    556       ENDIF
    557 
    558       END SUBROUTINE BRACWSA
    559 
    560 
    561 !*****************************************************************************
    562 !*     REAL FUNCTION WVCOND()
    563       REAL FUNCTION WVCOND(WSA,T,P,SAt)
    564 !*****************************************************************************
    565 !* Condensation de H2O selon WSA, T et P et H2SO4tot
    566 !*
    567 !* Adapted for VenusGCM A. Stolzenbach 07/2014
    568 !     INPUT:
    569 !     SAt  : VMR of total H2SO4
    570 !     WSA: aerosol H2SO4 weight fraction (fraction)
    571 !     T: temperature (K)
    572 !     P: pressure (Pa)
    573 !     OUTPUT:
    574 !       WVCOND : VMR H2O condense
    575 
    576 !      USE chemparam_mod
    577 
    578       IMPLICIT NONE
    579 
    580       REAL, INTENT(IN) :: SAt, WSA
    581       REAL, INTENT(IN) :: T, P
    582 
    583 !     working variables
    584       REAL SA, WV
    585       REAL  DND2,pstand,lpar,acidps
    586       REAL  x1, satpacid
    587       REAL , DIMENSION(2):: act
    588       REAL  CONCM
    589       REAL  NH2SO4
    590       REAL  H2OCOND, H2SO4COND
    591 
    592 
    593       CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS!
    594 
    595         NH2SO4=SAt*CONCM
    596 
    597       pstand=1.01325E+5 !Pa  1 atm pressure
    598 
    599         x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
    600 
    601         CALL zeleznik(x1,T,act)
    602 
    603 !pure acid satur vapor pressure
    604         lpar= -11.695+DLOG(pstand) ! Zeleznik
    605         acidps=1/360.15-1.0/T+0.38/545.
    606      & *(1.0+DLOG(360.15/T)-360.15/T)
    607         acidps = 10156.0*acidps +lpar
    608         acidps = DEXP(acidps)    !Pa
    609 
    610 !acid sat.vap.PP over mixture (flat surface):
    611         satpacid=act(2)*acidps ! Pa
    612 
    613 !       Conversion from Pa to N.D #/m3
    614         DND2=satpacid/(1.3806488E-23*T)
    615 
    616 !       H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
    617         IF (NH2SO4.GT.DND2) THEN
    618         H2SO4COND=NH2SO4-DND2
    619 !       calcul de H2O cond correspondant a H2SO4 cond
    620         H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
    621 
    622 !       Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30
    623         ELSE
    624         H2OCOND=1.0E-30*CONCM
    625         END IF
    626 
    627 !*****************************************************
    628 !     ATTENTION: Ici on ne prends pas en compte
    629 !                si H2O en defaut!
    630 !                On veut la situation theorique
    631 !                a l'equilibre
    632 !*****************************************************
    633 !       Test si H2O en defaut H2Ocond>H2O dispo
    634 !       IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN
    635 
    636 !       On peut alors condenser tout le H2O dispo
    637 !       H2OCOND=NH2O
    638 !       On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
    639 !       H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
    640 
    641 !      END IF
    642 
    643 !     Calcul de H2O condensŽe VMR
    644       WVCOND=H2OCOND/CONCM
    645 
    646       END FUNCTION WVCOND
    647 
    648 !*****************************************************************************
    649 !*     REAL FUNCTION IRFRMWV()
    650       REAL      FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,
    651      + WVTOT,SATOT,NROOT)
    652 !*****************************************************************************
    653 !* Iterative Root Finder Ridder's Method for Water Vapor calculus
    654 !* From Numerical Recipes
    655 !* Adapted for VenusGCM A. Stolzenbach 07/2014
    656 !*
    657 !* Les iterations sur [X1,X2] sont [WV1,WV2]
    658 !* la variable X est WV
    659 !* IRFRMWV sort en OUTPUT : WSALOC pour ITERWV=0 (ou WVEQ=0)
    660 
    661       IMPLICIT NONE
    662 
    663       REAL, INTENT(IN) ::    X1, X2
    664       REAL, INTENT(IN) ::    XACC
    665       INTEGER, INTENT(IN) :: MAXIT,NROOT
    666 
    667 !     LOCAL VARIABLES
    668       REAL :: XL, XH, XM, XNEW, X
    669       REAL :: WSALOC, WVEQ, WVLIQ
    670       REAL :: FL, FH, FM, FNEW
    671       REAL :: ANS, S, FSIGN
    672       INTEGER i
    673 
    674 !     External variables needed:
    675       REAL, INTENT(IN) :: TAIR,PAIR
    676       REAL, INTENT(IN) :: WVTOT,SATOT
    677       REAL, INTENT(IN) :: RADIUS
    678 
    679 
    680 !     Initialisation
    681       X=X1
    682       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
    683       FL=WVEQ
    684       X=X2
    685       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
    686       FH=WVEQ
    687 
    688 !     Test Bracketed values
    689       IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.
    690      &   ((FL.GT.0.).AND.(FH.LT.0.)))
    691      &  THEN
    692          XL=X1
    693          XH=X2
    694          ANS=-9.99e99
    695 
    696          DO i=1, MAXIT
    697             XM=0.5*(XL+XH)
    698             CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    699      &       TAIR,PAIR,RADIUS)
    700             FM=WVEQ
    701             S=SQRT(FM*FM-FL*FH)
    702 
    703             IF (S.EQ.0.0) THEN
    704                IRFRMWV=WSALOC
    705                RETURN
    706             ENDIF
    707 
    708             IF (FL.GT.FH) THEN
    709                FSIGN=1.0
    710             ELSE
    711                FSIGN=-1.0
    712             ENDIF
    713 
    714             XNEW=XM+(XM-XL)*(FSIGN*FM/S)
    715 
    716             IF (ABS(XNEW-ANS).LE.XACC) THEN
    717                IRFRMWV=WSALOC
    718                RETURN
    719             ENDIF
    720 
    721             ANS=XNEW
    722             CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    723      &       TAIR,PAIR,RADIUS)
    724             FNEW=WVEQ
    725 
    726             IF (FNEW.EQ.0.0) THEN
    727                IRFRMWV=WSALOC
    728                RETURN
    729             ENDIF
    730 
    731             IF (SIGN(FM, FNEW).NE.FM) THEN
    732                XL=XM
    733                FL=FM
    734                XH=ANS
    735                FH=FNEW
    736                ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
    737                   XH=ANS
    738                   FH=FNEW
    739                ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
    740                   XL=ANS
    741                   FL=FNEW
    742                ELSE
    743                   PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
    744                   PRINT*,'you shall not PAAAAAASS'
    745                   STOP
    746             ENDIF
    747          ENDDO
    748          PRINT*,'Paaaaas bien MAXIT atteint'
    749          PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
    750          PRINT*,'you shall not PAAAAAASS'
    751          XL=X1
    752          XH=X2
    753          ANS=-9.99e99
    754 
    755          DO i=1, MAXIT
    756             XM=0.5*(XL+XH)
    757             CALL ITERWV(XM,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    758      &       TAIR,PAIR,RADIUS)
    759             FM=WVEQ
    760             S=SQRT(FM*FM-FL*FH)
    761             IF (FL.GT.FH) THEN
    762                FSIGN=1.0
    763             ELSE
    764                FSIGN=-1.0
    765             ENDIF
    766 
    767             XNEW=XM+(XM-XL)*(FSIGN*FM/S)
    768 
    769             ANS=XNEW
    770             CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    771      &       TAIR,PAIR,RADIUS)
    772             FNEW=WVEQ
    773             PRINT*,'WVliq',WVLIQ,'WVtot',WVTOT,'WVeq',WVEQ
    774             PRINT*,'WSA',WSALOC,'SAtot',SATOT
    775             PRINT*,'T',TAIR,'P',PAIR
    776 
    777             IF (SIGN(FM, FNEW).NE.FM) THEN
    778                XL=XM
    779                FL=FM
    780                XH=ANS
    781                FH=FNEW
    782                ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
    783                   XH=ANS
    784                   FH=FNEW
    785                ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
    786                   XL=ANS
    787                   FL=FNEW
    788                ELSE
    789                   PRINT*,'PROBLEM IRFRMWV dans new_cloud_venus'
    790                   PRINT*,'you shall not PAAAAAASS TWIIICE???'
    791                   STOP
    792             ENDIF
    793          ENDDO
    794          STOP
    795       ELSE
    796          PRINT*,'IRFRMWV must be bracketed'
    797          PRINT*,'NROOT de BRACWV', NROOT
    798          IF (ABS(FL).LT.XACC) THEN
    799          PRINT*,'IRFRMWV FL == 0',FL
    800          PRINT*,'X1',X1,'X2',X2,'FH',FH
    801            CALL ITERWV(X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    802      &      TAIR,PAIR,RADIUS)
    803            IRFRMWV=WSALOC
    804          RETURN
    805          ENDIF
    806          IF (ABS(FH).LT.XACC) THEN
    807          PRINT*,'IRFRMWV FH == 0',FH
    808          PRINT*,'X1',X1,'X2',X2,'FL',FL
    809            CALL ITERWV(X2,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    810      &      TAIR,PAIR,RADIUS)
    811            IRFRMWV=WSALOC
    812          RETURN
    813          ENDIF
    814          IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN
    815            PRINT*,'STOP dans IRFRMWV avec rien == 0'
    816            PRINT*,'X1',X1,'X2',X2
    817            PRINT*,'Fcalc',FL,FH
    818            PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS
    819            STOP
    820          ENDIF
    821          IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN
    822            PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC'
    823            PRINT*,FL,FH
    824            STOP
    825          ENDIF
    826 
    827 
    828       ENDIF
    829 !  FIN Test Bracketed values
    830 
    831       END FUNCTION IRFRMWV
    832 
    833 !*****************************************************************************
    834 !*     REAL FUNCTION IRFRMSA()
    835       REAL      FUNCTION IRFRMSA(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,LPPWV,
    836      +               NB)
    837 !*****************************************************************************
    838 !* Iterative Root Finder Ridder's Method for Sulfuric Acid calculus
    839 !* From Numerical Recipes
    840 !* Adapted for VenusGCM A. Stolzenbach 07/2014
    841 !*
    842 !* Les iterations sur [X1,X2] sont [WSA1,WSA2]
    843 !* la variable X est WSA
    844 !* IRFRMSA sort en OUTPUT : WSA pour KEEQ=0
    845 
    846       IMPLICIT NONE
    847 
    848       REAL, INTENT(IN) ::    X1, X2
    849       REAL, INTENT(IN) ::    XACC
    850       INTEGER, INTENT(IN) :: MAXIT, NB
    851 
    852 !     LOCAL VARIABLES
    853       REAL XL, XH, XM, XNEW
    854       REAL Fl, FH, FM, FNEW
    855       REAL ANS, S, FSIGN
    856       INTEGER i
    857 
    858 !     External variables needed:
    859       REAL, INTENT(IN) :: TAIR,PAIR
    860       REAL, INTENT(IN) :: LPPWV
    861       REAL, INTENT(IN) :: RADIUS
    862 
    863 !     External functions needed:
    864       REAL KEEQ
    865 
    866 
    867 
    868 !     Initialisation
    869       FL=KEEQ(RADIUS,TAIR,X1,LPPWV)
    870       FH=KEEQ(RADIUS,TAIR,X2,LPPWV)
    871 
    872 !     Test Bracketed values
    873       IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.((FL.GT.0.).AND.(FH.LT.0.)))
    874      &  THEN
    875          XL=X1
    876          XH=X2
    877          ANS=-9.99e99
    878 
    879          DO i=1, MAXIT
    880             XM=0.5*(XL+XH)
    881             FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
    882             S=SQRT(FM*FM-FL*FH)
    883 
    884             IF (S.EQ.0.0) THEN
    885                IRFRMSA=ANS
    886                RETURN
    887             ENDIF
    888 
    889             IF (FL.GT.FH) THEN
    890                FSIGN=1.0
    891             ELSE
    892                FSIGN=-1.0
    893             ENDIF
    894 
    895             XNEW=XM+(XM-XL)*(FSIGN*FM/S)
    896 
    897             IF (ABS(XNEW-ANS).LE.XACC) THEN
    898                IRFRMSA=ANS
    899                RETURN
    900             ENDIF
    901 
    902             ANS=XNEW
    903             FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
    904 
    905             IF (FNEW.EQ.0.0) THEN
    906                IRFRMSA=ANS
    907                RETURN
    908             ENDIF
    909 
    910             IF (SIGN(FM, FNEW).NE.FM) THEN
    911                XL=XM
    912                FL=FM
    913                XH=ANS
    914                FH=FNEW
    915                ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
    916                   XH=ANS
    917                   FH=FNEW
    918                ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
    919                   XL=ANS
    920                   FL=FNEW
    921                ELSE
    922                   PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
    923                   PRINT*,'you shall not PAAAAAASS'
    924                   STOP
    925             ENDIF
    926          ENDDO
    927          PRINT*,'Paaaaas bien MAXIT atteint'
    928          PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
    929          PRINT*,'you shall not PAAAAAASS'
    930          XL=X1
    931          XH=X2
    932          PRINT*,'Borne XL',XL,'XH',XH
    933          ANS=-9.99e99
    934 
    935          DO i=1, MAXIT
    936             XM=0.5*(XL+XH)
    937             FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
    938             S=SQRT(FM*FM-FL*FH)
    939 
    940             IF (FL.GT.FH) THEN
    941                FSIGN=1.0
    942             ELSE
    943                FSIGN=-1.0
    944             ENDIF
    945 
    946             XNEW=XM+(XM-XL)*(FSIGN*FM/S)
    947 
    948             ANS=XNEW
    949             FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
    950             PRINT*,'KEEQ result',FNEW,'T',TAIR,'R',RADIUS
    951             IF (SIGN(FM, FNEW).NE.FM) THEN
    952                XL=XM
    953                FL=FM
    954                XH=ANS
    955                FH=FNEW
    956                ELSEIF (SIGN(FL, FNEW).NE.FL) THEN
    957                   XH=ANS
    958                   FH=FNEW
    959                ELSEIF (SIGN(FH, FNEW).NE.FH) THEN
    960                   XL=ANS
    961                   FL=FNEW
    962                ELSE
    963                   PRINT*,'PROBLEM IRFRMSA dans new_cloud_venus'
    964                   PRINT*,'you shall not PAAAAAASS'
    965                   STOP
    966             ENDIF
    967          ENDDO
    968          STOP
    969       ELSE
    970          PRINT*,'IRFRMSA must be bracketed'
    971          IF (FL.EQ.0.0) THEN
    972            PRINT*,'IRFRMSA FL == 0',Fl
    973            IRFRMSA=X1
    974            RETURN
    975          ENDIF
    976          IF (FH.EQ.0.0) THEN
    977            PRINT*,'IRFRMSA FH == 0',FH
    978            IRFRMSA=X2
    979            RETURN
    980          ENDIF
    981          IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN
    982            PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH
    983            PRINT*,'X1',X1,'X2',X2
    984            PRINT*,'Kind F', KIND(FL), KIND(FH)
    985            PRINT*,'Kind X', KIND(X1), KIND(X2)
    986            PRINT*,'Logical: ',(SIGN(FL,FH).NE.FL)
    987            PRINT*,'Logical: ',(SIGN(FH,FL).NE.FH)
    988            PRINT*,'nb root BRACWSA',NB
    989            STOP
    990          ENDIF
    991 
    992       ENDIF
    993 !  FIN Test Bracketed values
    994 
    995       END function IRFRMSA
    996 
    997 !*****************************************************************************
    998 !*     REAL FUNCTION KEEQ()
    999       REAL      FUNCTION KEEQ(RADIUS,TAIR,WSA,LPPWV)
    1000 !*****************************************************************************
    1001 !* Kelvin Equation EQuality
    1002 !* ln(PPWV_eq) - (2Mh2o sigma)/(R T r rho) - ln(ph2osa) = 0
    1003 !*
    1004 
    1005       IMPLICIT NONE
    1006 
    1007       REAL, INTENT(IN) :: RADIUS,TAIR,WSA,LPPWV
    1008 
    1009 !     Physical constants:
    1010       REAL   MH2O
    1011       REAL   RGAS
    1012       PARAMETER(
    1013 !       Molar weight of water (kg/mole)
    1014      +          MH2O=18.0153d-3,
    1015 !       Universal gas constant (J/(mole K))
    1016      +          RGAS=8.314462175d0)
    1017 !
    1018 !     External functions needed:
    1019       REAL   PWVSAS_GV,SIGMADROPLET,RHODROPLET
    1020 !     PWVSAS_GV:      Natural logaritm of water vapor pressure over
    1021 !                  sulfuric acid solution
    1022 !     SIGMADROPLET:       Surface tension of sulfuric acid solution
    1023 !     RHODROPLET:       Density of sulfuric acid solution
    1024 !
    1025 !     Auxiliary local variables:
    1026       REAL   C1
    1027 
    1028       PARAMETER(
    1029      +        C1=2.0d0*MH2O/RGAS)
    1030 
    1031 
    1032       KEEQ=LPPWV-C1*SIGMADROPLET(WSA,TAIR)/
    1033      &     (TAIR*RADIUS*RHODROPLET(WSA,TAIR))-
    1034      &     PWVSAS_GV(TAIR,WSA)
    1035 
    1036       END FUNCTION KEEQ
    1037 
    1038 *****************************************************************************
    1039 *     REAL FUNCTION PWVSAS_GV(TAIR,WSA)                                       
    1040       REAL FUNCTION PWVSAS_GV(TAIR,WSA)
    1041 *****************************************************************************
    1042 *
    1043 *     Natural logaritm of saturated water vapor pressure over plane
    1044 *     sulfuric acid solution.
    1045 *
    1046 *     Source: J.I.Gmitro & T.Vermeulen: A.I.Ch.E.J.  10,740,1964.
    1047 *             W.F.Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
    1048 *
    1049 *     The formula of Gmitro & Vermeulen for saturation pressure
    1050 *     is used:
    1051 *                 ln(p) = A ln(298/T) + B/T + C + DT
    1052 *     with values of A,B,C and D given by Gmitro & Vermeulen,
    1053 *     and calculated from partial molal properties given by Giauque et al.
    1054 *     
    1055 *     
    1056 *
    1057 *     Input:  TAIR: Temperature (K)
    1058 *             WSA:  Weight fraction of H2SO4  [0;1]
    1059 *     Output: Natural logaritm of water vapor pressure
    1060 *             over sulfuric acid solution   ( ln(Pa) )
    1061 *
    1062 *
    1063 *     External functions needed for calculation of partial molal
    1064 *     properties of pure components at 25 ! as function of W.
    1065       IMPLICIT NONE
    1066 
    1067       REAL :: CPH2O,ALH2O,FFH2O,LH2O
    1068 *     CPH2O:  Partial molal heat capacity of sulfuric acid solution.
    1069 *     ALH2O:  Temparature derivative of CPH2O
    1070 *     FFH2O:  Partial molal free energy of sulfuric acid solution.
    1071 *     LH2O:   Partial molal enthalpy of sulfuric acid
    1072 *
    1073 !
    1074 !
    1075       REAL, INTENT(IN) :: TAIR,WSA
    1076       REAL :: ADOT,BDOT,CDOT,DDOT
    1077       REAL :: RGAS,MMHGPA
    1078       REAL :: K1,K2
    1079       REAL :: A,B,C,D,CP,L,F,ALFA
    1080 !     Physical constants given by Gmitro & Vermeulen:
    1081       PARAMETER(
    1082      +        ADOT=-3.67340,
    1083      +        BDOT=-4143.5,
    1084      +        CDOT=10.24353,
    1085      +        DDOT=0.618943d-3)
    1086       PARAMETER(
    1087 !     Gas constant (cal/(deg mole)):
    1088      +        RGAS=1.98726,
    1089 !     Natural logarith of conversion factor between atm. and Pa:     
    1090      +     MMHGPA=11.52608845,
    1091      +     K1=298.15,
    1092      +     K2=K1*K1/2.0)
    1093 !
    1094       INTEGER :: KLO,KHI,K,I,NPOINT
    1095       PARAMETER(NPOINT=110)
    1096       REAL, DIMENSION(NPOINT) :: WTAB(NPOINT)
    1097       DATA (WTAB(I),I=1,NPOINT)/
    1098      +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
    1099      +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
    1100      +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
    1101      +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
    1102      +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
    1103      +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
    1104      +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
    1105      +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
    1106      +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
    1107      +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
    1108      +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
    1109      +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
    1110      +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
    1111      +0.99908,0.99927,0.99945,0.99963,0.99982,1.0000/
    1112 !
    1113       KLO=1
    1114       KHI=NPOINT
    1115  1    IF(KHI-KLO.GT.1) THEN
    1116           K=(KHI+KLO)/2
    1117           IF(WTAB(K).GT.MAX(WTAB(1),WSA)) THEN
    1118               KHI=K
    1119           ELSE
    1120               KLO=K
    1121           ENDIF
    1122           GOTO 1
    1123       ENDIF
    1124 !
    1125       CP=CPH2O(WSA,KHI,KLO)
    1126       F=-FFH2O(WSA,KHI,KLO)
    1127       L=-LH2O(WSA,KHI,KLO)
    1128       ALFA=ALH2O(WSA,KHI,KLO)
    1129 !
    1130       A=ADOT+(CP-K1*ALFA)/RGAS
    1131       B=BDOT+(L-K1*CP+K2*ALFA)/RGAS
    1132       C=CDOT+(CP+(F-L)/K1)/RGAS
    1133       D=DDOT-ALFA/(2.0d0*RGAS)
    1134 !
    1135 !     WRITE(*,*) 'TAIR= ',TAIR,'  WSA= ',WSA
    1136 !     WRITE(*,*) 'CPH2O(WSA)= ',CP
    1137 !     WRITE(*,*) 'ALFAH2O(WSA)= ',ALFA
    1138 !     WRITE(*,*) 'FFH2O(WSA)= ',F
    1139 !     WRITE(*,*) 'LH2O(WSA)= ',L
    1140 !
    1141       PWVSAS_GV=A*DLOG(K1/TAIR)+B/TAIR+C+D*TAIR+MMHGPA
    1142      
    1143       END FUNCTION PWVSAS_GV
    1144 *******************************************************************************
    1145 *     REAL FUNCTION CPH2O(W)
    1146       REAL FUNCTION CPH2O(W,khi_in,klo_in)
    1147 *******************************************************************************
    1148 *
    1149 *     Relative partial molal heat capacity of water (cal/(deg mole) in
    1150 *     sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
    1151 *     calculated by cubic spline fitting.
    1152 *
    1153 *     Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
    1154 *
    1155       IMPLICIT NONE
    1156 
    1157       INTEGER :: NPOINT,I
    1158       PARAMETER(NPOINT=109)
    1159       REAL, DIMENSION(NPOINT) :: WTAB(NPOINT),CPHTAB(NPOINT),
    1160      +              Y2(NPOINT),YWORK(NPOINT)
    1161       REAL, INTENT(IN):: W
    1162       INTEGER, INTENT(IN):: khi_in,klo_in
    1163       INTEGER :: khi,klo
    1164       REAL :: CPH
    1165       LOGICAL :: FIRST
    1166       DATA (WTAB(I),I=1,NPOINT)/
    1167      +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
    1168      +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
    1169      +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
    1170      +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
    1171      +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
    1172      +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
    1173      +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
    1174      +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
    1175      +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
    1176      +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
    1177      +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
    1178      +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
    1179      +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
    1180      +0.99908,0.99927,0.99945,0.99963,0.99982/
    1181       DATA (CPHTAB(I),I=1,NPOINT)/
    1182      + 17.996, 17.896, 17.875, 17.858, 17.840, 17.820, 17.800, 17.791,
    1183      + 17.783, 17.777, 17.771, 17.769, 17.806, 17.891, 18.057, 18.248,
    1184      + 18.429, 18.567, 18.613, 18.640, 18.660, 18.660, 18.642, 18.592,
    1185      + 18.544, 18.468, 18.348, 18.187, 17.995, 17.782, 17.562, 17.352,
    1186      + 17.162, 16.993, 16.829, 16.657, 16.581, 16.497, 16.405, 16.302,
    1187      + 16.186, 16.053, 15.901, 15.730, 15.540, 15.329, 15.101, 14.853,
    1188      + 14.586, 14.296, 13.980, 13.638, 13.274, 12.896, 12.507, 12.111,
    1189      + 11.911, 11.711, 11.514, 11.320, 11.130, 10.940, 10.760, 10.570,
    1190      + 10.390, 10.200, 10.000, 9.8400, 9.7600, 9.7900, 9.9500, 10.310,
    1191      + 10.950, 11.960, 13.370, 15.060, 16.860, 18.550, 20.000, 21.170,
    1192      + 22.030, 22.570, 22.800, 22.750, 22.420, 21.850, 21.120, 20.280,
    1193      + 19.360, 18.350, 17.220, 15.940, 14.490, 12.840, 10.800, 9.8000,
    1194      + 7.8000, 3.8000,0.20000,-5.4000,-7.0000,-8.8000,-10.900,-13.500,
    1195      +-17.000,-22.000,-29.000,-40.000,-59.000/
    1196       DATA FIRST/.TRUE./
    1197       SAVE FIRST,WTAB,CPHTAB,Y2
    1198 !
    1199       IF(FIRST) THEN
    1200           FIRST=.FALSE.
    1201           CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2)
    1202       ENDIF
    1203 
    1204       if(khi_in.GT.NPOINT) then
    1205          khi=NPOINT
    1206          klo=NPOINT-1
    1207       else
    1208          khi=khi_in
    1209          klo=klo_in
    1210       endif
    1211 
    1212       CALL SPLINT(WTAB(khi),WTAB(klo),CPHTAB(khi),CPHTAB(klo),
    1213      .              Y2(khi),Y2(klo),W,CPH)
    1214       CPH2O=CPH
    1215      
    1216       END FUNCTION CPH2O
    1217 !
    1218 *******************************************************************************
    1219       REAL FUNCTION FFH2O(W,khi,klo)
    1220 *     REAL FUNCTION FFH2O(W)
    1221 *******************************************************************************
    1222 *
    1223 *     Relative partial molal free energy water (cal/mole) in
    1224 *     sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
    1225 *     calculated by cubic spline fitting.
    1226 *
    1227 *     Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
    1228 *
    1229       IMPLICIT NONE
    1230 
    1231       INTEGER :: NPOINT,I
    1232       PARAMETER(NPOINT=110)
    1233       REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
    1234       REAL, INTENT(IN) :: W
    1235       INTEGER, INTENT(IN):: khi,klo
    1236       REAL :: FF
    1237       LOGICAL :: FIRST
    1238       DATA (WTAB(I),I=1,NPOINT)/
    1239      +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
    1240      +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
    1241      +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
    1242      +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
    1243      +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
    1244      +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
    1245      +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
    1246      +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
    1247      +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
    1248      +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
    1249      +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
    1250      +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
    1251      +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
    1252      +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
    1253       DATA (FFTAB(I),I=1,NPOINT)/
    1254      +0.00000, 22.840, 25.810, 29.250, 33.790, 39.970, 48.690, 54.560,
    1255      + 61.990, 71.790, 85.040, 103.70, 130.70, 145.20, 163.00, 184.50,
    1256      + 211.50, 245.60, 266.40, 290.10, 317.40, 349.00, 385.60, 428.40,
    1257      + 452.50, 478.80, 507.50, 538.80, 573.30, 611.60, 653.70, 700.50,
    1258      + 752.60, 810.60, 875.60, 948.60, 980.60, 1014.3, 1049.7, 1087.1,
    1259      + 1126.7, 1168.7, 1213.5, 1261.2, 1312.0, 1366.2, 1424.3, 1486.0,
    1260      + 1551.8, 1622.3, 1697.8, 1778.5, 1864.9, 1956.8, 2055.8, 2162.0,
    1261      + 2218.0, 2276.0, 2337.0, 2400.0, 2466.0, 2535.0, 2607.0, 2682.0,
    1262      + 2760.0, 2842.0, 2928.0, 3018.0, 3111.0, 3209.0, 3311.0, 3417.0,
    1263      + 3527.0, 3640.0, 3757.0, 3878.0, 4002.0, 4130.0, 4262.0, 4397.0,
    1264      + 4535.0, 4678.0, 4824.0, 4973.0, 5128.0, 5287.0, 5454.0, 5630.0,
    1265      + 5820.0, 6031.0, 6268.0, 6541.0, 6873.0, 7318.0, 8054.0, 8284.0,
    1266      + 8579.0, 8997.0, 9295.0, 9720.0, 9831.0, 9954.0, 10092., 10248.,
    1267      + 10423., 10618., 10838., 11099., 11460., 12014./
    1268       DATA FIRST/.TRUE./
    1269       SAVE FIRST,WTAB,FFTAB,Y2
    1270 !
    1271       IF(FIRST) THEN
    1272           FIRST=.FALSE.
    1273           CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2)
    1274       ENDIF
    1275 
    1276       CALL SPLINT(WTAB(khi),WTAB(klo),FFTAB(khi),FFTAB(klo),
    1277      .              Y2(khi),Y2(klo),W,FF)
    1278       FFH2O=FF
    1279      
    1280       END FUNCTION FFH2O
    1281 !
    1282 *******************************************************************************
    1283       REAL FUNCTION LH2O(W,khi,klo)
    1284 *     REAL FUNCTION LH2O(W)
    1285 *******************************************************************************
    1286 *
    1287 *     Relative partial molal heat content of water (cal/mole) in
    1288 *     sulfuric acid solution, as a function of H2SO4 weight fraction [0;1],
    1289 *     calculated by cubic spline fitting.
    1290 *
    1291 *     Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
    1292 *
    1293       IMPLICIT NONE
    1294 
    1295       INTEGER :: NPOINT,I
    1296       PARAMETER(NPOINT=110)
    1297       REAL, DIMENSION(NPOINT) ::  WTAB,LTAB,Y2,YWORK
    1298       REAL, INTENT(IN) :: W
    1299       INTEGER, INTENT(IN):: khi,klo
    1300       REAL :: L
    1301       LOGICAL :: FIRST
    1302       DATA (WTAB(I),I=1,NPOINT)/
    1303      +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
    1304      +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
    1305      +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
    1306      +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
    1307      +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
    1308      +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
    1309      +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
    1310      +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
    1311      +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
    1312      +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
    1313      +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
    1314      +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
    1315      +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
    1316      +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
    1317       DATA (LTAB(I),I=1,NPOINT)/
    1318      +0.00000, 5.2900, 6.1000, 7.1800, 8.7800, 11.210, 15.290, 18.680,
    1319      + 23.700, 31.180, 42.500, 59.900, 89.200, 106.70, 128.60, 156.00,
    1320      + 190.40, 233.80, 260.10, 290.00, 324.00, 362.50, 406.50, 456.10,
    1321      + 483.20, 512.40, 543.60, 577.40, 613.80, 653.50, 696.70, 744.50,
    1322      + 797.20, 855.80, 921.70, 995.70, 1028.1, 1062.3, 1098.3, 1136.4,
    1323      + 1176.7, 1219.3, 1264.7, 1313.0, 1364.3, 1418.9, 1477.3, 1539.9,
    1324      + 1607.2, 1679.7, 1757.9, 1842.7, 1934.8, 2035.4, 2145.5, 2267.0,
    1325      + 2332.0, 2401.0, 2473.0, 2550.0, 2631.0, 2716.0, 2807.0, 2904.0,
    1326      + 3007.0, 3118.0, 3238.0, 3367.0, 3507.0, 3657.0, 3821.0, 3997.0,
    1327      + 4186.0, 4387.0, 4599.0, 4819.0, 5039.0, 5258.0, 5476.0, 5694.0,
    1328      + 5906.0, 6103.0, 6275.0, 6434.0, 6592.0, 6743.0, 6880.0, 7008.0,
    1329      + 7133.0, 7255.0, 7376.0, 7497.0, 7618.0, 7739.0, 7855.0, 7876.0,
    1330      + 7905.0, 7985.0, 8110.0, 8415.0, 8515.0, 8655.0, 8835.0, 9125.0,
    1331      + 9575.0, 10325., 11575., 13500., 15200., 16125./
    1332       DATA FIRST/.TRUE./
    1333       SAVE FIRST,WTAB,LTAB,Y2
    1334 !
    1335       IF(FIRST) THEN
    1336           FIRST=.FALSE.
    1337           CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2)
    1338       ENDIF
    1339 
    1340       CALL SPLINT(WTAB(khi),WTAB(klo),LTAB(khi),LTAB(klo),
    1341      .              Y2(khi),Y2(klo),W,L)
    1342       LH2O=L
    1343      
    1344       END FUNCTION LH2O
    1345 *******************************************************************************
    1346       REAL FUNCTION ALH2O(W,khi_in,klo_in)
    1347 *     REAL FUNCTION ALH2O(W)
    1348 *******************************************************************************
    1349 *
    1350 *     Relative partial molal temperature derivative of heat capacity (water)
    1351 *     in sulfuric acid solution, (cal/deg**2), calculated by
    1352 *     cubic spline fitting.
    1353 *
    1354 *     Source: Giauque et al.: J. Amer. Chem. Soc. 82,62,1960.
    1355 *
    1356       IMPLICIT NONE
    1357 
    1358       INTEGER :: NPOINT,I
    1359       PARAMETER(NPOINT=96)
    1360       REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
    1361       REAL, INTENT(IN) :: W
    1362       INTEGER, INTENT(IN):: khi_in,klo_in
    1363       INTEGER :: khi,klo
    1364       REAL :: A
    1365       LOGICAL :: FIRST
    1366       DATA (WTAB(I),I=1,NPOINT)/
    1367      +0.29517,0.31209,
    1368      +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
    1369      +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
    1370      +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
    1371      +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
    1372      +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
    1373      +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
    1374      +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
    1375      +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
    1376      +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
    1377      +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
    1378      +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
    1379      +0.99908,0.99927,0.99945,0.99963,0.99982, 1.0000/
    1380       DATA (ATAB(I),I=1,NPOINT)/
    1381      + 0.0190, 0.0182, 0.0180, 0.0177, 0.0174, 0.0169, 0.0167, 0.0164,
    1382      + 0.0172, 0.0212, 0.0239, 0.0264, 0.0276, 0.0273, 0.0259, 0.0238,
    1383      + 0.0213, 0.0190, 0.0170, 0.0155, 0.0143, 0.0133, 0.0129, 0.0124,
    1384      + 0.0120, 0.0114, 0.0106, 0.0097, 0.0084, 0.0067, 0.0047, 0.0024,
    1385      +-0.0002,-0.0031,-0.0063,-0.0097,-0.0136,-0.0178,-0.0221,-0.0263,
    1386      +-0.0303,-0.0340,-0.0352,-0.0360,-0.0362,-0.0356,-0.0343,-0.0321,
    1387      +-0.0290,-0.0251,-0.0201,-0.0137,-0.0058, 0.0033, 0.0136, 0.0254,
    1388      + 0.0388, 0.0550, 0.0738, 0.0962, 0.1198, 0.1300, 0.1208, 0.0790,
    1389      + 0.0348, 0.0058,-0.0102,-0.0211,-0.0292,-0.0350,-0.0390,-0.0418,
    1390      +-0.0432,-0.0436,-0.0429,-0.0411,-0.0384,-0.0346,-0.0292,-0.0220,
    1391      +-0.0130,-0.0110,-0.0080,-0.0060,-0.0040,-0.0030,-0.0030,-0.0020,
    1392      +-0.0020,-0.0020,-0.0020,-0.0010,-0.0010, 0.0000, 0.0000, 0.0000/
    1393       DATA FIRST/.TRUE./
    1394       SAVE FIRST,WTAB,ATAB,Y2
    1395 !
    1396       IF(FIRST) THEN
    1397           FIRST=.FALSE.
    1398           CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2)
    1399       ENDIF
    1400 
    1401       if(klo_in.LT.15) then
    1402          khi=2
    1403          klo=1
    1404       else
    1405          khi=khi_in-14
    1406          klo=klo_in-14
    1407       endif
    1408 
    1409       CALL SPLINT(WTAB(khi),WTAB(klo),ATAB(khi),ATAB(klo),
    1410      .              Y2(khi),Y2(klo),W,A)
    1411       ALH2O=A
    1412      
    1413       END FUNCTION ALH2O
    1414 !******************************************************************************
    1415       SUBROUTINE SPLINE(X,Y,N,WORK,Y2)
    1416 !******************************************************************************
    1417 !     Routine to calculate 2.nd derivatives of tabulated function
    1418 !     Y(i)=Y(Xi), to be used for cubic spline calculation.
    1419 !
    1420       IMPLICIT NONE
    1421 
    1422       INTEGER, INTENT(IN) :: N
    1423       INTEGER :: I
    1424       REAL, DIMENSION(N), INTENT(IN) :: X,Y
    1425       REAL, DIMENSION(N), INTENT(OUT) :: Y2,WORK
    1426       REAL   SIG,P,QN,UN,YP1,YPN
    1427 
    1428 !AM Venus: Let's check the values
    1429 !      write(*,*) 'In spline, N ', N
    1430 
    1431       YP1=(Y(2)-Y(1))/(X(2)-X(1))
    1432       YPN=(Y(N)-Y(N-1))/(X(N)-X(N-1))
    1433       IF(YP1.GT.99.0E+30) THEN
    1434           Y2(1)=0.0
    1435           WORK(1)=0.0
    1436       ELSE
    1437           Y2(1)=-0.5d0
    1438           WORK(1)=(3.0d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
    1439       ENDIF
    1440       DO I=2,N-1
    1441 !         write(*,*) 'In spline, I ', I
    1442           SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
    1443           P=SIG*Y2(I-1)+2.0d0
    1444           Y2(I)=(SIG-1.0d0)/P
    1445           WORK(I)=(6.0d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
    1446      +             /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*WORK(I-1))/P
    1447       ENDDO
    1448       IF(YPN.GT.99.0E+30) THEN
    1449           QN=0.0
    1450           UN=0.0
    1451       ELSE
    1452           QN=0.5d0
    1453           UN=(3.0d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
    1454       ENDIF
    1455       Y2(N)=(UN-QN*WORK(N-1))/(QN*Y2(N-1)+1.0d0)
    1456       DO I=N-1,1,-1
    1457 !         write(*,*) 'In spline, I ', I
    1458           Y2(I)=Y2(I)*Y2(I+1)+WORK(I)
    1459       ENDDO
    1460 !
    1461       END SUBROUTINE SPLINE
    1462 
    1463 !******************************************************************************
    1464       SUBROUTINE SPLINT(XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo,X,Y)
    1465 !******************************************************************************
    1466 !     Cubic spline calculation
    1467 
    1468       IMPLICIT NONE
    1469 
    1470       REAL, INTENT(IN) :: XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo
    1471       REAL, INTENT(IN) :: X
    1472       REAL, INTENT(OUT) :: Y
    1473       REAL :: H,A,B
    1474 !
    1475       H=XAhi-XAlo
    1476       A=(XAhi-X)/H
    1477       B=(X-XAlo)/H
    1478       Y=A*YAlo+B*YAhi+
    1479      +        ((A**3-A)*Y2Alo+(B**3-B)*Y2Ahi)*(H**2)/6.0d0
    1480 !
    1481 
    1482       END SUBROUTINE SPLINT
    1483296!******************************************************************
    1484297      SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
    … …  
    1593406      END SUBROUTINE CALCM_SAT
    1594407
    1595        SUBROUTINE Zeleznik(x,T,act)
    1596 
    1597   !+++++++++++++++++++++++++++++++++++++++++++++++++++
    1598   !     Water and sulfuric acid activities in liquid
    1599   !     aqueous solutions.
    1600   !     Frank J. Zeleznik, Thermodynnamic properties
    1601   !     of the aqueous sulfuric acid system to 220K-350K,
    1602   !     mole fraction 0,...,1
    1603   !     J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991
    1604   !+++++++++++++++++++++++++++++++++++++++++++++++++++
    1605 
    1606          IMPLICIT NONE
    1607 
    1608          REAL, INTENT(IN) :: x,T
    1609          REAL :: activitya, activityw
    1610          REAL, INTENT(OUT), DIMENSION(2):: act
    1611 !         REAL x,T, activitya, activityw
    1612 !         REAL, DIMENSION(2):: act
    1613 
    1614 
    1615 !         write(*,*) 'x, T ', x, T
    1616 
    1617          act(2)=activitya(x,T)
    1618          act(1)=activityw(x,T)
    1619 
    1620 !         write(*,*) 'act ', act
    1621 
    1622        END SUBROUTINE Zeleznik
    1623 
    1624 !start of functions related to zeleznik activities
    1625 
    1626       REAL FUNCTION m111(T)
    1627 
    1628        REAL, INTENT(IN) :: T
    1629        m111=-23.524503387D0
    1630      &    +0.0406889449841D0*T
    1631      &    -0.151369362907D-4*T**2+2961.44445015D0/T
    1632      &    +0.492476973663D0*dlog(T)
    1633        END FUNCTION m111
    1634 
    1635       REAL FUNCTION m121(T)
    1636 
    1637        REAL, INTENT(IN) :: T
    1638        m121=1114.58541077D0-1.1833078936D0*T
    1639      &    -0.00209946114412D0*T**2-246749.842271D0/T
    1640      &    +34.1234558134D0*dlog(T)
    1641        END FUNCTION m121
    1642 
    1643        FUNCTION m221(T)
    1644 
    1645        REAL, INTENT(IN) :: T
    1646        m221=-80.1488100747D0-0.0116246143257D0*T
    1647      &    +0.606767928954D-5*T**2+3092.72150882D0/T
    1648      &    +12.7601667471D0*dlog(T)
    1649        END FUNCTION m221
    1650 
    1651       REAL FUNCTION m122(T)
    1652 
    1653        REAL, INTENT(IN) :: T
    1654        m122=888.711613784D0-2.50531359687D0*T
    1655      &    +0.000605638824061D0*T**2-196985.296431D0/T
    1656      &    +74.550064338D0*dlog(T)
    1657        END FUNCTION m122
    1658 
    1659       REAL FUNCTION e111(T)
    1660 
    1661        REAL, INTENT(IN) :: T
    1662        e111=2887.31663295D0-3.32602457749D0*T
    1663      &    -0.2820472833D-2*T**2-528216.112353D0/T
    1664      &    +0.68699743564D0*dlog(T)
    1665        END FUNCTION e111
    1666 
    1667       REAL FUNCTION e121(T)
    1668 
    1669        REAL, INTENT(IN) :: T
    1670        e121=-370.944593249D0-0.690310834523D0*T
    1671      &    +0.56345508422D-3*T**2-3822.52997064D0/T
    1672      &    +94.2682037574D0*dlog(T)
    1673        END FUNCTION e121
    1674 
    1675       REAL FUNCTION e211(T)
    1676 
    1677        REAL, INTENT(IN) :: T
    1678        e211=38.3025318809D0-0.0295997878789D0*T
    1679      &    +0.120999746782D-4*T**2-3246.97498999D0/T
    1680      &    -3.83566039532D0*dlog(T)
    1681        END FUNCTION e211
    1682 
    1683       REAL FUNCTION e221(T)
    1684 
    1685        REAL, INTENT(IN) :: T
    1686        e221=2324.76399402D0-0.141626921317D0*T
    1687      &    -0.00626760562881D0*T**2-450590.687961D0/T
    1688      &    -61.2339472744D0*dlog(T)
    1689        END FUNCTION e221
    1690 
    1691       REAL FUNCTION e122(T)
    1692 
    1693        REAL, INTENT(IN) :: T
    1694        e122=-1633.85547832D0-3.35344369968D0*T
    1695      &    +0.00710978119903D0*T**2+198200.003569D0/T
    1696      &    +246.693619189D0*dlog(T)
    1697        END FUNCTION e122
    1698 
    1699       REAL FUNCTION e212(T)
    1700 
    1701        REAL, INTENT(IN) :: T
    1702        e212=1273.75159848D0+1.03333898148D0*T
    1703      &    +0.00341400487633D0*T**2+195290.667051D0/T
    1704      &    -431.737442782D0*dlog(T)
    1705        END FUNCTION e212
    1706 
    1707       REAL FUNCTION lnAa(x1,T)
    1708 
    1709        REAL, INTENT(IN) :: T,x1
    1710        REAL ::
    1711      &          m111,m121,m221,m122
    1712      &            ,e111,e121,e211,e122,e212,e221
    1713        lnAa=-(
    1714      &    (2*m111(T)+e111(T)*(2*dlog(x1)+1))*x1
    1715      &    +(2*m121(T)+e211(T)*dlog(1-x1)+e121(T)*(dlog(x1)+1))*(1-x1)
    1716      &    -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1
    1717      &    -(2*m121(T)+e121(T)*(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1)
    1718      &    -(2*m122(T)+e122(T)*dlog(x1)
    1719      &           +e212(T)*dlog(1-x1))*(1-x1))*x1*(1-x1)
    1720      &    -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2
    1721      &    -x1*(1-x1)*(
    1722      &                  (6*m122(T)+e122(T)*(3*dlog(x1)+1)
    1723      &                          +e212(T)*(3*dlog(1-x1)+1)
    1724      &                   )*x1*(1-x1)
    1725      &                -(2*m122(T)+e122(T)*(dlog(x1)+1)
    1726      &                                   +e212(T)*dlog(1-x1)
    1727      &                    )*(1-x1))
    1728      &     )
    1729        END FUNCTION lnAa
    1730 
    1731       REAL FUNCTION lnAw(x1,T)
    1732 
    1733        REAL, INTENT(IN) :: T,x1
    1734        REAL ::
    1735      &          m111,m121,m221,m122
    1736      &            ,e111,e121,e211,e122,e212,e221
    1737        lnAw=-(
    1738      &  (2*m121(T)+e121(T)*dlog(x1)+e211(T)*(dlog(1-x1)+1))*x1
    1739      &  +(2*m221(T)+e221(T)*(2*dlog(1-x1)+1))*(1-x1)
    1740      &  -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1
    1741      & -(2*m121(T)+e121(T)*(dlog(x1)+1)
    1742      &            +e211(T)*(dlog(1-x1)+1))*x1*(1-x1)
    1743      &        -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2
    1744      &   +x1*(2*m122(T)+e122(T)*dlog(x1)+e212(T)*dlog(1-x1))*x1*(1-x1)
    1745      &  +x1*(1-x1)*((2*m122(T)+e122(T)*dlog(x1)
    1746      &                        +e212(T)*(dlog(1-x1)+1))*x1
    1747      &               -(6*m122(T)+e122(T)*(3*dlog(x1)+1)
    1748      &                        +e212(T)*(3*dlog(1-x1)+1))*(1-x1)*x1)
    1749      &     )
    1750        END FUNCTION lnAw
    1751 
    1752       REAL FUNCTION activitya(xal,T)
    1753 
    1754        REAL, INTENT(IN) :: T,xal
    1755        REAL :: lnAa
    1756 !       &          ,m111,m121,m221,m122 &
    1757 !       &            ,e111,e121,e211,e122,e212,e221
    1758 
    1759 !       write(*,*) 'in activitya ', xal, T
    1760        activitya=DEXP(lnAa(xal,T)-lnAa(1.D0-1.D-12,T))
    1761        END FUNCTION activitya
    1762 
    1763        FUNCTION activityw(xal,T)
    1764 
    1765        REAL, INTENT(IN) :: T,xal
    1766        REAL :: lnAw
    1767 
    1768        activityw=DEXP(lnAw(xal,T)-lnAw(1.D-12,T))
    1769        END FUNCTION activityw
    1770 
    1771 ! end of functions related to zeleznik activities
    1772 
    1773 
    1774 
    1775 
    1776       FUNCTION SIGMADROPLET(xmass,t)
    1777 ! calculates the surface tension of the liquid in J/m^2
    1778 ! xmass=mass fraction of h2so4, t in kelvins
    1779 ! about 230-323 K , x=0,...,1
    1780 !(valid down to the solid phase limit temp, which depends on molefraction)
    1781       IMPLICIT NONE
    1782       REAL :: SIGMADROPLET
    1783       REAL, INTENT(IN):: xmass, t
    1784       REAL :: a, b, t1, tc, xmole
    1785       REAL, PARAMETER :: Msa=98.078
    1786       REAL, PARAMETER :: Mwv=18.0153
    1787 
    1788        IF (t .LT. 305.15) THEN
    1789 !low temperature surface tension
    1790 ! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari
    1791 ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
    1792 ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
    1793 !and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631
    1794       a=0.11864+xmass*(-0.11651+xmass*(0.76852+xmass*
    1795      & (-2.40909+xmass*(2.95434-xmass*1.25852))))
    1796       b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass*
    1797      & (0.007611235+xmass*(-0.00937386+xmass*0.00389722))))
    1798       SIGMADROPLET=a+t*b
    1799       ELSE
    1800 
    1801       xmole = (xmass/Msa)*(1./((xmass/Msa)+(1.-xmass)/Mwv))
    1802 ! high temperature surface tension
    1803 !H. Vehkam‰ki and M. Kulmala and K.E. J. lehtinen, 2003,
    1804 !Modelling binary homogeneous nucleation of water-sulfuric acid vapours:
    1805 ! parameterisation for high temperature emissions, () Environ. Sci. Technol., 37, 3392-3398
    1806 
    1807       tc= 647.15*(1.0-xmole)*(1.0-xmole)+900.0*xmole*xmole+
    1808      & 3156.186*xmole*(1-xmole) !critical temperature
    1809       t1=1.0-t/tc
    1810       a= 0.2358+xmole*(-0.529+xmole*(4.073+xmole*(-12.6707+xmole*
    1811      & (15.3552+xmole*(-6.3138)))))
    1812       b=-0.14738+xmole*(0.6253+xmole*(-5.4808+xmole*(17.2366+xmole*
    1813      & (-21.0487+xmole*(8.719)))))
    1814       SIGMADROPLET=(a+b*t1)*t1**(1.256)
    1815       END IF
    1816 
    1817       RETURN
    1818       END FUNCTION SIGMADROPLET
    1819 
    1820       FUNCTION RHODROPLET(xmass,t)
    1821 !
    1822 ! calculates the density of the liquid in kg/m^3
    1823 ! xmass=mass fraction of h2so4, t in kelvins
    1824 ! Hanna Vehkam‰ki and Markku Kulmala and Ismo Napari
    1825 ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
    1826 ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
    1827 !and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631
    1828 
    1829 ! about 220-373 K , x=0,...,1
    1830 !(valid down to the solid phase limit temp, which depends on molefraction)
    1831 
    1832       IMPLICIT NONE
    1833       REAL :: RHODROPLET
    1834       REAL, INTENT(IN) :: xmass, t
    1835       REAL ::  a,b,c
    1836 
    1837 
    1838       a=0.7681724+xmass*(2.1847140+xmass*(7.1630022+xmass*
    1839      & (-44.31447+xmass*
    1840      & (88.75606+xmass*(-75.73729+xmass*23.43228)))))
    1841       b=1.808225e-3+xmass*(-9.294656e-3+xmass*(-0.03742148+
    1842      &  xmass*(0.2565321+xmass*(-0.5362872+xmass*
    1843      &  (0.4857736-xmass*0.1629592)))))
    1844       c=-3.478524e-6+xmass*(1.335867e-5+xmass*
    1845      & (5.195706e-5+xmass*(-3.717636e-4+xmass*
    1846      & (7.990811e-4+xmass*(-7.458060e-4+xmass*2.58139e-4)))))
    1847       RHODROPLET=a+t*(b+c*t) ! g/cm^3
    1848       RHODROPLET= RHODROPLET*1.0e3 !kg/m^3
    1849       RETURN
    1850       END FUNCTION RHODROPLET

  • trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F

    r1642 r1661  
    341341c Variables tendance sedimentation
    342342
    343       REAL :: d_tr_sed(klon,klev,2)
     343      REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
     344      REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
     345      REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
     346      REAL :: aer_flux(klev)
     347      REAL :: d_tr_sed(klon,klev,nqmax)
    344348      REAL :: d_tr_ssed(klon)
    345349c
    … …  
    640644c===============================================
    641645     
    642       if ((nlon .EQ. 1) .AND. ok_cloud) then
     646c MICROPHY SANS CHIMIE: seulement si full microphy (cl_scheme=2)
     647
     648      if (ok_chem.and..not.ok_cloud) then
     649        print*,"LA CHIMIE A BESOIN DE LA MICROPHYSIQUE"
     650        print*,"ok_cloud doit etre = a ok_chem"
     651      stop
     652      endif
     653      if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.1)) then
     654        print*,"cl_scheme=1 doesnot work without chemistry"
     655      stop
     656      endif
     657       if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
     658        print*,"Full microphysics without chemistry"
     659c indexation of microphys tracers
     660        CALL chemparam_ini()
     661      endif
     662   
     663c number of microphysical tracers
     664      nmicro=0
     665      if (ok_cloud .AND. (cl_scheme.eq.1)) nmicro=2
     666      if (ok_cloud .AND. (cl_scheme.eq.2)) nmicro=8
     667 
     668c CAS 1D POUR MICROPHYS Aurelien
     669      if ((nlon .EQ. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
    643670        PRINT*,'Open profile_cloud_parameters.csv'
    644671        OPEN(66,file='profile_cloud_parameters.csv',
    … …  
    646673      endif
    647674
    648       if ((nlon .EQ. 1) .AND. ok_sedim) then
     675      if ((nlon .EQ. 1) .AND. ok_sedim .AND. (cl_scheme.eq.1)) then
    649676        PRINT*,'Open profile_cloud_sedim.csv'
    650677        OPEN(77,file='profile_cloud_sedim.csv',
    … …  
    652679      endif
    653680           
    654       if ((nlon .GT. 1) .AND. (ok_chem.OR.ok_cloud)) then
     681c INIT PHOTOCHEMISTRY ! includes the indexation of microphys tracers
     682      if ((nlon .GT. 1) .AND. ok_chem) then
    655683c !!! DONC 3D !!!
    656684        CALL chemparam_ini()
    657685      endif
    658686         
    659       if ((nlon .GT. 1) .AND. ok_cloud) then
     687c INIT MICROPHYS SCHEME 1 (AURELIEN) 
     688      if ((nlon .GT. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
    660689c !!! DONC 3D !!!
    661690        CALL cloud_ini(nlon,nlev)
    … …  
    889918c        CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
    890919
     920C === SEDIMENTATION ===
     921
    891922         if (ok_sedim) then
     923
     924c !! Depend on cl_scheme !!
     925
     926          if (cl_scheme.eq.1) then
     927c         ================
    892928
    893929           CALL new_cloud_sedim(
    … …  
    899935     I                     t_seri,
    900936     I                 tr_seri,
    901      O                     d_tr_sed,
     937     O                     d_tr_sed(:,:,1:2),
    902938     O                     d_tr_ssed,
    903939     I                     nqmax,
    … …  
    933969     
    934970        Fsedim(:,klev+1) = 0.
     971       
     972          elseif (cl_scheme.eq.2) then
     973c         ====================
     974
     975           d_tr_sed(:,:,:) = 0.D0
     976
     977           DO i=1, klon
     978
     979c Mode 1
     980         m0_mode1(:,1)=tr_seri(i,:,i_m0_mode1drop)
     981         m0_mode1(:,2)=tr_seri(i,:,i_m0_mode1ccn)
     982         m3_mode1(:,1)=tr_seri(i,:,i_m3_mode1sa)
     983         m3_mode1(:,2)=tr_seri(i,:,i_m3_mode1w)
     984         m3_mode1(:,3)=tr_seri(i,:,i_m3_mode1ccn)
     985         
     986         call drop_sedimentation(dtime,klev,m0_mode1,m3_mode1,
     987     &        paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
     988     &        1,
     989     &        d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
     990
     991        d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
     992     &                              + d_drop_sed
     993        d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
     994     &                              + d_ccn_sed(:,1)
     995        d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
     996     &                              + d_ccn_sed(:,2)
     997        d_tr_sed(i,:,i_m3_mode1sa)  = d_tr_sed(i,:,i_m3_mode1sa)
     998     &                              + d_liq_sed(:,1)
     999        d_tr_sed(i,:,i_m3_mode1w)   = d_tr_sed(i,:,i_m3_mode1w)
     1000     &                              + d_liq_sed(:,2)
     1001
     1002c Mode 2
     1003         m0_mode2(:,1)=tr_seri(i,:,i_m0_mode2drop)
     1004         m0_mode2(:,2)=tr_seri(i,:,i_m0_mode2ccn)
     1005         m3_mode2(:,1)=tr_seri(i,:,i_m3_mode2sa)
     1006         m3_mode2(:,2)=tr_seri(i,:,i_m3_mode2w)
     1007         m3_mode2(:,3)=tr_seri(i,:,i_m3_mode2ccn)
     1008         
     1009         call drop_sedimentation(dtime,klev,m0_mode2,m3_mode2,
     1010     &        paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
     1011     &        2,
     1012     &        d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
     1013
     1014        d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
     1015     &                              + d_drop_sed
     1016        d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
     1017     &                              + d_ccn_sed(:,1)
     1018        d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
     1019     &                              + d_ccn_sed(:,2)
     1020        d_tr_sed(i,:,i_m3_mode2sa)  = d_tr_sed(i,:,i_m3_mode2sa)
     1021     &                              + d_liq_sed(:,1)
     1022        d_tr_sed(i,:,i_m3_mode2w)   = d_tr_sed(i,:,i_m3_mode2w)
     1023     &                              + d_liq_sed(:,2)
     1024
     1025c Aer
     1026        call aer_sedimentation(dtime,klev,tr_seri(i,:,i_m0_aer),
     1027     &        tr_seri(i,:,i_m3_aer),paprs(i,:),
     1028     &        zzlev(i,:),zzlay(i,:),t_seri(i,:),
     1029     &        d_tr_sed(i,:,i_m0_aer),d_tr_sed(i,:,i_m3_aer),aer_flux)
     1030
     1031           END DO
     1032         
     1033           DO iq=nqmax-nmicro+1,nqmax
     1034            tr_seri(:,:,iq)  = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)
     1035            d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq) / dtime
     1036           END DO
     1037
     1038        endif
     1039c         ====================
     1040
     1041C === FIN SEDIMENTATION ===
    9351042
    9361043         endif ! ok_sedim

  • trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F

    r1530 r1661  
    1414     I                    pdtphys,
    1515     I                    temp,
    16      I                    pplev,
     16     I                    pplay,
    1717     O                    trac)
    1818c======================================================================
    … …  
    2929     
    3030#include "clesphys.h"
    31 c#include "temps.h"
    32 c#include "paramet.h"
    33 c#include "comcstfi.h" !me permet de recuperer mugaz et d'autres constantes comme rad,pi etc
    3431#include "YOMCST.h"
    3532c======================================================================
    … …  
    6158      real  trac_sav(n_lon,n_lev,nqmax)
    6259      real  trac_sum(n_lon,n_lev)
    63       real  pplev(n_lon,n_lev)  ! pression pour le mileu de chaque couche (en Pa)
     60      real  pplay(n_lon,n_lev)  ! pression pour le mileu de chaque couche (en Pa)
    6461      real  lon_sun
    6562
    … …  
    9390         PRINT*,'n_lon 1D: ',n_lon
    9491         end if
    95                  
    96 c         if ((n_lon .GT. 1) .AND. ok_chem) then
    97 c !!! DONC 3D !!!       
    98 c           CALL chemparam_ini()
    99 c         endif
    100          
    101 c         if ((n_lon .GT. 1) .AND. ok_cloud) then
    102 c !!! DONC 3D !!!
    103 c         CALL cloud_ini(n_lon,n_lev)
    104 c         endif
    105            
     92
     93c=============================================================
     94c == CASE INIT PHOTOCHEM ==
     95c=============================================================
     96                             
    10697       IF (reinit_trac) THEN
    10798       PRINT*,'REINIT MIXING RATIO TRACEURS'
    10899
    109 c       =============================================================
    110 c                                       Passage de Rm à Rv
     100c                                       Passage de Rm a Rv
    111101c       =============================================================
    112102c       Necessaire si on reprend les start.nc qui sont en MMR
    113 
    114          DO iq=1,nqmax
     103c SEULEMENT LES GAZ
     104         DO iq=1,nqmax-nmicro
    115105          trac(:,:,iq)=trac(:,:,iq)*mmean(:,:)/M_tr(iq)
    116106         END DO
    117107c       =============================================================
    118108         
    119 c=============================================================
    120109c               Initialisation des profils traceurs en Rv
    121110c=============================================================
    … …  
    146135c     On a donc le CO2 qui est le restant d atmosphere Venus 
    147136         trac_sum(:,:)=0.0
    148         DO iq=2,nqmax
     137c SEULEMENT LES GAZ
     138        DO iq=2,nqmax-nmicro
    149139         trac_sum(:,:)= trac_sum(:,:) + trac(:,:,iq)
    150140        END DO
    … …  
    154144c=============================================================
    155145     
    156 c       =============================================================
    157 c                                       Passage de Rv à Rm
    158 c       =============================================================
    159          DO iq=1,nqmax
     146c                                       Passage de Rv a Rm
     147c       =============================================================
     148c SEULEMENT LES GAZ
     149         DO iq=1,nqmax-nmicro
    160150          trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
    161151         END DO
    … …  
    164154       ENDIF  !FIN REINIT TRAC
    165155
     156c=============================================================
     157c == CASE INIT GAZ when muphy without chemistry ==
     158c=============================================================
     159                             
     160       if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
     161
     162c                                       Passage de Rm a Rv
     163c       =============================================================
     164c       Necessaire si on reprend les start.nc qui sont en MMR
     165c SEULEMENT LES GAZ
     166         DO iq=1,nqmax-nmicro
     167          trac(:,:,iq)=trac(:,:,iq)*mmean(:,:)/M_tr(iq)
     168         END DO
     169c       =============================================================
     170         
     171        call vapors4muphy_ini(n_lon,n_lev,trac)
     172
     173c                                       Passage de Rv a Rm
     174c       =============================================================
     175c SEULEMENT LES GAZ
     176         DO iq=1,nqmax-nmicro
     177          trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
     178         END DO
     179c       =============================================================
     180   
     181      endif
    166182       
    167183c-------------
    … …  
    172188
    173189c       =============================================================
    174 c                                       Passage de Rm à Rv
    175 c       =============================================================
    176        DO iq=1,nqmax
    177          trac(:,:,iq)=MAX(trac(:,:,iq)*mmean(:,:)/M_tr(iq),1.E-30)
    178        END DO
    179 c       =============================================================
    180 
    181 
    182 c       =============================================================
    183 c                        Appel Microphysique (sans nucleation)
    184 c                  Volume Mixing Ratio
    185 c       =============================================================
    186 
    187          IF (ok_cloud) THEN
    188                
     190c                                       Passage de Rm a Rv
     191c                                        DES GAZ
     192c       =============================================================
     193      DO iq=1,nqmax-nmicro
     194        trac(:,:,iq)=MAX(trac(:,:,iq)*mmean(:,:)/M_tr(iq),1.E-30)
     195      END DO
     196c       =============================================================
     197
     198
     199c       =============================================================
     200c           Appel Microphysique (schema simple, these Aurelien Stolzenbach)
     201c       =============================================================
     202
     203      IF (ok_cloud .AND. cl_scheme.eq.1) THEN
     204
     205c         Passage de Rm a Rv pour les liq
     206         trac(:,:,i_h2so4liq)=MAX(trac(:,:,i_h2so4liq)
     207     &                          *mmean(:,:)/M_tr(i_h2so4liq),1.E-30)
     208         trac(:,:,i_h2oliq)=MAX(trac(:,:,i_h2oliq)
     209     &                          *mmean(:,:)/M_tr(i_h2oliq),1.E-30)
     210             
    189211c      PRINT*,'DEBUT CLOUD'     
    190212c     On remet tout le RM liq dans la partie gaz
    191 c     !!! On reforme un nuage à chaque fois !!!
     213c     !!! On reforme un nuage a chaque fois !!!
    192214
    193215      DO ilev=1, n_lev
    … …  
    203225                   
    204226      CALL new_cloud_venus(n_lev, n_lon,
    205      e temp,pplev,
     227     e temp,pplay,
    206228     e mrtwv,mrtsa,
    207229     e mrwv,mrsa)
    … …  
    210232c       Actualisation des mixing ratio liq et gaz
    211233c       =========================================
    212 c       Si la routine new_cloud_venus n'a pas actualisé mrwv et mrsa
     234c       Si la routine new_cloud_venus n'a pas actualise mrwv et mrsa
    213235c       on a alors bien mr=mrt pour sa et wv, donc les parties liq sont=0 hors du nuage
    214236c       ou si on ne condense pas
    … …  
    229251
    230252c       =============================================================
    231 c      PRINT*,'FIN CLOUD'
    232       ENDIF
     253c      PRINT*,'FIN CLOUD, scheme 1'
     254      ENDIF 
     255
     256c       =============================================================
     257c           Appel Microphysique (schema complet, these Sabrina Guilbon)
     258c       =============================================================
     259
     260      IF (ok_cloud .AND. cl_scheme.eq.2) THEN
     261
     262c       Boucle sur grille (n_lon) et niveaux (n_lev)
     263        DO ilon=1, n_lon       
     264         DO ilev=1, n_lev
     265
     266           CALL MAD_MUPHY(pdtphys,                               ! Timestep
     267     &  temp(ilon,ilev),pplay(ilon,ilev),                        ! Temperature and pressure
     268     &  trac(ilon,ilev,i_h2o),trac(ilon,ilev,i_h2so4),           ! Mixing ratio of SA and W
     269     &  trac(ilon,ilev,i_m0_aer),trac(ilon,ilev,i_m3_aer),       ! Moments of aerosols
     270     &  trac(ilon,ilev,i_m0_mode1drop),trac(ilon,ilev,i_m0_mode1ccn),  ! Moments of mode 1
     271     &  trac(ilon,ilev,i_m3_mode1sa),trac(ilon,ilev,i_m3_mode1w),      ! Moments of mode 1
     272     &  trac(ilon,ilev,i_m3_mode1ccn),                                 ! Moments of mode 1
     273     &  trac(ilon,ilev,i_m0_mode2drop),trac(ilon,ilev,i_m0_mode2ccn),  ! Moments of mode 2
     274     &  trac(ilon,ilev,i_m3_mode2sa),trac(ilon,ilev,i_m3_mode2w),      ! Moments of mode 2
     275     &  trac(ilon,ilev,i_m3_mode2ccn))                                 ! Moments of mode 2
     276
     277         ENDDO
     278        ENDDO
     279
     280
     281c       =============================================================
     282c      PRINT*,'FIN CLOUD, scheme 2'
     283      ENDIF 
     284
    233285           
    234286c=============================================================
    … …  
    239291c     Ici, la longitude au midi local se deplace vers l'Ouest
    240292c     c'est le sens terrestre
    241 c     pour Vénus on prend juste l'opposé de la longitude et on a la rotation
    242 c     de Vénus et donc le midi local qui se déplace vers l'Est
     293c     pour Venus on prend juste l'oppose de la longitude et on a la rotation
     294c     de Venus et donc le midi local qui se deplace vers l'Est
    243295     
    244296      lon_sun = (0.5 - gmtime) * 2.0 * RPI
    … …  
    246298      lat_local = lat * RPI/180.0E+0
    247299       
     300      IF (ok_chem) THEN
     301
     302c      PRINT*,'DEBUT CHEMISTRY'
    248303      DO ilon=1, n_lon
    249304
    … …  
    255310c      PRINT*,'sza_local :', sza_local     
    256311   
    257       IF (ok_chem) THEN
    258 c      PRINT*,'DEBUT CHEMISTRY'
    259312c       =============================================================
    260313c                                       Appel Photochimie
    261314c       =============================================================
    262 c     Pression en hPa => pplev/100.
     315c     Pression en hPa => pplay/100.
    263316       
    264317      CALL new_photochemistry_venus(n_lev, n_lon, pdtphys,
    265      e                         pplev(ilon,:)/100.,
     318     e                         pplay(ilon,:)/100.,
    266319     e                         temp(ilon,:),
    267320     e                         trac(ilon,:,:),
    … …  
    269322     e                         sza_local, nqmax)
    270323c       =============================================================
     324
     325      END DO
    271326c      PRINT*,'FIN CHEMISTRY'
    272327   
    273         END IF
     328      END IF
     329c=============================================================
     330
     331c       =============================================================
     332c                                       Passage de Rv a Rm
     333c       =============================================================
     334c                                        DES GAZ
     335      DO iq=1,nqmax-nmicro
     336                trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
    274337
    275338      END DO
    276 c       =============================================================
    277 c                                       Passage de Rv à Rm
    278 c       =============================================================
    279         DO iq=1,nqmax
    280 c               trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/RMD
    281                 trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
    282 
    283         END DO
     339
     340      IF (ok_cloud .AND. cl_scheme.eq.1) THEN
     341c         Passage de Rm a Rv pour les liq
     342         trac(:,:,i_h2so4liq)=trac(:,:,i_h2so4liq)*M_tr(i_h2so4liq)
     343     &                          /mmean(:,:)
     344         trac(:,:,i_h2oliq)=trac(:,:,i_h2oliq)*M_tr(i_h2oliq)
     345     &                          /mmean(:,:)
     346      ENDIF
     347     
    284348c       =============================================================   
    285349C      PRINT*,'FIN PHYTRAC'

  • trunk/LMDZ.VENUS/libf/phyvenus/write_histins.h

    r1572 r1661  
    6868       call histwrite_phy(nid_ins,.false.,"tops",itau_w,topsw)
    6969     
    70        IF (ok_cloud) THEN
     70       IF (ok_cloud.and.(cl_scheme.eq.1)) THEN
    7171       
    7272       IF (nb_mode.GE.1) THEN
    … …  
    109109     & rho_droplet)
    110110                ENDIF
    111        IF (ok_sedim) THEN
     111
     112       IF (ok_sedim.and.(cl_scheme.eq.1)) THEN
    112113     
    113114      call histwrite_phy(nid_ins,.false.,"d_tr_sed_H2SO4",itau_w,
    … …  
    118119      call histwrite_phy(nid_ins,.false.,"F_sedim",itau_w,Fsedim)
    119120                ENDIF             
     121
    120122      ENDIF !lev_histins.GE.2
    121123

  • trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h

    r1572 r1661  
    115115      call histwrite_phy(nid_mth,.false.,"dtcond",itau_w,d_t_conduc)
    116116      call histwrite_phy(nid_mth,.false.,"dumolvis",itau_w,d_u_molvis)
    117       call histwrite_phy(nid_mth,.false.,"dvmolvis",itau_w,d_v_molvis)
     117      call histwrite_phy(nid_mth,.false.,"dvmolvis",itau_w,-1.*d_v_molvis)
    118118
    119119c     call histwrite_phy(nid_mth,.false.,"dtec",itau_w,d_t_ec)
Note: See TracChangeset for help on using the changeset viewer.