Index: /trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 521)
+++ /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 522)
@@ -1,7 +1,9 @@
-       SUBROUTINE watercloud(ngrid,nlay, ptimestep, 
+       SUBROUTINE watercloud(ngrid,nlay,ptimestep, 
      &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
      &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
      &                rsedcloud,rhocloud)
+! to use  'getin'
+      USE ioipsl_getincom
       IMPLICIT NONE
 
@@ -13,8 +15,11 @@
 c    - An improved microphysical scheme (see improvedclouds.F)
 c
+c  There is a time loop specific to cloud formation and sedimentation
+c  due to timescales smaller than the GCM integration timestep.
+c
 c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 
-c           J.-B. Madeleine
+c           J.-B. Madeleine, Thomas Navarro
 c
-c  2004 - Oct. 2011
+c  2004 - 2012
 c=======================================================================
 
@@ -35,28 +40,28 @@
 
       INTEGER ngrid,nlay
-      integer nq                 ! nombre de traceurs 
+      INTEGER nq                 ! nombre de traceurs 
       REAL ptimestep             ! pas de temps physique (s)
       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      REAL pdpsrf(ngrid)         ! tendance surf pressure
+      REAL pdpsrf(ngrid)         ! tendence surf pressure
       REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
       REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
       REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
-      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
+      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
 
       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
-      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
-
-      REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
-      REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
-      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
+      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
+
+      REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
+      REAL tauscaling(ngridmx)   ! Convertion factor for dust amount
+      real rdust(ngridmx,nlay)   ! Dust geometric mean radius (m)
 
 c   Outputs:
 c   -------
 
-      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
+      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
       real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
-      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
-                                   !   a la chaleur latente
+      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
+                                   ! a la chaleur latente
 
       REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
@@ -64,11 +69,42 @@
       REAL nuice(ngrid,nlay)   ! Estimated effective variance
                                !   of the size distribution
-      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
-      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
+      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
+      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
 
 c   local:
 c   ------
 
-      INTEGER ig,l
+      ! for sedimentation
+      REAL zq(ngridmx,nlay,nqmx)    ! local value of tracers
+      real masse (ngridmx,nlay)     ! Layer mass (kg.m-2)
+      real epaisseur (ngridmx,nlay) ! Layer thickness (m)
+      real wq(ngridmx,nlay+1)       ! displaced tracer mass (kg.m-2)
+      
+      ! for ice radius computation
+      REAL Mo,No
+      REAl ccntyp
+      real beta      ! correction for the shape of the ice particles (cf. newsedim)
+      save beta
+      
+      ! for time loop
+      INTEGER microstep  ! time subsampling step variable
+      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
+      SAVE imicro
+      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
+      SAVE microtimestep
+      
+      ! tendency given by clouds (inside the micro loop)
+      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
+      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
+
+      ! global tendency (clouds+sedim+physics)
+      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
+      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
+      
+      REAL CBRT
+      EXTERNAL CBRT
+
+      
+      INTEGER iq,ig,l
       LOGICAL,SAVE :: firstcall=.true.
 
@@ -93,21 +129,298 @@
         write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
         write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
+                
+
+        write(*,*) "correction for the shape of the ice particles ?"
+        beta=0.75 ! default value
+        call getin("ice_shape",beta)
+        write(*,*) " ice_shape = ",beta
+        
+        write(*,*) "time subsampling for microphysic ?"
+        imicro = 1
+        call getin("imicro",imicro)
+        write(*,*)"imicro = ",imicro
+        
+        microtimestep = ptimestep/float(imicro)
+        write(*,*)"Physical timestep is",ptimestep 
+        write(*,*)"Microphysics timestep is",microtimestep 
 
         firstcall=.false.
       ENDIF ! of IF (firstcall)
-
-c     Main call to the different cloud schemes:
-      IF (microphys) THEN
-        CALL improvedclouds(ngrid,nlay,ptimestep,
-     &             pplev,pplay,pt,pdt,
-     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
+      
+c-----Initialization
+      subpdq(1:ngrid,1:nlay,1:nq) = 0
+      subpdt(1:ngrid,1:nlay)      = 0
+      pdqscloud(1:ngrid,1:nq)     = 0
+      zq(1:ngrid,1:nlay,1:nq)     = 0
+
+c-----Computing the different layer properties for clouds sedimentation      
+      do l=1,nlay
+        do ig=1, ngrid
+          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
+         enddo
+      enddo
+
+c------------------------------------------------------------------
+c Time subsampling for coupled microphysics and sedimentation
+c------------------------------------------------------------------
+      DO microstep=1,imicro 
+      
+c-------------------------------------------------------------------
+c   1.  Tendencies: 
+c------------------
+
+
+c------ Temperature tendency subpdt
+        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
+        ! If imicro=1 subpdt is the same as pdt
+        DO l=1,nlay
+          DO ig=1,ngrid
+             subpdt(ig,l) = subpdt(ig,l)
+     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
+          ENDDO
+        ENDDO
+c------ Traceurs tendencies subpdq
+c------ At each micro timestep we add pdq in order to have a stepped entry
+        IF (microphys) THEN
+          DO l=1,nlay
+            DO ig=1,ngrid
+              subpdq(ig,l,igcm_dust_mass) = 
+     &            subpdq(ig,l,igcm_dust_mass)
+     &          + pdq(ig,l,igcm_dust_mass)
+              subpdq(ig,l,igcm_dust_number) = 
+     &            subpdq(ig,l,igcm_dust_number)
+     &          + pdq(ig,l,igcm_dust_number)
+              subpdq(ig,l,igcm_ccn_mass) = 
+     &            subpdq(ig,l,igcm_ccn_mass)
+     &          + pdq(ig,l,igcm_ccn_mass)
+              subpdq(ig,l,igcm_ccn_number) = 
+     &            subpdq(ig,l,igcm_ccn_number)
+     &          + pdq(ig,l,igcm_ccn_number)
+            ENDDO
+          ENDDO
+        ENDIF
+        DO l=1,nlay
+          DO ig=1,ngrid
+            subpdq(ig,l,igcm_h2o_ice) = 
+     &          subpdq(ig,l,igcm_h2o_ice)
+     &        + pdq(ig,l,igcm_h2o_ice)
+            subpdq(ig,l,igcm_h2o_vap) = 
+     &          subpdq(ig,l,igcm_h2o_vap)
+     &        + pdq(ig,l,igcm_h2o_vap)
+          ENDDO
+        ENDDO
+        
+        
+c-------------------------------------------------------------------
+c   2.  Main call to the different cloud schemes:
+c------------------------------------------------
+        IF (microphys) THEN
+           CALL improvedclouds(ngrid,nlay,microtimestep,
+     &             pplev,pplay,pzlev,pt,subpdt,
+     &             pq,subpdq,subpdqcloud,subpdtcloud,
      &             nq,tauscaling,rdust,rice,nuice,
      &             rsedcloud,rhocloud)
-      ELSE
-        CALL simpleclouds(ngrid,nlay,ptimestep,
-     &             pplev,pplay,pzlev,pzlay,pt,pdt,
-     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
+        ELSE
+           CALL simpleclouds(ngrid,nlay,microtimestep,
+     &             pplev,pplay,pzlev,pzlay,pt,subpdt,
+     &             pq,subpdq,subpdqcloud,subpdtcloud,
      &             nq,tau,rice,nuice,rsedcloud)
-      ENDIF
+        ENDIF
+      
+c--------------------------------------------------------------------
+c    3.  CCN & ice sedimentation:
+c--------------------------------
+! Sedimentation is done here for water ice and its CCN only
+! callsedim in physics is done for dust (and others species if any)           
+                   
+        DO l=1,nlay
+         DO ig=1,ngrid
+          subpdt(ig,l) =
+     &      subpdt(ig,l) + subpdtcloud(ig,l)
+         ENDDO
+        ENDDO
+        
+c------- water ice update before sedimentation (radius is done in the cloud scheme routine)    
+         DO l=1,nlay
+          DO ig=1, ngrid
+          zq(ig,l,igcm_h2o_ice)= max(1e-30,
+     &       pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice)
+     &       + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep)
+!          zq(ig,l,igcm_h2o_vap)= max(1e-30,
+!     &       pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap)
+!     &       + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep)
+          ENDDO
+         ENDDO
+        
+      
+c------- CCN update before sedimentation  
+        IF (microphys) THEN
+         DO l=1,nlay
+          DO ig=1,ngrid
+           zq(ig,l,igcm_ccn_number)=
+     &       pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number)
+     &       + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
+           zq(ig,l,igcm_ccn_number)=  max(1e-30,
+     &       zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ?
+           zq(ig,l,igcm_ccn_mass)=
+     &       pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass)
+     &       + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
+           zq(ig,l,igcm_ccn_mass)=max(1e-30,
+     &       zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ?
+          ENDDO
+         ENDDO
+        ENDIF
+        
+
+        
+        IF (microphys) THEN
+        
+c------- CCN number sedimentation
+          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
+     &        rhocloud,zq(1,1,igcm_ccn_number),wq,beta)
+          do ig=1,ngrid
+            ! matters if one would like to know ccn surface deposition
+            pdqscloud(ig,igcm_ccn_number)=
+     &          pdqscloud(ig,igcm_ccn_number) + wq(ig,1)
+          enddo
+          
+c------- CCN mass sedimentation
+          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
+     &        rhocloud,zq(1,1,igcm_ccn_mass),wq,beta)
+          do ig=1,ngrid
+            ! matters if one would like to know ccn surface deposition
+            pdqscloud(ig,igcm_ccn_mass)=
+     &          pdqscloud(ig,igcm_ccn_mass) + wq(ig,1)
+          enddo
+          
+c------- Water ice sedimentation -- improved microphys
+          call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
+     &        rhocloud,zq(1,1,igcm_h2o_ice),wq,beta)
+        ELSE
+        
+c------- Water ice sedimentation -- simple microphys
+          call newsedim(ngrid,nlay,ngrid*nlay,1,
+     &        microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
+     &        rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta)
+     
+        ENDIF
+
+
+c------- Surface ice tendency update
+        DO ig=1,ngrid
+          pdqscloud(ig,igcm_h2o_ice)=
+     &          pdqscloud(ig,igcm_h2o_ice) + wq(ig,1)
+        ENDDO
+      
+
+c-------------------------------------------------------------------
+c   5.  Updating tendencies after sedimentation:
+c-----------------------------------------------
+
+        DO l=1,nlay
+         DO ig=1,ngrid
+           
+           subpdq(ig,l,igcm_h2o_ice) = 
+     &     (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice))
+     &      /microtimestep
+
+     
+           subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap)
+     &        +subpdqcloud(ig,l,igcm_h2o_vap)
+     
+         ENDDO
+        ENDDO
+      
+      
+        IF (microphys) then
+         DO l=1,nlay
+          DO ig=1,ngrid
+           subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
+     &        -pq(ig,l,igcm_ccn_number))/microtimestep
+           subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
+     &        -pq(ig,l,igcm_ccn_mass))/microtimestep
+          ENDDO
+         ENDDO
+        ENDIF
+        
+
+
+      ENDDO ! of DO microstep=1,imicro
+      
+c-------------------------------------------------------------------
+c   6.  Compute final tendencies after time loop:
+c------------------------------------------------
+
+c------ Whole temperature tendency pdtcloud
+       DO l=1,nlay
+         DO ig=1,ngrid
+             pdtcloud(ig,l) =
+     &         subpdt(ig,l)/imicro-pdt(ig,l)
+          ENDDO
+       ENDDO
+c------ Traceurs tendencies pdqcloud
+       DO iq=1,nq
+        DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro
+     &       - pdq(ig,l,iq)
+         ENDDO
+        ENDDO
+       ENDDO
+c------ Traceurs surface tendencies pdqscloud
+       DO iq=1,nq
+        DO ig=1,ngrid
+           pdqscloud(ig,iq) =
+     &         pdqscloud(ig,iq)/ptimestep
+        ENDDO
+       ENDDO
+       
+
+          
+c------Update the ice particle size "rice" for output or photochemistry.
+c------It is not used again in the water cycle.
+       IF(scavenging) THEN 
+        DO l=1, nlay
+         DO ig=1,ngrid
+          Mo = pq(ig,l,igcm_h2o_ice) 
+     &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
+     &         + (pq(ig,l,igcm_ccn_mass)
+     &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
+     &         *tauscaling(ig) + 1.e-30
+          No = (pq(ig,l,igcm_ccn_number)
+     &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
+     &         *tauscaling(ig) + 1.e-30
+          rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) +  
+     &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
+     &                   / Mo * rho_ice
+     &                  + (pq(ig,l,igcm_ccn_mass)
+     &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
+     &                   *tauscaling(ig)/ Mo * rho_dust
+          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
+          rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
+           if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8
+         ENDDO
+        ENDDO
+       ELSE
+        DO l=1,nlay
+          DO ig=1,ngrid
+            ccntyp = 
+     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
+            ccntyp = ccntyp /ccn_factor
+            rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
+     &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
+     &      +ccntyp*(4./3.)*pi*rdust(ig,l)**3.)
+     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
+          ENDDO
+        ENDDO
+       ENDIF ! of IF(scavenging)
+                
+!--------------------------------------------------------------
+!--------------------------------------------------------------
+      
 
 c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
@@ -122,4 +435,33 @@
       end do
 
+c=======================================================================
+
+!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice
+!!      if (photochem) then
+!!c        computation of dust and ice surface area (micron2/cm3)
+!!c        for heterogeneous chemistry
+!!
+!!         do l = 1,nlay
+!!            do ig = 1,ngrid
+!!c                                                                      
+!!c              npart: number density of ccn in #/cm3    
+!!c                                                     
+!!               npart(ig,l) = 1.e-6*ccn(ig,l)  
+!!     $                       *masse(ig,l)/epaisseur(ig,l) 
+!!c
+!!c              dust and ice surface area
+!!c                                                                  
+!!               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
+!!c                                                                  
+!!               if (rice(ig,l) .ge. rdust(ig,l)) then                   
+!!                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
+!!                  surfdust(ig,l) = 0.
+!!               else                                                    
+!!                  surfice(ig,l) = 0.                                 
+!!               end if                                             
+!!            end do      ! of do ig=1,ngrid
+!!         end do         ! of do l=1,nlay
+!!      end if            ! of photochem
+
       RETURN
       END
