Index: trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 2934)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 2966)
@@ -27,4 +27,5 @@
      &                      qperemin
       use dimradmars_mod, only: naerkind
+      use conc_mod, only: mmean
       use write_output_mod, only: write_output
       IMPLICIT NONE
@@ -45,4 +46,8 @@
 c
 c  2004 - 2012
+c
+c     2023: J. Naar, adding different subtimestep on each grid cell
+c          plus not doing microphysics if no water present
+c          plus simpleclouds no longer in the loop for microphysics
 c=======================================================================
 
@@ -105,4 +110,5 @@
 
       ! global tendency (clouds+physics)
+      ! JN : keeping this for simpleclouds scheme
       REAL sum_subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
       REAL sum_subpdt(ngrid,nlay)         ! cf. pdtcloud
@@ -137,4 +143,17 @@
       
       !$OMP THREADPRIVATE(flagcloud)
+ 
+! Scheme for adaptative timestep J. Naar 2023
+c      LOGICAL :: computed_micro(ngrid,nlay) ! Check if microphy was done in this cell
+      REAL :: computed_micro(ngrid,nlay) ! Check if microphy was done inthis cell (logical)
+      REAL :: zt_micro(ngrid,nlay) ! Temperature during microphysics (K)
+      REAL :: zq_micro(ngrid,nlay,nq) ! Tracers during microphysics (kg/kg) 
+      REAL :: zqsat_micro(ngrid,nlay) ! Theoritical q(h2o_vap) at saturation (kg/kg)
+      INTEGER :: zimicro(ngrid,nlay)    ! Number of ptimestep divisions
+      REAL :: zpotcond_inst(ngrid,nlay) ! condensable water at the beginning of physics (kg/kg)
+      REAL :: zpotcond_full(ngrid,nlay) ! condensable water with integrated tendancies (kg/kg)
+      REAL :: zpotcond(ngrid,nlay) ! maximal condensable water, used to
+compute adaptative subdivision of ptimestep
+
 
 c    ** un petit test de coherence
@@ -282,69 +301,56 @@
 c------------------------------------------------------------------ 
       rhocloud(1:ngrid,1:nlay) = rho_dust
-      DO microstep=1,imicro 
-      
+
+
+c     Initialisation of all the stuff JN2023
+c      computed_micro(1:ngrid,1:nlay)=.false.
+      computed_micro(1:ngrid,1:nlay)=0.
+      zt_micro(:,:)=pt(:,:)
+      zq_micro(:,:,:)=pq(:,:,:)
+      zq_micro(:,:,:)=pq(:,:,:)
+      call watersat(ngrid*nlay,zt_micro,pplay,zqsat_micro)
+      zpotcond_inst=zq_micro(:,:,igcm_h2o_vap) - zqsat_micro
+      call watersat(ngrid*nlay,zt_micro+pdt*ptimestep,pplay,zqsat_micro)
+      zpotcond_full=(zq_micro(:,:,igcm_h2o_vap)+
+     &             pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat_micro
+      zimicro(1:ngrid,1:nlay)=imicro
+      if (cloud_adapt_ts) then
+              call adapt_imicro(ptimestep,zpotcond(ig,l),
+     $                   zimicro(ig,l))
+      endif! (cloud_adapt_ts) then
+      DO l=1,nlay
+        DO ig=1,ngrid
+c         Start by computing the condensable water vapor amount
+          if (zpotcond_full(ig,l).gt.0.) then
+            zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
+          else if (zpotcond_full(ig,l).le.0.) then
+            zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
+          endif! (zpotcond_full.gt.0.) then
+          microtimestep=ptimestep/real(zimicro(ig,l))
+c         Check if microphysics is even needed, that is if enough action
+c         is happening water-wise
+          if ((pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
+     &      .gt.1e-22) .or. (abs(zpotcond(ig,l)).gt.1e-22)) then
+c         Eventuellement sortir simpleclouds de la boucle egalement
+          computed_micro(ig,l)=1.
+          DO microstep=1,zimicro(ig,l) 
+      
+
+c JN : incrementing after main microphysics scheme
+c Previously we were incrementing tendencies, we now
+c increment tracers and temperature directly
+c We are thus starting at the end of the first iteration
+c 
 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
-             sum_subpdt(ig,l) = sum_subpdt(ig,l)
-     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
-          ENDDO
-        ENDDO
-c------ Tracers tendencies subpdq are additionned
-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
-              sum_subpdq(ig,l,igcm_dust_mass) = 
-     &            sum_subpdq(ig,l,igcm_dust_mass)
-     &          + pdq(ig,l,igcm_dust_mass)
-              sum_subpdq(ig,l,igcm_dust_number) = 
-     &            sum_subpdq(ig,l,igcm_dust_number)
-     &          + pdq(ig,l,igcm_dust_number)
-              sum_subpdq(ig,l,igcm_ccn_mass) = 
-     &            sum_subpdq(ig,l,igcm_ccn_mass)
-     &          + pdq(ig,l,igcm_ccn_mass)
-              sum_subpdq(ig,l,igcm_ccn_number) = 
-     &            sum_subpdq(ig,l,igcm_ccn_number)
-     &          + pdq(ig,l,igcm_ccn_number)
-            ENDDO
-          ENDDO
-        ENDIF
-        DO l=1,nlay
-          DO ig=1,ngrid
-            sum_subpdq(ig,l,igcm_h2o_ice) = 
-     &          sum_subpdq(ig,l,igcm_h2o_ice)
-     &        + pdq(ig,l,igcm_h2o_ice)
-            sum_subpdq(ig,l,igcm_h2o_vap) = 
-     &          sum_subpdq(ig,l,igcm_h2o_vap)
-     &        + pdq(ig,l,igcm_h2o_vap)
-            IF (hdo) THEN
-            sum_subpdq(ig,l,igcm_hdo_ice) = 
-     &          sum_subpdq(ig,l,igcm_hdo_ice)
-     &        + pdq(ig,l,igcm_hdo_ice)
-            sum_subpdq(ig,l,igcm_hdo_vap) = 
-     &          sum_subpdq(ig,l,igcm_hdo_vap)
-     &        + pdq(ig,l,igcm_hdo_vap)
-            ENDIF !hdo
-          ENDDO
-        ENDDO      
-        
-c-------------------------------------------------------------------
-c   2.  Main call to the different cloud schemes:
+c   1.  Main call to the different cloud schemes:
 c------------------------------------------------
         IF (microphys) THEN
-           CALL improvedclouds(ngrid,nlay,microtimestep,
-     &             pplay,pteff,sum_subpdt, 
-     &             pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
-     &             nq,tauscaling)
+           CALL improvedclouds(microtimestep,
+     &          pplay(ig,l),zt_micro(ig,l), 
+     &          zq_micro(ig,l,:),subpdqcloud(ig,l,:),
+     &          subpdtcloud(ig,l),nq,tauscaling(ig),mmean(ig,l))
 
         ELSE
+c Simpleclouds should maybe be taken out and put in a specific loop ?
            CALL simpleclouds(ngrid,nlay,microtimestep,
      &             pplay,pzlay,pteff,sum_subpdt,
@@ -354,55 +360,51 @@
 
 c-------------------------------------------------------------------
-c   3.  Updating tendencies after cloud scheme:
+c   2.  Updating tracers and temperature after cloud scheme:
 c-----------------------------------------------
 
         IF (microphys) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              sum_subpdq(ig,l,igcm_dust_mass) =
-     &            sum_subpdq(ig,l,igcm_dust_mass)
-     &          + subpdqcloud(ig,l,igcm_dust_mass)
-              sum_subpdq(ig,l,igcm_dust_number) =
-     &            sum_subpdq(ig,l,igcm_dust_number)
-     &          + subpdqcloud(ig,l,igcm_dust_number)
-              sum_subpdq(ig,l,igcm_ccn_mass) =
-     &            sum_subpdq(ig,l,igcm_ccn_mass)
-     &          + subpdqcloud(ig,l,igcm_ccn_mass)
-              sum_subpdq(ig,l,igcm_ccn_number) =
-     &            sum_subpdq(ig,l,igcm_ccn_number)
-     &          + subpdqcloud(ig,l,igcm_ccn_number)
-            ENDDO
-          ENDDO
+              zq_micro(ig,l,igcm_dust_mass) =
+     &         zq_micro(ig,l,igcm_dust_mass)+(pdq(ig,l,igcm_dust_mass)
+     &         +subpdqcloud(ig,l,igcm_dust_mass))*microtimestep
+              zq_micro(ig,l,igcm_dust_number) =
+     &         zq_micro(ig,l,igcm_dust_number)
+     &         +(pdq(ig,l,igcm_dust_number)
+     &         + subpdqcloud(ig,l,igcm_dust_number))*microtimestep
+              zq_micro(ig,l,igcm_ccn_mass) =
+     &         zq_micro(ig,l,igcm_ccn_mass) +
+     &         (pdq(ig,l,igcm_ccn_mass)
+     &         +subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
+              zq_micro(ig,l,igcm_ccn_number) =
+     &          zq_micro(ig,l,igcm_ccn_number) +
+     &         (pdq(ig,l,igcm_ccn_number)
+     &          + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
         ENDIF
-        DO l=1,nlay
-          DO ig=1,ngrid
-            sum_subpdq(ig,l,igcm_h2o_ice) =
-     &          sum_subpdq(ig,l,igcm_h2o_ice)
-     &        + subpdqcloud(ig,l,igcm_h2o_ice)
-            sum_subpdq(ig,l,igcm_h2o_vap) =
-     &          sum_subpdq(ig,l,igcm_h2o_vap)
-     &        + subpdqcloud(ig,l,igcm_h2o_vap)
+            zq_micro(ig,l,igcm_h2o_ice) =
+     &       zq_micro(ig,l,igcm_h2o_ice)+
+     &         (pdq(ig,l,igcm_h2o_ice)
+     &        + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep
+            zq_micro(ig,l,igcm_h2o_vap) =
+     &       zq_micro(ig,l,igcm_h2o_vap)+
+     &         (pdq(ig,l,igcm_h2o_vap)
+     &        + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep
 
             IF (hdo) THEN
-            sum_subpdq(ig,l,igcm_hdo_ice) =
-     &          sum_subpdq(ig,l,igcm_hdo_ice)
-     &        + subpdqcloud(ig,l,igcm_hdo_ice)
-            sum_subpdq(ig,l,igcm_hdo_vap) =
-     &          sum_subpdq(ig,l,igcm_hdo_vap)
-     &        + subpdqcloud(ig,l,igcm_hdo_vap)
+            zq_micro(ig,l,igcm_hdo_ice) =
+     &       zq_micro(ig,l,igcm_hdo_ice)+
+     &         (pdq(ig,l,igcm_hdo_ice)
+     &        + subpdqcloud(ig,l,igcm_hdo_ice))*microtimestep
+            zq_micro(ig,l,igcm_hdo_vap) =
+     &       zq_micro(ig,l,igcm_hdo_vap)+
+     &         (pdq(ig,l,igcm_hdo_vap)
+     &        + subpdqcloud(ig,l,igcm_hdo_vap))*microtimestep
             ENDIF ! hdo
 
-          ENDDO
-        ENDDO
-        
+c  Could also set subpdtcloud to 0 if not activice to make it simpler
+            zt_micro(ig,l) = zt_micro(ig,l)+
+     &           pdt(ig,l)*microtimestep
         IF (activice) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              sum_subpdt(ig,l) =
-     &            sum_subpdt(ig,l) + subpdtcloud(ig,l)
-            ENDDO
-          ENDDO
+              zt_micro(ig,l) = zt_micro(ig,l)+
+     &           subpdtcloud(ig,l)*microtimestep
         ENDIF
-
 !      !! Example of how to use writediagmicrofi useful to 
 !      !! get outputs at each microphysical sub-timestep (better to be used in 1D)
@@ -411,8 +413,18 @@
 !     &      'subpdtcloud','K/s',1,subpdtcloud(:,:))     
  
-      ENDDO ! of DO microstep=1,imicro
-      
+          ENDDO ! of DO microstep=1,imicro
+          endif! (zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
+!     &      .gt.1e-22).or.(abs(zpotcond).gt.1e-22) then
+        ENDDO ! ig=1,ngrid
+      ENDDO ! l=1,nlay
+
+      
+c------ Useful outputs to check how it went
+      call write_output("computed_micro","computed_micro "//
+     &   "after microphysics","logical",computed_micro(:,:))
+      call write_output("zimicro","Used number of subtimestep "//
+     &   "in cloud microphysics","integer",real(zimicro(:,:)))
 c-------------------------------------------------------------------
-c   6.  Compute final tendencies after time loop:
+c   3.  Compute final tendencies after time loop:
 c------------------------------------------------
 
@@ -420,51 +432,59 @@
        DO l=1,nlay
          DO ig=1,ngrid
-             pdtcloud(ig,l) =
-     &         sum_subpdt(ig,l)/real(imicro)-pdt(ig,l)
-          ENDDO
+         pdtcloud(ig,l) = -pdt(ig,l)+
+     &        (zt_micro(ig,l)-pt(ig,l)) / ptimestep
+         ENDDO
        ENDDO
-       
+
 c------ Tracers tendencies pdqcloud
        DO l=1,nlay
          DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_h2o_ice) = 
-     &        sum_subpdq(ig,l,igcm_h2o_ice)/real(imicro) 
-     &       - pdq(ig,l,igcm_h2o_ice)
-            pdqcloud(ig,l,igcm_h2o_vap) = 
-     &        sum_subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
-     &       - pdq(ig,l,igcm_h2o_vap)
+            pdqcloud(ig,l,igcm_h2o_ice) =
+     &        -pdq(ig,l,igcm_h2o_ice)+
+     &       (zq_micro(ig,l,igcm_h2o_ice) -
+     &        pq(ig,l,igcm_h2o_ice)) / ptimestep
+            pdqcloud(ig,l,igcm_h2o_vap) =
+     &        -pdq(ig,l,igcm_h2o_vap)+
+     &       (zq_micro(ig,l,igcm_h2o_vap) -
+     &        pq(ig,l,igcm_h2o_vap)) / ptimestep
             IF (hdo) THEN
-            pdqcloud(ig,l,igcm_hdo_ice) = 
-     &        sum_subpdq(ig,l,igcm_hdo_ice)/real(imicro) 
-     &       - pdq(ig,l,igcm_hdo_ice)
-            pdqcloud(ig,l,igcm_hdo_vap) = 
-     &        sum_subpdq(ig,l,igcm_hdo_vap)/real(imicro) 
-     &       - pdq(ig,l,igcm_hdo_vap)
+            pdqcloud(ig,l,igcm_hdo_ice) =
+     &        -pdq(ig,l,igcm_hdo_ice)+
+     &       (zq_micro(ig,l,igcm_hdo_ice) -
+     &        pq(ig,l,igcm_hdo_ice)) / ptimestep
+            pdqcloud(ig,l,igcm_hdo_vap) =
+     &        -pdq(ig,l,igcm_hdo_vap)+
+     &       (zq_micro(ig,l,igcm_hdo_vap) -
+     &        pq(ig,l,igcm_hdo_vap)) / ptimestep
             ENDIF !hdo
          ENDDO
-       ENDDO
-       
+       ENDDO       
+
        IF(microphys) THEN
         DO l=1,nlay
          DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_ccn_mass) = 
-     &        sum_subpdq(ig,l,igcm_ccn_mass)/real(imicro) 
-     &       - pdq(ig,l,igcm_ccn_mass)
-            pdqcloud(ig,l,igcm_ccn_number) = 
-     &        sum_subpdq(ig,l,igcm_ccn_number)/real(imicro) 
-     &       - pdq(ig,l,igcm_ccn_number)
+            pdqcloud(ig,l,igcm_ccn_mass) =
+     &        -pdq(ig,l,igcm_ccn_mass)+
+     &       (zq_micro(ig,l,igcm_ccn_mass) -
+     &        pq(ig,l,igcm_ccn_mass)) / ptimestep
+            pdqcloud(ig,l,igcm_ccn_number) =
+     &        -pdq(ig,l,igcm_ccn_number)+
+     &       (zq_micro(ig,l,igcm_ccn_number) -
+     &        pq(ig,l,igcm_ccn_number)) / ptimestep
          ENDDO
         ENDDO
        ENDIF
-       
+
        IF(scavenging) THEN
         DO l=1,nlay
          DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_dust_mass) = 
-     &        sum_subpdq(ig,l,igcm_dust_mass)/real(imicro) 
-     &       - pdq(ig,l,igcm_dust_mass)
-            pdqcloud(ig,l,igcm_dust_number) = 
-     &        sum_subpdq(ig,l,igcm_dust_number)/real(imicro)
-     &       - pdq(ig,l,igcm_dust_number)
+            pdqcloud(ig,l,igcm_dust_mass) =
+     &        -pdq(ig,l,igcm_dust_mass)+
+     &       (zq_micro(ig,l,igcm_dust_mass) -
+     &        pq(ig,l,igcm_dust_mass)) / ptimestep
+            pdqcloud(ig,l,igcm_dust_number) =
+     &        -pdq(ig,l,igcm_dust_number)+
+     &       (zq_micro(ig,l,igcm_dust_number) -
+     &        pq(ig,l,igcm_dust_number)) / ptimestep
          ENDDO
         ENDDO
@@ -764,3 +784,21 @@
        end subroutine end_watercloud_mod
 
+
+
+       SUBROUTINE adapt_imicro(ptimestep,potcond,
+     $                     zimicro)
+
+c Pas de temps adaptatif pour les nuages
+
+      real,intent(in) :: ptimestep ! total duration of physics (sec)
+      real,intent(in) :: potcond ! total duration of physics (sec)
+      real :: alpha, beta ! total duration of physics (sec)
+      integer,intent(out) :: zimicro ! number of ptimestep division
+
+c       zimicro = ptimestep*alpha*potcond**beta
+       zimicro = 30
+
+       END SUBROUTINE adapt_imicro
+
+
       END MODULE watercloud_mod 
