SUBROUTINE new_cloud_sedim(n_lon,n_lev,ptimestep,
     &                pmidlay,pbndlay, 
     &                pt,
     &                pq, pdqsed,pdqs_sed,nq,F_sed)

      USE ioipsl
      USE infotrac
      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
               PRINT*,'STOP sedimentation on epuise tout le VMR present'
               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
               PRINT*,'STOP sedimentation on epuise tout le VMR present'
               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
*******************************************************************************