Index: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/mad_muphy.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/mad_muphy.F90	(revision 1666)
+++ trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/mad_muphy.F90	(revision 1667)
@@ -122,8 +122,8 @@
   IF (SEDIM .and. MERGE) THEN
      redge = (r1*r2)**0.5D0 ! Edge radius of mode 1
-     write(*,*)'The intersection size or virtual bonduary', redge
+!     write(*,*)'The intersection size or virtual bonduary', redge
 
      IF (r1 .EQ. redge) THEN
-        write(*,*)'Mode 1 merge to mode 2', redge, r1
+!        write(*,*)'Mode 1 merge to mode 2', redge, r1
         CALL MERGING(r1,N1,sigma1,sigma2,dM0_merge_m1,dM3_merge_m1, &
              &   dM0_merge_m2,dM3_merge_m2)
@@ -143,5 +143,5 @@
 
      IF (r2 .EQ. redge) THEN
-        write(*,*)'Mode 2 merge to mode 1',redge,r2
+!        write(*,*)'Mode 2 merge to mode 1',redge,r2
         CALL MERGING(r2,N2,sigma2,sigma1,dM0_merge_m2,dM3_merge_m2, &
              &   dM0_merge_m1,dM3_merge_m1)
@@ -325,6 +325,6 @@
   M3_m2(3) = M3_m2(3) + dM3_m2_ccn
   !------------------------------------!
-  write(*,*)'l.653: after up', M3_m1(1),M3_m1(2),M3_m1(3),dM3_m1_SA,dM3_m1_WV,dM3_m1_ccn
-  write(*,*)'l.654: after up', M0_m1(1),M0_m1(2),dM0_m1_drop
+!  write(*,*)'l.653: after up', M3_m1(1),M3_m1(2),M3_m1(3),dM3_m1_SA,dM3_m1_WV,dM3_m1_ccn
+!  write(*,*)'l.654: after up', M0_m1(1),M0_m1(2),dM0_m1_drop
 !-----------------------------------------------------------------------
   ! To controle/check the mass conservation !
@@ -375,5 +375,5 @@
         N1 = 0.D0
      ENDIF
-     write(*,*)'rmean not ZERO',r1,N1,M0_m1(1),M3_m1(1),WSA
+!     write(*,*)'rmean not ZERO',r1,N1,M0_m1(1),M3_m1(1),WSA
   ENDIF
   
@@ -408,5 +408,5 @@
   frac_ccn_m2 = M0_m2(2)/(M0_m2(1)+M0_m2(2))
 
-  write(*,*)'rmean',r1,'NTOT',N1, 'SH2SO4', SH2SO4
+!  write(*,*)'rmean',r1,'NTOT',N1, 'SH2SO4', SH2SO4
 
 
@@ -416,8 +416,8 @@
   IF (MERGE) THEN
      redge = (r1*r2)**0.5D0 ! Edge radius of mode 1
-     write(*,*)'The intersection size or virtual bonduary', redge
+!     write(*,*)'The intersection size or virtual bonduary', redge
 
      IF (r1 .EQ. redge) THEN
-        write(*,*)'Mode 1 merge to mode 2', redge, r1
+!        write(*,*)'Mode 1 merge to mode 2', redge, r1
         CALL MERGING(r1,N1,sigma1,sigma2,dM0_merge_m1,dM3_merge_m1, &
              &   dM0_merge_m2,dM3_merge_m2)
@@ -437,5 +437,5 @@
 
      IF (r2 .EQ. redge) THEN
-        write(*,*)'Mode 2 merge to mode 1',redge,r2
+!        write(*,*)'Mode 2 merge to mode 1',redge,r2
         CALL MERGING(r2,N2,sigma2,sigma1,dM0_merge_m2,dM3_merge_m2, &
              &   dM0_merge_m1,dM3_merge_m1)
@@ -476,12 +476,12 @@
 !-----------------------------------------------------------------------
   if (dt .eq. 1) then
-     write(*,*) 'Ntot, rc mode 1: ', N1,r1
+!     write(*,*) 'Ntot, rc mode 1: ', N1,r1
      ! Log-normal distribution function
      CALL build_radius_grid(nbin,rmin,rmax,vratio)
      CALL logdist(sigma1,N1,r1,rad_cld,n_g1)
-     DO i = 1, nbin
-        open(888,file='ND_t3600_nuc-het-flux_dt1_1step_sssat')
-        write(888,'(I10,2X,3(ES15.7,1X))') i, rad_cld(i), n_g1(i)
-     ENDDO
+!     DO i = 1, nbin
+!        open(888,file='ND_t3600_nuc-het-flux_dt1_1step_sssat')
+!        write(888,'(I10,2X,3(ES15.7,1X))') i, rad_cld(i), n_g1(i)
+!     ENDDO
   endif
 
@@ -507,5 +507,5 @@
 
   MFPAIR = FPLAIR(TAIR,PAIR) ! Mean free path of air molecules
-  write(*,*) 'MFPAIR',MFPAIR, TAIR, PAIR
+!  write(*,*) 'MFPAIR',MFPAIR, TAIR, PAIR
   Kn = MFPAIR/ri1             ! Knudsen number 
 
@@ -531,5 +531,5 @@
   real, intent(in) :: TAIR, PAIR, MRSA_loc
   real :: ROSAS, PSASAS_ZELE,  STSAS
-  write(*,*)'checkpoint change_wsa', ST, RHOSA, RHOsasat, PSASAS_ZELE(TAIR,WSA), TAIR, WSA
+!  write(*,*)'checkpoint change_wsa', ST, RHOSA, RHOsasat, PSASAS_ZELE(TAIR,WSA), TAIR, WSA
   ST = STSAS(TAIR,WSA)      ! Surface tension of sulfuric acid solution/vapor (N/m)
 
@@ -539,5 +539,5 @@
   SH2SO4 = (PAIR*MRSA_loc)/RHOsasat    ! Saturation Ratio
 
-  write(*,*)'checkpoint change_wsa', PSASAS_ZELE(TAIR,WSA)
+!  write(*,*)'checkpoint change_wsa', PSASAS_ZELE(TAIR,WSA)
 
 END SUBROUTINE change_wsa
Index: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F	(revision 1667)
+++ trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F	(revision 1667)
@@ -0,0 +1,344 @@
+      SUBROUTINE new_cloud_sedim(n_lon,n_lev,ptimestep,
+     &                pmidlay,pbndlay, 
+     &                pt,
+     &                pq, pdqsed,pdqs_sed,nq,F_sed)
+
+      USE ioipsl
+      USE dimphy
+      USE chemparam_mod
+      IMPLICIT NONE
+
+c=======================================================================
+c      
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+#include "YOMCST.h"     
+c#include "dimphys.h"
+c#include "comcstfi.h"
+c#include "tracer.h"
+c#include "callkeys.h"
+	
+c
+c   arguments:
+c   ----------
+
+      INTEGER n_lon                 ! number of horizontal grid points
+      INTEGER n_lev                 ! number of atmospheric layers
+      REAL ptimestep                ! physics time step (s)
+      REAL pmidlay(n_lon,n_lev)     ! pressure at middle layers (Pa)
+      REAL pt(n_lon,n_lev)          ! temperature at mid-layer (l)
+      REAL pbndlay(n_lon,n_lev+1)   ! pressure at layer boundaries
+
+c    Traceurs :
+      integer nq                    ! number of tracers
+      real pq(n_lon,n_lev,nq)       ! tracers (kg/kg)
+c      real pdqfi(n_lon,n_lev,nq)   ! tendency before sedimentation (kg/kg.s-1)
+      real pdqsed(n_lon,n_lev,2)    ! tendency due to sedimentation (kg/kg)
+      real pdqs_sed(n_lon)          ! surface density (Flux if /ptimestep) at surface due to sedimentation (kg.m-2)
+      
+c   local:
+c   ------
+      integer imode
+      integer ig
+      integer iq
+      integer l
+
+      real zlev(n_lon,n_lev+1)      ! altitude at layer boundaries
+      real zlay(n_lon,n_lev)        ! altitude at the midlle layer
+      real zqi_wv(n_lon,n_lev)      ! to locally store H2O tracer
+      real zqi_sa(n_lon,n_lev)      ! to locally store H2SO4 tracer
+      real m_lay (n_lon,n_lev)      ! Layer Pressure over gravity (Dp/g == kg.m-2)
+      real wq(n_lon,n_lev+1)        ! displaced tracer mass (kg.m-2)
+
+c    Physical constant
+c    ~~~~~~~~~~~~~~~~~
+c     Gas molecular viscosity (N.s.m-2)
+c      real,parameter :: visc=1.e-5       ! CO2
+	REAL :: VISCOSITY_CO2
+c     Effective gas molecular radius (m)
+      real,parameter :: molrad=2.2e-10   ! CO2
+      
+c     Cloud density (kg.m-3)
+c     ~~~~~~~~~~~~~~~~~~~~~~
+c      real, DIMENSION(n_lon,n_lev) ::  rho_droplet
+
+      REAL, DIMENSION(n_lon,n_lev+1) ::
+     + wgt_SA                         ! Fraction of H2SO4 in droplet local
+
+c     Stokes speed and sedimentation flux variable
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      REAL :: A1,A2,A3,A4,            ! coeff du DL du Flux de sedimentation
+     + D_stokes,                      ! coeff de la vitesse de Stokes
+     + Rp_DL,                         ! "Point" du DL
+     + l_mean,                        ! libre parcours moyen (m)
+     + a,b_exp,c                      ! coeff du calcul du Flux de sedimentation
+      REAL, DIMENSION(n_lon,n_lev+1) ::
+     + F_sed                          ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
+      
+      
+      REAL :: R_mode0                 ! Rayon mode 0 (m), rayon le plus frequent
+
+
+
+!      PRINT*,'RHO_DROPLET new_cloud_sedim.F'
+!      PRINT*,'rho_droplet',rho_droplet(16,21)
+!      PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
+
+c-----------------------------------------------------------------------
+c    1. Initialization
+c    -----------------
+    
+c     Updating the droplet mass mixing ratio with the partition H2O/H2SO4
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      do l=1,n_lev
+         do ig=1,n_lon
+         zqi_wv(ig,l) = pq(ig,l,i_h2oliq)
+         zqi_sa(ig,l) = pq(ig,l,i_h2so4liq)
+         wgt_SA(ig,l) = WH2SO4(ig,l)
+         enddo
+      enddo 
+
+c     Init F_sed
+      F_sed(:,:) = 0.0E+0
+
+c     Au niveau top+1 , tout égal a 0      
+      wgt_SA(:,n_lev+1) = 0.0E+0   
+
+c    Computing the different layer properties
+c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c    m_lay (kg.m-2)
+c    Ici g=8.87, conflit pour g entre #include "YOMCST.h" 
+c	et #include "comcstfi.h"
+
+      do  l=1,n_lev
+         do ig=1, n_lon
+         m_lay(ig,l)=(pbndlay(ig,l) - pbndlay(ig,l+1)) /8.87E+0 
+            IF (m_lay(ig,l).LE.0.0) THEN
+            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
+            PRINT*,'!!!!          m_lay <= 0        !!!!'
+            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
+            ENDIF
+         end do
+      end do
+	
+c         Computing sedimentation for droplet "tracer"
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c	    pbndlay(:,51)=0 (en parallèle c'est sûr), ne pas l'utiliser pour Fse
+	
+	DO imode=1, nbr_mode
+         DO l = cloudmin, cloudmax
+            DO ig=1,n_lon
+
+c     RD=1000.*RNAVO*RKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1     
+           D_stokes=((rho_droplet(ig,l)-pmidlay(ig,l)/(RD*pt(ig,l))))
+     &	*(2./9.)*(RG/VISCOSITY_CO2(pt(ig,l)))
+      
+           l_mean=(pt(ig,l)/pmidlay(ig,l))*
+     &	(0.707*R/(4.*RPI* molrad*molrad * RNAVO))
+      
+           R_mode0=R_MEDIAN(ig,l,imode)*
+     &     EXP(-LOG(STDDEV(ig,l,imode))**2.)
+              IF ((l_mean/(R_mode0)).GT.10.) THEN
+              Rp_DL=R_MEDIAN(ig,l,imode)*
+     &        EXP(3.*LOG(STDDEV(ig,l,imode))**2.)
+              ELSE
+              Rp_DL=R_MEDIAN(ig,l,imode)*
+     &        EXP(4.*LOG(STDDEV(ig,l,imode))**2.)
+              ENDIF
+		
+           a=1.246*l_mean
+	
+           c=0.87/l_mean
+	
+           b_exp=0.42*l_mean*EXP(-c*Rp_DL)
+	
+           A1=a+b_exp*(1.+c*Rp_DL
+     &	+0.5*(Rp_DL*c)**2
+     &	+1./6.*(Rp_DL*c)**3)
+     
+           A2=1.-b_exp*(c
+     &	+Rp_DL*c**2
+     &	+0.5*(Rp_DL**2)*(c**3))
+	
+           A3=0.5*b_exp*(c**2+Rp_DL*c**3)
+	
+           A4=-b_exp*1./6.*c**3
+
+c     Addition des Flux de tous les modes presents	
+       F_sed(ig,l)=F_sed(ig,l)+(rho_droplet(ig,l)*4./3.*RPI*
+     &  NBRTOT(ig,l,imode)*1.0E6*D_stokes*(
+     &  A1*R_MEDIAN(ig,l,imode)**4
+     &  *EXP(8.0*LOG(STDDEV(ig,l,imode))**2.)
+     &  +A2*R_MEDIAN(ig,l,imode)**5
+     &  *EXP(12.5*LOG(STDDEV(ig,l,imode))**2.)
+     &  +A3*R_MEDIAN(ig,l,imode)**6
+     &  *EXP(18.0*LOG(STDDEV(ig,l,imode))**2.)
+     &  +A4*R_MEDIAN(ig,l,imode)**7
+     &  *EXP(24.5*LOG(STDDEV(ig,l,imode))**2.)))
+      
+c      PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
+      
+      	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
+      	PRINT*,'==============================================='
+      	PRINT*,'WARNING On a epuise la couche', ig, l
+      	PRINT*,'On epuise pas une couche avec une espèce
+     & 	 minoritaire, c est pas bien maaaaaal'
+            PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
+      	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
+      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
+      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
+     & 	rho_droplet(ig,l)
+     		PRINT*,'Ntot',NBRTOT(ig,l,:)
+     		PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
+     		PRINT*,'K_MASS',K_MASS(ig,l,:)
+     		PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
+     		
+c     		ELSE
+c     		
+c     		PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+c      	PRINT*,'WARNING On a PAS epuise la couche', ig, l
+c      	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
+c      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
+c      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
+c     & 	rho_droplet(ig,l)(ig,l)
+c     		PRINT*,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)*1.0e6
+c     		PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l) 
+            STOP    		
+     	      ENDIF
+	
+           IF (F_sed(ig,l).LT.0.0d0) THEN
+	      PRINT*,"F_sed est négatif !!!"
+	      PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
+      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
+      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
+      	PRINT*,'Temp',pt(ig,l),'Rho',
+     & 	rho_droplet(ig,l)
+     		PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
+     &	NBRTOT(ig,l,imode)*1.0e6
+     		PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
+     &		R_MEDIAN(ig,l,imode)
+     		PRINT*,'A1',A1,'A2',A2
+     		PRINT*,'A3',A1,'A4',A2
+     		PRINT*,'D_stokes',D_stokes
+     		STOP
+           ENDIF
+           
+	      ENDDO
+	      
+c           ELSE           
+c           F_sed(:,l)=0.0d0           
+c           ENDIF
+           
+         ENDDO
+      ENDDO
+
+c     Passage du Flux au Flux pour un pas de temps (== kg.m-2)     
+      F_sed(:,:)=F_sed(:,:)*ptimestep
+
+
+c	VENUS: le flux à la surface est fixé à 0
+c     les conditions P/T en surface ne permettent pas la condensation
+      DO ig=1,n_lon
+      pdqs_sed(ig) = 0.0d0
+      ENDDO
+      
+c         Compute the final tendency:
+c         ---------------------------
+
+c     Partie H2SO4l
+c     ~~~~~~~~~~~~
+
+      DO l = 1, n_lev
+         DO ig=1,n_lon
+            zqi_sa(ig,l) = zqi_sa(ig,l) + (
+     &                         F_sed(ig,l+1)*wgt_SA(ig,l+1)
+     &                       - F_sed(ig,l)*wgt_SA(ig,l))
+     &                       / m_lay(ig,l)
+c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
+             IF (zqi_sa(ig,l).LT.0.0D0) THEN
+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
+c              de tps de la phys, microphys et chimie
+c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
+c              mais en fait on l epuise pour un pdt<pdtphys
+               zqi_sa(ig,l) = 0.0D0
+             ENDIF
+            pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)                       
+         ENDDO
+      ENDDO
+
+c     Partie H2Ol
+c     ~~~~~~~~~~~
+                     
+      DO l = 1, n_lev
+         DO ig=1,n_lon
+            zqi_wv(ig,l) = zqi_wv(ig,l) + (
+     &                         F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
+     &                       - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
+     &                       / m_lay(ig,l)
+c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
+             IF (zqi_wv(ig,l).LT.0.0D0) THEN
+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
+c              de tps de la phys, microphys et chimie
+c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
+c              mais en fait on l epuise pour un pdt<pdtphys
+               zqi_wv(ig,l) = 0.0D0
+             ENDIF
+            pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)                   
+         ENDDO
+      ENDDO
+
+c		Save output file in 1D model
+c		============================ 
+	        
+c      IF (n_lon .EQ. 1) THEN
+c      PRINT*,'Save output sedim' 	
+c      DO l = 1, n_lev
+c      	DO ig=1,n_lon
+c      	WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
+c     & 	(WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
+c     & 	(1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
+c     & 	pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
+c      ENDDO
+c      	ENDDO
+c      ENDIF   
+
+      RETURN
+      END
+
+*******************************************************************************
+      REAL FUNCTION VISCOSITY_CO2(temp)
+c	Aurélien Stolzenbach 2015
+c	Calcul de la viscosité dynamique du CO2 80°K -> 300°K
+c	Viscosité dynamique en Pa.s
+c	Source: Johnston & Grilly (1942)
+
+c	température en °K
+	REAL, INTENT(IN) :: temp
+	
+	REAL :: denom, numer
+	
+c	Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
+	
+	numer = 200.**(2.27/4.27)-0.435
+	denom = temp**(2.27/4.27)-0.435
+	
+	VISCOSITY_CO2 = (numer/denom)*1015.*(temp/200.)**(3./2.)
+
+c	convertion de Poises*1e7 -> Pa.s	
+	VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8	
+
+      END FUNCTION VISCOSITY_CO2
+*******************************************************************************
+
+
Index: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_venus.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_venus.F	(revision 1667)
+++ trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_venus.F	(revision 1667)
@@ -0,0 +1,406 @@
+!******************************************************************************
+!*     venus_SAS_composition SUBROUTINE
+!*     modified from
+!*     PROGRAM PSC_MODEL_E
+!*     by A. Määttänen
+!*     subroutine for LMDZ+photochemistry VENUS
+!*     by A. Stolzenbach
+!*
+!*     Input/Output files:
+!*     -------------------
+!*
+!----------------------------------------------------------------------------
+      SUBROUTINE new_cloud_venus(
+     + nblev, nblon,
+     + TT,PP,
+     + mrt_wv,mrt_sa,
+     + mr_wv,mr_sa)
+
+      USE chemparam_mod
+      IMPLICIT NONE
+
+#include "YOMCST.h"
+
+      INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
+      INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales
+
+!----------------------------------------------------------------------------
+!     Ambient air state variables:
+      REAL, INTENT(IN), DIMENSION(nblon,nblev) :: mrt_wv,mrt_sa,
+     +                                            TT,PP
+      REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: mr_wv,mr_sa
+!----------------------------------------------------------------------------
+      INTEGER :: ilon, ilev, imode
+!----------------------------------------------------------------------------
+!     Thermodynamic functions:
+      REAL :: ROSAS
+!----------------------------------------------------------------------------
+!     Auxilary variables:
+      REAL :: NH2SO4,NH2O
+      REAL :: H2SO4_liq,H2O_liq
+      REAL :: CONCM
+      REAL :: MCONDTOT
+      REAL :: RMODE
+      REAL :: WSAFLAG
+      REAL :: K_SAV
+!----------------------------------------------------------------------------
+!     Ridder's Method variables:
+      REAL :: WVMIN, WVMAX, WVACC
+
+      INTEGER :: NBROOT
+
+      INTEGER :: MAXITE
+      PARAMETER(MAXITE=20)
+
+      INTEGER :: NBRAC
+      PARAMETER(NBRAC=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:
+
+!AM Venus
+! These aerosols will then be given an equilibrium composition for the given size distribution
+
+  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
+  ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
+  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
+  !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
+
+!===========================================
+!     Debut boucle sur niveau et lat,lon
+!===========================================
+!     Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0
+      NBRTOT(:,:,:)=0.0E+0
+      WH2SO4(:,:)=0.0E+0
+      rho_droplet(:,:)=0.0E+0
+
+      DO ilev=cloudmin, cloudmax
+      DO ilon=1, nblon
+
+!       Boucle sur les modes
+        RMODE=0.0E+0
+        K_SAV = 0.0
+
+        DO imode=1, nbr_mode
+          IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
+!       RMODE est le rayon modal de la distribution en volume du mode le plus
+!       representatif pour la Mtot
+            RMODE=R_MEDIAN(ilon,ilev,imode)*
+     &      EXP(2.*(DLOG(STDDEV(ilon,ilev,imode))**2.))
+            K_SAV=K_MASS(ilon,ilev,imode)
+          ENDIF
+        ENDDO ! FIN boucle imode
+
+!       Initialisation des bornes pour WV
+        WVMIN=1.E-90
+        WVMAX=mrt_wv(ilon,ilev)
+
+!       Accuracy de WVeq
+        WVACC=WVMAX*1.0E-3
+
+!       BRACWV borne la fonction f(WV) - WV = 0
+!       de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
+!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
+!       Elle fait appel à la fct/ssrtine ITERWV()
+
+        CALL BRACWV(TT(ilon,ilev),PP(ilon,ilev),WVMIN,WVMAX,NBRAC,
+     &  RMODE,mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),FLAG,WSAFLAG,NBROOT)
+
+        SELECT CASE(FLAG)
+
+        CASE(1)
+!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
+!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
+!       Elle fait appel  la fct/ssrtine ITERWV()
+
+          WH2SO4(ilon,ilev)=IRFRMWV(TT(ilon,ilev),PP(ilon,ilev),
+     &    WVMIN,WVMAX,WVACC,MAXITE,RMODE,
+     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
+
+          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
+
+!          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
+!            PRINT*,'PROBLEM RHO_DROPLET'
+!            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
+!            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
+!            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
+!            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
+!            STOP
+!          ENDIF
+
+          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
+
+	      NH2SO4=mrt_sa(ilon,ilev)*CONCM
+	      NH2O=mrt_wv(ilon,ilev)*CONCM
+
+          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
+     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
+     &       H2SO4_liq,H2O_liq,MCONDTOT)
+
+!       Boucle sur les modes
+! *********************************************
+!	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)/(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
+!	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
+          H2O_liq=H2O_liq/CONCM
+          H2SO4_liq=H2SO4_liq/CONCM
+
+          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
+          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
+
+!          Problemes quand on a condense tout, on peut obtenir des -1e-24
+!		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_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
+
+
+
+        CASE(2)
+!       Cas NROOT=0 mais proche de 0
+
+          WH2SO4(ilon,ilev)=WSAFLAG
+
+          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(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
+!     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 ROSAS (cf test externe en Python), sinon, la valeur est trop
+!     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 bonne valeur (tests externes Python confirment)
+
+          IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
+     &      (WH2SO4(ilon,ilev).GT.0.1))THEN
+            PRINT*,'PROBLEM RHO_DROPLET'
+            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
+            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
+            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
+            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
+            STOP
+          ENDIF
+
+
+          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
+
+	    NH2SO4=mrt_sa(ilon,ilev)*CONCM
+	    NH2O=mrt_wv(ilon,ilev)*CONCM
+
+          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
+     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
+     &       H2SO4_liq,H2O_liq,MCONDTOT)
+
+!       Boucle sur les modes
+! *********************************************
+!	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)/(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
+!	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
+          H2O_liq=H2O_liq/CONCM
+          H2SO4_liq=H2SO4_liq/CONCM
+
+          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
+          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
+
+!          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_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
+
+        CASE(3)
+!         Cas 0 NROOT
+            mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
+            mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
+            rho_droplet(ilon,ilev)=0.0E+0
+            WH2SO4(ilon,ilev)=0.0E+0
+            DO imode=1, nbr_mode
+              NBRTOT(ilon,ilev,imode)=0.0E+0
+            ENDDO
+
+        END SELECT
+      ENDDO   !FIN boucle ilon
+      ENDDO   !FIN boucle ilev
+
+      END SUBROUTINE new_cloud_venus
+
+
+!******************************************************************
+      SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
+     + T,H2SO4COND,H2OCOND,RMTOT)
+
+!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
+!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
+!                                       SIZE DISTRIBTUION
+!     ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
+!    ---------------------------------------------------------------
+!     INPUT:
+!     H2SO4: #/m3 of total H2SO4
+!	H2O  : #/m3 of total H2O
+!     WSA: aerosol H2SO4 weight fraction (fraction)
+!     DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume)
+!       for total mass, almost same result with ro=1.67 gr/cm3
+!     RSTDEV: standard deviation of aerosol distribution (no unit)
+!     RADIUS: MEDIAN radius (m)
+!     T: temperature (K)
+!
+!     OUTPUT:
+!     RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
+!            mais rho_droplet et M_tot_distrib doivent etre de meme dimension
+!	H2OCOND
+!	H2SO4COND
+
+
+
+      IMPLICIT NONE
+
+      REAL, INTENT(IN) :: H2SO4, H2O, WSA
+      REAL, INTENT(IN) :: DENSO4, T
+      REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
+!     working variables
+      REAL :: RMH2S4
+      REAL :: DND2,pstand,lpar,acidps
+      REAL :: x1, satpacid
+      REAL , DIMENSION(2):: act
+!
+!     masse of an H2SO4 molecule (kg)
+      RMH2S4=98.078/(6.02214129E+26)
+
+      pstand=1.01325E+5 !Pa  1 atm pressure
+
+        x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
+
+        call zeleznik(x1,t,act)
+
+!pure acid satur vapor pressure
+        lpar= -11.695+DLOG(pstand) ! Zeleznik
+        acidps=1/360.15-1.0/t+0.38/545.
+     + *(1.0+DLOG(360.15/t)-360.15/t)
+        acidps = 10156.0*acidps +lpar
+        acidps = DEXP(acidps)    !Pa
+
+!acid sat.vap.PP over mixture (flat surface):
+        satpacid=act(2)*acidps ! Pa
+
+!       Conversion from Pa to N.D #/m3
+        DND2=satpacid/(1.3806488E-23*T)
+!       Conversion from N.D #/m3 TO #/cm3
+!        DND2=DND2*1.d-6
+
+!	H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
+	IF (H2SO4.GE.DND2) THEN
+	H2SO4COND=H2SO4-DND2
+!	calcul de H2O cond correspondant a H2SO4 cond
+	H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
+
+!     RMTOT: = Mass of H2SO4 satur per 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)
+
+!	Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
+	ELSE
+	H2SO4COND=0.0E+0
+	H2OCOND=0.0E+0
+	RMTOT=0.0E+0
+	END IF
+
+!	Test si H2O en defaut H2Ocond>H2O dispo
+	IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
+
+!     Si H2O en dfaut, on as pas le bon WSA!
+!     En effet, normalement, on a exactement le WSA correspondant a
+!     WVg + WVl = WVtot
+!     Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou
+!     quand T est grand (Lower Haze), le modle reprsente mal le WSA
+!     cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA
+
+!      PRINT*,'PROBLEM H2O EN DEFAUT'
+!      PRINT*,'H2OCOND',H2OCOND,'H2O',H2O
+!      PRINT*,'WSA',WSA,'RHO',DENSO4
+!      STOP
+
+
+!	On peut alors condenser tout le H2O dispo
+	H2OCOND=H2O
+!	On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
+	H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
+
+!     RMTOT: = Mass of H2SO4 satur per 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
+
Index: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/wsa_new.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/wsa_new.F90	(revision 1666)
+++ trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/wsa_new.F90	(revision 1667)
@@ -204,5 +204,5 @@
   WSAMIN = 0.1D0
   WSAMAX = 1.0D0
-  LPPWV = LOG(PAIR*WV)
+  LPPWV = DLOG(PAIR*WV)
 
   ! Appel Bracket de KEEQ         
@@ -293,5 +293,5 @@
         ! Si NROOT>1 on place la borne sup output  la borne min du calcul en i             
         IF (NROOT.GT.1) THEN
-           XB=(1.-LOG(DBLE(I+1))/LOG(DBLE(N)))*XMAX
+           XB=(1.-DLOG(DBLE(I+1))/DLOG(DBLE(N)))*XMAX
         ENDIF
 
@@ -314,5 +314,5 @@
 
      DO J=N-1,1,-1
-        X=(1.-LOG(DBLE(N-J))/LOG(DBLE(N)))*XMAX
+        X=(1.-DLOG(DBLE(N-J))/DLOG(DBLE(N)))*XMAX
         !             write(*,*) 'BRACWV, bf 4th ITERWV (cas 2) '
         CALL ITERWV(TAIR,PAIR,X,WVLIQ,WVEQ,WVTOT,WSAOUT,SATOT,RADIUS)
@@ -365,10 +365,10 @@
   DX=(XB-XA)/N
   X=XB
-  FP=KEEQ(RADIUS,X,LPPWVINP)
+  FP=KEEQ(TAIR,RADIUS,X,LPPWVINP)
 
   DO I=N,1,-1
      X=X-DX
 
-     FC=KEEQ(RADIUS,X,LPPWVINP)
+     FC=KEEQ(TAIR,RADIUS,X,LPPWVINP)
 
      IF ((FP*FC).LE.0.) THEN
@@ -752,5 +752,5 @@
 
   ! Physical constants:
-  REAL :: MH2O, KEEQ
+  REAL :: KEEQ
 
   !     External functions needed:
@@ -764,6 +764,5 @@
   REAL :: C1
 
-  MH2O= MWV
-  C1=2.0D0*MH2O/RGAS
+  C1=2.0D0*MWV/RGAS
 
   KEEQ=LPPWV-C1*STSAS(TAIR,WX)/(TAIR*RADIUS*ROSAS(TAIR,WX))- &
@@ -818,8 +817,8 @@
 
   !pure acid satur vapor pressure
-  lpar= -11.695D0 + LOG(pstand) ! Zeleznik
-  acidps = 1/360.15D0 - 1.0/T + 0.38D0/545.D0*(1.0+LOG(360.15D0/T)-360.15D0/T)
+  lpar= -11.695D0 + DLOG(pstand) ! Zeleznik
+  acidps = 1/360.15D0 - 1.0/T + 0.38D0/545.D0*(1.0+DLOG(360.15D0/T)-360.15D0/T)
   acidps = 10156.D0*acidps + lpar
-  acidps = EXP(acidps)    !Pa
+  acidps = DEXP(acidps)    !Pa
 
   !acid sat.vap.PP over mixture (flat surface):
Index: trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h	(revision 1666)
+++ trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h	(revision 1667)
@@ -225,24 +225,24 @@
      .                "ave(X)", zsto,zout)
 c
-         CALL histdef(nid_mth,"dtlwrNLTE", "LWNLTE radiation dT",
-     .                "K/s",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
-     .                "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_mth,"dtlwrLTE","LW_LTE radiation dT",
-     .                "K/s",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
-     .                "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_mth, "dteuv", "UV radiation dT", "K/s",
-     .                nbp_lon,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_mth, "dtcond", "Therm conduction", "K/s",
-     .                nbp_lon,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_mth, "dumolvis", "molec viscosity (u)"
-     .                ,"m/s2",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_mth, "dvmolvis", "molec viscosity (v)", 
-     .                "m/s2",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
-     .                "ave(X)", zsto,zout)
+c        CALL histdef(nid_mth,"dtlwrNLTE", "LWNLTE radiation dT",
+c    .                "K/s",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
+c    .                "ave(X)", zsto,zout)
+c
+c        CALL histdef(nid_mth,"dtlwrLTE","LW_LTE radiation dT",
+c    .                "K/s",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
+c    .                "ave(X)", zsto,zout)
+c
+c        CALL histdef(nid_mth, "dteuv", "UV radiation dT", "K/s",
+c    .                nbp_lon,jj_nb,nhori, klev,1,klev,nvert, 32,
+c    .                "ave(X)", zsto,zout)
+c        CALL histdef(nid_mth, "dtcond", "Therm conduction", "K/s",
+c    .                nbp_lon,jj_nb,nhori, klev,1,klev,nvert, 32,
+c    .                "ave(X)", zsto,zout)
+c        CALL histdef(nid_mth, "dumolvis", "molec viscosity (u)"
+c    .                ,"m/s2",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
+c    .                "ave(X)", zsto,zout)
+c        CALL histdef(nid_mth, "dvmolvis", "molec viscosity (v)", 
+c    .                "m/s2",nbp_lon,jj_nb,nhori,klev,1,klev,nvert,32,
+c    .                "ave(X)", zsto,zout)
 
 c
Index: trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F	(revision 1666)
+++ 	(revision )
@@ -1,344 +1,0 @@
-      SUBROUTINE new_cloud_sedim(n_lon,n_lev,ptimestep,
-     &                pmidlay,pbndlay, 
-     &                pt,
-     &                pq, pdqsed,pdqs_sed,nq,F_sed)
-
-      USE ioipsl
-      USE dimphy
-      USE chemparam_mod
-      IMPLICIT NONE
-
-c=======================================================================
-c      
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-#include "YOMCST.h"     
-c#include "dimphys.h"
-c#include "comcstfi.h"
-c#include "tracer.h"
-c#include "callkeys.h"
-	
-c
-c   arguments:
-c   ----------
-
-      INTEGER n_lon                 ! number of horizontal grid points
-      INTEGER n_lev                 ! number of atmospheric layers
-      REAL ptimestep                ! physics time step (s)
-      REAL pmidlay(n_lon,n_lev)     ! pressure at middle layers (Pa)
-      REAL pt(n_lon,n_lev)          ! temperature at mid-layer (l)
-      REAL pbndlay(n_lon,n_lev+1)   ! pressure at layer boundaries
-
-c    Traceurs :
-      integer nq                    ! number of tracers
-      real pq(n_lon,n_lev,nq)       ! tracers (kg/kg)
-c      real pdqfi(n_lon,n_lev,nq)   ! tendency before sedimentation (kg/kg.s-1)
-      real pdqsed(n_lon,n_lev,2)    ! tendency due to sedimentation (kg/kg)
-      real pdqs_sed(n_lon)          ! surface density (Flux if /ptimestep) at surface due to sedimentation (kg.m-2)
-      
-c   local:
-c   ------
-      integer imode
-      integer ig
-      integer iq
-      integer l
-
-      real zlev(n_lon,n_lev+1)      ! altitude at layer boundaries
-      real zlay(n_lon,n_lev)        ! altitude at the midlle layer
-      real zqi_wv(n_lon,n_lev)      ! to locally store H2O tracer
-      real zqi_sa(n_lon,n_lev)      ! to locally store H2SO4 tracer
-      real m_lay (n_lon,n_lev)      ! Layer Pressure over gravity (Dp/g == kg.m-2)
-      real wq(n_lon,n_lev+1)        ! displaced tracer mass (kg.m-2)
-
-c    Physical constant
-c    ~~~~~~~~~~~~~~~~~
-c     Gas molecular viscosity (N.s.m-2)
-c      real,parameter :: visc=1.e-5       ! CO2
-	REAL :: VISCOSITY_CO2
-c     Effective gas molecular radius (m)
-      real,parameter :: molrad=2.2e-10   ! CO2
-      
-c     Cloud density (kg.m-3)
-c     ~~~~~~~~~~~~~~~~~~~~~~
-c      real, DIMENSION(n_lon,n_lev) ::  rho_droplet
-
-      REAL, DIMENSION(n_lon,n_lev+1) ::
-     + wgt_SA                         ! Fraction of H2SO4 in droplet local
-
-c     Stokes speed and sedimentation flux variable
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      REAL :: A1,A2,A3,A4,            ! coeff du DL du Flux de sedimentation
-     + D_stokes,                      ! coeff de la vitesse de Stokes
-     + Rp_DL,                         ! "Point" du DL
-     + l_mean,                        ! libre parcours moyen (m)
-     + a,b_exp,c                      ! coeff du calcul du Flux de sedimentation
-      REAL, DIMENSION(n_lon,n_lev+1) ::
-     + F_sed                          ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
-      
-      
-      REAL :: R_mode0                 ! Rayon mode 0 (m), rayon le plus frequent
-
-
-
-!      PRINT*,'RHO_DROPLET new_cloud_sedim.F'
-!      PRINT*,'rho_droplet',rho_droplet(16,21)
-!      PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
-
-c-----------------------------------------------------------------------
-c    1. Initialization
-c    -----------------
-    
-c     Updating the droplet mass mixing ratio with the partition H2O/H2SO4
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      do l=1,n_lev
-         do ig=1,n_lon
-         zqi_wv(ig,l) = pq(ig,l,i_h2oliq)
-         zqi_sa(ig,l) = pq(ig,l,i_h2so4liq)
-         wgt_SA(ig,l) = WH2SO4(ig,l)
-         enddo
-      enddo 
-
-c     Init F_sed
-      F_sed(:,:) = 0.0E+0
-
-c     Au niveau top+1 , tout égal a 0      
-      wgt_SA(:,n_lev+1) = 0.0E+0   
-
-c    Computing the different layer properties
-c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c    m_lay (kg.m-2)
-c    Ici g=8.87, conflit pour g entre #include "YOMCST.h" 
-c	et #include "comcstfi.h"
-
-      do  l=1,n_lev
-         do ig=1, n_lon
-         m_lay(ig,l)=(pbndlay(ig,l) - pbndlay(ig,l+1)) /8.87E+0 
-            IF (m_lay(ig,l).LE.0.0) THEN
-            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
-            PRINT*,'!!!!          m_lay <= 0        !!!!'
-            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
-            ENDIF
-         end do
-      end do
-	
-c         Computing sedimentation for droplet "tracer"
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c	    pbndlay(:,51)=0 (en parallèle c'est sûr), ne pas l'utiliser pour Fse
-	
-	DO imode=1, nbr_mode
-         DO l = cloudmin, cloudmax
-            DO ig=1,n_lon
-
-c     RD=1000.*RNAVO*RKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1     
-           D_stokes=((rho_droplet(ig,l)-pmidlay(ig,l)/(RD*pt(ig,l))))
-     &	*(2./9.)*(RG/VISCOSITY_CO2(pt(ig,l)))
-      
-           l_mean=(pt(ig,l)/pmidlay(ig,l))*
-     &	(0.707*R/(4.*RPI* molrad*molrad * RNAVO))
-      
-           R_mode0=R_MEDIAN(ig,l,imode)*
-     &     EXP(-LOG(STDDEV(ig,l,imode))**2.)
-              IF ((l_mean/(R_mode0)).GT.10.) THEN
-              Rp_DL=R_MEDIAN(ig,l,imode)*
-     &        EXP(3.*LOG(STDDEV(ig,l,imode))**2.)
-              ELSE
-              Rp_DL=R_MEDIAN(ig,l,imode)*
-     &        EXP(4.*LOG(STDDEV(ig,l,imode))**2.)
-              ENDIF
-		
-           a=1.246*l_mean
-	
-           c=0.87/l_mean
-	
-           b_exp=0.42*l_mean*EXP(-c*Rp_DL)
-	
-           A1=a+b_exp*(1.+c*Rp_DL
-     &	+0.5*(Rp_DL*c)**2
-     &	+1./6.*(Rp_DL*c)**3)
-     
-           A2=1.-b_exp*(c
-     &	+Rp_DL*c**2
-     &	+0.5*(Rp_DL**2)*(c**3))
-	
-           A3=0.5*b_exp*(c**2+Rp_DL*c**3)
-	
-           A4=-b_exp*1./6.*c**3
-
-c     Addition des Flux de tous les modes presents	
-       F_sed(ig,l)=F_sed(ig,l)+(rho_droplet(ig,l)*4./3.*RPI*
-     &  NBRTOT(ig,l,imode)*1.0E6*D_stokes*(
-     &  A1*R_MEDIAN(ig,l,imode)**4
-     &  *EXP(8.0*LOG(STDDEV(ig,l,imode))**2.)
-     &  +A2*R_MEDIAN(ig,l,imode)**5
-     &  *EXP(12.5*LOG(STDDEV(ig,l,imode))**2.)
-     &  +A3*R_MEDIAN(ig,l,imode)**6
-     &  *EXP(18.0*LOG(STDDEV(ig,l,imode))**2.)
-     &  +A4*R_MEDIAN(ig,l,imode)**7
-     &  *EXP(24.5*LOG(STDDEV(ig,l,imode))**2.)))
-      
-c      PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
-      
-      	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
-      	PRINT*,'==============================================='
-      	PRINT*,'WARNING On a epuise la couche', ig, l
-      	PRINT*,'On epuise pas une couche avec une espèce
-     & 	 minoritaire, c est pas bien maaaaaal'
-            PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
-      	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
-      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
-      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
-     & 	rho_droplet(ig,l)
-     		PRINT*,'Ntot',NBRTOT(ig,l,:)
-     		PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
-     		PRINT*,'K_MASS',K_MASS(ig,l,:)
-     		PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
-     		
-c     		ELSE
-c     		
-c     		PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-c      	PRINT*,'WARNING On a PAS epuise la couche', ig, l
-c      	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
-c      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
-c      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
-c     & 	rho_droplet(ig,l)(ig,l)
-c     		PRINT*,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)*1.0e6
-c     		PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l) 
-            STOP    		
-     	      ENDIF
-	
-           IF (F_sed(ig,l).LT.0.0d0) THEN
-	      PRINT*,"F_sed est négatif !!!"
-	      PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
-      	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
-      	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
-      	PRINT*,'Temp',pt(ig,l),'Rho',
-     & 	rho_droplet(ig,l)
-     		PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
-     &	NBRTOT(ig,l,imode)*1.0e6
-     		PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
-     &		R_MEDIAN(ig,l,imode)
-     		PRINT*,'A1',A1,'A2',A2
-     		PRINT*,'A3',A1,'A4',A2
-     		PRINT*,'D_stokes',D_stokes
-     		STOP
-           ENDIF
-           
-	      ENDDO
-	      
-c           ELSE           
-c           F_sed(:,l)=0.0d0           
-c           ENDIF
-           
-         ENDDO
-      ENDDO
-
-c     Passage du Flux au Flux pour un pas de temps (== kg.m-2)     
-      F_sed(:,:)=F_sed(:,:)*ptimestep
-
-
-c	VENUS: le flux à la surface est fixé à 0
-c     les conditions P/T en surface ne permettent pas la condensation
-      DO ig=1,n_lon
-      pdqs_sed(ig) = 0.0d0
-      ENDDO
-      
-c         Compute the final tendency:
-c         ---------------------------
-
-c     Partie H2SO4l
-c     ~~~~~~~~~~~~
-
-      DO l = 1, n_lev
-         DO ig=1,n_lon
-            zqi_sa(ig,l) = zqi_sa(ig,l) + (
-     &                         F_sed(ig,l+1)*wgt_SA(ig,l+1)
-     &                       - F_sed(ig,l)*wgt_SA(ig,l))
-     &                       / m_lay(ig,l)
-c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
-             IF (zqi_sa(ig,l).LT.0.0D0) THEN
-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
-c              de tps de la phys, microphys et chimie
-c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
-c              mais en fait on l epuise pour un pdt<pdtphys
-               zqi_sa(ig,l) = 0.0D0
-             ENDIF
-            pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)                       
-         ENDDO
-      ENDDO
-
-c     Partie H2Ol
-c     ~~~~~~~~~~~
-                     
-      DO l = 1, n_lev
-         DO ig=1,n_lon
-            zqi_wv(ig,l) = zqi_wv(ig,l) + (
-     &                         F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
-     &                       - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
-     &                       / m_lay(ig,l)
-c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
-             IF (zqi_wv(ig,l).LT.0.0D0) THEN
-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
-c              de tps de la phys, microphys et chimie
-c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
-c              mais en fait on l epuise pour un pdt<pdtphys
-               zqi_wv(ig,l) = 0.0D0
-             ENDIF
-            pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)                   
-         ENDDO
-      ENDDO
-
-c		Save output file in 1D model
-c		============================ 
-	        
-c      IF (n_lon .EQ. 1) THEN
-c      PRINT*,'Save output sedim' 	
-c      DO l = 1, n_lev
-c      	DO ig=1,n_lon
-c      	WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
-c     & 	(WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
-c     & 	(1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
-c     & 	pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
-c      ENDDO
-c      	ENDDO
-c      ENDIF   
-
-      RETURN
-      END
-
-*******************************************************************************
-      REAL FUNCTION VISCOSITY_CO2(temp)
-c	Aurélien Stolzenbach 2015
-c	Calcul de la viscosité dynamique du CO2 80°K -> 300°K
-c	Viscosité dynamique en Pa.s
-c	Source: Johnston & Grilly (1942)
-
-c	température en °K
-	REAL, INTENT(IN) :: temp
-	
-	REAL :: denom, numer
-	
-c	Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
-	
-	numer = 200.**(2.27/4.27)-0.435
-	denom = temp**(2.27/4.27)-0.435
-	
-	VISCOSITY_CO2 = (numer/denom)*1015.*(temp/200.)**(3./2.)
-
-c	convertion de Poises*1e7 -> Pa.s	
-	VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8	
-
-      END FUNCTION VISCOSITY_CO2
-*******************************************************************************
-
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_venus.F	(revision 1666)
+++ 	(revision )
@@ -1,407 +1,0 @@
-!******************************************************************************
-!*     venus_SAS_composition SUBROUTINE
-!*     modified from
-!*     PROGRAM PSC_MODEL_E
-!*     by A. Määttänen
-!*     subroutine for LMDZ+photochemistry VENUS
-!*     by A. Stolzenbach
-!*
-!*     Input/Output files:
-!*     -------------------
-!*
-!----------------------------------------------------------------------------
-      SUBROUTINE new_cloud_venus(
-     + nblev, nblon,
-     + TT,PP,
-     + mrt_wv,mrt_sa,
-     + mr_wv,mr_sa)
-
-      USE chemparam_mod
-      IMPLICIT NONE
-
-#include "YOMCST.h"
-
-      INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
-      INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales
-
-!----------------------------------------------------------------------------
-!     Ambient air state variables:
-      REAL, INTENT(IN), DIMENSION(nblon,nblev) :: mrt_wv,mrt_sa,
-     +                                            TT,PP
-      REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: mr_wv,mr_sa
-!----------------------------------------------------------------------------
-      INTEGER :: ilon, ilev, imode
-!----------------------------------------------------------------------------
-!     Thermodynamic functions:
-      REAL :: ROSAS
-!----------------------------------------------------------------------------
-!     Auxilary variables:
-      REAL :: NH2SO4,NH2O
-      REAL :: H2SO4_liq,H2O_liq
-      REAL :: CONCM
-      REAL :: MCONDTOT
-      REAL :: RMODE
-      REAL :: WSAFLAG
-      REAL :: K_SAV
-!----------------------------------------------------------------------------
-!     Ridder's Method variables:
-      REAL :: WVMIN, WVMAX, WVACC
-
-      INTEGER :: NBROOT
-
-      INTEGER :: MAXITE
-      PARAMETER(MAXITE=20)
-
-      INTEGER :: NBRAC
-      PARAMETER(NBRAC=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:
-
-!AM Venus
-! These aerosols will then be given an equilibrium composition for the given size distribution
-
-  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
-  ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
-  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
-  !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
-
-!===========================================
-!     Debut boucle sur niveau et lat,lon
-!===========================================
-!     Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0
-      NBRTOT(:,:,:)=0.0E+0
-      WH2SO4(:,:)=0.0E+0
-      rho_droplet(:,:)=0.0E+0
-
-      DO ilev=cloudmin, cloudmax
-      DO ilon=1, nblon
-
-!       Boucle sur les modes
-        RMODE=0.0E+0
-        K_SAV = 0.0
-
-        DO imode=1, nbr_mode
-          IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
-!       RMODE est le rayon modal de la distribution en volume du mode le plus
-!       representatif pour la Mtot
-            RMODE=R_MEDIAN(ilon,ilev,imode)*
-     &      EXP(2.*(DLOG(STDDEV(ilon,ilev,imode))**2.))
-            K_SAV=K_MASS(ilon,ilev,imode)
-          ENDIF
-        ENDDO ! FIN boucle imode
-
-!       Initialisation des bornes pour WV
-        WVMIN=1.E-90
-        WVMAX=mrt_wv(ilon,ilev)
-
-!       Accuracy de WVeq
-        WVACC=WVMAX*1.0E-3
-
-!       BRACWV borne la fonction f(WV) - WV = 0
-!       de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
-!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
-!       Elle fait appel à la fct/ssrtine ITERWV()
-
-        CALL BRACWV(WVMIN,WVMAX,NBRAC,RMODE,
-     &  mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),TT(ilon,ilev),
-     &  PP(ilon,ilev),FLAG,WSAFLAG,NBROOT)
-
-        SELECT CASE(FLAG)
-
-        CASE(1)
-!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
-!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
-!       Elle fait appel  la fct/ssrtine ITERWV()
-
-          WH2SO4(ilon,ilev)=IRFRMWV(WVMIN,WVMAX,WVACC,MAXITE,RMODE,
-     &    TT(ilon,ilev),PP(ilon,ilev),
-     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
-
-          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
-
-!          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
-!            PRINT*,'PROBLEM RHO_DROPLET'
-!            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
-!            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
-!            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
-!            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
-!            STOP
-!          ENDIF
-
-          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
-
-	      NH2SO4=mrt_sa(ilon,ilev)*CONCM
-	      NH2O=mrt_wv(ilon,ilev)*CONCM
-
-          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
-     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
-     &       H2SO4_liq,H2O_liq,MCONDTOT)
-
-!       Boucle sur les modes
-! *********************************************
-!	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)/(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
-!	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
-          H2O_liq=H2O_liq/CONCM
-          H2SO4_liq=H2SO4_liq/CONCM
-
-          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
-          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
-
-!          Problemes quand on a condense tout, on peut obtenir des -1e-24
-!		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_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
-
-
-
-        CASE(2)
-!       Cas NROOT=0 mais proche de 0
-
-          WH2SO4(ilon,ilev)=WSAFLAG
-
-          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(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
-!     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 ROSAS (cf test externe en Python), sinon, la valeur est trop
-!     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 bonne valeur (tests externes Python confirment)
-
-          IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
-     &      (WH2SO4(ilon,ilev).GT.0.1))THEN
-            PRINT*,'PROBLEM RHO_DROPLET'
-            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
-            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
-            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
-            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
-            STOP
-          ENDIF
-
-
-          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
-
-	    NH2SO4=mrt_sa(ilon,ilev)*CONCM
-	    NH2O=mrt_wv(ilon,ilev)*CONCM
-
-          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
-     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
-     &       H2SO4_liq,H2O_liq,MCONDTOT)
-
-!       Boucle sur les modes
-! *********************************************
-!	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)/(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
-!	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
-          H2O_liq=H2O_liq/CONCM
-          H2SO4_liq=H2SO4_liq/CONCM
-
-          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
-          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
-
-!          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_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
-
-        CASE(3)
-!         Cas 0 NROOT
-            mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
-            mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
-            rho_droplet(ilon,ilev)=0.0E+0
-            WH2SO4(ilon,ilev)=0.0E+0
-            DO imode=1, nbr_mode
-              NBRTOT(ilon,ilev,imode)=0.0E+0
-            ENDDO
-
-        END SELECT
-      ENDDO   !FIN boucle ilon
-      ENDDO   !FIN boucle ilev
-
-      END SUBROUTINE new_cloud_venus
-
-
-!******************************************************************
-      SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
-     + T,H2SO4COND,H2OCOND,RMTOT)
-
-!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
-!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
-!                                       SIZE DISTRIBTUION
-!     ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
-!    ---------------------------------------------------------------
-!     INPUT:
-!     H2SO4: #/m3 of total H2SO4
-!	H2O  : #/m3 of total H2O
-!     WSA: aerosol H2SO4 weight fraction (fraction)
-!     DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume)
-!       for total mass, almost same result with ro=1.67 gr/cm3
-!     RSTDEV: standard deviation of aerosol distribution (no unit)
-!     RADIUS: MEDIAN radius (m)
-!     T: temperature (K)
-!
-!     OUTPUT:
-!     RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
-!            mais rho_droplet et M_tot_distrib doivent etre de meme dimension
-!	H2OCOND
-!	H2SO4COND
-
-
-
-      IMPLICIT NONE
-
-      REAL, INTENT(IN) :: H2SO4, H2O, WSA
-      REAL, INTENT(IN) :: DENSO4, T
-      REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
-!     working variables
-      REAL :: RMH2S4
-      REAL :: DND2,pstand,lpar,acidps
-      REAL :: x1, satpacid
-      REAL , DIMENSION(2):: act
-!
-!     masse of an H2SO4 molecule (kg)
-      RMH2S4=98.078/(6.02214129E+26)
-
-      pstand=1.01325E+5 !Pa  1 atm pressure
-
-        x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
-
-        call zeleznik(x1,t,act)
-
-!pure acid satur vapor pressure
-        lpar= -11.695+DLOG(pstand) ! Zeleznik
-        acidps=1/360.15-1.0/t+0.38/545.
-     + *(1.0+DLOG(360.15/t)-360.15/t)
-        acidps = 10156.0*acidps +lpar
-        acidps = DEXP(acidps)    !Pa
-
-!acid sat.vap.PP over mixture (flat surface):
-        satpacid=act(2)*acidps ! Pa
-
-!       Conversion from Pa to N.D #/m3
-        DND2=satpacid/(1.3806488E-23*T)
-!       Conversion from N.D #/m3 TO #/cm3
-!        DND2=DND2*1.d-6
-
-!	H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
-	IF (H2SO4.GE.DND2) THEN
-	H2SO4COND=H2SO4-DND2
-!	calcul de H2O cond correspondant a H2SO4 cond
-	H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
-
-!     RMTOT: = Mass of H2SO4 satur per 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)
-
-!	Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
-	ELSE
-	H2SO4COND=0.0E+0
-	H2OCOND=0.0E+0
-	RMTOT=0.0E+0
-	END IF
-
-!	Test si H2O en defaut H2Ocond>H2O dispo
-	IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
-
-!     Si H2O en dfaut, on as pas le bon WSA!
-!     En effet, normalement, on a exactement le WSA correspondant a
-!     WVg + WVl = WVtot
-!     Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou
-!     quand T est grand (Lower Haze), le modle reprsente mal le WSA
-!     cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA
-
-!      PRINT*,'PROBLEM H2O EN DEFAUT'
-!      PRINT*,'H2OCOND',H2OCOND,'H2O',H2O
-!      PRINT*,'WSA',WSA,'RHO',DENSO4
-!      STOP
-
-
-!	On peut alors condenser tout le H2O dispo
-	H2OCOND=H2O
-!	On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
-	H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
-
-!     RMTOT: = Mass of H2SO4 satur per 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
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h	(revision 1666)
+++ trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h	(revision 1667)
@@ -110,10 +110,10 @@
       call histwrite_phy(nid_mth,.false.,"dtlwr",itau_w,dtlw)
 c en K/s      
-      call histwrite_phy(nid_mth,.false.,"dtlwrNLTE",itau_w,d_t_nlte)
-      call histwrite_phy(nid_mth,.false.,"dtlwrLTE",itau_w,-1.*cool)
-      call histwrite_phy(nid_mth,.false.,"dteuv",itau_w,d_t_euv)
-      call histwrite_phy(nid_mth,.false.,"dtcond",itau_w,d_t_conduc)
-      call histwrite_phy(nid_mth,.false.,"dumolvis",itau_w,d_u_molvis)
-      call histwrite_phy(nid_mth,.false.,"dvmolvis",itau_w,-1.*d_v_molvis)
+c     call histwrite_phy(nid_mth,.false.,"dtlwrNLTE",itau_w,d_t_nlte)
+c     call histwrite_phy(nid_mth,.false.,"dtlwrLTE",itau_w,-1.*cool)
+c     call histwrite_phy(nid_mth,.false.,"dteuv",itau_w,d_t_euv)
+c     call histwrite_phy(nid_mth,.false.,"dtcond",itau_w,d_t_conduc)
+c     call histwrite_phy(nid_mth,.false.,"dumolvis",itau_w,d_u_molvis)
+c     call histwrite_phy(nid_mth,.false.,"dvmolvis",itau_w,-1.*d_v_molvis)
 
 c     call histwrite_phy(nid_mth,.false.,"dtec",itau_w,d_t_ec)
