Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8522)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8531)
@@ -28,10 +28,10 @@
       &                 fr1_i0, fr2_i0
    USE ice1D          ! thermodynamic sea-ice variables
-   USE icethd_dif     ! vertical diffusion
+   USE icethd_zdf     ! vertical diffusion
    USE icethd_dh      ! ice-snow growth and melt
    USE icethd_da      ! lateral melting
    USE icethd_sal     ! ice salinity
    USE icethd_ent     ! ice enthalpy redistribution
-   USE icethd_lac     ! lateral accretion
+   USE icethd_do      ! lateral accretion
    USE iceitd         ! remapping thickness distribution
    USE icetab         ! 1D <==> 2D transformation
@@ -51,4 +51,10 @@
    PUBLIC   ice_thd_init    ! called by ice_init
 
+   !!** namelist (namthd) **
+   LOGICAL ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F)
+   LOGICAL ::   ln_icedA         ! activate lateral melting param. (T) or not (F)
+   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
+   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -70,9 +76,11 @@
       !!               at the ice base, snow acc.,heat budget of the leads)
       !!             - selection of the icy points and put them in an array
-      !!             - call ice_thd_dif  for vertical heat diffusion
+      !!             - call ice_thd_zdf  for vertical heat diffusion
       !!             - call ice_thd_dh   for vertical ice growth and melt
       !!             - call ice_thd_ent  for enthalpy remapping
       !!             - call ice_thd_sal  for ice desalination
       !!             - call ice_thd_temp to  retrieve temperature from ice enthalpy
+      !!             - call ice_thd_lam  for extra lateral ice melt if monocat
+      !!             - call ice_thd_da   for lateral ice melt
       !!             - back to the geographic grid
       !!---------------------------------------------------------------------
@@ -223,10 +231,10 @@
             !
             IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
-                              CALL ice_thd_dif                             ! Ice/Snow Temperature profile
+                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
                               CALL ice_thd_ent( e_i_1d(1:nidx,:) )         ! Ice enthalpy remapping
             ENDIF
             !
-                              CALL ice_thd_sal                      ! --- Ice salinity --- !    
+                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
             !
                               CALL ice_thd_temp                     ! --- temperature update --- !
@@ -261,5 +269,5 @@
       IF( jpl > 1 )        CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- !
       !
-      IF( ln_icedO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- !
+      IF( ln_icedO )       CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
       !
       ! controls
@@ -523,81 +531,46 @@
       !!                   ***  ROUTINE ice_thd_init *** 
       !!                 
-      !! ** Purpose :   Physical constants and parameters linked to the ice 
-      !!              thermodynamics
+      !! ** Purpose :   Physical constants and parameters associated with
+      !!                ice thermodynamics
       !!
-      !! ** Method  :   Read the namice_thd namelist and check the ice-thermo
-      !!              parameter values called at the first timestep (nit000)
+      !! ** Method  :   Read the namthd namelist and check the parameters
+      !!                called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namicether
+      !! ** input   :   Namelist namthd
       !!-------------------------------------------------------------------
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_thd/ rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   &
-         &                 ln_icedH, rn_blow_s,                                                    &
-         &                 ln_icedA, rn_beta, rn_dmin,                                             &
-         &                 ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz,          &
-         &                 nn_iceflx
+      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )              ! Namelist namice_thd in reference namelist : Ice thermodynamics
-      READ  ( numnam_ice_ref, namice_thd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namice_thd in configuration namelist : Ice thermodynamics
-      READ  ( numnam_ice_cfg, namice_thd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_thd )
+      REWIND( numnam_ice_ref )              ! Namelist namthd in reference namelist : Ice thermodynamics
+      READ  ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd in configuration namelist : Ice thermodynamics
+      READ  ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd )
       !
       !
       IF(lwp) THEN                          ! control print
-         WRITE(numout,*) 'ice_thd_init : Ice Thermodynamics'
-         WRITE(numout,*) '~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_thd'
-         WRITE(numout,*) '   -- icethd_dif --'
-         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
-         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
-         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
-         WRITE(numout,*) '      thermal conductivity of the snow                        rn_cnd_s     = ', rn_cnd_s
-         WRITE(numout,*) '   -- icethd_dh --'
-         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH     = ', ln_icedH
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall          rn_blow_s    = ', rn_blow_s
-         WRITE(numout,*) '   -- icethd_da --'
-         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA     = ', ln_icedA
-         WRITE(numout,*) '      Coef. beta for lateral melting param.                   rn_beta      = ', rn_beta
-         WRITE(numout,*) '      Minimum floe diameter for lateral melting param.        rn_dmin      = ', rn_dmin
-         WRITE(numout,*) '   -- icethd_lac --'
-         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO     = ', ln_icedO
-         WRITE(numout,*) '      ice thickness for lateral accretion                     rn_hinew     = ', rn_hinew
-         WRITE(numout,*) '      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil
-         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom   rn_maxfraz   = ', rn_maxfraz
-         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil rn_vfraz     = ', rn_vfraz
-         WRITE(numout,*) '      Squeezing coefficient for collection of frazil          rn_Cfraz     = ', rn_Cfraz
-         WRITE(numout,*) '   -- icestp --'
-         WRITE(numout,*) '      Multicategory heat flux formulation                     nn_iceflx    = ', nn_iceflx
+         WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namthd:'
+         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH  = ', ln_icedH
+         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA  = ', ln_icedA
+         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO  = ', ln_icedO
+         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS  = ', ln_icedS
+     ENDIF
+      !
+                       CALL ice_thd_zdf_init   ! set ice heat diffusion parameters
+      IF( ln_icedA )   CALL ice_thd_da_init    ! set ice lateral melting parameters
+      IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters
+                       CALL ice_thd_sal_init   ! set ice salinity parameters
+      !
+      IF( ln_icedS .AND. nn_icesal == 1 ) THEN
+         ln_icedS = .FALSE.
+         CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)')
       ENDIF
-      !
-      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
-         CALL ctl_stop( 'ice_thd_init: choose one and only one formulation for thermal conductivity (ln_cndi_U64 or ln_cndi_P07)' )
-      ENDIF
-      !
-      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_init : rn_hinew should be >= rn_himin' )
-      !
-      IF(lwp) WRITE(numout,*)
-      SELECT CASE( nn_iceflx )         ! LIM3 Multi-category heat flux formulation
-      CASE( -1  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use per-category fluxes (nn_iceflx = -1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
-      CASE(  0  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_iceflx = 0) '
-      CASE(  1  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_iceflx = 1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
-      CASE(  2  )
-         IF(lwp) WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_iceflx = 2) '
-         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in forced mode cannot be 2' )
-      CASE DEFAULT
-         CALL ctl_stop( 'ice_thd_init: LIM3 option, nn_iceflx, should be between -1 and 2' )
-      END SELECT
       !
    END SUBROUTINE ice_thd_init
