Changeset 8752

Timestamp:

2017-11-20T13:54:32+01:00 (5 years ago)

Author:

dancopsey

Message:

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

Location:

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM

Files:

• CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg (modified) (5 diffs)
• CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg (modified) (1 diff)
• CONFIG/SHARED/field_def_nemo-lim.xml (modified) (3 diffs)
• CONFIG/SHARED/namelist_ice_lim3_ref (modified) (9 diffs)
• CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_ice_cfg (modified) (1 diff)
• NEMO/LIM_SRC_3/ice.F90 (modified) (3 diffs)
• NEMO/LIM_SRC_3/ice1D.F90 (deleted)
• NEMO/LIM_SRC_3/ice1d.F90 (copied) (copied from branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90)
• NEMO/LIM_SRC_3/icealb.F90 (modified) (7 diffs)
• NEMO/LIM_SRC_3/icectl.F90 (modified) (3 diffs)
• NEMO/LIM_SRC_3/icedyn.F90 (modified) (1 diff)
• NEMO/LIM_SRC_3/icedyn_adv_pra.F90 (modified) (8 diffs)
• NEMO/LIM_SRC_3/icedyn_adv_umx.F90 (modified) (1 diff)
• NEMO/LIM_SRC_3/icedyn_rdgrft.F90 (modified) (8 diffs)
• NEMO/LIM_SRC_3/icedyn_rhg.F90 (modified) (3 diffs)
• NEMO/LIM_SRC_3/icedyn_rhg_evp.F90 (modified) (18 diffs)
• NEMO/LIM_SRC_3/iceforcing.F90 (modified) (5 diffs)
• NEMO/LIM_SRC_3/iceistate.F90 (modified) (3 diffs)
• NEMO/LIM_SRC_3/iceitd.F90 (modified) (7 diffs)
• NEMO/LIM_SRC_3/icerst.F90 (modified) (4 diffs)
• NEMO/LIM_SRC_3/icestp.F90 (modified) (7 diffs)
• NEMO/LIM_SRC_3/icethd.F90 (modified) (12 diffs)
• NEMO/LIM_SRC_3/icethd_dh.F90 (modified) (7 diffs)
• NEMO/LIM_SRC_3/icethd_pnd.F90 (copied) (copied from branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90)
• NEMO/LIM_SRC_3/icethd_zdf.F90 (modified) (20 diffs)
• NEMO/LIM_SRC_3/iceupdate.F90 (modified) (3 diffs)
• NEMO/LIM_SRC_3/icevar.F90 (modified) (8 diffs)
• NEMO/LIM_SRC_3/icewri.F90 (modified) (5 diffs)
• NEMO/LIM_SRC_3/limmp.F90 (deleted)
• NEMO/OPA_SRC/BDY/bdyice.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/albedooce.F90 (deleted)
• NEMO/OPA_SRC/SBC/ocealb.F90 (copied) (copied from branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/ocealb.F90)
• NEMO/OPA_SRC/SBC/sbc_ice.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/SBC/sbcblk.F90 (modified) (11 diffs)
• NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (11 diffs)

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg

@@ -13 +13 @@
 !!             11 - Ice growth in open water           (namthd_do)
 !!             12 - Ice salinity                       (namthd_sal)
-!!             13 - Ice melt ponds                     (nammp)
+!!             13 - Ice melt ponds                     (namthd_pnd)
 !!             14 - Ice initialization                 (namini)
 !!             15 - Ice/snow albedos                   (namalb)
@@ -20 +20 @@
 !
 !------------------------------------------------------------------------------
-&nampar     !   Generic parameters
+&nampar         !   Generic parameters
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namitd     !   Ice discretization
+&namitd         !   Ice discretization
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namdyn     !   Ice dynamics
+&namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
 /
@@ -48 +48 @@
 /
 !------------------------------------------------------------------------------
-&namthd     !   Ice thermodynamics
+&namthd         !   Ice thermodynamics
 !------------------------------------------------------------------------------
 /
@@ -56 +56 @@
 /
 !------------------------------------------------------------------------------
-&namthd_da     !   Ice lateral melting
+&namthd_da      !   Ice lateral melting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namthd_do     !   Ice growth in open water
+&namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
 /
@@ -68 +68 @@
 /
 !------------------------------------------------------------------------------
-&nammp      !   Melt ponds
+&namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namini     !   Ice initialization
+&namini         !   Ice initialization
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namalb     !   albedo parameters
+&namalb         !   albedo parameters
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namdia     !   Diagnostics
+&namdia         !   Diagnostics
 !------------------------------------------------------------------------------
 /

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg

@@ -1 +1 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !! ESIM configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref
-!!              1 - Generic parameters                 (namice_run)
-!!              2 - Ice thickness discretization       (namice_itd)
-!!              3 - Ice dynamics                       (namice_dyn)
-!!              4 - Ice ridging/rafting                (namice_rdgrft)
-!!              5 - Ice rheology                       (namice_rhg)
-!!              6 - Ice advection                      (namice_adv)
-!!              7 - Ice thermodynamics                 (namice_thd)
-!!              8 - Ice salinity                       (namice_sal)
-!!              9 - Ice melt ponds                     (namice_mp)
-!!             10 - Ice initialization                 (namice_ini)
-!!             11 - Ice/snow albedos                   (namice_alb)
-!!             12 - Ice diagnostics                    (namice_dia)
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface forcing                (namforcing)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (namthd_pnd)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !
 !------------------------------------------------------------------------------
-&namice_run     !   Generic parameters
+&nampar         !   Generic parameters
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_itd     !   Ice discretization
+&namitd         !   Ice discretization
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_dyn     !   Ice dynamics
+&namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_rdgrft  !   Ice ridging/rafting
+&namdyn_rdgrft  !   Ice ridging/rafting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_rhg     !   Ice rheology
+&namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_adv     !   Ice advection
+&namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_thd     !   Ice thermodynamics
+&namforcing     !   Ice surface forcing
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_sal     !   Ice salinity
+&namthd         !   Ice thermodynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicemp      !   Melt ponds
+&namthd_zdf     !   Ice heat diffusion
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_ini     !   Ice initialization
+&namthd_da      !   Ice lateral melting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_alb     !   albedo parameters
+&namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_dia     !   Diagnostics
+&namthd_sal     !   Ice salinity
 !------------------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+/

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml

@@ -132 +132 @@
          <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            /> 
 
-         <field id="iceconc_cat"  long_name="Sea-ice area fractions in thickness categories"               unit=""   grid_ref="grid_T_3D_ncatice" />
+         <field id="iceconc_cat"  long_name="Sea-ice concentration in thickness categories"                unit=""   grid_ref="grid_T_3D_ncatice" />
          <field id="icethic_cat"  long_name="Sea-ice thickness in thickness categories"                    unit="m"  grid_ref="grid_T_3D_ncatice" />
          <field id="snowthic_cat" long_name="Snow thickness in thickness categories"                       unit="m"  grid_ref="grid_T_3D_ncatice" />
@@ -140 +140 @@
          <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degC"   grid_ref="grid_T_3D_ncatice" />
          <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degC"   grid_ref="grid_T_3D_ncatice" />
-         <field id="iceamp_cat"   long_name="Ice melt pond fraction for categories"                        unit="%"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="iceamp_cat"   long_name="Ice melt pond concentration for categories"                   unit="%"      grid_ref="grid_T_3D_ncatice" /> 
          <field id="icevmp_cat"   long_name="Ice melt pond volume for categories"                          unit="m"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="icehmp_cat"   long_name="Ice melt pond thickness for categories"                       unit="m"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="iceafp_cat"   long_name="Ice melt pond fraction for categories"                        unit="m"      grid_ref="grid_T_3D_ncatice" /> 
 
          <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     />
@@ -186 +188 @@
          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   />
          <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   />
+         <field id="vfxpnd"       long_name="melt pond water flux to ocean"                                unit="kg/m2/s"   />
 
          <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   />

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

@@ -13 +13 @@
 !!             11 - Ice growth in open water           (namthd_do)
 !!             12 - Ice salinity                       (namthd_sal)
-!!             13 - Ice melt ponds                     (nammp)
+!!             13 - Ice melt ponds                     (namthd_pnd)
 !!             14 - Ice initialization                 (namini)
 !!             15 - Ice/snow albedos                   (namalb)
@@ -20 +20 @@
 !
 !------------------------------------------------------------------------------
-&nampar     !   Generic parameters
+&nampar         !   Generic parameters
 !------------------------------------------------------------------------------
    jpl              =   5             !  number of ice  categories
@@ -39 +39 @@
 /
 !------------------------------------------------------------------------------
-&namitd     !   Ice discretization
-!------------------------------------------------------------------------------
-   rn_himean        =   2.0           !  expected domain-average ice thickness (m)
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05)
+      rn_himean     =   2.0           !  expected domain-average ice thickness (m)
+   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m)
+      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0  
    rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping
 /
 !------------------------------------------------------------------------------
-&namdyn     !   Ice dynamics
+&namdyn         !   Ice dynamics
 !------------------------------------------------------------------------------
    ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
@@ -90 +93 @@
 !------------------------------------------------------------------------------
    ln_rhg_EVP       = .true.          !  EVP rheology
-      rn_creepl     =   1.0e-12       !     creep limit [1/s]
+      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017)
+      rn_creepl     =   2.0e-9        !     creep limit [1/s]
       rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
       nn_nevp       = 120             !     number of EVP subcycles                             
@@ -118 +122 @@
                                       !     = 2  Redistribute a single flux over categories
                                       !          ==> coupled mode only
-/
-!------------------------------------------------------------------------------
-&namthd     !   Ice thermodynamics
+   nice_jules       = 1               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
+/
+!------------------------------------------------------------------------------
+&namthd         !   Ice thermodynamics
 !------------------------------------------------------------------------------
    ln_icedH         = .true.          !  activate ice thickness change from growing/melting (T) or not (F)
    ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F)
    ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
-   ln_icedS         = .true.          !  activate gravity drainage and flushing (T) or not (F)
+   ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
 /
 !------------------------------------------------------------------------------
@@ -136 +141 @@
                                       !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
    rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
-   ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F)
 /
 !------------------------------------------------------------------------------
@@ -150 +154 @@
 /
 !------------------------------------------------------------------------------
-&namthd_do     !   Ice growth in open water
+&namthd_do      !   Ice growth in open water
 !------------------------------------------------------------------------------
    rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
@@ -174 +178 @@
 /
 !------------------------------------------------------------------------------
-&nammp      !   Melt ponds
-!------------------------------------------------------------------------------
-   ln_pnd           = .false.         !  active melt ponds
-   ln_pnd_rad       = .false.         !  active melt ponds radiative coupling
-   ln_pnd_fw        = .false.         !  active melt ponds freshwater coupling
-   nn_pnd_scheme    =   0             !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic
-   rn_apnd          =   0.2           !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0)
-   rn_hpnd          =   0.05          !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0)
-/
-!------------------------------------------------------------------------------
-&namini     !   Ice initialization
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+      ln_pnd_fwb    = .false.         !     melt ponds store freshwater or not
+   ln_pnd_CST       = .false.         !  activate constant melt ponds
+      rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0
+      rn_hpnd       =   0.05          !     prescribed pond depth, at Tsu=0
+   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
@@ -209 +213 @@
 /
 !------------------------------------------------------------------------------
-&namalb     !   albedo parameters
-!------------------------------------------------------------------------------
-   nn_ice_alb       =    1            !  parameterization of ice/snow albedo
-                                      !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
-                                      !     1: "home made" based on Brandt et al. (JClim 2005)
-                                      !                         and Grenfell & Perovich (JGR 2004), giving cloud-sky albedo
-                                      !     2: as 1 with melt ponds
-   rn_alb_sdry      =   0.85          !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
-   rn_alb_smlt      =   0.75          !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
-   rn_alb_idry      =   0.60          !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
-   rn_alb_imlt      =   0.50          !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
-   rn_alb_dpnd      =   0.27          !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' )
-/
-!------------------------------------------------------------------------------
-&namdia     !   Diagnostics
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+   !                                  !                          !  obs range (cloud-sky)
+   rn_alb_sdry      =   0.85          !  dry snow albedo         :  0.85 -- 0.87
+   rn_alb_smlt      =   0.75          !  melting snow albedo     :  0.72 -- 0.82
+   rn_alb_idry      =   0.60          !  dry ice albedo          :  0.54 -- 0.65
+   rn_alb_imlt      =   0.50          !  bare puddled ice albedo :  0.49 -- 0.58
+   rn_alb_dpnd      =   0.27          !  ponded ice albedo       :  0.10 -- 0.30 
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
 !------------------------------------------------------------------------------
    ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F)

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_ice_cfg

@@ -1 +1 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! LIM3 namelist:  
-!!              1 - Generic parameters                 (namicerun)
-!!              2 - Diagnostics                        (namicediag)
-!!              3 - Ice initialization                 (namiceini)
-!!              4 - Ice discretization                 (namiceitd)
-!!              5 - Ice dynamics and transport         (namicedyn)
-!!              6 - Ice thermodynamics                 (namicethd)
-!!              7 - Ice salinity                       (namicesal)
-!!              8 - Ice mechanical redistribution      (namiceitdme)
-!!              9 - Ice/snow albedos                   (namicealb)
+!! ESIM namelist:  
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface forcing                (namforcing)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (namthd_pnd)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !
 !------------------------------------------------------------------------------
-&namicerun     !   Generic parameters
+&nampar         !   Generic parameters
 !------------------------------------------------------------------------------
-   jpl              =    1          !  number of ice  categories
-   nlay_i           =    1          !  number of ice  layers
-   ln_limthd        =  .false.       !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   ln_limdyn        =  .true.       !  ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   nn_limdyn        =   0           !     (ln_limdyn=T) switch for ice dynamics
-                                    !      2: total
-                                    !      1: advection only (no diffusion, no ridging/rafting)
-                                    !      0: advection only (as 1 but with prescribed velocity, bypass rheology)
-   rn_uice          =   0.5    !     (nn_limdyn=0) ice u-velocity
-   rn_vice          =   0.0    !     (nn_limdyn=0) ice v-velocity
+   jpl              =   1             !  number of ice  categories
+   nlay_i           =   1             !  number of ice  layers
+   ln_icedyn        = .true.          !  ice dynamics (T) or not (F)
+   ln_icethd        = .false.         !  ice thermo   (T) or not (F)
 /
 !------------------------------------------------------------------------------
-&namicediag    !   Diagnostics
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping
+/
+!------------------------------------------------------------------------------
+&namdyn         !   Ice dynamics
+!------------------------------------------------------------------------------
+   ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
+   ln_dynADV        = .true.          !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+      rn_uice       =   0.5           !        prescribed ice u-velocity
+      rn_vice       =   0.            !        prescribed ice v-velocity
+/
+!------------------------------------------------------------------------------
+&namdyn_rdgrft  !   Ice ridging/rafting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namiceini     !   Ice initialization
-!------------------------------------------------------------------------------
-                  ! -- limistate -- !
-   ln_limini      = .true.         !  activate ice initialization (T) or not (F)
-   ln_limini_file = .true.         !  netcdf file provided for initialization (T) or not (F)
-   cn_dir="./"
-   sn_hti = 'initice'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-   sn_hts = 'initice'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-   sn_ati = 'initice'    , -12 ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-   sn_tsu = 'initice'    , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-   sn_tmi = 'initice'    , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-   sn_smi = 'initice'    , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
-/
-!------------------------------------------------------------------------------
-&namiceitd     !   Ice discretization
+&namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicedyn     !   Ice dynamics and transport
+&namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
-                  ! -- limtrp & limadv -- !
-   nn_limadv      =    0            !  choose the advection scheme (-1=Prather ; 0=Ultimate-Macho)
-   nn_limadv_ord  =    5            !  choose the order of the advection scheme (if nn_limadv=0)
 /
 !------------------------------------------------------------------------------
-&namicethd     !   Ice thermodynamics
+&namforcing     !   Ice surface forcing
 !------------------------------------------------------------------------------
-                  ! -- limthd_dh -- !
-   ln_limdH       = .true.          !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-                  ! -- limthd_da -- !
-   ln_limdA       = .true.          !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-                 ! -- limthd_lac -- !
-   ln_limdO       = .true.          !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   rn_hnewice     = 0.02            !  thickness for new ice formation in open water (m)
-                  ! -- limitd_th -- !
-   rn_himin       = 0.01            !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice
 /
 !------------------------------------------------------------------------------
-&namicesal     !   Ice salinity
+&namthd         !   Ice thermodynamics
 !------------------------------------------------------------------------------
-                 ! -- limthd_sal -- !
-   ln_limdS       = .true.          !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
 /
 !------------------------------------------------------------------------------
-&namiceitdme   !   Ice mechanical redistribution (ridging and rafting)
+&namthd_zdf     !   Ice heat diffusion
 !------------------------------------------------------------------------------
-                  ! -- limitd_me -- !
-   ln_ridging     =   .true.        !  ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   ln_rafting     =   .true.        !  rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
 /
-!-----------------------------------------------------------------------
-&namicealb     !   albedo parameters
-!-----------------------------------------------------------------------
+!------------------------------------------------------------------------------
+&namthd_da      !   Ice lateral melting
+!------------------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&namthd_do      !   Ice growth in open water
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+/

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

@@ -178 +178 @@
    !
    !                                     !!** ice-rheology namelist (namrhg) **
+   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
@@ -195 +196 @@
    !                                      !   = 2  Redistribute a single flux over categories
 
+   INTEGER , PUBLIC            ::   nice_jules            ! Choice of jules coupling 
+   !                                                      ! Associated indices
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   ! no Jules coupling (ice thermodynamics forced via qsr and qns)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   ! emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   ! active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
+
    !                                     !!** ice-salinity namelist (namthd_sal) **
    INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model
@@ -204 +211 @@
    REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]
 
-   ! MV MP 2016
-   !                                     !!** melt pond namelist (nammp)
-   LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not
-   LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not
-   LOGICAL , PUBLIC ::   ln_pnd_fw        !: ponds active wrt meltwater or not
-   INTEGER , PUBLIC ::   nn_pnd_scheme    !: type of melt pond scheme:   =0 prescribed, =1 empirical, =2 topographic
-   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1), only if nn_pnd_scheme = 0
-   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1), only if nn_pnd_scheme = 0
-   ! END MV MP 2016
+   !                                     !!** namelist namthd_pnd
+   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012
+   LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater
+   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
+   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
+   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1)
+   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo
+
    !                                     !!** ice-diagnostics namelist (namdia) **
    LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

@@ -36 +36 @@
    !
    ! ** albedo namelist (namalb)
-   INTEGER  ::   nn_ice_alb       ! type of albedo scheme: 0: Shine & Henderson-Sellers (JGR 1985)
-   !                                      !                         1: "home made" based on Brandt et al. (JClim 2005)
-   !                                      !                                             and Grenfell & Perovich (JGR 2004)
-   !                                      !                         2: Same as 1 but with melt ponds
    REAL(wp) ::   rn_alb_sdry      ! dry snow albedo
    REAL(wp) ::   rn_alb_smlt      ! melting snow albedo
@@ -53 +49 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -59 +55 @@
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
       !!       
-      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
-      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
-      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
-      !!                                                                           and Grenfell & Perovich (JGR 2004)
-      !!                  2: fractional surface-based formulation of scheme 1 (NEW)
-      !!                Description of scheme 1:
+      !! ** Method  :   The scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                      and Grenfell & Perovich (JGR 2004)
       !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
       !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
@@ -76 +68 @@
       !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
       !!
-      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                     compilation of values from literature (reference overcast sky values)
+      !!                        rn_alb_sdry = 0.85      ! dry snow
+      !!                        rn_alb_smlt = 0.75      ! melting snow
+      !!                        rn_alb_idry = 0.60      ! bare frozen ice
+      !!                        rn_alb_imlt = 0.50      ! bare puddled ice albedo
+      !!                        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
+      !!                                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
+      !!                     Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+      !!                        rn_alb_sdry = 0.85      ! dry snow
+      !!                        rn_alb_smlt = 0.72      ! melting snow
+      !!                        rn_alb_idry = 0.65      ! bare frozen ice
+      !!                     Brandt et al 2005 (East Antarctica)
+      !!                        rn_alb_sdry = 0.87      ! dry snow
+      !!                        rn_alb_smlt = 0.82      ! melting snow
+      !!                        rn_alb_idry = 0.54      ! bare frozen ice
+      !!
+      !! ** Note    :   The old parameterization from Shine & Henderson-Sellers (not here anymore) presented several misconstructions
       !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
       !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
@@ -87 +95 @@
       !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice       !  ice surface temperature (Kelvin)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_su        !  ice surface temperature (Kelvin)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth
+      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
-      LOGICAL , INTENT(in   )                   ::   ld_pnd       !  melt ponds radiatively active or not
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs    !  albedo of ice under clear    sky
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os    !  albedo of ice under overcast sky
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky
       !
       INTEGER  ::   ji, jj, jl                ! dummy loop indices
-      REAL(wp) ::   zswitch, z1_c1, z1_c2     ! local scalar
-      REAL(wp) ::   z1_href_pnd               !   -      -                 
-      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
-      !
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free)
-!! MV MP
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd        ! ponded sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice        ! bare sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw        ! snow-covered sea ice albedo
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw      ! relative snow fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice      ! relative ice fraction
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd      ! relative ice fraction (effective)
+      REAL(wp) ::   z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar
+      REAL(wp) ::   z1_href_pnd               ! inverse of the characteristic length scale (Lecomte et al. 2015)
+      REAL(wp) ::   zalb_pnd, zafrac_pnd      ! ponded sea ice albedo & relative pound fraction
+      REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction
+      REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('icealb')
       !
-      !-----------------------------------------------------
-      !  Snow-free albedo (no ice thickness dependence yet)
-      !-----------------------------------------------------
-      !
-      ! Part common to nn_ice_alb = 0, 1, 2
-      !
-      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
-         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   ;   zalb(:,:,:) = rn_alb_imlt
-         ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry
-         END WHERE
-      ENDIF
-
-      SELECT CASE ( nn_ice_alb )       ! select a parameterization
-      !
-      !              !------------------------------------------
-      CASE( 0 )      !  Shine and Henderson-Sellers (1985)     !   (based on clear sky values)
-         !           !------------------------------------------
-         !
-         !                       ! Thickness-dependent bare ice albedo
-         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
-         ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
-         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              &
-            &                                                                  - 0.8608 * ph_ice * ph_ice     &
-            &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice
-         ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
-            z1_href_pnd = 1. / 0.05    ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
-            !                          ! Snow-free ice albedo is a function of pond fraction
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   
-               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
-            END WHERE
-         ENDIF 
-         !
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !                    ! Freezing snow
-                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > ppc2
-                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - ppc1 ) ) )
-                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * pp1_c1  )   &
-                     &    +           zswitch   * rn_alb_sdry  
-                     !
-                  !                    ! Melting snow
-                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > ppc2
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - ppc2 ) )
-                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * pp1_c2 )   &
-                     &      +         zswitch   *   rn_alb_smlt 
-                     !
-                  !                    ! Snow albedo
-                  zswitch =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
-                  zalb_st =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-                  !
-                  !                    ! Surface albedo
-                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
-                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
-                  !
-               END DO
+      z1_href_pnd = 0.05
+      z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+      z1_c2 = 1. / 0.05
+      z1_c3 = 1. / 0.02
+      z1_c4 = 1. / 0.03
+      !
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !                       !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time)
+               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN
+                  zafrac_snw = 0._wp
+                  IF( ld_pnd_alb ) THEN
+                     zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                  ELSE
+                     zafrac_pnd = 0._wp
+                  ENDIF
+                  zafrac_ice = 1._wp - zafrac_pnd
+               ELSE
+                  zafrac_snw = 1._wp      ! Snow fully "shades" melt ponds and ice
+                  zafrac_pnd = 0._wp
+                  zafrac_ice = 0._wp
+               ENDIF
+               !
+               !                       !--- Bare ice albedo (for hi > 150cm)
+               IF( ld_pnd_alb ) THEN
+                  zalb_ice = rn_alb_idry
+               ELSE
+                  IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN  ;   zalb_ice = rn_alb_imlt
+                  ELSE                                                               ;   zalb_ice = rn_alb_idry   ;   ENDIF
+               ENDIF
+               !                       !--- Bare ice albedo (for hi < 150cm)
+               IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN      ! 5cm < hi < 150cm
+                  zalb_ice = zalb_ice    + ( 0.18 - zalb_ice   ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice(ji,jj,jl)) )
+               ELSEIF( ph_ice(ji,jj,jl) <= 0.05 ) THEN                               ! 0cm < hi < 5cm
+                  zalb_ice = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl)
+               ENDIF
+               !
+               !                       !--- Snow-covered ice albedo (freezing, melting cases)
+               IF( pt_su(ji,jj,jl) < rt0_snow ) THEN
+                  zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 )
+               ELSE
+                  zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c4 )
+               ENDIF
+               !                       !--- Ponded ice albedo
+               IF( ld_pnd_alb ) THEN
+                  zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd ) 
+               ELSE
+                  zalb_pnd = rn_alb_dpnd
+               ENDIF
+               !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
+               palb_os(ji,jj,jl) = zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice
+               !
             END DO
          END DO
-         !
-         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + ppcloud       ! Oberhuber correction for overcast sky
-         !
-         !
-         !           !------------------------------------------
-      CASE( 1 )      !  New parameterization (2016)            !    (based on overcast sky values)
-         !           !------------------------------------------
-         !
-         !                    compilation of values from literature (reference overcast sky values)
-         !                          rn_alb_sdry = 0.85      ! dry snow
-         !                          rn_alb_smlt = 0.75      ! melting snow
-         !                          rn_alb_idry = 0.60      ! bare frozen ice
-         !                          rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015)
-         !                                                  ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich
-         !                    Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
-         !                          rn_alb_sdry = 0.85      ! dry snow
-         !                          rn_alb_smlt = 0.72      ! melting snow
-         !                          rn_alb_idry = 0.65      ! bare frozen ice
-         !                    Brandt et al 2005 (East Antarctica)
-         !                          rn_alb_sdry = 0.87      ! dry snow
-         !                          rn_alb_smlt = 0.82      ! melting snow
-         !                          rn_alb_idry = 0.54      ! bare frozen ice
-         ! 
-         !              !--- Computation of snow-free ice albedo 
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
-         z1_c2 = 1. / 0.05
-
-         !              !--- Accounting for the ice-thickness dependency
-         WHERE     ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice) )
-         ELSE WHERE                                        ;   zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo
-            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
-            !
-            !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )
-               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd
-            END WHERE
-         ENDIF 
-         !
-         !              !--- Overcast sky surface albedo (accounting for snow, ice melt ponds)
-         z1_c1 = 1. / 0.02
-         z1_c2 = 1. / 0.03
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! Snow depth dependence of snow albedo
-                  zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
-                  zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
-
-                  ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow)
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
-                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-
-                  ! Surface albedo   
-                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
-                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
-
-              END DO
-            END DO
-         END DO
-         !
-         !              !--- Clear sky surface albedo
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
-         !
-         !
-         !           !------------------------------------------
-      CASE( 2 )      !  Fractional surface-based parameterization (2017)
-         !           !------------------------------------------
-         !              MV: I propose this formulation that is more elegant, and more easy to expand towards
-         !              varying pond and snow fraction.
-         !              Formulation 1 has issues to handle ponds and snow together that can't easily be fixed. 
-         !              This one handles it better, I believe.
-         !
-         !                 !==  Snow, bare ice and ponded ice fractions  ==!
-         !
-         !                       ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category
-         !
-         !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
-         zafrac_pnd = 0._wp
-         IF ( ld_pnd ) THEN
-            WHERE( ph_snw == 0._wp )   zafrac_pnd = pafrac_pnd   ! Snow fully "shades" melt ponds
-         ENDIF         
-         !
-         !                       !--- Specific snow fraction (for now, prescribed)
-         WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1.
-         ELSEWHERE                 ;   zafrac_snw = 0.
-         END WHERE
-         !
-         !                       !--- Specific ice fraction 
-         zafrac_ice = 1. - zafrac_snw - zafrac_pnd
-         !
-         !                 !==  Snow-covered, pond-covered, and bare ice albedos  ==!
-         !
-         !                       !--- Bare ice albedo
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
-         z1_c2 = 1. / 0.05
-         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_ice = zalb
-         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  &
-            &                                                                       ( LOG(1.5) - LOG(ph_ice) )
-         ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
-         END WHERE
-         !
-         z1_c1 = 1. / 0.02       !--- Snow-covered ice albedo (freezing, melting cases)
-         z1_c2 = 1. / 0.03
-         !
-         WHERE( pt_ice < rt0_snow )   ;   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 )
-         ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 )
-         END WHERE
-         !
-         IF ( ld_pnd ) THEN      !--- Depth-dependent ponded ice albedo
-            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015)
-            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd * z1_href_pnd ) 
-         ELSE
-            zalb_pnd  = rn_alb_dpnd
-         ENDIF
-         !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions
-         pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 )
-
-      END SELECT
+      END DO
+      !
+      palb_cs(:,:,:) = palb_os(:,:,:) - ( - 0.1010 * palb_os(:,:,:) * palb_os(:,:,:) + 0.1933 * palb_os(:,:,:) - 0.0148 )
       !
       IF( nn_timing == 1 )   CALL timing_stop('icealb')
@@ -317 +188 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namalb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
+      NAMELIST/namalb/ rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
       !!----------------------------------------------------------------------
       !
@@ -334 +205 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namalb:'
-         WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
          WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
          WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

@@ -76 +76 @@
       IF( icount == 0 ) THEN
          !                          ! water flux
-         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    &
+         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
+            &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
+            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
+            &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  &
             &                  ) * e1e2t(:,:) ) * zconv
          !
@@ -96 +97 @@
 
          pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
-            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv
+            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv
 
       ELSEIF( icount == 1 ) THEN
 
          ! water flux
-         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
+         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
+            &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
+            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
+            &                wfx_ice_sub(:,:) + wfx_spr(:,:)  &
             &              ) * e1e2t(:,:) ) * zconv - pdiag_fv
 
@@ -117 +119 @@
  
          ! outputs
-         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  &
+         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_v ) * r1_rdtice - zfv ) * rday
 
-         zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  &
+         zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_s ) * r1_rdtice + zfs ) * rday
 

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

@@ -172 +172 @@
             END DO
          END DO
-      END DO
-            
-      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
-         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      ENDIF
+      END DO 
+      !                                           !-- correct pond fraction to avoid a_ip > a_i
+      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
       !
    END SUBROUTINE Hbig

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90

@@ -120 +120 @@
             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
          END DO
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
@@ -167 +167 @@
                      &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
@@ -220 +220 @@
                      &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
@@ -249 +249 @@
          END DO
          ! MV MP 2016
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
             pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
@@ -755 +755 @@
                   sxyage(:,:,jl)= z2d(:,:)
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   DO jl = 1, jpl 
                      WRITE(zchar,'(I2.2)') jl
@@ -833 +833 @@
                syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
                sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
                   syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
@@ -854 +854 @@
             syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
             sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-            IF ( nn_pnd_scheme > 0 ) THEN
+            IF ( ln_pnd_H12 ) THEN
                sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
                syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
@@ -989 +989 @@
             END DO
          END DO
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             DO jl = 1, jpl 
                WRITE(zchar,'(I2.2)') jl

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90

@@ -124 +124 @@
             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
-            IF ( nn_pnd_scheme > 0 ) THEN
+            IF ( ln_pnd_H12 ) THEN
                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

@@ -253 +253 @@
          CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-            CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         ENDIF
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -270 +268 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-         ENDIF
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
          
          !-----------------------------------------------------------------------------!
@@ -320 +316 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-            CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         ENDIF
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -337 +331 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-         ENDIF
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
 
       ENDIF ! npti > 0
@@ -651 +643 @@
 
                !MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   aprdg1     = a_ip_2d(ji,jl1) * afrdg
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
@@ -679 +671 @@
                ! Put the melt pond water into the ocean
                !------------------------------------------            
-               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
+               IF ( ln_pnd_fwb ) THEN
                   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
                      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
@@ -702 +694 @@
                oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
                ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
@@ -766 +758 @@
                      &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) )
                   ! MV MP 2016
-                  IF ( nn_pnd_scheme > 0 ) THEN
+                  IF ( ln_pnd_H12 ) THEN
                      v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
                         &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90

@@ -78 +78 @@
       CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic !
          !                             !------------------------!
-         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
       END SELECT
@@ -107 +107 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
       !!-------------------------------------------------------------------
       !
@@ -125 +125 @@
          WRITE(numout,*) '   Namelist namdyn_rhg:'
          WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP
+         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP
          WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl
          WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90

@@ -53 +53 @@
 CONTAINS
 
-   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice, pshear_i, pdivu_i, pdelta_i )
+   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
@@ -95 +95 @@
       !!              5) Diagnostics including charge ellipse
       !!
-      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
-      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
-      !!              but this solution appears very unstable (see Kimmritz et al 2016)
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
       !!
       !! References : Hunke and Dukowicz, JPO97
       !!              Bouillon et al., Ocean Modelling 2009
       !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! time step
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   u_ice, v_ice 
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i, pdivu_i, pdelta_i 
       !!
@@ -114 +115 @@
       REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
       REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
-      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
+      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                       ! alpha coef from Bouillon 2009 or Kimmritz 2017
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
@@ -126 +127 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
-      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
@@ -231 +234 @@
 
       ! alpha parameters (Bouillon 2009)
-      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
-      zalph2 = zalph1 * z1_ecc2
-
-      ! alpha and beta parameters (Bouillon 2013)
-      !!zalph1 = 40.
-      !!zalph2 = 40.
-      !!zbeta  = 3000.
-      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
-
-      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ENDIF
+         
       ! Initialise stress tensor 
       zs1 (:,:) = pstress1_i (:,:) 
@@ -304 +303 @@
             zmf(ji,jj)      = zm1 * ff_t(ji,jj)
 
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
             ! m/dt
             zmU_t(ji,jj)    = zmassU * z1_dtevp
             zmV_t(ji,jj)    = zmassV * z1_dtevp
-
+            
             ! Drag ice-atm.
             zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
@@ -326 +328 @@
          END DO
       END DO
-      CALL lbc_lnk( zmf, 'T', 1. )
+      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
       !
       !------------------------------------------------------------------------------!
@@ -380 +382 @@
                ! P/delta at T points
                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+                  zalph2   = zalph1
+                  z1_alph2 = z1_alph1
+               ENDIF
                
                ! stress at T points
@@ -392 +404 @@
             DO ji = 1, jpim1
 
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  zbeta(ji,jj) = zalph2
+               ENDIF
+               
                ! P/delta at F points
                zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
@@ -411 +430 @@
                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
                   &                  ) * r1_e1e2u(ji,jj)
-                  !
-                  !                !--- V points
+               !
+               !                !--- V points
                zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
                   &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
@@ -419 +438 @@
                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
                   &                  ) * r1_e1e2v(ji,jj)
-                  !
-                  !                !--- u_ice at V point
+               !
+               !                !--- u_ice at V point
                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
-                  !
-                  !                !--- v_ice at U point
+               !
+               !                !--- v_ice at U point
                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
-               !
             END DO
          END DO
          !
          ! --- Computation of ice velocity --- !
-         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
          !  Bouillon et al. 2009 (eq 34-35) => stable
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
@@ -438 +456 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                     !                 !--- tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
-                     !                 !--- Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-                     !
-                     !                 !--- tau_bottom/v_ice
+                  !
+                  !                 !--- tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-                     !
-                     !                 !--- Coriolis at V-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
-                     !
-                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-                     !
-                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                     !
-                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -466 +491 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                     !
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
-                  !
+                  ENDIF
                END DO
             END DO
@@ -483 +503 @@
             !
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                               
-                  ! tau_io/(u_oce - u_ice)
+               DO ji = fs_2, fs_jpim1          
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-                  
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -514 +539 @@
                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -532 +554 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-
-                  ! tau_io/(u_oce - u_ice)
+                  !                 !--- tau_io/(u_oce - u_ice)
                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
+                  !
+                  !                 !--- tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
                   zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin 
+                     &            ) * zmaskU(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
+                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                      &            ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+                  ENDIF
                END DO
             END DO
@@ -578 +602 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-         
-                  ! tau_io/(v_oce - v_ice)
+                  !                 !--- tau_io/(v_oce - v_ice)
                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress
                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
-                  ! tau_bottom/v_ice
+                  !
+                  !                 !--- tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
                   ztauB = - tau_icebfr(ji,jj) / zvel
-                  
-                  ! Coriolis at V-points (energy conserving formulation)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
-
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
+                  !
+                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                  
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
+                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
@@ -608 +637 @@
                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                      &            ) * zmaskV(ji,jj)
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
+                  ENDIF
                END DO
             END DO

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

@@ -99 +99 @@
       !!
       !! ** Action  : It provides the following fields used in sea ice model:
-      !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
       !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
       !!                sprecip                                  = solid precipitation                           [Kg/m2/s]
@@ -130 +129 @@
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
@@ -142 +141 @@
          !
       CASE( jp_blk )                !--- bulk formulation
-                                CALL blk_ice_flx( t_su, alb_ice )    ! 
-         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
+                                CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    ! 
+         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_iceflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
       CASE ( jp_purecpl )           !--- coupled formulation
-                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
+                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
          IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
@@ -268 +267 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx
+      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx, nice_jules
       !!-------------------------------------------------------------------
       !
@@ -285 +284 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namforcing:'
-         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio    = ', rn_cio
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s = ', rn_blow_s
-         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx = ', nn_iceflx
+         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio     = ', rn_cio
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s  = ', rn_blow_s
+         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx  = ', nn_iceflx
+         WRITE(numout,*) '      Jules coupling (0=no, 1=emulated, 2=active)      nice_jules = ', nice_jules
       ENDIF
       !

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

@@ -94 +94 @@
       REAL(wp) ::   ztmelts, zdh
       INTEGER  ::   i_hemis, i_fill, jl0
-      REAL(wp) ::   zarg, zV, zconv, zdv
+      REAL(wp) ::   zarg, zV, zconv, zdv, zfac
       INTEGER , DIMENSION(4)           ::   itest
       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
@@ -314 +314 @@
 
          ! for constant salinity in time
-         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+         IF( nn_icesal /= 2 )  THEN
             CALL ice_var_salprof
             sv_i = s_i * v_i
@@ -358 +358 @@
          tn_ice (:,:,:) = t_su (:,:,:)
 
-         ! MV MP 2016
          ! Melt pond volume and fraction
-         IF ( ln_pnd ) THEN
-            DO jl = 1, jpl
-               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
-               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
-               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
-               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
-            END DO
-         ELSE
-            a_ip(:,:,:)      = 0._wp
-            v_ip(:,:,:)      = 0._wp
-            a_ip_frac(:,:,:) = 0._wp
-            h_ip     (:,:,:) = 0._wp
-         ENDIF
-         ! END MV MP 2016
+         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
+         ELSE                                     ;   zfac = 0._wp   ;   ENDIF 
+         DO jl = 1, jpl
+            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac
+            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac
+         END DO
+         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
+         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
 
       ELSE ! if ln_iceini=false

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

@@ -36 +36 @@
    PUBLIC   ice_itd_reb   ! called in icecor
 
+   INTEGER ::              nice_catbnd   ! choice of the type of ice category function
+   !                                       ! associated indices:
+   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function
+   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user
+   !
    ! ** namelist (namitd) **
+   LOGICAL  ::   ln_cat_hfn  ! ice categories are defined by function like rn_himean**(-0.05)
    REAL(wp) ::   rn_himean   ! mean thickness of the domain
-
+   LOGICAL  ::   ln_cat_usr  ! ice categories are defined by rn_catbnd
+   REAL(wp), DIMENSION(0:100) ::   rn_catbnd   ! ice categories bounds
+   !
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2017)
@@ -293 +301 @@
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
-               ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
-                  a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
-               ENDIF
-               ! END MV MP 2016
+               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
                h_i_1d(ji) = rn_himin
             ENDIF
@@ -457 +461 @@
                zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
                !  
-               ! MV MP 2016 
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   !                                                  ! Pond fraction
                   ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
@@ -468 +471 @@
                   v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
                ENDIF
-               ! END MV MP 2016
                !
             ENDIF   ! jl1 >0
@@ -643 +645 @@
       !! ** input   :   Namelist namitd
       !!-------------------------------------------------------------------
-      INTEGER  ::   jl    ! dummy loop index
-      INTEGER  ::   ios   ! Local integer output status for namelist read
+      INTEGER  ::   jl            ! dummy loop index
+      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
       REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
-      !!
-      NAMELIST/namitd/ rn_himean, rn_himin
+      !
+      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
       !!------------------------------------------------------------------
       !
@@ -664 +666 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namitd: '
-         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
-         WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin 
+         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
+         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
+         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
+         WRITE(numout,*) '         ice categories boundaries (m)                                  rn_catbnd  = ', rn_catbnd 
+         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
       ENDIF
       !
@@ -671 +676 @@
       !  Thickness categories boundaries  !
       !-----------------------------------!
-      !
-      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function)
-      zhmax  = 3._wp * rn_himean
-      DO jl = 1, jpl
-         znum = jpl * ( zhmax+1 )**zalpha
-         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
-         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
-      END DO
+      !                             !== set the choice of ice categories ==!
+      ioptio = 0 
+      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
+      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
+      !
+      SELECT CASE( nice_catbnd )
+      !                                !------------------------!
+      CASE( np_cathfn )                ! h^(-alpha) function
+         !                             !------------------------!
+         zalpha = 0.05_wp
+         zhmax  = 3._wp * rn_himean
+         DO jl = 1, jpl
+            znum = jpl * ( zhmax+1 )**zalpha
+            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
+            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
+         END DO
+         !                             !------------------------!
+      CASE( np_catusr )                ! user defined
+         !                             !------------------------!
+         DO jl = 0, jpl
+            hi_max(jl) = rn_catbnd(jl)
+         END DO
+         !
+      END SELECT
       !
       DO jl = 1, jpl                ! mean thickness by category

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90

@@ -147 +147 @@
       END DO
 
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'a_ip'//'_htc'//zchar
-            z2d(:,:) = a_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
-            znam = 'v_ip'//'_htc'//zchar
-            z2d(:,:) = v_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
-         END DO
-      ENDIF
-      ! END MV MP 2016
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I2.2)') jl
+         znam = 'a_ip'//'_htc'//zchar
+         z2d(:,:) = a_ip(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
+         znam = 'v_ip'//'_htc'//zchar
+         z2d(:,:) = v_ip(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
+      END DO
 
       DO jl = 1, jpl 
@@ -193 +189 @@
       !! ** purpose  :   read restart file
       !!----------------------------------------------------------------------
-      INTEGER  :: jk, jl
-      REAL(wp) :: zfice, ziter
+      INTEGER  ::   jk, jl
+      INTEGER  ::   id1            ! local integer
+      REAL(wp) ::   zfice, ziter
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d
       CHARACTER(len=25) ::   znam
@@ -248 +245 @@
       END DO
 
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
+      id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. )
+      IF( id1 > 0 ) THEN                       ! fields exist (melt ponds)
          DO jl = 1, jpl 
             WRITE(zchar,'(I2.2)') jl
@@ -259 +256 @@
             v_ip(:,:,jl) = z2d(:,:)
          END DO
-      ENDIF
-      ! END MV MP 2016
+      ELSE                                     ! start from rest
+         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it'
+         a_ip(:,:,:) = 0._wp
+         v_ip(:,:,:) = 0._wp
+      ENDIF
 
       DO jl = 1, jpl 

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

@@ -35 +35 @@
    USE icedyn         ! sea-ice: dynamics
    USE icethd         ! sea-ice: thermodynamics
-   USE limmp          ! sea-ice: melt ponds
    USE icecor         ! sea-ice: corrections
    USE iceupdate      ! sea-ice: sea surface boundary condition update
@@ -154 +153 @@
          !------------------------------------------------------!
          ! It provides the following fields used in sea ice model:
-         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%]
          !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
          !    sprecip              = solid precipitation                           [Kg/m2/s]
@@ -170 +168 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-         IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds
          !
          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
@@ -238 +235 @@
       CALL ice_itd_init                ! ice thickness distribution initialization
       !
-      CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here)
-      !
+      IF( ln_icethd ) THEN
+         CALL ice_thd_init             ! set ice thermodynics parameters (clem: important to call it first for melt ponds)
+      ENDIF   
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
@@ -255 +253 @@
          CALL ice_dyn_init             ! set ice dynamics parameters
       ENDIF
-      !
-      IF( ln_icethd ) THEN
-         CALL ice_thd_init             ! set ice thermodynics parameters
-      ENDIF   
       !
       CALL ice_update_init             ! ice surface boundary condition
@@ -325 +319 @@
          IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
       ENDIF
-      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
-         CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
-      ENDIF
+!     IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
+!        CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
+!     ENDIF
       !
       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
@@ -398 +392 @@
       wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
       wfx_snw_sni(:,:) = 0._wp 
-      ! MV MP 2016
       wfx_pnd(:,:) = 0._wp
-      ! END MV MP 2016
 
       hfx_thd(:,:) = 0._wp   ;

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

@@ -26 +26 @@
    USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
    USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
-      &                 fr1_i0, fr2_i0
+      &                 qml_ice, qcn_ice, qsr_ice_tr
    USE ice1D          ! sea-ice: thermodynamics variables
    USE icethd_zdf     ! sea-ice: vertical heat diffusion
@@ -34 +34 @@
    USE icethd_ent     ! sea-ice: enthalpy redistribution
    USE icethd_do      ! sea-ice: growth in open water
+   USE icethd_pnd     ! sea-ice: melt ponds
    USE iceitd         ! sea-ice: remapping thickness distribution
    USE icetab         ! sea-ice: 1D <==> 2D transformation
@@ -86 +87 @@
       !!             - call ice_thd_rem  for remapping thickness distribution
       !!             - call ice_thd_do   for ice growth in leads
-      !!---------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt    ! number of iteration
       !
@@ -230 +231 @@
             s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here) 
             dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp
-            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp
+            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
             !
             IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
                               CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
+                              CALL ice_thd_pnd                             ! Melt ponds formation
                               CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
             ENDIF
@@ -362 +364 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
          CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)   )
-         END DO
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
+            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
+         END DO
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        )
          CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d  (1:npti), qsr_ice (:,:,kl) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d   (1:npti), fr1_i0           )
-         CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d   (1:npti), fr2_i0           )
          CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d  (1:npti), qns_ice (:,:,kl) )
          CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d  (1:npti), ftr_ice (:,:,kl) )
@@ -388 +391 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d    (1:npti), fhtur            )
          CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d     (1:npti), fhld             )
+         
+         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d   (1:npti), qml_ice      (:,:,kl)  )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d   (1:npti), qcn_ice      (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_tr_1d(1:npti), qsr_ice_tr   (:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   )
@@ -403 +410 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          )
@@ -454 +462 @@
          !
          ! Change thickness to volume (replaces routine ice_var_eqv2glo)
-         v_i_1d(1:npti)  = h_i_1d(1:npti) * a_i_1d(1:npti)
-         v_s_1d(1:npti)  = h_s_1d(1:npti) * a_i_1d(1:npti)
-         sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti)
+         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
+         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
+         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
+         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
          
          CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
          CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
-         END DO
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
+            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
+         END DO
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
@@ -488 +500 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )
@@ -523 +536 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
          !
       END SELECT
@@ -530 +544 @@
 
    SUBROUTINE ice_thd_init
-      !!-----------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_thd_init *** 
       !!                 
@@ -570 +584 @@
       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
+                       CALL ice_thd_pnd_init   ! set melt ponds parameters
       !
    END SUBROUTINE ice_thd_init

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

@@ -130 +130 @@
       !------------------------------------------------------------------------------!
       !
-      DO ji = 1, npti
-         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+      SELECT CASE ( nice_jules )
+
+         CASE ( np_jules_ACTIVE )
+
+            DO ji = 1, npti
+               zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+            END DO
+
+         CASE ( np_jules_EMULE )
+
+            DO ji = 1, npti
+               zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
+               qml_ice_1d(ji)   =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+               zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )
+            END DO
+
+         CASE ( np_jules_OFF ) 
+
+            DO ji = 1, npti
+               zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) ) 
+               zdum             =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+               zq_su(ji)        =   MAX( 0._wp, zdum             * rdt_ice )
+            END DO
+
+      END SELECT
+         
+      DO ji = 1, npti
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
-
-         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -148 +171 @@
             ! Contribution to mass flux
             wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
+            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji)
             ! updates
-            h_s_1d(ji)   = 0._wp
+            h_s_1d(ji)    = 0._wp
             e_s_1d (ji,1) = 0._wp
             t_s_1d (ji,1) = rt0
@@ -211 +235 @@
          ! snow melting only = water into the ocean (then without snow precip), >0
          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          ! updates available heat + precipitations after melting
          zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
@@ -233 +258 @@
             hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk)
             ! updates available heat + thickness
-            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
+            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
          END DO
@@ -571 +597 @@
          zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
          dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
-         h_s_1d   (ji)  = h_s_1d(ji)   + zdeltah(ji,1)
+         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1)
         
          zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
@@ -577 +603 @@
          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
-         wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          !    
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
@@ -611 +638 @@
 
          ! Case constant salinity in time: virtual salt flux to keep salinity constant
-         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+         IF( nn_icesal /= 2 )  THEN
             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean 

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

@@ -34 +34 @@
    LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
    LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
-   REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K]
    REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F)
-
+   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K]
+
+   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation
+   !                                           ! associated indices:
+   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999)
+
+   INTEGER , PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns
+   INTEGER , PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns
+   INTEGER , PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd
+   
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2017)
@@ -45 +52 @@
 CONTAINS
 
-   SUBROUTINE ice_thd_zdf
+      SUBROUTINE ice_thd_zdf
+      
+      !!
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_thd_zdf  ***
       !! ** Purpose :
+      !!           This chooses between the appropriate routine for the 
+      !!           computation of vertical diffusion
+      !!
+      !!-------------------------------------------------------------------
+      !!
+   
+      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver
+      
+                                       !-------------      
+         CASE( np_BL99 )               ! BL99 solver
+                                       !-------------
+         
+            IF ( nice_jules == np_jules_OFF ) THEN   ! No Jules coupler     
+
+               CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )
+                                       
+            ENDIF
+            
+            IF ( nice_jules == np_jules_EMULE ) THEN   ! Jules coupler is emulated
+            
+               CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )
+               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
+                                                   
+            ENDIF
+
+            IF ( nice_jules == np_jules_ACTIVE ) THEN   ! Jules coupler is emulated
+            
+               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
+                                                   
+            ENDIF
+            
+      END SELECT
+      
+   END SUBROUTINE ice_thd_zdf
+   
+   
+   
+   SUBROUTINE ice_thd_zdf_BL99(k_jules)
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
+      !! ** Purpose :
       !!           This routine determines the time evolution of snow and sea-ice 
-      !!           temperature profiles.
+      !!           temperature profiles, using the original Bitz and Lipscomb (1999) algorithm
+      !!
       !! ** Method  :
       !!           This is done by solving the heat equation diffusion with
@@ -83 +134 @@
       !!           total ice/snow thickness : h_i_1d, h_s_1d
       !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND)
+      
       INTEGER ::   ji, jk         ! spatial loop index
       INTEGER ::   jm             ! current reference number of equation
@@ -101 +154 @@
       REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
       REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   z1_hsu
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zcpi                      ! Ice specific heat
@@ -111 +163 @@
       REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
       REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
-      
+      
+      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
@@ -136 +187 @@
       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
+
+      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor
       
       ! Mono-category
@@ -166 +219 @@
       ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
       END WHERE
-      !
-      ! temperatures
-      ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration
+
+      ! Store initial temperatures and non solar heat fluxes
+      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
+
+         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1
+         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value
+         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux 
+         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
+
+         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
+
+      ENDIF   
+
       ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
       ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
-      t_su_1d(1:npti)   = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! necessary
-      !
+
       !-------------
       ! 2) Radiation
       !-------------
-      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
-      DO ji = 1, npti
-         ! --- Computation of i0 --- !
-         ! i0 describes the fraction of solar radiation which does not contribute
-         ! to the surface energy budget but rather penetrates inside the ice.
-         ! We assume that no radiation is transmitted through the snow
-         ! If there is no no snow
-         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
-         ! zftrice = io.qsr_ice       is below the surface 
-         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
-         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
-         zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )     
-         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
-
-         ! --- Solar radiation absorbed / transmitted at the surface --- !
-         !     Derivative of the non solar flux
-         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
-         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
-      END DO
-
       ! --- Transmission/absorption of solar radiation in the ice --- !
-      zradtr_s(1:npti,0) = zftrice(1:npti)
+!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value
+!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1
+
+!     DO ji = 1, npti
+
+!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )     
+
+!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque
+
+!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation 
+
+!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji)
+!              zftrice(ji) = qsr_ice_tr_1d(ji)
+!              i0(ji) = zfrqsr_tr_i
+
+!     END DO
+
+      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)
       DO jk = 1, nlay_s
          DO ji = 1, npti
@@ -209 +267 @@
       END DO
          
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )
       DO jk = 1, nlay_i 
          DO ji = 1, npti
@@ -330 +388 @@
             END DO
          END DO
-         !
-         !----------------------------
-         ! 6) surface flux computation
-         !----------------------------
-         IF ( ln_dqns_i ) THEN 
-            DO ji = 1, npti
-               ! update of the non solar flux according to the update in T_su
+         
+         !
+         !----------------------------------------!
+         !                                        !
+         !       JULES IF (OFF or SND MODE)       !
+         !                                        !
+         !----------------------------------------!
+         !
+         
+         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
+         
+            !
+            ! In OFF mode the original BL99 temperature computation is used
+            ! (with qsr_ice, qns_ice and dqns_ice as inputs)
+            !
+            ! In SND mode, the computation is required to compute the conduction fluxes
+            !
+         
+            !----------------------------
+            ! 6) surface flux computation
+            !----------------------------
+
+            DO ji = 1, npti
+            ! update of the non solar flux according to the update in T_su
                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
             END DO
-         ENDIF
-
-         DO ji = 1, npti
-            zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
-         END DO
-         !
-         !----------------------------
-         ! 7) tridiagonal system terms
-         !----------------------------
-         !!layer denotes the number of the layer in the snow or in the ice
-         !!jm denotes the reference number of the equation in the tridiagonal
-         !!system, terms of tridiagonal system are indexed as following :
-         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
-         !!ice interior terms (top equation has the same form as the others)
-         ztrid   (1:npti,:,:) = 0._wp
-         zindterm(1:npti,:)   = 0._wp
-         zindtbis(1:npti,:)   = 0._wp
-         zdiagbis(1:npti,:)   = 0._wp
-
-         DO jm = nlay_s + 2, nlay_s + nlay_i 
+
+            DO ji = 1, npti
+               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar
+            END DO
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            !!layer denotes the number of the layer in the snow or in the ice
+            !!jm denotes the reference number of the equation in the tridiagonal
+            !!system, terms of tridiagonal system are indexed as following :
+            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            !!ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
             DO ji = 1, npti
                jk = jm - nlay_s - 1
@@ -367 +441 @@
                zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
             END DO
-         ENDDO
-
-         jm =  nlay_s + nlay_i + 1
-         DO ji = 1, npti
+            ENDDO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
             !!ice bottom term
             ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
@@ -377 +451 @@
             zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
                &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
-         ENDDO
-
-
-         DO ji = 1, npti
+            ENDDO
+
+
+            DO ji = 1, npti
             !                               !---------------------!
-            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
+            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells !
                !                            !---------------------!
                ! snow interior terms (bottom equation has the same form as the others)
                DO jm = 3, nlay_s + 1
-                  jk = jm - 1
-                  ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-                  ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                  ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                  zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+               jk = jm - 1
+               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
+               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
                END DO
 
                ! case of only one layer in the ice (ice equation is altered)
                IF ( nlay_i == 1 ) THEN
-                  ztrid(ji,nlay_s+2,3)  = 0.0
-                  zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+               ztrid(ji,nlay_s+2,3)  = 0.0
+               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
                ENDIF
 
                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
 
-                  jm_min(ji) = 1
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! surface equation
-                  ztrid(ji,1,1)  = 0.0
-                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
-                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
-
-                  ! first layer of snow equation
-                  ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
-                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
+               jm_min(ji) = 1
+               jm_max(ji) = nlay_i + nlay_s + 1
+
+               ! surface equation
+               ztrid(ji,1,1)  = 0.0
+               ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
+               ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
+               zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+
+               ! first layer of snow equation
+               ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
+               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+               ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
+               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
 
                ELSE                            !--  case 2 : surface is melting
-                  !
-                  jm_min(ji) = 2
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! first layer of snow equation
-                  ztrid(ji,2,1)  = 0.0
-                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
-                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+               !
+               jm_min(ji) = 2
+               jm_max(ji) = nlay_i + nlay_s + 1
+
+               ! first layer of snow equation
+               ztrid(ji,2,1)  = 0.0
+               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
+               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
+               &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
             !                               !---------------------!
@@ -433 +507 @@
                !
                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
-                  !
-                  jm_min(ji) = nlay_s + 1
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! surface equation   
-                  ztrid(ji,jm_min(ji),1)  = 0.0
-                  ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
-                  ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
-                  zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
-
-                  ! first layer of ice equation
-                  ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
-                  ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
-                  ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)  
-                  zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
-
-                  ! case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,jm_min(ji),1)    = 0.0
-                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
-                     ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
-                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
-                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,jm_min(ji)+1,3)  = 0.0
-                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
-                        &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
-                  ENDIF
+               !
+               jm_min(ji) = nlay_s + 1
+               jm_max(ji) = nlay_i + nlay_s + 1
+
+               ! surface equation   
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
+               ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
+               zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
+
+               ! first layer of ice equation
+               ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
+               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
+               ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)  
+               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
+
+               ! case of only one layer in the ice (surface & ice equations are altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,jm_min(ji),1)    = 0.0
+               ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
+               ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
+               ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
+               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+               ztrid(ji,jm_min(ji)+1,3)  = 0.0
+               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
+               &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
+               ENDIF
 
                ELSE                            !--  case 2 : surface is melting
 
-                  jm_min(ji)    =  nlay_s + 2
-                  jm_max(ji)    =  nlay_i + nlay_s + 1
-
-                  ! first layer of ice equation
-                  ztrid(ji,jm_min(ji),1)  = 0.0
-                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
-                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
-                     &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
-
-                  ! case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,jm_min(ji),1)  = 0.0
-                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,jm_min(ji),3)  = 0.0
-                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
-                        &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
-                  ENDIF
+               jm_min(ji)    =  nlay_s + 2
+               jm_max(ji)    =  nlay_i + nlay_s + 1
+
+               ! first layer of ice equation
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
+               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
+               &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
+
+               ! case of only one layer in the ice (surface & ice equations are altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+               ztrid(ji,jm_min(ji),3)  = 0.0
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
+               &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
+               ENDIF
 
                ENDIF
@@ -488 +562 @@
             zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
             !
-         END DO
-         !
-         !------------------------------
-         ! 8) tridiagonal system solving
-         !------------------------------
-         ! Solve the tridiagonal system with Gauss elimination method.
-         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-         jm_maxt = 0
-         jm_mint = nlay_i+5
-         DO ji = 1, npti
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
             jm_mint = MIN(jm_min(ji),jm_mint)
             jm_maxt = MAX(jm_max(ji),jm_maxt)
-         END DO
-
-         DO jk = jm_mint+1, jm_maxt
+            END DO
+
+            DO jk = jm_mint+1, jm_maxt
             DO ji = 1, npti
                jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
@@ -508 +582 @@
                zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
             END DO
-         END DO
-
-         DO ji = 1, npti
+            END DO
+
+            DO ji = 1, npti
             ! ice temperatures
             t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-         END DO
-
-         DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
             DO ji = 1, npti
                jk    =  jm - nlay_s - 1
                t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
             END DO
-         END DO
-
-         DO ji = 1, npti
+            END DO
+
+            DO ji = 1, npti
             ! snow temperatures      
             IF( h_s_1d(ji) > 0._wp ) THEN
                t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-                  &                / zdiagbis(ji,nlay_s+1)
+               &                / zdiagbis(ji,nlay_s+1)
             ENDIF
             ! surface temperature
@@ -532 +606 @@
             IF( t_su_1d(ji) < rt0 ) THEN
                t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
-            ENDIF
-         END DO
-         !
-         !--------------------------------------------------------------
-         ! 9) Has the scheme converged ?, end of the iterative procedure
-         !--------------------------------------------------------------
-         ! check that nowhere it has started to melt
-         ! zdti_max is a measure of error, it has to be under zdti_bnd
-         zdti_max = 0._wp
-         DO ji = 1, npti
+               &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+            ENDIF
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged ?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+            DO ji = 1, npti
             t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
             zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
-         END DO
-
-         DO jk  =  1, nlay_s
+            END DO
+
+            DO jk  =  1, nlay_s
             DO ji = 1, npti
                t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
                zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
             END DO
-         END DO
-
-         DO jk  =  1, nlay_i
+            END DO
+
+            DO jk  =  1, nlay_i
             DO ji = 1, npti
                ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
@@ -560 +634 @@
                zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
             END DO
-         END DO
-
-         ! Compute spatial maximum over all errors
-         ! note that this could be optimized substantially by iterating only the non-converging points
-         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
-
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
+         !
+         !----------------------------------------!
+         !                                        !
+         !       JULES IF (RCV MODE)              !
+         !                                        !
+         !----------------------------------------!
+         !
+
+         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode
+         
+            !
+            ! In RCV mode, we use a modified BL99 solver 
+            ! with conduction flux (qcn_ice) as forcing term
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            !!layer denotes the number of the layer in the snow or in the ice
+            !!jm denotes the reference number of the equation in the tridiagonal
+            !!system, terms of tridiagonal system are indexed as following :
+            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            !!ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
+            DO ji = 1, npti
+               jk = jm - nlay_s - 1
+               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
+               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+            END DO
+            ENDDO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
+            !!ice bottom term
+            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
+            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
+            ztrid(ji,jm,3)  = 0.0
+            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
+            ENDDO
+
+
+            DO ji = 1, npti
+            !                              !---------------------!
+            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
+               !                           !---------------------!
+               ! snow interior terms (bottom equation has the same form as the others)
+               DO jm = 3, nlay_s + 1
+               jk = jm - 1
+               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
+               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+               END DO
+
+               ! case of only one layer in the ice (ice equation is altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,nlay_s+2,3)  = 0.0
+               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+               ENDIF
+
+               jm_min(ji) = 2
+               jm_max(ji) = nlay_i + nlay_s + 1
+
+               ! first layer of snow equation
+               ztrid(ji,2,1)  = 0.0
+               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
+               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
+               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
+               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
+
+            !                               !---------------------!
+            ELSE                            ! cells without snow  !
+            !                               !---------------------!
+
+               jm_min(ji)    =  nlay_s + 2
+               jm_max(ji)    =  nlay_i + nlay_s + 1
+
+               ! first layer of ice equation
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)  
+               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
+               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+
+               ! case of only one layer in the ice (surface & ice equations are altered)
+               IF ( nlay_i == 1 ) THEN
+               ztrid(ji,jm_min(ji),1)  = 0.0
+               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
+               ztrid(ji,jm_min(ji),3)  = 0.0
+               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    &
+               &                    + qcn_ice_1d(ji) )
+
+               ENDIF
+
+            ENDIF
+            !
+            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
+            !
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
+            jm_mint = MIN(jm_min(ji),jm_mint)
+            jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+
+            DO jk = jm_mint+1, jm_maxt
+            DO ji = 1, npti
+               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
+               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
+               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
+            END DO
+            END DO
+
+            DO ji = 1, npti
+            ! ice temperatures
+            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+            DO ji = 1, npti
+               jk    =  jm - nlay_s - 1
+               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+            END DO
+            END DO
+
+            DO ji = 1, npti
+            ! snow temperatures      
+            IF( h_s_1d(ji) > 0._wp ) THEN
+               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
+               &                / zdiagbis(ji,nlay_s+1)
+            ENDIF
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged ?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+
+            DO jk  =  1, nlay_s
+            DO ji = 1, npti
+               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+            END DO
+            END DO
+
+            DO jk  =  1, nlay_i
+            DO ji = 1, npti
+               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+            END DO
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
+
+         ENDIF ! k_jules
+         
       END DO  ! End of the do while iterative procedure
-
+      
       IF( ln_icectl .AND. lwp ) THEN
          WRITE(numout,*) ' zdti_max : ', zdti_max
          WRITE(numout,*) ' iconv    : ', iconv
       ENDIF
+      
       !
       !-----------------------------
       ! 10) Fluxes at the interfaces
       !-----------------------------
+      !
+      ! --- update conduction fluxes
+      !
       DO ji = 1, npti
-         !                                ! surface ice conduction flux
+      !                                ! surface ice conduction flux
          fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
             &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
-         !                                ! bottom ice conduction flux
+      !                                ! bottom ice conduction flux
          fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
-
-      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
-      CALL ice_thd_enmelt
-
-      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      IF ( ln_dqns_i ) THEN
+      
+      !
+      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
+      !
+      DO ji = 1, npti
+            IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN 
+                   ! OFF or SND mode
+               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
+            ELSE   ! RCV mode
+               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 
+            ENDIF
+      END DO
+      
+      !
+      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
+      !
+
+      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF
+      
+         CALL ice_thd_enmelt       
+
+         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
          DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
+            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
+               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
+               
+            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN
+               
+                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
+                     zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
+                  ELSE                          ! case T_su = 0degC
+                     zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+                  ENDIF
+            
+            ELSE ! RCV CASE
+            
+               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+            
+            ENDIF
+            
+            ! total heat sink to be sent to the ocean
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+
+            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
+   
+         END DO 
+           
+         !
+         ! --- SIMIP diagnostics
+         !
+
+         DO ji = 1, npti
+            !--- Conduction fluxes (positive downwards)
+            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
+   
+            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
+               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
+                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
+            ELSE
+               t_si_1d(ji) = t_su_1d(ji)
+            ENDIF
          END DO
-      END IF
-
-      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
-      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
-      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-      DO ji = 1, npti
-         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
-            &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
-
-         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
-         ELSE                          ! case T_su = 0degC
-            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
-         ENDIF
-         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
-
-         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
-         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
-
-      END DO   
-
-      ! --- Diagnostics SIMIP --- !
-      DO ji = 1, npti
-         !--- Conduction fluxes (positive downwards)
-         diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-         diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
-
-         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
-            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
-         ELSE
-            t_si_1d(ji) = t_su_1d(ji)
-         ENDIF
-      END DO
-      !
-   END SUBROUTINE ice_thd_zdf
+     
+      ENDIF
+      !
+      !---------------------------------------------------------------------------------------
+      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
+      !---------------------------------------------------------------------------------------
+      !
+      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
+      
+         ! Restore temperatures to their initial values
+         t_s_1d(1:npti,:)        = ztsold (1:npti,:)
+         t_i_1d(1:npti,:)        = ztiold (1:npti,:)
+         qcn_ice_1d(1:npti)      = fc_su(1:npti)
+         
+      ENDIF
+      
+   END SUBROUTINE ice_thd_zdf_BL99
+
 
 
@@ -663 +954 @@
 
 
+
    SUBROUTINE ice_thd_zdf_init
       !!-----------------------------------------------------------------------
@@ -675 +967 @@
       !! ** input   :   Namelist namthd_zdf
       !!-------------------------------------------------------------------
-      INTEGER  ::   ios   ! Local integer output status for namelist read
+      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i
+      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i
       !!-------------------------------------------------------------------
       !
@@ -694 +986 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namthd_zdf:'
-         WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99
+         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99
          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,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
          WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
      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_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
       ENDIF
+
+      !                             !== set the choice of ice vertical thermodynamic formulation ==!
+      ioptio = 0 
+      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties)
+      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF
+      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' )
       !
    END SUBROUTINE ice_thd_zdf_init
@@ -711 +1009 @@
    !!----------------------------------------------------------------------
    !!   Default option       Dummy Module             No ESIM sea-ice model
-   !!----------------------------------------------------------------------
+   !!---------------------------------------------------------------------
 #endif
 

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

@@ -164 +164 @@
             ! mass flux from ice/ocean
             wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
-               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
-
-            IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
+               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
 
             ! add the snow melt water to snow mass flux to the ocean
@@ -199 +197 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
@@ -246 +244 @@
       CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt 
       CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt 
-      IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
+      IF ( ln_pnd_fwb )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
 
       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

@@ -97 +97 @@
       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
 
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN                     ! Melt pond
-         at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
-         vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
-      ENDIF
-      ! END MP 2016
+      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
+      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
 
       ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction  
@@ -167 +163 @@
 !!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
 
-      !-------------------------------------------------------
-      ! Ice thickness, snow thickness, ice salinity, ice age
-      !-------------------------------------------------------
+      !---------------------------------------------------------------
+      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
+      !---------------------------------------------------------------
       !                                            !--- inverse of the ice area
       WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
@@ -178 +174 @@
       ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
       END WHERE
-      !
-      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
+      !                                           !--- ice thickness
+      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
 
       zhmax    =          hi_max(jpl)
       z1_zhmax =  1._wp / hi_max(jpl)               
-      WHERE( h_i(:,:,jpl) > zhmax )               !--- bound h_i by hi_max (i.e. 99 m) with associated update of ice area
+      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
          h_i  (:,:,jpl) = zhmax
          a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
-         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as h_i > hi_max
+         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
       END WHERE
-
-      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)    !--- snow thickness
-      
-      o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age
-
-      IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere)
+      !                                           !--- snow thickness
+      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
+      !                                           !--- ice age      
+      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
+      !                                           !--- pond fraction and thickness      
+      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
+      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
+      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
+      END WHERE
+      !
+      !                                           !---  salinity (with a minimum value imposed everywhere)     
+      IF( nn_icesal == 2 ) THEN
          WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
          ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
          END WHERE
       ENDIF
-
-      CALL ice_var_salprof      ! salinity profile
+      CALL ice_var_salprof   ! salinity profile
 
       !-------------------
@@ -242 +243 @@
       vt_s (:,:) = SUM( v_s, dim=3 )
       at_i (:,:) = SUM( a_i, dim=3 )
-
-! MV MP 2016
-! probably should resum for melt ponds ???
-
       !
    END SUBROUTINE ice_var_glo2eqv
@@ -258 +255 @@
       !!-------------------------------------------------------------------
       !
-      v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:)
-      v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:)
-      sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:)
+      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
+      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
+      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
+      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
       !
    END SUBROUTINE ice_var_eqv2glo
@@ -478 +476 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
                hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
+               IF( ln_pnd_fwb ) THEN 
+                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
+               ENDIF
                !-----------------------------------------------------------------
                ! Zap snow energy 
@@ -498 +499 @@
                h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
 
-            END DO
-         END DO
-
-         IF( ln_pnd ) THEN 
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( ln_pnd_fw )   &
-                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
-                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
-                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         
+               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
+               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+
+            END DO
+         END DO
+        
       END DO 
 
@@ -520 +513 @@
       ! open water = 1 if at_i=0
       WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
+
       !
    END SUBROUTINE ice_var_zapsmall

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

@@ -128 +128 @@
       CALL iom_put( "snowvol"     , vt_s  * zswi      )        ! snow volume
       
-      IF ( ln_pnd ) THEN
-         CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
-         CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
-      ENDIF
+      CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
+      CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
 
       !----------------------------------
@@ -145 +143 @@
       IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume
 
-      IF ( ln_pnd ) THEN
-         IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories
-         IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories
-         IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories
-         IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories
-      ENDIF
+      IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories
+      IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories
+      IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories
+      IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories
 
       !--------------------------------
@@ -316 +312 @@
       CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN
-         CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   &
-      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-         CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   &
-      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      ENDIF
-      ! END MV MP 2016
-
+      CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
@@ -370 +360 @@
       CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
 
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN
-         CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) )
-         CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) )
-      ENDIF
-      ! END MV MP 2016
+      CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) )
 
       CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
@@ -384 +370 @@
       CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
-      IF ( ln_pnd ) &
-         CALL histwrite( kid, "vfxpnd", kt, wfx_pnd     , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxpnd", kt, wfx_pnd        , jpi*jpj, (/1/) )
 
       CALL histwrite( kid, "sithicat", kt, h_i         , jpi*jpj*jpl, (/1/) )    

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90

@@ -19 +19 @@
    USE eosbn2          ! equation of state
    USE oce             ! ocean dynamics and tracers variables
-   USE ice             ! LIM_3 ice variables
-   USE icevar
-   USE icectl
+   USE ice             ! sea-ice: variables
+   USE icevar          ! sea-ice: operations
+   USE iceitd          ! sea-ice: rebining
+   USE icectl          ! sea-ice: control prints
    USE par_oce         ! ocean parameters
    USE dom_oce         ! ocean space and time domain variables 
@@ -71 +72 @@
       END DO
       !
-                        CALL ice_var_zapsmall
-                        CALL ice_var_agg(1)
+      CALL ice_var_zapsmall
+      CALL ice_var_agg(1)
       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
       !
@@ -248 +249 @@
       END DO !jl
       !
+      ! --- In case categories are out of bounds, do a remapping --- !
+      !     i.e. inputs have not the same ice thickness distribution 
+      !          (set by rn_himean) than the regional simulation
+      IF( jpl > 1 )   CALL ice_itd_reb( kt )
       !      
       IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs')

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -48 +48 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice_tr     !: solar flux transmitted below the ice surface      [W/m2]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
@@ -123 +125 @@
       ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
          &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
-         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   &
+         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  &
          &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   &
-         &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
          &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
          &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
@@ -139 +141 @@
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
                 STAT= ierr(2) )
-      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
-         &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , &
+      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , tn_ice (jpi,jpj,1)    , &
+         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , &
          &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , &
          &                     STAT= ierr(3) )      
@@ -168 +170 @@
    REAL            , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-]
    INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

@@ -17 +17 @@
    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
+   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle) 
    !!----------------------------------------------------------------------
 
@@ -42 +43 @@
    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
    USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
@@ -64 +66 @@
    PUBLIC   blk_ice_tau   ! routine called in icestp module
    PUBLIC   blk_ice_flx   ! routine called in icestp module
-#endif
+   PUBLIC   blk_ice_qcn   ! routine called in icestp module
+#endif 
 
 !!Lolo: should ultimately be moved in the module with all physical constants ?
@@ -697 +700 @@
 
 
-   SUBROUTINE blk_ice_flx( ptsu, palb )
+   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_ice_flx  ***
@@ -710 +713 @@
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
       !!
@@ -716 +721 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
+      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice
@@ -722 +729 @@
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3)
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa
+
       !!---------------------------------------------------------------------
       !
@@ -830 +838 @@
       CALL wrk_dealloc( jpi,jpj,   zevap, zsnw )
 
-      !--------------------------------------------------------------------
-      ! FRACTIONs of net shortwave radiation which is not absorbed in the
-      ! thin surface layer and penetrates inside the ice cover
-      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
-      !
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-      !
-      !
+      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
+      !
+      ! former coding was
+      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+
+      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
+
+      qsr_ice_tr(:,:,:) = 0._wp
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
+               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+
+               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN 
+                  zfrqsr_tr_i  = zfr1            + zfac * zfr2
+               ELSE 
+                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
+               ENDIF
+
+               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation 
+
+            END DO
+         END DO
+      END DO
+                                                                                      
+
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
@@ -854 +885 @@
 
    END SUBROUTINE blk_ice_flx
+   
+   
+
+   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi )
+   
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_qcn  ***
+      !!
+      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
+      !!                to force sea ice / snow thermodynamics
+      !!                in the case JULES coupler is emulated
+      !!                
+      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
+      !!                following the 0-layer Semtner (1976) approach
+      !!
+      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
+      !!              - qcn_ice : surface inner conduction flux (W/m2)
+      !!
+      !!---------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option
+
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature
+
+      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness
+      
+      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices
+      INTEGER                                    ::   iter                      !
+      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors
+      REAL(wp)                                   ::   zkeff_h, ztsu             !
+      REAL(wp)                                   ::   zqc, zqnet                !
+      REAL(wp)                                   ::   zhe, zqa0                 !
+           
+      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations
+      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
+
+      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor
+      
+      !!---------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn')
+      !
+      CALL wrk_alloc( jpi,jpj,jpl,   zgfac )
+
+      ! -------------------------------------!
+      !      I   Enhanced conduction factor  !
+      ! -------------------------------------!
+      !
+      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3)
+      ! Fichefet and Morales Maqueda, JGR 1997
+      !
+      zgfac(:,:,:) = 1._wp
+      
+      SELECT CASE ( k_monocat ) 
+      
+         CASE ( 1 , 3 )
+      
+            zfac    =   1._wp /  ( rn_cnd_s + rcdic )
+            zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon
+            zfac3   =   2._wp / zepsilon
+      
+            DO jl = 1, jpl                
+               DO jj = 1 , jpj
+                  DO ji = 1, jpi
+                     !                                ! Effective thickness
+                     zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac
+                     !                                ! Enhanced conduction factor
+                     IF( zhe >=  zfac2 )  &
+                        zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) )
+                  END DO
+               END DO
+            END DO
+      
+      END SELECT
+      
+      ! -------------------------------------------------------------!
+      !      II   Surface temperature and conduction flux            !
+      ! -------------------------------------------------------------!
+
+      zfac   =   rcdic * rn_cnd_s 
+                                      ! ========================== !
+      DO jl = 1, jpl                  !  Loop over ice categories  !
+         !                            ! ========================== !
+         DO jj = 1 , jpj
+            DO ji = 1, jpi
+               !                      ! Effective conductivity of the snow-ice system divided by thickness
+               zkeff_h   =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) )
+                                      ! Store initial surface temperature
+               ztsu      =   ptsu(ji,jj,jl)
+                                      ! Net initial atmospheric heat flux
+               zqa0      =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)
+               
+               DO iter = 1, nit       ! --- Iteration loop
+               
+                   !                  ! Conduction heat flux through snow-ice system (>0 downwards)
+                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) )
+                   !                  ! Surface energy budget
+                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc
+                   !                  ! Temperature update
+                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )
+               
+               END DO
+               
+               ptsu(ji,jj,jl)    = MIN( rt0, ztsu )
+               
+               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
+               
+            END DO
+         END DO
+         
+      END DO 
+      
+      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac )
+      
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn')
+
+   END SUBROUTINE blk_ice_qcn
    
 #endif
@@ -1093 +1243 @@
       zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
       !
-!!      DO jj = 2, jpjm1
-!!         DO ji = fs_2, fs_jpim1
-      DO jj = 1, jpj
-         DO ji = 1, jpi
+      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
+         DO ji = fs_2, fs_jpim1
             ! Virtual potential temperature [K]
             zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
@@ -1140 +1288 @@
          END DO
       END DO
-!!      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
+      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. )
       !
    END SUBROUTINE Cdn10_Lupkes2015

branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -32 +32 @@
    USE geo2ocean      ! 
    USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
-   USE albedooce      ! 
+   USE ocealb         ! 
    USE eosbn2         ! 
    USE sbcrnf, ONLY : l_rnfcpl
@@ -178 +178 @@
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
 
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2] 
@@ -202 +202 @@
       ierr(:) = 0
       !
-      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
+      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
       
 #if ! defined key_lim3 && ! defined key_cice
@@ -737 +737 @@
       !     2. receiving mixed oce-ice solar radiation 
       IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
-         CALL albedo_oce( zaos, zacs )
+         CALL oce_alb( zaos, zacs )
          ! Due to lack of information on nebulosity : mean clear/overcast sky
-         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
+         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
       ENDIF
 
@@ -1530 +1530 @@
    
 
-   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist )
+   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1585 +1585 @@
       REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl         ! dummy loop index
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      !
+      INTEGER ::   ji,jj,jl         ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i
       !!----------------------------------------------------------------------
       !
@@ -1598 +1600 @@
       CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i )
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1890 +1892 @@
 !       ( see OASIS3 user guide, 5th edition, p39 )
          zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( albedo_oce_mix(:,:  ) * ziceld(:,:)       &
+            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
             &                     + palbi         (:,:,1) * picefr(:,:) ) )
       END SELECT
@@ -1951 +1953 @@
       END SELECT
 
-      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
-      ! Used for LIM3
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+      ! --- Transmitted shortwave radiation (W/m2) --- !
+      
+      IF ( nice_jules == 0 ) THEN
+               
+         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
+     
+         ! former coding was    
+         ! Coupled case: since cloud cover is not received from atmosphere 
+         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+         
+         ! to retrieve that coding, we needed to access h_i & h_s from here
+         ! we could even retrieve cloud fraction from the coupler
+               
+         DO jl = 1, jpl
+            DO jj = 1 , jpj
+               DO ji = 1, jpi
+            
+                  !--- surface transmission parameter (Grenfell Maykut 77) --- !
+                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 
+               
+                  ! --- influence of snow and thin ice --- !
+                  IF ( phs(ji,jj,jl) >= 0.0_wp ) zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
+                  IF ( phi(ji,jj,jl) <= 0.1_wp ) zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
+               END DO
+            END DO
+         END DO
+         
+         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation 
+               
+      ENDIF
+      
+      IF ( nice_jules == 2 ) THEN
+      
+         ! here we must receive the qsr_ice_tr array from the coupler
+         ! for now just assume zero
+         
+         qsr_ice_tr(:,:,:) = 0.0_wp
+      
+      ENDIF
+
+
 
       CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
       CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i )
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
@@ -2057 +2096 @@
                    ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
                 ELSEWHERE
-                   ztmp1(:,:) = albedo_oce_mix(:,:)
+                   ztmp1(:,:) = alb_oce_mix(:,:)
                 END WHERE
              CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
@@ -2085 +2124 @@
 
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
-         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
+         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
          DO jl=1,jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)