Index: /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F	(revision 1638)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F	(revision 1639)
@@ -262,6 +262,6 @@
 c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
              IF (zqi_sa(ig,l).LT.0.0D0) THEN
-               PRINT*,'STOP sedimentation on epuise tout le VMR present'
-               PRINT*,'couche',ig,'level',l
+c              PRINT*,'STOP sedimentation on epuise tout le VMR present'
+c              PRINT*,'couche',ig,'level',l
 c               STOP
 c              Ce n est pas juste mais il faudrait alors adapter les pas
@@ -286,6 +286,6 @@
 c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
              IF (zqi_wv(ig,l).LT.0.0D0) THEN
-               PRINT*,'STOP sedimentation on epuise tout le VMR present'
-               PRINT*,'couche',ig,'level',l
+c              PRINT*,'STOP sedimentation on epuise tout le VMR present'
+c              PRINT*,'couche',ig,'level',l
 c               STOP
 c              Ce n est pas juste mais il faudrait alors adapter les pas
Index: /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F	(revision 1638)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F	(revision 1639)
@@ -19,10 +19,10 @@
       USE chemparam_mod
       IMPLICIT NONE
-      
-#include "YOMCST.h" 
-      
+
+#include "YOMCST.h"
+
       INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
       INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales
-      
+
 !----------------------------------------------------------------------------
 !     Ambient air state variables:
@@ -34,5 +34,5 @@
 !----------------------------------------------------------------------------
 !     Thermodynamic functions:
-      REAL :: RHODROPLET 
+      REAL :: RHODROPLET
 !----------------------------------------------------------------------------
 !     Auxilary variables:
@@ -47,22 +47,29 @@
 !     Ridder's Method variables:
       REAL :: WVMIN, WVMAX, WVACC
-      
+
       INTEGER :: NBROOT
-      
+
       INTEGER :: MAXITE
       PARAMETER(MAXITE=20)
 
       INTEGER :: NBRAC
-      PARAMETER(NBRAC=20)
-      
+      PARAMETER(NBRAC=5)
+
       INTEGER :: FLAG
 !----------------------------------------------------------------------------
-
+c     Ratio radius shell model du mode 3
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c	Ce ratio correspond aux mesures effectuées par J. Cimino (1982), Icarus
+c     Fixer ce parametre a 0 revient a une gouttelette pure en liquide acide sulfurique
+      REAL, PARAMETER :: qrad = 0.97
+      REAL :: qmass
+c	masse volumique du coeur (kg.m-3)
+      REAL, PARAMETER :: rho_core = 2500.0
 !----------------------------------------------------------------------------
 !     External functions needed:
       REAL :: IRFRMWV
-!---------------------------------------------------------------------------- 
-
-           
+!----------------------------------------------------------------------------
+
+
 ! >>> Program starts here:
 
@@ -70,7 +77,7 @@
 ! These aerosols will then be given an equilibrium composition for the given size distribution
 
-  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari 
+  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
   ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
-  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric 
+  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
   !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
 
@@ -82,12 +89,12 @@
       WH2SO4(:,:)=0.0E+0
       rho_droplet(:,:)=0.0E+0
-                  
+
       DO ilev=cloudmin, cloudmax
-      DO ilon=1, nblon 
-          
+      DO ilon=1, nblon
+
 !       Boucle sur les modes
         RMODE=0.0E+0
         K_SAV = 0.0
-        
+
         DO imode=1, nbr_mode
           IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
@@ -106,27 +113,27 @@
 !       Accuracy de WVeq
         WVACC=WVMAX*1.0E-3
-                      
+
 !       BRACWV borne la fonction f(WV) - WV = 0
 !       de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
-!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0 
+!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
 !       Elle fait appel à la fct/ssrtine ITERWV()
-      
+
         CALL BRACWV(WVMIN,WVMAX,NBRAC,RMODE,
      &  mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),TT(ilon,ilev),
      &  PP(ilon,ilev),FLAG,WSAFLAG,NBROOT)
-        
+
         SELECT CASE(FLAG)
-             
-        CASE(1) 
-!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)         
-!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0 
+
+        CASE(1)
+!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
+!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
 !       Elle fait appel  la fct/ssrtine ITERWV()
-           
+
           WH2SO4(ilon,ilev)=IRFRMWV(WVMIN,WVMAX,WVACC,MAXITE,RMODE,
      &    TT(ilon,ilev),PP(ilon,ilev),
-     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)             
+     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
 
           rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
-     &    TT(ilon,ilev)) 
+     &    TT(ilon,ilev))
 
 !          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
@@ -139,9 +146,9 @@
 !            STOP
 !          ENDIF
-                
+
           CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
 
-	    NH2SO4=mrt_sa(ilon,ilev)*CONCM
-	    NH2O=mrt_wv(ilon,ilev)*CONCM
+	      NH2SO4=mrt_sa(ilon,ilev)*CONCM
+	      NH2O=mrt_wv(ilon,ilev)*CONCM
 
           CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
@@ -150,14 +157,31 @@
 
 !       Boucle sur les modes
-          DO imode=1, nbr_mode
-            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN       
+! *********************************************
+!	AVEC MODE 3 type J. Cimino 1982, Icarus
+! *********************************************
+!	Mode 1 et 2 avec les Kmass modifies
+          DO imode=1, nbr_mode-1
+!           calcul qmass
+	      qmass = (rho_core*qrad**3)/
+     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
+
+            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
               NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
-     &        K_MASS(ilon,ilev,imode)*MCONDTOT*
-     &        EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
+     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
+     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
      &        (R_MEDIAN(ilon,ilev,imode)**3.)
             ELSE
               NBRTOT(ilon,ilev,imode)=0.0E+0
-            ENDIF     
-          ENDDO   
+            ENDIF
+          ENDDO
+!	Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
+          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
+            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
+     &      K_MASS(ilon,ilev,3)*MCONDTOT*
+     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
+     &      (R_MEDIAN(ilon,ilev,3)**3.)
+          ELSE
+              NBRTOT(ilon,ilev,3)=0.0E+0
+          ENDIF
 
 !       Passage de #/m3 en VMR
@@ -167,8 +191,8 @@
           mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
           mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
-        
+
 !          Problemes quand on a condense tout, on peut obtenir des -1e-24
 !		aprs la soustraction et conversion de ND  VMR
-          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30          
+          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
           IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
 
@@ -176,22 +200,22 @@
 
         CASE(2)
-!       Cas NROOT=0 mais proche de 0 
+!       Cas NROOT=0 mais proche de 0
 
           WH2SO4(ilon,ilev)=WSAFLAG
-          
+
           rho_droplet(ilon,ilev)=RHODROPLET(WH2SO4(ilon,ilev),
-     &    TT(ilon,ilev))         
+     &    TT(ilon,ilev))
 
 !     ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
 !     ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
-!     incohrentes avec TT et WH2SO4 (a priori lorsque NTOT=0) 
+!     incoherentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
 !     Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
 !     donne par RHODROPLET (cf test externe en Python), sinon, la valeur est trop
-!     basse (de l'ordre de 1000 kg/m3) et correspond parfois  la valeur avec
-!     WSA=0.1 (pas totalement sur) 
+!     basse (de l'ordre de 1000 kg/m3) et correspond parfois a la valeur avec
+!     WSA=0.1 (pas totalement sur)
 !     En tous cas, incoherent avec ce qui est attendue pour le WSA et T donnee
-!     La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a 
+!     La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
 !     la bonne valeur (tests externes Python confirment)
-      
+
           IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
      &      (WH2SO4(ilon,ilev).GT.0.1))THEN
@@ -203,7 +227,7 @@
             PRINT*,'FLAG',FLAG,'NROOT',NBROOT
             STOP
-          ENDIF 
- 
-          
+          ENDIF
+
+
           CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
 
@@ -216,14 +240,32 @@
 
 !       Boucle sur les modes
-          DO imode=1, nbr_mode
-            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN       
+! *********************************************
+!	AVEC MODE 3 type J. Cimino 1982, Icarus
+! *********************************************
+!	Mode 1 et 2 avec alcul coeff*Kmass
+          DO imode=1, nbr_mode-1
+!           calcul qmass
+	      qmass = (rho_core*qrad**3)/
+     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
+
+            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
               NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
-     &        K_MASS(ilon,ilev,imode)*MCONDTOT*
-     &        EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
+     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
+     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
      &        (R_MEDIAN(ilon,ilev,imode)**3.)
             ELSE
               NBRTOT(ilon,ilev,imode)=0.0E+0
-            ENDIF     
-          ENDDO   
+            ENDIF
+          ENDDO
+!	Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
+          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
+            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
+     &      K_MASS(ilon,ilev,3)*MCONDTOT*
+     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
+     &      (R_MEDIAN(ilon,ilev,3)**3.)
+          ELSE
+              NBRTOT(ilon,ilev,3)=0.0E+0
+          ENDIF
+
 
 !       Passage de #/m3 en VMR
@@ -233,12 +275,12 @@
           mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
           mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
-        
+
 !          Problmes quand on a condense tout, on peut obtenir des -1e-24
 !		aprs la soustraction et conversion de ND  VMR
-          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30          
+          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
           IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
-          
+
         CASE(3)
-!         Cas 0 NROOT  
+!         Cas 0 NROOT
             mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
             mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
@@ -246,62 +288,62 @@
             WH2SO4(ilon,ilev)=0.0E+0
             DO imode=1, nbr_mode
-              NBRTOT(ilon,ilev,imode)=0.0E+0    
-            ENDDO   
-
-        END SELECT 
-      ENDDO   !FIN boucle ilon 
-      ENDDO   !FIN boucle ilev 
-             
+              NBRTOT(ilon,ilev,imode)=0.0E+0
+            ENDDO
+
+        END SELECT
+      ENDDO   !FIN boucle ilon
+      ENDDO   !FIN boucle ilev
+
       END SUBROUTINE new_cloud_venus
-      
-      
+
+
 !*****************************************************************************
-!*    SUBROUTINE ITERWV()                                         
+!*    SUBROUTINE ITERWV()
       SUBROUTINE ITERWV(WV,WVLIQ,WVEQOUT,WVTOT,WSAOUT,SATOT,
      + TAIR,PAIR,RADIUS)
 !*****************************************************************************
-!* Cette routine est la solution par itration afin de trouver WSA pour un WV,
-!* et donc LPPWV, donn. Ce qui nous donne egalement le WV correspondant au 
-!* WSA solution   
-!* For VenusGCM by A. Stolzenbach 07/2014 
+!* Cette routine est la solution par iteration afin de trouver WSA pour un WV,
+!* et donc LPPWV, donne. Ce qui nous donne egalement le WV correspondant au
+!* WSA solution
+!* For VenusGCM by A. Stolzenbach 07/2014
 !* OUTPUT: WVEQ et WSAOUT
 
       IMPLICIT NONE
-      REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS 
-      
+      REAL, INTENT(IN) :: WV, WVTOT, SATOT, TAIR, PAIR, RADIUS
+
       REAl, INTENT(OUT) :: WVEQOUT, WSAOUT, WVLIQ
-      
+
       REAL :: WSAMIN, WSAMAX, WSAACC
-      PARAMETER(WSAACC=0.001)
-      
+      PARAMETER(WSAACC=0.01)
+
       REAL :: LPPWV
-      
+
       INTEGER :: MAXITSA, NBRACSA, NBROOT
       PARAMETER(MAXITSA=20)
-      PARAMETER(NBRACSA=20)
-      
+      PARAMETER(NBRACSA=5)
+
       LOGICAl :: FLAG1,FLAG2
 
-!     External Function      
+!     External Function
       REAl :: IRFRMSA, WVCOND
-      
+
       IF (RADIUS.LT.1E-30) THEN
         PRINT*,'RMODE == 0 FLAG 3'
         STOP
       ENDIF
-!     Initialisation WSA=[0.1,1.0]      
+!     Initialisation WSA=[0.1,1.0]
       WSAMIN=0.1
       WSAMAX=1.0
-           
+
       LPPWV=DLOG(PAIR*WV)
-           
-!     Appel Bracket de KEEQ         
+
+!     Appel Bracket de KEEQ
       CALL BRACWSA(WSAMIN,WSAMAX,NBRACSA,RADIUS,TAIR,
      &     LPPWV,FLAG1,FLAG2,NBROOT)
-           
-      IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN          
+
+      IF ((.NOT.FLAG1).AND.(.NOT.FLAG2).AND.(NBROOT.EQ.1)) THEN
 !          Appel Ridder's Method
 
-        WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA, 
+        WSAOUT=IRFRMSA(WSAMIN,WSAMAX,WSAACC,MAXITSA,
      &  RADIUS,TAIR,PAIR,LPPWV,NBROOT)
 !              IF (WSAOUT.EQ.1.0) WSAOUT=0.999999
@@ -313,27 +355,27 @@
         IF (FLAG2.AND.(.NOT.FLAG1)) WSAOUT=WSAMIN
         IF (FLAG1.AND.FLAG2) THEN
-          PRINT*,'FLAGs BARCWSA tous TRUE'
+          PRINT*,'FLAGs BRACWSA tous TRUE'
           STOP
         ENDIF
       ENDIF
-           
-          
+
+
 !     WVEQ output correspondant a WVliq lie a WSA calcule
       WVLIQ=WVCOND(WSAOUT,TAIR,PAIR,SATOT)
       WVEQOUT=(WVLIQ+WV)/WVTOT-1.0
-                           
+
       END SUBROUTINE ITERWV
 
 
 !*****************************************************************************
-!*    SUBROUTINE BRACWV()                                         
+!*    SUBROUTINE BRACWV()
       SUBROUTINE BRACWV(XA,XB,N,RADIUS,WVTOT,SATOT,TAIR,PAIR,
      +           FLAGWV,WSAFLAG,NROOT)
 !*****************************************************************************
 !* Bracket de ITERWV
-!* From Numerical Recipes     
-!* Adapted for VenusGCM A. Stolzenbach 07/2014 
+!* From Numerical Recipes
+!* Adapted for VenusGCM A. Stolzenbach 07/2014
 !* X est WVinput
-!* OUTPUT: XA et XB      
+!* OUTPUT: XA et XB
 
       IMPLICIT NONE
@@ -341,15 +383,15 @@
       REAL, INTENT(IN) :: WVTOT,SATOT,RADIUS,TAIR,PAIR
       INTEGER, INTENT(IN) :: N
-      
+
       REAL, INTENT(INOUT) :: XA,XB
       REAL, INTENT(OUT) :: WSAFLAG
-      
+
       INTEGER :: I,J
-      
+
       INTEGER, INTENT(OUT) :: NROOT
-      
+
       REAL :: FP, FC, X, WVEQ, WVLIQ, WSAOUT
       REAL :: XMAX,XMIN,WVEQACC
-      
+
       INTEGER, INTENT(OUT) :: FLAGWV
 
@@ -358,28 +400,30 @@
 !     s'approche de 0 mais ne l'atteint pas.
       WVEQACC=1.0E-3
-       
+
       FLAGWV=1
 
       NROOT=0
 
-      X=XA
+!     25/11/2016
+!     On change ordre on va du max au min
+      X=XB
       XMAX=XB
       XMIN=XA
 
 !     CAS 1 On borne la fonction (WVEQ=0)
-      
+
       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
       FP=WVEQ
-      
-      DO I=1,N-1
+
+      DO I=N-1,1,-1
          X=(1.-DLOG(REAL(N-I))/DLOG(REAL(N)))*XMAX
          CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,TAIR,PAIR,RADIUS)
          FC=WVEQ
 
-          IF ((FP*FC).LT.0.0) THEN  
+          IF ((FP*FC).LT.0.0) THEN
             NROOT=NROOT+1
-!           Si NROOT>1 on place la borne sup output  la borne min du calcul en i             
+!           Si NROOT>1 on place la borne sup output  la borne min du calcul en i
             IF (NROOT.GT.1) THEN
-               XB=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
+               XB=(1.-DLOG(REAL(I+1))/DLOG(REAL(N)))*XMAX
             ENDIF
 
@@ -387,22 +431,22 @@
               XA=XMIN
             ELSE
-              XA=(1.-DLOG(REAL(N-I+1))/DLOG(REAL(N)))*XMAX
+              XA=X
             ENDIF
-            XB=X
-         ENDIF         
+            RETURN
+         ENDIF
          FP=FC
       ENDDO
 
-!     CAS 2 on refait la boucle pour tester si WVEQ est proche de 0 
+!     CAS 2 on refait la boucle pour tester si WVEQ est proche de 0
 !     avec le seuil WVEQACC
       IF (NROOT.EQ.0) THEN
-         X=XMIN
+         X=XMAX
          CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
      +   TAIR,PAIR,RADIUS)
-          DO J=1,N-1
+          DO J=N-1,1,-1
              X=(1.-DLOG(REAL(N-J))/DLOG(REAL(N)))*XMAX
              CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,
      +       TAIR,PAIR,RADIUS)
-             
+
              IF (ABS(WVEQ).LE.WVEQACC) THEN
                 WSAFLAG=WSAOUT
@@ -412,5 +456,5 @@
           ENDDO
 
-!     CAS 3 Pas de borne, WVEQ jamais proche de 0          
+!     CAS 3 Pas de borne, WVEQ jamais proche de 0
           FLAGWV=3
           RETURN
@@ -418,14 +462,14 @@
 
       END SUBROUTINE BRACWV
-      
+
 !*****************************************************************************
-!*    SUBROUTINE BRACWSA()                                         
+!*    SUBROUTINE BRACWSA()
       SUBROUTINE BRACWSA(XA,XB,N,RADIUS,TAIR,LPPWVINP,FLAGH,FLAGL,
      +           NROOT)
 !*****************************************************************************
 !* Bracket de KEEQ
-!* From Numerical Recipes     
-!* Adapted for VenusGCM A. Stolzenbach 07/2014 
-     
+!* From Numerical Recipes
+!* Adapted for VenusGCM A. Stolzenbach 07/2014
+
       IMPLICIT NONE
 
@@ -433,36 +477,36 @@
 !     External functions needed:
       REAl KEEQ
-!---------------------------------------------------------------------------- 
-
-      REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP 
+!----------------------------------------------------------------------------
+
+      REAL, INTENT(IN) :: RADIUS,TAIR,LPPWVINP
       INTEGER, INTENT(IN) :: N
-      
+
       REAL, INTENT(INOUT) :: XA,XB
-      
+
       INTEGER, INTENT(OUT) ::  NROOT
 
       INTEGER :: I, J
-      
+
       REAL :: DX, FP, FC, X
-      
+
       LOGICAL, INTENT(OUT) :: FLAGH,FLAGL
-    
-      
+
+
       FLAGL=.FALSE.
-      FLAGH=.FALSE.    
+      FLAGH=.FALSE.
       NROOT=0
       DX=(XB-XA)/N
-      X=XA
+      X=XB
       FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
-      
-      DO I=1,N
-         X=X+DX
+
+      DO I=N,1,-1
+         X=X-DX
          FC=KEEQ(RADIUS,TAIR,X,LPPWVINP)
-         
+
          IF ((FP*FC).LE.0.) THEN
             NROOT=NROOT+1
-            XA=X-DX
-            XB=X
-!            RETURN
+            XA=X
+            XB=X+DX
+            RETURN
 !            IF (NROOT.GT.1) THEN
 !               PRINT*,'On a plus d1 intervalle KEEQ=0'
@@ -479,8 +523,8 @@
 !            ENDIF
          ENDIF
-         
+
          FP=FC
       ENDDO
-      
+
       IF (NROOT.EQ.0) THEN
 !         PRINT*,'On a 0 intervalle KEEQ=0'
@@ -492,5 +536,5 @@
 !         PRINT*,'NBRAC',N
 !         STOP
-         
+
 !         X=XA
 !         FP=KEEQ(RADIUS,TAIR,X,LPPWVINP)
@@ -503,5 +547,5 @@
 
 
-!        Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]        
+!        Test determine la tendance globale KEEQ sur [WSAMIN,WSAMAX]
          IF ((ABS(KEEQ(RADIUS,TAIR,XA,LPPWVINP))-
      &    ABS(KEEQ(RADIUS,TAIR,XB,LPPWVINP))).GT.0.0) FLAGH=.TRUE.
@@ -511,10 +555,10 @@
 !        STOP
       ENDIF
-      
+
       END SUBROUTINE BRACWSA
-          
-           
+
+
 !*****************************************************************************
-!*     REAL FUNCTION WVCOND()                                         
+!*     REAL FUNCTION WVCOND()
       REAL FUNCTION WVCOND(WSA,T,P,SAt)
 !*****************************************************************************
@@ -527,9 +571,9 @@
 !     T: temperature (K)
 !     P: pressure (Pa)
-!     OUTPUT: 
+!     OUTPUT:
 !	WVCOND : VMR H2O condense
 
 !      USE chemparam_mod
-      
+
       IMPLICIT NONE
 
@@ -545,10 +589,10 @@
       REAL  NH2SO4
       REAL  H2OCOND, H2SO4COND
-   
+
 
       CONCM= (P)/(1.3806488E-23*T) !air number density, molec/m3? CHECK UNITS!
 
 	NH2SO4=SAt*CONCM
-      
+
       pstand=1.01325E+5 !Pa  1 atm pressure
 
@@ -565,9 +609,9 @@
 
 !acid sat.vap.PP over mixture (flat surface):
-        satpacid=act(2)*acidps ! Pa 
+        satpacid=act(2)*acidps ! Pa
 
 !       Conversion from Pa to N.D #/m3
         DND2=satpacid/(1.3806488E-23*T)
-	        
+
 !	H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
 	IF (NH2SO4.GT.DND2) THEN
@@ -576,5 +620,5 @@
 	H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
 
-!	Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30 
+!	Si on a H2SO4<H2SO4sat on ne condense rien, VMR = 1.0E-30
 	ELSE
 	H2OCOND=1.0E-30*CONCM
@@ -582,26 +626,26 @@
 
 !*****************************************************
-!     ATTENTION: Ici on ne prends pas en compte 
+!     ATTENTION: Ici on ne prends pas en compte
 !                si H2O en defaut!
 !                On veut la situation theorique
 !                a l'equilibre
-!*****************************************************	     	
+!*****************************************************
 !	Test si H2O en defaut H2Ocond>H2O dispo
 !	IF ((H2OCOND.GT.NH2O).AND.(NH2SO4.GE.DND2)) THEN
-	
+
 !	On peut alors condenser tout le H2O dispo
 !	H2OCOND=NH2O
 !	On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
 !	H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
-	    
+
 !      END IF
 
-!     Calcul de H2O condense VMR          
+!     Calcul de H2O condense VMR
       WVCOND=H2OCOND/CONCM
-      
+
       END FUNCTION WVCOND
 
 !*****************************************************************************
-!*     REAL FUNCTION IRFRMWV()                                         
+!*     REAL FUNCTION IRFRMWV()
       REAL	FUNCTION IRFRMWV(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,
      + WVTOT,SATOT,NROOT)
@@ -616,9 +660,9 @@
 
       IMPLICIT NONE
-      
+
       REAL, INTENT(IN) ::    X1, X2
-      REAL, INTENT(IN) ::    XACC 
+      REAL, INTENT(IN) ::    XACC
       INTEGER, INTENT(IN) :: MAXIT,NROOT
-      
+
 !     LOCAL VARIABLES
       REAL :: XL, XH, XM, XNEW, X
@@ -627,5 +671,5 @@
       REAL :: ANS, S, FSIGN
       INTEGER i
-      
+
 !     External variables needed:
       REAL, INTENT(IN) :: TAIR,PAIR
@@ -633,14 +677,14 @@
       REAL, INTENT(IN) :: RADIUS
 
-      
+
 !     Initialisation
       X=X1
-      CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS) 
+      CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
       FL=WVEQ
       X=X2
       CALL ITERWV(X,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,TAIR,PAIR,RADIUS)
       FH=WVEQ
-      
-!     Test Bracketed values 
+
+!     Test Bracketed values
       IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.
      &   ((FL.GT.0.).AND.(FH.LT.0.)))
@@ -649,5 +693,5 @@
          XH=X2
          ANS=-9.99e99
-         
+
          DO i=1, MAXIT
             XM=0.5*(XL+XH)
@@ -656,10 +700,10 @@
             FM=WVEQ
             S=SQRT(FM*FM-FL*FH)
-            
+
             IF (S.EQ.0.0) THEN
                IRFRMWV=WSALOC
                RETURN
-            ENDIF    
-             
+            ENDIF
+
             IF (FL.GT.FH) THEN
                FSIGN=1.0
@@ -667,22 +711,22 @@
                FSIGN=-1.0
             ENDIF
-            
-            XNEW=XM+(XM-XL)*(FSIGN*FM/S) 
-            
+
+            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
+
             IF (ABS(XNEW-ANS).LE.XACC) THEN
                IRFRMWV=WSALOC
                RETURN
-            ENDIF    
-            
+            ENDIF
+
             ANS=XNEW
             CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
      &       TAIR,PAIR,RADIUS)
             FNEW=WVEQ
-            
+
             IF (FNEW.EQ.0.0) THEN
                IRFRMWV=WSALOC
                RETURN
             ENDIF
-            
+
             IF (SIGN(FM, FNEW).NE.FM) THEN
                XL=XM
@@ -708,5 +752,5 @@
          XH=X2
          ANS=-9.99e99
-         
+
          DO i=1, MAXIT
             XM=0.5*(XL+XH)
@@ -720,7 +764,7 @@
                FSIGN=-1.0
             ENDIF
-            
-            XNEW=XM+(XM-XL)*(FSIGN*FM/S)     
-            
+
+            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
+
             ANS=XNEW
             CALL ITERWV(ANS,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
@@ -753,5 +797,5 @@
          PRINT*,'NROOT de BRACWV', NROOT
          IF (ABS(FL).LT.XACC) THEN
-         PRINT*,'IRFRMWV FL == 0',FL 
+         PRINT*,'IRFRMWV FL == 0',FL
          PRINT*,'X1',X1,'X2',X2,'FH',FH
            CALL ITERWV(X1,WVLIQ,WVEQ,WVTOT,WSALOC,SATOT,
@@ -768,25 +812,25 @@
          RETURN
          ENDIF
-         IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN 
+         IF ((ABS(FL).GT.XACC).AND.(ABS(FH).GT.XACC)) THEN
            PRINT*,'STOP dans IRFRMWV avec rien == 0'
            PRINT*,'X1',X1,'X2',X2
            PRINT*,'Fcalc',FL,FH
            PRINT*,'T',TAIR,'P',PAIR,'R',RADIUS
-           STOP   
+           STOP
          ENDIF
-         IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN 
+         IF ((ABS(FL).LT.XACC).AND.(ABS(FH).LT.XACC)) THEN
            PRINT*,'STOP dans IRFRMWV Trop de solution < WVACC'
            PRINT*,FL,FH
-           STOP   
+           STOP
          ENDIF
-         
-         
+
+
       ENDIF
 !  FIN Test Bracketed values
-       
+
       END FUNCTION IRFRMWV
-                 
+
 !*****************************************************************************
-!*     REAL FUNCTION IRFRMSA()                                         
+!*     REAL FUNCTION IRFRMSA()
       REAL	FUNCTION IRFRMSA(X1,X2,XACC,MAXIT,RADIUS,TAIR,PAIR,LPPWV,
      +               NB)
@@ -801,9 +845,9 @@
 
       IMPLICIT NONE
-      
+
       REAL, INTENT(IN) ::    X1, X2
-      REAL, INTENT(IN) ::    XACC 
+      REAL, INTENT(IN) ::    XACC
       INTEGER, INTENT(IN) :: MAXIT, NB
-      
+
 !     LOCAL VARIABLES
       REAL XL, XH, XM, XNEW
@@ -811,20 +855,20 @@
       REAL ANS, S, FSIGN
       INTEGER i
-      
+
 !     External variables needed:
       REAL, INTENT(IN) :: TAIR,PAIR
       REAL, INTENT(IN) :: LPPWV
       REAL, INTENT(IN) :: RADIUS
-      
+
 !     External functions needed:
       REAL KEEQ
 
 
-      
+
 !     Initialisation
       FL=KEEQ(RADIUS,TAIR,X1,LPPWV)
       FH=KEEQ(RADIUS,TAIR,X2,LPPWV)
-      
-!     Test Bracketed values 
+
+!     Test Bracketed values
       IF (((FL.LT.0.).AND.(FH.GT.0.)).OR.((FL.GT.0.).AND.(FH.LT.0.)))
      &  THEN
@@ -832,15 +876,15 @@
          XH=X2
          ANS=-9.99e99
-         
+
          DO i=1, MAXIT
             XM=0.5*(XL+XH)
             FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
             S=SQRT(FM*FM-FL*FH)
-            
+
             IF (S.EQ.0.0) THEN
                IRFRMSA=ANS
                RETURN
-            ENDIF    
-             
+            ENDIF
+
             IF (FL.GT.FH) THEN
                FSIGN=1.0
@@ -848,20 +892,20 @@
                FSIGN=-1.0
             ENDIF
-            
-            XNEW=XM+(XM-XL)*(FSIGN*FM/S) 
-            
+
+            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
+
             IF (ABS(XNEW-ANS).LE.XACC) THEN
                IRFRMSA=ANS
                RETURN
-            ENDIF    
-            
+            ENDIF
+
             ANS=XNEW
             FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
-            
+
             IF (FNEW.EQ.0.0) THEN
                IRFRMSA=ANS
                RETURN
             ENDIF
-            
+
             IF (SIGN(FM, FNEW).NE.FM) THEN
                XL=XM
@@ -888,10 +932,10 @@
          PRINT*,'Borne XL',XL,'XH',XH
          ANS=-9.99e99
-         
+
          DO i=1, MAXIT
             XM=0.5*(XL+XH)
             FM=KEEQ(RADIUS,TAIR,XM,LPPWV)
             S=SQRT(FM*FM-FL*FH)
- 
+
             IF (FL.GT.FH) THEN
                FSIGN=1.0
@@ -899,7 +943,7 @@
                FSIGN=-1.0
             ENDIF
-            
-            XNEW=XM+(XM-XL)*(FSIGN*FM/S)  
-            
+
+            XNEW=XM+(XM-XL)*(FSIGN*FM/S)
+
             ANS=XNEW
             FNEW=KEEQ(RADIUS,TAIR,ANS,LPPWV)
@@ -937,5 +981,5 @@
          IF ((FL.NE.0.).AND.(FH.NE.0.)) THEN
            PRINT*,'IRFRMSA FH and FL neq 0: ', FL, FH
-           PRINT*,'X1',X1,'X2',X2 
+           PRINT*,'X1',X1,'X2',X2
            PRINT*,'Kind F', KIND(FL), KIND(FH)
            PRINT*,'Kind X', KIND(X1), KIND(X2)
@@ -944,13 +988,13 @@
            PRINT*,'nb root BRACWSA',NB
            STOP
-         ENDIF   
-      
+         ENDIF
+
       ENDIF
 !  FIN Test Bracketed values
-       
+
       END function IRFRMSA
-      
+
 !*****************************************************************************
-!*     REAL FUNCTION KEEQ()                                         
+!*     REAL FUNCTION KEEQ()
       REAL	FUNCTION KEEQ(RADIUS,TAIR,WSA,LPPWV)
 !*****************************************************************************
@@ -985,11 +1029,11 @@
      +        C1=2.0d0*MH2O/RGAS)
 
-      
+
       KEEQ=LPPWV-C1*SIGMADROPLET(WSA,TAIR)/
      &     (TAIR*RADIUS*RHODROPLET(WSA,TAIR))-
      &     PWVSAS_GV(TAIR,WSA)
-      
+
       END FUNCTION KEEQ
-      
+
 *****************************************************************************
 *     REAL FUNCTION PWVSAS_GV(TAIR,WSA)                                        
@@ -1048,9 +1092,39 @@
      +     K2=K1*K1/2.0)
 !
-!
-      CP=CPH2O(WSA)
-      F=-FFH2O(WSA)
-      L=-LH2O(WSA)
-      ALFA=ALH2O(WSA)
+      INTEGER :: KLO,KHI,K,I,NPOINT
+      PARAMETER(NPOINT=110)
+      REAL, DIMENSION(NPOINT) :: WTAB(NPOINT)
+      DATA (WTAB(I),I=1,NPOINT)/
+     +0.00000,0.08932,0.09819,0.10792,0.11980,0.13461,0.15360,0.16525,
+     +0.17882,0.19482,0.21397,0.23728,0.26629,0.27999,0.29517,0.31209,
+     +0.33107,0.35251,0.36430,0.37691,0.39043,0.40495,0.42059,0.43749,
+     +0.44646,0.45580,0.46555,0.47572,0.48634,0.49745,0.50908,0.52126,
+     +0.53405,0.54747,0.56159,0.57646,0.58263,0.58893,0.59537,0.60195,
+     +0.60868,0.61557,0.62261,0.62981,0.63718,0.64472,0.65245,0.66037,
+     +0.66847,0.67678,0.68530,0.69404,0.70300,0.71220,0.72164,0.73133,
+     +0.73628,0.74129,0.74637,0.75152,0.75675,0.76204,0.76741,0.77286,
+     +0.77839,0.78399,0.78968,0.79545,0.80130,0.80724,0.81327,0.81939,
+     +0.82560,0.83191,0.83832,0.84482,0.85143,0.85814,0.86495,0.87188,
+     +0.87892,0.88607,0.89334,0.90073,0.90824,0.91588,0.92365,0.93156,
+     +0.93959,0.94777,0.95610,0.96457,0.97319,0.98196,0.99090,0.99270,
+     +0.99452,0.99634,0.99725,0.99817,0.99835,0.99853,0.99872,0.99890,
+     +0.99908,0.99927,0.99945,0.99963,0.99982,1.0000/
+!
+      KLO=1
+      KHI=NPOINT
+ 1    IF(KHI-KLO.GT.1) THEN
+          K=(KHI+KLO)/2
+          IF(WTAB(K).GT.MAX(WTAB(1),WSA)) THEN
+              KHI=K
+          ELSE
+              KLO=K
+          ENDIF
+          GOTO 1
+      ENDIF
+!
+      CP=CPH2O(WSA,KHI,KLO)
+      F=-FFH2O(WSA,KHI,KLO)
+      L=-LH2O(WSA,KHI,KLO)
+      ALFA=ALH2O(WSA,KHI,KLO)
 !
       A=ADOT+(CP-K1*ALFA)/RGAS
@@ -1070,5 +1144,5 @@
 *******************************************************************************
 *     REAL FUNCTION CPH2O(W)
-      REAL FUNCTION CPH2O(W)
+      REAL FUNCTION CPH2O(W,khi_in,klo_in)
 *******************************************************************************
 *
@@ -1086,4 +1160,6 @@
      +              Y2(NPOINT),YWORK(NPOINT)
       REAL, INTENT(IN):: W
+      INTEGER, INTENT(IN):: khi_in,klo_in
+      INTEGER :: khi,klo
       REAL :: CPH
       LOGICAL :: FIRST
@@ -1125,5 +1201,15 @@
           CALL SPLINE(WTAB,CPHTAB,NPOINT,YWORK,Y2)
       ENDIF
-      CALL SPLINT(WTAB,CPHTAB,Y2,NPOINT,W,CPH)
+
+      if(khi_in.GT.NPOINT) then
+         khi=NPOINT
+         klo=NPOINT-1
+      else
+         khi=khi_in
+         klo=klo_in
+      endif
+
+      CALL SPLINT(WTAB(khi),WTAB(klo),CPHTAB(khi),CPHTAB(klo),
+     .              Y2(khi),Y2(klo),W,CPH)
       CPH2O=CPH
       
@@ -1131,5 +1217,5 @@
 !
 *******************************************************************************
-      REAL FUNCTION FFH2O(W)
+      REAL FUNCTION FFH2O(W,khi,klo)
 *     REAL FUNCTION FFH2O(W)
 *******************************************************************************
@@ -1147,4 +1233,5 @@
       REAL, DIMENSION(NPOINT) :: WTAB,FFTAB,Y2,YWORK
       REAL, INTENT(IN) :: W
+      INTEGER, INTENT(IN):: khi,klo
       REAL :: FF
       LOGICAL :: FIRST
@@ -1186,5 +1273,7 @@
           CALL SPLINE(WTAB,FFTAB,NPOINT,YWORK,Y2)
       ENDIF
-      CALL SPLINT(WTAB,FFTAB,Y2,NPOINT,W,FF)
+
+      CALL SPLINT(WTAB(khi),WTAB(klo),FFTAB(khi),FFTAB(klo),
+     .              Y2(khi),Y2(klo),W,FF)
       FFH2O=FF
       
@@ -1192,5 +1281,5 @@
 !
 *******************************************************************************
-      REAL FUNCTION LH2O(W)
+      REAL FUNCTION LH2O(W,khi,klo)
 *     REAL FUNCTION LH2O(W)
 *******************************************************************************
@@ -1208,4 +1297,5 @@
       REAL, DIMENSION(NPOINT) ::  WTAB,LTAB,Y2,YWORK
       REAL, INTENT(IN) :: W
+      INTEGER, INTENT(IN):: khi,klo
       REAL :: L
       LOGICAL :: FIRST
@@ -1247,10 +1337,12 @@
           CALL SPLINE(WTAB,LTAB,NPOINT,YWORK,Y2)
       ENDIF
-      CALL SPLINT(WTAB,LTAB,Y2,NPOINT,W,L)
+
+      CALL SPLINT(WTAB(khi),WTAB(klo),LTAB(khi),LTAB(klo),
+     .              Y2(khi),Y2(klo),W,L)
       LH2O=L
       
       END FUNCTION LH2O
 *******************************************************************************
-      REAL FUNCTION ALH2O(W)
+      REAL FUNCTION ALH2O(W,khi_in,klo_in)
 *     REAL FUNCTION ALH2O(W)
 *******************************************************************************
@@ -1268,4 +1360,6 @@
       REAL, DIMENSION(NPOINT) :: WTAB,ATAB,Y2,YWORK
       REAL, INTENT(IN) :: W
+      INTEGER, INTENT(IN):: khi_in,klo_in
+      INTEGER :: khi,klo
       REAL :: A
       LOGICAL :: FIRST
@@ -1304,5 +1398,15 @@
           CALL SPLINE(WTAB,ATAB,NPOINT,YWORK,Y2)
       ENDIF
-      CALL SPLINT(WTAB,ATAB,Y2,NPOINT,MAX(WTAB(1),W),A)
+
+      if(klo_in.LT.15) then
+         khi=2
+         klo=1
+      else
+         khi=khi_in-14
+         klo=klo_in-14
+      endif
+
+      CALL SPLINT(WTAB(khi),WTAB(klo),ATAB(khi),ATAB(klo),
+     .              Y2(khi),Y2(klo),W,A)
       ALH2O=A
       
@@ -1358,5 +1462,5 @@
 
 !******************************************************************************
-      SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
+      SUBROUTINE SPLINT(XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo,X,Y)
 !******************************************************************************
 !     Cubic spline calculation
@@ -1364,27 +1468,14 @@
       IMPLICIT NONE
 
-      INTEGER, INTENT(IN) :: N
-      INTEGER :: KLO,KHI,K
-      REAL, INTENT(IN), DIMENSION(N) :: XA,YA,Y2A
+      REAL, INTENT(IN) :: XAhi,XAlo,YAhi,YAlo,Y2Ahi,Y2Alo
       REAL, INTENT(IN) :: X
       REAL, INTENT(OUT) :: Y
       REAL :: H,A,B
 !
-      KLO=1
-      KHI=N
- 1    IF(KHI-KLO.GT.1) THEN
-          K=(KHI+KLO)/2
-          IF(XA(K).GT.X) THEN
-              KHI=K
-          ELSE
-              KLO=K
-          ENDIF
-          GOTO 1
-      ENDIF
-      H=XA(KHI)-XA(KLO)
-      A=(XA(KHI)-X)/H
-      B=(X-XA(KLO))/H
-      Y=A*YA(KLO)+B*YA(KHI)+
-     +        ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0d0
+      H=XAhi-XAlo
+      A=(XAhi-X)/H
+      B=(X-XAlo)/H
+      Y=A*YAlo+B*YAhi+
+     +        ((A**3-A)*Y2Alo+(B**3-B)*Y2Ahi)*(H**2)/6.0d0
 !
 
@@ -1394,6 +1485,6 @@
      + T,H2SO4COND,H2OCOND,RMTOT)
 
-!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION) 
-!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL 
+!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
+!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
 !                                       SIZE DISTRIBTUION
 !     ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
@@ -1409,12 +1500,12 @@
 !     T: temperature (K)
 !
-!     OUTPUT: 
+!     OUTPUT:
 !     RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
-!            mais rho_droplet et M_tot_distrib doivent tre de meme dimension
+!            mais rho_droplet et M_tot_distrib doivent etre de meme dimension
 !	H2OCOND
 !	H2SO4COND
 
 
-      
+
       IMPLICIT NONE
 
@@ -1430,5 +1521,5 @@
 !     masse of an H2SO4 molecule (kg)
       RMH2S4=98.078/(6.02214129E+26)
-      
+
       pstand=1.01325E+5 !Pa  1 atm pressure
 
@@ -1445,5 +1536,5 @@
 
 !acid sat.vap.PP over mixture (flat surface):
-        satpacid=act(2)*acidps ! Pa 
+        satpacid=act(2)*acidps ! Pa
 
 !       Conversion from Pa to N.D #/m3
@@ -1451,5 +1542,5 @@
 !       Conversion from N.D #/m3 TO #/cm3
 !        DND2=DND2*1.d-6
-	        
+
 !	H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
 	IF (H2SO4.GE.DND2) THEN
@@ -1458,7 +1549,8 @@
 	H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
 
-!     RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3 
+!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
+!     RMTOT = Mtot/rho_droplet
 !	RMTOT=M_distrib/rho_droplet
-	
+
 	RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
 
@@ -1469,5 +1561,5 @@
 	RMTOT=0.0E+0
 	END IF
-	     	
+
 !	Test si H2O en defaut H2Ocond>H2O dispo
 	IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
@@ -1484,19 +1576,19 @@
 !      PRINT*,'WSA',WSA,'RHO',DENSO4
 !      STOP
-      
-	
+
+
 !	On peut alors condenser tout le H2O dispo
 	H2OCOND=H2O
 !	On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
 	H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
-	
-!     RMTOT: = Mass of H2SO4 satur per cm3 of air/ Mass of sulfuric acid part of droplet solution per cm3 
-!	RMTOT=Volume of aerosol cm3 /cm3 of air
-!	Volume of aerosol/cm3 air
-	
+
+!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
+!	RMTOT=Volume of aerosol m3 /m3 of air
+!	Volume of aerosol/m3 air
+
 	RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
-     
+
       END IF
-      		
+
       END SUBROUTINE CALCM_SAT
 
@@ -1510,5 +1602,5 @@
   !     mole fraction 0,...,1
   !     J. Phys. Chem. Ref. Data, Vol. 20, No. 6,PP.1157, 1991
-  !+++++++++++++++++++++++++++++++++++++++++++++++++++ 
+  !+++++++++++++++++++++++++++++++++++++++++++++++++++
 
          IMPLICIT NONE
@@ -1523,5 +1615,5 @@
 !         write(*,*) 'x, T ', x, T
 
-         act(2)=activitya(x,T) 
+         act(2)=activitya(x,T)
          act(1)=activityw(x,T)
 
@@ -1535,7 +1627,7 @@
 
        REAL, INTENT(IN) :: T
-       m111=-23.524503387D0 
-     &    +0.0406889449841D0*T 
-     &    -0.151369362907D-4*T**2+2961.44445015D0/T 
+       m111=-23.524503387D0
+     &    +0.0406889449841D0*T
+     &    -0.151369362907D-4*T**2+2961.44445015D0/T
      &    +0.492476973663D0*dlog(T)
        END FUNCTION m111
@@ -1544,6 +1636,6 @@
 
        REAL, INTENT(IN) :: T
-       m121=1114.58541077D0-1.1833078936D0*T 
-     &    -0.00209946114412D0*T**2-246749.842271D0/T 
+       m121=1114.58541077D0-1.1833078936D0*T
+     &    -0.00209946114412D0*T**2-246749.842271D0/T
      &    +34.1234558134D0*dlog(T)
        END FUNCTION m121
@@ -1552,6 +1644,6 @@
 
        REAL, INTENT(IN) :: T
-       m221=-80.1488100747D0-0.0116246143257D0*T 
-     &    +0.606767928954D-5*T**2+3092.72150882D0/T 
+       m221=-80.1488100747D0-0.0116246143257D0*T
+     &    +0.606767928954D-5*T**2+3092.72150882D0/T
      &    +12.7601667471D0*dlog(T)
        END FUNCTION m221
@@ -1560,6 +1652,6 @@
 
        REAL, INTENT(IN) :: T
-       m122=888.711613784D0-2.50531359687D0*T 
-     &    +0.000605638824061D0*T**2-196985.296431D0/T 
+       m122=888.711613784D0-2.50531359687D0*T
+     &    +0.000605638824061D0*T**2-196985.296431D0/T
      &    +74.550064338D0*dlog(T)
        END FUNCTION m122
@@ -1568,6 +1660,6 @@
 
        REAL, INTENT(IN) :: T
-       e111=2887.31663295D0-3.32602457749D0*T 
-     &    -0.2820472833D-2*T**2-528216.112353D0/T 
+       e111=2887.31663295D0-3.32602457749D0*T
+     &    -0.2820472833D-2*T**2-528216.112353D0/T
      &    +0.68699743564D0*dlog(T)
        END FUNCTION e111
@@ -1576,6 +1668,6 @@
 
        REAL, INTENT(IN) :: T
-       e121=-370.944593249D0-0.690310834523D0*T 
-     &    +0.56345508422D-3*T**2-3822.52997064D0/T 
+       e121=-370.944593249D0-0.690310834523D0*T
+     &    +0.56345508422D-3*T**2-3822.52997064D0/T
      &    +94.2682037574D0*dlog(T)
        END FUNCTION e121
@@ -1584,6 +1676,6 @@
 
        REAL, INTENT(IN) :: T
-       e211=38.3025318809D0-0.0295997878789D0*T 
-     &    +0.120999746782D-4*T**2-3246.97498999D0/T 
+       e211=38.3025318809D0-0.0295997878789D0*T
+     &    +0.120999746782D-4*T**2-3246.97498999D0/T
      &    -3.83566039532D0*dlog(T)
        END FUNCTION e211
@@ -1592,6 +1684,6 @@
 
        REAL, INTENT(IN) :: T
-       e221=2324.76399402D0-0.141626921317D0*T 
-     &    -0.00626760562881D0*T**2-450590.687961D0/T 
+       e221=2324.76399402D0-0.141626921317D0*T
+     &    -0.00626760562881D0*T**2-450590.687961D0/T
      &    -61.2339472744D0*dlog(T)
        END FUNCTION e221
@@ -1600,6 +1692,6 @@
 
        REAL, INTENT(IN) :: T
-       e122=-1633.85547832D0-3.35344369968D0*T 
-     &    +0.00710978119903D0*T**2+198200.003569D0/T 
+       e122=-1633.85547832D0-3.35344369968D0*T
+     &    +0.00710978119903D0*T**2+198200.003569D0/T
      &    +246.693619189D0*dlog(T)
        END FUNCTION e122
@@ -1608,6 +1700,6 @@
 
        REAL, INTENT(IN) :: T
-       e212=1273.75159848D0+1.03333898148D0*T 
-     &    +0.00341400487633D0*T**2+195290.667051D0/T 
+       e212=1273.75159848D0+1.03333898148D0*T
+     &    +0.00341400487633D0*T**2+195290.667051D0/T
      &    -431.737442782D0*dlog(T)
        END FUNCTION e212
@@ -1616,22 +1708,22 @@
 
        REAL, INTENT(IN) :: T,x1
-       REAL :: 
-     &          m111,m121,m221,m122 
+       REAL ::
+     &          m111,m121,m221,m122
      &            ,e111,e121,e211,e122,e212,e221
-       lnAa=-( 
-     &    (2*m111(T)+e111(T)*(2*dlog(x1)+1))*x1 
-     &    +(2*m121(T)+e211(T)*dlog(1-x1)+e121(T)*(dlog(x1)+1))*(1-x1) 
-     &    -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 
-     &    -(2*m121(T)+e121(T)*(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1) 
-     &    -(2*m122(T)+e122(T)*dlog(x1) 
-     &           +e212(T)*dlog(1-x1))*(1-x1))*x1*(1-x1) 
-     &    -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 
-     &    -x1*(1-x1)*( 
-     &                  (6*m122(T)+e122(T)*(3*dlog(x1)+1) 
-     &                          +e212(T)*(3*dlog(1-x1)+1) 
-     &                   )*x1*(1-x1) 
-     &                -(2*m122(T)+e122(T)*(dlog(x1)+1) 
-     &                                   +e212(T)*dlog(1-x1) 
-     &                    )*(1-x1)) 
+       lnAa=-(
+     &    (2*m111(T)+e111(T)*(2*dlog(x1)+1))*x1
+     &    +(2*m121(T)+e211(T)*dlog(1-x1)+e121(T)*(dlog(x1)+1))*(1-x1)
+     &    -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1
+     &    -(2*m121(T)+e121(T)*(dlog(x1)+1)+e211(T)*(dlog(1-x1)+1)
+     &    -(2*m122(T)+e122(T)*dlog(x1)
+     &           +e212(T)*dlog(1-x1))*(1-x1))*x1*(1-x1)
+     &    -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2
+     &    -x1*(1-x1)*(
+     &                  (6*m122(T)+e122(T)*(3*dlog(x1)+1)
+     &                          +e212(T)*(3*dlog(1-x1)+1)
+     &                   )*x1*(1-x1)
+     &                -(2*m122(T)+e122(T)*(dlog(x1)+1)
+     &                                   +e212(T)*dlog(1-x1)
+     &                    )*(1-x1))
      &     )
        END FUNCTION lnAa
@@ -1640,19 +1732,19 @@
 
        REAL, INTENT(IN) :: T,x1
-       REAL :: 
-     &          m111,m121,m221,m122 
+       REAL ::
+     &          m111,m121,m221,m122
      &            ,e111,e121,e211,e122,e212,e221
-       lnAw=-( 
-     &  (2*m121(T)+e121(T)*dlog(x1)+e211(T)*(dlog(1-x1)+1))*x1 
-     &  +(2*m221(T)+e221(T)*(2*dlog(1-x1)+1))*(1-x1) 
-     &  -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1 
-     & -(2*m121(T)+e121(T)*(dlog(x1)+1) 
-     &            +e211(T)*(dlog(1-x1)+1))*x1*(1-x1) 
-     &        -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2 
-     &   +x1*(2*m122(T)+e122(T)*dlog(x1)+e212(T)*dlog(1-x1))*x1*(1-x1) 
-     &  +x1*(1-x1)*((2*m122(T)+e122(T)*dlog(x1) 
-     &                        +e212(T)*(dlog(1-x1)+1))*x1 
-     &               -(6*m122(T)+e122(T)*(3*dlog(x1)+1) 
-     &                        +e212(T)*(3*dlog(1-x1)+1))*(1-x1)*x1) 
+       lnAw=-(
+     &  (2*m121(T)+e121(T)*dlog(x1)+e211(T)*(dlog(1-x1)+1))*x1
+     &  +(2*m221(T)+e221(T)*(2*dlog(1-x1)+1))*(1-x1)
+     &  -(m111(T)+e111(T)*(dlog(x1)+1))*x1*x1
+     & -(2*m121(T)+e121(T)*(dlog(x1)+1)
+     &            +e211(T)*(dlog(1-x1)+1))*x1*(1-x1)
+     &        -(m221(T)+e221(T)*(dlog(1-x1)+1))*(1-x1)**2
+     &   +x1*(2*m122(T)+e122(T)*dlog(x1)+e212(T)*dlog(1-x1))*x1*(1-x1)
+     &  +x1*(1-x1)*((2*m122(T)+e122(T)*dlog(x1)
+     &                        +e212(T)*(dlog(1-x1)+1))*x1
+     &               -(6*m122(T)+e122(T)*(3*dlog(x1)+1)
+     &                        +e212(T)*(3*dlog(1-x1)+1))*(1-x1)*x1)
      &     )
        END FUNCTION lnAw
@@ -1671,5 +1763,5 @@
        FUNCTION activityw(xal,T)
 
-       REAL, INTENT(IN) :: T,xal 
+       REAL, INTENT(IN) :: T,xal
        REAL :: lnAw
 
@@ -1702,5 +1794,5 @@
       a=0.11864+xmass*(-0.11651+xmass*(0.76852+xmass*
      & (-2.40909+xmass*(2.95434-xmass*1.25852))))
-      b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass* 
+      b=-0.00015709+xmass*(0.00040102+xmass*(-0.00239950+xmass*
      & (0.007611235+xmass*(-0.00937386+xmass*0.00389722))))
       SIGMADROPLET=a+t*b
@@ -1745,5 +1837,5 @@
 
       a=0.7681724+xmass*(2.1847140+xmass*(7.1630022+xmass*
-     & (-44.31447+xmass* 
+     & (-44.31447+xmass*
      & (88.75606+xmass*(-75.73729+xmass*23.43228)))))
       b=1.808225e-3+xmass*(-9.294656e-3+xmass*(-0.03742148+
@@ -1757,6 +1849,2 @@
       RETURN
       END FUNCTION RHODROPLET
-
-
-
-
Index: /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1638)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1639)
@@ -25,4 +25,6 @@
 c With extension NLTE (G. Gilli, 2014)
 
+c Ksi matrices latitudinaly interpolated  (I. Garate-Lopez, 2016)
+
 c======================================================================
       use dimphy
@@ -74,9 +76,15 @@
       integer mat0,lat,ips,isza,ips0,isza0
       real    factp,factz,ksi
-
+c ------- for lat-interp ----------------
+	integer mat0A, mat0B, latA, latB, kasua
+	integer	ipsA, ipsB, iszaA, iszaB, ips0A, ips0B, isza0A, isza0B
+	real	lat_deg, latA_deg, latB_deg
+	real	factlat, k1, k2, k3, k4
+c --------------------------------------
       logical firstcall
       data    firstcall/.true./
       save    firstcall
       
+cERROR		! For checking if the file it's being read
 c-------------------------------------------
 c  Initialisations
@@ -99,4 +107,5 @@
       ENDDO
 
+
 c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
       DO j = 1, klon
@@ -110,8 +119,14 @@
         zcool(k) = 0.0
        ENDDO
+c       zheat(1:klev)=0.0	!Explicit loop (no change in performance)
+c       zcool(1:klev)=0.0
+
        DO k = 1, klev+1
         ZFLNET(k) = 0.0
         ZFSNET(k) = 0.0
        ENDDO
+c       ZFLNET(1:klev+1)=0.0
+c       ZFSNET(1:klev+1)=0.0
+
        ztopsw = 0.0
        ztoplw = 0.0
@@ -120,7 +135,7 @@
        zsollwdown = 0.0
       
-         zfract = fract(j)
-         zrmu0 = rmu0(j)
- 
+       zfract = fract(j)
+       zrmu0 = rmu0(j)
+
        DO k = 1, klev+1
          PPB(k) = paprs(j,k)/1.e5
@@ -132,5 +147,5 @@
        ENDDO
        pt0(klev+1) = 0.
- 
+       
        DO k = 0,klev+1
        DO i = 0,klev+1
@@ -139,5 +154,5 @@
        ENDDO
        ENDDO
-       
+
 c======================================================================
 c Getting psi and deltapsi
@@ -160,26 +175,94 @@
 c finding the right matrixes
 c --------------------------
-       mat0 = 0
-       if (nlatve.eq.1) then   ! clouds are taken as uniform
-         lat=1
-       else
-         if (abs(latitude_deg(j)).le.50.) then
-	   lat=1
-	 elseif (abs(latitude_deg(j)).le.60.) then
-	   lat=2
-	 elseif (abs(latitude_deg(j)).le.70.) then
-	   lat=3
-	 elseif (abs(latitude_deg(j)).le.80.) then
-	   lat=4
-	 else
-	   lat=5
-	 endif
-       endif
-       
+
+	mat0  = 0
+	mat0A = 0
+	mat0B = 0
+
+c    Latitude
+c    --------
+	lat  = 0
+	latA = 0
+	latB = 0
+
+c	write(*,*) 'nlatve:', nlatve
+
+	lat_deg = abs(latitude_deg(j))
+
+c	if (nlatve.eq.1) then   ! clouds are taken as uniform
+	if ((nlatve.eq.1).or.(lat_deg.le.25.)) then
+	    lat  = 1
+        elseif (lat_deg.le.50.) then
+	    lat  = 1
+	    latA = 1
+	    latB = 2
+	    latA_deg = 25.0
+	    latB_deg = 55.0
+        elseif (lat_deg.le.55.) then
+	    lat  = 2
+	    latA = 1
+	    latB = 2
+	    latA_deg = 25.0
+	    latB_deg = 55.0
+        elseif (lat_deg.le.60.) then
+	    lat  = 2
+	    latA = 2
+	    latB = 3
+	    latA_deg = 55.0
+	    latB_deg = 65.0
+        elseif (lat_deg.le.65.) then
+	    lat  = 3
+	    latA = 2
+	    latB = 3
+	    latA_deg = 55.0
+	    latB_deg = 65.0
+        elseif (lat_deg.le.70.) then
+	    lat  = 3
+	    latA = 3
+	    latB = 4
+	    latA_deg = 65.0
+	    latB_deg = 75.0
+        elseif (lat_deg.le.75.) then
+	    lat  = 4
+	    latA = 3
+	    latB = 4
+	    latA_deg = 65.0
+	    latB_deg = 75.0
+        elseif (lat_deg.le.80.) then
+	    lat  = 4
+	    latA = 4
+	    latB = 5
+	    latA_deg = 75.0
+	    latB_deg = 85.0
+        elseif (lat_deg.le.85.) then
+            lat = 5
+	    latA = 4
+	    latB = 5
+	    latA_deg = 75.0
+	    latB_deg = 85.0
+        else
+            lat = 5
+	endif
+
+c        write(*,*) 'Lat',lat,'LatA',latA,'LatB',latB
+
+	factlat = 0
+	if (latA.gt.0) then
+	  factlat = (lat_deg - latA_deg) / (latB_deg - latA_deg)
+	endif
+
+c       write (*,*) 'Factor de correccion:', factlat
+
+
+c    Pressure at Surface
+c    -------------------    
+
        ips0=0
-       if (nbpsve(lat).gt.1) then
+       ips0A=0
+       ips0B=0
+       if (nbpsve(lat).gt.1) then            ! Interpolation on ps
        do ips=1,nbpsve(lat)-1
          if (  (psurfve(ips,lat).ge.paprs(j,1))
-     .    .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
+     .  .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
               ips0  = ips
 c             print*,'ig=',j,'  ips0=',ips
@@ -192,9 +275,52 @@
         ips0=1
        endif
+
+       if (latA.eq.lat) then
+         ips0A=ips0
+       else
+         if (latA.gt.0) then
+           if (nbpsve(latA).gt.1) then
+	     do ipsA=1,nbpsve(latA)-1
+	       if (  (psurfve(ipsA,latA).ge.paprs(j,1))
+     .        .and.(psurfve(ipsA+1,latA).lt.paprs(j,1)) ) then
+                ips0A = ipsA
+                exit
+	       endif
+	     enddo
+	   else            ! Only one ps, no interpolation
+	     ips0A=1
+	   endif     ! nbpsve(latA).gt.1
+         endif       ! latA.gt.0   (if latA=0 ips0A is not used, so it doesn't matter)
+       endif         ! latA.eq.lat
+
+       if (latB.eq.lat) then
+         ips0B=ips0
+       else
+         if (latB.gt.0) then
+           if (nbpsve(latB).gt.1) then
+	     do ipsB=1,nbpsve(latB)-1
+	       if (  (psurfve(ipsB,latB).ge.paprs(j,1))
+     .        .and.(psurfve(ipsB+1,latB).lt.paprs(j,1)) ) then
+                 ips0B = ipsB
+                 exit
+	       endif
+	     enddo
+	   else
+	     ips0B=1
+	   endif     ! nbpsve(latB).gt.1
+         endif       ! latB.gt.0   (if latB=0 ips0B is not used, so it doesn't matter)
+       endif         ! latB.eq.lat
+
+
+c    Solar Zenith Angle
+c    ------------------
+
        isza0=0
+       isza0A=0
+       isza0B=0
        if (nbszave(lat).gt.1) then
         do isza=1,nbszave(lat)-1
          if (  (szave(isza,lat).ge.zrmu0)
-     .    .and.(szave(isza+1,lat).lt.zrmu0) ) then
+     .  .and.(szave(isza+1,lat).lt.zrmu0) ) then
               isza0  = isza
 c             print*,'ig=',j,'  isza0=',isza
@@ -207,4 +333,44 @@
         isza0=-99
        endif
+
+
+       if (latA.eq.lat) then
+         isza0A=isza0
+       else
+         if (latA.gt.0) then
+           if (nbszave(latA).gt.1) then
+	     do iszaA=1,nbszave(latA)-1
+	       if (  (szave(iszaA,latA).ge.zrmu0)
+     .        .and.(szave(iszaA+1,latA).lt.zrmu0) ) then
+                 isza0A = iszaA
+                 exit
+	       endif
+	     enddo
+	   else            ! Only one sza, no interpolation
+	     isza0A=-99
+	   endif     ! nbszave(latA).gt.1
+         endif       ! latA.gt.0   (if latA=0 isza0A is not used, so it doesn't matter)
+       endif         ! latA.eq.lat
+
+       if (latB.eq.lat) then
+         isza0B=isza0
+       else
+         if (latB.gt.0) then
+           if (nbszave(latB).gt.1) then
+	     do iszaB=1,nbszave(latB)-1
+	       if (  (szave(iszaB,latB).ge.zrmu0)
+     .        .and.(szave(iszaB+1,latB).lt.zrmu0) ) then
+                 isza0B = iszaB
+                 exit
+	       endif
+	     enddo
+	   else            ! Only one sza, no interpolation
+	     isza0B=-99
+	   endif     ! nbszave(latB).gt.1
+         endif       ! latB.gt.0   (if latB=0 isza0B is not used, so it doesn't matter)
+       endif         ! latB.eq.lat
+
+c        write(*,*) 'nbszave', nbszave(lat),'nbpsve(lat)',nbpsve(lat)
+
        
 c -------- Probleme aux bords
@@ -215,6 +381,6 @@
      .               /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
        endif
-       if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1))) 
-     .   then
+       if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1)))
+     . then
               ips0  = nbpsve(lat)-1
               print*,'Extrapolation! ig=',j,'  ips0=',ips0
@@ -230,32 +396,136 @@
          stop
        endif
-       
+
+
+
        if (isza0.eq.-99) then
-         mat0 = indexve(lat)+ips0
+         mat0  = indexve(lat) +ips0
+         mat0A = indexve(latA)+ips0A
+         mat0B = indexve(latB)+ips0B
        else
-         mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0
+         mat0  = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
+         mat0A = indexve(latA)+(isza0A-1)*nbpsve(latA)+ips0A
+         mat0B = indexve(latB)+(isza0B-1)*nbpsve(latB)+ips0B
        endif
-       
+
+c        write(*,*) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
+   
 c interpolation of ksi and computation of psi,deltapsi
 c ----------------------------------------------------
-       do band=1,nnuve
-        do k=0,klev+1
-         do i=0,klev+1
-          if (isza0.eq.-99) then
-           ksi = ksive(i,k,band,mat0)*(1-factp)
-     .          +ksive(i,k,band,mat0+1)*factp
-          else
-           ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
-     .          +ksive(i,k,band,mat0+1)*factp  *(1-factz)
-     .          +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
-     .          +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
-          endif
-     
-          psi(i,k) = psi(i,k) + 
-     .               RPI*ksi*(bplck(i,band)-bplck(k,band))
-          deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
-         enddo
-        enddo
-       enddo
+
+       if (isza0.eq.-99) then
+         if (latA.gt.0) then		! Not being in the two extremal bins
+
+           do band=1,nnuve
+             do k=0,klev+1
+               do i=k+1,klev+1
+	         k1 = ksive(i,k,band,mat0A)*(1-factlat)
+     .              + ksive(i,k,band,mat0B)*factlat
+	         k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
+     .              + ksive(i,k,band,mat0B+1)*factlat
+                 ksi = k1*(1-factp) + k2*factp
+                 psi(i,k) = psi(i,k) + 
+     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
+c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
+c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
+c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
+
+	         kasua=1
+               enddo
+             enddo
+           enddo
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 psi(k,i) = -psi(i,k) 
+               enddo
+             enddo
+
+         else 		! latA=0 --> extremal bins
+
+           do band=1,nnuve
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 ksi = ksive(i,k,band,mat0)*(1-factp)
+     .               + ksive(i,k,band,mat0+1)*factp
+                 psi(i,k) = psi(i,k) + 
+     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
+c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
+c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
+c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
+
+	         kasua=2
+               enddo
+             enddo
+           enddo
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 psi(k,i) = -psi(i,k) 
+               enddo
+             enddo
+
+         endif 		! latA.gt.0
+
+       else		! isza0=!-99
+
+         if (latA.gt.0) then		! Not being in the two extremal bins
+
+           do band=1,nnuve
+             do k=0,klev+1
+               do i=k+1,klev+1
+	         k1 = ksive(i,k,band,mat0A)*(1-factlat)
+     .              + ksive(i,k,band,mat0B)*factlat
+	         k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
+     .              + ksive(i,k,band,mat0B+1)*factlat
+	         k3 = ksive(i,k,band,mat0A+nbpsve(latA))*(1-factlat)
+     .              + ksive(i,k,band,mat0B+nbpsve(latB))*factlat
+	         k4 = ksive(i,k,band,mat0A+nbpsve(latA)+1)*(1-factlat)
+     .              + ksive(i,k,band,mat0B+nbpsve(latB)+1)*factlat
+	         ksi = ( k1*(1-factp) + k2*factp )*(1-factz)
+     .              + ( k3*(1-factp) + k4*factp )*factz
+                 psi(i,k) = psi(i,k) + 
+     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
+c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
+c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
+c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
+
+	        kasua=3
+               enddo
+             enddo
+           enddo
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 psi(k,i) = -psi(i,k) 
+               enddo
+             enddo
+
+         else 		! latA=0 --> extremal bins
+
+           do band=1,nnuve
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
+     .               + ksive(i,k,band,mat0+1)*factp  *(1-factz)
+     .               + ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
+     .               + ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
+                 psi(i,k) = psi(i,k) + 
+     .                      RPI*ksi*(bplck(i,band)-bplck(k,band))
+c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
+c                 deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
+c                 deltapsi(k,i) = deltapsi(k,i) + RPI*ksi*zdblay(k,band)
+
+	         kasua=4
+               enddo
+             enddo
+           enddo
+             do k=0,klev+1
+               do i=k+1,klev+1
+                 psi(k,i) = -psi(i,k) 
+               enddo
+             enddo
+
+         endif 		! latA.gt.0
+       endif		! isza0.eq.-99
+
+c	write(*,*) 'Kasua:', kasua
 
 c======================================================================
@@ -358,3 +628,2 @@
       RETURN
       END
-
