MODULE improvedclouds_mod

IMPLICIT NONE

CONTAINS

!=======================================================================

SUBROUTINE improvedclouds(ngrid,nlay,ptimestep,pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro,zt,zq)

use updaterad, only: updaterice_micro, updaterccn
use watersat_mod, only: watersat
use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, &
                      igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, &
                      igcm_hdo_ice,qperemin
use conc_mod, only: mmean
use comcstfi_h, only: pi, cpp
use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp
use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To
use nuclea_mod, only: nuclea
use sig_h2o_mod, only: sig_h2o
use growthrate_mod, only: growthrate
use write_output_mod, only: write_output
use callkeys_mod, only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac

implicit none

!------------------------------------------------------------------
!  This routine is used to form clouds when a parcel of the GCM is
!    saturated. It includes the ability to have supersaturation, a
!    computation of the nucleation rates, growthrates and the
!    scavenging of dust particles by clouds.
!  It is worth noting that the amount of dust is computed using the
!    dust optical depth computed in aeropacity.F. That's why
!    the variable called "tauscaling" is used to convert
!    pq(dust_mass) and pq(dust_number), which are relative
!    quantities, to absolute and realistic quantities stored in zq.
!    This has to be done to convert the inputs into absolute
!    values, but also to convert the outputs back into relative
!    values which are then used by the sedimentation and advection
!    schemes.

!  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
!           (October 2011)
!           T. Navarro, debug,correction, new scheme (October-April 2011)
!           A. Spiga, optimization (February 2012)
!           J. Naar, adaptative subtimestep now done here (June 2023)
!------------------------------------------------------------------

!------------------------------------------------------------------
!     Inputs/outputs:
INTEGER, INTENT(IN) :: ngrid,nlay
INTEGER, INTENT(IN) :: nq ! nombre de traceurs
REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa)            
REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K)
REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s)
REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg)
REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s)
REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility)

REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg) 
REAL, dimension(ngrid,nlay), INTENT(OUT) :: zt ! temperature post microphy (K)

!------------------------------------------------------------------
!     Local variables:
LOGICAL, SAVE ::  firstcall = .true.
!$OMP THREADPRIVATE(firstcall)
REAL*8 :: derf ! Error function
!external derf
INTEGER :: ig,l,i
REAL, dimension(ngrid,nlay) :: zqsat ! saturation
REAL :: lw ! Latent heat of sublimation (J.kg-1) 
REAL :: cste
REAL :: dMice ! mass of condensed ice
REAL :: dMicetot ! JN
REAL :: dMice_hdo ! mass of condensed HDO ice
REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient
! REAL :: sumcheck
REAL*8 :: ph2o ! Water vapor partial pressure (Pa)
REAL*8 :: satu ! Water vapor saturation ratio over ice
REAL*8 :: Mo,No, Rn, Rm, dev2, n_derf, m_derf
REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin
REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin
REAL :: dN, dM, seq
REAL, dimension(nbin_cld) :: rate ! nucleation rate
REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004)
REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3)
REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m)
REAL :: res ! Resistance growth
REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient

!     Parameters of the size discretization used by the microphysical scheme
DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
DOUBLE PRECISION :: vrat_cld ! Volume ratio
DOUBLE PRECISION, dimension(nbin_cld+1), SAVE :: rb_cld ! boundary values of each rad_cld bin (m)
!$OMP THREADPRIVATE(rb_cld)
DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! width of each rad_cld bin (m)
DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! particle volume for each bin (m3)
REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions
!$OMP THREADPRIVATE(sigma_ice)

!----------------------------------      
! JN : used in subtimestep
REAL :: microtimestep! subdivision of ptimestep (s)
REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat
REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers 
! REAL, dimension(ngrid,nlay,nq) :: zt0 ! local initial value of temperature
INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep
INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls
REAL, dimension(ngrid,nlay) :: zpotcond_inst ! condensable water at the beginning of physics (kg/kg)
REAL, dimension(ngrid,nlay) :: zpotcond_full ! condensable water with integrated tendancies (kg/kg)
REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two)
REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two)
REAL :: spenttime ! time spent in while loop
REAL :: zdq ! used to compute adaptative timestep
LOGICAL :: ending_ts ! Condition to end while loop

!----------------------------------      
! TESTS
INTEGER :: countcells
LOGICAL, SAVE :: test_flag    ! flag for test/debuging outputs
!$OMP THREADPRIVATE(test_flag)
REAL, dimension(ngrid,nlay) :: satubf,satuaf, res_out

!------------------------------------------------------------------

! AS: firstcall OK absolute
IF (firstcall) THEN
!=============================================================
! 0. Definition of the size grid
!=============================================================
! rad_cld is the primary radius grid used for microphysics computation.
! The grid spacing is computed assuming a constant volume ratio
! between two consecutive bins; i.e. vrat_cld.
! vrat_cld is determined from the boundary values of the size grid: 
! rmin_cld and rmax_cld.
! The rb_cld array contains the boundary values of each rad_cld bin.
! dr_cld is the width of each rad_cld bin.

  ! Volume ratio between two adjacent bins
  ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
  ! vrat_cld = exp(vrat_cld)
  vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
  vrat_cld = exp(vrat_cld)
  write(*,*) "vrat_cld", vrat_cld

  rb_cld(1)  = rbmin_cld
  rad_cld(1) = rmin_cld
  vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
  ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld

  do i=1,nbin_cld-1
    rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
    vol_cld(i+1)  = vol_cld(i) * vrat_cld
  enddo
      
  do i=1,nbin_cld
    rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i)
    dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
  enddo
  rb_cld(nbin_cld+1) = rbmax_cld
  dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)

  print*, ' '
  print*,'Microphysics: size bin information:'
  print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
  print*,'-----------------------------------'
  do i=1,nbin_cld
    write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i)
  enddo
  write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
  print*,'-----------------------------------'

  ! we save that so that it is not computed at each timestep and gridpoint
  ! rb_cld(i) = log(rb_cld(i))  
  rb_cld(:) = log(rb_cld(:))

  ! Contact parameter of water ice on dust ( m=cos(theta) )
  !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  !  mteta is initialized in conf_phys
  write(*,*) 'water_param contact parameter:', mteta

  ! Volume of a water molecule (m3)
  vo1 = mh2o / dble(rho_ice)
  ! Variance of the ice and CCN distributions
  sigma_ice = sqrt(log(1.+nuice_sed))
      
  write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
  write(*,*) 'nuice for sedimentation:', nuice_sed
  write(*,*) 'Volume of a water molecule:', vo1

  test_flag = .false.
  firstcall=.false.
ENDIF

!=============================================================
! 1. Initialisation
!=============================================================

res_out(:,:) = 0
rice(:,:) = 1.e-8

!     Initialize the temperature, tracers
zt(:,:)=pt(:,:)
zq(:,:,:)=pq(:,:,:)
subpdtcloud(:,:)=0

WHERE( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30

!=============================================================
! 2. Compute saturation
!=============================================================

dev2 = 1. / ( sqrt(2.) * sigma_ice )

! Compute the condensable amount of water vapor or the sublimable water ice (if negative)
! Attention ici pdt*ptimestep peut etre severe
call watersat(ngrid*nlay,max(1.,zt+pdt*ptimestep),pplay,zqsat) ! DIFF: "temp+trend" at least 1
zpotcond_full=(zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
call watersat(ngrid*nlay,zt,pplay,zqsat)
zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat
zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
! zpotcond is the most strict criterion between the two 
DO l=1,nlay
  DO ig=1,ngrid
    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
    ! zpotcond(ig,l)=1.e-2
  ENDDO !ig=1,ngrid
ENDDO !l=1,nlay
      
countcells = 0
zimicro(:,:)=imicro
count_micro(:,:)=0

! Main loop over the GCM's grid
DO l=1,nlay
  DO ig=1,ngrid
    ! Subtimestep : here we go
    IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l))
    spenttime = 0.
    dMicetot= 0. ! DIFF: dMicetot=new var
    ending_ts=.false.
    DO while (.not.ending_ts)
      call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
      zqsat(ig,l)=zqsat_tmp(1)
      ! Get the partial pressure of water vapor and its saturation ratio
      ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
      satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
      microtimestep=ptimestep/real(zimicro(ig,l))
      ! if (satu.lt.1.0) then
      !   microtimestep=ptimestep/30.
      !   write(*,*) "saturation correction" !JN
      ! endif
      
      ! Initialize tracers for scavenging + hdo computations (JN)
      zq0(ig,l,:)=zq(ig,l,:)

      ! Check if we are integrating over ptimestep
      if (spenttime+microtimestep.ge.ptimestep) then
        microtimestep=ptimestep-spenttime
        ! If so : last call !
        ending_ts=.true.
      endif! (spenttime+microtimestep.ge.ptimestep) then

!=============================================================
! 3. Nucleation
!=============================================================

      IF ( satu .ge. 1. ) THEN ! if there is condensation
        call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))

        ! Expand the dust moments into a binned distribution
        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
        Rn = rdust(ig,l)
        Rn = -log(Rn) 
        Rm = Rn - 3. * sigma_ice*sigma_ice  
        n_derf = derf( (rb_cld(1)+Rn) *dev2)
        m_derf = derf( (rb_cld(1)+Rm) *dev2)
        do i = 1, nbin_cld
          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
        enddo
      
        ! sumcheck = 0
        ! do i = 1, nbin_cld
        !   sumcheck = sumcheck + n_aer(i)
        ! enddo
        ! sumcheck = abs(sumcheck/No - 1)
        ! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
        !   print*, "WARNING, No sumcheck PROBLEM"
        !   print*, "sumcheck, No",sumcheck, No
        !   print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
        !   print*, "Dust binned distribution", n_aer
        ! endif
        !        
        ! sumcheck = 0
        ! do i = 1, nbin_cld
        !   sumcheck = sumcheck + m_aer(i)
        ! enddo
        ! sumcheck = abs(sumcheck/Mo - 1)
        ! if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
        !   print*, "WARNING, Mo sumcheck PROBLEM"
        !   print*, "sumcheck, Mo",sumcheck, Mo
        !   print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
        !   print*, "Dust binned distribution", m_aer
        ! endif

        ! Get the rates of nucleation
        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
      
        dN = 0.
        dM = 0.
        do i = 1, nbin_cld
          dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
          dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
        enddo

        ! Update Dust particles
        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10) 
        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
        ! Update CCNs
        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
      
      ENDIF ! of is satu >1

!=============================================================
! 4. Ice growth: scheme for radius evolution
!=============================================================
! We trigger crystal growth if and only if there is at least one nuclei (N>1).
! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
! to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
      IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth 
        call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l))

        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30

        ! saturation at equilibrium
        ! rice should not be too small, otherwise seq value is not valid
        seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
      
        ! get resistance growth
        call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)

        res_out(ig,l) = res

        ! implicit scheme of mass growth
        ! cste here must be computed at each step
        cste = 4*pi*rho_ice*microtimestep

        dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)

        ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
        ! dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice) - pdq(ig,l,igcm_h2o_ice)*microtimestep) ! DIFF: take trend into account
        ! dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
        dMice = min(dMice,zq(ig,l,igcm_h2o_vap) + pdq(ig,l,igcm_h2o_vap)*microtimestep)

        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
        zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice

        countcells = countcells + 1 
      
        ! latent heat release
        lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
      
        ! DIFF: trend is enforce in a range, stabilize the scheme ?
        if (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
          subpdtcloud(ig,l)=5./microtimestep
        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
        if (subpdtcloud(ig,l)*microtimestep.lt.-5.0) then
            subpdtcloud(ig,l)=-5./microtimestep
        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then

        ! Special case of the isotope of water HDO    
        if (hdo) then
          ! condensation
          if (dMice.gt.0.0) then
            ! do we use fractionation?
            if (hdofrac) then
              ! Calculation of the HDO vapor coefficient
              Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) * kbz * zt(ig,l) / ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) * sqrt(1.+mhdo/mco2) )
              ! Calculation of the fractionnation coefficient at equilibrium
              alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
              ! alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
              ! Calculation of the 'real' fractionnation coefficient
              alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
              ! alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
            else
              alpha_c(ig,l) = 1.d0
            endif
            if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
              dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
            else
              dMice_hdo=0.
            endif
          !! sublimation
          else
            if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
              dMice_hdo=dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
            else
              dMice_hdo=0.
            endif
          endif !if (dMice.gt.0.0)

          dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
          dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))

          zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
          zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo

        endif ! if (hdo)        

!=============================================================
! 5. Dust cores released, tendancies, latent heat, etc ...
!=============================================================

        ! If all the ice particles sublimate, all the condensation
        ! nuclei are released:
        if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
          ! Water
          zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
          zq(ig,l,igcm_h2o_ice) = 0.
          if (hdo) then
            zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
            zq(ig,l,igcm_hdo_ice) = 0.
          endif
          ! Dust particles
          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
          ! CCNs
          zq(ig,l,igcm_ccn_mass) = 0.
          zq(ig,l,igcm_ccn_number) = 0.
        endif
      
      ELSE
        ! Initialization of dMice when it's not computed 
        dMice=0 ! no condensation/sublimation to account for
        subpdtcloud(ig,l)=0 ! no condensation/sublimation to account for
      ENDIF !of if Nccn>1
      
      ! No more getting tendency : we increment tracers & temp directly 

      ! Increment tracers & temp for following microtimestep
      ! Dust :
      ! Special treatment of dust_mass / number for scavenging ?
      IF (scavenging) THEN
        zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+pdq(ig,l,igcm_dust_mass)*microtimestep
        zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+pdq(ig,l,igcm_dust_number)*microtimestep
      ELSE
        zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
        zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
      ENDIF !(scavenging) THEN
      zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)*microtimestep
      zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)*microtimestep

      ! Water :
      zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*microtimestep
      zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*microtimestep

      ! HDO (if computed) :
      IF (hdo) THEN
        zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*microtimestep
        zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*microtimestep
      ENDIF ! hdo

      ! Could also set subpdtcloud to 0 if not activice to make it simpler or change name of the flag
      IF (.not.activice) subpdtcloud(ig,l)=0.

      ! Temperature :
      zt(ig,l) = zt(ig,l)+(pdt(ig,l)+subpdtcloud(ig,l))*microtimestep

      ! Prevent negative tracers ! JN
      WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30

      IF (cloud_adapt_ts) then
        ! Estimation of how much is actually condensing/subliming
         dMicetot=dMicetot+abs(dMice) ! DIFF: accumulation of absolute dMice
        ! dMicetot=dMicetot+dMice
        ! write(*,*) "dMicetot= ", dMicetot
        ! write(*,*) "we are in (l,ig):", l,ig !JN
        IF (spenttime.ne.0) then
          zdq=(dMicetot/spenttime)!*(ptimestep-spenttime)
        ELSE
          !Initialization for spenttime=0
          zdq=zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
        ENDIF
        ! Threshold with powerlaw (sanity check)
        ! zdq=min(abs(zdq),
        ! & abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
        zdq=abs(zdq)
        call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
      ENDIF! (cloud_adapt_ts) then
      ! Increment time spent in here
      spenttime=spenttime+microtimestep
      count_micro(ig,l)=count_micro(ig,l)+1
    ENDDO ! while (.not. ending_ts)
  ENDDO ! of ig loop
ENDDO ! of nlayer loop

!------ Useful outputs to check how it went
call write_output("zpotcond_inst","zpotcond_inst microphysics","(kg/kg)",zpotcond_inst(:,:))
call write_output("zpotcond_full","zpotcond_full microphysics","(kg/kg)",zpotcond_full(:,:))
call write_output("zpotcond","zpotcond microphysics","(kg/kg)",zpotcond(:,:))
call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:))

!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
IF (test_flag) then
  ! error2d(:) = 0.
  DO l=1,nlay
    DO ig=1,ngrid
      ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
      satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
      satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l) 
    ENDDO
  ENDDO

  print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
ENDIF ! endif test_flag

END SUBROUTINE improvedclouds

!=============================================================
! The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
! It is an analytical solution to the ice radius growth equation, 
! with the approximation of a constant 'reduced' cunningham correction factor 
! (lambda in growthrate.F) taken at radius req instead of rice    
!=============================================================

! subroutine phi(rice,req,coeff1,coeff2,time)
!      
! implicit none
!      
! ! inputs
! real rice ! ice radius
! real req  ! ice radius at equilibirum
! real coeff1  ! coeff for the log
! real coeff2  ! coeff for the arctan
!
! ! output      
! real time
!      
! !local
! real var
!      
!  1.73205 is sqrt(3)
!      
! var = max(
! &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
!            
! time = 
! &   coeff1 * 
! &   log( var )
! & + coeff2 * 1.73205 *
! &   atan( (2*rice+req) / (1.73205*req) )
!      
! return
! end

!=======================================================================

SUBROUTINE adapt_imicro(ptimestep,potcond,zimicro)

! Adaptative timestep for water ice clouds.
! Works using a powerlaw to compute the minimal duration of subtimestep
! (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
! Then, we use the instantaneous vap (ice) gradient extrapolated to the
! rest of duration to increase subtimestep duration, for computing
! efficiency.

real,intent(in) :: ptimestep ! total duration of physics (sec)
real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
real :: alpha, beta ! Coefficients for power law
real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
integer,intent(out) :: zimicro ! number of ptimestep division

! Default ptimestep : defstep (7.5 mins)
defstep=88775.*5./960.
coef=ptimestep/defstep
! Conservative coefficients :
! alpha=1.81846504e+11
! beta=1.54550140e+00
alpha=1.88282793e+05 ! DIFF: new power law coefficients
beta=4.57764370e-01
zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,50.),7000.))
zimicro=2*zimicro ! DIFF: prefiction times two, allow to complete year

END SUBROUTINE adapt_imicro

END MODULE improvedclouds_mod