MODULE ice
   !!======================================================================
   !!                        ***  MODULE ice  ***
   !! LIM-3 Sea Ice physics:  diagnostics variables of ice defined in memory
   !!=====================================================================
   !! History :  3.0  ! 2008-03  (M. Vancoppenolle) original code LIM-3
   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
   !!----------------------------------------------------------------------
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3'                                      LIM-3 sea-ice model
   !!----------------------------------------------------------------------
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! MPP library

   IMPLICIT NONE
   PRIVATE

   PUBLIC    ice_alloc  ! called by icestp.F90

   !!======================================================================
   !! LIM3 by the use of sweat, agile fingers and sometimes brain juice, 
   !!  was developed in Louvain-la-Neuve by : 
   !!    * Martin Vancoppenolle (UCL-ASTR, Belgium)
   !!    * Sylvain Bouillon (UCL-ASTR, Belgium)
   !!    * Miguel Angel Morales Maqueda (NOC-L, UK)
   !! 
   !! Based on extremely valuable earlier work by
   !!    * Thierry Fichefet
   !!    * Hugues Goosse
   !!
   !! The following persons also contributed to the code in various ways
   !!    * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
   !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) 
   !!      and Elisabeth Hunke (LANL), USA.
   !! 
   !! For more info, the interested user is kindly invited to consult the following references
   !!    For model description and validation :
   !!    * Vancoppenolle et al., Ocean Modelling, 2008a.
   !!    * Vancoppenolle et al., Ocean Modelling, 2008b.
   !!    For a specific description of EVP :
   !!    * Bouillon et al., Ocean Modelling 2009.
   !!
   !!    Or the reference manual, that should be available by 2011
   !!======================================================================
   !!                                                                     |
   !!              I C E   S T A T E   V A R I A B L E S                  |
   !!                                                                     |
   !! Introduction :                                                      |
   !! --------------                                                      |
   !! Every ice-covered grid cell is characterized by a series of state   |
   !! variables. To account for unresolved spatial variability in ice     |
   !! thickness, the ice cover in divided in ice thickness categories.    |
   !!                                                                     |
   !! Sea ice state variables depend on the ice thickness category        |
   !!                                                                     |
   !! Those variables are divided into two groups                         |
   !! * Extensive (or global) variables.                                  |
   !!   These are the variables that are transported by all means         |
   !! * Intensive (or equivalent) variables.                              |
   !!   These are the variables that are either physically more           |
   !!   meaningful and/or used in ice thermodynamics                      |
   !!                                                                     |
   !! Routines in icevar.F90 perform conversions                          |
   !!  - ice_var_glo2eqv  : from global to equivalent variables           |
   !!  - ice_var_eqv2glo  : from equivalent to global variables           |
   !!                                                                     |
   !! For various purposes, the sea ice state variables have sometimes    |
   !! to be aggregated over all ice thickness categories. This operation  |
   !! is done in :                                                        |
   !!  - ice_var_agg                                                      |
   !!                                                                     |
   !! in icestp.F90, the routines that compute the changes in the ice     |
   !! state variables are called                                          |
   !! - ice_rhg : ice dynamics                                            |
   !! - ice_adv : ice advection                                           |
   !! - ice_rdgrft : ice ridging and rafting                              |
   !! - ice_thd : ice halo-thermodynamics and creation of new ice         |
   !! - ice_itd : thermodynamic changes in ice thickness distribution     |
   !!                                                                     |
   !! See the associated routines for more information                    |
   !!                                                                     |
   !! List of ice state variables :                                       |
   !! -----------------------------                                       |
   !!                                                                     |
   !!-------------|-------------|---------------------------------|-------|
   !!   name in   |   name in   |              meaning            | units |
   !! 2D routines | 1D routines |                                 |       |
   !!-------------|-------------|---------------------------------|-------|
   !!                                                                     |
   !! ******************************************************************* |
   !! ***         Dynamical variables (prognostic)                    *** |
   !! ******************************************************************* |
   !!                                                                     |
   !! u_ice       |      -      |    Comp. U of the ice velocity  | m/s   |
   !! v_ice       |      -      |    Comp. V of the ice velocity  | m/s   |
   !!                                                                     |
   !! ******************************************************************* |
   !! ***         Category dependent state variables (prognostic)     *** |
   !! ******************************************************************* |
   !!                                                                     |
   !! ** Global variables                                                 |
   !!-------------|-------------|---------------------------------|-------|
   !! a_i         | a_i_1d      |    Ice concentration            |       |
   !! v_i         |      -      |    Ice volume per unit area     | m     |
   !! v_s         |      -      |    Snow volume per unit area    | m     |
   !! smv_i       |      -      |    Sea ice salt content         | ppt.m |
   !! oa_i        !      -      !    Sea ice areal age content    | s     |
   !! e_i         !      -      !    Ice enthalpy                 | J/m2  | 
   !!      -      ! e_i_1d      !    Ice enthalpy per unit vol.   | J/m3  | 
   !! e_s         !      -      !    Snow enthalpy                | J/m2  | 
   !!      -      ! e_s_1d      !    Snow enthalpy per unit vol.  | J/m3  | 
   !!                                                                     |
   !!-------------|-------------|---------------------------------|-------|
   !!                                                                     |
   !! ** Equivalent variables                                             |
   !!-------------|-------------|---------------------------------|-------|
   !!                                                                     |
   !! ht_i        | ht_i_1d     |    Ice thickness                | m     |
   !! ht_s        ! ht_s_1d     |    Snow depth                   | m     |
   !! sm_i        ! sm_i_1d     |    Sea ice bulk salinity        ! ppt   |
   !! s_i         ! s_i_1d      |    Sea ice salinity profile     ! ppt   |
   !! o_i         !      -      |    Sea ice Age                  ! s     |
   !! t_i         ! t_i_1d      |    Sea ice temperature          ! K     |
   !! t_s         ! t_s_1d      |    Snow temperature             ! K     |
   !! t_su        ! t_su_1d     |    Sea ice surface temperature  ! K     |
   !!                                                                     |
   !! notes: the ice model only sees a bulk (i.e., vertically averaged)   |
   !!        salinity, except in thermodynamic computations, for which    |
   !!        the salinity profile is computed as a function of bulk       |
   !!        salinity                                                     |
   !!                                                                     |
   !!        the sea ice surface temperature is not associated to any     |
   !!        heat content. Therefore, it is not a state variable and      |
   !!        does not have to be advected. Nevertheless, it has to be     |
   !!        computed to determine whether the ice is melting or not      |
   !!                                                                     |
   !! ******************************************************************* |
   !! ***         Category-summed state variables (diagnostic)        *** |
   !! ******************************************************************* |
   !! at_i        | at_i_1d     |    Total ice concentration      |       |
   !! vt_i        |      -      |    Total ice vol. per unit area | m     |
   !! vt_s        |      -      |    Total snow vol. per unit ar. | m     |
   !! smt_i       |      -      |    Mean sea ice salinity        | ppt   |
   !! tm_i        |      -      |    Mean sea ice temperature     | K     |
   !! et_i        !      -      !    Total ice enthalpy           | J/m2  | 
   !! et_s        !      -      !    Total snow enthalpy          | J/m2  | 
   !! bv_i        !      -      !    relative brine volume        | ???   | 
   !!=====================================================================

   !!--------------------------------------------------------------------------
   !! * Share Module variables
   !!--------------------------------------------------------------------------
   !                                     !!** ice-generic parameters namelist (namice_run) **
   INTEGER           , PUBLIC ::   jpl             !: number of ice  categories 
   INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers 
   INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers 
   LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F)
   LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F)
   REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere
   REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere
   CHARACTER(len=32) , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
   CHARACTER(len=32) , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory

   !                                     !!** ice-itd namelist (namice_itd) **
   REAL(wp), PUBLIC ::   rn_himin         !: minimum ice thickness
   
   !                                     !!** ice-dynamics namelist (namice_dyn) **
   LOGICAL,  PUBLIC ::   ln_dynFULL       !:  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
   LOGICAL,  PUBLIC ::   ln_dynRHGADV     !:  dyn.: no ridge/raft & no corrections  (rheology + advection)
   LOGICAL,  PUBLIC ::   ln_dynADV        !:  dyn.: only advection w prescribed vel.(rn_uvice + advection)
   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice
   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice) 
   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
   !
   !                                     !!** ice-rdige/raft namelist (namice_rdgrft) **
   LOGICAL , PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79)
   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79
   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
   LOGICAL , PUBLIC ::   ln_str_R75       !: ice strength parameterization (Rothrock75)
   REAL(wp), PUBLIC ::   rn_perdg         !: ridging work divided by pot. energy change in ridging
   !
   !                                     !!** ice-rheology namelist (namice_rhg) **
   LOGICAL , PUBLIC ::   ln_rhg_EVP       !: EVP rheology
   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
   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
   !
   !                                     !!** ice-advection namelist (namice_adv) **
   LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme
   LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme
   INTEGER , PUBLIC ::   nn_UMx           !: order of the UMx advection scheme   
   !
   !                                     !!** ice-thermodynamics namelist (namice_thd) **
                                          ! -- icethd_dif -- !
   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
   LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964)
   LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007)
   LOGICAL , PUBLIC ::   ln_dqns_i        !: change non-solar surface flux with changing surface temperature (T) or not (F)
   INTEGER , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1-4) or not (0)
   REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]
                                          ! -- icethd_dh -- !
   LOGICAL , PUBLIC ::   ln_icedH         !: activate ice thickness change from growing/melting (T) or not (F)
   REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice
                                          ! -- icethd_da -- !
   LOGICAL , PUBLIC ::   ln_icedA         !: activate lateral melting param. (T) or not (F)
   REAL(wp), PUBLIC ::   rn_beta          !: coef. beta for lateral melting param.
   REAL(wp), PUBLIC ::   rn_dmin          !: minimum floe diameter for lateral melting param.
                                          ! -- icethd_lac -- !
   LOGICAL , PUBLIC ::   ln_icedO         !: activate ice growth in open-water (T) or not (F)
   REAL(wp), PUBLIC ::   rn_hinew         !: thickness for new ice formation (m)
   LOGICAL , PUBLIC ::   ln_frazil        !: use of frazil ice collection as function of wind (T) or not (F)
   REAL(wp), PUBLIC ::   rn_maxfraz       !: maximum portion of frazil ice collecting at the ice bottom
   REAL(wp), PUBLIC ::   rn_vfraz         !: threshold drift speed for collection of bottom frazil ice
   REAL(wp), PUBLIC ::   rn_Cfraz         !: squeezing coefficient for collection of bottom frazil ice
                                          ! -- icethd -- !
   INTEGER , PUBLIC ::   nn_iceflx        !: Redistribute heat flux over ice categories
   !                                      !   =-1  Do nothing (needs N(cat) fluxes)
   !                                      !   = 0  Average N(cat) fluxes then apply the average over the N(cat) ice 
   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice
   !                                      !                                   using T-ice and albedo sensitivity
   !                                      !   = 2  Redistribute a single flux over categories

   !                                     !!** ice-salinity namelist (namice_sal) **
   LOGICAL , PUBLIC ::   ln_icedS         !: activate gravity drainage and flushing (T) or not (F)
   INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model
   !                                      ! 1 - constant salinity in both space and time
   !                                      ! 2 - prognostic salinity (s(z,t))
   !                                      ! 3 - salinity profile, constant in time
   REAL(wp), PUBLIC ::   rn_icesal        !: bulk salinity (ppt) in case of constant salinity
   REAL(wp), PUBLIC ::   rn_sal_gd        !: restoring salinity for gravity drainage [PSU]
   REAL(wp), PUBLIC ::   rn_time_gd       !: restoring time constant for gravity drainage (= 20 days) [s]
   REAL(wp), PUBLIC ::   rn_sal_fl        !: restoring salinity for flushing [PSU]
   REAL(wp), PUBLIC ::   rn_time_fl       !: restoring time constant for gravity drainage (= 10 days) [s]
   REAL(wp), PUBLIC ::   rn_simax         !: maximum ice salinity [PSU]
   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]

   ! MV MP 2016
   !                                     !!** melt pond namelist (namicemp)
   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
   !                                     !!** ice-diagnostics namelist (namice_dia) **
   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
   LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F)
   LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F)
   INTEGER , PUBLIC ::   iiceprt          !: debug i-point
   INTEGER , PUBLIC ::   jiceprt          !: debug j-point

   !                                     !!** some other parameters 
   INTEGER , PUBLIC ::   kt_ice           !: iteration number
   REAL(wp), PUBLIC ::   rdt_ice          !: ice time step
   REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice
   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i
   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s 
   REAL(wp), PUBLIC ::   rswitch          !: switch for the presence of ice (1) or not (0)
   REAL(wp), PUBLIC ::   rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft   !: conservation diagnostics
   REAL(wp), PUBLIC, PARAMETER ::   epsi06 = 1.e-06_wp  !: small number 
   REAL(wp), PUBLIC, PARAMETER ::   epsi10 = 1.e-10_wp  !: small number 
   REAL(wp), PUBLIC, PARAMETER ::   epsi20 = 1.e-20_wp  !: small number 


   !                                     !!** define arrays
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce,v_oce !: surface ocean velocity used in ice dynamics
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hicol       !: ice collection thickness accreted in leads
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength    !: ice strength
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   delta_i     !: ice rheology elta factor (Flato & Hibler 95) [s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   divu_i      !: Divergence of the velocity field [s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field [s-1]
   !
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qlead       !: heat balance of the lead (or of the open ocean)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhtur       !: net downward heat flux from the ice to the ocean
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange   [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sni !: snow ice growth component of wfx_snw [kg.m-2.s-1]
   ! MV MP 2016
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1]
   ! END MV MP 2016
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice  [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice    [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: snow sublimation           [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: ice sublimation            [kg.m-2.s-1]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw    [kg.m-2.s-1]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: ice-ocean mass exchange                   [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni     !: snow ice growth component of wfx_ice      [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw     !: lateral ice growth component of wfx_ice   [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog     !: bottom ice growth component of wfx_ice    [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn     !: dynamical ice growth component of wfx_ice [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom     !: bottom melt component of wfx_ice          [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum     !: surface melt component of wfx_ice         [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_lam     !: lateral melt component of wfx_ice         [kg.m-2.s-1]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res     !: residual component of wfx_ice             [kg.m-2.s-1]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_lam     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sum     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sni     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_opw     !: salt flux due to ice growth/melt                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bri     !: salt flux due to brine rejection                      [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_dyn     !: salt flux due to porous ridged ice formation          [PSU/m2/s]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth        [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt          [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt         [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion             [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations   [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation [kg.m-2.s-1]
   
   ! heat flux associated with ice-atmosphere mass exchange
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2]

   ! heat flux associated with ice-ocean mass exchange
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (icethd_dh) [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from ridging                      [W.m-2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d     !: maximum ice concentration 2d array
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice        !: transmitted solar radiation under ice

   !!--------------------------------------------------------------------------
   !! * Ice global state variables
   !!--------------------------------------------------------------------------
   !! Variables defined for each ice category
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ht_i      !: Ice thickness (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i       !: Ice fractional areas (concentration)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_i       !: Ice volume per unit area (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_s       !: Snow volume per unit area(m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ht_s      !: Snow thickness (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_su      !: Sea-Ice Surface Temperature (K)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sm_i      !: Sea-Ice Bulk salinity (ppt)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   smv_i     !: Sea-Ice Bulk salinity times volume per area (ppt.m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (s)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (s)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i      !: brine volume

   !! Variables summed over all categories, or associated to all the ice in a single grid cell
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity (m/s)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area (m)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i        !: =1-at_i ; total open water fractional area
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s  !: ice and snow total heat content
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i        !: mean sea ice salinity averaged over all categories [PSU]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_i        !: mean ice  thickness over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tau_icebfr   !: ice friction with bathy (landfast param activated)

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s      !: Snow temperatures [K]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s      !: Snow ...      
      
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i      !: ice temperatures          [K]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i      !: ice thermal contents    [J/m2]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   s_i      !: ice salinities          [PSU]

   ! MV MP 2016
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond fraction per grid cell area
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area [m]
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond volume per ice area
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond thickness [m]

   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area [m]
   ! END MV MP 2016

   !!--------------------------------------------------------------------------
   !! * Moments for advection
   !!--------------------------------------------------------------------------
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   !: open water in sea ice
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   !: ice thickness 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    !: snow thickness
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     !: lead fraction
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    !: snow thermal content
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   !: ice salinity
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   !: ice age
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     !: ice layers heat content
   ! MV MP 2016
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    !:  melt pond fraction
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    !:  melt pond volume
   ! END MV MP 2016

   !!--------------------------------------------------------------------------
   !! * Old values of global variables
   !!--------------------------------------------------------------------------
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, ht_s_b, ht_i_b  !: snow and ice volumes/thickness
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, smv_i_b, oa_i_b        !:
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                         !: snow heat content
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                         !: ice temperatures
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b              !: ice velocity
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                        !: ice concentration (total)
            
   !!--------------------------------------------------------------------------
   !! * Ice thickness distribution variables
   !!--------------------------------------------------------------------------
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max         !: Boundary of ice thickness categories in thickness space
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean        !: Mean ice thickness in catgories 
   !
   !!--------------------------------------------------------------------------
   !! * Ice diagnostics
   !!--------------------------------------------------------------------------
   ! thd refers to changes induced by thermodynamics
   ! trp   ''         ''     ''       advection (transport of ice)
   !
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vi   !: transport of ice volume
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vs   !: transport of snw volume
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_smv  !: transport of salt content
   !
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat     !: snw/ice heat content variation   [W/m2] 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_smvi     !: ice salt content variation   [] 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s] 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s] 

   !
   !!--------------------------------------------------------------------------
   !! * SIMIP extra diagnostics
   !!--------------------------------------------------------------------------
   ! Extra sea ice diagnostics to address the data request
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_si          !: Temperature at Snow-ice interface (K) 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K) 
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)

   !
   !!----------------------------------------------------------------------
   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   FUNCTION ice_alloc()
      !!-----------------------------------------------------------------
      !!               *** Routine ice_alloc ***
      !!-----------------------------------------------------------------
      INTEGER :: ice_alloc
      !
      INTEGER :: ierr(18), ii
      !!-----------------------------------------------------------------

      ierr(:) = 0

      ! What could be one huge allocate statement is broken-up to try to
      ! stay within Fortran's max-line length limit.
      ii = 1
      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , hicol    (jpi,jpj) ,                        &
         &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  &
         &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) )

      ii = ii + 1
      ALLOCATE( t_bo   (jpi,jpj) , wfx_snw_sni(jpi,jpj) ,                                                &
         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
         &      wfx_ice(jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
         ! MV MP 2016
         &      wfx_pnd(jpi,jpj) ,                                                              &
         ! END MV MP 2016
         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                         &
         &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                           &
         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  &
         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld   (jpi,jpj) ,                        &
         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     &
         &      hfx_opw(jpi,jpj) , hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,     &
         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj)        ,  STAT=ierr(ii) )

      ! * Ice global state variables
      ii = ii + 1
      ALLOCATE( ftr_ice(jpi,jpj,jpl) ,                                                 &
         &      ht_i   (jpi,jpj,jpl) , a_i   (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,     &
         &      v_s    (jpi,jpj,jpl) , ht_s  (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,     &
         &      sm_i   (jpi,jpj,jpl) , smv_i (jpi,jpj,jpl) , o_i   (jpi,jpj,jpl) ,     &
         &      oa_i   (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                       &
         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     &
         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                              , STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , s_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )

      ! MV MP 2016
      ii = ii + 1
      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , &
         &      h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
      ii = ii + 1
      ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
      ! END MV MP 2016

      ! * Moments for advection
      ii = ii + 1
      ALLOCATE( sxopw(jpi,jpj) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj) , STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
         &      STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
         &      sxc0 (jpi,jpj,jpl) , syc0 (jpi,jpj,jpl) , sxxc0 (jpi,jpj,jpl) , syyc0 (jpi,jpj,jpl) , sxyc0 (jpi,jpj,jpl) ,   &
         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   &
         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   &
         &      STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) ,     &
         &      syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl)                            , STAT=ierr(ii) )

      ! MV MP 2016
      ii = ii + 1
      ALLOCATE( sxap(jpi,jpj,jpl) , syap(jpi,jpj,jpl) , sxxap(jpi,jpj,jpl) , syyap(jpi,jpj,jpl) , sxyap(jpi,jpj,jpl) ,   &
         &      sxvp(jpi,jpj,jpl) , syvp(jpi,jpj,jpl) , sxxvp(jpi,jpj,jpl) , syyvp(jpi,jpj,jpl) , sxyvp(jpi,jpj,jpl) ,   &
         &      STAT = ierr(ii) )
      ! END MV MP 2016

      ! * Old values of global variables
      ii = ii + 1
      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b  (jpi,jpj,jpl) , ht_s_b(jpi,jpj,jpl)        , ht_i_b(jpi,jpj,jpl)        ,   &
         &      a_i_b  (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b (jpi,jpj,nlay_i,jpl) , e_s_b (jpi,jpj,nlay_s,jpl) ,   &
         &      oa_i_b (jpi,jpj,jpl)                                                     , STAT=ierr(ii) )
      ii = ii + 1
      ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
      
      ! * Ice thickness distribution variables
      ii = ii + 1
      ALLOCATE( hi_max(0:jpl), hi_mean(jpl),  STAT=ierr(ii) )

      ! * Ice diagnostics
      ii = ii + 1
      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),   & 
         &      diag_trp_es(jpi,jpj) , diag_trp_smv(jpi,jpj) , diag_heat  (jpi,jpj),   &
         &      diag_smvi  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) )

      ! * SIMIP diagnostics
      ii = ii + 1
      ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
                diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
                STAT = ierr(ii) )

      ice_alloc = MAXVAL( ierr(:) )
      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc: failed to allocate arrays.')
      !
   END FUNCTION ice_alloc

#else
   !!----------------------------------------------------------------------
   !!   Default option         Empty module            NO LIM sea-ice model
   !!----------------------------------------------------------------------
#endif

   !!======================================================================
END MODULE ice