Index: trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90
===================================================================
--- trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90	(revision 1926)
+++ trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90	(revision 3083)
@@ -1,21 +1,21 @@
-! Copyright 2017 Université de Reims Champagne-Ardenne 
+! Copyright 2017 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 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
 ! 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 license and that you accept its terms.
@@ -44,17 +44,17 @@
   IMPLICIT NONE
 
-  CONTAINS 
-    
+  CONTAINS
+
   SUBROUTINE mmp_initialize(dt,p_prod,tx_prod,rc_prod,rplanet,g0, air_rad,air_mmol,clouds,cfgpath)
     !! Initialize global parameters of the model.
-    !! 
+    !!
     !! The function initializes all the global parameters of the model from direct input.
-    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their 
-    !! default values are suitable for production runs.  
+    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their
+    !! default values are suitable for production runs.
     !! @note
     !! If the subroutine fails to initialize parameters, the run is aborted.
     !!
     !! @warning
-    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it 
+    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
     !! initializes global variable that are not thread private.
     !!
@@ -87,7 +87,7 @@
 
     INTEGER                                           :: coag_choice
-    REAL(kind=mm_wp)                                  :: fiad_max, fiad_min,df,rm,rho_aer
-    LOGICAL                                           :: w_h_prod, w_h_sed, w_h_coag, w_c_sed, w_c_nucond, &
-                                                         no_fiadero, fwsed_m0, fwsed_m3
+    REAL(kind=mm_wp)                                  :: fiad_max,fiad_min,df,rm,rho_aer
+    LOGICAL                                           :: w_h_prod,w_h_sed,w_h_coag,w_c_sed,w_c_nucond, &
+                                                         no_fiadero,fwsed_m0,fwsed_m3
     TYPE(error)                                       :: err
     INTEGER                                           :: i
@@ -96,4 +96,6 @@
     CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
     REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE       :: tmp
+    REAL(kind=mm_wp)                                  :: m0as_min,rcs_min,m0af_min,rcf_min,m0n_min
+    LOGICAL                                           :: wdebug
 
     w_h_prod    = .true.
@@ -104,20 +106,26 @@
     fwsed_m0    = .true.
     fwsed_m3    = .false.
-    no_fiadero  = .false.
+    no_fiadero  = .true.
     fiad_min    = 0.1_mm_wp
     fiad_max    = 10._mm_wp
     coag_choice = 7
+    wdebug      = .false.
+    m0as_min    = 1e-10_mm_wp
+    rcs_min     = 1e-9_mm_wp
+    m0af_min    = 1e-10_mm_wp
+    rcf_min     = 1e-9_mm_wp
+    m0n_min     = 1e-10_mm_wp
 
     WRITE(*,'(a)') "##### MMP_GCM SPEAKING #####"
-    WRITE(*,'(a)') "I will initialize ze microphysics model in moments YAMMS"
+    WRITE(*,'(a)') "I will initialize the microphysics model in moments YAMMS"
     WRITE(*,'(a)') "On error I will simply abort the program. Stay near your computer !"
     WRITE(*,*)
     WRITE(*,'(a)') "Reading muphys configuration file ("//trim(cfgpath)//")..."
-    err = cfg_read_config(cparser,TRIM(cfgpath),.true.) 
+    err = cfg_read_config(cparser,TRIM(cfgpath),.true.)
     IF (err /= 0) THEN
       ! RETURN AN ERROR !!
       call abort_program(err)
     ENDIF
-   
+
     ! YAMMS internal parameters:
     err = mm_check_opt(cfg_get_value(cparser,"rm",rm),rm,50e-9_mm_wp,mm_log)
@@ -126,26 +134,21 @@
     ! the following parameters are primarily used to test and debug YAMMS.
     ! They are set in an optional configuration file and default to suitable values for production runs.
-    err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)    ,w_h_prod   ,.true.   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)  ,w_h_sed    ,.true.   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)   ,w_h_coag   ,.true.   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed),w_c_sed    ,clouds   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond) ,w_c_nucond ,clouds   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)            ,fwsed_m0   ,.true.   ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)            ,fwsed_m3   ,.false.  ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)       ,no_fiadero ,.false.  ,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)  ,fiad_min   ,0.1_mm_wp,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)  ,fiad_max   ,10._mm_wp,mm_log)
-    err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7,mm_log)
-
-    ! optic look-up table file path.
-    mmp_optic_file = ''
-    opt_file = ''
-    err = mm_check_opt(cfg_get_value(cparser,"optics/optic_file",opt_file),opt_file,'',mm_log)
-    IF (err /= 0) THEN
-      WRITE(*,'(a)') "Warning: I was unable to retrieve the path of the optic look-up table file:"
-      WRITE(*,'(a)') "  The GCM may abort if it uses YAMMS optical properties calculation module !"
-    ELSE
-      mmp_optic_file = TRIM(opt_file)
-    ENDIF
+    err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)          ,w_h_prod   ,.true.     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)        ,w_h_sed    ,.true.     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)         ,w_h_coag   ,.true.     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed)      ,w_c_sed    ,clouds     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond)       ,w_c_nucond ,clouds     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)                  ,fwsed_m0   ,.true.     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)                  ,fwsed_m3   ,.false.    ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)             ,no_fiadero ,.true.     ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)        ,fiad_min   ,0.1_mm_wp  ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)        ,fiad_max   ,10._mm_wp  ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7          ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"m0as_min",m0as_min)                 ,m0as_min   ,1e-10_mm_wp,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"rcs_min",rcs_min)                   ,rcs_min    ,1e-9_mm_wp ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"m0af_min",m0af_min)                 ,m0af_min   ,1e-10_mm_wp,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"rcf_min",rcf_min)                   ,rcf_min    ,rm         ,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"m0n_min",m0n_min)                   ,m0n_min    ,1e-10_mm_wp,mm_log)
+    err = mm_check_opt(cfg_get_value(cparser,"debug",wdebug)                      ,wdebug     ,.false.    ,mm_log)
 
     ! Retrieve clouds species configuration file
@@ -156,66 +159,5 @@
     ENDIF
 
-    ! YAMMS initialization.
-    err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
-                           air_rad,air_mmol,coag_choice,clouds,spcpath,  &
-                           w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
-                           w_c_sed,fwsed_m0,fwsed_m3, &
-                           no_fiadero,fiad_min,fiad_max)
-    IF (err /= 0) call abort_program(err)
-
-    ! Extra initialization (needed for YAMMS method interfaces)
-    err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
-    IF (err/=0) call abort_program(err)
-    err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
-    IF (err/=0) call abort_program(err)
-
-    ! initialize transfert probabilities look-up tables
-    IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
-      err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
-      IF (err /= 0) call abort_program(err)
-
-      IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
-        call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1)) 
-      ENDIF
-      IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
-        call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1)) 
-      ENDIF
-      IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
-        call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
-      ENDIF
-      IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
-        call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
-      ENDIF
-    ENDIF
-    ! initialize mean electric correction look-up tables
-    IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
-      err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
-      IF (err /= 0) call abort_program(err)
-
-      IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
-        call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
-      ELSE
-        mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x) 
-        mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
-        mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
-        mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
-      ENDIF
-      IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
-        call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
-      ELSE
-        mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x) 
-        mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
-        mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
-        mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
-      ENDIF
-      IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
-        call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
-      ELSE
-        mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x) 
-        mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
-        mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
-        mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
-      ENDIF
-    ENDIF
+    ! Setup alpha function: THEY ARE REQUIRED IN YAMMS global initialization !
     ! spherical mode inter-moments function parameters
     IF (.NOT.cfg_has_section(cparser,'alpha_s')) call abort_program(error("Cannot find [alpha_s] section",-1))
@@ -227,11 +169,67 @@
     IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))
 
-    ! get size-distribution laws parameters
-    IF (.NOT.cfg_has_section(cparser,'dndr_s')) call abort_program(error("Cannot find [dndr_s] section",-2))
-    err = read_nprm(cparser,'dndr_s',mmp_pns)
-    IF (err /= 0) call abort_program(error("dndr_s: "//TRIM(err%msg),-2))
-    IF (.NOT.cfg_has_section(cparser,'dndr_f')) call abort_program(error("Cannot find [dndr_f] section",-2))
-    err = read_nprm(cparser,'dndr_f',mmp_pnf)
-    IF (err /= 0) call abort_program(error("dndr_f: "//TRIM(err%msg),-2))
+    ! YAMMS initialization.
+    err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
+                           air_rad,air_mmol,coag_choice,clouds,spcpath,  &
+                           w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
+                           w_c_sed,fwsed_m0,fwsed_m3, &
+                           no_fiadero,fiad_min,fiad_max, &
+                           m0as_min,rcs_min,m0af_min,rcf_min,m0n_min,wdebug)
+    IF (err /= 0) call abort_program(err)
+
+    ! Extra initialization (needed for YAMMS method interfaces)
+    err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
+    IF (err/=0) call abort_program(err)
+    err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
+    IF (err/=0) call abort_program(err)
+
+    ! initialize transfert probabilities look-up tables
+    IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
+      err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
+      IF (err /= 0) call abort_program(err)
+
+      IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
+        call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1))
+      ENDIF
+      IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
+        call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1))
+      ENDIF
+      IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
+        call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
+      ENDIF
+      IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
+        call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
+      ENDIF
+    ENDIF
+    ! initialize mean electric correction look-up tables
+    IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
+      err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
+      IF (err /= 0) call abort_program(err)
+
+      IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
+        call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
+      ELSE
+        mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x)
+        mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
+        mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
+        mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
+      ENDIF
+      IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
+        call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
+      ELSE
+        mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x)
+        mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
+        mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
+        mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
+      ENDIF
+      IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
+        call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
+      ELSE
+        mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x)
+        mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
+        mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
+        mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
+      ENDIF
+    ENDIF
 
     ! btk coefficients
@@ -249,10 +247,16 @@
     WRITE(*,'(a,L2)')     "electric_charging    : ", mmp_w_qe
     call mm_dump_parameters()
-      
+
+    IF (clouds) THEN
+      DO i=1, size(mm_xESPS)
+        print*, TRIM(mm_xESPS(i)%name), " fmol2fmas = ", mm_xESPS(i)%fmol2fmas
+      ENDDO
+    ENDIF
+
   END SUBROUTINE mmp_initialize
 
   FUNCTION read_aprm(parser,sec,pp) RESULT(err)
-    !! Read and store [[mmp_gcm(module):aprm(type)]] parameters. 
-    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser 
+    !! Read and store [[mmp_gcm(module):aprm(type)]] parameters.
+    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser
     CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
     TYPE(aprm), INTENT(out)      :: pp     !! [[mmp_gcm(module):aprm(type)]] object that stores the parameters values.
@@ -266,20 +270,4 @@
   END FUNCTION read_aprm
 
-  FUNCTION read_nprm(parser,sec,pp) RESULT(err)
-    !! Read and store [[mmp_gcm(module):nprm(type)]] parameters.
-    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser 
-    CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
-    TYPE(nprm), INTENT(out)      :: pp     !! [[mmp_gcm(module):nprm(type)]] object that stores the parameters values.
-    TYPE(error) :: err                     !! Error status of the function.
-    err = cfg_get_value(parser,TRIM(sec)//'/rc',pp%rc) ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'/a0',pp%a0) ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'/c',pp%c)   ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'/a',pp%a)   ; IF (err /= 0) RETURN
-    err = cfg_get_value(parser,TRIM(sec)//'/b',pp%b)   ; IF (err /= 0) RETURN
-    IF (SIZE(pp%a) /= SIZE(pp%b)) &
-      err = error("Inconsistent number of coefficients (a and b must have the same size)",3)
-    RETURN
-  END FUNCTION read_nprm
-
 END MODULE MMP_GCM
 
