Changeset 1639


Ignore:
Timestamp:
Dec 6, 2016, 11:50:13 AM (8 years ago)
Author:
slebonnois
Message:

SL, VENUS: optimisation of computer time for photochemistry and radiative transfer

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r1621 r1639  
    262262c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
    263263             IF (zqi_sa(ig,l).LT.0.0D0) THEN
    264                PRINT*,'STOP sedimentation on epuise tout le VMR present'
    265                PRINT*,'couche',ig,'level',l
     264c              PRINT*,'STOP sedimentation on epuise tout le VMR present'
     265c              PRINT*,'couche',ig,'level',l
    266266c               STOP
    267267c              Ce n est pas juste mais il faudrait alors adapter les pas
    … …  
    286286c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
    287287             IF (zqi_wv(ig,l).LT.0.0D0) THEN
    288                PRINT*,'STOP sedimentation on epuise tout le VMR present'
    289                PRINT*,'couche',ig,'level',l
     288c              PRINT*,'STOP sedimentation on epuise tout le VMR present'
     289c              PRINT*,'couche',ig,'level',l
    290290c               STOP
    291291c              Ce n est pas juste mais il faudrait alors adapter les pas

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

    r1442 r1639  
    1919      USE chemparam_mod
    2020      IMPLICIT NONE
    21      
    22 #include "YOMCST.h" 
    23      
     21
     22#include "YOMCST.h"
     23
    2424      INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
    2525      INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales
    26      
     26
    2727!----------------------------------------------------------------------------
    2828!     Ambient air state variables:
    … …  
    3434!----------------------------------------------------------------------------
    3535!     Thermodynamic functions:
    36       REAL :: RHODROPLET 
     36      REAL :: RHODROPLET
    3737!----------------------------------------------------------------------------
    3838!     Auxilary variables:
    … …  
    4747!     Ridder's Method variables:
    4848      REAL :: WVMIN, WVMAX, WVACC
    49      
     49
    5050      INTEGER :: NBROOT
    51      
     51
    5252      INTEGER :: MAXITE
    5353      PARAMETER(MAXITE=20)
    5454
    5555      INTEGER :: NBRAC
    56       PARAMETER(NBRAC=20)
    57      
     56      PARAMETER(NBRAC=5)
     57
    5858      INTEGER :: FLAG
    5959!----------------------------------------------------------------------------
    60 
     60c     Ratio radius shell model du mode 3
     61c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     62c       Ce ratio correspond aux mesures effectuées par J. Cimino (1982), Icarus
     63c     Fixer ce parametre a 0 revient a une gouttelette pure en liquide acide sulfurique
     64      REAL, PARAMETER :: qrad = 0.97
     65      REAL :: qmass
     66c       masse volumique du coeur (kg.m-3)
     67      REAL, PARAMETER :: rho_core = 2500.0
    6168!----------------------------------------------------------------------------
    6269!     External functions needed:
    6370      REAL :: IRFRMWV
    64 !---------------------------------------------------------------------------- 
    65 
    66            
     71!----------------------------------------------------------------------------
     72
     73
    6774! >>> Program starts here:
    6875
    … …  
    7077! These aerosols will then be given an equilibrium composition for the given size distribution
    7178
    72   ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari 
     79  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
    7380  ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
    74   ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric 
     81  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
    7582  !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
    7683
    … …  
    8289      WH2SO4(:,:)=0.0E+0
    8390      rho_droplet(:,:)=0.0E+0
    84                  
     91
    8592      DO ilev=cloudmin, cloudmax
    86       DO ilon=1, nblon 
    87          
     93      DO ilon=1, nblon
     94
    8895!       Boucle sur les modes
    8996        RMODE=0.0E+0
    9097        K_SAV = 0.0
    91        
     98
    9299        DO imode=1, nbr_mode
    93100          IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
    … …  
    106113!       Accuracy de WVeq
    107114        WVACC=WVMAX*1.0E-3
    108                      
     115
    109116!       BRACWV borne la fonction f(WV) - WV = 0
    110117!       de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
    111 !       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0 
     118!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
    112119!       Elle fait appel à la fct/ssrtine ITERWV()
    113      
     120
    114121        CALL BRACWV(WVMIN,WVMAX,NBRAC,RMODE,
    115122     &  mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),TT(ilon,ilev),
    116123     &  PP(ilon,ilev),FLAG,WSAFLAG,NBROOT)
    117        
     124
    118125        SELECT CASE(FLAG)
    119              
    120         CASE(1) 
    121 !         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)         
    122 !       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0 
     126
     127        CASE(1)
     128!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
     129!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
    123130!       Elle fait appel ˆ la fct/ssrtine ITERWV()
    124            
     131
    125132          WH2SO4(ilon,ilev)=IRFRMWV(WVMIN,WVMAX,WVACC,MAXITE,RMODE,
    126133     &    TT(ilon,ilev),PP(ilon,ilev),
    127      &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)             
     134     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
    128135
    129136          rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
    130      &    TT(ilon,ilev)) 
     137     &    TT(ilon,ilev))
    131138
    132139!          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
    … …  
    139146!            STOP
    140147!          ENDIF
    141                
     148
    142149          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
    143150
    144             NH2SO4=mrt_sa(ilon,ilev)*CONCM
    145             NH2O=mrt_wv(ilon,ilev)*CONCM
     151              NH2SO4=mrt_sa(ilon,ilev)*CONCM
     152              NH2O=mrt_wv(ilon,ilev)*CONCM
    146153
    147154          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
    … …  
    150157
    151158!       Boucle sur les modes
    152           DO imode=1, nbr_mode
    153             IF (K_MASS(ilon,ilev,imode).GT.0.) THEN       
     159! *********************************************
     160!       AVEC MODE 3 type J. Cimino 1982, Icarus
     161! *********************************************
     162!       Mode 1 et 2 avec les Kmass modifies
     163          DO imode=1, nbr_mode-1
     164!           calcul qmass
     165              qmass = (rho_core*qrad**3)/
     166     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
     167
     168            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
    154169              NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
    155      &        K_MASS(ilon,ilev,imode)*MCONDTOT*
    156      &        EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
     170     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
     171     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
    157172     &        (R_MEDIAN(ilon,ilev,imode)**3.)
    158173            ELSE
    159174              NBRTOT(ilon,ilev,imode)=0.0E+0
    160             ENDIF     
    161           ENDDO   
     175            ENDIF
     176          ENDDO
     177!       Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
     178          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
     179            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
     180     &      K_MASS(ilon,ilev,3)*MCONDTOT*
     181     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
     182     &      (R_MEDIAN(ilon,ilev,3)**3.)
     183          ELSE
     184              NBRTOT(ilon,ilev,3)=0.0E+0
     185          ENDIF
    162186
    163187!       Passage de #/m3 en VMR
    … …  
    167191          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
    168192          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
    169        
     193
    170194!          Problemes quand on a condense tout, on peut obtenir des -1e-24
    171195!               aprs la soustraction et conversion de ND ˆ VMR
    172           IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30         
     196          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
    173197          IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
    174198
    … …  
    176200
    177201        CASE(2)
    178 !       Cas NROOT=0 mais proche de 0 
     202!       Cas NROOT=0 mais proche de 0
    179203
    180204          WH2SO4(ilon,ilev)=WSAFLAG
    181          
     205
    182206          rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
    183      &    TT(ilon,ilev))         
     207     &    TT(ilon,ilev))
    184208
    185209!     ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
    186210!     ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
    187 !     incohŽrentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
     211!     incoherentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
    188212!     Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
    189213!     donne par RHODROPLET (cf test externe en Python), sinon, la valeur est trop
    190 !     basse (de l'ordre de 1000 kg/m3) et correspond parfois ˆ la valeur avec
    191 !     WSA=0.1 (pas totalement sur) 
     214!     basse (de l'ordre de 1000 kg/m3) et correspond parfois a la valeur avec
     215!     WSA=0.1 (pas totalement sur)
    192216!     En tous cas, incoherent avec ce qui est attendue pour le WSA et T donneeŽ
    193 !     La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a 
     217!     La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
    194218!     la bonne valeur (tests externes Python confirment)
    195      
     219
    196220          IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
    197221     &      (WH2SO4(ilon,ilev).GT.0.1))THEN
    … …  
    203227            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
    204228            STOP
    205           ENDIF 
    206  
    207          
     229          ENDIF
     230
     231
    208232          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
    209233
    … …  
    216240
    217241!       Boucle sur les modes
    218           DO imode=1, nbr_mode
    219             IF (K_MASS(ilon,ilev,imode).GT.0.) THEN       
     242! *********************************************
     243!       AVEC MODE 3 type J. Cimino 1982, Icarus
     244! *********************************************
     245!       Mode 1 et 2 avec alcul coeff*Kmass
     246          DO imode=1, nbr_mode-1
     247!           calcul qmass
     248              qmass = (rho_core*qrad**3)/
     249     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
     250
     251            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
    220252              NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
    221      &        K_MASS(ilon,ilev,imode)*MCONDTOT*
    222      &        EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
     253     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
     254     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
    223255     &        (R_MEDIAN(ilon,ilev,imode)**3.)
    224256            ELSE
    225257              NBRTOT(ilon,ilev,imode)=0.0E+0
    226             ENDIF     
    227           ENDDO   
     258            ENDIF
     259          ENDDO
     260!       Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
     261          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
     262            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
     263     &      K_MASS(ilon,ilev,3)*MCONDTOT*
     264     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
     265     &      (R_MEDIAN(ilon,ilev,3)**3.)
     266          ELSE
     267              NBRTOT(ilon,ilev,3)=0.0E+0
     268          ENDIF
     269
    228270
    229271!       Passage de #/m3 en VMR
    … …  
    233275          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
    234276          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
    235        
     277
    236278!          Problmes quand on a condense tout, on peut obtenir des -1e-24
    237279!               aprs la soustraction et conversion de ND ˆ VMR
    238           IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30         
     280          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
    239281          IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
    240          
     282
    241283        CASE(3)
    242 !         Cas 0 NROOT 
     284!         Cas 0 NROOT
    243285            mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
    244286            mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
    … …  
    246288            WH2SO4(ilon,ilev)=0.0E+0
    247289            DO imode=1, nbr_mode
    248               NBRTOT(ilon,ilev,imode)=0.0E+0   
    249             ENDDO   
    250 
    251         END SELECT 
    252       ENDDO   !FIN boucle ilon 
    253       ENDDO   !FIN boucle ilev 
    254              
     290              NBRTOT(ilon,ilev,imode)=0.0E+0
     291            ENDDO
     292
     293        END SELECT
     294      ENDDO   !FIN boucle ilon
     295      ENDDO   !FIN boucle ilev
     296
    255297      END SUBROUTINE new_cloud_venus
    256      
    257      
     298
     299
    258300!*****************************************************************************
    259 !*    SUBROUTINE ITERWV()                                         
     301!*    SUBROUTINE ITERWV()
    260302      SUBROUTINE ITERWV(WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,
    261303     + TAIR,PAIR,RADIUS)
    262304!*****************************************************************************
    263 !* Cette routine est la solution par itŽration afin de trouver WSA pour un WV,
    264 !* et donc LPPWV, donnŽ. Ce qui nous donne egalement le WV correspondant au
    265 !* WSA solution   
    266 !* For VenusGCM by A. Stolzenbach 07/2014 
     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
    267309!* OUTPUT: WVEQ et WSAOUT
    268310
    269311      IMPLICIT NONE
    270       REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS 
    271      
     312      REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS
     313
    272314      REAl, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ
    273      
     315
    274316      REAL :: WSAMIN, WSAMAX, WSAACC
    275       PARAMETER(WSAACC=0.001)
    276      
     317      PARAMETER(WSAACC=0.01)
     318
    277319      REAL :: LPPWV
    278      
     320
    279321      INTEGER :: MAXITSA, NBRACSA, NBROOT
    280322      PARAMETER(MAXITSA=20)
    281       PARAMETER(NBRACSA=20)
    282      
     323      PARAMETER(NBRACSA=5)
     324
    283325      LOGICAl :: FLAG1,FLAG2
    284326
    285 !     External Function     
     327!     External Function
    286328      REAl :: IRFRMSA, WVCOND
    287      
     329
    288330      IF (RADIUS.LT.1E-30) THEN
    289331        PRINT*,'RMODE == 0 FLAG 3'
    290332        STOP
    291333      ENDIF
    292 !     Initialisation WSA=[0.1,1.0]     
     334!     Initialisation WSA=[0.1,1.0]
    293335      WSAMIN=0.1
    294336      WSAMAX=1.0
    295            
     337
    296338      LPPWV=DLOG(PAIR*WV)
    297            
    298 !     Appel Bracket de KEEQ         
     339
     340!     Appel Bracket de KEEQ
    299341      CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR,
    300342     &     LPPWV,FLAG1,FLAG2,NBROOT)
    301            
    302       IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN         
     343
     344      IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN
    303345!          Appel Ridder's Method
    304346
    305         WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA, 
     347        WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA,
    306348     &  RADIUS,TAIR,PAIR,LPPWV,NBROOT)
    307349!              IF (WSAOUT.EQ.1.0) WSAOUT=0.999999
    … …  
    313355        IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT=WSAMIN
    314356        IF (FLAG1.AND.FLAG2) THEN
    315           PRINT*,'FLAGs BARCWSA tous TRUE'
     357          PRINT*,'FLAGs BRACWSA tous TRUE'
    316358          STOP
    317359        ENDIF
    318360      ENDIF
    319            
    320          
     361
     362
    321363!     WVEQ output correspondant a WVliq lie a WSA calcule
    322364      WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT)
    323365      WVEQOUT=(WVLIQ+WV)/WVTOT-1.0
    324                            
     366
    325367      END SUBROUTINE ITERWV
    326368
    327369
    328370!*****************************************************************************
    329 !*    SUBROUTINE BRACWV()                                         
     371!*    SUBROUTINE BRACWV()
    330372      SUBROUTINE BRACWV(XA,XB,N,RADIUS,WVTOT,SATOT,TAIR,PAIR,
    331373     +           FLAGWV,WSAFLAG,NROOT)
    332374!*****************************************************************************
    333375!* Bracket de ITERWV
    334 !* From Numerical Recipes     
    335 !* Adapted for VenusGCM A. Stolzenbach 07/2014 
     376!* From Numerical Recipes
     377!* Adapted for VenusGCM A. Stolzenbach 07/2014
    336378!* X est WVinput
    337 !* OUTPUT: XA et XB     
     379!* OUTPUT: XA et XB
    338380
    339381      IMPLICIT NONE
    … …  
    341383      REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR,PAIR
    342384      INTEGER, INTENT(IN) :: N
    343      
     385
    344386      REAL, INTENT(INOUT) :: XA,XB
    345387      REAL, INTENT(OUT) :: WSAFLAG
    346      
     388
    347389      INTEGER :: I,J
    348      
     390
    349391      INTEGER, INTENT(OUT) :: NROOT
    350      
     392
    351393      REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT
    352394      REAL :: XMAX,XMIN,WVEQACC
    353      
     395
    354396      INTEGER, INTENT(OUT) :: FLAGWV
    355397
    … …  
    358400!     s'approche de 0 mais ne l'atteint pas.
    359401      WVEQACC=1.0E-3
    360        
     402
    361403      FLAGWV=1
    362404
    363405      NROOT=0
    364406
    365       X=XA
     407!     25/11/2016
     408!     On change ordre on va du max au min
     409      X=XB
    366410      XMAX=XB
    367411      XMIN=XA
    368412
    369413!     CAS 1 On borne la fonction (WVEQ=0)
    370      
     414
    371415      CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
    372416      FP=WVEQ
    373      
    374       DO I=1,N-1
     417
     418      DO I=N-1,1,-1
    375419         X=(1.-DLOG(REAL(N-I))/DLOG(REAL(N)))*XMAX
    376420         CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
    377421         FC=WVEQ
    378422
    379           IF ((FP*FC).LT.0.0) THEN 
     423          IF ((FP*FC).LT.0.0) THEN
    380424            NROOT=NROOT+1
    381 !           Si NROOT>1 on place la borne sup output ˆ la borne min du calcul en i             
     425!           Si NROOT>1 on place la borne sup output ˆ la borne min du calcul en i
    382426            IF (NROOT.GT.1) THEN
    383                XB=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
     427               XB=(1.-DLOG(REAL(I+1))/DLOG(REAL(N)))*XMAX
    384428            ENDIF
    385429
    … …  
    387431              XA=XMIN
    388432            ELSE
    389               XA=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
     433              XA=X
    390434            ENDIF
    391             XB=X
    392          ENDIF         
     435            RETURN
     436         ENDIF
    393437         FP=FC
    394438      ENDDO
    395439
    396 !     CAS 2 on refait la boucle pour tester si WVEQ est proche de 0 
     440!     CAS 2 on refait la boucle pour tester si WVEQ est proche de 0
    397441!     avec le seuil WVEQACC
    398442      IF (NROOT.EQ.0) THEN
    399          X=XMIN
     443         X=XMAX
    400444         CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
    401445     +   TAIR,PAIR,RADIUS)
    402           DO J=1,N-1
     446          DO J=N-1,1,-1
    403447             X=(1.-DLOG(REAL(N-J))/DLOG(REAL(N)))*XMAX
    404448             CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
    405449     +       TAIR,PAIR,RADIUS)
    406              
     450
    407451             IF (ABS(WVEQ).LE.WVEQACC) THEN
    408452                WSAFLAG=WSAOUT
    … …  
    412456          ENDDO
    413457
    414 !     CAS 3 Pas de borne, WVEQ jamais proche de 0         
     458!     CAS 3 Pas de borne, WVEQ jamais proche de 0
    415459          FLAGWV=3
    416460          RETURN
    … …  
    418462
    419463      END SUBROUTINE BRACWV
    420      
     464
    421465!*****************************************************************************
    422 !*    SUBROUTINE BRACWSA()                                         
     466!*    SUBROUTINE BRACWSA()
    423467      SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL,
    424468     +           NROOT)
    425469!*****************************************************************************
    426470!* Bracket de KEEQ
    427 !* From Numerical Recipes     
    428 !* Adapted for VenusGCM A. Stolzenbach 07/2014 
    429      
     471!* From Numerical Recipes
     472!* Adapted for VenusGCM A. Stolzenbach 07/2014
     473
    430474      IMPLICIT NONE
    431475
    … …  
    433477!     External functions needed:
    434478      REAl KEEQ
    435 !---------------------------------------------------------------------------- 
    436 
    437       REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP 
     479!----------------------------------------------------------------------------
     480
     481      REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP
    438482      INTEGER, INTENT(IN) :: N
    439      
     483
    440484      REAL, INTENT(INOUT) :: XA,XB
    441      
     485
    442486      INTEGER, INTENT(OUT) ::  NROOT
    443487
    444488      INTEGER :: I, J
    445      
     489
    446490      REAL :: DX, FP, FC, X
    447      
     491
    448492      LOGICAL, INTENT(OUT) :: FLAGH,FLAGL
    449    
    450      
     493
     494
    451495      FLAGL=.FALSE.
    452       FLAGH=.FALSE.   
     496      FLAGH=.FALSE.
    453497      NROOT=0
    454498      DX=(XB-XA)/N
    455       X=XA
     499      X=XB
    456500      FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    457      
    458       DO I=1,N
    459          X=X+DX
     501
     502      DO I=N,1,-1
     503         X=X-DX
    460504         FC=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    461          
     505
    462506         IF ((FP*FC).LE.0.) THEN
    463507            NROOT=NROOT+1
    464             XA=X-DX
    465             XB=X
    466 !            RETURN
     508            XA=X
     509            XB=X+DX
     510            RETURN
    467511!            IF (NROOT.GT.1) THEN
    468512!               PRINT*,'On a plus d1 intervalle KEEQ=0'
    … …  
    479523!            ENDIF
    480524         ENDIF
    481          
     525
    482526         FP=FC
    483527      ENDDO
    484      
     528
    485529      IF (NROOT.EQ.0) THEN
    486530!         PRINT*,'On a 0 intervalle KEEQ=0'
    … …  
    492536!         PRINT*,'NBRAC',N
    493537!         STOP
    494          
     538
    495539!         X=XA
    496540!         FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
    … …  
    503547
    504548
    505 !        Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]       
     549!        Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]
    506550         IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
    507551     &    ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE.
    … …  
    511555!        STOP
    512556      ENDIF
    513      
     557
    514558      END SUBROUTINE BRACWSA
    515          
    516            
     559
     560
    517561!*****************************************************************************
    518 !*     REAL FUNCTION WVCOND()                                         
     562!*     REAL FUNCTION WVCOND()
    519563      REAL FUNCTION WVCOND(WSA,T,P,SAt)
    520564!*****************************************************************************
    … …  
    527571!     T: temperature (K)
    528572!     P: pressure (Pa)
    529 !     OUTPUT: 
     573!     OUTPUT:
    530574!       WVCOND : VMR H2O condense
    531575
    532576!      USE chemparam_mod
    533      
     577
    534578      IMPLICIT NONE
    535579
    … …  
    545589      REAL  NH2SO4
    546590      REAL  H2OCOND, H2SO4COND
    547    
     591
    548592
    549593      CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS!
    550594
    551595        NH2SO4=SAt*CONCM
    552      
     596
    553597      pstand=1.01325E+5 !Pa  1 atm pressure
    554598
    … …  
    565609
    566610!acid sat.vap.PP over mixture (flat surface):
    567         satpacid=act(2)*acidps ! Pa 
     611        satpacid=act(2)*acidps ! Pa
    568612
    569613!       Conversion from Pa to N.D #/m3
    570614        DND2=satpacid/(1.3806488E-23*T)
    571                
     615
    572616!       H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
    573617        IF (NH2SO4.GT.DND2) THEN
    … …  
    576620        H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
    577621
    578 !       Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30 
     622!       Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30
    579623        ELSE
    580624        H2OCOND=1.0E-30*CONCM
    … …  
    582626
    583627!*****************************************************
    584 !     ATTENTION: Ici on ne prends pas en compte 
     628!     ATTENTION: Ici on ne prends pas en compte
    585629!                si H2O en defaut!
    586630!                On veut la situation theorique
    587631!                a l'equilibre
    588 !*****************************************************         
     632!*****************************************************
    589633!       Test si H2O en defaut H2Ocond>H2O dispo
    590634!       IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN
    591        
     635
    592636!       On peut alors condenser tout le H2O dispo
    593637!       H2OCOND=NH2O
    594638!       On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
    595639!       H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
    596            
     640
    597641!      END IF
    598642
    599 !     Calcul de H2O condensŽe VMR         
     643!     Calcul de H2O condensŽe VMR
    600644      WVCOND=H2OCOND/CONCM
    601      
     645
    602646      END FUNCTION WVCOND
    603647
    604648!*****************************************************************************
    605 !*     REAL FUNCTION IRFRMWV()                                         
     649!*     REAL FUNCTION IRFRMWV()
    606650      REAL      FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,
    607651     + WVTOT,SATOT,NROOT)
    … …  
    616660
    617661      IMPLICIT NONE
    618      
     662
    619663      REAL, INTENT(IN) ::    X1, X2
    620       REAL, INTENT(IN) ::    XACC 
     664      REAL, INTENT(IN) ::    XACC
    621665      INTEGER, INTENT(IN) :: MAXIT,NROOT
    622      
     666
    623667!     LOCAL VARIABLES
    624668      REAL :: XL, XH, XM, XNEW, X
    … …  
    627671      REAL :: ANS, S, FSIGN
    628672      INTEGER i
    629      
     673
    630674!     External variables needed:
    631675      REAL, INTENT(IN) :: TAIR,PAIR
    … …  
    633677      REAL, INTENT(IN) :: RADIUS
    634678
    635      
     679
    636680!     Initialisation
    637681      X=X1
    638       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS) 
     682      CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
    639683      FL=WVEQ
    640684      X=X2
    641685      CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
    642686      FH=WVEQ
    643      
    644 !     Test Bracketed values 
     687
     688!     Test Bracketed values
    645689      IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.
    646690     &   ((FL.GT.0.).AND.(FH.LT.0.)))
    … …  
    649693         XH=X2
    650694         ANS=-9.99e99
    651          
     695
    652696         DO i=1, MAXIT
    653697            XM=0.5*(XL+XH)
    … …  
    656700            FM=WVEQ
    657701            S=SQRT(FM*FM-FL*FH)
    658            
     702
    659703            IF (S.EQ.0.0) THEN
    660704               IRFRMWV=WSALOC
    661705               RETURN
    662             ENDIF   
    663              
     706            ENDIF
     707
    664708            IF (FL.GT.FH) THEN
    665709               FSIGN=1.0
    … …  
    667711               FSIGN=-1.0
    668712            ENDIF
    669            
    670             XNEW=XM+(XM-XL)*(FSIGN*FM/S) 
    671            
     713
     714            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
     715
    672716            IF (ABS(XNEW-ANS).LE.XACC) THEN
    673717               IRFRMWV=WSALOC
    674718               RETURN
    675             ENDIF   
    676            
     719            ENDIF
     720
    677721            ANS=XNEW
    678722            CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    679723     &       TAIR,PAIR,RADIUS)
    680724            FNEW=WVEQ
    681            
     725
    682726            IF (FNEW.EQ.0.0) THEN
    683727               IRFRMWV=WSALOC
    684728               RETURN
    685729            ENDIF
    686            
     730
    687731            IF (SIGN(FM, FNEW).NE.FM) THEN
    688732               XL=XM
    … …  
    708752         XH=X2
    709753         ANS=-9.99e99
    710          
     754
    711755         DO i=1, MAXIT
    712756            XM=0.5*(XL+XH)
    … …  
    720764               FSIGN=-1.0
    721765            ENDIF
    722            
    723             XNEW=XM+(XM-XL)*(FSIGN*FM/S)     
    724            
     766
     767            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
     768
    725769            ANS=XNEW
    726770            CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    … …  
    753797         PRINT*,'NROOT de BRACWV', NROOT
    754798         IF (ABS(FL).LT.XACC) THEN
    755          PRINT*,'IRFRMWV FL == 0',FL 
     799         PRINT*,'IRFRMWV FL == 0',FL
    756800         PRINT*,'X1',X1,'X2',X2,'FH',FH
    757801           CALL ITERWV(X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
    … …  
    768812         RETURN
    769813         ENDIF
    770          IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN 
     814         IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN
    771815           PRINT*,'STOP dans IRFRMWV avec rien == 0'
    772816           PRINT*,'X1',X1,'X2',X2
    773817           PRINT*,'Fcalc',FL,FH
    774818           PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS
    775            STOP   
     819           STOP
    776820         ENDIF
    777          IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN 
     821         IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN
    778822           PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC'
    779823           PRINT*,FL,FH
    780            STOP   
     824           STOP
    781825         ENDIF
    782          
    783          
     826
     827
    784828      ENDIF
    785829!  FIN Test Bracketed values
    786        
     830
    787831      END FUNCTION IRFRMWV
    788                  
     832
    789833!*****************************************************************************
    790 !*     REAL FUNCTION IRFRMSA()                                         
     834!*     REAL FUNCTION IRFRMSA()
    791835      REAL      FUNCTION IRFRMSA(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,LPPWV,
    792836     +               NB)
    … …  
    801845
    802846      IMPLICIT NONE
    803      
     847
    804848      REAL, INTENT(IN) ::    X1, X2
    805       REAL, INTENT(IN) ::    XACC 
     849      REAL, INTENT(IN) ::    XACC
    806850      INTEGER, INTENT(IN) :: MAXIT, NB
    807      
     851
    808852!     LOCAL VARIABLES
    809853      REAL XL, XH, XM, XNEW
    … …  
    811855      REAL ANS, S, FSIGN
    812856      INTEGER i
    813      
     857
    814858!     External variables needed:
    815859      REAL, INTENT(IN) :: TAIR,PAIR
    816860      REAL, INTENT(IN) :: LPPWV
    817861      REAL, INTENT(IN) :: RADIUS
    818      
     862
    819863!     External functions needed:
    820864      REAL KEEQ
    821865
    822866
    823      
     867
    824868!     Initialisation
    825869      FL=KEEQ(RADIUS,TAIR,X1,LPPWV)
    826870      FH=KEEQ(RADIUS,TAIR,X2,LPPWV)
    827      
    828 !     Test Bracketed values 
     871
     872!     Test Bracketed values
    829873      IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.((FL.GT.0.).AND.(FH.LT.0.)))
    830874     &  THEN
    … …  
    832876         XH=X2
    833877         ANS=-9.99e99
    834          
     878
    835879         DO i=1, MAXIT
    836880            XM=0.5*(XL+XH)
    837881            FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
    838882            S=SQRT(FM*FM-FL*FH)
    839            
     883
    840884            IF (S.EQ.0.0) THEN
    841885               IRFRMSA=ANS
    842886               RETURN
    843             ENDIF   
    844              
     887            ENDIF
     888
    845889            IF (FL.GT.FH) THEN
    846890               FSIGN=1.0
    … …  
    848892               FSIGN=-1.0
    849893            ENDIF
    850            
    851             XNEW=XM+(XM-XL)*(FSIGN*FM/S) 
    852            
     894
     895            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
     896
    853897            IF (ABS(XNEW-ANS).LE.XACC) THEN
    854898               IRFRMSA=ANS
    855899               RETURN
    856             ENDIF   
    857            
     900            ENDIF
     901
    858902            ANS=XNEW
    859903            FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
    860            
     904
    861905            IF (FNEW.EQ.0.0) THEN
    862906               IRFRMSA=ANS
    863907               RETURN
    864908            ENDIF
    865            
     909
    866910            IF (SIGN(FM, FNEW).NE.FM) THEN
    867911               XL=XM
    … …  
    888932         PRINT*,'Borne XL',XL,'XH',XH
    889933         ANS=-9.99e99
    890          
     934
    891935         DO i=1, MAXIT
    892936            XM=0.5*(XL+XH)
    893937            FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
    894938            S=SQRT(FM*FM-FL*FH)
    895  
     939
    896940            IF (FL.GT.FH) THEN
    897941               FSIGN=1.0
    … …  
    899943               FSIGN=-1.0
    900944            ENDIF
    901            
    902             XNEW=XM+(XM-XL)*(FSIGN*FM/S) 
    903            
     945
     946            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
     947
    904948            ANS=XNEW
    905949            FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
    … …  
    937981         IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN
    938982           PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH
    939            PRINT*,'X1',X1,'X2',X2 
     983           PRINT*,'X1',X1,'X2',X2
    940984           PRINT*,'Kind F', KIND(FL), KIND(FH)
    941985           PRINT*,'Kind X', KIND(X1), KIND(X2)
    … …  
    944988           PRINT*,'nb root BRACWSA',NB
    945989           STOP
    946          ENDIF   
    947      
     990         ENDIF
     991
    948992      ENDIF
    949993!  FIN Test Bracketed values
    950        
     994
    951995      END function IRFRMSA
    952      
     996
    953997!*****************************************************************************
    954 !*     REAL FUNCTION KEEQ()                                         
     998!*     REAL FUNCTION KEEQ()
    955999      REAL      FUNCTION KEEQ(RADIUS,TAIR,WSA,LPPWV)
    9561000!*****************************************************************************
    … …  
    9851029     +        C1=2.0d0*MH2O/RGAS)
    9861030
    987      
     1031
    9881032      KEEQ=LPPWV-C1*SIGMADROPLET(WSA,TAIR)/
    9891033     &     (TAIR*RADIUS*RHODROPLET(WSA,TAIR))-
    9901034     &     PWVSAS_GV(TAIR,WSA)
    991      
     1035
    9921036      END FUNCTION KEEQ
    993      
     1037
    9941038*****************************************************************************
    9951039*     REAL FUNCTION PWVSAS_GV(TAIR,WSA)                                       
    … …  
    10481092     +     K2=K1*K1/2.0)
    10491093!
    1050 !
    1051       CP=CPH2O(WSA)
    1052       F=-FFH2O(WSA)
    1053       L=-LH2O(WSA)
    1054       ALFA=ALH2O(WSA)
     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)
    10551129!
    10561130      A=ADOT+(CP-K1*ALFA)/RGAS
    … …  
    10701144*******************************************************************************
    10711145*     REAL FUNCTION CPH2O(W)
    1072       REAL FUNCTION CPH2O(W)
     1146      REAL FUNCTION CPH2O(W,khi_in,klo_in)
    10731147*******************************************************************************
    10741148*
    … …  
    10861160     +              Y2(NPOINT),YWORK(NPOINT)
    10871161      REAL, INTENT(IN):: W
     1162      INTEGER, INTENT(IN):: khi_in,klo_in
     1163      INTEGER :: khi,klo
    10881164      REAL :: CPH
    10891165      LOGICAL :: FIRST
    … …  
    11251201          CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2)
    11261202      ENDIF
    1127       CALL SPLINT(WTAB,CPHTAB,Y2,NPOINT,W,CPH)
     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)
    11281214      CPH2O=CPH
    11291215     
    … …  
    11311217!
    11321218*******************************************************************************
    1133       REAL FUNCTION FFH2O(W)
     1219      REAL FUNCTION FFH2O(W,khi,klo)
    11341220*     REAL FUNCTION FFH2O(W)
    11351221*******************************************************************************
    … …  
    11471233      REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
    11481234      REAL, INTENT(IN) :: W
     1235      INTEGER, INTENT(IN):: khi,klo
    11491236      REAL :: FF
    11501237      LOGICAL :: FIRST
    … …  
    11861273          CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2)
    11871274      ENDIF
    1188       CALL SPLINT(WTAB,FFTAB,Y2,NPOINT,W,FF)
     1275
     1276      CALL SPLINT(WTAB(khi),WTAB(klo),FFTAB(khi),FFTAB(klo),
     1277     .              Y2(khi),Y2(klo),W,FF)
    11891278      FFH2O=FF
    11901279     
    … …  
    11921281!
    11931282*******************************************************************************
    1194       REAL FUNCTION LH2O(W)
     1283      REAL FUNCTION LH2O(W,khi,klo)
    11951284*     REAL FUNCTION LH2O(W)
    11961285*******************************************************************************
    … …  
    12081297      REAL, DIMENSION(NPOINT) ::  WTAB,LTAB,Y2,YWORK
    12091298      REAL, INTENT(IN) :: W
     1299      INTEGER, INTENT(IN):: khi,klo
    12101300      REAL :: L
    12111301      LOGICAL :: FIRST
    … …  
    12471337          CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2)
    12481338      ENDIF
    1249       CALL SPLINT(WTAB,LTAB,Y2,NPOINT,W,L)
     1339
     1340      CALL SPLINT(WTAB(khi),WTAB(klo),LTAB(khi),LTAB(klo),
     1341     .              Y2(khi),Y2(klo),W,L)
    12501342      LH2O=L
    12511343     
    12521344      END FUNCTION LH2O
    12531345*******************************************************************************
    1254       REAL FUNCTION ALH2O(W)
     1346      REAL FUNCTION ALH2O(W,khi_in,klo_in)
    12551347*     REAL FUNCTION ALH2O(W)
    12561348*******************************************************************************
    … …  
    12681360      REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
    12691361      REAL, INTENT(IN) :: W
     1362      INTEGER, INTENT(IN):: khi_in,klo_in
     1363      INTEGER :: khi,klo
    12701364      REAL :: A
    12711365      LOGICAL :: FIRST
    … …  
    13041398          CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2)
    13051399      ENDIF
    1306       CALL SPLINT(WTAB,ATAB,Y2,NPOINT,MAX(WTAB(1),W),A)
     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)
    13071411      ALH2O=A
    13081412     
    … …  
    13581462
    13591463!******************************************************************************
    1360       SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
     1464      SUBROUTINE SPLINT(XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo,X,Y)
    13611465!******************************************************************************
    13621466!     Cubic spline calculation
    … …  
    13641468      IMPLICIT NONE
    13651469
    1366       INTEGER, INTENT(IN) :: N
    1367       INTEGER :: KLO,KHI,K
    1368       REAL, INTENT(IN), DIMENSION(N) :: XA,YA,Y2A
     1470      REAL, INTENT(IN) :: XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo
    13691471      REAL, INTENT(IN) :: X
    13701472      REAL, INTENT(OUT) :: Y
    13711473      REAL :: H,A,B
    13721474!
    1373       KLO=1
    1374       KHI=N
    1375  1    IF(KHI-KLO.GT.1) THEN
    1376           K=(KHI+KLO)/2
    1377           IF(XA(K).GT.X) THEN
    1378               KHI=K
    1379           ELSE
    1380               KLO=K
    1381           ENDIF
    1382           GOTO 1
    1383       ENDIF
    1384       H=XA(KHI)-XA(KLO)
    1385       A=(XA(KHI)-X)/H
    1386       B=(X-XA(KLO))/H
    1387       Y=A*YA(KLO)+B*YA(KHI)+
    1388      +        ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0d0
     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
    13891480!
    13901481
    … …  
    13941485     + T,H2SO4COND,H2OCOND,RMTOT)
    13951486
    1396 !     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION) 
    1397 !     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL 
     1487!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
     1488!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
    13981489!                                       SIZE DISTRIBTUION
    13991490!     ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
    … …  
    14091500!     T: temperature (K)
    14101501!
    1411 !     OUTPUT: 
     1502!     OUTPUT:
    14121503!     RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
    1413 !            mais rho_droplet et M_tot_distrib doivent tre de meme dimension
     1504!            mais rho_droplet et M_tot_distrib doivent etre de meme dimension
    14141505!       H2OCOND
    14151506!       H2SO4COND
    14161507
    14171508
    1418      
     1509
    14191510      IMPLICIT NONE
    14201511
    … …  
    14301521!     masse of an H2SO4 molecule (kg)
    14311522      RMH2S4=98.078/(6.02214129E+26)
    1432      
     1523
    14331524      pstand=1.01325E+5 !Pa  1 atm pressure
    14341525
    … …  
    14451536
    14461537!acid sat.vap.PP over mixture (flat surface):
    1447         satpacid=act(2)*acidps ! Pa 
     1538        satpacid=act(2)*acidps ! Pa
    14481539
    14491540!       Conversion from Pa to N.D #/m3
    … …  
    14511542!       Conversion from N.D #/m3 TO #/cm3
    14521543!        DND2=DND2*1.d-6
    1453                
     1544
    14541545!       H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
    14551546        IF (H2SO4.GE.DND2) THEN
    … …  
    14581549        H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
    14591550
    1460 !     RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3
     1551!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
     1552!     RMTOT = Mtot/rho_droplet
    14611553!       RMTOT=M_distrib/rho_droplet
    1462        
     1554
    14631555        RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
    14641556
    … …  
    14691561        RMTOT=0.0E+0
    14701562        END IF
    1471                
     1563
    14721564!       Test si H2O en defaut H2Ocond>H2O dispo
    14731565        IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
    … …  
    14841576!      PRINT*,'WSA',WSA,'RHO',DENSO4
    14851577!      STOP
    1486      
    1487        
     1578
     1579
    14881580!       On peut alors condenser tout le H2O dispo
    14891581        H2OCOND=H2O
    14901582!       On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
    14911583        H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
    1492        
    1493 !     RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3
    1494 !       RMTOT=Volume of aerosol cm3 /cm3 of air
    1495 !       Volume of aerosol/cm3 air
    1496        
     1584
     1585!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
     1586!       RMTOT=Volume of aerosol m3 /m3 of air
     1587!       Volume of aerosol/m3 air
     1588
    14971589        RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
    1498      
     1590
    14991591      END IF
    1500                
     1592
    15011593      END SUBROUTINE CALCM_SAT
    15021594
    … …  
    15101602  !     mole fraction 0,...,1
    15111603  !     J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991
    1512   !+++++++++++++++++++++++++++++++++++++++++++++++++++ 
     1604  !+++++++++++++++++++++++++++++++++++++++++++++++++++
    15131605
    15141606         IMPLICIT NONE
    … …  
    15231615!         write(*,*) 'x, T ', x, T
    15241616
    1525          act(2)=activitya(x,T) 
     1617         act(2)=activitya(x,T)
    15261618         act(1)=activityw(x,T)
    15271619
    … …  
    15351627
    15361628       REAL, INTENT(IN) :: T
    1537        m111=-23.524503387D0 
    1538      &    +0.0406889449841D0*T 
    1539      &    -0.151369362907D-4*T**2+2961.44445015D0/T 
     1629       m111=-23.524503387D0
     1630     &    +0.0406889449841D0*T
     1631     &    -0.151369362907D-4*T**2+2961.44445015D0/T
    15401632     &    +0.492476973663D0*dlog(T)
    15411633       END FUNCTION m111
    … …  
    15441636
    15451637       REAL, INTENT(IN) :: T
    1546        m121=1114.58541077D0-1.1833078936D0*T 
    1547      &    -0.00209946114412D0*T**2-246749.842271D0/T 
     1638       m121=1114.58541077D0-1.1833078936D0*T
     1639     &    -0.00209946114412D0*T**2-246749.842271D0/T
    15481640     &    +34.1234558134D0*dlog(T)
    15491641       END FUNCTION m121
    … …  
    15521644
    15531645       REAL, INTENT(IN) :: T
    1554        m221=-80.1488100747D0-0.0116246143257D0*T 
    1555      &    +0.606767928954D-5*T**2+3092.72150882D0/T 
     1646       m221=-80.1488100747D0-0.0116246143257D0*T
     1647     &    +0.606767928954D-5*T**2+3092.72150882D0/T
    15561648     &    +12.7601667471D0*dlog(T)
    15571649       END FUNCTION m221
    … …  
    15601652
    15611653       REAL, INTENT(IN) :: T
    1562        m122=888.711613784D0-2.50531359687D0*T 
    1563      &    +0.000605638824061D0*T**2-196985.296431D0/T 
     1654       m122=888.711613784D0-2.50531359687D0*T
     1655     &    +0.000605638824061D0*T**2-196985.296431D0/T
    15641656     &    +74.550064338D0*dlog(T)
    15651657       END FUNCTION m122
    … …  
    15681660
    15691661       REAL, INTENT(IN) :: T
    1570        e111=2887.31663295D0-3.32602457749D0*T 
    1571      &    -0.2820472833D-2*T**2-528216.112353D0/T 
     1662       e111=2887.31663295D0-3.32602457749D0*T
     1663     &    -0.2820472833D-2*T**2-528216.112353D0/T
    15721664     &    +0.68699743564D0*dlog(T)
    15731665       END FUNCTION e111
    … …  
    15761668
    15771669       REAL, INTENT(IN) :: T
    1578        e121=-370.944593249D0-0.690310834523D0*T 
    1579      &    +0.56345508422D-3*T**2-3822.52997064D0/T 
     1670       e121=-370.944593249D0-0.690310834523D0*T
     1671     &    +0.56345508422D-3*T**2-3822.52997064D0/T
    15801672     &    +94.2682037574D0*dlog(T)
    15811673       END FUNCTION e121
    … …  
    15841676
    15851677       REAL, INTENT(IN) :: T
    1586        e211=38.3025318809D0-0.0295997878789D0*T 
    1587      &    +0.120999746782D-4*T**2-3246.97498999D0/T 
     1678       e211=38.3025318809D0-0.0295997878789D0*T
     1679     &    +0.120999746782D-4*T**2-3246.97498999D0/T
    15881680     &    -3.83566039532D0*dlog(T)
    15891681       END FUNCTION e211
    … …  
    15921684
    15931685       REAL, INTENT(IN) :: T
    1594        e221=2324.76399402D0-0.141626921317D0*T 
    1595      &    -0.00626760562881D0*T**2-450590.687961D0/T 
     1686       e221=2324.76399402D0-0.141626921317D0*T
     1687     &    -0.00626760562881D0*T**2-450590.687961D0/T
    15961688     &    -61.2339472744D0*dlog(T)
    15971689       END FUNCTION e221
    … …  
    16001692
    16011693       REAL, INTENT(IN) :: T
    1602        e122=-1633.85547832D0-3.35344369968D0*T 
    1603      &    +0.00710978119903D0*T**2+198200.003569D0/T 
     1694       e122=-1633.85547832D0-3.35344369968D0*T
     1695     &    +0.00710978119903D0*T**2+198200.003569D0/T
    16041696     &    +246.693619189D0*dlog(T)
    16051697       END FUNCTION e122
    … …  
    16081700
    16091701       REAL, INTENT(IN) :: T
    1610        e212=1273.75159848D0+1.03333898148D0*T 
    1611      &    +0.00341400487633D0*T**2+195290.667051D0/T 
     1702       e212=1273.75159848D0+1.03333898148D0*T
     1703     &    +0.00341400487633D0*T**2+195290.667051D0/T
    16121704     &    -431.737442782D0*dlog(T)
    16131705       END FUNCTION e212
    … …  
    16161708
    16171709       REAL, INTENT(IN) :: T,x1
    1618        REAL :: 
    1619      &          m111,m121,m221,m122 
     1710       REAL ::
     1711     &          m111,m121,m221,m122
    16201712     &            ,e111,e121,e211,e122,e212,e221
    1621        lnAa=-( 
    1622      &    (2*m111(T)+e111(T)*(2*dlog(x1)+1))*x1 
    1623      &    +(2*m121(T)+e211(T)*dlog(1-x1)+e121(T)*(dlog(x1)+1))*(1-x1) 
    1624      &    -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 
    1625      &    -(2*m121(T)+e121(T)*(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1) 
    1626      &    -(2*m122(T)+e122(T)*dlog(x1) 
    1627      &           +e212(T)*dlog(1-x1))*(1-x1))*x1*(1-x1) 
    1628      &    -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 
    1629      &    -x1*(1-x1)*( 
    1630      &                  (6*m122(T)+e122(T)*(3*dlog(x1)+1) 
    1631      &                          +e212(T)*(3*dlog(1-x1)+1) 
    1632      &                   )*x1*(1-x1) 
    1633      &                -(2*m122(T)+e122(T)*(dlog(x1)+1) 
    1634      &                                   +e212(T)*dlog(1-x1) 
    1635      &                    )*(1-x1)) 
     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))
    16361728     &     )
    16371729       END FUNCTION lnAa
    … …  
    16401732
    16411733       REAL, INTENT(IN) :: T,x1
    1642        REAL :: 
    1643      &          m111,m121,m221,m122 
     1734       REAL ::
     1735     &          m111,m121,m221,m122
    16441736     &            ,e111,e121,e211,e122,e212,e221
    1645        lnAw=-( 
    1646      &  (2*m121(T)+e121(T)*dlog(x1)+e211(T)*(dlog(1-x1)+1))*x1 
    1647      &  +(2*m221(T)+e221(T)*(2*dlog(1-x1)+1))*(1-x1) 
    1648      &  -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 
    1649      & -(2*m121(T)+e121(T)*(dlog(x1)+1) 
    1650      &            +e211(T)*(dlog(1-x1)+1))*x1*(1-x1) 
    1651      &        -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 
    1652      &   +x1*(2*m122(T)+e122(T)*dlog(x1)+e212(T)*dlog(1-x1))*x1*(1-x1) 
    1653      &  +x1*(1-x1)*((2*m122(T)+e122(T)*dlog(x1) 
    1654      &                        +e212(T)*(dlog(1-x1)+1))*x1 
    1655      &               -(6*m122(T)+e122(T)*(3*dlog(x1)+1) 
    1656      &                        +e212(T)*(3*dlog(1-x1)+1))*(1-x1)*x1) 
     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)
    16571749     &     )
    16581750       END FUNCTION lnAw
    … …  
    16711763       FUNCTION activityw(xal,T)
    16721764
    1673        REAL, INTENT(IN) :: T,xal 
     1765       REAL, INTENT(IN) :: T,xal
    16741766       REAL :: lnAw
    16751767
    … …  
    17021794      a=0.11864+xmass*(-0.11651+xmass*(0.76852+xmass*
    17031795     & (-2.40909+xmass*(2.95434-xmass*1.25852))))
    1704       b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass* 
     1796      b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass*
    17051797     & (0.007611235+xmass*(-0.00937386+xmass*0.00389722))))
    17061798      SIGMADROPLET=a+t*b
    … …  
    17451837
    17461838      a=0.7681724+xmass*(2.1847140+xmass*(7.1630022+xmass*
    1747      & (-44.31447+xmass* 
     1839     & (-44.31447+xmass*
    17481840     & (88.75606+xmass*(-75.73729+xmass*23.43228)))))
    17491841      b=1.808225e-3+xmass*(-9.294656e-3+xmass*(-0.03742148+
    … …  
    17571849      RETURN
    17581850      END FUNCTION RHODROPLET
    1759 
    1760 
    1761 
    1762 

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

    r1591 r1639  
    2525c With extension NLTE (G. Gilli, 2014)
    2626
     27c Ksi matrices latitudinaly interpolated  (I. Garate-Lopez, 2016)
     28
    2729c======================================================================
    2830      use dimphy
    … …  
    7476      integer mat0,lat,ips,isza,ips0,isza0
    7577      real    factp,factz,ksi
    76 
     78c ------- for lat-interp ----------------
     79        integer mat0A, mat0B, latA, latB, kasua
     80        integer ipsA, ipsB, iszaA, iszaB, ips0A, ips0B, isza0A, isza0B
     81        real    lat_deg, latA_deg, latB_deg
     82        real    factlat, k1, k2, k3, k4
     83c --------------------------------------
    7784      logical firstcall
    7885      data    firstcall/.true./
    7986      save    firstcall
    8087     
     88cERROR          ! For checking if the file it's being read
    8189c-------------------------------------------
    8290c  Initialisations
    … …  
    99107      ENDDO
    100108
     109
    101110c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
    102111      DO j = 1, klon
    … …  
    110119        zcool(k) = 0.0
    111120       ENDDO
     121c       zheat(1:klev)=0.0       !Explicit loop (no change in performance)
     122c       zcool(1:klev)=0.0
     123
    112124       DO k = 1, klev+1
    113125        ZFLNET(k) = 0.0
    114126        ZFSNET(k) = 0.0
    115127       ENDDO
     128c       ZFLNET(1:klev+1)=0.0
     129c       ZFSNET(1:klev+1)=0.0
     130
    116131       ztopsw = 0.0
    117132       ztoplw = 0.0
    … …  
    120135       zsollwdown = 0.0
    121136     
    122          zfract = fract(j)
    123          zrmu0 = rmu0(j)
    124  
     137       zfract = fract(j)
     138       zrmu0 = rmu0(j)
     139
    125140       DO k = 1, klev+1
    126141         PPB(k) = paprs(j,k)/1.e5
    … …  
    132147       ENDDO
    133148       pt0(klev+1) = 0.
    134  
     149      
    135150       DO k = 0,klev+1
    136151       DO i = 0,klev+1
    … …  
    139154       ENDDO
    140155       ENDDO
    141        
     156
    142157c======================================================================
    143158c Getting psi and deltapsi
    … …  
    160175c finding the right matrixes
    161176c --------------------------
    162        mat0 = 0
    163        if (nlatve.eq.1) then   ! clouds are taken as uniform
    164          lat=1
    165        else
    166          if (abs(latitude_deg(j)).le.50.) then
    167            lat=1
    168          elseif (abs(latitude_deg(j)).le.60.) then
    169            lat=2
    170          elseif (abs(latitude_deg(j)).le.70.) then
    171            lat=3
    172          elseif (abs(latitude_deg(j)).le.80.) then
    173            lat=4
    174          else
    175            lat=5
    176          endif
    177        endif
    178        
     177
     178        mat0  = 0
     179        mat0A = 0
     180        mat0B = 0
     181
     182c    Latitude
     183c    --------
     184        lat  = 0
     185        latA = 0
     186        latB = 0
     187
     188c       write(*,*) 'nlatve:', nlatve
     189
     190        lat_deg = abs(latitude_deg(j))
     191
     192c       if (nlatve.eq.1) then   ! clouds are taken as uniform
     193        if ((nlatve.eq.1).or.(lat_deg.le.25.)) then
     194            lat  = 1
     195        elseif (lat_deg.le.50.) then
     196            lat  = 1
     197            latA = 1
     198            latB = 2
     199            latA_deg = 25.0
     200            latB_deg = 55.0
     201        elseif (lat_deg.le.55.) then
     202            lat  = 2
     203            latA = 1
     204            latB = 2
     205            latA_deg = 25.0
     206            latB_deg = 55.0
     207        elseif (lat_deg.le.60.) then
     208            lat  = 2
     209            latA = 2
     210            latB = 3
     211            latA_deg = 55.0
     212            latB_deg = 65.0
     213        elseif (lat_deg.le.65.) then
     214            lat  = 3
     215            latA = 2
     216            latB = 3
     217            latA_deg = 55.0
     218            latB_deg = 65.0
     219        elseif (lat_deg.le.70.) then
     220            lat  = 3
     221            latA = 3
     222            latB = 4
     223            latA_deg = 65.0
     224            latB_deg = 75.0
     225        elseif (lat_deg.le.75.) then
     226            lat  = 4
     227            latA = 3
     228            latB = 4
     229            latA_deg = 65.0
     230            latB_deg = 75.0
     231        elseif (lat_deg.le.80.) then
     232            lat  = 4
     233            latA = 4
     234            latB = 5
     235            latA_deg = 75.0
     236            latB_deg = 85.0
     237        elseif (lat_deg.le.85.) then
     238            lat = 5
     239            latA = 4
     240            latB = 5
     241            latA_deg = 75.0
     242            latB_deg = 85.0
     243        else
     244            lat = 5
     245        endif
     246
     247c        write(*,*) 'Lat',lat,'LatA',latA,'LatB',latB
     248
     249        factlat = 0
     250        if (latA.gt.0) then
     251          factlat = (lat_deg - latA_deg) / (latB_deg - latA_deg)
     252        endif
     253
     254c       write (*,*) 'Factor de correccion:', factlat
     255
     256
     257c    Pressure at Surface
     258c    -------------------   
     259
    179260       ips0=0
    180        if (nbpsve(lat).gt.1) then
     261       ips0A=0
     262       ips0B=0
     263       if (nbpsve(lat).gt.1) then            ! Interpolation on ps
    181264       do ips=1,nbpsve(lat)-1
    182265         if (  (psurfve(ips,lat).ge.paprs(j,1))
    183      .    .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
     266     .  .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
    184267              ips0  = ips
    185268c             print*,'ig=',j,'  ips0=',ips
    … …  
    192275        ips0=1
    193276       endif
     277
     278       if (latA.eq.lat) then
     279         ips0A=ips0
     280       else
     281         if (latA.gt.0) then
     282           if (nbpsve(latA).gt.1) then
     283             do ipsA=1,nbpsve(latA)-1
     284               if (  (psurfve(ipsA,latA).ge.paprs(j,1))
     285     .        .and.(psurfve(ipsA+1,latA).lt.paprs(j,1)) ) then
     286                ips0A = ipsA
     287                exit
     288               endif
     289             enddo
     290           else            ! Only one ps, no interpolation
     291             ips0A=1
     292           endif     ! nbpsve(latA).gt.1
     293         endif       ! latA.gt.0   (if latA=0 ips0A is not used, so it doesn't matter)
     294       endif         ! latA.eq.lat
     295
     296       if (latB.eq.lat) then
     297         ips0B=ips0
     298       else
     299         if (latB.gt.0) then
     300           if (nbpsve(latB).gt.1) then
     301             do ipsB=1,nbpsve(latB)-1
     302               if (  (psurfve(ipsB,latB).ge.paprs(j,1))
     303     .        .and.(psurfve(ipsB+1,latB).lt.paprs(j,1)) ) then
     304                 ips0B = ipsB
     305                 exit
     306               endif
     307             enddo
     308           else
     309             ips0B=1
     310           endif     ! nbpsve(latB).gt.1
     311         endif       ! latB.gt.0   (if latB=0 ips0B is not used, so it doesn't matter)
     312       endif         ! latB.eq.lat
     313
     314
     315c    Solar Zenith Angle
     316c    ------------------
     317
    194318       isza0=0
     319       isza0A=0
     320       isza0B=0
    195321       if (nbszave(lat).gt.1) then
    196322        do isza=1,nbszave(lat)-1
    197323         if (  (szave(isza,lat).ge.zrmu0)
    198      .    .and.(szave(isza+1,lat).lt.zrmu0) ) then
     324     .  .and.(szave(isza+1,lat).lt.zrmu0) ) then
    199325              isza0  = isza
    200326c             print*,'ig=',j,'  isza0=',isza
    … …  
    207333        isza0=-99
    208334       endif
     335
     336
     337       if (latA.eq.lat) then
     338         isza0A=isza0
     339       else
     340         if (latA.gt.0) then
     341           if (nbszave(latA).gt.1) then
     342             do iszaA=1,nbszave(latA)-1
     343               if (  (szave(iszaA,latA).ge.zrmu0)
     344     .        .and.(szave(iszaA+1,latA).lt.zrmu0) ) then
     345                 isza0A = iszaA
     346                 exit
     347               endif
     348             enddo
     349           else            ! Only one sza, no interpolation
     350             isza0A=-99
     351           endif     ! nbszave(latA).gt.1
     352         endif       ! latA.gt.0   (if latA=0 isza0A is not used, so it doesn't matter)
     353       endif         ! latA.eq.lat
     354
     355       if (latB.eq.lat) then
     356         isza0B=isza0
     357       else
     358         if (latB.gt.0) then
     359           if (nbszave(latB).gt.1) then
     360             do iszaB=1,nbszave(latB)-1
     361               if (  (szave(iszaB,latB).ge.zrmu0)
     362     .        .and.(szave(iszaB+1,latB).lt.zrmu0) ) then
     363                 isza0B = iszaB
     364                 exit
     365               endif
     366             enddo
     367           else            ! Only one sza, no interpolation
     368             isza0B=-99
     369           endif     ! nbszave(latB).gt.1
     370         endif       ! latB.gt.0   (if latB=0 isza0B is not used, so it doesn't matter)
     371       endif         ! latB.eq.lat
     372
     373c        write(*,*) 'nbszave', nbszave(lat),'nbpsve(lat)',nbpsve(lat)
     374
    209375       
    210376c -------- Probleme aux bords
    … …  
    215381     .               /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
    216382       endif
    217        if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1))) 
    218      .   then
     383       if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1)))
     384     . then
    219385              ips0  = nbpsve(lat)-1
    220386              print*,'Extrapolation! ig=',j,'  ips0=',ips0
    … …  
    230396         stop
    231397       endif
    232        
     398
     399
     400
    233401       if (isza0.eq.-99) then
    234          mat0 = indexve(lat)+ips0
     402         mat0  = indexve(lat) +ips0
     403         mat0A = indexve(latA)+ips0A
     404         mat0B = indexve(latB)+ips0B
    235405       else
    236          mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0
     406         mat0  = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
     407         mat0A = indexve(latA)+(isza0A-1)*nbpsve(latA)+ips0A
     408         mat0B = indexve(latB)+(isza0B-1)*nbpsve(latB)+ips0B
    237409       endif
    238        
     410
     411c        write(*,*) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
     412   
    239413c interpolation of ksi and computation of psi,deltapsi
    240414c ----------------------------------------------------
    241        do band=1,nnuve
    242         do k=0,klev+1
    243          do i=0,klev+1
    244           if (isza0.eq.-99) then
    245            ksi = ksive(i,k,band,mat0)*(1-factp)
    246      .          +ksive(i,k,band,mat0+1)*factp
    247           else
    248            ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
    249      .          +ksive(i,k,band,mat0+1)*factp  *(1-factz)
    250      .          +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
    251      .          +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
    252           endif
    253      
    254           psi(i,k) = psi(i,k) +
    255      .               RPI*ksi*(bplck(i,band)-bplck(k,band))
    256           deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
    257          enddo
    258         enddo
    259        enddo
     415
     416       if (isza0.eq.-99) then
     417         if (latA.gt.0) then            ! Not being in the two extremal bins
     418
     419           do band=1,nnuve
     420             do k=0,klev+1
     421               do i=k+1,klev+1
     422                 k1 = ksive(i,k,band,mat0A)*(1-factlat)
     423     .              + ksive(i,k,band,mat0B)*factlat
     424                 k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
     425     .              + ksive(i,k,band,mat0B+1)*factlat
     426                 ksi = k1*(1-factp) + k2*factp
     427                 psi(i,k) = psi(i,k) +
     428     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
     429c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
     430c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
     431c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
     432
     433                 kasua=1
     434               enddo
     435             enddo
     436           enddo
     437             do k=0,klev+1
     438               do i=k+1,klev+1
     439                 psi(k,i) = -psi(i,k)
     440               enddo
     441             enddo
     442
     443         else           ! latA=0 --> extremal bins
     444
     445           do band=1,nnuve
     446             do k=0,klev+1
     447               do i=k+1,klev+1
     448                 ksi = ksive(i,k,band,mat0)*(1-factp)
     449     .               + ksive(i,k,band,mat0+1)*factp
     450                 psi(i,k) = psi(i,k) +
     451     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
     452c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
     453c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
     454c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
     455
     456                 kasua=2
     457               enddo
     458             enddo
     459           enddo
     460             do k=0,klev+1
     461               do i=k+1,klev+1
     462                 psi(k,i) = -psi(i,k)
     463               enddo
     464             enddo
     465
     466         endif          ! latA.gt.0
     467
     468       else             ! isza0=!-99
     469
     470         if (latA.gt.0) then            ! Not being in the two extremal bins
     471
     472           do band=1,nnuve
     473             do k=0,klev+1
     474               do i=k+1,klev+1
     475                 k1 = ksive(i,k,band,mat0A)*(1-factlat)
     476     .              + ksive(i,k,band,mat0B)*factlat
     477                 k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
     478     .              + ksive(i,k,band,mat0B+1)*factlat
     479                 k3 = ksive(i,k,band,mat0A+nbpsve(latA))*(1-factlat)
     480     .              + ksive(i,k,band,mat0B+nbpsve(latB))*factlat
     481                 k4 = ksive(i,k,band,mat0A+nbpsve(latA)+1)*(1-factlat)
     482     .              + ksive(i,k,band,mat0B+nbpsve(latB)+1)*factlat
     483                 ksi = ( k1*(1-factp) + k2*factp )*(1-factz)
     484     .              + ( k3*(1-factp) + k4*factp )*factz
     485                 psi(i,k) = psi(i,k) +
     486     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
     487c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
     488c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
     489c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
     490
     491                kasua=3
     492               enddo
     493             enddo
     494           enddo
     495             do k=0,klev+1
     496               do i=k+1,klev+1
     497                 psi(k,i) = -psi(i,k)
     498               enddo
     499             enddo
     500
     501         else           ! latA=0 --> extremal bins
     502
     503           do band=1,nnuve
     504             do k=0,klev+1
     505               do i=k+1,klev+1
     506                 ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
     507     .               + ksive(i,k,band,mat0+1)*factp  *(1-factz)
     508     .               + ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
     509     .               + ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
     510                 psi(i,k) = psi(i,k) +
     511     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
     512c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
     513c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
     514c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
     515
     516                 kasua=4
     517               enddo
     518             enddo
     519           enddo
     520             do k=0,klev+1
     521               do i=k+1,klev+1
     522                 psi(k,i) = -psi(i,k)
     523               enddo
     524             enddo
     525
     526         endif          ! latA.gt.0
     527       endif            ! isza0.eq.-99
     528
     529c       write(*,*) 'Kasua:', kasua
    260530
    261531c======================================================================
    … …  
    358628      RETURN
    359629      END
    360 
Note: See TracChangeset for help on using the changeset viewer.