Index: trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90
===================================================================
--- trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90	(revision 2109)
+++ trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90	(revision 3083)
@@ -1,21 +1,21 @@
-! Copyright 2013-2015,2017 Université de Reims Champagne-Ardenne 
+! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
 ! Contributor: J. Burgalat (GSMA, URCA)
 ! email of the author : jeremie.burgalat@univ-reims.fr
-! 
+!
 ! This software is a computer program whose purpose is to compute
 ! microphysics processes using a two-moments scheme.
-! 
+!
 ! This library is governed by the CeCILL-B license under French law and
-! abiding by the rules of distribution of free software.  You can  use, 
+! abiding by the rules of distribution of free software.  You can  use,
 ! modify and/ or redistribute the software under the terms of the CeCILL-B
 ! license as circulated by CEA, CNRS and INRIA at the following URL
-! "http://www.cecill.info". 
-! 
+! "http://www.cecill.info".
+!
 ! As a counterpart to the access to the source code and  rights to copy,
 ! modify and redistribute granted by the license, users are provided only
 ! with a limited warranty  and the software's author,  the holder of the
 ! economic rights,  and the successive licensors  have only  limited
-! liability. 
-! 
+! liability.
+!
 ! In this respect, the user's attention is drawn to the risks associated
 ! with loading,  using,  modifying and/or developing or reproducing the
@@ -25,8 +25,8 @@
 ! professionals having in-depth computer knowledge. Users are therefore
 ! encouraged to load and test the software's suitability as regards their
-! requirements in conditions enabling the security of their systems and/or 
-! data to be ensured and,  more generally, to use and operate it in the 
-! same conditions as regards security. 
-! 
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and,  more generally, to use and operate it in the
+! same conditions as regards security.
+!
 ! The fact that you are presently reading this means that you have had
 ! knowledge of the CeCILL-B license and that you accept its terms.
@@ -35,5 +35,5 @@
 !! summary: Haze microphysics module.
 !! author: J. Burgalat
-!! date: 2013-2015,2017
+!! date: 2013-2015,2017,2022
 MODULE MM_HAZE
   !! Haze microphysics module.
@@ -46,19 +46,19 @@
   !!
   !! @note
-  !! The production function is specific to Titan, where aerosols are created above the detached 
-  !! haze layer. No other source is taken into account.  This process is controled by two parameters, 
-  !! the pressure level of production and the production rate. Then both M0 and M3 of the aerosols 
+  !! The production function is specific to Titan, where aerosols are created above the detached
+  !! haze layer. No other source is taken into account.  This process is controled by two parameters,
+  !! the pressure level of production and the production rate. Then both M0 and M3 of the aerosols
   !! distribution are updated in the production zone by addition of the production rate along a
   !! gaussian shape.
   !!
   !! @note
-  !! The interface methods always uses the global variables defined in [[mm_globals(module)]] when 
+  !! The interface methods always uses the global variables defined in [[mm_globals(module)]] when
   !! values (any kind, temperature, pressure, moments...) over the vertical grid are required.
   !!
   !! @warning
-  !! The tendencies returned by the method are always defined over the vertical grid from __TOP__ 
+  !! The tendencies returned by the method are always defined over the vertical grid from __TOP__
   !! to __GROUND__.
   !!
-  !! @todo 
+  !! @todo
   !! Modify tests on tendencies vectors to get sure that allocation is done:
   !! Currently, we assume the compiler handles automatic allocation of arrays.
@@ -72,7 +72,7 @@
 
   PUBLIC :: mm_haze_microphysics, mm_haze_coagulation, mm_haze_sedimentation, &
-            mm_haze_production
-
-  CONTAINS
+    mm_haze_production
+
+CONTAINS
 
   !============================================================================
@@ -83,15 +83,15 @@
     !! Get the evolution of moments tracers through haze microphysics processes.
     !!
-    !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies 
-    !! of moments tracers for coagulation, sedimentation and production processes for the 
+    !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies
+    !! of moments tracers for coagulation, sedimentation and production processes for the
     !! atmospheric column.
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s
-      !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s
-      !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f
-      !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f
-      !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as
@@ -107,5 +107,5 @@
     zdm3af(1:mm_nla) = 0._mm_wp
 
-    IF (mm_w_haze_coag) THEN 
+    IF (mm_w_haze_coag) THEN
       ! Calls coagulation
       call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f)
@@ -120,5 +120,5 @@
 
       ! Updates tendencies
-      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 
+      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as
       dm0a_f=dm0a_f+zdm0af ; dm3a_f=dm3a_f+zdm3af
     ENDIF
@@ -127,5 +127,5 @@
       call mm_haze_production(zdm0as,zdm3as)
       ! We only produce spherical aerosols
-      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as 
+      dm0a_s=dm0a_s+zdm0as ; dm3a_s=dm3a_s+zdm3as
     ENDIF
 
@@ -137,8 +137,8 @@
   ! COAGULATION PROCESS RELATED METHODS
   !============================================================================
- 
+
   SUBROUTINE mm_haze_coagulation(dM0s,dM3s,dM0f,dM3f)
     !! Get the evolution of the aerosols moments vertical column due to coagulation process.
-    !! 
+    !!
     !! This is main method of the coagulation process:
     !!
@@ -149,13 +149,13 @@
     !! 5. Finally computes the tendencies of the moments.
     !!
-    !! All arguments are assumed vectors of __N__ elements where __N__ is the total number of 
+    !! All arguments are assumed vectors of __N__ elements where __N__ is the total number of
     !! vertical __layers__.
     !!
     !! @note
-    !! The method uses directly the global variables related to the vertical atmospheric structure 
-    !! stored in [[mm_globals(module)]]. Consequently they must be updated before calling the subroutine. 
+    !! The method uses directly the global variables related to the vertical atmospheric structure
+    !! stored in [[mm_globals(module)]]. Consequently they must be updated before calling the subroutine.
     !!
     !! @bug
-    !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm), 
+    !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm),
     !! a floating point exception occured (i.e. a NaN) as we perform a division by zero
     !!
@@ -163,13 +163,13 @@
     !! Get rid of the fu\*\*\*\* STOP statement...
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s
-      !! Tendency of the 0th order moment of the spherical size-distribution over a time step (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the spherical size-distribution over a time step (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3s
-      !! Tendency of the 3rd order moment of the spherical size-distribution (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the spherical size-distribution (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0f
-      !! Tendency of the 0th order moment of the fractal size-distribution over a time step (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the fractal size-distribution over a time step (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3f
-      !! Tendency of the 3rd order moment of the fractal size-distribution over a time step (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the fractal size-distribution over a time step (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c_kco,c_kfm,c_slf,tmp, &
-                                                   kco,kfm,pco,pfm,mq
+      kco,kfm,pco,pfm,mq
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a_ss,a_sf,b_ss,b_ff,c_ss,c_sf
     INTEGER                                     :: i
@@ -180,12 +180,12 @@
     ! Alloctes local arrays
     ALLOCATE(kco(mm_nla),kfm(mm_nla),c_slf(mm_nla), &
-             c_kco(mm_nla),c_kfm(mm_nla),mq(mm_nla), &
-             pco(mm_nla),pfm(mm_nla))
+      c_kco(mm_nla),c_kfm(mm_nla),mq(mm_nla), &
+      pco(mm_nla),pfm(mm_nla))
     ALLOCATE(a_ss(mm_nla),a_sf(mm_nla), &
-             b_ss(mm_nla),b_ff(mm_nla), &
-             c_ss(mm_nla),c_sf(mm_nla))
-             
+      b_ss(mm_nla),b_ff(mm_nla), &
+      c_ss(mm_nla),c_sf(mm_nla))
+
     a_ss(:) = 0._mm_wp ; a_sf(:) = 0._mm_wp
-    b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp 
+    b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp
     c_ss(:) = 0._mm_wp ; c_sf(:) = 0._mm_wp
 
@@ -193,5 +193,5 @@
     c_kco(:) = mm_get_kco(mm_temp) ; c_kfm(:) = mm_get_kfm(mm_temp)
     ! get slf (slip-flow factor)
-    c_slf(:) = mm_akn * mm_lambda_g(mm_temp,mm_play) 
+    c_slf(:) = mm_akn * mm_lambda_g(mm_temp,mm_play)
 
     DO i=1,mm_nla
@@ -202,6 +202,6 @@
         pfm(i) = mm_ps2s(mm_rcs(i),0,1,mm_temp(i),mm_play(i))
         ! (A_SS_CO x A_SS_FM) / (A_SS_CO + A_SS_FM)
-        kco(i) = g0ssco(mm_rcs(i),c_slf(i),c_kco(i)) 
-        kfm(i) = g0ssfm(mm_rcs(i),c_kfm(i)) 
+        kco(i) = g0ssco(mm_rcs(i),c_slf(i),c_kco(i))
+        kfm(i) = g0ssfm(mm_rcs(i),c_kfm(i))
         IF (kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp) /=0) THEN
           a_ss(i) = (kco(i)*(pco(i)-2._mm_wp)*kfm(i)*(pfm(i)-2._mm_wp))/(kco(i)*(pco(i)-2._mm_wp)+kfm(i)*(pfm(i)-2._mm_wp))
@@ -267,5 +267,5 @@
       ENDIF
     ENDDO
-   
+
     DEALLOCATE(kco,kfm,c_kco,c_kfm,pco,pfm,c_slf)
 
@@ -302,5 +302,5 @@
     !! Get &gamma; pre-factor for the 0th order moment with SS interactions in the continuous flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -320,5 +320,5 @@
     !! Get &gamma; pre-factor for the 0th order moment with SF interactions in the continuous flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ or __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -329,11 +329,11 @@
     REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, e, rcff
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
-    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 
-    ! computes mm_alpha coefficients
-    a1=mm_alpha_s(1._mm_wp)   ; a2=mm_alpha_f(-e)   ; a3=mm_alpha_f(e) 
+    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra
+    ! computes mm_alpha coefficients
+    a1=mm_alpha_s(1._mm_wp)   ; a2=mm_alpha_f(-e)   ; a3=mm_alpha_f(e)
     a4=mm_alpha_s(-1._mm_wp) ; a5=mm_alpha_s(-2._mm_wp) ; a6=mm_alpha_f(-2._mm_wp*e)
     ! Computes gamma pre-factor
     res = c_kco*( 2._mm_wp + a1*a2*rcs/rcff + a4*a3*rcff/rcs + c_slf*( a4/rcs + &
-                         a2/rcff + a5*a3*rcff/rcs**2 + a1*a6*rcs/rcff**2))
+      a2/rcff + a5*a3*rcff/rcs**2 + a1*a6*rcs/rcff**2))
     RETURN
   END FUNCTION g0sfco
@@ -342,5 +342,5 @@
     !! Get &gamma; pre-factor for the 0th order moment with FF interactions in the continuous flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
@@ -361,5 +361,5 @@
     !! Get &gamma; pre-factor for the 3rd order moment with SS interactions in the continuous flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -370,10 +370,10 @@
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(2._mm_wp)  ; a3=mm_alpha_s(1._mm_wp) 
+    a1=mm_alpha_s(3._mm_wp) ; a2=mm_alpha_s(2._mm_wp)  ; a3=mm_alpha_s(1._mm_wp)
     a4=mm_alpha_s(4._mm_wp) ; a5=mm_alpha_s(-1._mm_wp) ; a6=mm_alpha_s(-2._mm_wp)
 
     ! Computes gamma pre-factor
-    res = (2._mm_wp*a1 + a2*a3 + a4*a5 + c_slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))* & 
-          c_kco/(a1**2*rcs**3)
+    res = (2._mm_wp*a1 + a2*a3 + a4*a5 + c_slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2))* &
+      c_kco/(a1**2*rcs**3)
     RETURN
   END FUNCTION g3ssco
@@ -382,5 +382,5 @@
     !! Get &gamma; pre-factor for the 3rd order moment with SF interactions in the continuous flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ or __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -392,12 +392,12 @@
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra 
-    a1=mm_alpha_s(3._mm_wp)    ; a2=mm_alpha_s(4._mm_wp) ; a3=mm_alpha_f(-e) 
+    e=3._mm_wp/mm_df ; rcff = rcf**e*mm_rb2ra
+    a1=mm_alpha_s(3._mm_wp)    ; a2=mm_alpha_s(4._mm_wp) ; a3=mm_alpha_f(-e)
     a4=mm_alpha_s(2._mm_wp)    ; a5=mm_alpha_f(e)    ; a6=mm_alpha_s(1._mm_wp)
     a7=mm_alpha_f(-2._mm_wp*e) ; a8=mm_alpha_f(3._mm_wp)
     ! Computes gamma pre-factor
     res = (2._mm_wp*a1*rcs**3 + a2*rcs**4*a3/rcff + a4*rcs**2*a5*rcff + &
-                c_slf *( a4*rcs**2 + a1*rcs**3*a3/rcff + a6*rcs*a5*rcff + &
-                a2*rcs**4*a7/rcff**2))* c_kco/(a1*a8*(rcs*rcf)**3)
+      c_slf *( a4*rcs**2 + a1*rcs**3*a3/rcff + a6*rcs*a5*rcff + &
+      a2*rcs**4*a7/rcff**2))* c_kco/(a1*a8*(rcs*rcf)**3)
     RETURN
   END FUNCTION g3sfco
@@ -406,5 +406,5 @@
     !! Get &gamma; pre-factor for the 0th order moment with SS interactions in the Free Molecular flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -422,7 +422,7 @@
 
   ELEMENTAL FUNCTION g0sffm(rcs, rcf, c_kfm) RESULT(res)
-    !> Get &gamma; pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime. 
-    !!
-    !! @note 
+    !> Get &gamma; pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime.
+    !!
+    !! @note
     !! If __rcs__ or __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -435,38 +435,38 @@
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    e1 = 3._mm_wp/mm_df 
-    e2 = 6._mm_wp/mm_df 
+    e1 = 3._mm_wp/mm_df
+    e2 = 6._mm_wp/mm_df
     e3 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
     e4 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
 
-    rcff1 = mm_rb2ra * rcf**e1 
+    rcff1 = mm_rb2ra * rcf**e1
     rcff2 = rcff1**2
-    rcff3 = mm_rb2ra * rcf**e3 
+    rcff3 = mm_rb2ra * rcf**e3
     rcff4 = mm_rb2ra**2 * rcf**e4
 
-    a1=mm_alpha_s(0.5_mm_wp)  ; a2=mm_alpha_s(-0.5_mm_wp) ; a3=mm_alpha_f(e1) 
+    a1=mm_alpha_s(0.5_mm_wp)  ; a2=mm_alpha_s(-0.5_mm_wp) ; a3=mm_alpha_f(e1)
     a4=mm_alpha_s(-1.5_mm_wp) ; a5=mm_alpha_f(e2)      ; a6=mm_alpha_s(2._mm_wp)
-    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(1._mm_wp)   ; a9=mm_alpha_f(e3) 
+    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(1._mm_wp)   ; a9=mm_alpha_f(e3)
     a10=mm_alpha_f(e4)
 
     ! Computes gamma pre-factor
     res = (a1*rcs**0.5_mm_wp + 2._mm_wp*rcff1*a2*a3/rcs**0.5_mm_wp + a4*a5*rcff2/rcs**1.5_mm_wp + &
-           a6*a7*rcs**2/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs*rcff3 + a10*rcff4 &
-          )*mm_get_btk(4,0)*c_kfm
+      a6*a7*rcs**2/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs*rcff3 + a10*rcff4 &
+      )*mm_get_btk(4,0)*c_kfm
     RETURN
   END FUNCTION g0sffm
 
   ELEMENTAL FUNCTION g0fffm(rcf, c_kfm) RESULT(res)
-    !! Get &gamma; pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime. 
-    !!
-    !! @note 
+    !! Get &gamma; pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime.
+    !!
+    !! @note
     !! If __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcf   !! Characteristic radius of the fractal size-distribution.
     REAL(kind=mm_wp), INTENT(in) :: c_kfm !! Thermodynamic free molecular flow regime pre-factor.
-    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor. 
+    REAL(kind=mm_wp) :: res               !! &gamma; coagulation kernel pre-factor.
     REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, e1, e2, e3, rcff
     res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    e1=3._mm_wp/mm_df ; e2=(6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 
+    e1=3._mm_wp/mm_df ; e2=(6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
     e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
     rcff=mm_rb2ra**2*rcf**e3
@@ -489,11 +489,11 @@
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(1._mm_wp)  ; a3=mm_alpha_s(2.5_mm_wp) 
+    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(1._mm_wp)  ; a3=mm_alpha_s(2.5_mm_wp)
     a4=mm_alpha_s(2._mm_wp)   ; a5=mm_alpha_s(1.5_mm_wp) ; a6=mm_alpha_s(5._mm_wp)
-    a7=mm_alpha_s(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)  ; a9=mm_alpha_s(-0.5_mm_wp) 
+    a7=mm_alpha_s(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)  ; a9=mm_alpha_s(-0.5_mm_wp)
     a10=mm_alpha_s(3._mm_wp) ; a11=mm_alpha_s(0.5_mm_wp)
     ! Computes gamma pre-factor
     res = (a1 + 2._mm_wp*a2*a3 + a4*a5 + a6*a7 + 2._mm_wp*a8*a9 + a10*a11) &
-          *mm_get_btk(1,3)*c_kfm/(a10**2*rcs**2.5_mm_wp)
+      *mm_get_btk(1,3)*c_kfm/(a10**2*rcs**2.5_mm_wp)
     RETURN
   END FUNCTION g3ssfm
@@ -502,5 +502,5 @@
     !! Get &gamma; pre-factor for the 3rd order moment with SF interactions in the Free Molecular flow regime.
     !!
-    !! @note 
+    !! @note
     !! If __rcs__ or __rcf__ is 0, the function returns 0.
     REAL(kind=mm_wp), INTENT(in) :: rcs   !! Characteristic radius of the spherical size-distribution.
@@ -512,16 +512,16 @@
     res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN
     ! computes mm_alpha coefficients
-    e1=3._mm_wp/mm_df 
-    e2=(6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) 
+    e1=3._mm_wp/mm_df
+    e2=(6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
     e3=(12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df)
     rcff1=mm_rb2ra*rcf**e1 ; rcff2=mm_rb2ra*rcf**e2 ; rcff3=mm_rb2ra**2*rcf**e3
-    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(2.5_mm_wp)   ; a3=mm_alpha_f(e1) 
-    a4=mm_alpha_s(1.5_mm_wp)  ; a5=mm_alpha_f(2._mm_wp*e1) ; a6=mm_alpha_s(5._mm_wp) 
-    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)    ; a9=mm_alpha_f(e2) 
+    a1=mm_alpha_s(3.5_mm_wp)  ; a2=mm_alpha_s(2.5_mm_wp)   ; a3=mm_alpha_f(e1)
+    a4=mm_alpha_s(1.5_mm_wp)  ; a5=mm_alpha_f(2._mm_wp*e1) ; a6=mm_alpha_s(5._mm_wp)
+    a7=mm_alpha_f(-1.5_mm_wp) ; a8=mm_alpha_s(4._mm_wp)    ; a9=mm_alpha_f(e2)
     a10=mm_alpha_s(3._mm_wp)  ; a11=mm_alpha_f(e3)         ; a12=mm_alpha_f(3._mm_wp)
     ! Computes gamma pre-factor
     res = (a1*rcs**3.5_mm_wp + 2._mm_wp*a2*a3*rcs**2.5_mm_wp*rcff1 + a4*a5*rcs**1.5_mm_wp*rcff1**2 + &
-          a6*a7*rcs**5/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs**4*rcff2 + & 
-          a10*a11*rcs**3*rcff3)*mm_get_btk(4,3)*c_kfm/(a10*a12*(rcs*rcf)**3)
+      a6*a7*rcs**5/rcf**1.5_mm_wp + 2._mm_wp*a8*a9*rcs**4*rcff2 + &
+      a10*a11*rcs**3*rcff3)*mm_get_btk(4,3)*c_kfm/(a10*a12*(rcs*rcf)**3)
     RETURN
   END FUNCTION g3sffm
@@ -530,21 +530,20 @@
   ! SEDIMENTATION PROCESS RELATED METHODS
   !============================================================================
-  
+
   SUBROUTINE mm_haze_sedimentation(dm0s,dm3s,dm0f,dm3f)
     !! Interface to sedimentation algorithm.
     !!
     !! The subroutine computes the evolution of each moment of the aerosols tracers
-    !! through sedimentation process and returns their tendencies for a timestep. 
+    !! through sedimentation process and returns their tendencies for a timestep.
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s
-      !! Tendency of the 0th order moment of the spherical mode (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the spherical mode (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s
-      !! Tendency of the 3rd order moment of the spherical mode (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the spherical mode (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0f
-      !! Tendency of the 0th order moment of the fractal mode (\(m^{-3}\)).
+    !! Tendency of the 0th order moment of the fractal mode (\(m^{-3}\)).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3f
-      !! Tendency of the 3rd order moment of the fractal mode (\(m^{3}.m^{-3}\)).
+    !! Tendency of the 3rd order moment of the fractal mode (\(m^{3}.m^{-3}\)).
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: ft,fdcor,wth
-    REAL(kind=mm_wp)                            :: m,n,p 
-    REAL(kind=mm_wp), PARAMETER                 :: fac = 4._mm_wp/3._mm_wp * mm_pi 
+    REAL(kind=mm_wp), PARAMETER                 :: fac = 4._mm_wp/3._mm_wp * mm_pi
 
     ALLOCATE(ft(mm_nle),wth(mm_nle),fdcor(mm_nle))
@@ -624,5 +623,5 @@
     !! Compute the tendency of the moment through sedimentation process.
     !!
-    !! 
+    !!
     !! The method computes the time evolution of the \(k^{th}\) order moment through sedimentation:
     !!
@@ -630,6 +629,6 @@
     !!
     !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method).
-    !! 
-    !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells 
+    !!
+    !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells
     !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm
     !! originally implemented in the LMDZ-Titan 2D GCM.
@@ -638,5 +637,5 @@
     REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: ft  !! Downward sedimentation flux  (effective velocity of the moment).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of \(k^{th}\) order moment (in \(m^{k}.m^{-3}\)).
-    INTEGER                                 :: i 
+    INTEGER                                 :: i
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko
     ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla))
@@ -648,5 +647,5 @@
       bs(1:mm_nla) = -(ft(2:mm_nle)+mm_dzlay(1:mm_nla)/dt)
       cs(1:mm_nla) = -mm_dzlay(1:mm_nla)/dt*mk(1:mm_nla)
-      ! (Tri)diagonal matrix inversion 
+      ! (Tri)diagonal matrix inversion
       mko(1) = cs(1)/bs(1)
       DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO
@@ -660,6 +659,6 @@
       ! interior
       mko(2:mm_nla-1)=(bs(2:mm_nla-1)*mk(1:mm_nla-2) + &
-                       cs(2:mm_nla-1)*mk(2:mm_nla-1)   &
-                      )/as(2:mm_nla-1)
+        cs(2:mm_nla-1)*mk(2:mm_nla-1)   &
+        )/as(2:mm_nla-1)
     ENDIF
     dmk = mko - mk
@@ -670,12 +669,12 @@
   SUBROUTINE get_weff(mk,k,df,rc,dt,afun,wth,corf)
     !! Get the effective settling velocity for aerosols moments.
-    !! 
-    !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol 
-    !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following 
+    !!
+    !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol
+    !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following
     !! equation:
-    !! 
-    !! $$ 
+    !!
+    !! $$
     !! \begin{eqnarray*}
-    !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr 
+    !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr
     !!                     == M_{k}  \times v^{eff}_{M_{k}} \\
     !!     v^{eff}_{M_{k} &=& \dfrac{2 \rho g r_{c}^{\dfrac{3D_{f}-3}{D_{f}}}}
@@ -686,24 +685,24 @@
     !! $$
     !!
-    !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm 
+    !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm
     !! as defined in \cite{toon1988b}.
     !!
     !! @warning
-    !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will 
+    !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will
     !! be computed.
     REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: mk
-      !! Moment of order __k__ (\(m^{k}.m^{-3}\)) at each layer.
+    !! Moment of order __k__ (\(m^{k}.m^{-3}\)) at each layer.
     REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla)            :: rc
-      !! Characteristic radius associated to the moment at each layer.
+    !! Characteristic radius associated to the moment at each layer.
     REAL(kind=mm_wp), INTENT(in)                               :: k
-      !! The order of the moment.
+    !! The order of the moment.
     REAL(kind=mm_wp), INTENT(in)                               :: df
-      !! Fractal dimension of the aersols.
+    !! Fractal dimension of the aersols.
     REAL(kind=mm_wp), INTENT(in)                               :: dt
-      !! Time step (s).
+    !! Time step (s).
     REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle)           :: wth
-      !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)).
-    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf 
-      !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__.
+    !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)).
+    REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle), OPTIONAL :: corf
+    !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__.
     INTERFACE
       FUNCTION afun(order)
@@ -713,5 +712,5 @@
         REAL(kind=mm_wp) :: afun              !! Alpha value.
       END FUNCTION afun
-    END INTERFACE 
+    END INTERFACE
     INTEGER                             :: i
     REAL(kind=mm_wp)                    :: af1,af2,ar1,ar2
@@ -720,22 +719,27 @@
     REAL(kind=mm_wp), DIMENSION(mm_nle) :: zcorf
     ! ------------------
-    
+
     wth(:) = 0._mm_wp ; zcorf(:) = 1._mm_wp
 
     ar1 = (3._mm_wp*df -3._mm_wp)/df    ; ar2 = (3._mm_wp*df -6._mm_wp)/df
-    af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df 
+    af1 = (df*(k+3._mm_wp)-3._mm_wp)/df ; af2 = (df*(k+3._mm_wp)-6._mm_wp)/df
     rb2ra = mm_rm**((df-3._mm_wp)/df)
     DO i=2,mm_nla
-      IF (rc(i-1) <= 0._mm_wp) CYCLE 
+      IF (rc(i-1) <= 0._mm_wp) CYCLE
       dzb = (mm_dzlay(i)+mm_dzlay(i-1))/2._mm_wp
       csto = 2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))/(9._mm_wp*mm_eta_g(mm_btemp(i)))
       cslf = mm_akn * mm_lambda_g(mm_btemp(i),mm_plev(i))
-      wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2)) 
+      wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2))
+
+      ! >>> [TEMPO : BBT]
+      !wth(i) = wth(i) * (2574e3 / (2574e3+mm_zlev(i)))**4
+      ! <<< [TEMPO : BBT]
+      
       ! now correct velocity to reduce numerical diffusion
       IF (.NOT.mm_no_fiadero_w) THEN
         IF (mk(i) <= 0._mm_wp) THEN
-          ratio = mm_fiadero_max 
+          ratio = mm_fiadero_max
         ELSE
-          ratio = MAX(MIN(mk(i-1)/mk(i),mm_fiadero_max),mm_fiadero_min) 
+          ratio = MAX(MIN(mk(i-1)/mk(i),mm_fiadero_max),mm_fiadero_min)
         ENDIF
         ! apply correction
@@ -760,10 +764,10 @@
   SUBROUTINE mm_haze_production(dm0s,dm3s)
     !! Compute the production of aerosols moments.
-    !! 
-    !! The method computes the tendencies of M0 and M3 for the spherical mode through production process. 
-    !! Production values are distributed along a normal law in altitude, centered  at 
+    !!
+    !! The method computes the tendencies of M0 and M3 for the spherical mode through production process.
+    !! Production values are distributed along a normal law in altitude, centered  at
     !! [[mm_globals(module):mm_p_prod(variable)]] pressure level with a fixed sigma of 20km.
     !!
-    !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation a spherical 
+    !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation a spherical
     !! characteristic radius set to [[mm_globals(module):mm_rc_prod(variable)]].
     !!
@@ -773,5 +777,5 @@
     REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (\(m^{3}.m^{-3}\)).
     INTEGER                     :: i
-    REAL(kind=mm_wp)            :: zprod,cprod,timefact 
+    REAL(kind=mm_wp)            :: zprod,cprod,timefact
     REAL(kind=mm_wp), PARAMETER :: sigz  = 20e3_mm_wp, &
                                    fnorm = 1._mm_wp/(dsqrt(2._mm_wp*mm_pi)*sigz), &
@@ -793,5 +797,5 @@
     dm3s(:)= mm_tx_prod *0.75_mm_wp/mm_pi *mm_dt / mm_rhoaer / 2._mm_wp / mm_dzlev(1:mm_nla) * &
              (erf((mm_zlev(1:mm_nla)-zprod)/znorm) - &
-             erf((mm_zlev(2:mm_nla+1)-zprod)/znorm)) 
+             erf((mm_zlev(2:mm_nla+1)-zprod)/znorm))
     dm0s(:) = dm3s(:)/(mm_rc_prod**3*mm_alpha_s(3._mm_wp))
 
@@ -803,5 +807,4 @@
     ENDIF
 
-
   END SUBROUTINE mm_haze_production
 
