Changeset 834

Timestamp:

2008-03-07T18:11:35+01:00 (16 years ago)

Author:

ctlod

Message:

Clean comments and useless lines, see ticket:#72

Location:

trunk/NEMO/LIM_SRC_3

Files:

• dom_ice.F90 (modified) (1 diff)
• ice.F90 (modified) (11 diffs)
• iceini.F90 (modified) (12 diffs)
• icestp.F90 (modified) (3 diffs)
• limadv.F90 (modified) (1 diff)
• limcons.F90 (modified) (4 diffs)
• limdia.F90 (modified) (14 diffs)
• limdyn.F90 (modified) (12 diffs)
• limflx.F90 (modified) (4 diffs)
• limhdf.F90 (modified) (1 diff)
• limicepoints.F90 (deleted)
• limistate.F90 (modified) (3 diffs)
• limitd_me.F90 (modified) (67 diffs)
• limitd_th.F90 (modified) (33 diffs)
• limmsh.F90 (modified) (6 diffs)
• limrhg.F90 (modified) (34 diffs)
• limrst.F90 (modified) (3 diffs)
• limtab.F90 (modified) (2 diffs)
• limthd.F90 (modified) (12 diffs)
• limthd_dh.F90 (modified) (43 diffs)
• limthd_dif.F90 (modified) (16 diffs)
• limthd_ent.F90 (modified) (35 diffs)
• limthd_lab.F90 (deleted)
• limthd_lac.F90 (modified) (17 diffs)
• limthd_sal.F90 (modified) (12 diffs)
• limtrp.F90 (modified) (7 diffs)
• limupdate.F90 (modified) (8 diffs)
• limvar.F90 (modified) (11 diffs)
• limwri.F90 (modified) (5 diffs)
• thd_ice.F90 (modified) (3 diffs)

trunk/NEMO/LIM_SRC_3/dom_ice.F90

@@ -6 +6 @@
    !! History :
    !!   2.0  !  03-08  (C. Ethe)  Free form and module
+   !!   3.0  !  03-08  (M. Vancoppenolle), Ice cats, salinity and EVP-C
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !!   LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005)
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/dom_ice.F90,v 1.4 2005/03/27 18:34:40 opalod Exp $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt

trunk/NEMO/LIM_SRC_3/ice.F90

@@ -1 +1 @@
 MODULE ice
-   !!======================================================================
-   !!                        ***  MODULE ice  ***
-   !! Sea Ice physics:  diagnostics variables of ice defined in memory
-   !!=====================================================================
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3' :                                   LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !! History :
    !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !!   3.0  !  08-03  (M. Vancoppenolle) : LIM3 !
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !!  LIM 3.0, UCL-LOCEAN-IPSL (2005)
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/ice.F90,v 1.4 2005/03/27 18:34:41 opalod Exp $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
@@ -20 +17 @@
    IMPLICIT NONE
    PRIVATE
+   !!
+   !!======================================================================
+   !!                        ***  MODULE ice  ***
+   !!
+   !!                           **************
+   !!                           * L I M  3.0 *
+   !!                           **************
+   !!
+   !!                         ''in ice we trust''
+   !!
+   !!                   This module contains the sea ice 
+   !!                 diagnostics variables of ice defined 
+   !!                             in memory
+   !!
+   !!======================================================================
+   !!
+   !! 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 (POL, 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 
+   !!      and Rachid Benshila (LOCEAN-IPSL, France)
+   !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
+   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) 
+   !!      and Elisabeth Hunke (LANL), USA.
+   !! 
+   !! (c) UCL-ASTR, 2005-2008
+   !!
+   !! 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., in prep for 2008.
+   !!
+   !!    Or the reference manual, that should be available by 2009
+   !!
+   !!======================================================================
+   !!                                                                     |
+   !!            *****************************************                |
+   !!            *                                       *                |
+   !! ************ 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 limvar.F90 perform conversions                          |
+   !!  - lim_var_glo2eqv  : from global to equivalent variables           |
+   !!  - lim_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 :                                                        |
+   !!  - lim_var_agg                                                      |
+   !!                                                                     |
+   !! in icestp.F90, the routines that compute the changes in the ice     |
+   !! state variables are called                                          |
+   !! - lim_dyn : ice dynamics                                            |
+   !! - lim_trp : ice transport                                           |
+   !! - lim_itd_me : mechanical redistribution (ridging and rafting)      |
+   !! - lim_thd : ice halo-thermodynamics                                 |
+   !! - lim_itd_th : thermodynamic changes in ice thickness distribution  |
+   !!                and creation of new ice                              |
+   !!                                                                     |
+   !! 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_b       |    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    | day   |
+   !! e_i         !      -      !    Ice enthalpy                 | 10^9 J| 
+   !!      -      ! q_i_b       !    Ice enthalpy per unit vol.   | J/m3  | 
+   !! e_s         !      -      !    Snow enthalpy                | 10^9 J| 
+   !!      -      ! q_s_b       !    Snow enthalpy per unit vol.  | J/m3  | 
+   !!                                                                     |
+   !!-------------|-------------|---------------------------------|-------|
+   !!                                                                     |
+   !! ** Equivalent variables                                             |
+   !!                                                                     |
+   !!-------------|-------------|---------------------------------|-------|
+   !!                                                                     |
+   !! ht_i        | ht_i_b      |    Ice thickness                | m     |
+   !! ht_s        ! ht_s_b      |    Snow depth                   | m     |
+   !! sm_i        ! sm_i_b      |    Sea ice bulk salinity        ! ppt   |
+   !! s_i         ! s_i_b       |    Sea ice salinity profile     ! ppt   |
+   !! o_i         !      -      |    Sea ice Age                  ! days  |
+   !! t_i         ! t_i_b       |    Sea ice temperature          ! K     |
+   !! t_s         ! t_s_b       |    Snow temperature             ! K     |
+   !! t_su        ! t_su_b      |    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_b      |    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     |
+   !! ot_i        !      -      !    Sea ice areal age content    | day   |
+   !! et_i        !      -      !    Total ice enthalpy           | 10^9 J| 
+   !! et_s        !      -      !    Total snow enthalpy          | 10^9 J| 
+   !! bv_i        !      -      !    Mean relative brine volume   | ???   | 
+   !!                                                                     |
+   !!                                                                     |
+   !!=====================================================================
 
    LOGICAL, PUBLIC ::    &
       con_i = .false.           ! switch for conservation test
 
+   !!--------------------------------------------------------------------------
    !! * Share Module variables
-   INTEGER , PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) **
-      nbiter = 1      ,  &  !: number of sub-time steps for relaxation
-      nbitdr = 250    ,  &  !: maximum number of iterations for relaxation
-      nevp   = 400    ,  &  !: number of iterations for subcycling
-      nlay_i = 5            !: number of layers in the ice
-
-   REAL(wp), PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) **
-      epsd   = 1.0e-20,  &  !: tolerance parameter for dynamic
-      alpha  = 0.5    ,  &  !: coefficient for semi-implicit coriolis
-      dm     = 0.6e+03,  &  !: diffusion constant for dynamics
-      om     = 0.5    ,  &  !: relaxation constant
-      resl   = 5.0e-05,  &  !: maximum value for the residual of relaxation
-      cw     = 5.0e-03,  &  !: drag coefficient for oceanic stress
-      angvg  = 0.0    ,  &  !: turning angle for oceanic stress
-      pstar  = 1.0e+04,  &  !: first bulk-rheology parameter
-      c_rhg  = 20.0   ,  &  !: second bulk-rhelogy parameter
-      etamn  = 0.0e+07,  &  !: minimun value for viscosity
-      creepl = 2.0e-08,  &  !: creep limit
+   !!--------------------------------------------------------------------------
+   INTEGER , PUBLIC ::   &     !!: ** ice-dynamic namelist (namicedyn) **
+      nbiter = 1      ,  &     !: number of sub-time steps for relaxation
+      nbitdr = 250    ,  &     !: maximum number of iterations for relaxation
+      nevp   = 400    ,  &     !: number of iterations for subcycling
+      nlay_i = 5               !: number of layers in the ice
+
+   REAL(wp), PUBLIC ::   &     !!: ** ice-dynamic namelist (namicedyn) **
+      epsd   = 1.0e-20,  &     !: tolerance parameter for dynamic
+      alpha  = 0.5    ,  &     !: coefficient for semi-implicit coriolis
+      dm     = 0.6e+03,  &     !: diffusion constant for dynamics
+      om     = 0.5    ,  &     !: relaxation constant
+      resl   = 5.0e-05,  &     !: maximum value for the residual of relaxation
+      cw     = 5.0e-03,  &     !: drag coefficient for oceanic stress
+      angvg  = 0.0    ,  &     !: turning angle for oceanic stress
+      pstar  = 1.0e+04,  &     !: determines ice strength (N/M), Hibler JPO79
+      c_rhg  = 20.0   ,  &     !: determines changes in ice strength
+      etamn  = 0.0e+07,  &     !: minimun value for viscosity : has to be 0
+      creepl = 2.0e-08,  &     !: creep limit : has to be under 1.0e-9
       ecc    = 2.0    ,  &  !: eccentricity of the elliptical yield curve
       ahi0   = 350.e0 ,  &  !: sea-ice hor. eddy diffusivity coeff. (m2/s)
-      ! EVP1
       telast = 2880.0 ,  &  !: timescale for elastic waves (s) !SB
       alphaevp = 1.0  ,  &  !: coeficient of the internal stresses !SB
@@ -51 +214 @@
 
    REAL(wp), PUBLIC ::   & !!: ** ice-salinity namelist (namicesal) **
-      s_i_max  = 20.0 ,  &  !: maximum ice salinity
-      s_i_min  =  0.1 ,  &  !: minimum ice salinity
+      s_i_max  = 20.0 ,  &  !: maximum ice salinity (ppt)
+      s_i_min  =  0.1 ,  &  !: minimum ice salinity (ppt)
       s_i_0    =  3.5 ,  &  !: 1st sal. value for the computation of sal .prof.
+                            !: (ppt)
       s_i_1    =  4.5 ,  &  !: 2nd sal. value for the computation of sal .prof.
+                            !: (ppt)
       sal_G    = 5.00 ,  &  !: restoring salinity for gravity drainage
+                            !: (ppt)
       sal_F    = 2.50 ,  &  !: restoring salinity for flushing
-      time_G   = 1.728e+06,&!: restoring time constant for gravity drainage 
-      time_F   = 1.728e+06,&!: restoring time constant for gravity drainage 
-      bulk_sal = 4.0       !!: bulk salinity in case of constant salinity
+                            !: (ppt)
+      time_G   = 1.728e+06,&!: restoring time constant for gravity drainage
+                            !: (= 20 days, in s)
+      time_F   = 8.640e+05,&!: restoring time constant for gravity drainage 
+                            !: (= 10 days, in s)
+      bulk_sal = 4.0        !: bulk salinity (ppt) in case of constant salinity
 
    INTEGER , PUBLIC ::   & !!: ** ice-salinity namelist (namicesal) **
-      num_sal  = 1    ,  &  !:
+      num_sal  = 1    ,  &  !: salinity configuration used in the model
+                            !: 1 - s constant in space and time
+                            !: 2 - prognostic salinity (s(z,t))
+                            !: 3 - salinity profile, constant in time
+                            !: 4 - salinity variations affect only ice
+                            !      thermodynamics
       sal_prof = 1    ,  &  !: salinity profile or not 
-      thcon_i_swi = 1       !: thermal conductivity of Untersteiner (1964) or
-                            !: Pringle et al (2007) 
+      thcon_i_swi = 1       !: thermal conductivity of Untersteiner (1964) (1) or
+                            !: Pringle et al (2007) (2)
 
    REAL(wp), PUBLIC ::   & !!: ** ice-mechanical redistribution namelist (namiceitdme)
@@ -74 +248 @@
       Gstar = 0.15    ,  & !!: fractional area of young ice contributing to ridging
       astar = 0.05    ,  & !!: equivalent of G* for an exponential participation function
-      Hstar = 100.0   ,  & !!: fractional area of young ice contributing to ridging
-      hparmeter = 0.75,  & !!: threshold thickness for rafting / ridging 
-      Craft = 5.0     ,  & !!: coefficient for the hyperbolic tangent in rafting
-      ridge_por = 0.0 ,  & !!: initial porosity of ridges
-      sal_max_ridge = 15.0, & !!: maximum ridged ice salinity
+      Hstar = 100.0   ,  & !!: thickness that determines the maximal thickness of ridged
+                           !!: ice
+      hparmeter = 0.75,  & !!: threshold thickness (m) for rafting / ridging 
+      Craft = 5.0     ,  & !!: coefficient for smoothness of the hyperbolic tangent in rafting
+      ridge_por = 0.0 ,  & !!: initial porosity of ridges (0.3 regular value)
+      sal_max_ridge = 15.0, & !!: maximum ridged ice salinity (ppt)
       betas    = 1.0      , & !:: coef. for partitioning of snowfall between leads and sea ice
       kappa_i  = 1.0      , & !!: coefficient for the extinction of radiation
-      nconv_i_thd = 50    , & !!: maximal number of iterations in heat diffusion
-      maxer_i_thd = 1.0e-4    !!: maximal tolerated error for heat diffusion
+                              !!: Grenfell et al. (2006) (m-1)
+      nconv_i_thd = 50    , & !!: maximal number of iterations for heat diffusion
+      maxer_i_thd = 1.0e-4    !!: maximal tolerated error (C) for heat diffusion
 
    INTEGER , PUBLIC ::   & !!: ** ice-mechanical redistribution namelist (namiceitdme)
       ridge_scheme_swi = 0, & !!: scheme used for ice ridging
       raftswi          = 1, & !!: rafting of ice or not                        
-      partfun_swi      = 1, & !!: participation function TH75 (0) or Letal07 (1) 
-      transfun_swi     = 0, & !!: transfer function of H80 (0) or Letal07 (1)
+      partfun_swi      = 1, & !!: participation function Thorndike et al. JGR75 (0) 
+                              !!: or Lipscomb et al. JGR07 (1) 
+      transfun_swi     = 0, & !!: transfer function of Hibler, MWR80 (0) 
+                              !!: or Lipscomb et al., 2007 (1)
       brinstren_swi    = 0    !!: use brine volume to diminish ice strength
 
@@ -162 +340 @@
 
    INTEGER, PUBLIC, DIMENSION(jpi, jpj, jpl) ::          &   !::
-      patho_case ! number of the pathological case (if any)
-
-!---------------------------------------------------------------------------------------------------
-!  ) Ice global state variables
-!---------------------------------------------------------------------------------------------------
+      patho_case ! number of the pathological case (if any, of course)
+
+   !!--------------------------------------------------------------------------
+   !! * Ice global state variables
+   !!--------------------------------------------------------------------------
+   !! Variables defined for each ice category
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl)        ::   &  !:
-      ht_i   ,   &  !: Ice thickness
-      a_i    ,   &  !: Ice fractional areas
-      v_i    ,   &  !: Ice volume per unit area
-      v_s    ,   &  !: Snow volume per unit area
-      ht_s   ,   &  !: Snow thickness
-      t_su   ,   &  !: Sea-Ice Surface Temperature [K]
-      sm_i   ,   &  !: Sea-Ice Bulk salinity [ppt]
-      smv_i  ,   &  !: Sea-Ice Bulk salinity [ppt] times volume
+      ht_i   ,   &  !: Ice thickness (m)
+      a_i    ,   &  !: Ice fractional areas (concentration)
+      v_i    ,   &  !: Ice volume per unit area (m)
+      v_s    ,   &  !: Snow volume per unit area(m)
+      ht_s   ,   &  !: Snow thickness (m)
+      t_su   ,   &  !: Sea-Ice Surface Temperature (K)
+      sm_i   ,   &  !: Sea-Ice Bulk salinity (ppt)
+      smv_i  ,   &  !: Sea-Ice Bulk salinity times volume per area (ppt.m)
+                    !: this is an extensive variable that has to be transported
       o_i    ,   &  !: Sea-Ice Age (days)
-      ov_i   ,   &  !: Sea-Ice Age times volume / area (days.m)
+      ov_i   ,   &  !: Sea-Ice Age times volume per area (days.m)
       oa_i          !: Sea-Ice Age times ice area (days)
 
+   !! Variables summed over all categories, or associated to 
+   !! all the ice in a single grid cell
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
       u_ice, v_ice,   &  !: two components of the ice velocity (m/s)
       tio_u, tio_v,   &  !: two components of the ice-ocean stress (N/m2)
-      vt_i        ,   &  !: ice total volume
-      vt_s        ,   &  !: snow total volume
-      at_i        ,   &  !: ice total fractional area
-      ato_i       ,   &  !: total open water fractional area
+      vt_i        ,   &  !: ice total volume per unit area (m)
+      vt_s        ,   &  !: snow total volume per unit area (m)
+      at_i        ,   &  !: ice total fractional area (ice concentration)
+      ato_i       ,   &  !: total open water fractional area (1-at_i)
+      et_i        ,   &  !: total ice heat content
+      et_s        ,   &  !: total snow heat content
       ot_i        ,   &  !: mean age over all categories
-      et_i        ,   &  !: total ice heat content
-      bv_i        ,   &  !: total ice heat content
-      tm_i        ,   &  !: mean ice temperature
-      et_s        ,   &  !: total snow heat content
-      smt_i
+      tm_i        ,   &  !: mean ice temperature over all categories
+      bv_i        ,   &  !: brine volume averaged over all categories
+      smt_i              !: mean sea ice salinity averaged over all categories
 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpm) ::   &  !:
@@ -216 +398 @@
       reslum        !: Relative absorption of solar radiation in each ocean level
 
-!---------------------------------------------------------------------------------------------------
-!  ) Moments for Advection
-!---------------------------------------------------------------------------------------------------
-
+   !!--------------------------------------------------------------------------
+   !! * Moments for advection
+   !!--------------------------------------------------------------------------
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  !:
          sxopw, syopw, sxxopw, syyopw, sxyopw          !: open water in sea ice
@@ -234 +415 @@
          sxe ,  sye ,  sxxe ,  syye ,  sxye            !: ice layers heat content
 
-!---------------------------------------------------------------------------------------------------
-!  ) Old values of global variables
-!---------------------------------------------------------------------------------------------------
+   !!--------------------------------------------------------------------------
+   !! * Old values of global variables
+   !!--------------------------------------------------------------------------
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: 
       old_v_s, old_v_i,                          &  !: snow and ice volumes
@@ -247 +428 @@
       old_u_ice, old_v_ice
 
-!---------------------------------------------------------------------------------------------------
-!  ) Values of increments of global variables
-!---------------------------------------------------------------------------------------------------
-! thd refers to changes induced by thermodynamics
-! trp   ''         ''     ''       advection (transport of ice)
-
+   !!--------------------------------------------------------------------------
+   !! * Increment of global variables
+   !!--------------------------------------------------------------------------
+   ! thd refers to changes induced by thermodynamics
+   ! trp   ''         ''     ''       advection (transport of ice)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: 
       d_a_i_thd  , d_a_i_trp ,                   &  !: icefractions                  
@@ -271 +451 @@
       d_u_ice_dyn, d_v_ice_dyn
 
-!---------------------------------------------------------------------------------------------------
-!  ) Ice thickness distribution variables
-!---------------------------------------------------------------------------------------------------
-
+   !!--------------------------------------------------------------------------
+   !! * Ice thickness distribution variables
+   !!--------------------------------------------------------------------------
+   ! REMOVE
    INTEGER(wp), PUBLIC, DIMENSION(jpl)                ::   &  !:
       ice_types      !: Vector making the connection between types and categories
 
    INTEGER(wp), PUBLIC, DIMENSION(jpm,2)              ::   &  !:
-      ice_cat_bounds !: Matrix containing the integer upper and lower boundaries of ice thickness categories
-
+      ice_cat_bounds !: Matrix containing the integer upper and 
+                     !: lower boundaries of ice thickness categories
+
+   ! REMOVE
    INTEGER(wp), PUBLIC, DIMENSION(jpm)                ::   &  !:
       ice_ncat_types !: Vector containing the number of thickness categories in each ice type
@@ -290 +472 @@
       hi_mean         !: Mean ice thickness in catgories 
 
+   ! REMOVE
    REAL(wp), PUBLIC, DIMENSION(0:jpl,jpm)         ::   &  !:
-      hi_max_typ     !: Boundary of ice thickness categories in thickness space (same but specific for each ice type)
-
-!---------------------------------------------------------------------------------------------------
-!  ) Ice diagnostics                     
-!---------------------------------------------------------------------------------------------------
+      hi_max_typ     !: Boundary of ice thickness categories 
+                     !:in thickness space (same but specific for each ice type)
+
+   !!--------------------------------------------------------------------------
+   !! * Ice diagnostics
+   !!--------------------------------------------------------------------------
+   !! Check if everything down here is necessary
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: volume of ice formed in the leads
       v_newice
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: thermodynamic growth rates 
       dv_dt_thd,  &
-      izero, fstroc, fhbricat ! to remove
-
+      izero, fstroc, fhbricat
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   & 
       diag_sni_gr,                           & ! snow ice growth 
@@ -309 +493 @@
       diag_bot_me,                           & ! vertical bottom melt 
       diag_sur_me                              ! vertical surface melt
-
    INTEGER , PUBLIC ::   &                      !: indexes of the debugging
       jiindex,           &                      !  point

trunk/NEMO/LIM_SRC_3/iceini.F90

@@ -22 +22 @@
    USE par_ice
    USE limvar
-   USE limicepoints
 
    IMPLICIT NONE
@@ -51 +50 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008) 
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/iceini.F90,v 1.4 2005/03/27 18:34:41 opalod Exp $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
@@ -64 +64 @@
       !!
       !! History :
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form and modules
+      !!   2.0  !  02-08  (G. Madec)  F90: Free form and modules
+      !!   3.0  !  08-03  (M. Vancop) ITD, salinity, EVP-C
       !!----------------------------------------------------------------------
-
-      INTEGER :: ji,jj,jk,jl, index
 
       ! Open the namelist file 
@@ -73 +72 @@
 
       CALL ice_run                    !  read in namelist some run parameters
-
-      CALL lim_icepoints              !  define specific checking ice points
                  
       ! Louvain la Neuve Ice model
@@ -87 +84 @@
       CALL lim_msh                    ! ice mesh initialization
      
-      CALL lim_itd_ini
+      CALL lim_itd_ini                ! initialize the ice thickness
+                                      ! distribution
       ! Initial sea-ice state
       IF( .NOT.ln_rstart ) THEN
          numit = 0
-         ! martin ajoute
          numit = nit000 - 1
          CALL lim_istate              ! start from rest: sea-ice deduced from sst
-         CALL lim_var_agg(1)
-         CALL lim_var_glo2eqv
+         CALL lim_var_agg(1)          ! aggregate category variables in
+                                      ! bulk variables
+         CALL lim_var_glo2eqv         ! convert global variables in equivalent
+                                      ! variables
       ELSE
          CALL lim_rst_read( numit )   ! start from a restart file
-         ! martin ajoute
          numit = nit000 - 1
-         CALL lim_var_agg(1)
-         CALL lim_var_glo2eqv
+         CALL lim_var_agg(1)          ! aggregate ice variables
+         CALL lim_var_glo2eqv         ! convert global var in equivalent variables
       ENDIF
 
@@ -130 +128 @@
       !! history :
       !!   2.0  !  03-08 (C. Ethe)  Original code
+      !!   3.0  !  08-03 (M. Vancop) LIM3
       !!-------------------------------------------------------------------
       NAMELIST/namicerun/ ln_limdyn, acrit, hsndif, hicdif, cai, cao, ln_nicep
@@ -158 +157 @@
         !!            Initializes the ice thickness distribution
         !! ** Method  :
-        !!            Very simple
+        !!            Very simple. Currently there are no ice types in the
+        !!            model...
         !!
         !! ** Arguments :
@@ -172 +172 @@
         !! ** History :
         !!           (12-2005) Martin Vancoppenolle
-        !!           Rien n'est jamais acquis
-        !!           A l'homme ni sa force
-        !!           Ni sa faiblesse ni son coeur
-        !!           Et quand il croit
-        !!           Ouvrir ses bras son ombre
-        !!           Est celle d'une croix
         !!
         !!------------------------------------------------------------------
@@ -183 +177 @@
 
        !! * Local variables
-       INTEGER ::   ji,       &   ! spatial dummy loop index
-                    jl,       &   ! ice category dummy loop index
+       INTEGER ::   jl,       &   ! ice category dummy loop index
                     jm            ! ice types    dummy loop index
 
@@ -198 +191 @@
 
 !!-- End of declarations
-!!----------------------------------------------------------------------------------------------
-
-!--------------------------------------------------------------------------------------------------!
-! 1) Ice thickness distribution parameters initialization                                          !
-!--------------------------------------------------------------------------------------------------!
-
-      !- Types boundaries in integer thickness space
-      !----------------------------------------------
-      !- Type 1 (undeformed ice) Boundaries
+!!------------------------------------------------------------------------------
+
+!------------------------------------------------------------------------------!
+! 1) Ice thickness distribution parameters initialization    
+!------------------------------------------------------------------------------!
+
+      !- Types boundaries (integer)
+      !----------------------------
       ice_cat_bounds(1,1) = 1
       ice_cat_bounds(1,2) = jpl
-!     !- Type 2 (ridged ice  ) Boundaries
-!     ice_cat_bounds(2,1) = 4
-!     ice_cat_bounds(2,2) = 5
-!     !- Type 3 (rafted ice  ) Boundaries
-!     ice_cat_bounds(3,1) = 6
-!     ice_cat_bounds(3,2) = 6
 
       !- Number of ice thickness categories in each ice type
@@ -223 +209 @@
       !- Make the correspondence between thickness categories and ice types
       !---------------------------------------------------------------------
-      !- ice_types = 1 -> undeformed ice
-      !- ice_types = 2 -> ridged   ice
-      !- ice_types = 3 -> rafted   ice
       DO jm = 1, jpm       !over types
          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories
@@ -242 +225 @@
       WRITE(numout,*)
 
-      !- Thickness categories boundaries in thickness space
-      !---------------------------------------------------------------------
-
-      ! NEW DEFINITION
-
-         hi_max(:) = 0.0
-         hi_max_typ(:,:) = 0.0
-
-! new categories
-!     !- Type 3 - rafted ice
-!        hi_max(ice_cat_bounds(3,2)) = 5.0
-
-!     !- Type 2 - ridged ice
-!        zc1 = 3./REAL(ice_cat_bounds(2,2)-ice_cat_bounds(2,1)+1)
-!        zc2 = 10.0*zc1
-!        zc3 = 3.0
-
-!        DO jl = ice_cat_bounds(2,1),ice_cat_bounds(2,2)
-!           zx1 = REAL(jl-ice_cat_bounds(2,1)+1) / REAL(ice_cat_bounds(2,2)-ice_cat_bounds(2,1)+1)
-!           hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1.0 + TANH ( zc3 * (zx1 - 1.0 ) ) )
-!        END DO
-
-!        ! force the thickness into the ice categories
-!        hi_max(ice_cat_bounds(2,1)) = 15.0
-!        hi_max(ice_cat_bounds(2,2)) = 25.0
+      !- Thickness categories boundaries 
+      !----------------------------------
+      hi_max(:) = 0.0
+      hi_max_typ(:,:) = 0.0
 
       !- Type 1 - undeformed ice
-         zc1 = 3./REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1)
-         zc2 = 10.0*zc1
-         zc3 = 3.0
-
-         DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-            zx1 = REAL(jl-1) / REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1)
-            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1.0 + TANH ( zc3 * (zx1 - 1.0 ) ) )
-         END DO
+      zc1 = 3./REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1)
+      zc2 = 10.0*zc1
+      zc3 = 3.0
+
+      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
+         zx1 = REAL(jl-1) / REAL(ice_cat_bounds(1,2)-ice_cat_bounds(1,1)+1)
+         hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1.0 + TANH ( zc3 * (zx1 - 1.0 ) ) )
+      END DO
 
       !- Fill in the hi_max_typ vector, useful in other circumstances
-         DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-            hi_max_typ(jl,1) = hi_max(jl)
-         END DO
-!        DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2)
-!           hi_max_typ(jl-ice_cat_bounds(2,1)+1,2) = hi_max(jl)
-!        END DO
-!        DO jl = ice_cat_bounds(3,1), ice_cat_bounds(3,2)
-!           hi_max_typ(jl-ice_cat_bounds(3,1)+1,3) = hi_max(jl)
-!        END DO
-
-         WRITE(numout,*) ' Thickness category boundaries independently of ice type '
-         WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
-
-         WRITE(numout,*) ' Thickness category boundaries inside ice types '
-         DO jm = 1, jpm
-            WRITE(numout,*) ' Type number ', jm
-            WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)
-         END DO
-
-         DO jl = 1, jpl
-            hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.0
-         END DO
-
-         tn_ice(:,:,:) = t_su(:,:,:)
+      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
+         hi_max_typ(jl,1) = hi_max(jl)
+      END DO
+
+      WRITE(numout,*) ' Thickness category boundaries independently of ice type '
+      WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
+
+      WRITE(numout,*) ' Thickness category boundaries inside ice types '
+      DO jm = 1, jpm
+         WRITE(numout,*) ' Type number ', jm
+         WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)
+      END DO
+
+      DO jl = 1, jpl
+         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.0
+      END DO
+
+      tn_ice(:,:,:) = t_su(:,:,:)
 
     END SUBROUTINE lim_itd_ini

trunk/NEMO/LIM_SRC_3/icestp.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   Lim sea-ice model
+   !!   'key_lim3' :                                   LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   ice_stp       : sea-ice model time-stepping
@@ -27 +27 @@
    USE limdyn
    USE limtrp
-   USE limicepoints
    USE limthd
    USE limitd_th
@@ -102 +101 @@
       !+++++
       index = 12
-      jiindex = arc_sp_grid(index,1)
-      jjindex = arc_sp_grid(index,2)
       jiindex = 44
       jjindex = 140

trunk/NEMO/LIM_SRC_3/limadv.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
+   !!   'key_lim3'                                     LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_adv_x  : advection of sea ice on x axis

trunk/NEMO/LIM_SRC_3/limcons.F90

@@ -1 +1 @@
 MODULE limcons
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                   LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!
    !!======================================================================
    !!                     ***  MODULE  limcons  ***
@@ -10 +15 @@
    !! 
    !!======================================================================
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
-   !!----------------------------------------------------------------------
    !!    lim_cons   :   checks whether energy/mass are conserved 
    !!----------------------------------------------------------------------
+   !!
    !! * Modules used
 
@@ -33 +36 @@
 
    !! * Module variables
-
-   !!----------------------------------------------------------------------
-   !!   LIM 2.1 , UCL-ASTR-LODYC-IPSL  (2004)
+   !!----------------------------------------------------------------------
+   !!   LIM 3.0 , UCL-ASTR-LODYC-IPSL  (2008)
    !!----------------------------------------------------------------------
 
@@ -222 +224 @@
    END SUBROUTINE lim_cons_check
 
-!==============================================================================
-
-   END MODULE limcons
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module            NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+   !!======================================================================
+END MODULE limcons

trunk/NEMO/LIM_SRC_3/limdia.F90

@@ -10 +10 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3' :                                   LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_dia      : computation of the time evolution of keys var.
@@ -95 +95 @@
       !! * Local variables
        INTEGER  ::   jv,ji,jj,jl ! dummy loop indices
-       INTEGER  ::   nv          ! indice of variable 
        REAL(wp), DIMENSION(jpinfmx) ::  & 
           vinfor           ! temporary working space 
@@ -138 +137 @@
                 vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
-!               vinfor(11) = vinfor(11) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !undef ice volume
-!               vinfor(13) = vinfor(13) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !def ice volume
                 vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
                 vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
+                ! the computation of this diagnostic is not reliable
                 vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
-                                                        v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
+                                                        v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 
                 vinfor(53) = vinfor(53) + fsalt(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux
                 vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux
                 vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux
-!slightly modified this
                 vinfor(59) = vinfor(59) + sst_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST
                 vinfor(61) = vinfor(61) + sss_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
@@ -157 +154 @@
                 vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
-!               vinfor(79) = vinfor(79) + v_i(ji,jj,6)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(79) = 0.0
                 vinfor(81) = vinfor(81) + fmass(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
-!               vinfor(83) = vinfor(83) + ht_i(ji,jj,5)*aire(ji,jj) / 1.0e12 ! mass flux
              ENDIF
           END DO
@@ -173 +168 @@
        END DO
 
-!      DO jl = ice_cat_bounds(2,1), ice_cat_bounds(3,2)
-!         DO jj = njeq, jpjm1
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               vinfor(13) = vinfor(13) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !def ice volume (inc. rafted ice)
-!            END DO
-!         END DO
-!      END DO
        vinfor(13) = 0.0
 
        vinfor(15) = vinfor(15) / vinfor(7) ! these have to be divided by total ice volume to have the
        vinfor(29) = vinfor(29) / vinfor(7) ! right value
-!      vinfor(31) = vinfor(31) / vinfor(7)
        vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) )
        vinfor(67) = vinfor(67) / vinfor(7)
@@ -192 +179 @@
        vinfor(57) = vinfor(57) / vinfor(5) ! 
        vinfor(79) = vinfor(79) / vinfor(5) !
-!      vinfor(83) = vinfor(83) / vinfor(5) !
 
        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) !
@@ -273 +259 @@
        !! Fram strait in ORCA2 = 5 points
        !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i
-!      jj = 137
        jj = 136 ! C grid
        vinfor(83) = 0.0
        vinfor(84) = 0.0
-!      WRITE(numout,*) ' Fram Strait Export '
        DO ji = 134, 138
-!         vinfor(83) = vinfor(83) - ( v_ice(ji,jj) + v_ice(ji+1,jj) ) / 2.0 * & 
-!                                     e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12
-!         vinfor(84) = vinfor(84) - ( v_ice(ji,jj) + v_ice(ji+1,jj) ) / 2.0 * & 
-!                                     e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12
-! C grid
           vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 
                                       e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12
           vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 
                                       e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12
-!         WRITE(numout,*)
-!         WRITE(numout,*) ' ji , jj : ', ji, jj
-!         WRITE(numout,*) ' v_ice ji  -jj  : ', v_ice(ji,jj)
-!         WRITE(numout,*) ' v_ice ji+1-jj  : ', v_ice(ji+1,jj)
-!         WRITE(numout,*) ' e1t(ji,jj) : ', e1t(ji,jj)
-!         WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj)
-!         WRITE(numout,*) ' rdt_ice  : ', rdt_ice
-!         WRITE(numout,*) ' area export : ', vinfor(83)
-!         WRITE(numout,*) ' vol  export : ', vinfor(84)
        END DO
 
@@ -311 +281 @@
                 vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
-!               vinfor(12) = vinfor(12) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !undef ice volume
-!               vinfor(14) = vinfor(14) + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !def ice volume
                 vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
                 vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
-!               vinfor(32) = vinfor(32) + u_ice(ji,jj)*u_ice(ji,jj) + v_ice(ji,jj)*v_ice(ji,jj) !ice vel
+                ! this diagnostic is not well computed (weighted by vol instead
+                ! of area)
                 vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
                                                         v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
@@ -324 +293 @@
                 vinfor(62) = vinfor(62) + sss_io(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
                 vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature
-                vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12       ! ice heat content
+                vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy
                 vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume
@@ -330 +299 @@
                 vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
                 vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
-!               vinfor(80) = vinfor(80) + v_i(ji,jj,6)*aire(ji,jj) / 1.0e12 ! mass flux
                 vinfor(80) = 0.0
                 vinfor(82) = vinfor(82) + fmass(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
-!               vinfor(84) = vinfor(84) + ht_i(ji,jj,5)*aire(ji,jj) / 1.0e12 ! mass flux
              ENDIF
           END DO
@@ -346 +313 @@
        END DO
 
-!      DO jl = ice_cat_bounds(2,1), ice_cat_bounds(3,2)
-!         DO jj = 2, njeqm1
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               vinfor(14) = vinfor(14) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !def ice volume (inc. rafted ice)
-!            END DO
-!         END DO
-!      END DO
        vinfor(14) = 0.0
 
        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8))) 
-       vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by total ice volume to have the
-       vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! right value
+       vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol
+       vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! 
        vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) )
-       vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! right value
+       vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! 
 
        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6))) 
-       vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by total ice extent to have the
-       vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! right value 
+       vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt
+       vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! 
        vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) ! 
        vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) !
@@ -408 +368 @@
        END DO
        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) !
-       vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! these have to be divided by total ice extent to have the
+       vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt
        !! 2.2) Diagnostics dependent on age
        !!------------------------------------
@@ -472 +432 @@
        !! history :
        !!  8.5  ! 03-08 (C. Ethe) original code
+       !!  9.0  ! 08-03 (M. Vancoppenolle) LIM3
        !!-------------------------------------------------------------------
        NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
@@ -479 +440 @@
           &          ndeb ,     &
           &          irecl
-
-       INTEGER  ::   nv            ! indice of variable 
 
        REAL(wp) ::   zxx0, zxx1    ! temporary scalars

trunk/NEMO/LIM_SRC_3/limdyn.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3' :                                 LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!    lim_dyn      : computes ice velocities
@@ -36 +36 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
@@ -58 +58 @@
       !!   1.0  !  01-04   (LIM)  Original code
       !!   2.0  !  02-08   (C. Ethe, G. Madec)  F90, mpp
-      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) LIM3, EVP, C-grid
+      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) 
+      !!                   LIM3, EVP, C-grid
       !!------------------------------------------------------------------------------------
       !! * Local variables
-      INTEGER ::   ji, jj, jl        ! dummy loop indices
+      INTEGER ::   ji, jj            ! dummy loop indices
       INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
-! nemo modif
-!        jhemis                      ! jhemis = 1 (NH) ; jhemis = -1 (SH)
+
       REAL(wp) ::   &
          ztairx, ztairy,          &  ! tempory scalars
@@ -70 +70 @@
          ztglx , ztgly ,          &
          zt11, zt12, zt21, zt22 , &
-         zustm, zsfrld, zsfrldm4, &
+         zustm,                   &
          zsfrldmx2, zsfrldmy2,    &
          zu_ice, zv_ice, ztair2
-! nemo modif
+
       REAL(wp),DIMENSION(jpj) ::   &
            zind,                     &  ! i-averaged indicator of sea-ice
@@ -86 +86 @@
       IF ( ln_limdyn ) THEN
 
-         ! Mean ice and snow thicknesses.          
-
+         ! ocean velocity
          u_oce(:,:)  = u_io(:,:) * tmu(:,:)
          v_oce(:,:)  = v_io(:,:) * tmv(:,:)
@@ -94 +93 @@
          old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
 
-         !                                         ! Rheology (ice dynamics)
-         !                                         ! ========
+         ! Rheology (ice dynamics)
+         ! ========
 
          !  Define the j-limits where ice rheology is computed
@@ -105 +104 @@
             IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
             CALL lim_rhg( i_j1, i_jpj )
-
          ELSE                                 ! optimization of the computational area
 
@@ -124 +122 @@
                i_j1 = MAX( 1, i_j1-1 )
                IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-
                CALL lim_rhg( i_j1, i_jpj )
 
@@ -169 +166 @@
          ENDIF
 
-         !                                         ! Ice-Ocean stress
-         !                                         ! ================
+         ! Ice-Ocean stress
+         ! ================
          DO jj = 2, jpjm1
-!           jhemis = SIGN(1, jj - jeq )
-!           zsang  = REAL(jhemis) * sangvg
-!nemo new version modif
             zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
 
@@ -193 +187 @@
                ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
 
-!              ! test on oscillations
-!              ztairx = 0.0
-!              ztairy = 0.0
-
                zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
                zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
@@ -208 +198 @@
                ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
                ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
-
-!              ! test on oscillations
-!              ztglx  = 0.0
-!              ztgly  = 0.0
+!
+!              ! IMPORTANT
+!              ! these lignes are bound to prevent numerical oscillations
+!              ! in the ice-ocean stress
+!              ! They are physically ill-based. There is a cleaner solution
+!              ! to try (remember discussion in Paris Gurvan)
+!
                ztglx   = ztglx * exp( - zmod / 0.5 )
                ztgly   = ztglx * exp( - zmod / 0.5 )
@@ -281 +274 @@
         CALL prt_ctl(tab2d_1=at_i  , clinfo1=' lim_dyn  : at_i  :')
       ENDIF
-
-!OPA9_DYN LIM2.1 2005
-!END OPA9_DYN LIM2.1 2005
 
    END SUBROUTINE lim_dyn

trunk/NEMO/LIM_SRC_3/limflx.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_flx  : flux at the ice / ocean interface
@@ -21 +21 @@
    USE in_out_manager
    USE albedo
-   USE limicepoints
    USE par_ice
    USE prtctl           ! Print control
@@ -69 +68 @@
       !! * Modules used
       !! * Local variables
-      INTEGER ::   ji, jj, jl         ! dummy loop indices
+      INTEGER ::   ji, jj             ! dummy loop indices
 
       INTEGER ::   &
@@ -80 +79 @@
 !!         zfcm2  ,              &  !  non solar heat fluxes
            zfold  ,            &
-           zdfseqv,            &    ! increment to the ice salt flux
 #if defined key_lim_fdd   
          zfons,                &  ! salt exchanges at the ice/ocean interface

trunk/NEMO/LIM_SRC_3/limhdf.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_hdf  : diffusion trend on sea-ice variable

trunk/NEMO/LIM_SRC_3/limistate.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3' :                                    LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_istate      :  Initialisation of diagnostics ice variables
@@ -21 +21 @@
    USE dom_ice
    USE ice
-   USE limicepoints
    USE lbclnk
 
@@ -72 +71 @@
 
       !! * Local variables
-      INTEGER  ::   ji, jj, jk, jl, jm         ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
 
       REAL(wp) ::   zidto,    &  ! temporary scalar
-         zs0, ztf, zbin, zai, zvi, zvs, zmi, zms, &
+         zs0, ztf, zbin, &
          zeps6, zeps, ztmelts, &
          epsi06

trunk/NEMO/LIM_SRC_3/limitd_me.F90

@@ -2 +2 @@
 
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                    LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!
    !!======================================================================
    !!                       ***  MODULE limitd_me ***
@@ -7 +11 @@
    !!                     computation of changes in g(h)      
    !!======================================================================
-
+   !!
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -17 +21 @@
    USE ice_oce          ! ice variables
    USE thd_ice
-   USE limicepoints
    USE limistate
    USE in_out_manager
    USE ice
    USE par_ice
-   USE limthd_lab
    USE limthd_lac
    USE limvar
@@ -74 +76 @@
       REAL(wp) :: &                               
          Cp 
-
 !
 !-----------------------------------------------------------------------
-!     Ridging diagnostic arrays for history files
+! Ridging diagnostic arrays for history files
 !-----------------------------------------------------------------------
 !
@@ -88 +89 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM-@ 3.0,  UCL-ASTR (2005-2006)
+   !!   LIM-@ 3.0,  UCL-ASTR (2008)
    !!   (c) UCL-ASTR and Martin Vancoppenolle
    !!----------------------------------------------------------------------
@@ -114 +115 @@
         !! ** External : 
         !!
+        !! ** Steps :
+        !!  1) Thickness categories boundaries, ice / o.w. concentrations
+        !!     Ridge preparation
+        !!  2) Dynamical inputs (closing rate, divu_adv, opning)
+        !!  3) Ridging iteration
+        !!  4) Ridging diagnostics
+        !!  5) Heat, salt and freshwater fluxes
+        !!  6) Compute increments of tate variables and come back to old values
+        !!
         !! ** References : There are a lot of references and can be difficult / 
         !!                 boring to read
@@ -138 +148 @@
         !!
         !! ** History :
-        !!           This routine was freely inpired from CICE
-        !!           authors William H. Lipscomb, LANL
-        !!                   Elizabeth C. Hunke (LANL)
-        !!                the routines were available online 
-        !!                so, big thanks to the authors and to LANL
-        !!           The modified features are the ridging/rafting algorithm
+        !!           This routine is based on CICE code
+        !!           and authors William H. Lipscomb,
+        !!           and Elizabeth C. Hunke, LANL
+        !!           are gratefully acknowledged
         !!
         !!           (02-2006) Martin Vancoppenolle, UCL-ASTR 
-        !!           Voila, c'est la derniere etape du codage!!! J'en ai  
-        !!           plein le cul mais je m'en sors mieux qu'avant!!! 
-        !!           je ne sais pas quand ca va marcher
-        !!           O toi, futur utilisateur, si tu pestes sur ce code
-        !!           Pardonne moi, sache que j'ai bien peste dessus aussi
         !!
         !!--------------------------------------------------------------------!
         !! * Arguments
 
-       !! * Local variables
-       INTEGER ::   ji,       &   ! spatial dummy loop index
-                    jj,       &   ! spatial dummy loop index
-                    jk,       &   ! vertical layering dummy loop index
-                    jl,       &   ! ice category dummy loop index
-                    jm,       &   ! ice types    dummy loop index
-                    index,    &   ! indexes for ice points debugging
-                    jbnd1,    &
-                    jbnd2,    &  
-                    niter,    &   ! iteration counter
-                    nitermax = 20 ! max number of ridging iterations (it used to
-!                    be 20
-
-       REAL(wp)  ::           &  ! constant values
-          zeps      =  1.0e-10, &
-          epsi10    =  1.0e-10, &
-          epsi06    =  1.0e-6
-
-       REAL(wp)  ::           &  ! constant values for ice enthalpy
-          zindb                
-
-       REAL(wp), DIMENSION(jpi,jpj) :: &
-          closing_net,        &  ! net rate at which area is removed    (1/s)
-                                 ! (ridging ice area - area of new ridges) / dt
-          divu_adv   ,        &  ! divu as implied by transport scheme  (1/s)
-          opning     ,        &  ! rate of opening due to divergence/shear
-          closing_gross,      &  ! rate at which area removed, not counting
-                                 ! area of new ridges
-          msnow_mlt  ,        &  ! mass of snow added to ocean (kg m-2)
-          esnow_mlt              ! energy needed to melt snow in ocean (J m-2)
-
-       REAL(wp) ::            &
-          w1,                 &  ! temporary variable
-          tmpfac,             &  ! factor by which opening/closing rates are cut
-          dti                    ! 1 / dt
-
-       LOGICAL   ::           &
-         iterate_ridging,     &  ! if true, repeat the ridging
-         asum_error              ! flag for asum .ne. 1
-
-       REAL(wp) ::  &         
-         big = 1.0e8
-
-      REAL (wp), DIMENSION(jpi,jpj) :: &  ! 
-              vt_i_init, vt_i_final,   &  !  ice volume summed over categories
-              vt_s_init, vt_s_final,   &  !  snow volume summed over categories
-              et_i_init, et_i_final,   &  !  ice energy summed over categories
-              et_s_init, et_s_final       !  snow energy summed over categories
-
-      CHARACTER (len = 15) :: fieldid
+        !! * Local variables
+        INTEGER ::   ji,       &   ! spatial dummy loop index
+                     jj,       &   ! spatial dummy loop index
+                     jk,       &   ! vertical layering dummy loop index
+                     jl,       &   ! ice category dummy loop index
+                     niter,    &   ! iteration counter
+                     nitermax = 20 ! max number of ridging iterations 
+
+        REAL(wp)  ::             &  ! constant values
+           zeps      =  1.0e-10, &
+           epsi10    =  1.0e-10, &
+           epsi06    =  1.0e-6
+
+        REAL(wp), DIMENSION(jpi,jpj) :: &
+           closing_net,        &  ! net rate at which area is removed    (1/s)
+                                  ! (ridging ice area - area of new ridges) / dt
+           divu_adv   ,        &  ! divu as implied by transport scheme  (1/s)
+           opning     ,        &  ! rate of opening due to divergence/shear
+           closing_gross,      &  ! rate at which area removed, not counting
+                                  ! area of new ridges
+           msnow_mlt  ,        &  ! mass of snow added to ocean (kg m-2)
+           esnow_mlt              ! energy needed to melt snow in ocean (J m-2)
+
+        REAL(wp) ::            &
+           w1,                 &  ! temporary variable
+           tmpfac,             &  ! factor by which opening/closing rates are cut
+           dti                    ! 1 / dt
+
+        LOGICAL   ::           &
+           iterate_ridging,    &  ! if true, repeat the ridging
+           asum_error              ! flag for asum .ne. 1
+
+        REAL(wp) ::  &         
+           big = 1.0e8
+
+        REAL (wp), DIMENSION(jpi,jpj) :: &  ! 
+           vt_i_init, vt_i_final       !  ice volume summed over categories
+
+        CHARACTER (len = 15) :: fieldid
 
 !!-- End of declarations
@@ -223 +215 @@
 ! Set hi_max(ncat) to a big value to ensure that all ridged ice 
 ! is thinner than hi_max(ncat).
+
       Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0    ! proport const for PE
-! BACK TO LIMA_MEC
-!     hi_max(ice_cat_bounds(2,2)) = big
-
-      CALL lim_itd_me_ridgeprep
-
-!     ! conservation check
-!     CALL lim_column_sum (jpl,   v_i, vt_i_init)
+      CALL lim_itd_me_ridgeprep ! prepare ridging
+
+      ! conservation check
+      IF ( con_i) CALL lim_column_sum (jpl,   v_i, vt_i_init)
 
 ! Initialize arrays.
@@ -288 +278 @@
       END DO
 
-!+++++ [
-      WRITE(numout,*) ' 2. Dynamical inputs '
-      WRITE(numout,*) ' closing_net ', closing_net(jiindex,jjindex)
-      WRITE(numout,*) ' divu_adv    ', divu_adv   (jiindex,jjindex)
-      WRITE(numout,*) ' opning      ', opning     (jiindex,jjindex)
-!+++++ ]
-
-
 !-----------------------------------------------------------------------------!
 ! 3) Ridging iteration
@@ -304 +286 @@
 
       DO WHILE ( iterate_ridging .AND. niter < nitermax )
-
-!+++++ [
-!     WRITE(numout,*)
-      WRITE(numout,*) ' 3. Ridging iteration : ', niter
-      WRITE(numout,*) ' ---------------------- '
-!     WRITE(numout,*)
-!     WRITE(numout,*) 'a_i ', a_i(jiindex, jjindex,1:jpl)
-!     WRITE(numout,*) 'v_i ', v_i(jiindex, jjindex,1:jpl)
-!+++++ ]
-
-! 3.1 ice / o.w. concentrationts, participation and transfer functions
 
       DO jj = 1, jpj
@@ -342 +313 @@
       END DO ! jj
 
-!+++++ [
-      WRITE(numout,*) ' correction 1'
-      WRITE(numout,*) ' closing_gross   ', closing_gross(jiindex,jjindex)
-      WRITE(numout,*) ' opning          ', opning       (jiindex,jjindex)
-!+++++ [
-
 ! correction to closing rate / opening if excessive ice removal
 !---------------------------------------------------------------
@@ -356 +321 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-!              IF (numit.GE.565) WRITE(numout,*) 'Stop Hier : ',ji,jj,jl
                IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
                   w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
@@ -369 +333 @@
       END DO !jl
 
-!+++++ [
-      WRITE(numout,*) ' correction 2'
-      WRITE(numout,*) ' closing_gross   ', closing_gross(jiindex,jjindex)
-      WRITE(numout,*) ' opning          ', opning(jiindex,jjindex) 
-!+++++ ]
-
 ! 3.3 Redistribute area, volume, and energy.
 !-----------------------------------------------------------------------------!
@@ -385 +343 @@
 
       CALL lim_itd_me_asumr
-
-!+++++ [
-      WRITE(numout,*)
-      WRITE(numout,*) ' back from asumr'
-      WRITE(numout,*) ' asum, ato_i, at_i : ', asum(jiindex,jjindex), &
-      ato_i(jiindex,jjindex), at_i(jiindex, jjindex)
-      WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,1:jpl)
-      WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,1:jpl)
-!+++++ ]
 
 ! 3.5 Do we keep on iterating ???
@@ -457 +406 @@
 !-----------------------------------------------------------------------------!
          ! fresh water source for ocean
-         fmmec(ji,jj)      = fmmec(ji,jj)      + msnow_mlt(ji,jj)*dti  !in kg / m2 /s !!! very important
-!        fresh_hist(ji,jj) = fresh_hist(ji,jj) + msnow_mlt(ji,jj)*dti
+         fmmec(ji,jj)      = fmmec(ji,jj)      + msnow_mlt(ji,jj)*dti  
       
          ! heat sink for ocean
          fhmec(ji,jj)      = fhmec(ji,jj)      + esnow_mlt(ji,jj)*dti
-!        fhnet_hist(ji,jj) = fhnet_hist(ji,jj) + esnow_mlt(ji,jj)*dti
 
          END DO
@@ -486 +433 @@
 
       ! Conservation check
-!     CALL lim_column_sum (jpl,   v_i, vt_i_final)
-!     fieldid = ' v_i : limitd_me '
-!     CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
+      IF ( con_i ) THEN
+         CALL lim_column_sum (jpl,   v_i, vt_i_final)
+         fieldid = ' v_i : limitd_me '
+         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
+      ENDIF
 
 !-----------------------------------------------------------------------------!
@@ -494 +443 @@
 !-----------------------------------------------------------------------------!
 
-     !+++++
-! BACK TO LIMA_MEC
-!     hi_max(ice_cat_bounds(2,2)) = 20.0
-     ! only for debugging
-     CALL lim_var_glo2eqv
-
-      WRITE(numout,*) ' End of lim_itd_me ****** '
-      WRITE(numout,*) ' u_ice      : ', u_ice(jiindex,jjindex)
-      WRITE(numout,*) ' v_ice      : ', v_ice(jiindex,jjindex)
-
-      DO jl = 1, jpl
-         WRITE(numout,*) '* - category number ', jl
-         WRITE(numout,*) ' a_i        : ', a_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' ht_i       : ', ht_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_i        : ', v_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_s        : ', v_s(jiindex,jjindex,jl)
-         WRITE(numout,*) ' e_s        : ', e_s(jiindex,jjindex,1,jl)/1.0e9
-         WRITE(numout,*) ' e_i1       : ', e_i(jiindex,jjindex,1,jl)/1.0e9
-         WRITE(numout,*) ' e_i2       : ', e_i(jiindex,jjindex,2,jl)/1.0e9
-         WRITE(numout,*) ' smv_i      : ', smv_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' sm_i       : ', sm_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' oa_i       : ', oa_i (jiindex,jjindex,jl)
-         WRITE(numout,*) ' o_i        : ', o_i (jiindex,jjindex,jl)
-         WRITE(numout,*) ' t_snow     : ', t_s(jiindex,jjindex,1,jl)
-         WRITE(numout,*) ' t_su       : ', t_su(jiindex,jjindex,jl)
-         WRITE(numout,*) ' t_i1       : ', t_i(jiindex,jjindex,1,jl)
-         WRITE(numout,*) ' t_i2       : ', t_i(jiindex,jjindex,2,jl)
-         WRITE(numout,*)
-      END DO
-
-!        DO jl = 1, jpl
-!           DO jj = 1, jpj
-!              DO ji = 1, jpi
-!                 IF (old_a_i(ji,jj,jl).eq.0.00) THEN
-!                     o_i(ji,jj,jl) = MAX( MIN( rdt_ice*float(numit)/86400.0, o_i(ji,jj,jl) ), 0.0 )
-!                 ENDIF
-!              END DO
-!           END DO
-!        END DO
-     !+++++
-
-!     ! Zero out categories with very small areas
+      CALL lim_var_glo2eqv
       CALL lim_itd_me_zapsmall
 
@@ -615 +523 @@
          END DO
       END DO
-
-
          
    END SUBROUTINE lim_itd_me
@@ -624 +530 @@
    SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6)
 
-        !!-----------------------------------------------------------------------------
+        !!----------------------------------------------------------------------
         !!                ***  ROUTINE lim_itd_me_icestrength ***
         !! ** Purpose :
@@ -634 +540 @@
         !! dissipated per unit area removed from the ice pack under compression,
         !! and assumed proportional to the change in potential energy caused
-        !! by ridging.  
+        !! by ridging. Note that only Hibler's formulation is stable and that
+        !! ice strength has to be smoothed
         !!
         !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
@@ -640 +547 @@
         !! ** External : 
         !!
-        !! ** References : There are a lot of references and they are difficult to read
-        !!                 so the best thing to do is to accept what is here
-        !!
-        !! ** History :
-        !!           This routine was freely inpired from CICE
-        !!           authors William H. Lipscomb, LANL
-        !!                   Elizabeth C. Hunke (LANL)
-        !!                the routines were available online 
-        !!                so, big thanks to the authors and to LANL
-        !!           The modified features are the ridging/rafting algorithm
-        !!
-        !!           (02-2006) Martin Vancoppenolle, UCL-ASTR 
-        !!           Voila, c'est la derniere etape du codage!!! J'en ai  
-        !!           plein le cul mais je m'en sors mieux qu'avant!!! 
-        !!------------------------------------------------------------------
+        !! ** References :
+        !!                
+        !!----------------------------------------------------------------------
         !! * Arguments
  
@@ -674 +569 @@
       REAL(wp), DIMENSION(jpi,jpj) :: &
          zworka              !: temporary array used here
-
-      WRITE(numout,*) ' lim_itd_me_icestrength : compute ice strength        '
-      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~ '
 
 !------------------------------------------------------------------------------!
@@ -707 +599 @@
                      hi * hi
 
-!                    strength(ji,jj) = strength(ji,jj)                         &
-!                    + athorn(ji,jj,jl)/krdg(ji,jj,jl)                         &
-!                    * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3)     &
-!                    / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))                      
                      !--------------------------
                      ! PE gain from rafting ice
@@ -720 +608 @@
                      ! PE gain from ridging ice
                      !----------------------------
-                     ! anything from ridge porosity ?
                      strength(ji,jj) = strength(ji,jj)                         &
                      + aridge(ji,jj,jl)/krdg(ji,jj,jl)                         &
@@ -753 +640 @@
          END DO                 ! i
 
-         WRITE(numout,*) ' Pstar    : ', Pstar
-         WRITE(numout,*) ' vt_i     : ', vt_i(jiindex,jjindex)
-         WRITE(numout,*) ' C_rhg    : ', C_rhg
-         WRITE(numout,*) ' at_i     : ', at_i(jiindex,jjindex)
-
          ksmooth = 1
 
@@ -766 +648 @@
 ! 5) Impact of brine volume
 !------------------------------------------------------------------------------!
+! CAN BE REMOVED
 !
       IF ( brinstren_swi .EQ. 1 ) THEN
-
-         WRITE(numout,*) ' mean brine volume : ', bv_i(jiindex,jjindex)
 
          DO jj = 1, jpj
@@ -778 +659 @@
                   zdummy = 0.0
                ENDIF
-!              zdummy = 0.0
-!              strength(ji,jj) = strength(ji,jj) * ( 1.0 - zdummy )
-! Timco and O'brien, CRST, 1994, Flexural strength equation for sea ice
                strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
             END DO              ! j
          END DO                 ! i
 
-         WRITE(numout,*) ' strength after brine volume correction : ', &
-                           strength(jiindex,jjindex)
-         
       ENDIF
 
@@ -861 +736 @@
       ENDIF ! ksmooth
       
-      WRITE(numout,*) ' Ice strength : ', strength(jiindex,jjindex)
-
       ! Boundary conditions
       CALL lbc_lnk( strength, 'T', 1. )
@@ -885 +758 @@
         !! ** External : 
         !!
-        !! ** History :
-        !!           This routine was freely inspired from CICE
-        !!           authors William H. Lipscomb, LANL
-        !!                   Elizabeth C. Hunke (LANL)
-        !!                the routines were available online 
-        !!                so, big thanks to the authors and to LANL
-        !!           The modified features are the ridging/rafting algorithm
-        !!
-        !!           (02-2006) Martin Vancoppenolle, UCL-ASTR 
-        !!           J'espere que nos lignes sont courbes
-        !!           Mais j'ai bien peur qu'elles ne soient droites
         !!---------------------------------------------------------------------!
         !! * Arguments
@@ -901 +763 @@
       INTEGER :: &
          ji,jj,  &          ! horizontal indices
-         jl,jl1, &          ! thickness category index
-         jm,     &          ! ice type index
-         index,  &
+         jl,     &          ! thickness category index
          krdg_index         ! which participation function using
 
@@ -913 +773 @@
          Gsum             ! Gsum(n) = sum of areas in categories 0 to n
 
-      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::    &
-         Gsumra           ! Gsum(n) = sum of areas in categories 0 to n
-
-      REAL(wp), DIMENSION(jpi,jpj,jpm)    ::    &
-         ai_typ           ! Surface of  respective ice types
-
       REAL(wp) ::    &
          hi,         &    ! ice thickness for each cat (m)
@@ -927 +781 @@
 
       REAL(wp)            ::     &
-         zindb, zdummy,          &
-         zweight,                &
+         zdummy,                 &
          epsi06 = 1.0e-6
 
-      INTEGER, DIMENSION(jpi,jpj,-1:jpl) :: &
-         index_cumul       ! temporary array used here
 !------------------------------------------------------------------------------!
 
@@ -958 +809 @@
       ! because of divergence during transport
 
-!+++++
-! 0.0 Compute surfaces of relative ice types
-! BACK TO LIMA_MEC
-!     ai_typ(:,:,:) = 0.0
-!     DO jm = 1, jpm
-!        DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-!           DO jj = 1, jpj
-!              DO ji = 1, jpi
-!                 ai_typ(ji,jj,jm) = ai_typ(ji,jj,jm) + a_i(ji,jj,jl) 
-!              END DO
-!           END DO
-!        END DO
-!     END DO
-!     WRITE(numout,*) ' ai_typ ', ai_typ(jiindex,jjindex,1:3)
-! BACK TO LIMA_MEC
-
       ! Compute cumulative thickness distribution function
       ! Compute the cumulative thickness distribution function Gsum,
       ! where Gsum(n) is the fractional area in categories 0 to n.
-
       ! initial value (in h = 0) equals open water area
 
@@ -1006 +840 @@
 
       ! Normalize the cumulative distribution to 1
-      !former line
-!     zworka(:,:) = 1.0 / Gsum(:,:,jpl)
-      WRITE(numout,*) ' AVANT : numit '
-      !new lines
       DO jj = 1, jpj 
          DO ji = 1, jpi
-!           IF (numit.GE.46365) THEN 
-            IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN
-               WRITE(numout,*) ' ji ,jj : ',ji,jj, 'Gsum:', &
-                                                 Gsum(ji,jj,jpl)
-            ENDIF
             zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl)
          END DO
       END DO
 
-      WRITE(numout,*) ' APRES '
       DO jl = 0, jpl
          DO jj = 1, jpj 
             DO ji = 1, jpi
-            IF ( (ji .EQ. jiindex) .AND. ( jj.EQ.jjindex ) ) THEN
-               WRITE(numout,*) ' ji ,jj : ',ji,jj,jl, 'Gsum:', &
-                                                 Gsum(ji,jj,jl)
-               WRITE(numout,*) ' zworka ', zworka(ji,jj)
-            ENDIF
-
                Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj)
             END DO 
          END DO
       END DO
-
-! BACK TO LIMA_MEC
-! For rafted ice
-!     jm = 3
-!     DO jj = 1, jpj
-!        DO ji = 1, jpi
-!           Gsumra(ji,jj,ice_cat_bounds(jm,1)-1) = 0.0
-!!          zworka(ji,jj)  = 1.0 / ai_typ(ji,jj,jm) ! used to normalize Gsum
-!        END DO
-!     END DO
-
-!     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-!        DO jj = 1, jpj 
-!           DO ji = 1, jpi
-!              Gsumra(ji,jj,jl) = Gsumra(ji,jj,jl-1) + a_i(ji,jj,jl)*zworka(ji,jj)
-!           END DO
-!        END DO
-!     END DO
-
-! +++++ [
-      WRITE(numout,*) 
-      WRITE(numout,*) ' Gsum : ', Gsum(jiindex,jjindex,-1:jpl)
-! +++++ ]
-
 
 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
@@ -1098 +892 @@
       ! precompute exponential terms using Gsum as a work array
       zdummy = 1.0 / (1.0-EXP(-astari))
-      WRITE(numout,*) ' zdummy : ', zdummy
 
       DO jl = -1, jpl
@@ -1119 +912 @@
       ENDIF ! krdg_index
 
-      WRITE(numout,*) ' astari : ', astari
-      WRITE(numout,*) ' athorn : ', athorn(jiindex,jjindex,0:jpl)
-
-! BACK TO LIMA_MEC
-!     jm = 3
-!     DO jl = ice_cat_bounds(3,1), ice_cat_bounds(3,2) ! only undeformed ice participates
-!        DO jj = 1, jpj 
-!           DO ji = 1, jpi
-!           IF (Gsumra(ji,jj,jl) < Gstar) THEN
-!              athorn(ji,jj,jl) = Gstari * (Gsumra(ji,jj,jl)-Gsumra(ji,jj,jl-1)) * &
-!                   (2.0 - (Gsumra(ji,jj,jl-1)+Gsumra(ji,jj,jl))*Gstari)
-!           ELSEIF (Gsumra(ji,jj,jl-1) < Gstar) THEN
-!              athorn(ji,jj,jl) = Gstari * (Gstar-Gsumra(ji,jj,jl-1)) *  &
-!                   (2.0 - (Gsumra(ji,jj,jl-1)+Gstar)*Gstari)
-!           ELSE
-!              athorn(ji,jj,jl) = 0.0
-!           ENDIF
-!           END DO ! ji
-!        END DO ! jj
-!     END DO ! jl 
-! BACK TO LIMA_MEC
-
-!!!
-! BACK TO LIMA_MEC
-      ! Weighting
-      ! Undeformed ice
-!     DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-!        DO jj = 1, jpj
-!           DO ji = 1, jpi
-!              zindb            = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )
-!              zweight          = zindb*ai_typ(ji,jj,1) /  &
-!                                 MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3), epsi06 )
-!              athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)  
-!           END DO
-!        END DO
-!     END DO
-
-      ! Rafted ice
-!     jm = 3
-!     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-!        DO jj = 1, jpj
-!           DO ji = 1, jpi
-!              zindb              = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi06 ) )
-!              zweight            = zindb*ai_typ(ji,jj,jm) /  & 
-!                                   MAX( ai_typ(ji,jj,1) + ai_typ(ji,jj,3) , epsi06)
-!              athorn(ji,jj,jl) = zweight*athorn(ji,jj,jl)  
-!           END DO
-!        END DO
-!     END DO
-      
-!+++++
-! BACK TO LIMA_MEC
-!!!
-      
       ! Ridging and rafting ice participation functions
       IF ( raftswi .EQ. 1 ) THEN
@@ -1224 +963 @@
 
       ENDIF
-!+++++
-      WRITE(numout,*)
-      WRITE(numout,*) ' aridge : ', aridge(jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' araft  : ', araft (jiindex, jjindex, 1:jpl)
-!+++++
 
 !-----------------------------------------------------------------
@@ -1290 +1024 @@
             DO ji = 1, jpi
                aksum(ji,jj)    = aksum(ji,jj)                          &
-!                      + athorn(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))
                        + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))    &
                        + araft (ji,jj,jl) * (1.0 - 1.0/kraft)
@@ -1297 +1030 @@
       END DO
 
-      WRITE(numout,*) ' krdg  : ', krdg(jiindex,jjindex,1:jpl)
-      WRITE(numout,*) ' hraft : ', hraft(jiindex,jjindex,1:jpl)
-      WRITE(numout,*) ' aksum : ', aksum(jiindex,jjindex)
-
       END SUBROUTINE lim_itd_me_ridgeprep
 
@@ -1324 +1053 @@
         !! ** External : 
         !!
-        !! ** History :
-        !!           This routine was freely inpired from CICE
-        !!           authors William H. Lipscomb, LANL
-        !!                   Elizabeth C. Hunke (LANL)
-        !!                the routines were available online 
-        !!                so, big thanks to the authors and to LANL
-        !!           The modified features are the ridging/rafting algorithm
-        !!
-        !!           (02-2006) Martin Vancoppenolle, UCL-ASTR 
-        !!     Dans 5 mois je suis sensé avoir fini ma thèse alors j'ai interet a cavaler
-        !!     sinon ca va mal se passer
-        !!     J'ai trop chaud, je suis fatigué
-        !!     Mon ame se desseche t elle ??
-        !!     Suis je en train de changer ???
-        !!     Oublie-je l'essentiel de mon histoire ??
-        !!
 
       REAL (wp), DIMENSION(jpi,jpj), INTENT(IN)   :: &
@@ -1355 +1068 @@
          jk,           &   ! ice layer index
          ij,           &   ! horizontal index, combines i and j loops
-         icells,       &   ! number of cells with aicen > puny
-         index             ! index for debug ice points
+         icells            ! number of cells with aicen > puny
 
       INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
@@ -1370 +1082 @@
          vsnon_init,            &  ! snow volume before ridging
          esnon_init,            &  ! snow energy before ridging
-         smv_i_init,             &  ! ice salinity before ridging
-         oa_i_init                  ! ice age before ridging
+         smv_i_init,            &  ! ice salinity before ridging
+         oa_i_init                 ! ice age before ridging
 
       REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: &
@@ -1380 +1092 @@
          ardg1      , &     ! area of ice ridged
          ardg2      , &     ! area of new ridges
-         virdg      , &     ! ice volume of ridging ice
          vsrdg      , &     ! snow volume of ridging ice
          esrdg      , &     ! snow energy of ridging ice
-         smrdg      , &     ! salinity of ridging ice
-         oirdg      , &     ! volumic age content of ridging ice
          oirdg1     , &     ! areal age content of ridged ice
          oirdg2     , &     ! areal age content of ridging ice
@@ -1397 +1106 @@
          srdg1      , &     ! sal*volume of ice ridged
          srdg2      , &     ! sal*volume of new ridges
-         smsw       , &     ! sal*volume of water trapped into ridges
-         totvrdg1   , &     ! total volume of ice ridged
-         totvrdg2           ! total volume of new ridges
+         smsw               ! sal*volume of water trapped into ridges
 
       REAL(wp), DIMENSION(jpi,jpj) ::           &
@@ -1409 +1116 @@
          esrft      , &     ! snow energy of rafting ice
          smrft      , &     ! salinity of rafting ice
-         oirft      , &     ! volumic age content of rafting ice
          oirft1     , &     ! areal age content of rafted ice
          oirft2             ! areal age content of rafting ice
 
       REAL(wp), DIMENSION(jpi,jpj,jkmax) ::    &
-         eirdg      , &     ! ice energy of ridging ice
          eirft      , &     ! ice energy of rafting ice
          erdg1      , &     ! enth*volume of ice ridged
@@ -1477 +1182 @@
          END DO !jj
       END DO !ji
-
-      !+++++
-!     WRITE(numout,*) ' 1. Change in open water area '
-!     WRITE(numout,*) ' closing_gross', closing_gross(jiindex, jjindex)
-!     WRITE(numout,*) ' opning       ', opning       (jiindex, jjindex)
-!     WRITE(numout,*) ' ato_i        ', ato_i        (jiindex, jjindex)
-!     WRITE(numout,*) 
-      !+++++ 
 
       ! if negative open water area alert it
@@ -1527 +1224 @@
       END DO !jl
 
-      !+++++
-      WRITE(numout,*) ' State variables before ridging is computed '
-      WRITE(numout,*) ' a_i         ', a_i(jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' v_i         ', v_i(jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' v_s         ', v_s(jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' e_s         ', e_s(jiindex, jjindex, 1, 1:jpl)
-      WRITE(numout,*) ' smv_i       ', smv_i(jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' oa_i        ', oa_i (jiindex, jjindex, 1:jpl)
-      WRITE(numout,*) ' e_i 1       ', e_i (jiindex, jjindex, 1, 1:jpl)
-      WRITE(numout,*) ' e_i 2       ', e_i (jiindex, jjindex, 2, 1:jpl)
-      WRITE(numout,*)
-!     !+++++
-
+!
 !-----------------------------------------------------------------
 ! 3) Pump everything from ice which is being ridged / rafted
@@ -1552 +1237 @@
       !------------------------------------------------
 
-      !+++++
-!     WRITE(numout,*)
-!     WRITE(numout,*) ' *** jl1 ***  ', jl1
-      !+++++
          icells = 0
          DO jj = 1, jpj
@@ -1571 +1252 @@
          large_afrft = .false.
 
-!        !+++++ 
-!        WRITE(numout,*) ' 3.1 ridging cells '
-!        WRITE(numout,*) ' icells     ', icells
-!        WRITE(numout,*)
-         !+++++
-
          DO ij = 1, icells
             ji = indxi(ij)
@@ -1594 +1269 @@
             oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
             oirft2(ji,jj)= oirft1(ji,jj) / kraft
-
-            !+++++
-!           IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN
-!              WRITE(numout,*) ' 3.2 area of ridging ice and of new ridge '
-!              WRITE(numout,*) ' aridge    ', aridge(ji,jj,jl1)
-!              WRITE(numout,*) ' araft     ', araft (ji,jj,jl1)
-!              WRITE(numout,*) ' krdg      ', krdg  (ji,jj,jl1)
-!              WRITE(numout,*) ' kraft     ', kraft
-!              WRITE(numout,*) ' clogross  ', closing_gross(ji,jj)
-!              WRITE(numout,*) ' ardg1 : ', ardg1(jiindex,jjindex)
-!              WRITE(numout,*) ' arft1 : ', arft1(jiindex,jjindex)
-!              WRITE(numout,*) ' ardg2 : ', ardg2(jiindex,jjindex)
-!              WRITE(numout,*) ' arft2 : ', arft2(jiindex,jjindex)
-!              WRITE(numout,*)
-!           ENDIF
-            !+++++
 
       !---------------------------------------------------------------
@@ -1629 +1288 @@
             ENDIF
 
-!           !+++++
-!           IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN
-!              WRITE(numout,*) ' 3.3 Ridging / rafting fractions '
-!              WRITE(numout,*) ' afrac : ', afrac(ji,jj)
-!              WRITE(numout,*) ' afrft : ', afrft(ji,jj)
-!              WRITE(numout,*)
-!           ENDIF
-!           !+++++ 
-
       !--------------------------------------------------------------------------
       ! 3.4) Subtract area, volume, and energy from ridging 
       !     / rafting category n1.
       !--------------------------------------------------------------------------
-!        jl1 looping 1-jpl
-!           ij looping 1-icells
-
-            ! ridging volumes, heat contents ...
-            ! virdg is the volume of ridged ice, at the end of ridging process
-            ! i.e., it includes 30% porosity
-!           virdg(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 
             vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) /             &
                            ( 1.0 + ridge_por )
@@ -1659 +1302 @@
                             ( 1. + ridge_por )
             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
-!           oirdg(ji,jj) = oa_i_init(ji,jj,jl1)  * afrac(ji,jj) * &
-!                          ( 1.0 - ridge_por )
 
             ! rafting volumes, heat contents ...
@@ -1667 +1308 @@
             esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
             smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
-!           oirft(ji,jj) = oa_i_init(ji,jj,jl1)  * afrft(ji,jj) 
 
             ! substract everything
@@ -1683 +1323 @@
             ! Salinity
             !-------------
-            ! salinity times volume of the seawater into ridges
-!           WRITE(numout,*) ' sss_io : ', sss_io(ji,jj)
-!           WRITE(numout,*) ' vsw    : ', vsw(ji,jj)
-!           WRITE(numout,*) ' ridge_por : ', ridge_por
-            
             smsw(ji,jj)       = sss_io(ji,jj) * vsw(ji,jj) * ridge_por 
 
@@ -1704 +1339 @@
             srdg2(ji,jj)      = sm_newridge * vrdg2(ji,jj) 
 
-
       !------------------------------------            
-      ! 3.5 Increment ridging diagnostics
+      ! 3.6 Increment ridging diagnostics
       !------------------------------------            
 
@@ -1720 +1354 @@
 
       !------------------------------------------            
-      ! 3.6 Put the snow somewhere in the ocean
+      ! 3.7 Put the snow somewhere in the ocean
       !------------------------------------------            
 
@@ -1742 +1376 @@
 
       !-----------------------------------------------------------------
-      ! 3.7) Compute quantities used to apportion ice among categories
+      ! 3.8 Compute quantities used to apportion ice among categories
       ! in the n2 loop below
       !-----------------------------------------------------------------
@@ -1754 +1388 @@
 
 
-         !++++++++
-         IF ( vrdg1(ji,jj) .GT. 0.0 ) THEN
-            IF ( ( vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por ) .GT. epsi10 ) ) & 
-            THEN 
-            WRITE(numout,*) ' ALERTE in ridging routine  '
-            WRITE(numout,*) ' difference : ', vrdg2(ji,jj) / vrdg1(ji,jj) - ( 1. + ridge_por )
-            WRITE(numout,*) ' vrdg2/vrdg1              ', vrdg2(ji,jj) / &
-            vrdg1(ji,jj)
-            WRITE(numout,*) ' should be equal to 1+p : ', 1 + ridge_por
-            WRITE(numout,*) ' vrdg1 : ', vrdg1(ji,jj)
-            WRITE(numout,*) ' vrdg2 : ', vrdg2(ji,jj)
-            ENDIF
-         ENDIF
-
-         !++++++++
-         !
-
          END DO                 ! ij
 
-!        !+++++
-!        !+++++
-
-!        !+++++
-!        WRITE(numout,*)
-!        WRITE(numout,*) ' Mass of salt of ridging ice 1 : ', srdg1(jiindex,jjindex)
-!        WRITE(numout,*) ' Mass of salt of ridging ice 2 : ', srdg2(jiindex,jjindex)
-!        WRITE(numout,*) ' Salinity of ridging ice     : ', smrdg(jiindex,jjindex) / & 
-!                                                           virdg(jiindex,jjindex) / (1.0-ridge_por)
-!        WRITE(numout,*) ' Mass of salt after ridging  : ', sm_i(ji,jj,jl1)
-!        WRITE(numout,*) ' Sal after removal of ridg.ice ', sm_i(ji,jj,jl1) / &
-!                                                           vicen_init(ji,jj,jl1)
-         !+++++
-
-
       !--------------------------------------------------------------------
-      ! 3.8) Compute ridging ice enthalpy, remove it from ridging ice and
+      ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
       !      compute ridged ice enthalpy 
       !--------------------------------------------------------------------
@@ -1838 +1440 @@
                   ji = indxi(ij)
                   jj = indxj(ij)
-            ! ridging diagnostics
-!           WRITE(numout,*) ' debut ', eice_init(ji,jj)
-!           WRITE(numout,*) ' erdg2 ', erdg2(ji,jj,jk)  
-!           WRITE(numout,*) ' erdg1 ', erdg1(ji,jj,jk)  
-!           WRITE(numout,*) ' ersw  ', ersw(ji,jj,jk)  
-!           WRITE(numout,*) ' vsw   ', vsw(ji,jj)  
                   eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - &
                   erdg1(ji,jj,jk)
@@ -1849 +1445 @@
             END DO !jk
          ENDIF
-
-!+++++ [
-!        WRITE(numout,*)
-!        WRITE(numout,*) ' ridging ice enthalpy ', jl1
-!        WRITE(numout,*) ' eird1     ', erdg1(jiindex,jjindex,1:nlay_i)
-!        WRITE(numout,*) ' eird2 old ', erdg2(jiindex,jjindex,1:nlay_i)
-!        WRITE(numout,*) ' eirft     ', eirft(jiindex,jjindex,1:2)
-!        WRITE(numout,*) ' e_i       ', e_i  (jiindex,jjindex,1:nlay_i,jl1)
-!+++++ ]
 
          IF (large_afrac) THEN  ! there is a bug
@@ -1884 +1471 @@
          ENDIF                  ! large_afrft
 
-!+++++ [
-!           IF ( ( ji.EQ.jiindex ) .AND. (jj.EQ.jjindex) ) THEN
-!              WRITE(numout,*) ' 3.4 Substract area, volume and energy '
-!              WRITE(numout,*) ' virdg : ', virdg(ji,jj)
-!              WRITE(numout,*) ' vsrdg     ', vsrdg(ji,jj)
-!              WRITE(numout,*) ' esrdg     ', esrdg(ji,jj)
-!              WRITE(numout,*) ' smrdg     ', smrdg(ji,jj)
-!              WRITE(numout,*) ' oirdg     ', oirdg(ji,jj)
-!              WRITE(numout,*)
-!              WRITE(numout,*) ' virft : ', virft(ji,jj)
-!              WRITE(numout,*) ' vsrft     ', vsrft(ji,jj)
-!              WRITE(numout,*) ' esrft     ', esrft(ji,jj)
-!              WRITE(numout,*) ' smrft     ', smrft(ji,jj)
-!              WRITE(numout,*) ' oirft     ', oirft(ji,jj)
-!              WRITE(numout,*)  
-
-!              WRITE(numout,*) ' a_i   : ', a_i(ji,jj,jl1)
-!              WRITE(numout,*) ' v_i   : ', v_i(ji,jj,jl1)
-!              WRITE(numout,*) ' v_s       ', v_s(ji,jj,jl1)
-!              WRITE(numout,*) ' e_s       ', e_s(ji,jj,1,jl1)
-!              WRITE(numout,*) ' smv_i      ', smv_i(ji,jj,jl1)
-!              WRITE(numout,*) ' oa_i       ', oa_i(ji,jj,jl1)
-!           ENDIF
-!+++++ ]
-
 !-------------------------------------------------------------------------------
 ! 4) Add area, volume, and energy of new ridge to each category jl2
@@ -1946 +1508 @@
 
             END DO ! ij
-!           WRITE(numout,*) ' new ridges being built from ', jl1
-!           WRITE(numout,*) ' added to ', jl2
-!           WRITE(numout,*) ' a_i   : ', a_i(jiindex,jjindex,jl2)
-!           WRITE(numout,*) ' v_i   : ', v_i(jiindex,jjindex,jl2)
 
             ! Transfer ice energy to category jl2 by ridging
@@ -1964 +1522 @@
          END DO                 ! jl2 (new ridges)            
 
-!        WRITE(numout,*) ' new ridges built '
-!        WRITE(numout,*) ' a_i   : ', a_i(jiindex,jjindex,jl1)
-!        WRITE(numout,*) ' v_i   : ', v_i(jiindex,jjindex,jl1)
-!        WRITE(numout,*) ' e_i   : ', e_i(jiindex,jjindex,:,jl1)
-
          DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
-               ! jl2  = ice_cat_bounds(3,1) ! rafted category
 
             DO ij = 1, icells
@@ -1990 +1542 @@
                   oa_i(ji,jj,jl2)  = oa_i(ji,jj,jl2)                  &
                                    + oirft2(ji,jj)    
-
-!                 IF ( (ji.EQ.jiindex) .AND. (jj.EQ.jjindex) ) THEN
-!                         WRITE(numout,*) ' raft ice from jl1 ', jl1 , &
-!                                         ' to jl2 :', jl2
-!                         WRITE(numout,*) ' hraft (jl1) ', hraft(ji,jj,jl1)
-!                         WRITE(numout,*) ' hi_maxs ', hi_max(jl2-1), hi_max(jl2)
-!                 ENDIF
-
                ENDIF ! hraft
 
@@ -2017 +1561 @@
          END DO ! jl2
 
-!        WRITE(numout,*) ' new rafted ice built '
-!        WRITE(numout,*) ' a_i   : ', a_i(jiindex,jjindex,jl1)
-!        WRITE(numout,*) ' v_i   : ', v_i(jiindex,jjindex,jl1)
-!        WRITE(numout,*) ' e_i   : ', e_i(jiindex,jjindex,:,jl1)
-
       END DO ! jl1 (deforming categories)
 
@@ -2060 +1599 @@
 
       INTEGER :: ji, jj, jl
-
-!     WRITE(numout,*) ' lim_itd_me_asumr : sum the ice / open water fractions'
-!     WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
 
       !-----------------------------------------------------------------
@@ -2172 +1708 @@
            xtmp      ! temporary variable
 
-      WRITE(numout,*) ' lim_itd_me_zapsmall '
-      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~ '
-      
       DO jl = 1, jpl
 

trunk/NEMO/LIM_SRC_3/limitd_th.F90

@@ -1 +1 @@
 MODULE limitd_th
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                   LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limitd_th ***
@@ -15 +18 @@
    USE ice_oce          ! ice variables
    USE thd_ice
-   USE limicepoints
    USE limistate
    USE in_out_manager
    USE ice
    USE par_ice
-   USE limthd_lab
    USE limthd_lac
    USE limvar
@@ -44 +45 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM-@ 4.0,  UCL-ASTR (2005)
+   !!   LIM3.0,  UCL-ASTR (2008)
    !!   (c) UCL-ASTR and Martin Vancoppenolle
    !!----------------------------------------------------------------------
@@ -73 +74 @@
         !! ** History :
         !!           (12-2005) Martin Vancoppenolle 
-        !!            Au moment ou j'ecris ces lignes, je ne me rends pas 
-        !!             compte du boulot que j'entame... un truc de malate comme
-        !!             on dit ici chez les Belches
         !!
         !!------------------------------------------------------------------
@@ -81 +79 @@
 
        !! * Local variables
-       INTEGER ::   ji,       &   ! spatial dummy loop index
-                    jj,       &   ! spatial dummy loop index
-                    jk,       &   ! vertical layering dummy loop index
-                    jl,       &   ! ice category dummy loop index
-                    jm,       &   ! ice types    dummy loop index
-                    index,    &
+       INTEGER ::   jm,       &   ! ice types    dummy loop index
                     jbnd1,    &
                     jbnd2
@@ -94 +87 @@
           epsi10    =  1.0e-10
 
-       REAL(wp)  ::           &  ! constant values for ice enthalpy
-          zindb                
-
 !!-- End of declarations
 !!----------------------------------------------------------------------------------------------
 
-       IF(lwp) THEN
+       IF (lwp) THEN
           WRITE(numout,*)
-          WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'
+          WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution'
           WRITE(numout,*) '~~~~~~~~~~~'
        ENDIF
@@ -120 +110 @@
       CALL lim_var_agg(1)
 
-!+++++
-      WRITE(numout,*) ' From limitd_th : ' 
-      WRITE(numout,*) ' at_i       : ', at_i(jiindex,jjindex)
-      WRITE(numout,*) ' vt_i       : ', vt_i(jiindex,jjindex)
-      WRITE(numout,*) ' vt_s       : ', vt_s(jiindex,jjindex)
-      DO jl = 1, jpl
-         WRITE(numout,*) '* - category number ', jl
-         WRITE(numout,*) ' a_i        : ', a_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_i        : ', v_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' ht_i       : ', ht_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_s        : ', v_s(jiindex,jjindex,jl)
-         WRITE(numout,*) ' ht_s       : ', ht_s(jiindex,jjindex,jl)
-         WRITE(numout,*) ' e_s        : ', e_s(jiindex,jjindex,1,jl)/1.0e9
-         WRITE(numout,*) ' e_i        : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9
-         WRITE(numout,*) ' t_su       : ', t_su(jiindex,jjindex,jl)
-         WRITE(numout,*) ' t_snow     : ', t_s(jiindex,jjindex,1,jl)
-         WRITE(numout,*) ' t_i        : ', t_i(jiindex,jjindex,1:nlay_i,jl)
-         WRITE(numout,*) ' smv_i      : ', smv_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' oa_i       : ', oa_i(jiindex,jjindex,jl)
-         WRITE(numout,*)
-      END DO
-!+++++
-
-!------------------------------------------------------------------------------|
-!  2) Melt ice laterally.
-!------------------------------------------------------------------------------|
-!     DO jm = 1, jpm
-!        CALL lim_thd_lab(ice_cat_bounds(jm,1),ice_cat_bounds(jm,2))
-!     END DO
-!     CALL lim_thd_lab
-
 !------------------------------------------------------------------------------|
 !  3) Add frazil ice growing in leads.
@@ -158 +117 @@
       CALL lim_var_glo2eqv ! only for info
 
-!+++++
-      WRITE(numout,*) ' limthd_lac, new values ***** '
-      DO jl = 1, jpl
-         WRITE(numout,*) '* - category number ', jl
-         WRITE(numout,*) ' a_i        : ', a_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' ht_i       : ', ht_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_i        : ', v_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' v_s        : ', v_s(jiindex,jjindex,jl)
-         WRITE(numout,*) ' e_i        : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9
-         WRITE(numout,*) ' smv_i      : ', smv_i(jiindex,jjindex,jl)
-         WRITE(numout,*) ' t_su       : ', t_su(jiindex,jjindex,jl)
-         WRITE(numout,*) ' t_snow     : ', t_s(jiindex,jjindex,1,jl)
-         WRITE(numout,*) ' t_i        : ', t_i(jiindex,jjindex,1:nlay_i,jl)
-         WRITE(numout,*)
-      END DO
-!+++++
 !----------------------------------------------------------------------------------------
 !  4) Computation of trend terms and get back to old values      
@@ -230 +173 @@
         !!           (06-2006) Adaptation to include salt, age and types
         !!           (04-2007) Mass conservation checked
-        !!
-        !!                     Je suis d'humeur massacrante aujourd'hui, tout 
-        !!                     le monde m'embete et m'empeche de coder
-        !!
-        !!                     Muere lentamente
-        !!                     quien evita una pasion y su remolino de emociones
-        !!                     justamente estas que regresan el brillo a los ojos
-        !!                     y restauran los corazones destrozados
-        !!
         !!------------------------------------------------------------------
         !! * Arguments
@@ -250 +184 @@
        INTEGER ::   ji,       &   ! spatial dummy loop index
                     jj,       &   ! spatial dummy loop index
-                    jk,       &   ! vertical layering dummy loop index
                     jl,       &   ! ice category dummy loop index
-                    index,    &   ! for ice points
                     zji, zjj, &   ! dummy indices used when changing coordinates
                     nd            ! used for thickness categories
@@ -307 +239 @@
        
        REAL(wp)  ::           &  ! constant values for ice enthalpy
-          zdummy, zdummy2,    &
           zslope                 ! used to compute local thermodynamic "speeds"
 
@@ -337 +268 @@
 !! 1) Compute thickness and changes in each ice category
 !!----------------------------------------------------------------------------------------------
-       IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'
-         WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '                 klbnd :       ', klbnd
-         WRITE(numout,*) '                 kubnd :       ', kubnd
-         WRITE(numout,*) '                 ntyp  :       ', ntyp 
+       IF (lwp) THEN
+       WRITE(numout,*)
+       WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution'
+       WRITE(numout,*) '~~~~~~~~~~~~~~~'
+       WRITE(numout,*) ' klbnd :       ', klbnd
+       WRITE(numout,*) ' kubnd :       ', kubnd
+       WRITE(numout,*) ' ntyp  :       ', ntyp 
        ENDIF
-
-! +++++ [
-!      index = 1
-!      jiindex = arc_sp_grid(index,1)
-!      jjindex = arc_sp_grid(index,2)
-!      WRITE(numout,*) '*', arc_sp_acro(index), ' ', arc_sp_name(index)
-!      WRITE(numout,*)
-!      WRITE(numout,*) ' a_i        : ', a_i(jiindex,jjindex,klbnd:kubnd)
-!      WRITE(numout,*) ' ht_i       : ', ht_i(jiindex,jjindex,klbnd:kubnd)
-!      WRITE(numout,*) ' v_i        : ', v_i(jiindex,jjindex,klbnd:kubnd)
-! +++++ ]
 
        zdhice(:,:,:) = 0.0
@@ -361 +281 @@
           DO jj = 1, jpj
              DO ji = 1, jpi
-
                 zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes
                 ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb
@@ -369 +288 @@
                    zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 
                 ENDIF
-
              END DO
           END DO
        END DO
-!+++++
-!      WRITE(numout,*) ' 1 *** '
-!      WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd
-!      WRITE(numout,*) 'ht_i   ', ht_i   (jiindex,jjindex,klbnd:kubnd)
-!      WRITE(numout,*) 'zht_i_o', zht_i_o(jiindex,jjindex,klbnd:kubnd)
-!      WRITE(numout,*) 'zdhice ', zdhice (jiindex,jjindex,klbnd:kubnd)
-!+++++
 
 !-----------------------------------------------------------------------------------------------
@@ -392 +303 @@
          END DO
       END DO
-
-!+++++
-!     WRITE(numout,*) ' 2 *** '
-!     WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd
-!     WRITE(numout,*) 'a_i    ', a_i   (jiindex,jjindex,klbnd:kubnd)
-!     WRITE(numout,*) 'at_i   ', at_i  (jiindex,jjindex)
-!+++++
 
 !-----------------------------------------------------------------------------------------------
@@ -478 +382 @@
          ! ji
       END DO !jl
-!+++++
-!     WRITE(numout,*) ' 4 *** '
-!     WRITE(numout,*) ' klbnd -> kubnd - 1 ', klbnd, kubnd - 1
-!     WRITE(numout,*) ' hi_max ', hi_max(klbnd:kubnd-1)
-!     WRITE(numout,*) ' zhbnew ', zhbnew(jiindex,jjindex,klbnd:kubnd-1)
-!+++++
             
 !-----------------------------------------------------------------------------------------------
@@ -521 +419 @@
      END DO !jj
 
-!+++++
-!    WRITE(numout,*) ' 6 *** '
-!    WRITE(numout,*) ' klbnd -1, kubnd ', klbnd - 1, kubnd 
-!    WRITE(numout,*) ' zhb0   ', zhb0(jiindex,jjindex)
-!    WRITE(numout,*) ' zhb1   ', zhb1(jiindex,jjindex)
-!    WRITE(numout,*) ' zhbnew klbnd-1 ', zhbnew(jiindex,jjindex,klbnd-1)
-!    WRITE(numout,*) ' zhbnew kubnd   ', zhbnew(jiindex,jjindex,klbnd)
-!+++++
-
 !-----------------------------------------------------------------------------------------------
 !  7) Compute g(h) 
@@ -537 +426 @@
                           g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
                           hR(:,:,klbnd), zremap_flag)
-
-!+++++
-!    WRITE(numout,*) ' 7a *** klbnd ', klbnd
-!    WRITE(numout,*) ' g0(klbnd) ', g0(jiindex,jjindex,klbnd)
-!    WRITE(numout,*) ' g1(klbnd) ', g1(jiindex,jjindex,klbnd)
-!    WRITE(numout,*) ' hL(klbnd) ', hL(jiindex,jjindex,klbnd)
-!    WRITE(numout,*) ' hR(klbnd) ', hR(jiindex,jjindex,klbnd)
-!+++++
 
      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd)
@@ -601 +482 @@
      END DO
 
-!+++++
-!    WRITE(numout,*) ' 7b *** klbnd->kubnd ', klbnd, kubnd
-!    WRITE(numout,*) ' g0 ', g0(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' g1 ', g1(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' hL ', hL(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' hR ', hR(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*)
-!    WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' a_i  ', a_i (jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' v_i  ', v_i (jiindex,jjindex,klbnd:kubnd)
-!+++++
-
 !-----------------------------------------------------------------------------------------------
 !  8) Compute area and volume to be shifted across each boundary
@@ -669 +538 @@
      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
 
-!    WRITE(numout,*) ' 9 *** '
-!    WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' a_i  ', a_i (jiindex,jjindex,klbnd:kubnd)
-!    WRITE(numout,*) ' v_i  ', v_i (jiindex,jjindex,klbnd:kubnd)
-
 !!----------------------------------------------------------------------------------------------
 !! 10) Make sure ht_i >= minimum ice thickness hi_min
@@ -741 +605 @@
         !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
         !!          (01-2006) Martin Vancoppenolle 
-        !!          Au moment ou j'ecris ces lignes, je ne me rends pas 
-        !!          compte du boulot que j'entame... un truc de malate comme
-        !!          on dit ici chez les Belches
-        !!          This routine was inspired from CICE (W.H. Lipscomb, E. Hunke, C. M. Bitz)
-        !!          the sea ice model of LANL, Los Alamos, USA. 
-        !!          Thanks to these guys and their team to put their routines online
         !!
         !!------------------------------------------------------------------
@@ -783 +641 @@
 !!-- End of declarations
 !!----------------------------------------------------------------------------------------------
-      
-!     WRITE(numout,*) ' lim_itd_fitline : linearly fitting the function g(h) '
-!     WRITE(numout,*) ' ~~~~~~~~~~~~~~~   in category number ', num_cat
 
       DO jj = 1, jpj
@@ -853 +708 @@
         !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
         !!          (01-2006) Martin Vancoppenolle 
-        !!          Et quand fatigues de s'etre souvenus
-        !!          Nos souvenirs fatigues ne seront plus que des loques...
-        !!          "!!! Il a attrape la ratatinette, l'epouvantable ratatinette"
-        !!
-        !!          This routine was largely inspired from CICE 
-        !!          (W.H. Lipscomb, E. Hunke, C. M. Bitz)
-        !!          the sea ice model of LANL, Los Alamos, USA. 
-        !!          Merci a eux et a leur equipe de mettre leurs routines en ligne
-        !!          
-        !!          J'ajoute aujourd'hui : le secret de ma vie pour l'instant
-        !!          c'est de trouver la balance entre mon ego et l'absence d'ego
-        !!          balance entre respect des autres et respect de moi
-        !!          trouver comment realiser la compréhension de l'autre
-        !!          sans déclencher ma propre aliénation 
         !!
         !!------------------------------------------------------------------
@@ -887 +728 @@
          jl1,           &  ! donor category
          jk,            &  ! ice layer index
-         zji, zjj,      &  ! indices when changing from 2D-1D is done
-         index             ! for ice points referencing
+         zji, zjj          ! indices when changing from 2D-1D is done
 
       REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
-         zaTsfn,        &
-         zqsnow
+         zaTsfn
 
       REAL(wp), DIMENSION(jpi,jpj) :: &
@@ -902 +741 @@
          zdeice      ,  &  ! ice energy transferred
          zdsm_vice      ,  &  ! ice salinity times volume transferred
-         zsm_v1      ,  &  ! ice salinity times volume
-         zsm_v2      ,  &  ! ice salinity times volume
          zdo_aice      ,  &  ! ice age times volume transferred
-         zo_v1       ,  &  ! ice age times volume 
-         zo_v2       ,  &  ! ice age times volume 
          zdaTsf      ,  &  ! aicen*Tsfcn transferred
          zindsn      ,  &  ! snow or not
@@ -928 +763 @@
 
 !!-- End of declarations
-!      WRITE(numout,*) ' lim_itd_shiftice : shifting ice between categories '
-!      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
+
 !----------------------------------------------------------------------------------------------
 ! 1) Define a variable equal to a_i*T_su
@@ -1024 +858 @@
       END DO !jl
 
-      !-----------------------------------------------------------------
-      ! error messages
-      !-----------------------------------------------------------------
-
-!        if (daice_negative) then
-!           do j = jlo,jhi
-!           do i = ilo,ihi
-!              if (donor(i,j,n) > 0 .and. daice(i,j,n) <= -puny) then
-!                 write(nu_diag,*) my_task,':',i,j,
-!    &                 'ITD Neg daice =',daice(i,j,n),' boundary',n
-!                 call abort_ice ('ice: ITD Neg daice')
-!              endif
-!           enddo
-!           enddo
-!        endif
-
-!        if (dvice_negative) then
-!           do j = jlo,jhi
-!           do i = ilo,ihi
-!              if (donor(i,j,n) > 0 .and. dvice(i,j,n) <= -puny) then
-!                 write(nu_diag,*) my_task,':',i,j,
-!    &                 'ITD Neg dvice =',dvice(i,j,n),' boundary',n
-!                 call abort_ice ('ice: ITD Neg dvice')
-!              endif
-!           enddo
-!           enddo
-!        endif
-
-!        if (daice_greater_aicen) then
-!           do j = jlo,jhi
-!           do i = ilo,ihi
-!              if (donor(i,j,n) > 0) then
-!                 n1 = donor(i,j,n)
-!                 if (daice(i,j,n) >= aicen(i,j,n1)+puny) then
-!                    write(nu_diag,*) my_task,':',i,j,
-!    &                    'ITD daice > aicen, cat',n1
-!                    write(nu_diag,*) my_task,':',i,j,
-!    &                    'daice =', daice(i,j,n),
-!    &                    'aicen =', aicen(i,j,n1)
-!                    call abort_ice ('ice: ITD daice > aicen')
-!                 endif
-!              endif
-!           enddo
-!           enddo
-!        endif
-
-!        if (dvice_greater_vicen) then
-!           do j = jlo,jhi
-!           do i = ilo,ihi
-!              if (donor(i,j,n) > 0) then
-!                 n1 = donor(i,j,n)
-!                 if (dvice(i,j,n) >= vicen(i,j,n1)+puny) then
-!                    write(nu_diag,*) my_task,':',i,j,
-!    &                    'ITD dvice > vicen, cat',n1
-!                    write(nu_diag,*) my_task,':',i,j,
-!    &                    'dvice =', dvice(i,j,n),
-!    &                    'vicen =', vicen(i,j,n1)
-!                    call abort_ice ('ice: ITD dvice > vicen')
-!                 endif
-!              endif
-!           enddo
-!           enddo
-!        endif
-
-!     enddo                     ! boundaries 1 to ncat-1
-
 !-------------------------------------------------------------------------------
 ! 3) Transfer volume and energy between categories
@@ -1210 +978 @@
                t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
                zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes
-!              t_s(ji,jj,1,jl) =  - zQsnow(ji,jj,jl) / MAX(v_s(ji,jj,jl),zeps) * zindsn + rtt
             ELSE
                ht_i(ji,jj,jl)  = 0.0
                t_su(ji,jj,jl)  = rtt
-!              t_s(ji,jj,1,jl) = rtt
             ENDIF
          END DO                 ! ji
@@ -1246 +1012 @@
         !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
         !!          (01-2006) Martin Vancoppenolle (adaptation)
-        !!
-        !!          Quel etre encore que celui-ci! Le Jugement Dernier sera la
-        !!          avant qu'il vous fasse jamais une avance sur votre mois,
-        !!          Seigneur! Tu peux supplier, te mettre en quatre,
-        !!          meme si tu es dans la misere, il ne te donnera rien,
-        !!          le vieux demon! Et quant on pense que, chez lui,
-        !!          sa cuisiniere lui donne des gifles! Je ne vois pas l'interet
-        !!          qu'il y a a travailler dans un ministere. Cela ne rapporte
-        !!          absolument rien.
-        !!
-        !!          This routine was largely inspired from CICE 
-        !!          (W.H. Lipscomb, E. Hunke, C. M. Bitz)
-        !!          the sea ice model of LANL, Los Alamos, USA. 
         !!
         !!------------------------------------------------------------------
@@ -1283 +1036 @@
       REAL(wp)  ::           &  ! constant values
          zeps      =  1.0e-10, &
-         epsi10    =  1.0e-10, &
-         zindb
+         epsi10    =  1.0e-10
 
       REAL (wp), DIMENSION(jpi,jpj) :: &  ! 
@@ -1294 +1046 @@
 !!-- End of declarations
 !------------------------------------------------------------------------------
-!     WRITE(numout,*) ' lim_itd_th_reb '
-!     WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
-!     WRITE(numout,*) ' Ice Type no    ', ntyp
-!     WRITE(numout,*) ' bounds of categories ', klbnd, kubnd
 
 !     ! conservation check
-!     CALL lim_column_sum (jpl,   v_i, vt_i_init)
-!     CALL lim_column_sum (jpl,   v_s, vt_s_init)
+      IF ( con_i ) THEN
+         CALL lim_column_sum (jpl,   v_i, vt_i_init)
+         CALL lim_column_sum (jpl,   v_s, vt_s_init)
+      ENDIF
 
 !
@@ -1307 +1057 @@
 ! 1) Compute ice thickness.
 !------------------------------------------------------------------------------
-! nothing to do
       DO jl = klbnd, kubnd
          DO jj = 1, jpj
@@ -1447 +1196 @@
 !------------------------------------------------------------------------------
 
-!   CALL lim_column_sum (jpl,   v_i, vt_i_final)
-!   fieldid = ' v_i : limitd_reb '
-!   CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-
-!   CALL lim_column_sum (jpl,   v_s, vt_s_final)
-!   fieldid = ' v_s : limitd_reb '
-!   CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
+    IF ( con_i ) THEN
+       CALL lim_column_sum (jpl,   v_i, vt_i_final)
+       fieldid = ' v_i : limitd_reb '
+       CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
+
+       CALL lim_column_sum (jpl,   v_s, vt_s_final)
+       fieldid = ' v_s : limitd_reb '
+       CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
+    ENDIF
 
     END SUBROUTINE lim_itd_th_reb

trunk/NEMO/LIM_SRC_3/limmsh.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_msh   : definition of the ice mesh
@@ -24 +24 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008) 
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limmsh.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
@@ -75 +75 @@
       fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
  
-!i    DO jj = 1, jpj
-!i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line
-!!ii     write(numout,*) jj, zind(jj)
-!i    END DO
-
       IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
          l_jeq = .TRUE.
@@ -118 +113 @@
       tmv(:,:) = 0.0e0 ! CGrid EVP
 !!i
-      
-      
       ! metric coefficients for sea ice dynamic
       !----------------------------------------
@@ -169 +162 @@
                &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
 
-! better written but change the last digit and thus solver in less than 100 timestep
-!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   &
-!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1) 
-
-!           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   &
-!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
-
-!!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
             zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
             zusden2 = zusden * 2.0 
@@ -231 +216 @@
       tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
 
-!i here we can use umask with a i and j shift of -1,-1
       tmu(:,1) = 0.e0
       tmu(1,:) = 0.e0
       tmv(:,1) = 0.e0
       tmv(1,:) = 0.e0
-!     DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask
-!        DO ji = 2, jpim1
-! bug
+
       DO jj = 1, jpj - 1
          DO ji = 2 , jpi - 1

trunk/NEMO/LIM_SRC_3/limrhg.F90

@@ -2 +2 @@
    !!======================================================================
    !!                     ***  MODULE  limrhg  ***
-   !!   Ice rheology :  performs sea ice rheology
+   !!   Ice rheology : sea ice rheology
    !!======================================================================
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_rhg   : computes ice velocities
@@ -36 +36 @@
       rone    = 1.e0
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrhg.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
@@ -44 +44 @@
 
    SUBROUTINE lim_rhg( k_j1, k_jpj )
+
       !!-------------------------------------------------------------------
-      !!                 ***  SUBROUTINR lim_rhg  ***
-      !!
-      !! ** purpose :   determines sea ice drift from wind stress, ice-ocean
+      !!                 ***  SUBROUTINE lim_rhg  ***
+      !!                          EVP-C-grid
+      !!
+      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
       !!  stress and sea-surface slope. Ice-ice interaction is described by 
-      !!  a non-linear elasto-viscous-plastic law including shear strength 
-      !!  and a bulk rheology (Hunke and Dukowicz, 2002).   
-      !!
-      !!  Grid B:
-      !!  the routine calculates 4 estimates of the divergence, tension, 
-      !!  and shear on each corner of a given scalar grid box. Likewise,
-      !!  four estimates of the components of the stress tensor are 
-      !!  calculated on each corner. The internal forces on a corner are 
-      !!  calculated as averages of the four stress tensor contributions.
-      !!
-      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity
-      !!
-      !! ** NOTE : for the ice dynamics, the labeling convention is 
-      !!  that the velocity point (i,j) is located on the southwest corner
-      !!  of a scalar (i,j) grid box. This choice is somewhat unfortunate,
-      !!  as most models locate the velocity point (i,j) on the northeast
-      !!  corner...
-      !!
-      !! ** MORE IMPORTANT NOTES : related to the note above, it is crucial 
-      !!  to make sure that all variables and coefficients are located on
-      !!  right points of the grid. Verify everywhere!
+      !!  a non-linear elasto-viscous-plastic (EVP) law including shear 
+      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
+      !!
+      !!  The points in the C-grid look like this, dear reader
+      !!
+      !!                              (ji,jj)
+      !!                                 |
+      !!                                 |
+      !!                      (ji-1,jj)  |  (ji,jj)
+      !!                             ---------   
+      !!                            |         |
+      !!                            | (ji,jj) |------(ji,jj)
+      !!                            |         |
+      !!                             ---------   
+      !!                     (ji-1,jj-1)     (ji,jj-1)
+      !!
+      !! ** Inputs  : - wind forcing (stress), oceanic currents
+      !!                ice total volume (vt_i) per unit area
+      !!                snow total volume (vt_s) per unit area
+      !!
+      !! ** Action  : - compute u_ice, v_ice : the components of the 
+      !!                sea-ice velocity vector
+      !!              - compute delta_i, shear_i, divu_i, which are inputs
+      !!                of the ice thickness distribution
+      !!
+      !! ** Steps   : 1) Compute ice snow mass, ice strength 
+      !!              2) Compute wind, oceanic stresses, mass terms and
+      !!                 coriolis terms of the momentum equation
+      !!              3) Solve the momentum equation (iterative procedure)
+      !!              4) Prevent high velocities if the ice is thin
+      !!              5) Recompute invariants of the strain rate tensor
+      !!                 which are inputs of the ITD, store stress
+      !!                 for the next time step
+      !!              6) Control prints of residual (convergence)
+      !!                 and charge ellipse.
+      !!                 The user should make sure that the parameters
+      !!                 nevp, telast and creepl maintain stress state
+      !!                 on the charge ellipse for plastic flow
+      !!                 e.g. in the Canadian Archipelago
+      !!
+      !! ** References : Hunke and Dukowicz, JPO97
+      !!                 Bouillon et al., 08, in prep (update this when
+      !!                 published)
+      !!                 Vancoppenolle et al., OM08
       !!
       !! History :
-      !!   1.0  !  07-03  (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) 
-      !!                  EVP, C grid, LIM3
+      !!   1.0  !  07-03  (M.A. Morales Maqueda, S. Bouillon)
+      !!   2.0  !  08-03  M. Vancoppenolle : LIM3
       !!
       !!-------------------------------------------------------------------
@@ -79 +104 @@
       !
       INTEGER, INTENT(in) :: &
-         k_j1 ,  &            !: southern j-index for ice computation
-         k_jpj                !: northern j-index for ice computation
+         k_j1 ,                      & !: southern j-index for ice computation
+         k_jpj                         !: northern j-index for ice computation
 
       ! * Local variables
-      INTEGER ::   ji, jj, jl !: dummy loop indices
+      INTEGER ::   ji, jj              !: dummy loop indices
 
       INTEGER  :: &
-         iim1, ijm1, iip1 , ijp1   , & ! temporary integers
-         iter, jter                    !    "          "
+         iter, jter                    !: temporary integers
 
       CHARACTER (len=50) ::   charout
 
       REAL(wp) :: &
-         zt11  , zt12  , zt21  , zt22  ,   &  ! temporary scalars
-         ztagnx, ztagny , delta               !    "         "
+         zt11, zt12, zt21, zt22,     & !: temporary scalars
+         ztagnx, ztagny,             & !: wind stress on U/V points                       
+         delta                         !
 
       REAL(wp) :: &
-         za, zstms, zsang, zmask
-
-      REAL(wp),DIMENSION(jpi,jpj) ::   &
-         zpresh, zpreshc, zfrld1, zfrld2, zmass1, zmass2, zcorl1, zcorl2, &
-         za1ct, za2ct, zc1, zusw ,                      &
-         u_oce1, v_oce1, u_oce2, v_oce2, u_ice2, v_ice1
+         za,                         & !:
+         zstms,                      & !: temporary scalar for ice strength
+         zsang,                      & !: temporary scalar for coriolis term
+         zmask                         !: mask for the computation of ice mass
+
+      REAL(wp),DIMENSION(jpi,jpj) :: &
+         zpresh        ,             & !: temporary array for ice strength
+         zpreshc       ,             & !: Ice strength on grid cell corners (zpreshc)
+         zfrld1, zfrld2,             & !: lead fraction on U/V points                                    
+         zmass1, zmass2,             & !: ice/snow mass on U/V points                                    
+         zcorl1, zcorl2,             & !: coriolis parameter on U/V points
+         za1ct, za2ct  ,             & !: temporary arrays
+         zc1           ,             & !: ice mass
+         zusw          ,             & !: temporary weight for the computation
+                                       !: of ice strength
+         u_oce1, v_oce1,             & !: ocean u/v component on U points                           
+         u_oce2, v_oce2,             & !: ocean u/v component on V points
+         u_ice2,                     & !: ice u component on V point
+         v_ice1                        !: ice v component on U point
 
       REAL(wp) :: &
-         dtevp,dtotel,ecc2,z0,zr,zcca,zccb,denr, &
-         zu_ice2,zv_ice1, zddc, zdtc, zdst, zdsshx, zdsshy, &
-         sigma1, sigma2
-
-      REAL(wp),DIMENSION(jpi,jpj) ::   &
-         zf1, zf2
-
-      REAL(wp),DIMENSION(jpi,jpj) ::   &
-         zdd,                & !: Divergence [
-         zdt,                & !:
-         zds,                & !:
-         deltat,             &
-         deltac,             &
-         zs1,                &
-         zs2,                &
-         zs12
-
-      REAL :: &
-         zumax                 !: Maximal tolerated ice velocity
-
-      REAL(wp) ::            &
-         zresm                 !: Maximal error on ice velocity
-
-      REAL(wp),DIMENSION(jpi,jpj) ::   &
-         zu_ice           ,  & !: Ice velocity on previous time step
-         zv_ice           ,  &
-         zresr                 !: Local error on velocity
+         dtevp,                      & !: time step for subcycling
+         dtotel,                     & !:
+         ecc2,                       & !: square of yield ellipse eccenticity
+         z0,                         & !: temporary scalar
+         zr,                         & !: temporary scalar
+         zcca, zccb,                 & !: temporary scalars
+         zu_ice2,                    & !: 
+         zv_ice1,                    & !:
+         zddc, zdtc,                 & !: temporary array for delta on corners
+         zdst,                       & !: temporary array for delta on centre
+         zdsshx, zdsshy,             & !: term for the gradient of ocean surface
+         sigma1, sigma2                !: internal ice stress
+
+      REAL(wp),DIMENSION(jpi,jpj) :: &
+         zf1, zf2                      !: arrays for internal stresses
+
+      REAL(wp),DIMENSION(jpi,jpj) :: &
+         zdd, zdt,                   & !: Divergence and tension at centre of grid cells
+         zds,                        & !: Shear on northeast corner of grid cells
+         deltat,                     & !: Delta at centre of grid cells
+         deltac,                     & !: Delta on corners
+         zs1, zs2,                   & !: Diagonal stress tensor components zs1 and zs2 
+         zs12                          !: Non-diagonal stress tensor component zs12
 
       REAL(wp) :: &
-         zindb            ,  & !: ice (1) or not (0)      
-         zdummy                !: ice computation
+         zresm            ,          & !: Maximal error on ice velocity
+         zindb            ,          & !: ice (1) or not (0)      
+         zdummy                        !: dummy argument
+
+      REAL(wp),DIMENSION(jpi,jpj) :: &
+         zu_ice           ,          & !: Ice velocity on previous time step
+         zv_ice           ,          &
+         zresr                         !: Local error on velocity
 
       INTEGER :: &
          stress_mean_swi = 1
-
-      !!----------------------------------------------------------------------
-         
-!
-!------------------------------------------------------------------------------!
-! 0) Ice-Snow mass (zc1), ice strength (zpresh)                                !
-!------------------------------------------------------------------------------!
-!
-      ! Maximal tolerated velocity
-!     zumax = 1.0
-
+!
+!------------------------------------------------------------------------------!
+! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                !
+!------------------------------------------------------------------------------!
+!
       ! Put every vector to 0
       zpresh(:,:)  = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:)  = 0.0
@@ -154 +188 @@
       za1ct(:,:)   = 0.0 ; za2ct(:,:)   = 0.0 
       zc1(:,:)     = 0.0 ; zusw(:,:)    = 0.0
-
       u_oce1(:,:)  = 0.0 ; v_oce1(:,:)  = 0.0 ; u_oce2(:,:)  = 0.0
       v_oce2(:,:)  = 0.0 ; u_ice2(:,:)  = 0.0 ; v_ice1(:,:)  = 0.0
-
-      zf1(:,:) = 0.0     ; zf2(:,:) = 0.0
-
+      zf1(:,:)     = 0.0 ; zf2(:,:) = 0.0
       zdd(:,:)     = 0.0 ; zdt(:,:)     = 0.0 ; zds(:,:)     = 0.0
       deltat(:,:)  = 0.0 ; deltac(:,:)  = 0.0 
-
-!
-!------------------------------------------------------------------------------!
-! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                !
-!------------------------------------------------------------------------------!
-!
-      ! compute ice strength
+      zpresh(:,:)  = 0.0 
+
+      ! Ice strength on T-points
       CALL lim_itd_me_icestrength(ridge_scheme_swi)
 
-      zpresh(:,:) = 0.0 
-
+      ! Ice mass and temp variables
       DO jj = k_j1 , k_jpj-1
          DO ji = 1 , jpi
@@ -178 +204 @@
             zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) / 2.
             ! tmi = 1 where ther is ice or on land
+            ! Claude sees a bug, next line : tmi(ji,jj)
             tmi           = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - &
                                     epsd ) ) ) * tms(ji,jj)
@@ -186 +213 @@
       CALL lbc_lnk( zpresh(:,:), 'T',  1. )
 
-      ! Ice strength (zpreshc) on grid cell corners (needed for    
-      ! calculation of shear stress 
+! Claude sees a bug
+!     CALL lbc_lnk( tmi(:,:), 'T',  1. )
+
+      ! Ice strength on grid cell corners (zpreshc)
+      ! needed for calculation of shear stress 
       DO jj = k_j1+1, k_jpj-1
          DO ji = 2, jpim1
-             zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2)           &
-                &            + tms(ji+1,jj  ) * wght(ji+1,jj+1,2,1) + tms(ji,jj  ) * wght(ji+1,jj+1,1,1)
+             zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
+                &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
+                &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
+                &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
              zusw(ji,jj)    = 1.0 / MAX( zstms, epsd )
-             zpreshc(ji,jj)  = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2)  &
-                &              + zpresh(ji+1,jj   ) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj  ) * wght(ji+1,jj+1,1,1) ) &
-                &              * zusw(ji,jj)
+             zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
+                &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
+                &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
+                &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
+                &             ) * zusw(ji,jj)
          END DO
       END DO
@@ -233 +267 @@
 
             ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = (zt12*zc1(ji,jj)+zt11*zc1(ji+1,jj))/(zt11+zt12+epsd)
-            zmass2(ji,jj) = (zt22*zc1(ji,jj)+zt21*zc1(ji,jj+1))/(zt21+zt22+epsd)
-
-            zcorl1(ji,jj)  = zmass1(ji,jj)*(e1t(ji+1,jj)*fcor(ji,jj)+e1t(ji,jj)*fcor(ji+1,jj))/(e1t(ji,jj)+e1t(ji+1,jj)+epsd)
-            zcorl2(ji,jj)  = zmass2(ji,jj)*(e2t(ji,jj+1)*fcor(ji,jj)+e2t(ji,jj)*fcor(ji,jj+1))/(e2t(ji,jj+1)+e2t(ji,jj)+epsd)
-            !
-            ! Reminder: since this is a C grid, the location of ocean currents 
-            ! calculated by OPA should coincide with ice model vector points:
-            ! no need for interpolation once the definition of u_io and v_io
-            ! has been changed in module "icestp". Arrays u_oce1 and v_oce2 could
-            ! then be replaced by u_oce and v_oce, respectively.
+            zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)
+            zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)
+            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + &
+                                              e1t(ji,jj)*fcor(ji+1,jj) ) &
+                                           / (e1t(ji,jj) + e1t(ji+1,jj) + epsd )
+            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + &
+                                              e2t(ji,jj)*fcor(ji,jj+1) ) &
+                                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
             !
             u_oce1(ji,jj)  = u_io(ji,jj)
-
-            ! SB modif because ocean has no slip boundary condition
+            v_oce2(ji,jj)  = v_io(ji,jj)
+
+            ! Ocean has no slip boundary condition
             v_oce1(ji,jj)  = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj)    &
                 &                 +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) &
                 &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)  
 
-            u_oce2(ji,jj)  = 0.5*( (u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)    &
+            u_oce2(ji,jj)  = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)     &
                 &                 +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) &
-                &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-    
-            v_oce2(ji,jj)  = v_io(ji,jj)
+                &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
 
             ! Wind stress.
@@ -262 +292 @@
             ztagny = ( 1. - zfrld2(ji,jj) ) * gtauy(ji,jj)
 
-            ! Computation of the velocity field taking into account the ice-ice interaction.                                 
+            ! Computation of the velocity field taking into account the ice internal interaction.
             ! Terms that are independent of the velocity field.
-            ! Reminder: does zcorl*(-v_oce,u_oce) represent the surface pressure gradient? Why is it still
-            !           calculated in this way? An ocean model with a free surface should provide the gradient
-            !           of surface elevation directly...
 
             ! SB On utilise maintenant le gradient de la pente de l'ocean
@@ -272 +299 @@
             !    zdsshx =  (ssh_io(ji+1,jj) - ssh_io(ji,jj))/e1u(ji,jj)
             !    zdsshy =  (ssh_io(ji,jj+1) - ssh_io(ji,jj))/e2v(ji,jj)  
+
             zdsshx = 0.0
             zdsshy = 0.0 
@@ -283 +311 @@
 !
 !------------------------------------------------------------------------------!
-! 3) Solution of the momentum equation
+! 3) Solution of the momentum equation, iterative procedure
 !------------------------------------------------------------------------------!
 !
@@ -319 +347 @@
          END DO      
 
-           DO jj = k_j1+1, k_jpj-1
-              DO ji = 2, jpim1
+         DO jj = k_j1+1, k_jpj-1
+            DO ji = 2, jpim1
 
          !  
@@ -327 +355 @@
          !- zds(:,:): shear on northeast corner of grid cells
                  !
-                 !- IMPORTANT REMINDER: note that, the way these terms are coded, there are many repeated 
-                 !                      calculations. Speed could be improved by regrouping terms. For
+                 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 
+                 !                      there are many repeated calculations. 
+                 !                      Speed could be improved by regrouping terms. For
                  !                      the moment, however, the stress is on clarity of coding to avoid
-                 !                      bugs (mamm).
-                 !- ALSO: arrays zdd, zdt, zds and delta could be removed in the future to minimise memory demand.
+                 !                      bugs (Martin, for Miguel).
+                 !
+                 !- ALSO: arrays zdd, zdt, zds and delta could 
+                 !  be removed in the future to minimise memory demand.
+                 !
                  !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
                  !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
@@ -337 +369 @@
                  !
                  !
-                 !          (ji,jj)
-                 !             |
-                 !             |
-                 !  (ji-1,jj)  |  (ji,jj)
-                 !         ---------   
-                 !        |         |
-                 !        | (ji,jj) |------(ji,jj)
-                 !        |         |
-                 !         ---------   
-                 !(ji-1,jj-1)     (ji,jj-1)
-                 ! 
-
                  zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
                     &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
@@ -366 +386 @@
 
                  !
-                 ! SB modif because ocean has no slip boundary condition 
                  zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
                     &            -u_ice(ji,jj)/e1u(ji,jj)                    &
@@ -379 +398 @@
 
  
-                 ! SB modif because need to compute zddc and zdtc correctly 
                  v_ice1(ji,jj)  = 0.5*( (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)) &
@@ -406 +424 @@
                     &         / area(ji,jj)
 
-! Old lines
-!                deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 
-!    &                                 ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & 
-!    &                               ) + creepl
-! New lines
-! Regularisation of dP/dx
                  delta = SQRT( zdd(ji,jj)*zdd(ji,jj) +                        & 
      &                       ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
                  deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 +                    &
                                  (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
-! End of new lines
 
                  !-Calculate stress tensor components zs1 and zs2 
                  !-at centre of grid cells (see section 3.5 of CICE user's guide).
-! Old lines
-!                zs1(ji,jj) = ( zs1(ji,jj)       &
-!                   &          - dtotel*((1.0-alphaevp)*zs1(ji,jj)+(1.0-zdd(ji,jj)/deltat(ji,jj))*zpresh(ji,jj)))       &
-!                   &        /(1.0+alphaevp*dtotel)
-! New lines
                  zs1(ji,jj) = ( zs1(ji,jj) &
                     &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
-                    &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) &       
+                    &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
+                                 * zpresh(ji,jj) ) )                          &       
                     &        / ( 1.0 + alphaevp * dtotel )
-! End of new lines
-
-                 zs2(ji,jj) = ( zs2(ji,jj)                                                                              &
-                    &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj)-zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)))        &
-                    &        /(1.0+alphaevp*ecc2*dtotel)
+
+                 zs2(ji,jj) = ( zs2(ji,jj)   &
+                    &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
+                                 zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
+                    &        / ( 1.0 + alphaevp*ecc2*dtotel )
 
               END DO
@@ -443 +450 @@
            DO jj = k_j1+1, k_jpj-1
               DO ji = 2, jpim1
-
                  !- Calculate Delta on corners
-                 
                  zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
                     &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
@@ -455 +460 @@
                     &        / ( e1f(ji,jj) * e2f(ji,jj) )
 
-
                  zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
                     &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
@@ -465 +469 @@
                     &        / ( e1f(ji,jj) * e2f(ji,jj) )
 
-                 deltac(ji,jj) = sqrt(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
+                 deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
 
                  !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-
-                 zs12(ji,jj) = ( zs12(ji,jj)                                                                            &
-                    &           -dtotel*((1.0-alphaevp)*ecc2*zs12(ji,jj)-zds(ji,jj)/(2.0*deltac(ji,jj))*zpreshc(ji,jj)))&
-                    &         /(1.0+alphaevp*ecc2*dtotel)
-
-              END DO
-           END DO
+                 zs12(ji,jj) = ( zs12(ji,jj)      &
+                    &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
+                    &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
+                    &         / ( 1.0 + alphaevp*ecc2*dtotel )
+
+              END DO ! ji
+           END DO ! jj
 
            CALL lbc_lnk( zs12(:,:), 'F', 1. )
@@ -481 +485 @@
            DO jj = k_j1+1, k_jpj-1
               DO ji = 2, jpim1
-
-     !
-     !- contribution of zs1, zs2 and zs12 to zf1
-     !
-     !          (ji,jj)
-     !             |
-     !             |
-     !  (ji-1,jj)  |  (ji,jj)
-     !         ---------   
-     !        |         |
-     !        | (ji,jj) |------(ji,jj)
-     !        |         |
-     !         ---------   
-     !(ji-1,jj-1)     (ji,jj-1)
-     ! 
-
-                zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)                                         &
-                   &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)           &
-                   &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)     &
-                   &             )                                                                              &
-                   &            /(e1u(ji,jj)*e2u(ji,jj))
-
-     ! 
-     !contribution of zs1, zs2 and zs12 to zf2
-     !
-     !          (ji,jj)
-     !             |
-     !             |
-     !  (ji-1,jj)  |  (ji,jj)
-     !         ---------   
-     !        |         |
-     !        | (ji,jj) |------(ji,jj)
-     !        |         |
-     !         ---------   
-     !(ji-1,jj-1)     (ji,jj-1)
-     !
-
-                zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj)                                         &
-                   &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2-zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj)           &
-                   &              +2.0*(zs12(ji,jj)*e2f(ji,jj)**2-zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj)     &
-                   &             )                                                                              &
-                   &            /(e1v(ji,jj)*e2v(ji,jj))
-
+                !- contribution of zs1, zs2 and zs12 to zf1
+                zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
+                   &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
+                   &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
+                   &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
+                ! contribution of zs1, zs2 and zs12 to zf2
+                zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
+                   &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
+                   &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
+                                    zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
+                   &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
               END DO
            END DO
@@ -531 +503 @@
          ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
          !
-           IF (mod(iter,2).eq.0) THEN 
+           IF (MOD(iter,2).eq.0) THEN 
 
               DO jj = k_j1+1, k_jpj-1
                  DO ji = 2, jpim1
-         !
-         !          (ji,jj)
-         !             |
-         !             |
-         !  (ji-1,jj)  |  (ji,jj)
-         !         ---------   
-         !        |         |
-         !        | (ji,jj) |------(ji,jj)
-         !        |         |
-         !         ---------   
-         !(ji-1,jj-1)     (ji,jj-1)
-         !
-                    zmask        = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj)
+                    zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
                     zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
                     z0           = zmass1(ji,jj)/dtevp
@@ -555 +515 @@
                       &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   &
                       &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-                    za           = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
-                    zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
+                    za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
+                                              (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
+                    zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
+                                   za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
                     zcca         = z0+za*cangvg
                     zccb         = zcorl1(ji,jj)+za*zsang
                     u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
-!                   u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) )
 
                  END DO
@@ -569 +530 @@
               DO jj = k_j1+1, k_jpj-1
                  DO ji = 2, jpim1
-         !
-         !          (ji,jj)
-         !             |
-         !             |
-         !  (ji-1,jj)  |  (ji,jj)
-         !         ---------   
-         !        |         |
-         !        | (ji,jj) |------(ji,jj)
-         !        |         |
-         !         ---------   
-         !(ji-1,jj-1)     (ji,jj-1)
-         !
+
                     zmask        = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj)
                     zsang        = sign(1.0,fcor(ji,jj))*sangvg
@@ -588 +538 @@
                 &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))   &
                 &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-                    za           = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                    zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
+                    za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
+                                   (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
+                    zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
+                                   za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
                     zcca         = z0+za*cangvg
                     zccb         = zcorl2(ji,jj)+za*zsang
                     v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
-!                   v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )
 
                  END DO
@@ -600 +551 @@
               CALL lbc_lnk( v_ice(:,:), 'V', -1. )
 
-           ELSE 
+         ELSE 
               DO jj = k_j1+1, k_jpj-1
                  DO ji = 2, jpim1
-         !
-         !          (ji,jj)
-         !             |
-         !             |
-         !  (ji-1,jj)  |  (ji,jj)
-         !         ---------   
-         !        |         |
-         !        | (ji,jj) |------(ji,jj)
-         !        |         |
-         !         ---------   
-         !(ji-1,jj-1)     (ji,jj-1)
-         !
                     zmask        = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj)
                     zsang        = sign(1.0,fcor(ji,jj))*sangvg
@@ -623 +562 @@
                        &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
 
-                    za           = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                    zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
+                    za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
+                                   (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
+                    zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
+                                   za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
                     zcca         = z0+za*cangvg
                     zccb         = zcorl2(ji,jj)+za*zsang
                     v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
-!                   v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )
 
                  END DO
@@ -637 +577 @@
               DO jj = k_j1+1, k_jpj-1
                  DO ji = 2, jpim1
-         !
-         !          (ji,jj)
-         !             |
-         !             |
-         !  (ji-1,jj)  |  (ji,jj)
-         !         ---------   
-         !        |         |
-         !        | (ji,jj) |------(ji,jj)
-         !        |         |
-         !         ---------   
-         ! (ji-1,jj-1)     (ji,jj-1)
-         !
-
-                    zmask        = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj)
-                    zsang        = sign(1.0,fcor(ji,jj))*sangvg
+                    zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
+                    zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
                     z0           = zmass1(ji,jj)/dtevp
                     ! SB modif because ocean has no slip boundary condition
@@ -658 +585 @@
                        &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
     
-                    za           = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
-                    zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
+                    za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
+                                   (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
+                    zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
+                                   za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
                     zcca         = z0+za*cangvg
                     zccb         = zcorl1(ji,jj)+za*zsang
                     u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
-!                   u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) )
-
-                 END DO
-              END DO
+                 END DO ! ji
+              END DO ! jj
 
               CALL lbc_lnk( u_ice(:,:), 'U', -1. )
 
-           ENDIF
-
-           !---  Convergence test.
-           DO jj = k_j1+1 , k_jpj-1
-              zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
-                            ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-           END DO
-           zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
-
-           !-------------------------------------------
-           ! Compute dynamical inputs of the itd model
-           !-------------------------------------------
-           ! Mean over all iterations
-
-      !    IF ( stress_mean_swi .EQ. 1 ) THEN
-      !       DO jj = k_j1+1, k_jpj-1
-      !          DO ji = 2, jpim1
-      !             divu_i(ji,jj) = divu_i(ji,jj) + zdd(ji,jj) / nevp
-      !             delta_i(ji,jj) = delta_i(ji,jj) + deltat (ji,jj) / nevp 
-      !             shear_i(ji,jj) = shear_i(ji,jj) + zds (ji,jj) / nevp 
-      !          END DO
-      !       END DO
-      !    ENDIF
+      ENDIF 
+
+      !---  Convergence test.
+      DO jj = k_j1+1 , k_jpj-1
+         zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
+                       ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+      END DO
+      zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
 
       !                                                   ! ==================== !
@@ -698 +610 @@
       !                                                   ! ==================== !
 
-      !-------------------------
-      ! Prevent high velocities
-      !-------------------------
+!
+!------------------------------------------------------------------------------!
+! 4) Prevent ice velocities when the ice is thin
+!------------------------------------------------------------------------------!
+!
       ! If the ice thickness is below 1cm then ice velocity should equal the
-      ! ocean velocity
+      ! ocean velocity, 
+      ! This prevents high velocity when ice is thin
       DO jj = k_j1+1, k_jpj-1
          DO ji = 2, jpim1
@@ -713 +628 @@
          END DO
       END DO
-
+!
+!------------------------------------------------------------------------------!
+! 5) Compute stress rain invariants and store strain tensor
+!------------------------------------------------------------------------------!
+!
+      ! Compute delta, shear and div, inputs for mechanical redistribution 
       DO jj = k_j1+1, k_jpj-1
          DO ji = 2, jpim1
-            ! Recompute stress tensor invariants
             !- zdd(:,:), zdt(:,:): divergence and tension at centre 
-            !  of grid cells
             !- zds(:,:): shear on northeast corner of grid cells
-
             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
             zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
@@ -782 +699 @@
       END DO !ji
 
-      ! dynamical invariants
+! MV This should be removed as it is the same as previous lines
+!     ! dynamical invariants
 !     IF ( stress_mean_swi .EQ. 0 ) THEN
+! the following should not be necessary
          DO jj = k_j1+1, k_jpj-1
             DO ji = 2, jpim1
@@ -797 +716 @@
       CALL lbc_lnk( shear_i(:,:), 'F', 1. )
 
-      IF(ln_ctl) THEN
-         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
-         CALL prt_ctl_info(charout)
-         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
-      ENDIF
-
-      ! Stress tensor
+
+      ! Store the stress tensor for next time step
+      ! this accelerates convergence and improves stability
       IF ( stress_mean_swi .EQ. 1 ) THEN
          DO jj = k_j1+1, k_jpj-1
@@ -814 +729 @@
       ENDIF
 
-      !Ice internal stresses
+      ! Lateral boundary condition
       CALL lbc_lnk( stress1_i(:,:), 'T', 1. )
       CALL lbc_lnk( stress2_i(:,:), 'T', 1. )
       CALL lbc_lnk( stress12_i(:,:), 'F', 1. )
  
-      !------------------------
-      ! Write charge ellipse
-      !------------------------
-
+!
+!------------------------------------------------------------------------------!
+! 6) Control prints of residual and charge ellipse
+!------------------------------------------------------------------------------!
+!
+      ! print the residual for convergence
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, iter
+         CALL prt_ctl_info(charout)
+         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
+      ENDIF
+
+      ! print charge ellipse
+      ! This can be desactivated once the user is sure that the stress state
+      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
       IF(ln_ctl) THEN
          CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)

trunk/NEMO/LIM_SRC_3/limrst.F90

@@ -6 +6 @@
 #if defined key_lim3
 !!----------------------------------------------------------------------
-!!   'key_lim3' :                                   LIM sea-ice model
+!!   'key_lim3' :                                    LIM3 sea-ice model
 !!----------------------------------------------------------------------
 !!   lim_rst_write   : write of the restart file 
@@ -19 +19 @@
 USE ice_oce         ! ice variables
 USE daymod
-USE limicepoints
 !USE limvar
 
@@ -65 +64 @@
    llbon
    INTEGER :: &
-   ji, jj, jk ,jl, jm, index
+   ji, jj, jk , jl
    INTEGER :: &
    inumwrs, it0, itime

trunk/NEMO/LIM_SRC_3/limtab.F90

@@ -5 +5 @@
    !!======================================================================
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   tab_2d_1d  : 2-D to 1-D
@@ -20 +22 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008) 
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limtab.F90,v 1.2 2005/03/27 18:34:42 opalod Exp $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 

trunk/NEMO/LIM_SRC_3/limthd.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_thd      : thermodynamic of sea ice
@@ -22 +22 @@
    USE dom_ice         ! LIM sea-ice domain
    USE iceini
-   USE limicepoints
    USE limthd_dif
    USE limthd_dh
@@ -84 +83 @@
       !! * Local variables
       INTEGER  ::  ji, jj, jk, jl,  &
-                   zjid , zjjd ,    &
                    zji  , zjj,      &   ! dummy loop indices
                    nbpb ,           &   ! nb of icy pts for thermo. cal.
@@ -95 +93 @@
       REAL(wp) ::   &
          zinda              ,  &   ! switch for test. the val. of concen.
-         zindb, zindg       ,  &   ! switches for test. the val of arg
-         za , zh, zthsnice  ,  &
+         zindb,                &   ! switches for test. the val of arg
+         zthsnice           ,  &
          zfric_u            ,  &   ! friction velocity 
          zfnsol             ,  &   ! total non solar heat
@@ -104 +102 @@
 
       REAL(wp) ::   &
-         ztmelts            ,  &   ! ice layers melting point
-         zaaa               ,  &   ! numeric value for 2nd degree polynomial
-         zbbb               ,  &   !    "      "    "   "    "        "
-         zccc               ,  &   !    "      "    "   "    "        "
-         zdiscrim           ,  &   !    "      "    "   "    "        "
-         zq                 ,  &   ! sea ice temporary energy of melting
-         zqold              ,  &   ! test variable for ice heat content
          zareamin
          
@@ -224 +215 @@
 
 ! still need to be updated : fdtcn !!!!
-!           !-- Lead heat budget (part 1, next one is in limthd_dh, then limthd_lab
+!           !-- Lead heat budget (part 1, next one is in limthd_dh
 !           !-- qldif -- (or qldif_1d in 1d routines)
             zfontn         = sprecip(ji,jj) * lfus              !   energy of melting
@@ -510 +501 @@
       REAL(wp) :: &
          zdes,    &          ! snow heat content increment (dummy)
-         ztmelts, &          ! layer melting point (dummy)
          zeps                ! very small value (1.e-10)
 
@@ -560 +550 @@
          jl                            !: category number
 
-      REAL(wp), DIMENSION(jpl) ::   &  !: ! goes to trash
-         zdummy
-
       REAL(wp)                 ::   &  !: ! goes to trash
          meance,                    &  !: mean conservation error
@@ -573 +560 @@
       INTEGER  :: &
          ji,jj,jk,                  &  !: loop indices
-         zji, zjj, zjiindex
+         zji, zjj
       !!---------------------------------------------------------------------
 
@@ -745 +732 @@
          jl                            !: category number
 
-      REAL(wp), DIMENSION(jpl) ::   &  !: ! goes to trash
-         zdummy
-
       REAL(wp)                 ::   &  !: ! goes to trash
          meance,                    &  !: mean conservation error
-         max_cons_err,              &  !: maximum tolerated conservation error
-         max_surf_err                  !: maximum tolerated surface error
+         max_cons_err                  !: maximum tolerated conservation error
 
       INTEGER ::                    &
@@ -758 +741 @@
       INTEGER  :: &
          ji,jj,jk,                  &  !: loop indices
-         zji, zjj, zjiindex
+         zji, zjj
 
       !!---------------------------------------------------------------------
@@ -974 +957 @@
       ENDIF
             
-      uscomi = 1.0 / ( 1.0 - amax )   ! inverse of minimum lead fraction
       rcdsn = hakdif * rcdsn 
       rcdic = hakdif * rcdic

trunk/NEMO/LIM_SRC_3/limthd_dh.F90

@@ -1 +1 @@
 MODULE limthd_dh
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_dh ***
@@ -19 +22 @@
    USE ice
    USE par_ice
-   USE limicepoints
       
    IMPLICIT NONE
@@ -36 +38 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2005)
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !!----------------------------------------------------------------------
 
@@ -70 +72 @@
        !! ** References : Bitz and Lipscomb, JGR 99
        !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
-       !!                 Vancoppenolle Fichefet and Bitz, GRL 2005
+       !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
+       !!                 Vancoppenolle et al., OM08
        !!
        !! ** History  : 
        !!   original code    01-04 (LIM)
        !!   original routine
-       !!               (05-2003) Martin Vancoppenolle, Louvain-La-Neuve, Belgium
-       !!               (06/07-2005) Martin Vancoppenolle for the 3D version
+       !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
+       !!               (06/07-2005) 3D version
+       !!               (03-2008)    Clean code
        !!
        !!------------------------------------------------------------------
@@ -90 +94 @@
           jk        ,         &  !: ice layer index
           isnow     ,         &  !: switch for presence (1) or absence (0) of snow
-          index     ,         &  !: ice point index (limicepoints)
           zji       ,         &  !: 2D corresponding indices to ji
           zjj       ,         &
@@ -107 +110 @@
           zhnfi     ,         &
           zihg      ,         &
-          zdhgm     ,         &
           zdhnm     ,         &
           zhnnew    ,         &
-          zdfseqv   ,         &
           zeps = 1.0e-13,     &
           zhisn     ,         &
@@ -116 +117 @@
                                  !: entrapment
           zds       ,         &  !: increment of bottom ice salinity
-          zoldsinew ,         &  !: dummy argument for error diagnosis
           zcoeff    ,         &  !: dummy argument for snowfall partitioning
                                  !: over ice and leads
-          zjiindex  ,         &  !: zdhs_temp, &
           zsm_snowice,        &  !: snow-ice salinity
           zswi1     ,         &  !: switch for computation of bottom salinity
@@ -159 +158 @@
 
        REAL(wp), DIMENSION(jpij,jkmax) :: &
-          zqt_i_lay  ,        &  !: total ice heat content
-          zqt_s_lay              !: total snow heat content
+          zqt_i_lay              !: total ice heat content
 
        ! Heat conservation 
@@ -181 +179 @@
       WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation'
       WRITE(numout,*) '~~~~~~~~~'
+
       zfsalt_melt(:)  = 0.0
       ftotal_fin(:)   = 0.0
@@ -211 +210 @@
       dsm_i_se_1d(:) = 0.0     
       dsm_i_si_1d(:) = 0.0     
-!     ! heat conservation
-!     IF ( con_i ) THEN
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' zf_surf   : ', z_f_surf(jiindex_1d)
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' qnsr_ice  : ', qnsr_ice_1d(jiindex_1d)
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' qsr_ice   : ', qsr_ice_1d(jiindex_1d)
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' fc_su     : ', fc_su(jiindex_1d)
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' i0        : ', i0(jiindex_1d)
-!     ENDIF
-
+!
 !------------------------------------------------------------------------------!
 !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
 !------------------------------------------------------------------------------!
-
-      !-----------------
+!
       ! Layer thickness
-      !-----------------
       DO ji = kideb,kiut
          zh_i(ji) = ht_i_b(ji) / nlay_i
@@ -232 +221 @@
       END DO
 
-      !--------------------------------------
-      ! Snow energy of melting, heat content
-      !--------------------------------------
+      ! Total enthalpy of the snow
       zqt_s(:) = 0.0
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            ! this line could be removed
-!           q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
-            ! heat conservation
             zqt_s(ji)    =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
          END DO
       END DO
 
-!     ! heat conservation
-!     qt_s_in(:,jl) = zqt_s(:)
-!     IF (jiindex_1d.GT.0) &
-!     WRITE(numout,*) 'jl : ', jl, ' zqts    : ', zqt_s(jiindex_1d) / &
-!                      rdt_ice
-
-!     IF (jiindex_1d.GT.0) THEN
-!     WRITE(numout,*) ' Vertical profile : '
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     WRITE(numout,*) ' qm0    : ', q_s_b(jiindex_1d,1:nlay_s) * ht_s_b(jiindex_1d) / &
-!     rdt_ice
-!     WRITE(numout,*) ' q_s_b  : ', q_s_b(jiindex_1d,1)
-!     WRITE(numout,*) ' t_s_b  : ', t_s_b(jiindex_1d,1)
-!     WRITE(numout,*) ' zthick0: ', ht_s_b(jiindex_1d)
-!     ENDIF
-
-      !-------------------------------------
-      ! Ice energy of melting, heat content
-      !-------------------------------------
+      ! Total enthalpy of the ice
       zqt_i(:) = 0.0
       DO jk = 1, nlay_i
          DO ji = kideb,kiut
-            ! Melting point in K
-!           ztmelts          =  - tmut * s_i_b(ji,jk) + rtt 
-            ! this line could be removed
-!           q_i_b(ji,jk)     =  rhoic *                                       &
-!                             ( cpic    * ( ztmelts - t_i_b(ji,jk) )          &
-!                             + lfus * ( 1.0 - (ztmelts-rtt) /                &
-!                               MIN( ( t_i_b(ji,jk) - rtt ) , - zeps ) )      &
-!                               - rcp      * ( ztmelts-rtt  ) )
             zqt_i(ji)        =  zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
             zqt_i_lay(ji,jk) =  q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
          END DO
       END DO
-!     IF (jiindex_1d.GT.0) &
-!     WRITE(numout,*) 'jl : ', jl, ' zqti    : ', zqt_i(jiindex_1d) / &
-!                      rdt_ice
-
-!     IF (jiindex_1d.GT.0) THEN
-!     WRITE(numout,*) ' --- Vertical profile : '
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' q_i_b  : ', q_i_b(jiindex_1d, 1:nlay_i)
-!     WRITE(numout,*) ' t_i_b  : ', t_i_b(jiindex_1d, 1:nlay_i)
-!     WRITE(numout,*) ' qm0    : ', zqt_i_lay(jiindex_1d,1:nlay_i) / rdt_ice
-!     WRITE(numout,*) ' zthick0: ', ht_i_b(jiindex_1d) / nlay_i
-!     ENDIF
-
 !
 !------------------------------------------------------------------------------|
@@ -298 +242 @@
 !------------------------------------------------------------------------------|
 !
-      !---------------------
-      ! Snow precips / melt
-      !---------------------
+      !-------------------------
+      ! 3.1 Snow precips / melt
+      !-------------------------
       ! Snow accumulation in one thermodynamic time step
       ! snowfall is partitionned between leads and ice
@@ -306 +250 @@
       ! but because of the winds, more snow falls on leads than on sea ice
       ! and a greater fraction (1-at_i)^beta of the total mass of snow 
-      ! (beta < 1) falls in leads
+      ! (beta < 1) falls in leads.
       ! In reality, beta depends on wind speed, 
       ! and should decrease with increasing wind speed but here, it is 
-      ! considered as a constant
-      ! an average value is 0.66
+      ! considered as a constant. an average value is 0.66
       ! Martin Vancoppenolle, December 2006
 
@@ -322 +265 @@
       ! Melt of fallen snow
       DO ji = kideb, kiut
-
          ! tatm_ice is now in K
          zqprec(ji)     =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )  
@@ -334 +276 @@
          qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
          zqt_s(ji)      =  zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji)
-         
-      END DO
-
-!     IF (jiindex_1d .GT. 0) THEN
-!        WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)
-!        WRITE(numout,*) ' zqprec    : ', zqprec(jiindex_1d)*( zdh_s_pre(jiindex_1d) + zdh_s_mel(jiindex_1d) ) / rdt_ice
-!        WRITE(numout,*) ' tatm      : ', tatm_ice_1d(jiindex_1d)
-!        WRITE(numout,*) 'jl : ', jl, ' zqts    : ', zqt_s(jiindex_1d) / &
-!                         rdt_ice
-!     ENDIF
-
-      ! updated total snow heat content
+      END DO
+
+      ! Update total snow heat content
       zqt_s(ji)         =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
-
-!     IF (jiindex_1D .GT. 0) &
-!     WRITE(numout,*) 'jl : ', jl, ' zqts    : ', zqt_s(jiindex_1d) / &
-!                      rdt_ice
 
       ! Snow melt due to surface heat imbalance
@@ -358 +287 @@
             zqfont_su(ji)  =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * &
                               q_s_b(ji,jk) 
-!           IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &
-!           zqfont_su(jiindex_1d)
             zdeltah(ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
             zdh_s_mel(ji)  =  zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow    
          END DO
       END DO
-
-!     !+++++
-!     IF (jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' dhs_pre :', zdh_s_pre(jiindex_1d)
-!     WRITE(numout,*) ' dhs_mel :', zdh_s_mel(jiindex_1d)
-!     WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)
-!     ENDIF
-!     !+++++
 
       ! Apply snow melt to snow depth
@@ -389 +308 @@
       END DO ! ji
 
-!     !+++++
-!     IF (jiindex_1d .GT. 1) WRITE(numout,*) ' dhs_tot   :', dh_s_tot(jiindex_1d)
-!     IF (jiindex_1d .GT. 1) WRITE(numout,*) ' zqfont_su :', zqfont_su(jiindex_1d)
-!     !+++++
-
-      !----------------------
-      ! Surface ice ablation 
-      !----------------------
+      !--------------------------
+      ! 3.2 Surface ice ablation 
+      !--------------------------
       DO ji = kideb, kiut 
-!        IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &
-!        zqfont_su(jiindex_1d)
          dh_i_surf(ji) =  0.0
-         ! Heat conservation test
+         ! For heat conservation test
          z_f_surf(ji)  =  zqfont_su(ji) / rdt_ice ! heat conservation test
          zdq_i(ji)     =  0.0
@@ -408 +320 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut 
+            ! melt of layer jk
             zdeltah(ji,jk)      = - zqfont_su(ji) / q_i_b(ji,jk)
+            ! recompute heat available
             zqfont_su(ji)       = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) *  &
                                   q_i_b(ji,jk) 
+            ! melt of layer jk cannot be higher than its thickness
             zdeltah(ji,jk)      = MAX( zdeltah(ji,jk) , - zh_i(ji) )
+            ! update surface melt
             dh_i_surf(ji)       = dh_i_surf(ji) + zdeltah(ji,jk) 
+            ! for energy conservation
             zdq_i(ji)           = zdq_i(ji) + zdeltah(ji,jk) *                &
                                         q_i_b(ji,jk) / rdt_ice
-            ! Contribution to salt flux -> to change
+            ! contribution to ice-ocean salt flux 
             zji                 = MOD( npb(ji) - 1, jpi ) + 1
             zjj                 = ( npb(ji) - 1 ) / jpi + 1
             zfsalt_melt(ji)     = zfsalt_melt(ji) +                           &
-!                                 ( sss_io(zji,zjj) - s_i_b(ji,jk) ) *        &
                                   ( sss_io(zji,zjj) - sm_i_b(ji)   ) *        &
                                   a_i_b(ji) *                                 &
                                   MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
-!+++++++++++
-!           IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN
-!              WRITE(numout,*) ' surf abl '
-!              WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) / ( sss_io(zji,zjj)+epsi16)
-!              WRITE(numout,*) ' sss_io      : ', sss_io(zji,zjj)
-!           ENDIF
-!+++++++++++
          END DO ! ji
       END DO ! jk
@@ -446 +355 @@
          ENDIF
 
-         IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!.GE. 1.0e-3 ) .AND. &
-!             ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN
+         IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
             WRITE(numout,*) ' ALERTE heat loss for surface melt '
             WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
@@ -471 +379 @@
       ENDIF ! con_i
 
-      !------------------
-      ! Snow sublimation
-      !------------------
-!     !++++++
-!     zqt_dummy(:) = 0.0
-!     DO jk = 1, nlay_s
-!        DO ji = kideb,kiut
-!           q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
-!           ! heat conservation
-!           zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * old_ht_s_b(ji) / nlay_s
-!        END DO
-!     END DO
-
-!     IF (jiindex_1d .GT. 0) &
-!     WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &
-!                     rdt_ice
-      !++++++
-!     !++++++
-!     zqt_dummy(:) = 0.0
-!     DO jk = 1, nlay_s
-!        DO ji = kideb,kiut
-!           q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
-!           ! heat conservation
-!           zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
-!        END DO
-!     END DO
-!     IF (jiindex_1d .GT. 0) &
-!     WRITE(numout,*) 'jl : ', jl, ' New zqts : ', zqt_dummy(jiindex_1d) / &
-!                     rdt_ice
-!     !++++++
+      !----------------------
+      ! 3.3 Snow sublimation
+      !----------------------
 
       DO ji = kideb, kiut
@@ -510 +391 @@
          ht_s_b(ji)      =  MAX( zzero , zdhcf )
          ! we recompute dh_s_tot 
-         ! which may be different in case of total snow melt
          dh_s_tot(ji)    =  ht_s_b(ji) - zhsold(ji)
-!        IF ( ht_s_b(ji) .LE. 0.0 ) THEN
-!           dh_s_tot(ji) =  MAX( 0.0 , dh_s_tot(ji) )
-!        ENDIF 
-
          qt_s_in(ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
       END DO  !ji
 
-!     IF ( jiindex_1d .GT. 0 ) THEN
-!        WRITE(numout,*) ' dh_s_sub : ', zdh_s_sub(jiindex_1d)
-!        WRITE(numout,*) ' dhs_tot  : ', dh_s_tot(jiindex_1d)
-!     ENDIF
-
-      !++++++
       zqt_dummy(:) = 0.0
       DO jk = 1, nlay_s
@@ -533 +403 @@
          END DO
       END DO
-!     IF (jiindex_1d .GT. 0) &
-!     WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &
-!                     rdt_ice
-!     !++++++
-
-!     IF (jiindex_1d .GT. 0) THEN
-!        WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
-!        WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
-!     ENDIF
 
       DO jk = 1, nlay_s !n 
@@ -553 +414 @@
       END DO 
 
-!     IF (jiindex_1d .GT. 0) THEN
-!        WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
-!        WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
-!     ENDIF
-
 !
 !------------------------------------------------------------------------------!
@@ -567 +423 @@
       ! the inner conductive flux  (fc_bo_i), from the open water heat flux 
       ! (qlbbqb) and the turbulent ocean flux (fbif). 
-      ! fc_bo_i is positive downwards
-      ! fbif and qlbbq are positive to the ice 
-
-      !---------------------------------------------
-      ! Basal growth - salinity not varying in time 
-      !---------------------------------------------
-      IF ( num_sal .NE. 2 ) THEN
+      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 
+
+      !-----------------------------------------------------
+      ! 4.1 Basal growth - (a) salinity not varying in time 
+      !-----------------------------------------------------
+      IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN
          DO ji = kideb, kiut
             IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
@@ -590 +445 @@
                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
                                       qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
-!              IF ( ji .EQ. jiindex_1d ) THEN
-!                 WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
-!                 WRITE(numout,*) ' ztmelts : ', ztmelts
-!                 WRITE(numout,*) ' t_bo    : ', t_bo_b(ji)
-!                 WRITE(numout,*) ' q_i_b   : ', q_i_b(ji,nlay_i+1)
-!                 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!              ENDIF
             ENDIF ! heat budget
          END DO ! ji
       ENDIF ! num_sal
 
-!     IF ( jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
-!     WRITE(numout,*) ' q_i_b     : ', q_i_b(jiindex_1d,nlay_i+1)
-!     WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(jiindex_1d)
-!     WRITE(numout,*) ' fbif_1d   : ', fbif_1d(jiindex_1d)
-!     WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(jiindex_1d)
-!     ENDIF
-
-      !-----------------------------------------
-      ! Basal growth - salinity varying in time 
-      !-----------------------------------------
+      !-------------------------------------------------
+      ! 4.1 Basal growth - (b) salinity varying in time 
+      !-------------------------------------------------
       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
          ! the growth rate (dh_i_bott) is function of the new ice
@@ -617 +457 @@
          ! salinity (snewice). snewice depends on dh_i_bott
          ! it converges quickly, so, no problem
+         ! See Vancoppenolle et al., OM08 for more info on this
 
          ! Initial value (tested 1D, can be anything between 1 and 20)
@@ -636 +477 @@
                                                  ( t_bo_b(ji) - rtt ) )       &
                                       - rcp * ( ztmelts-rtt ) )
-!              !+++++++
-!              IF (ji.EQ.jiindex_1d) THEN
-!                 WRITE(numout,*) ' ztmelts            : ', ztmelts
-!                 WRITE(numout,*) ' t_bo_b             : ', t_bo_b(ji)
-!                 WRITE(numout,*) ' q_i_b(ji,nlay_i+1) : ', q_i_b(ji,nlay_i+1)
-!              ENDIF
-!              !+++++++
                   ! Bottom growth rate = - F*dt / q
                   dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
@@ -650 +484 @@
                   ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
                   ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
-!                 zgrr   = MAX( dh_i_bott(ji) / rdt_ice, 0.0 )
                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) )
                   zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
                   zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
                   zswi1  = 1. - zswi2 * zswi12 
-!                 IF (ji.EQ.jiindex_1d) THEN
-!                    WRITE(numout,*) ' zgrr  : ', zgrr
-!                    WRITE(numout,*) ' zswi2 : ', zswi2, ' zswi12 :', zswi12, ' zswi1 :', zswi1
-!                 ENDIF
                   zfracs = zswi1  * 0.12 +  &
                            zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
@@ -673 +502 @@
             IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
                ! New ice salinity must not exceed 15 psu
-               zoldsinew   = s_i_new(ji)
                s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
                ! Metling point in K
@@ -684 +512 @@
                                    - rcp * ( ztmelts-rtt ) )
                ! Basal growth rate = - F*dt / q
-               ! sometimes, growth is very high because of very high bottom cond
-               ! flux... this is dangerous
                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
                                       qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
-      ! new lines
-      !-----------------
-      ! Salinity update
-      !-----------------
+               ! Salinity update
                ! entrapment during bottom growth
                dsm_i_se_1d(ji) =  ( s_i_new(ji)*dh_i_bott(ji) +              & 
@@ -701 +524 @@
       ENDIF ! num_sal
 
-!     IF ( jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
-!     WRITE(numout,*) ' s_i_new   : ', s_i_new(jiindex_1d)
-!     ENDIF
-
-      !------------
-      ! Basal melt
-      !------------
-      !++++++
+      !----------------
+      ! 4.2 Basal melt
+      !----------------
       meance_dh = 0.0
       numce_dh = 0
       innermelt(:) = 0
-      !++++++
 
       DO ji = kideb, kiut
@@ -729 +545 @@
       END DO
 
-!     IF ( jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
-!     WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(jiindex_1d)
-!     WRITE(numout,*) ' fbif      : ', fbif_1d(jiindex_1d)
-!     WRITE(numout,*) ' qlbbq     : ', qlbbq_1d(jiindex_1d)
-!     WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
-!     ENDIF
-
       DO jk = nlay_i, 1, -1
          DO ji = kideb, kiut
             IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
                ztmelts             =   - tmut * s_i_b(ji,jk) + rtt 
-               ! sometimes internal temperature gt melting point 
-               ! the whole layer is removed
-!              IF (ji.EQ.jiindex_1d) THEN
-!                 WRITE(numout,*) ' jk : ', jk
-!                 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!                 WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)
-!                 WRITE(numout,*) ' t_i_b     : ', t_i_b(ji,jk)
-!                 WRITE(numout,*) ' q_i_b     : ', q_i_b(ji,jk)
-!              ENDIF
                IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN
                   zdeltah(ji,jk)  = - zh_i(ji)
-      !           zqfont_bo(ji)   = zqfont_bo(ji)
                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
-                  !++++++
                   innermelt(ji)   = 1
-                  !++++++
-!                 IF (ji.EQ.jiindex_1d) THEN
-!                    WRITE(numout,*) ' Inner melt: ', innermelt(ji)
-!                    WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!                    WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!                 ENDIF
                ELSE  ! normal ablation
                   zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
@@ -774 +565 @@
                   zjj             = ( npb(ji) - 1 ) / jpi + 1
                   zfsalt_melt(ji) = zfsalt_melt(ji) +                         &
-!                                  ( sss_io(zji,zjj) - s_i_b(ji,jk) ) *       &
                                    ( sss_io(zji,zjj) - sm_i_b(ji)   ) *       &
                                    a_i_b(ji) * &
                                    MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
-!                 IF (ji.EQ.jiindex_1d) THEN
-!                    WRITE(numout,*) ' No inner melt '
-!                    WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!                    WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)
-!                    WRITE(numout,*) ' t_i_b     : ', t_i_b(ji,jk)
-!                    WRITE(numout,*) ' q_i_b     : ', q_i_b(ji,jk)
-!                    WRITE(numout,*)
-!                 ENDIF
                ENDIF
-!           IF ( ji .EQ. jiindex_1d ) THEN
-!              WRITE(numout,*) ' basal melt '
-!              WRITE(numout,*) ' sss_io      : ', sss_io(zji,zjj)
-!              WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji)
-!              WRITE(numout,*) ' zfsalt_melt good units : ', zfsalt_melt(ji) / ( sss_io(zji,zjj) + epsi16 )
-!           ENDIF
-            !+++++
             ENDIF
          END DO ! ji
       END DO ! jk
 
+      !-------------------
+      ! Conservation test
+      !-------------------
       IF ( con_i ) THEN
       DO ji = kideb, kiut
          IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
-            !-------------------
-            ! Conservation test
-            !-------------------
             IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
                numce_dh = numce_dh + 1
                meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
             ENDIF
-            IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!.GE. 1.0e-3 ) .AND. &
-!                ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN
+            IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
                 WRITE(numout,*) ' ALERTE heat loss for basal  melt '
                 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
@@ -834 +608 @@
       ENDIF ! con_i
 
-!     IF ( jiindex_1d .GT. 0 ) THEN
-!        WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
-!        WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
-!     ENDIF
 !
 !------------------------------------------------------------------------------!
@@ -843 +613 @@
 !------------------------------------------------------------------------------!
 !
-      !------------------------------------------
-      ! Excessive ablation in a 1-category model
-      !------------------------------------------
+      !----------------------------------------------
+      ! 5.1 Excessive ablation in a 1-category model
+      !----------------------------------------------
 
       DO ji = kideb, kiut
@@ -872 +642 @@
       END DO
 
-!     IF (jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     ENDIF
-
-      !--------------------------
-      ! If too much ice melts
-      !--------------------------
+      !-----------------------------------
+      ! 5.2 More than available ice melts
+      !-----------------------------------
       ! then heat applied minus heat content at previous time step
       ! should equal heat remaining 
@@ -890 +655 @@
          zhgnew(ji) =  zihgnew * zhgnew(ji) ! ice thickness is put to 0
          ! remaining heat
-         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) ! &
-!                       + zihgnew * MAX ( zfdt_init(ji) - zqt_i(ji) - zqt_s(ji), 0.0 )
-
-!        !++++++++++
-!        IF (ji.EQ.jiindex_1d) THEN
-!           WRITE(numout,*) ' zfdt_init : ', zfdt_init(ji) / rdt_ice
-!           WRITE(numout,*) ' zqt_i     : ', zqt_i(ji) / rdt_ice
-!           WRITE(numout,*) ' 1. zqt_s  : ', zqt_s(ji) / rdt_ice
-!           WRITE(numout,*) ' zqt       : ', ( zqt_i(ji) + zqt_s(ji) ) / rdt_ice
-!           WRITE(numout,*) ' zhgnew    : ', zhgnew(ji)
-!           WRITE(numout,*) ' zihgnew   : ', zihgnew
-!           WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
-!           WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-!           WRITE(numout,*) ' ht_s_b    : ', ht_s_b(ji)
-!           WRITE(numout,*) ' 1. zfdt_final: ', zfdt_final(ji) / rdt_ice
-!           WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) / rdt_ice
-!           WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
-!        ENDIF
-!        !++++++++++
+         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 
 
          ! If snow remains, energy is used to melt snow
@@ -923 +670 @@
          zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
 
-!        IF (ji.EQ.jiindex_1d) THEN
-!           WRITE(numout,*) ' 2. zfdt_final   : ', zfdt_final(ji) / rdt_ice
-!           WRITE(numout,*) ' ht_s_b       : ', ht_s_b(ji)
-!           WRITE(numout,*) ' zhgnew       : ', zhgnew(ji)
-!           WRITE(numout,*) ' 2. zqt_s     : ', zqt_s(ji) / rdt_ice
-!           WRITE(numout,*) ' zdhnm        : ', zdhnm
-!           WRITE(numout,*) 
-!        ENDIF
-
          ! Mass variations of ice and snow
          !---------------------------------
@@ -942 +680 @@
          ! Remaining heat to the ocean 
          !---------------------------------
-         ! check the switches here
          ! focea is in W.m-2 * dt
-
          focea(ji)  = - zfdt_final(ji) / rdt_ice
 
-!        IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN
-!                WRITE(numout,*) ' focea : ', focea(ji)
-!        ENDIF
-
-      END DO
-
-!     IF (jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     ENDIF
+      END DO
 
       ftotal_fin (:) = zfdt_final(:)  / rdt_ice
-!     IF (jiindex_1d .GT. 1 ) WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)
 
       !---------------------------
@@ -986 +712 @@
                          focea(ji) * a_i_b(ji) * rdt_ice
 
-!        ! astarojna
          zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
          ! equals 0 if ht_i = 0, 1 if ht_i gt 0
-         ! fstbif_1d est cumulatif merde
          fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
-         ! here there is a bug
          qldif_1d(ji)  = qldif_1d(ji)                                         &
                        + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 
@@ -1011 +734 @@
             t_i_b(ji,jk)  =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
             q_i_b(ji,jk)  =  zihgnew * q_i_b(ji,jk)
-!           IF (ji.eq.jiindex_1d) THEN
-!                   WRITE(numout,*) ' jk    : ', jk
-!                   WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)
-!                   WRITE(numout,*) ' t_i_b  : ', t_i_b(jiindex_1d, 1:nlay_i)
-!                   WRITE(numout,*) ' zihgnew : ', zihgnew, ' zhgnew : ', &
-!                   zhgnew(ji)
-!           ENDIF
          END DO
       END DO  ! ji
 
-!     IF (jiindex_1d.GT.0) THEN
-!     WRITE(numout,*) ' --- Vertical profile : '
-!     WRITE(numout,*) ' q_i_b  : ', q_i_b(jiindex_1d, 1:nlay_i)
-!     ENDIF
-
       DO ji = kideb, kiut
          ht_i_b(ji) = zhgnew(ji)
       END DO  ! ji
-
-!     IF (jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     ENDIF
 !
 !------------------------------------------------------------------------------|
@@ -1065 +771 @@
          zji                 = MOD( npb(ji) - 1, jpi ) + 1
          zjj                 = ( npb(ji) - 1 ) / jpi + 1
-!        zsm_snowice  = MIN ( s_i_max, ( rhoic - rhosn ) / rhoic *            &
-!                       sss_io(zji,zjj) )
 
          zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
@@ -1115 +819 @@
       END DO !ji
 
-!     IF (jiindex_1d .GT. 1 ) THEN
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     ENDIF
-
     END SUBROUTINE lim_thd_dh
 #else

trunk/NEMO/LIM_SRC_3/limthd_dif.F90

@@ -1 +1 @@
 MODULE limthd_dif
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_dif ***
@@ -14 +17 @@
    USE thd_ice
    USE iceini
-   USE limicepoints
    USE limistate
    USE in_out_manager
@@ -34 +36 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2005)
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !!   (c) UCL-ASTR and Martin Vancoppenolle
    !!----------------------------------------------------------------------
@@ -88 +90 @@
         !! ** History :
         !!           (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
-        !!           (06-2005) Martin Vancoppenolle still!!! for the 3d version...
+        !!           (06-2005) Martin Vancoppenolle, 3d version
         !!           (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
         !!           (04-2007) Energy conservation tested by M. Vancoppenolle
@@ -106 +108 @@
                     layer,    &   ! vertical dummy loop index 
                     nconv,    &   ! number of iterations in iterative procedure
-!                   nconvmax = 50, &! maximum number of iterations for temperature computations (50)
-                    index,    &   !
                     minnumeqmin, & !
                     maxnumeqmax
@@ -158 +158 @@
           zfsw     ,   & !solar radiation absorbed at the surface
           zf       ,   & ! surface flux function
-          dzf      ,   & ! derivative of the surface flux function
-          zksn           !effective snow conductivity
+          dzf            ! derivative of the surface flux function
 
        REAL(wp)  ::           &  ! constant values
@@ -168 +167 @@
           zbeta     =  0.117,     & !: for thermal conductivity (could be 0.13)
           zraext_s  =  1.0e08,    & !: extinction coefficient of radiation in the snow
-!         zraext_i  =  1.5,       & !extinction coefficient of radiation in the ice
-!         zraext_i  =  1.0,       & !extinction coefficient of radiation in the ice
-                                    ! MY 0.87, FY 1.2 (Grenfell, Perovich and Light
-                                    ! 07)
-!         zerrmax   =  1.0e-4 ,   & !: error criterion for convergence of the iterative process
-          zfc_int           ,     & !: conductive heat flux at the interface
           zkimin    =  0.10 ,     & !: minimum ice thermal conductivity
           zht_smin  =  1.0e-4       !: minimum snow depth
 
         REAL(wp)  ::          &  ! local variables 
-          zheshth,            &  ! = zhe(ji) / thth
           ztmelt_i,           &  ! ice melting temperature
-          zerritmax,          &  ! current maximal error on temperature 
-          zexp
+          zerritmax              ! current maximal error on temperature 
 
         REAL(wp), DIMENSION(jpij)  :: &
           zerrit,             &  ! current error on temperature 
           zdifcase,           &  ! case of the equation resolution (1->4)
-          zghe,               &  ! correction factor of the thermal conductivity
-          zhe,                &  ! effective thickness for compu. of equ. thermal conductivity
           zftrice,            &  ! solar radiation transmitted through the ice
-          zihic, zihe, zhsu
-
-        CHARACTER (len=50) ::   charout
+          zihic, zhsu
+
 !!-- End of declarations
 !!----------------------------------------------------------------------------------------------
 
-       ! new namelist parameters
-!      maxer_i_thd  = maxer_i_thd
-!      thcon_i_swi = thcon_i_swi
-!      kappa_i = kappa_i
-!      nconv_i_thd = nconv_i_thd
-
        IF(lwp) WRITE(numout,*)'lim_thd_dif : Heat diffusion in sea ice for cat :', jl
 
-       thcon_i_swi = 1 ! Pringle et al. (2007)
 !
 !------------------------------------------------------------------------------!
@@ -339 +320 @@
                         zradtr_i(ji,nlay_i)
        END DO
+       ! +++++
+
        DO layer = 1, nlay_i
           DO ji = kideb , kiut
@@ -345 +328 @@
        END DO
 
-       ! +++++
        
 !
@@ -379 +361 @@
        DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd))
 
-!      ! +++++ 
-!         DO ji = kideb , kiut
-!         zji                 = MOD( npb(ji) - 1, jpi ) + 1
-!         zjj                 = ( npb(ji) - 1 ) / jpi + 1
-!         IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN
-!            WRITE(numout,*) ' Iteration : ', nconv
-!            WRITE(numout,*) ' t_su     : ', t_su_b(ji)
-!            WRITE(numout,*) ' t_s      : ', t_s_b(ji,1)
-!            WRITE(numout,*) ' t_i      : ', t_i_b(ji,:)
-!            WRITE(numout,*) ' difcase  : ', zdifcase(ji)
-!         ENDIF
-!         END DO
-!      ! +++++ 
-
        nconv   =  nconv+1
 
@@ -410 +378 @@
 
        IF ( thcon_i_swi .EQ. 1 ) THEN
-       ! Pringle et al formula included, should be tested
+       ! Pringle et al formula included,
        ! 2.11 + 0.09 S/T - 0.011.T
        DO ji = kideb , kiut
@@ -489 +457 @@
           zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(zeps,zh_i(ji))
           zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji))
-
        !-- Interface
           zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, &
@@ -616 +583 @@
                    numeqmin(ji)    =  1
                    numeqmax(ji)    =  nlay_i + nlay_s + 1
+
                    !!surface equation
                    ztrid(ji,1,1) = 0.0
                    ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
                    ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
-                   ! bug again ???
-!                  zindterm(ji,1) = dzf(ji)*ztsuold(ji)   - zf(ji)
                    zindterm(ji,1) = dzf(ji)*t_su_b(ji)   - zf(ji)
 
@@ -853 +819 @@
        END DO
 
+       !-------------------------!
+       ! Heat conservation       !
+       !-------------------------!
        IF ( con_i ) THEN
-       !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-       ! ++++ heat conservation ++++ '
-       ! Internal fluxes ... goes to trash!!!
 
        DO ji = kideb, kiut
@@ -884 +850 @@
           zji                 = MOD( npb(ji) - 1, jpi ) + 1
           zjj                 = ( npb(ji) - 1 ) / jpi + 1
-!         IF ( ( zji.EQ.jiindex ) .AND. ( zjj.EQ.jjindex) ) THEN
-!            WRITE(numout,*) ' layer    : ', layer
-!            WRITE(numout,*) ' fc_i     : ', fc_i(ji,layer)
-!            WRITE(numout,*) ' zkappa_i : ', zkappa_i(ji,layer)
-!            WRITE(numout,*) ' ztcond_i : ', ztcond_i(ji,layer)
-!            WRITE(numout,*) ' zspeche_i: ', zspeche_i(ji,layer)
-!            WRITE(numout,*) ' ti2      : ', t_i_b(ji,layer+1)
-!            WRITE(numout,*) ' ti1      : ', t_i_b(ji,layer)
-!         ENDIF
        END DO
        END DO
@@ -903 +860 @@
 
        ENDIF
-
-       !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
     END SUBROUTINE lim_thd_dif

trunk/NEMO/LIM_SRC_3/limthd_ent.F90

@@ -1 +1 @@
 MODULE limthd_ent
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_ent   ***
@@ -7 +10 @@
    !!                       after vertical growth/decay
    !!======================================================================
-
-   !!----------------------------------------------------------------------
-   !!   lim_thd_ent : ice temperature redistribution 
+   !!   lim_thd_ent : ice redistribution of enthalpy
 
    !! * Modules used
@@ -22 +23 @@
    USE limistate
    USE ice
-!  USE limthd
    USE limvar
    USE par_ice
-   USE limicepoints
 
    IMPLICIT NONE
@@ -42 +41 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2005)
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -59 +58 @@
       !! ** Method  : linear conservative remapping
       !!           
-      !! ** Steps : 1) Switches
-      !!            2) Snow redistribution
-      !!            3) Ice enthalpy redistribution
-      !!            4) Ice salinity 
-      !!            5) Inverse temperature computation from enthalpy
+      !! ** Steps : 1) Grid
+      !!            2) Switches
+      !!            3) Snow redistribution
+      !!            4) Ice enthalpy redistribution
+      !!            5) Ice salinity, recover temperature
       !!
       !! ** Arguments
@@ -76 +75 @@
       !!               (07-2005) Martin for 3d adapatation
       !!               (11-2006) Vectorized by Xavier Fettweis (ASTR)
+      !!               (03-2008) Energy conservation and clean code
       !! * Arguments
 
@@ -84 +84 @@
 
       INTEGER ::            &
-         ji,jk,jk_a     ,   &  !  dummy loop indices
-         zji, zjj ,         &  ! dummy indices
+         ji,jk          ,   &  !  dummy loop indices
+         zji, zjj       ,   &  !  dummy indices
          ntop0          ,   &  !  old layer top index
          nbot1          ,   &  !  new layer bottom index
@@ -140 +140 @@
       ! Energy conservation
       REAL(wp), DIMENSION(jpij) :: &
-         zqti_in, zqts_in, zqt_in, &
-         zqti_fin, zqts_fin, zqt_fin
+         zqti_in, zqts_in,         &
+         zqti_fin, zqts_fin
+
+!------------------------------------------------------------------------------|
 
       zeps   = 1.0d-20
       zeps6  = 1.0d-06
-
-      IF(lwp) THEN
-         WRITE(numout,*) 
-         WRITE(numout,*) 'lim_thd_ent : Redistribution of energy in the ice ' 
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) ' kideb : ', kideb
-         WRITE(numout,*) ' kiut  : ', kiut 
-      ENDIF
-!
-!------------------------------------------------------------------------------|
-!  1) Number of layers / Thickness of the layers                               |
-!------------------------------------------------------------------------------|
-!
-      ! Number of layers depends on the ice thickness
-      ! 1 layer if ice thinner than 30 cm
-      ! 2 layers if ice thinner than 60 cm
-      ! 4 layers otherwise
-
-      nlays0        = nlay_s
-      nlayi0        = nlay_i
-
-      DO ji = kideb, kiut
-         zh_i(ji)   = old_ht_i_b(ji) / nlay_i 
-         zh_s(ji)   = old_ht_s_b(ji) / nlay_s
-      ENDDO
-
       zthick0(:,:) = 0.0
       zm0(:,:)     = 0.0
@@ -178 +154 @@
       z_i(:,:)     = 0.0
       z_s(:,:)     = 0.0
+
+!
+!------------------------------------------------------------------------------|
+!  1) Grid                                                                     |
+!------------------------------------------------------------------------------|
+!
+      nlays0        = nlay_s
+      nlayi0        = nlay_i
+
+      DO ji = kideb, kiut
+         zh_i(ji)   = old_ht_i_b(ji) / nlay_i 
+         zh_s(ji)   = old_ht_s_b(ji) / nlay_s
+      ENDDO
+
 !
 !------------------------------------------------------------------------------|
@@ -294 +284 @@
       ENDDO
 
-!     IF (jiindex_1d .GT. 0 ) THEN
-!     WRITE(numout,*) ' icsuind : ', icsuind(jiindex_1d)
-!     WRITE(numout,*) ' icsuswi : ', icsuswi(jiindex_1d)
-!     WRITE(numout,*) ' icboind : ', icboind(jiindex_1d)
-!     WRITE(numout,*) ' icboswi : ', icboswi(jiindex_1d)
-!     WRITE(numout,*) ' snicind : ', snicind(jiindex_1d)
-!     WRITE(numout,*) ' snicswi : ', snicswi(jiindex_1d)
-!     ENDIF
 !
 !------------------------------------------------------------------------------|
@@ -361 +343 @@
       ENDDO
 
-!     IF (jiindex_1d .GT. 0 ) THEN
-!        WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice
-!        WRITE(numout,*) ' q_s_b  : ', q_s_b(jiindex_1d,1)
-!        WRITE(numout,*) ' t_s_b  : ', t_s_b(jiindex_1d,1)
-!        WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1)
-!     ENDIF
-
-!     WRITE(numout,*) ' TEST : '
-!     WRITE(numout,*) ' maxnbot0 : ', maxnbot0
-!     WRITE(numout,*) ' kideb : ', kideb
-!     WRITE(numout,*) ' kiut  : ', kiut 
       DO jk = 2, maxnbot0
         DO ji = kideb, kiut
@@ -381 +352 @@
           zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, zeps - ht_s_b(ji) ) )
           zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch
-!         IF (ji.EQ.jiindex_1d) THEN
-!                 WRITE(numout,*) ' limsum ', limsum
-!                 WRITE(numout,*) ' snswi  : ', snswi(ji)
-!                 WRITE(numout,*) ' snind  : ', snind(ji)
-!                 WRITE(numout,*) ' t_s_b  : ', t_s_b(ji,1)
-!                 WRITE(numout,*) ' q_s_b  : ', rhosn * ( cpic * ( rtt - &
-!                                               t_s_b(ji,limsum) ) + lfus )
-!                 WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,jk)
-!         ENDIF
         END DO ! jk
       ENDDO ! ji
-
-      ! Heat conservation
-!     IF (jiindex_1d .GT. 0 ) THEN
-!        WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice
-!        WRITE(numout,*) ' q_s_b  : ', q_s_b(jiindex_1d,1)
-!        WRITE(numout,*) ' t_s_b  : ', t_s_b(jiindex_1d,1)
-!     ENDIF
 
       !------------------------------------------------
@@ -423 +378 @@
 
       !------------------
-      ! NEW SNOW PROFILE
+      ! new snow profile
       !------------------
 
@@ -460 +415 @@
         END DO
       ENDDO
-
-!     IF (jiindex_1d .GT. 0) THEN
-!        WRITE(numout,*) ' qm0      : ', qm0(jiindex_1d,0:maxnbot0) / rdt_ice
-!        WRITE(numout,*) ' maxnbot0 : ', maxnbot0
-!        WRITE(numout,*) ' zm0      : ', zm0(jiindex_1d,0:maxnbot0)
-!        WRITE(numout,*) ' z_s      : ', z_s(jiindex_1d,0:1)
-!        WRITE(numout,*) ' zhl0     : ', zhl0(jiindex_1d,0:maxnbot0)
-!        WRITE(numout,*) ' zhl0     : ', zhl0(jiindex_1d,0:maxnbot0)
-!     ENDIF
 
       !----------------
@@ -480 +426 @@
               q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) &
                                    * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps))
-!              IF (ji.EQ.jiindex_1d) THEN
-!                      WRITE(numout,*) ' layer0, layer1 ', layer0, layer1
-!                      WRITE(numout,*) ' q_s_b : ', q_s_b(ji,layer1) / rdt_ice
-!                      WRITE(numout,*) ' nbot0 : ', nbot0(ji)                 
-!                      WRITE(numout,*)
-!              ENDIF
            END DO
         END DO
@@ -516 +456 @@
       ENDIF
 
-!     WRITE(numout,*) ' ht_s_b   : ', ht_s_b(jiindex_1d)
-!     WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(jiindex_1d)
-!     WRITE(numout,*) ' dh_snowi : ', dh_snowice(jiindex_1d)
-!     WRITE(numout,*) ' zqts_fin : ', zqts_fin(jiindex_1d) / rdt_ice
-!     WRITE(numout,*) ' difference : ', ( zqts_fin(jiindex_1d) - &
-!                                       zqts_in(jiindex_1d) ) / rdt_ice
-
-!     IF (jiindex_1d .GT. 1) THEN
-!     WRITE(numout,*) ' Vertical profile : '
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     WRITE(numout,*) ' nbot0  : ', nbot0(jiindex_1d)
-!     WRITE(numout,*) ' qm0    : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &
-!     rdt_ice
-!     WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1:nbot0(jiindex_1d))
-!     WRITE(numout,*) ' zm0    : ', zm0(jiindex_1d,0:nbot0(jiindex_1d))
-!     ENDIF
-
       !---------------------
       ! Recover heat content
@@ -541 +464 @@
          END DO !ji
       ENDDO !jk  
-!     IF (jiindex_1d .GT. 1) THEN
-!     WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
-!     WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
-!     ENDIF
 
       !---------------------
       ! Recover temperature
       !---------------------
-      ! is it really necessary ???? since we don't use snow temperature anymore
-      ! it is now since, lim_thd_glohec uses t_s_b and not q_s_b!!!
-      ! it should not be necessary anymore
-
       zfac1 = 1. / ( rhosn * cpic )
       zfac2 = lfus / cpic  
@@ -563 +478 @@
         END DO
       ENDDO
-!     IF (jiindex_1d .GT. 1) THEN
-!     WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
-!     WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
-!     ENDIF
 !
 !------------------------------------------------------------------------------|
@@ -587 +498 @@
                          ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) ,    &
                          nlay_i + 2 ) )
-!       IF (ji.EQ.jiindex_1D) THEN
-!           WRITE(numout,*) ' nbot0 : ', nbot0(ji)
-!           WRITE(numout,*) ' nlayi0 : ', nlayi0
-!           WRITE(numout,*) ' icboind: ', icboind(ji)
-!           WRITE(numout,*) ' icsuind: ', icsuind(ji)
-!           WRITE(numout,*) ' icsuswi: ', icsuswi(ji)
-!           WRITE(numout,*) ' snicswi: ', snicswi(ji)
-!           WRITE(numout,*) ' jkmax  : ', jkmax
-!           WRITE(numout,*) ' term   : ', &
-!           nlayi0+(1-icboind(ji))+(1-icsuind(ji))*icsuswi(ji)+snicswi(ji)
-!       ENDIF
          ! maximum reference number of the bottommost layer over all domain
          maxnbot0     =  MAX( maxnbot0 , nbot0(ji) )
@@ -612 +512 @@
             ! the ice layer number goes from 1 to nlay_i
             ! limsum is the real ice layer number corresponding to present jk
-            limsum    =  ( (icsuswi(ji)*(icsuind(ji)+jk-1) + (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji)
+            limsum    =  ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 
+                           (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji)
             zm0(ji,jk)=  icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) &
                       +  limsum * zh_i(ji)
@@ -647 +548 @@
                       MIN((t_i_b(ji,limsum)-rtt),-zeps) ) - rcp*(ztmelts-rtt) ) &
                       * zthick0(ji,jk)
-
-!           IF (ji .EQ. jiindex_1d ) THEN
-!              WRITE(numout,*) 'jk : ', jk, ' t_i_b : ', t_i_b(ji,limsum), &
-!              ' s_i_b : ', s_i_b(ji,limsum)
-!              WRITE(numout,*) ' q_i_b  : ', &
-!              rhoic*(cpic*(ztmelts-t_i_b(ji,limsum)) + lfus *( &
-!              1.0-(ztmelts-rtt)/ &
-!                                   MIN((t_i_b(ji,limsum)-rtt),-zeps) ) - &
-!                                   rcp*(ztmelts-rtt) )
-!              WRITE(numout,*) ' zthick0: ', zthick0(ji,jk)
-!              WRITE(numout,*) ' limsum : ', limsum
-!              WRITE(numout,*) ' icsuswi: ', icsuswi(ji), ' icsuind : ', &
-!              icsuind(ji)
-!              WRITE(numout,*) ' snicswi: ', snicswi(ji)
-!              WRITE(numout,*) ' FAC1   : ', snicswi(ji)*(jk-1)
-!              WRITE(numout,*) ' FAC2   : ', icsuswi(ji)*(jk-1+icsuind(ji))
-!              WRITE(numout,*) ' FAC3   : ', (1-icsuswi(ji))*(1-snicswi(ji))*jk
-!           ENDIF
          END DO
       ENDDO
-!     IF (jiindex_1d .GT.0) WRITE(numout,*) ' qm0    : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &
-!     rdt_ice
-
-!     ! Heat conservation
-!     WRITE(numout,*) ' zqts_in : ', zqts_in(jiindex_1d) / rdt_ice
 
       !----------------------------
@@ -696 +574 @@
       ! Snow ice layer heat content
       !-----------------------------
-
       DO ji = kideb, kiut
          ! energy of the flooding seawater
@@ -704 +581 @@
          qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic 
 
-         ! The heat of the flooding seawater is removed from the lead
-! default line
-!        qldif_1d(ji)   = qldif_1d(ji) - zqsnic * a_i_b(ji)
-! LIM 2 line
-!         qldif_1d(ji)  = qldif_1d(ji) + zdrmh * ( 1.0 - alphs * ( rhosn/rhoic ) ) * xlic * ( 1.0 - frld_1d(ji) )
-! debug line : test a plus (should also be positive for the ice)
          qldif_1d(ji)   = qldif_1d(ji) + zqsnic * a_i_b(ji)
-!        qldif_1d(ji) = qldif_1d(ji) + rhosn * cpic *( t_bo_b(ji)-t_s_b(ji,1) ) * & 
-!                                      a_i_b(ji)*dh_snowice(ji)
-         !++++++
-!        IF ( ji.EQ.jiindex_1D ) THEN
-!           WRITE(numout,*) ' zqsnic : ', zqsnic / rdt_ice
-!           WRITE(numout,*) ' zqsnow : ', zqsnow(ji) / rdt_ice
-!        ENDIF
-!        !++++++
 
          ! enthalpy of the newly formed snow-ice layer
@@ -725 +588 @@
          qm0(ji,1)      =  snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1)
 
-!        !++++++
-!        IF (ji.EQ.jiindex_1D) THEN
-!           WRITE(numout,*) ' snicswi : ', snicswi(ji) 
-!           WRITE(numout,*) ' zqsnic  : ', zqsnic / rdt_ice
-!           WRITE(numout,*) ' qldif_1d : ', qldif_1d(ji)
-!           WRITE(numout,*) ' dhsnowice :', dh_snowice(ji)
-!           WRITE(numout,*) ' t_bo_b    :', t_bo_b(ji)   
-!           WRITE(numout,*) ' t_s_b     :', t_s_b(ji,1)   
-!           WRITE(numout,*) ' limsum    :', limsum   
-!           WRITE(numout,*) ' a_i_b     :', a_i_b(ji)
-!        ENDIF
-!        !++++++
-
       ENDDO ! ji
-
-!     IF (jiindex_1D.GT.1) THEN
-!     WRITE(numout,*) ' Vertical profile : '
-!     WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' nbot0  : ', nbot0(jiindex_1d)
-!     WRITE(numout,*) ' qm0    : ', qm0(jiindex_1d,1:nbot0(jiindex_1d)) / &
-!     rdt_ice
-!     WRITE(numout,*) ' zthick0: ', zthick0(jiindex_1d,1:nbot0(jiindex_1d))
-!     WRITE(numout,*) ' zm0    : ', zm0(jiindex_1d,0:nbot0(jiindex_1d))
-!     ENDIF
 
       DO jk = ntop0, maxnbot0
@@ -757 +596 @@
                        * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-zeps6+zeps) ) &
                        * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + zeps ) )
-!          IF (ji.EQ.jiindex_1D) THEN
-!                  WRITE(numout,*) ' jk  : ', jk
-!                  WRITE(numout,*) ' qm0 : ', qm0(ji,jk)
-!                  WRITE(numout,*) ' nbot0:', nbot0(ji)
-!                  WRITE(numout,*) ' fac 1 : ', MAX( 0.0 , SIGN( 1.0, nbot0(ji) &
-!                                               - jk + zeps ) )
-!          ENDIF
-        END DO
-      ENDDO
-!     IF (jiindex_1D .GT. 0 ) THEN
-!        WRITE(numout,*) ' qm0     : ', qm0(jiindex_1d,1:jkmax) / rdt_ice
-!        WRITE(numout,*) ' zqti_in : ', zqti_in(jiindex_1d) / rdt_ice
-!     ENDIF
-      ! Heat conservation
-!     WRITE(numout,*) ' zm0     : ', zm0(jiindex_1d,:) 
-!     WRITE(numout,*) ' zthick0 : ', zthick0(jiindex_1d,:) 
-!     WRITE(numout,*) ' qm0     : ', qm0(jiindex_1d,0:jkmax) / rdt_ice
+        END DO
+      ENDDO
 
       !-------------
@@ -803 +627 @@
       ENDDO
 
-      ! What is the difference with zthick0 ???
       !--thicknesses of the layers
       DO layer0 = ntop0, maxnbot0
@@ -810 +633 @@
         END DO
       ENDDO
-
-
-      ! Use less ??
-!     DO layer1 = ntop1, nbot1
-!       DO ji = kideb, kiut
-!           IF (ji.eq.jiindex_1d) THEN
-!                   WRITE(numout,*) ' jk    : ', layer1
-!                   WRITE(numout,*) ' q_i_b : ', q_i_b(ji,layer1)
-!           ENDIF
-!           q_i_b(ji,layer1)  =  q_i_b(ji,layer1) * (1.0 - MAX(0.0,SIGN(1.0,ht_i_b(ji)-zeps6+zeps)))
-!           IF (ji.eq.jiindex_1d) THEN
-!                   WRITE(numout,*) ' ht_i_b  ', ht_i_b(ji)
-!                   WRITE(numout,*) ' factor2 ', ht_i_b(ji)-zeps6+zeps
-!                   WRITE(numout,*) ' factor .... ',  (1.0 - &
-!                   MAX(0.0,SIGN(1.0,ht_i_b(ji)-zeps6+zeps)))
-!                   WRITE(numout,*) ' q_i_b : ', q_i_b(ji,layer1)
-!           ENDIF
-!       END DO
-!     ENDDO
 
       !------------------------
@@ -848 +652 @@
       ENDDO
 
-      
-
- !    WRITE(numout,*) ' s_i       : ', s_i_b(jiindex_1d,1:nlay_i)
-      
       !-------------------------
       ! Heat conservation check
@@ -861 +661 @@
          END DO
       END DO
-!     IF ( jiindex_1d .GT. 0 ) &
-!     WRITE(numout,*) ' zqti_fin : ', zqti_fin(jiindex_1d) / rdt_ice
 !
       DO ji = kideb, kiut
@@ -882 +680 @@
          ENDIF
       END DO
-!
-!     WRITE(numout,*) ' old_ht_i   : ', old_ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' ht_i_b     : ', ht_i_b(jiindex_1d)
-!     WRITE(numout,*) ' icsuswi    : ', icsuswi(jiindex_1d)
-!     WRITE(numout,*) ' icboswi    : ', icboswi(jiindex_1d)
-!     WRITE(numout,*) ' snicswi    : ', snicswi(jiindex_1d)
-!     WRITE(numout,*) ' dh_i_bott  : ', dh_i_bott(jiindex_1d)
-!     WRITE(numout,*) ' dh_i_surf  : ', dh_i_surf(jiindex_1d)
-!     WRITE(numout,*) ' dh_snowice : ', dh_snowice(jiindex_1d)
-!     WRITE(numout,*) ' zqti_fin   : ', zqti_fin(jiindex_1d) / rdt_ice
-!     WRITE(numout,*) ' difference : ', ( zqti_fin(jiindex_1d) - &
-!                                       zqti_in(jiindex_1d) ) / rdt_ice
 
       !----------------------
@@ -904 +690 @@
       ENDDO !jk  
 
-      !++++++++++++++++
       ! Heat conservation
       zqti_fin(:) = 0.0
@@ -912 +697 @@
          END DO
       END DO
-!     IF (jiindex_1d .GT. 0 ) &
-!     WRITE(numout,*) ' zqti_fin : ', zqti_fin(jiindex_1d) / rdt_ice
-!     WRITE(numout,*) ' q_i_b    : ', q_i_b(jiindex_1d,1:nlay_i) 
-      !++++++++++++++++
-      !-------------------------------------------------
-      ! Update salinity, bottom and snow ice entrapment 
-      !-------------------------------------------------
+
+!
+!------------------------------------------------------------------------------|
+!  5) Update salinity and recover temperature                                  |
+!------------------------------------------------------------------------------|
+!
+      ! Update salinity (basal entrapment, snow ice formation)
       DO ji = kideb, kiut
          sm_i_b(ji) = sm_i_b(ji)                                &
@@ -924 +709 @@
       END DO !ji
 
-      !---------------------
       ! Recover temperature
-      !---------------------
       DO jk = 1, nlay_i
 
@@ -944 +727 @@
       END DO !jk
 
-! Fin de la routine Enth
       END SUBROUTINE lim_thd_ent
 
@@ -957 +739 @@
 #endif
  END MODULE limthd_ent
-

trunk/NEMO/LIM_SRC_3/limthd_lac.F90

@@ -1 +1 @@
 MODULE limthd_lac
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_lac   ***
@@ -20 +23 @@
    USE iceini
    USE limtab
-   USE limicepoints
    USE taumod
    USE blk_oce
@@ -43 +45 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008) 
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limthd_lac.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
@@ -79 +81 @@
       !!
       !! History :
-      !!   1.0  !  01-04 (T. Fichefet, M. A. Morales Maqueda, H. Goosse) 
-      !!        !        original code
-      !!   2.0  !  02-08 (C. Ethe, G. Madec) F90, mpp
       !!   3.0  !  12-05 (M. Vancoppenolle)  Thorough rewrite of the routine
       !!                                     Salinity variations in sea ice, 
       !!                                     Multi-layer code
-      !!   3.1  !  01-06 (M. Vancoppenolle)  Ice thickness distribution
+      !!   3.1  !  01-06 (M. Vancoppenolle)  ITD
       !!   3.2  !  04-07 (M. Vancoppenolle)  Mass and energy conservation tested
       !!------------------------------------------------------------------------
@@ -185 +184 @@
 
       REAL(wp) ::            &
-         zde          ,      &       ! :increment of energy in category jl
-         zde_old      ,      &       ! :dummy variable for energy conservation
-         zde_new      ,      &       ! :dummy variable for energy conservation
-         zde_nice     ,      &       ! :dummy variable for energy conservation
-         zde_diff                    ! :dummy variable for energy conservation
+         zde                         ! :increment of energy in category jl
 
       CHARACTER (len = 15) :: fieldid
@@ -199 +194 @@
       vt_i_init(:,:) = 0.0
       vt_s_init(:,:) = 0.0
-      IF(lwp) THEN
-         WRITE(numout,*) 'lim_thd_lac : Creating new ice'
-         WRITE(numout,*) '~~~~~~~~~~~'
+      zeps6   = 1.0e-6
+
+!------------------------------------------------------------------------------!
+! 1) Conservation check and changes in each ice category
+!------------------------------------------------------------------------------!
+      IF ( con_i ) THEN
+         CALL lim_column_sum (jpl, v_i, vt_i_init)
+         CALL lim_column_sum (jpl, v_s, vt_s_init)
+         CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
+         CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
       ENDIF
-      zeps6   = 1.0e-6
-
-      !++++++++++
-      WRITE(numout,*) ' Step 0 '
-      WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,1:jpl)
-      WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,1:jpl)
-      DO jk = 1, nlay_i
-        WRITE(numout,*) ' e_i : ', jk, e_i(jiindex,jjindex,jk,1:jpl)
-      END DO
-      !++++++++++
-
-!------------------------------------------------------------------------------!
-! 1) Conservation check and changes in each ice category
-!------------------------------------------------------------------------------!
-      CALL lim_column_sum (jpl, v_i, vt_i_init)
-      CALL lim_column_sum (jpl, v_s, vt_s_init)
-      CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
-      CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
 
 !------------------------------------------------------------------------------|
@@ -497 +481 @@
 ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
 !------------------------------------------------------------------------------!
+
         !----------------------
         ! Thickness of new ice
@@ -511 +496 @@
               zjj                  = ( npac(ji) - 1 ) / jpi + 1
               WRITE(numout,*) ' collection thickness  <= 0 ', zh_newice(ji), ji, zji, zjj
-!             WRITE(numout,*) ' LATITUDE ', gphit(zji,zjj), ' LONGITUDE ', &
-!                                           glamt(zji,zjj)
-!             WRITE(numout,*) ' zh_newice ', zh_newice(ji)
-!             WRITE(numout,*) ' ji ', ji
-!             WRITE(numout,*) ' a_i ', a_i(zji,zjj,1:jpl)
-!             WRITE(numout,*) ' v_i ', v_i(zji,zjj,1:jpl)
            ENDIF
         END DO
@@ -554 +533 @@
                                            / ( t_bo_b(ji) - rtt ) )           &
                                          - rcp * ( ztmelts-rtt ) )
-        !+++++++++++
-        IF ( ji .EQ. jiindex_1D ) THEN
-        WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)
-        ENDIF
-        !+++++++++++
            ze_newice(ji)     =   MAX( ze_newice(ji) , 0.0 ) +                 &
                                  MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) )   & 
                                  * rhoic * lfus
-        !+++++++++++
-        IF ( ji .EQ. jiindex_1D ) THEN
-        WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)
-        ENDIF
-        !+++++++++++
         END DO ! ji
-        !+++++++++++
-        IF ( jiindex_1d .GT. 0) WRITE(numout,*) ' ze_newice : ', ze_newice(jiindex_1d)
-        !+++++++++++
-
         !----------------
         ! Age of new ice
@@ -585 +550 @@
            zqbgow(ji)        = qldif_1d(ji) - qcmif_1d(ji) !<0
         END DO ! ji
-        !+++++++++++
-        IF (jiindex_1d .GT. 0 ) THEN
-        WRITE(numout,*) ' qldif  : ', qldif_1d(jiindex_1d)
-        WRITE(numout,*) ' qcmif  : ', qcmif_1d(jiindex_1d)
-        WRITE(numout,*) ' zqbgow : ', zqbgow(jiindex_1d)
-        ENDIF
-        !+++++++++++
 
         !-------------------
@@ -605 +563 @@
            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
         END DO
-        !+++++++++++
-        IF (jiindex_1d .GT. 0) THEN
-        WRITE(numout,*) ' zv_newice : ', zv_newice(jiindex_1d)
-        ENDIF
-        !+++++++++++
 
         !------------------------------------
@@ -663 +616 @@
               za_newice(ji)  = za_newice(ji) - zda_res(ji)
               zv_newice(ji)  = zv_newice(ji) - zdv_res(ji)
-              IF ( ji .EQ. jiindex_1D ) THEN
-                      WRITE(numout,*) ' zv_newice : ', zv_newice(ji)
-                      WRITE(numout,*) ' zdv_res   : ', zdv_res  (ji)
-              ENDIF
            ELSE
               zda_res(ji) = 0.0
@@ -678 +627 @@
         zat_i_ac(:) = 0.0
 
-        WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)
         DO jl = 1, jpl
            DO ji = 1, nbpac
@@ -691 +639 @@
            END DO ! ji
         END DO ! jl
-        WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)
              
-!       !++++++++++++++++
-        DO ji = 1, nbpac
-!          !+++++
-           IF (zat_i_ac(ji).gt.1.0) THEN
-              zji            = MOD( npac(ji) - 1, jpi ) + 1
-              zjj            = ( npac(ji) - 1 ) / jpi + 1
-              WRITE(numout,*) ' *** ERROR MESSAGE *** '
-              WRITE(numout,*) ' something MUST be wrong '
-              WRITE(numout,*) ' at_i : ', zat_i_ac(ji)
-              WRITE(numout,*) ' is wrong '
-              WRITE(numout,*) ' point : ', zji, zjj
-           ENDIF
-        END DO ! ji
-        !++++++++++++++++
-
         !----------------------------------
         ! Heat content - lateral accretion
@@ -849 +781 @@
               zindb            = 1.0 - MAX( 0.0 , SIGN( 1.0 , -               &
                                  za_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes
-              !zo_i_ac(ji,jl)  = zv_old(ji,jl) * zo_i_ac(ji,jl) / &
-              !MAX(zv_i_ac(ji,jl),zeps) * zindb   
               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) /            &
                                  MAX( za_i_ac(ji,jl) , zeps ) * zindb   
@@ -940 +870 @@
       END DO
      
-      !++++
-      WRITE(numout,*) 'lim_thd_lac  : Salt flux diagnostic '
-      WRITE(numout,*) '~~~~~~~~~~~~'
-      WRITE(numout,*) ' *** Salt fluxes at bottom interface ***'
-      WRITE(numout,*) ' fseqv      : ', fseqv(jiindex,jjindex)
-      WRITE(numout,*)
-      !++++
-
 !------------------------------------------------------------------------------|
 ! 10) Conservation check and changes in each ice category
 !------------------------------------------------------------------------------|
 
+      IF ( con_i ) THEN 
       CALL lim_column_sum (jpl,   v_i, vt_i_final)
       fieldid = 'v_i, limthd_lac'
@@ -972 +895 @@
       WRITE(numout,*) ' et_i_init : ', et_i_init(jiindex,jjindex)
       WRITE(numout,*) ' et_i_final: ', et_i_final(jiindex,jjindex)
+
+      ENDIF
 
    END SUBROUTINE lim_thd_lac

trunk/NEMO/LIM_SRC_3/limthd_sal.F90

@@ -1 +1 @@
 MODULE limthd_sal
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limthd_sal ***
@@ -20 +23 @@
    USE limvar
    USE par_ice
-   USE limicepoints
 
    IMPLICIT NONE
@@ -60 +62 @@
       !!
       !! ** History  : 
+      !!
+      !! "Je ne suis pour l'instant qu'a 80% de ma condition, mais c'est 
+      !!  les 30% qui restent qui seront les plus difficiles"
+      !!                                           E. Mpenza
+      !!
       !!-------------------------------------------------------------------
       !! History :
@@ -65 +72 @@
       !!   3.0  !  05-12 Routine rewritten for the 3-D version
       !!---------------------------------------------------------------------
+      !!
       !! * Local variables
       INTEGER, INTENT(in) :: &
@@ -74 +82 @@
 
       REAL(wp) ::            &
-!        ds_i_gd,          &     !: salt loss by gravity drainage
-!        ds_i_seg,         &     !: salt gain by segregation of seawater at the bottom
-!        ds_i_flush,       &     !: summer salt loss by flushing
-         zargtemp,           &   !: used to compute Scharzacher profile 
          zsold,              &   !: old salinity
-!        ds_i_snice,       &     !: increase of salinity due to snow ice
          zeps=1.0e-06   ,    &   !: very small
          iflush         ,    &   !: flushing (1) or not (0)
@@ -87 +90 @@
          i_ice_switch   ,    &   !: ice thickness above a certain treshold or not
          ztmelts        ,    &   !: freezing point of sea ice
-         zs_snowice     ,    &   !: salinity of forming snow ice
          zaaa           ,    &   !: dummy factor
          zbbb           ,    &   !: dummy factor
@@ -95 +97 @@
       REAL(wp), DIMENSION(jpij) ::          &
          ze_init        ,    &   !initial total enthalpy
-         ze_end         ,    &   !final total enthalpy
          zhiold         ,    &
          zsiold
@@ -157 +158 @@
             END DO ! ji
          END DO ! jk
-
- !       IF (jiindex_1d .GT. 0 ) THEN
- !       WRITE(numout,*) ' s_i_b   : ', s_i_b(jiindex_1d,1:nlay_i)
- !       WRITE(numout,*) ' t_i_b   : ', t_i_b(jiindex_1d,1:nlay_i)
- !       WRITE(numout,*) ' q_i_b   : ', q_i_b(jiindex_1d,1:nlay_i)
- !       ENDIF
 
          DO ji = kideb, kiut
@@ -215 +210 @@
          CALL lim_var_salprof1d(kideb,kiut)
 
-         !--------------------------------------
-         ! Energy of melting for each ice layer
-         !--------------------------------------
-         ! q=q(S,T) J.m-3
-!        DO jk = 1, nlay_i
-!           DO ji = kideb, kiut
-!           ! New heat content after desalination
-!              ztmelts      = - tmut * s_i_b(ji,jk) + rtt ! Tfi in K 
-!              q_i_b(ji,jk) = rhoic *                                         &
-!                           ( cpic * ( ztmelts-t_i_b(ji,jk) )                 &
-!                           + lfus * ( 1.0 - (ztmelts-rtt) /                  &
-!                             MIN( ( t_i_b(ji,jk) - rtt ) , -zeps ) )         &
-!                           - rcp  * ( ztmelts-rtt ) )
-!           END DO ! ji
-!        END DO
-
-         !--------------------
-         ! Total heat content
-         !--------------------
-!        ze_end(:) = 0.0
-!        DO jk = 1, nlay_i
-!           DO ji = kideb, kiut
-!              ze_end(ji)  = ze_end(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
-!           END DO 
-!        END DO
-            
-!        IF (jiindex_1d .GT. 0 ) THEN
-!           WRITE(numout,*) ' ze_init : ', ze_init(jiindex_1d) / rdt_ice
-!           WRITE(numout,*) ' ze_end  : ', ze_end (jiindex_1d) /  rdt_ice
-!           WRITE(numout,*) ' s_i_b   : ', s_i_b(jiindex_1d,1:nlay_i)
-!           WRITE(numout,*) ' t_i_b   : ', t_i_b(jiindex_1d,1:nlay_i)
-!           WRITE(numout,*) ' q_i_b   : ', q_i_b(jiindex_1d,1:nlay_i)
-!        ENDIF
-
          !----------------------------
          ! Heat flux - brine drainage
          !----------------------------
-
-!        WRITE(numout,*) ' jiindex, jjindex : ', jiindex, jjindex
-!        WRITE(numout,*) ' jiindex_1d : ', jiindex_1d
-!        WRITE(numout,*) ' kideb, kiut: ', kideb, kiut
 
          DO ji = kideb, kiut
@@ -265 +222 @@
             iaccrbo      =  MAX( 0.0 , SIGN ( 1.0 , dh_i_bott(ji) ) )
 
-         !  IF ( ( iflush .EQ. 1 )  .OR. ( iaccrbo .EQ. 1 ) ) THEN
-!           fhbri_1d(ji) = fhbri_1d(ji) + ( ze_end(ji) - ze_init(ji) ) * &
-!                          a_i_b(ji) / at_i_b(ji) / rdt_ice
-                           ! ze_end GT ze_init... ?
-         !  ENDIF 
             fhbri_1d(ji) = 0.0
-
-            !++++++
-!           zji                 = MOD( npb(ji) - 1, jpi ) + 1
-!           zjj                 = ( npb(ji) - 1 ) / jpi + 1
-!           fhbricat(zji,zjj,jl)  =   ( ze_end(ji) - ze_init(ji) ) / &
-!                                   rdt_ice
-            !++++++
-!           ftotal_fin(ji) = ftotal_fin(ji) + fhbricat(zji,zjj,jl) 
-!           IF ( ji .EQ. jiindex_1d) WRITE(numout,*) ' fhbri : ', &
-!                                  fhbricat(zji,zjj,jl) * rdt_ice
-!           IF ( ji .EQ. jiindex_1d ) &
-!              WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)
-!           !+++++
-
-         END DO ! ji
-
-!        IF (jiindex_1d .GT. 0 ) THEN
-!           WRITE(numout,*) ' Category : ', jl
-!           WRITE(numout,*) ' jiindex_1d : ', jiindex_1d
-!           WRITE(numout,*) ' fhbricat : ', fhbricat(jiindex,jjindex,jl)
-!           WRITE(numout,*) ' fhbri_1d : ', fhbri_1d(jiindex_1d)
-!           WRITE(numout,*) ' ze_end   : ', ze_end(jiindex_1d)
-!           WRITE(numout,*) ' ze_init  : ', ze_init(jiindex_1d)
-!           WRITE(numout,*) ' a_i_b    : ', a_i_b(jiindex_1d)
-!           WRITE(numout,*) ' at_i_b   : ', at_i_b(jiindex_1d)
-!           WRITE(numout,*) ' rdt_ice  : ', rdt_ice
-!        ENDIF
+         END DO ! ji
 
          !----------------------------
@@ -377 +303 @@
             ENDIF
 
-      !!--multiyear sea-ice salinity
-      !!-- ca va pas alors on vire bordel de djeuille
-!            if (ht.gt.360) then
-!           do jk = 1, nlay_i
-!           s_i_b(ji,jk) = 1.58 + 0.18*ht_i_b(ji)
-!           end do
-!        endif 
             DO jk = 1, nlay_i
                s_i_b(ji,jk)   = sm_i_b(ji)
@@ -415 +334 @@
          END DO ! ji
       ENDIF
-
-!+++++
-!           IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN
-!              WRITE(numout,*) ' *** FS_EQV *** '
-!              WRITE(numout,*) ' limthd_sal '
-!              WRITE(numout,*) ' fseqv      ', fseqv_1d(ji)
-!              WRITE(numout,*) ' sss_io     ', sss_io(zji,zjj)
-!              WRITE(numout,*) ' s_i_new    ', s_i_new(ji)
-!              WRITE(numout,*) ' dh_i_bott  ', dh_i_bott(ji)
-!              WRITE(numout,*) ' a_i_b      ', a_i_b(ji)
-!              WRITE(numout,*) ' *** FS_BRI *** '
-!              WRITE(numout,*) ' fsbri      ', fsbri_1d(ji)
-!              WRITE(numout,*) ' ht_i       ', ht_i_b(ji)
-!              WRITE(numout,*) ' dsm_i_gd   ', dsm_i_gd_1d(ji)
-!              WRITE(numout,*) ' dsm_i_fl   ', dsm_i_fl_1d(ji)
-!           ENDIF
-
-!+++++
 
 !-- End of salinity computations

trunk/NEMO/LIM_SRC_3/limtrp.F90

@@ -6 +6 @@
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_trp      : advection/diffusion process of sea ice
@@ -26 +26 @@
    USE lib_mpp
    USE par_ice
-   USE limicepoints
 
    IMPLICIT NONE
@@ -50 +49 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
+   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limtrp.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
@@ -75 +74 @@
       !!---------------------------------------------------------------------
       !! * Local Variables
-      INTEGER  ::   ji, jj, jk, jl, index, layer, &  ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, layer, &  ! dummy loop indices
                     initad           ! number of sub-timestep for the advection
       INTEGER  ::   ji_maxu, ji_maxv, jj_maxu, jj_maxv
@@ -81 +80 @@
       REAL(wp) ::  &                              
          zindb  ,  &
-         zacrith, &
-         zindsn , &
-         zindic , &
-         zusvosn, &
-         zusvoic, &
-         zignm  , &
-         zindhe , &
-         zvbord , &
-         zcfl   , &
-         zusnit , &
-         zrtt, ztsn, ztic1, ztic2, zsal, zage, &
-         zq, zbigval, ze1, ze2, ze, &
+         zindsn ,  &
+         zindic ,  &
+         zusvosn,  &
+         zusvoic,  &
+         zvbord ,  &
+         zcfl   ,  &
+         zusnit ,  &
+         zrtt, zsal, zage, &
+         zbigval, ze, &
          zmaxu, zmaxv
 
@@ -101 +97 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl):: &  ! temporary workspace
          zs0ice, zs0sn, zs0a   ,         &
-         zs0c0 , zs0c1 , zs0c2 ,         &
+         zs0c0 ,                         &
          zs0sm , zs0oi
 
@@ -554 +550 @@
          END DO ! jl
 
-!-----------------------------------------------------------------------------!
-! 6) Correct age 
-!-----------------------------------------------------------------------------!
-
-!        DO jl = 1, jpl
-!           DO jj = 1, jpj
-!              DO ji = 1, jpi
-!                 IF (old_a_i(ji,jj,jl).eq.0.00) THEN
-!                     o_i(ji,jj,jl) = MAX( MIN( rdt_ice*float(numit)/86400.0, o_i(ji,jj,jl) ), 0.0 )
-!                 ENDIF
-!              END DO
-!           END DO
-!        END DO
-
-
       ENDIF
 

trunk/NEMO/LIM_SRC_3/limupdate.F90

@@ -9 +9 @@
    !!      from advection and ice thermodynamics 
    !!======================================================================
+#if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM sea-ice model
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    !!    lim_update   : computes update of sea-ice global variables 
@@ -42 +43 @@
    USE thd_ice         ! LIM thermodynamic sea-ice variables
    USE par_ice
-   USE limicepoints
    USE limitd_th
    USE limvar
@@ -80 +80 @@
       INTEGER ::      &
           ji, jj,     & ! geographical indices
-          jk, jl, jm, & ! layer, category and type indices
-          index, num_ex_cat, zji, zjj  ! dummy loop indices
+          jk, jl, jm    ! layer, category and type indices
       INTEGER ::      &
           jbnd1, jbnd2
@@ -96 +95 @@
          rzero  = 0.e0    ,  &
          rone   = 1.e0    ,  &
-         z_slope_s        ,  &    ! slope of the salinity profile
          zhimax                   ! maximum thickness tolerated for advection of
                                   ! in an ice-free cell
       REAL(wp) ::            &  ! dummy switches and arguments
-         zindhe, zconc, zindb, zindsn, zindic, zacrith, zignm,  &
-         zrtt, ztsn, zindg, zthsnice, za, zh, zdvres, zviold,ztic, &
-         zbigvalue, zvsold, zvold, z_dv_ex, z_da_ex, zamax, zdv, &
-         z_prescr_hi, zat_i_old, z_a_i, zalpha, zs_inf, zsmin, zsmax, &
+         zindb, zindsn, zindic, zacrith,  &
+         zrtt, zindg, zh, zdvres, zviold,                       &
+         zbigvalue, zvsold, z_da_ex, zamax,                     &
+         z_prescr_hi, zat_i_old,                                &
          ztmelts, ze_s
-
-      REAL(wp), DIMENSION(nlay_i) ::   &                                
-         z_de, zs_zero
 
       REAL(wp), DIMENSION(jpl) :: z_da_i, z_dv_i
@@ -113 +108 @@
       LOGICAL, DIMENSION(jpi,jpj,jpl) ::  &
          internal_melt
+
       INTEGER ::      &
          ind_im, layer      ! indices for internal melt
@@ -734 +730 @@
                  ( 1.0 - zindsn ) * 0.0
 
-!             ztsn            = MIN( MAX( zrtt, t_s(ji,jj,1,jl) ) , rtt )
-!             t_s(ji,jj,1,jl) = zindsn*ztsn + (1.0 - zindsn)*t_bo(ji,jj)
-
            END DO ! ji
         END DO ! jj
@@ -864 +857 @@
 
            ! 1) Count the number of existing categories
-!          num_ex_cat = 0
            DO jl  = 1, jpl
               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) ) 
@@ -1005 +997 @@
 
    END SUBROUTINE lim_update
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty Module               No sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_update     ! Empty routine
+   END SUBROUTINE lim_update
+
+#endif
+
 END MODULE limupdate

trunk/NEMO/LIM_SRC_3/limvar.F90

@@ -1 +1 @@
 MODULE limvar
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                       ***  MODULE limvar ***
@@ -38 +41 @@
    USE ice_oce          ! ice variables
    USE thd_ice
-   USE limicepoints
    USE in_out_manager
    USE ice
@@ -62 +64 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM-@ 4.0,  UCL-ASTR (2006)
+   !!   LIM3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
    !!   (c) UCL-ASTR and Martin Vancoppenolle
    !!----------------------------------------------------------------------
@@ -91 +93 @@
         !! ** History :
         !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
-        !!                     Take it easy man
-        !!                     Life is just a simple game, between ups / and downs \ :@)
         !!
         !! note : you could add an argument when you need only at_i, vt_i
@@ -103 +103 @@
                      jj,       &   ! spatial dummy loop index
                      jk,       &   ! vertical layering dummy loop index
-                     jl,       &   ! ice category dummy loop index
-                     jm,       &   ! ice types    dummy loop index
-                     index
+                     jl            ! ice category dummy loop index
  
         REAL ::      zeps, epsi16, zinda, epsi06
@@ -116 +114 @@
        epsi16 = 1.0e-16
        epsi06 = 1.0e-6
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_var_agg : Aggregate ice variables '
-         WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
 
        !------------------
@@ -228 +220 @@
                      jj,       &   ! spatial dummy loop index
                      jk,       &   ! vertical layering dummy loop index
-                     jl,       &   ! ice category dummy loop index
-                     jm,       &   ! ice types    dummy loop index
-                     index
+                     jl            ! ice category dummy loop index
 
         REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, &
@@ -242 +232 @@
 !!-- End of declarations
 !!------------------------------------------------------------------------------
-!     IF (lwp) THEN
-!     WRITE(numout,*)
-!     WRITE(numout,*) 'lim_var_glo2eqv : transform global variables into equivalent ice variables '
-!     WRITE(numout,*) '~~~~~~~~~~~~~~~'
-!     ENDIF
 
       !-------------------------------------------------------
@@ -377 +362 @@
         !!
         !!------------------------------------------------------------------
-        !! * Arguments
-
-        !! * Local variables
-        INTEGER ::   ji,       &   ! spatial dummy loop index
-                     jj,       &   ! spatial dummy loop index
-                     jk,       &   ! vertical layering dummy loop index
-                     jl,       &   ! ice category dummy loop index
-                     jm,       &   ! ice types    dummy loop index
-                     index
-
-        REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, &
-                ztmelts
-
-        REAL :: zeps
-
-        !Ice specific heat, massive latent heat
-        zeps    = 1.0e-10
-
-!!-- End of declarations
-!!------------------------------------------------------------------------------
-
-!      IF (lwp) THEN
-!      WRITE(numout,*)
-!      WRITE(numout,*) 'lim_var_eqv2glo : transform equivalent variables into global ice variables '
-!      WRITE(numout,*) '~~~~~~~~~~~~~~~'
-!      ENDIF
 
        v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
@@ -408 +367 @@
        smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
        oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
-
-!      DO jl = 1, jpl
-!        DO jk = 1, nlay_i
-!          DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               ! heat content per unit volume
-!               e_i(ji,jj,jk,jl) = zidto * rhoic * &
-!               (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-!               + lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
-!               - rcp      * ( ztmelts - rtt ) &
-!               )
-
-!               ! Correct dimensions to avoid big values
-!                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
-
-!               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-!                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
-!                                    area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
-!                                    FLOAT(nlay_i)
-
-!           END DO
-!         END DO
-!       END DO
-!     END DO
 
     END SUBROUTINE lim_var_eqv2glo
@@ -787 +722 @@
 
 #endif
- END MODULE limvar
+END MODULE limvar

trunk/NEMO/LIM_SRC_3/limwri.F90

@@ -1 +1 @@
 MODULE limwri
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
    !!======================================================================
    !!                     ***  MODULE  limwri  ***
@@ -8 +12 @@
    !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limwri.F90,v 1.4 2005/03/27 18:34:42 opalod Exp $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_wri      : write of the diagnostics variables in ouput file 
@@ -29 +29 @@
    USE iceini
    USE lbclnk
-   USE limicepoints
    USE par_ice
    USE limvar
@@ -98 +97 @@
          
       REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zfield, zfielda
+         zfield
 
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::  &
@@ -107 +106 @@
       CHARACTER(len = 40)  :: &
          clhstnam, clop, &
-         clhstnama, clopa
+         clhstnama
 
       INTEGER , SAVE ::      &
          nice, nhorid, ndim, niter, ndepid
       INTEGER , SAVE ::      &
-         nicea, nhorida, ndima, nitera, ndepida, ndimitd
+         nicea, nhorida, nitera, ndimitd
       INTEGER , DIMENSION( jpij ) , SAVE ::  &
          ndex51

trunk/NEMO/LIM_SRC_3/thd_ice.F90

@@ -17 +17 @@
    PRIVATE
 
+   !!---------------------------
    !! * Share Module variables
+   !!---------------------------
    REAL(wp) , PUBLIC ::   & !!! ** ice-thermo namelist (namicethd) **
-      hmelt   = -0.15  ,  &  !: maximum melting at the bottom
-      hicmin  = 0.2    ,  &  !: ice th. corr. to max. ener. in brine pocket
+      hmelt   = -0.15  ,  &  !: maximum melting at the bottom; active only for
+                             !: one category
+      hicmin  = 0.2    ,  &  !: (REMOVE)
       hiclim  = 0.05   ,  &  !: minimum ice thickness
       amax    = 0.999  ,  &  !: maximum lead fraction
-      sbeta   = 1.0    ,  &  !: numerical scheme for diffusion in ice 
-      parlat  = 0.0    ,  &  !: percent. of energy used for lateral ablation
-      hakspl  = 0.5    ,  &  !: slope of distr. for Hakkinen-Mellro's lat. melt
-      hibspl  = 0.5    ,  &  !: slope of distribution for Hibler's lat. melt
-      exld    = 2.0    ,  &  !: exponent for leads-closure rate
-      hakdif  = 1.0    ,  &  !: coefficient for diffusions of ice and snow
-      thth    = 0.2    ,  &  !: thick. for comp. of eq. thermal conduct
+      sbeta   = 1.0    ,  &  !: numerical scheme for diffusion in ice  (REMOVE)
+      parlat  = 0.0    ,  &  !: (REMOVE)
+      hakspl  = 0.5    ,  &  !: (REMOVE)
+      hibspl  = 0.5    ,  &  !: (REMOVE)
+      exld    = 2.0    ,  &  !: (REMOVE)
+      hakdif  = 1.0    ,  &  !: (REMOVE)
+      thth    = 0.2    ,  &  !: (REMOVE)
       hnzst   = 0.1    ,  &  !: thick. of the surf. layer in temp. comp.
       parsub  = 1.0    ,  &  !: switch for snow sublimation or not
@@ -41 +44 @@
       hiccrit = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m)
 
-   REAL(wp) , PUBLIC ::   &  !:
-      uscomi                 !: inverse of minimum lead fraction
+   !!-----------------------------
+   !! * Share 1D Module variables
+   !!-----------------------------
+   !: In ice thermodynamics, to spare memory, the vectors are folded
+   !: from 1D to 2D vectors. The following variables, with ending _1d (or _b)
+   !: are the variables corresponding to 2d vectors
 
    INTEGER , PUBLIC, DIMENSION(jpij) ::   &  !:
       npb     ,   &   !: number of points where computations has to be done
-      npac            !: correspondance between the points
+      npac            !: correspondance between the points (lateral accretion)
 
    REAL(wp), PUBLIC, DIMENSION(jpij) ::   &  !: 
@@ -123 +130 @@
       q_s_b                !:    Snow enthalpy per unit volume
 
-  ! goes to trash, just for conservation checks
-   REAL(wp), PUBLIC, DIMENSION(jpij,jpl)    ::   &  !: ! goes to trash
+  ! Clean the following ...
+  ! These variables are coded for conservation checks
+   REAL(wp), PUBLIC, DIMENSION(jpij,jpl)    ::   &  !
       qt_i_in   ,           &  !: ice energy summed over categories (initial)
       qt_i_fin  ,           &  !: ice energy summed over categories (final)