New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 2601 – NEMO

Changeset 2601


Ignore:
Timestamp:
2011-02-20T16:30:58+01:00 (13 years ago)
Author:
gm
Message:

dynamic mem: #785 ; LIM-3 case: changes required for compilation

Location:
branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90

    r2590 r2601  
    66   !! History :  3.0  ! 2003-08  (M. Vancoppenolle)  LIM-3 
    77   !!---------------------------------------------------------------------- 
    8    USE par_ice 
     8   USE par_ice        ! LIM-3 parameter 
     9   USE in_out_manager ! I/O manager 
    910 
    1011   IMPLICIT NONE 
     
    3132   !! $Id$ 
    3233   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    33    !!====================================================================== 
     34   !!---------------------------------------------------------------------- 
    3435CONTAINS 
    3536 
     
    4041      INTEGER :: dom_ice_alloc 
    4142      !!------------------------------------------------------------------- 
    42  
    43       ALLOCATE(fs2cor(jpi,jpj), fcor(jpi,jpj), & 
    44                covrai(jpi,jpj), area(jpi,jpj), & 
    45                tms(jpi,jpj)   , tmi(jpi,jpj) , & 
    46                tmu(jpi,jpj)   , tmv(jpi,jpj) , & 
    47                tmf(jpi,jpj)   ,                & 
    48                wght(jpi,jpj,2,2), Stat = dom_ice_alloc) 
    49  
    50       IF(dom_ice_alloc /= 0)THEN 
    51          CALL ctl_warn('dom_ice_alloc: failed to allocate arrays.') 
     43      ! 
     44      ALLOCATE( fs2cor(jpi,jpj) , fcor(jpi,jpj) ,      & 
     45         &      covrai(jpi,jpj) , area(jpi,jpj) ,      & 
     46         &      tms   (jpi,jpj) , tmi (jpi,jpj) ,      & 
     47         &      tmu   (jpi,jpj) , tmv (jpi,jpj) ,      & 
     48         &      tmf   (jpi,jpj) ,                      & 
     49         &      wght(jpi,jpj,2,2)               , STAT = dom_ice_alloc ) 
     50      ! 
     51      IF( dom_ice_alloc /= 0 ) THEN 
     52         CALL ctl_warn( 'dom_ice_alloc: failed to allocate arrays.' ) 
    5253      END IF 
    53  
     54      ! 
    5455   END FUNCTION dom_ice_alloc 
    5556 
     57   !!====================================================================== 
    5658END MODULE dom_ice 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r2528 r2601  
    44   !! LIM-3 Sea Ice physics:  diagnostics variables of ice defined in memory 
    55   !!===================================================================== 
    6    !! History :  3.0  ! 2008-03  (M. Vancoppenolle) : original code LIM-3 
     6   !! History :  3.0  ! 2008-03  (M. Vancoppenolle) original code LIM-3 
     7   !!            4.0  ! 3011-02  (G. Madec) dynamical allocation 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim3 
     
    1112   !!---------------------------------------------------------------------- 
    1213   USE par_ice          ! LIM sea-ice parameters 
     14   USE in_out_manager   ! I/O manager 
    1315 
    1416   IMPLICIT NONE 
    1517   PRIVATE 
     18 
     19   PUBLIC    ice_alloc  !  Called in iceini.F90 
    1620 
    1721   !!====================================================================== 
    1822   !! LIM3 by the use of sweat, agile fingers and sometimes brain juice,  
    1923   !!  was developed in Louvain-la-Neuve by :  
    20    !! 
    2124   !!    * Martin Vancoppenolle (UCL-ASTR, Belgium) 
    2225   !!    * Sylvain Bouillon (UCL-ASTR, Belgium) 
    23    !!    * Miguel Angel Morales Maqueda (POL, UK) 
     26   !!    * Miguel Angel Morales Maqueda (NOC-L, UK) 
    2427   !!  
    2528   !! Based on extremely valuable earlier work by 
    26    !! 
    2729   !!    * Thierry Fichefet 
    2830   !!    * Hugues Goosse 
    2931   !! 
    3032   !! The following persons also contributed to the code in various ways 
    31    !! 
    32    !!    * Gurvan Madec, Claude Talandier, Christian Ethe  
    33    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     33   !!    * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France) 
    3434   !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany) 
    3535   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa)  
    3636   !!      and Elisabeth Hunke (LANL), USA. 
    3737   !!  
    38    !! (c) UCL-ASTR, 2005-2008 
    39    !! 
    40    !! For more info, the interested user is kindly invited to consult the 
    41    !! following references 
     38   !! For more info, the interested user is kindly invited to consult the following references 
    4239   !!    For model description and validation : 
    4340   !!    * Vancoppenolle et al., Ocean Modelling, 2008a. 
    4441   !!    * Vancoppenolle et al., Ocean Modelling, 2008b. 
     42   !!    For a specific description of EVP : 
     43   !!    * Bouillon et al., Ocean Modelling 2009. 
    4544   !! 
    46    !!    For a specific description of EVP : 
    47    !!    * Bouillon et al., in prep for 2008. 
    48    !! 
    49    !!    Or the reference manual, that should be available by 2010 
     45   !!    Or the reference manual, that should be available by 2011 
    5046   !!====================================================================== 
    5147   !!                                                                     | 
     
    168164   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step 
    169165 
    170    INTEGER , PUBLIC ::   &     !!: ** ice-dynamic namelist (namicedyn) ** 
    171       nbiter = 1      ,  &     !: number of sub-time steps for relaxation 
    172       nbitdr = 250    ,  &     !: maximum number of iterations for relaxation 
    173       nevp   = 400    ,  &     !: number of iterations for subcycling 
    174       nlay_i = 5               !: number of layers in the ice 
    175  
    176    REAL(wp), PUBLIC ::   &     !!: ** ice-dynamic namelist (namicedyn) ** 
    177       epsd   = 1.0e-20,  &     !: tolerance parameter for dynamic 
    178       alpha  = 0.5    ,  &     !: coefficient for semi-implicit coriolis 
    179       dm     = 0.6e+03,  &     !: diffusion constant for dynamics 
    180       om     = 0.5    ,  &     !: relaxation constant 
    181       resl   = 5.0e-05,  &     !: maximum value for the residual of relaxation 
    182       cw     = 5.0e-03,  &     !: drag coefficient for oceanic stress 
    183       angvg  = 0.0    ,  &     !: turning angle for oceanic stress 
    184       pstar  = 1.0e+04,  &     !: determines ice strength (N/M), Hibler JPO79 
    185       c_rhg  = 20.0   ,  &     !: determines changes in ice strength 
    186       etamn  = 0.0e+07,  &     !: minimun value for viscosity : has to be 0 
    187       creepl = 2.0e-08,  &     !: creep limit : has to be under 1.0e-9 
    188       ecc    = 2.0    ,  &  !: eccentricity of the elliptical yield curve 
    189       ahi0   = 350.e0 ,  &  !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    190       telast = 2880.0 ,  &  !: timescale for elastic waves (s) !SB 
    191       alphaevp = 1.0  ,  &  !: coeficient of the internal stresses !SB 
    192       unit_fac = 1.0e9      !: conversion factor for ice / snow enthalpy 
    193  
    194    REAL(wp), PUBLIC ::   & !!: ** ice-salinity namelist (namicesal) ** 
    195       s_i_max  = 20.0 ,  &  !: maximum ice salinity (ppt) 
    196       s_i_min  =  0.1 ,  &  !: minimum ice salinity (ppt) 
    197       s_i_0    =  3.5 ,  &  !: 1st sal. value for the computation of sal .prof. 
    198                                 !: (ppt) 
    199       s_i_1    =  4.5 ,  &  !: 2nd sal. value for the computation of sal .prof. 
    200                                 !: (ppt) 
    201       sal_G    = 5.00 ,  &  !: restoring salinity for gravity drainage 
    202                                 !: (ppt) 
    203       sal_F    = 2.50 ,  &  !: restoring salinity for flushing 
    204                                 !: (ppt) 
    205       time_G   = 1.728e+06,&!: restoring time constant for gravity drainage 
    206                                 !: (= 20 days, in s) 
    207       time_F   = 8.640e+05,&!: restoring time constant for gravity drainage  
    208                                 !: (= 10 days, in s) 
    209       bulk_sal = 4.0        !: bulk salinity (ppt) in case of constant salinity 
    210  
    211    INTEGER , PUBLIC ::   & !!: ** ice-salinity namelist (namicesal) ** 
    212       num_sal  = 1    ,  &  !: salinity configuration used in the model 
    213                                 !: 1 - s constant in space and time 
    214                                 !: 2 - prognostic salinity (s(z,t)) 
    215                                 !: 3 - salinity profile, constant in time 
    216                                 !: 4 - salinity variations affect only ice 
    217                                 !      thermodynamics 
    218       sal_prof = 1    ,  &  !: salinity profile or not  
    219       thcon_i_swi = 1       !: thermal conductivity of Untersteiner (1964) (1) or 
    220    !: Pringle et al (2007) (2) 
    221  
    222    REAL(wp), PUBLIC ::   & !!: ** ice-mechanical redistribution namelist (namiceitdme) 
    223       Cs = 0.25       ,  & !!: fraction of shearing energy contributing to ridging             
    224       Cf = 17.0       ,  & !!: ratio of ridging work to PE loss 
    225       fsnowrdg = 0.5  ,  & !!: fractional snow loss to the ocean during ridging 
    226       fsnowrft = 0.5  ,  & !!: fractional snow loss to the ocean during ridging 
    227       Gstar = 0.15    ,  & !!: fractional area of young ice contributing to ridging 
    228       astar = 0.05    ,  & !!: equivalent of G* for an exponential participation function 
    229       Hstar = 100.0   ,  & !!: thickness that determines the maximal thickness of ridged 
    230                                 !!: ice 
    231       hparmeter = 0.75,  & !!: threshold thickness (m) for rafting / ridging  
    232       Craft = 5.0     ,  & !!: coefficient for smoothness of the hyperbolic tangent in rafting 
    233       ridge_por = 0.0 ,  & !!: initial porosity of ridges (0.3 regular value) 
    234       sal_max_ridge = 15.0, & !!: maximum ridged ice salinity (ppt) 
    235       betas    = 1.0      , & !:: coef. for partitioning of snowfall between leads and sea ice 
    236       kappa_i  = 1.0      , & !!: coefficient for the extinction of radiation 
    237                                 !!: Grenfell et al. (2006) (m-1) 
    238       nconv_i_thd = 50    , & !!: maximal number of iterations for heat diffusion 
    239       maxer_i_thd = 1.0e-4    !!: maximal tolerated error (C) for heat diffusion 
    240  
    241    INTEGER , PUBLIC ::   & !!: ** ice-mechanical redistribution namelist (namiceitdme) 
    242       ridge_scheme_swi = 0, & !!: scheme used for ice ridging 
    243       raftswi          = 1, & !!: rafting of ice or not                         
    244       partfun_swi      = 1, & !!: participation function Thorndike et al. JGR75 (0)  
    245                                 !!: or Lipscomb et al. JGR07 (1)  
    246       transfun_swi     = 0, & !!: transfer function of Hibler, MWR80 (0)  
    247                                 !!: or Lipscomb et al., 2007 (1) 
    248       brinstren_swi    = 0    !!: use brine volume to diminish ice strength 
    249  
    250    REAL(wp), PUBLIC ::   &  !: 
    251       usecc2          ,  &  !:  = 1.0 / ( ecc * ecc ) 
    252       rhoco           ,  &  !: = rau0 * cw 
    253       sangvg, cangvg  ,  &  !: sin and cos of the turning angle for ocean stress 
    254       pstarh                !: pstar / 2.0 
    255  
    256    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !: 
    257       u_oce, v_oce    ,  &  !: surface ocean velocity used in ice dynamics 
    258       ahiu , ahiv     ,  &  !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    259       pahu , pahv     ,  &  !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    260       ust2s, hicol    ,  &  !: friction velocity, ice collection thickness accreted in leads 
    261       strength        ,  &  !: ice strength 
    262       strp1, strp2    ,  &  !: strength at previous time steps 
    263       stress1_i       ,  &  !: first stress tensor element 
    264       stress2_i       ,  &  !: second stress tensor element 
    265       stress12_i      ,  &  !: diagonal stress tensor element 
    266       delta_i         ,  &  !: Delta factor for the ice rheology (see Flato and Hibler 95) [s-1] -> limrhg.F90 
    267       divu_i          ,  &  !: Divergence of the velocity field [s-1] -> limrhg.F90 
    268       shear_i               !: Shear of the velocity field [s-1] -> limrhg.F90 
    269  
    270    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    271       firic  ,   &  !: IR flux over the ice (only used for outputs) 
    272       fcsic  ,   &  !: Sensible heat flux over the ice (only used for outputs) 
    273       fleic  ,   &  !: Latent heat flux over the ice (only used for outputs) 
    274       qlatic ,   &  !: latent flux 
    275       rdvosif,   &  !: Variation of volume at surface (only used for outputs) 
    276       rdvobif,   &  !: Variation of ice volume at the bottom ice (only used for outputs) 
    277       fdvolif,   &  !: Total variation of ice volume (only used for outputs) 
    278       rdvonif,   &  !: Lateral Variation of ice volume (only used for outputs) 
    279       sist   ,   &  !: Average Sea-Ice Surface Temperature (Kelvin) 
    280       icethi ,   &  !: total ice thickness (for all categories) (only used for outputs) 
    281       t_bo   ,   &  !: Sea-Ice bottom temperature (Kelvin)       
    282       hicifp ,   &  !: Ice production/melting 
    283                                 !obsolete... can be removed 
    284       frld   ,   &  !: Leads fraction = 1-a/totalarea REFERS TO LEAD FRACTION everywhere 
    285                                 !: except in the OUTPUTS!!!! 
    286       pfrld  ,   &  !: Leads fraction at previous time   
    287       phicif ,   &  !: Old ice thickness 
    288       fbif   ,   &  !: Heat flux at the ice base 
    289       rdmsnif,   &  !: Variation of snow mass 
    290       rdmicif,   &  !: Variation of ice mass 
    291       qldif  ,   &  !: heat balance of the lead (or of the open ocean) 
    292       qcmif  ,   &  !: Energy needed to bring the ocean surface layer until its freezing  
    293       fdtcn  ,   &  !: net downward heat flux from the ice to the ocean 
    294       qdtcn  ,   &  !: energy from the ice to the ocean 
    295       fstric ,   &  !: transmitted solar radiation under ice 
    296       fscmbq ,   &  !: associated with lead chipotage with solar flux 
    297       ffltbif,   &  !: Array linked with the max heat contained in brine pockets (?) 
    298       fsbbq  ,   &  !: Also linked with the solar flux below the ice (?) 
    299       qfvbq  ,   &  !: Array used to store energy in case of toral lateral ablation (?) 
    300       dmgwi  ,   &  !: Variation of the mass of snow ice 
    301       fsalt_res, &  !: Residual salt flux due to correction of ice thickness 
    302       fsbri  ,   &  !: Salt flux due to brine rejection 
    303       fsalt_rpo, &  !: Salt flux associated with porous ridged ice formation 
    304       fheat_rpo, &  !: Heat flux associated with porous ridged ice formation 
    305       fhbri  ,   &  !: heat flux due to brine rejection 
    306       fmmec  ,   &  !: Mass flux due to snow loss during compression 
    307       fseqv  ,   &  !: Equivalent salt flux due to ice growth/melt 
    308       fheat_res, &  !: Residual heat flux due to correction of ice thickness 
    309       fhmec         !: Heat flux due to snow loss during compression 
     166   !                                    !!** ice-dynamic namelist (namicedyn) ** 
     167   INTEGER , PUBLIC ::   nbiter = 1      !: number of sub-time steps for relaxation 
     168   INTEGER , PUBLIC ::   nbitdr = 250    !: maximum number of iterations for relaxation 
     169   INTEGER , PUBLIC ::   nevp   = 400    !: number of iterations for subcycling 
     170   INTEGER , PUBLIC ::   nlay_i = 5      !: number of layers in the ice 
     171 
     172   !                                          !!** ice-dynamic namelist (namicedyn) ** 
     173   REAL(wp), PUBLIC ::   epsd   = 1.0e-20_wp   !: tolerance parameter for dynamic 
     174   REAL(wp), PUBLIC ::   alpha  = 0.5_wp       !: coefficient for semi-implicit coriolis 
     175   REAL(wp), PUBLIC ::   dm     = 0.6e+03_wp   !: diffusion constant for dynamics 
     176   REAL(wp), PUBLIC ::   om     = 0.5_wp       !: relaxation constant 
     177   REAL(wp), PUBLIC ::   resl   = 5.0e-05_wp   !: maximum value for the residual of relaxation 
     178   REAL(wp), PUBLIC ::   cw     = 5.0e-03_wp   !: drag coefficient for oceanic stress 
     179   REAL(wp), PUBLIC ::   angvg  = 0._wp        !: turning angle for oceanic stress 
     180   REAL(wp), PUBLIC ::   pstar  = 1.0e+04_wp   !: determines ice strength (N/M), Hibler JPO79 
     181   REAL(wp), PUBLIC ::   c_rhg  = 20._wp       !: determines changes in ice strength 
     182   REAL(wp), PUBLIC ::   etamn  = 0.0e+07_wp   !: minimun value for viscosity : has to be 0 
     183   REAL(wp), PUBLIC ::   creepl = 2.0e-08_wp   !: creep limit : has to be under 1.0e-9 
     184   REAL(wp), PUBLIC ::   ecc    = 2._wp        !: eccentricity of the elliptical yield curve 
     185   REAL(wp), PUBLIC ::   ahi0   = 350._wp      !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
     186   REAL(wp), PUBLIC ::   telast = 2880._wp     !: timescale for elastic waves (s) !SB 
     187   REAL(wp), PUBLIC ::   alphaevp = 1._wp      !: coeficient of the internal stresses !SB 
     188   REAL(wp), PUBLIC ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy 
     189 
     190   !                                              !!** ice-salinity namelist (namicesal) ** 
     191   REAL(wp), PUBLIC ::   s_i_max  = 20.0_wp        !: maximum ice salinity [PSU] 
     192   REAL(wp), PUBLIC ::   s_i_min  =  0.1_wp        !: minimum ice salinity [PSU] 
     193   REAL(wp), PUBLIC ::   s_i_0    =  3.5_wp        !: 1st sal. value for the computation of sal .prof. [PSU] 
     194   REAL(wp), PUBLIC ::   s_i_1    =  4.5_wp        !: 2nd sal. value for the computation of sal .prof. [PSU] 
     195   REAL(wp), PUBLIC ::   sal_G    =  5.0_wp        !: restoring salinity for gravity drainage [PSU] 
     196   REAL(wp), PUBLIC ::   sal_F    =  2.5_wp        !: restoring salinity for flushing [PSU] 
     197   REAL(wp), PUBLIC ::   time_G   =  1.728e+06_wp  !: restoring time constant for gravity drainage (= 20 days) [s] 
     198   REAL(wp), PUBLIC ::   time_F   =  8.640e+05_wp  !: restoring time constant for gravity drainage (= 10 days) [s] 
     199   REAL(wp), PUBLIC ::   bulk_sal =  4.0_wp        !: bulk salinity (ppt) in case of constant salinity 
     200 
     201   !                                      !!** ice-salinity namelist (namicesal) ** 
     202   INTEGER , PUBLIC ::   num_sal     = 1   !: salinity configuration used in the model 
     203   !                                       ! 1 - s constant in space and time 
     204   !                                       ! 2 - prognostic salinity (s(z,t)) 
     205   !                                       ! 3 - salinity profile, constant in time 
     206   !                                       ! 4 - salinity variations affect only ice thermodynamics 
     207   INTEGER , PUBLIC ::   sal_prof    = 1   !: salinity profile or not  
     208   INTEGER , PUBLIC ::   thcon_i_swi = 1   !: thermal conductivity: =1 Untersteiner (1964) ; =2 Pringle et al (2007) 
     209 
     210   !                                              !!** ice-mechanical redistribution namelist (namiceitdme) 
     211   REAL(wp), PUBLIC ::   Cs        = 0.25_wp       !: fraction of shearing energy contributing to ridging             
     212   REAL(wp), PUBLIC ::   Cf        = 17.0_wp       !: ratio of ridging work to PE loss 
     213   REAL(wp), PUBLIC ::   fsnowrdg  = 0.5_wp        !: fractional snow loss to the ocean during ridging 
     214   REAL(wp), PUBLIC ::   fsnowrft  = 0.5_wp        !: fractional snow loss to the ocean during ridging 
     215   REAL(wp), PUBLIC ::   Gstar     = 0.15_wp       !: fractional area of young ice contributing to ridging 
     216   REAL(wp), PUBLIC ::   astar     = 0.05_wp       !: equivalent of G* for an exponential participation function 
     217   REAL(wp), PUBLIC ::   Hstar     = 100.0_wp      !: thickness that determines the maximal thickness of ridged ice 
     218   REAL(wp), PUBLIC ::   hparmeter = 0.75_wp       !: threshold thickness (m) for rafting / ridging  
     219   REAL(wp), PUBLIC ::   Craft     = 5.0_wp        !: coefficient for smoothness of the hyperbolic tangent in rafting 
     220   REAL(wp), PUBLIC ::   ridge_por = 0.0_wp        !: initial porosity of ridges (0.3 regular value) 
     221   REAL(wp), PUBLIC ::   sal_max_ridge = 15.0_wp   !: maximum ridged ice salinity (ppt) 
     222   REAL(wp), PUBLIC ::   betas    = 1.0_wp         !: coef. for partitioning of snowfall between leads and sea ice 
     223   REAL(wp), PUBLIC ::   kappa_i  = 1.0_wp         !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
     224   REAL(wp), PUBLIC ::   nconv_i_thd = 50_wp       !: maximal number of iterations for heat diffusion 
     225   REAL(wp), PUBLIC ::   maxer_i_thd = 1.0e-4_wp   !: maximal tolerated error (C) for heat diffusion 
     226 
     227   !                                           !!** ice-mechanical redistribution namelist (namiceitdme) 
     228   INTEGER , PUBLIC ::   ridge_scheme_swi = 0   !: scheme used for ice ridging 
     229   INTEGER , PUBLIC ::   raftswi          = 1   !: rafting of ice or not                         
     230   INTEGER , PUBLIC ::   partfun_swi      = 1   !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 
     231   INTEGER , PUBLIC ::   transfun_swi     = 0   !: transfer function: =0 Hibler 1980, =1 Lipscomb et al. 2007 
     232   INTEGER , PUBLIC ::   brinstren_swi    = 0   !: use brine volume to diminish ice strength 
     233 
     234   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( ecc * ecc ) 
     235   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cw 
     236   REAL(wp), PUBLIC ::   sangvg, cangvg   !: sin and cos of the turning angle for ocean stress 
     237   REAL(wp), PUBLIC ::   pstarh           !: pstar / 2.0 
     238 
     239   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce, v_oce   !: surface ocean velocity used in ice dynamics 
     240   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahiu , ahiv    !: hor. diffusivity coeff. at U- and V-points [m2/s] 
     241   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pahu , pahv    !: ice hor. eddy diffusivity coef. at U- and V-points 
     242   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ust2s, hicol   !: friction velocity, ice collection thickness accreted in leads 
     243   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strp1, strp2   !: strength at previous time steps 
     244   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength       !: ice strength 
     245   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element 
     246   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   delta_i        !: ice rheology elta factor (Flato & Hibler 95) [s-1] 
     247   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   divu_i         !: Divergence of the velocity field [s-1] 
     248   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i        !: Shear of the velocity field [s-1] 
     249   ! 
     250   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   firic       !: IR flux over the ice (diag only) 
     251   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fcsic       !: Sensible heat flux over the ice (diag only) 
     252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fleic       !: Latent heat flux over the ice (diag only) 
     253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qlatic      !: latent flux 
     254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdvosif     !: Variation of volume at surface (diag only) 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdvobif     !: Variation of ice volume at the bottom ice (diag only) 
     256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fdvolif     !: Total variation of ice volume (diag only) 
     257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdvonif     !: Lateral Variation of ice volume (diag only) 
     258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin] 
     259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icethi      !: total ice thickness (for all categories) (diag only) 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
     261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hicifp      !: Ice production/melting==>!obsolete... can be removed 
     262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction 
     263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pfrld       !: Leads fraction at previous time   
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   phicif      !: Old ice thickness 
     265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif        !: Heat flux at the ice base 
     266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif     !: Variation of snow mass 
     267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif     !: Variation of ice mass 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif       !: heat balance of the lead (or of the open ocean) 
     269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif       !: Energy needed to bring the ocean to freezing  
     270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fdtcn       !: net downward heat flux from the ice to the ocean 
     271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qdtcn       !: energy from the ice to the ocean 
     272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fstric      !: transmitted solar radiation under ice 
     273   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fscmbq      !: associated with lead chipotage with solar flux 
     274   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ffltbif     !: related to max heat contained in brine pockets (?) 
     275   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbbq       !: Also linked with the solar flux below the ice (?) 
     276   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qfvbq       !: store energy in case of total lateral ablation (?) 
     277   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dmgwi       !: Variation of the mass of snow ice 
     278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness 
     279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection 
     280   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation 
     281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation 
     282   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: Mass flux due to snow loss during compression 
     284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: Equivalent salt flux due to ice growth/melt 
     285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhmec       !: Heat flux due to snow loss during compression 
     286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: Residual heat flux due to correction of ice thickness 
    310287 
    311288   ! temporary arrays for dummy version of the code 
    312    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    313       dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 
    314  
    315    INTEGER, PUBLIC, DIMENSION(jpi, jpj, jpl) ::          &   !:: 
    316       patho_case ! number of the pathological case (if any, of course) 
     289   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 
     290 
     291   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   patho_case ! number of the pathological case (if any, of course) 
    317292 
    318293   !!-------------------------------------------------------------------------- 
     
    320295   !!-------------------------------------------------------------------------- 
    321296   !! Variables defined for each ice category 
    322    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl)        ::   &  !: 
    323       ht_i   ,   &  !: Ice thickness (m) 
    324       a_i    ,   &  !: Ice fractional areas (concentration) 
    325       v_i    ,   &  !: Ice volume per unit area (m) 
    326       v_s    ,   &  !: Snow volume per unit area(m) 
    327       ht_s   ,   &  !: Snow thickness (m) 
    328       t_su   ,   &  !: Sea-Ice Surface Temperature (K) 
    329       sm_i   ,   &  !: Sea-Ice Bulk salinity (ppt) 
    330       smv_i  ,   &  !: Sea-Ice Bulk salinity times volume per area (ppt.m) 
    331                                 !: this is an extensive variable that has to be transported 
    332       o_i    ,   &  !: Sea-Ice Age (days) 
    333       ov_i   ,   &  !: Sea-Ice Age times volume per area (days.m) 
    334       oa_i          !: Sea-Ice Age times ice area (days) 
    335  
    336    !! Variables summed over all categories, or associated to  
    337    !! all the ice in a single grid cell 
    338    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    339       u_ice, v_ice,   &  !: two components of the ice velocity (m/s) 
    340       tio_u, tio_v,   &  !: two components of the ice-ocean stress (N/m2) 
    341       vt_i        ,   &  !: ice total volume per unit area (m) 
    342       vt_s        ,   &  !: snow total volume per unit area (m) 
    343       at_i        ,   &  !: ice total fractional area (ice concentration) 
    344       ato_i       ,   &  !: total open water fractional area (1-at_i) 
    345       et_i        ,   &  !: total ice heat content 
    346       et_s        ,   &  !: total snow heat content 
    347       ot_i        ,   &  !: mean age over all categories 
    348       tm_i        ,   &  !: mean ice temperature over all categories 
    349       bv_i        ,   &  !: brine volume averaged over all categories 
    350       smt_i              !: mean sea ice salinity averaged over all categories 
    351  
    352    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpm) ::   &  !: 
    353       at_i_typ    ,   &  !: total area contained in each ice type 
    354       vt_i_typ           !: total volume contained in each ice type 
    355  
    356    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,nlay_s,jpl) :: & !: 
    357       t_s,            &  !: Snow temperatures (K) 
    358       e_s 
    359  
    360    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: ! go to trash 
    361       e_i_cat 
    362  
    363    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) :: & !: 
    364       t_i,            &  !: Ice temperatures     [ Kelvins     ] 
    365       e_i,            &  !: Ice thermal contents [ Joules*10^9 ] 
    366       s_i                !: Ice salinities 
     297   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ht_i    !: Ice thickness (m) 
     298   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i     !: Ice fractional areas (concentration) 
     299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_i     !: Ice volume per unit area (m) 
     300   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_s     !: Snow volume per unit area(m) 
     301   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ht_s    !: Snow thickness (m) 
     302   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_su    !: Sea-Ice Surface Temperature (K) 
     303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sm_i    !: Sea-Ice Bulk salinity (ppt) 
     304   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   smv_i   !: Sea-Ice Bulk salinity times volume per area (ppt.m) 
     305   !                                                                  !  this is an extensive variable that has to be transported 
     306   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days) 
     307   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i    !: Sea-Ice Age times volume per area (days.m) 
     308   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days) 
     309 
     310   !! Variables summed over all categories, or associated to all the ice in a single grid cell 
     311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice   !: two components of the ice velocity (m/s) 
     312   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v   !: two components of the ice-ocean stress (N/m2) 
     313   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i           !: ice total volume per unit area (m) 
     314   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_s           !: snow total volume per unit area (m) 
     315   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i           !: ice total fractional area (ice concentration) 
     316   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i          !: total open water fractional area (1-at_i) 
     317   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i           !: total ice heat content 
     318   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_s           !: total snow heat content 
     319   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ot_i           !: mean age over all categories 
     320   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i           !: mean ice temperature over all categories 
     321   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bv_i           !: brine volume averaged over all categories 
     322   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories 
     323 
     324   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   at_i_typ   !: total area contained in each ice type 
     325   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   vt_i_typ   !: total volume contained in each ice type 
     326 
     327   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s    !: Snow temperatures (K) 
     328   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s    !: Snow ...       
     329 
     330   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_i_cat   !: ! go to trash 
     331       
     332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i   !: Ice temperatures     [ Kelvins     ] 
     333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i   !: Ice thermal contents [ Joules*10^9 ] 
     334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   s_i   !: Ice salinities 
    367335 
    368336   !!-------------------------------------------------------------------------- 
    369337   !! * Moments for advection 
    370338   !!-------------------------------------------------------------------------- 
    371    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  !: 
    372       sxopw, syopw, sxxopw, syyopw, sxyopw          !: open water in sea ice 
    373  
    374    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: 
    375       sxice, syice, sxxice, syyice, sxyice,      &  !: ice thickness moments for advection 
    376       sxsn,  sysn,  sxxsn,  syysn,  sxysn,       &  !: snow thickness 
    377       sxa,   sya,   sxxa,   syya,   sxya,        &  !: lead fraction 
    378       sxc0,  syc0,  sxxc0,  syyc0,  sxyc0,       &  !: snow thermal content 
    379       sxsal, sysal, sxxsal, syysal, sxysal,      &  !: ice salinity 
    380       sxage, syage, sxxage, syyage, sxyage          !: ice age 
    381  
    382    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) ::   &  !: 
    383       sxe ,  sye ,  sxxe ,  syye ,  sxye            !: ice layers heat content 
     339   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   !: open water in sea ice 
     340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   !: ice thickness  
     341   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: snow thickness 
     342   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa,   sya,   sxxa,   syya,   sxya     !: lead fraction 
     343   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxc0,  syc0,  sxxc0,  syyc0,  sxyc0    !: snow thermal content 
     344   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   !: ice salinity 
     345   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   !: ice age 
     346   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe ,  sye ,  sxxe ,  syye ,  sxye     !: ice layers heat content 
    384347 
    385348   !!-------------------------------------------------------------------------- 
    386349   !! * Old values of global variables 
    387350   !!-------------------------------------------------------------------------- 
    388    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !:  
    389       old_v_s, old_v_i,                          &  !: snow and ice volumes 
    390       old_a_i, old_smv_i, old_oa_i 
    391    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,nlay_s,jpl) ::   &  !: 
    392       old_e_s                                       !: snow heat content 
    393    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) ::   &  !: 
    394       old_e_i                          !: ice temperatures 
    395    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: ice velocity (gv6 and gv7) 
    396       old_u_ice, old_v_ice 
     351   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   old_v_s, old_v_i               !: snow and ice volumes 
     352   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   old_a_i, old_smv_i, old_oa_i   !: ??? 
     353   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   old_e_s                        !: snow heat content 
     354   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   old_e_i                        !: ice temperatures 
     355   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   old_u_ice, old_v_ice           !: ice velocity (gv6 and gv7) 
     356       
    397357 
    398358   !!-------------------------------------------------------------------------- 
     
    401361   ! thd refers to changes induced by thermodynamics 
    402362   ! trp   ''         ''     ''       advection (transport of ice) 
    403    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !:  
    404       d_a_i_thd  , d_a_i_trp ,                   &  !: icefractions                   
    405       d_v_s_thd  , d_v_s_trp,                    &  !: snow volume 
    406       d_v_i_thd  , d_v_i_trp,                    &  !: ice  volume 
    407       d_smv_i_thd, d_smv_i_trp,                  &      
    408       d_sm_i_fl  , d_sm_i_gd ,                   & 
    409       d_sm_i_se  , d_sm_i_si , d_sm_i_la ,       & 
    410       d_oa_i_thd , d_oa_i_trp, s_i_newice 
    411  
    412    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,nlay_s,jpl) ::  &  !: 
    413       d_e_s_thd, d_e_s_trp 
    414  
    415    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) ::   &  !: 
    416       d_e_i_thd, d_e_i_trp 
    417  
    418    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::       &  !: ice velocity  
    419       d_u_ice_dyn, d_v_ice_dyn 
    420  
     363   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_a_i_thd  , d_a_i_trp                 !: icefractions                   
     364   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_s_thd  , d_v_s_trp                 !: snow volume 
     365   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_i_thd  , d_v_i_trp                 !: ice  volume 
     366   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_smv_i_thd, d_smv_i_trp               !:      
     367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_fl  , d_sm_i_gd                 !: 
     368   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_se  , d_sm_i_si  , d_sm_i_la    !: 
     369   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_oa_i_thd , d_oa_i_trp , s_i_newice   !: 
     370 
     371   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_s_thd  , d_e_s_trp     !: 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_i_thd  , d_e_i_trp     !: 
     373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   d_u_ice_dyn, d_v_ice_dyn   !: ice velocity  
     374       
    421375   !!-------------------------------------------------------------------------- 
    422376   !! * Ice thickness distribution variables 
    423377   !!-------------------------------------------------------------------------- 
    424378   ! REMOVE 
    425    INTEGER, PUBLIC, DIMENSION(jpl)                ::   &  !: 
    426       ice_types      !: Vector making the connection between types and categories 
    427  
    428    INTEGER, PUBLIC, DIMENSION(jpm,2)              ::   &  !: 
    429       ice_cat_bounds !: Matrix containing the integer upper and  
    430    !: lower boundaries of ice thickness categories 
    431  
     379   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)    ::   ice_types      !: Vector connecting types and categories 
     380   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   ice_cat_bounds !: Matrix containing the integer upper and  
     381   !                                                                       !  lower boundaries of ice thickness categories 
    432382   ! REMOVE 
    433    INTEGER, PUBLIC, DIMENSION(jpm)                ::   &  !: 
    434       ice_ncat_types !: Vector containing the number of thickness categories in each ice type 
    435  
    436    REAL(wp), PUBLIC, DIMENSION(0:jpl)             ::   &  !: 
    437       hi_max          !: Boundary of ice thickness categories in thickness space 
    438  
    439    REAL(wp), PUBLIC, DIMENSION(jpl)               ::   &  !: 
    440       hi_mean         !: Mean ice thickness in catgories  
    441  
     383   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)    ::   ice_ncat_types !: nb of thickness categories in each ice type 
     384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max         !: Boundary of ice thickness categories in thickness space 
     385   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean        !: Mean ice thickness in catgories  
    442386   ! REMOVE 
    443    REAL(wp), PUBLIC, DIMENSION(0:jpl,jpm)         ::   &  !: 
    444       hi_max_typ     !: Boundary of ice thickness categories  
    445    !:in thickness space (same but specific for each ice type) 
     387   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_max_typ     !: Boundary of ice thickness categories  
     388   !                                                                       !  in thickness space (same but specific for each ice type) 
    446389 
    447390   !!-------------------------------------------------------------------------- 
    448391   !! * Ice Run 
    449392   !!-------------------------------------------------------------------------- 
    450    !! Namelist namicerun read in iceini 
     393   !                                                                     !!: ** Namelist namicerun read in iceini ** 
    451394   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
    452395   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
     
    463406   !!-------------------------------------------------------------------------- 
    464407   !! Check if everything down here is necessary 
    465    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: volume of ice formed in the leads 
    466       v_newice 
    467    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: thermodynamic growth rates  
    468       dv_dt_thd,  & 
    469       izero, fstroc, fhbricat 
    470    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  
    471       diag_sni_gr,                           & ! snow ice growth  
    472       diag_lat_gr,                           & ! lateral ice growth  
    473       diag_bot_gr,                           & ! bottom ice growth  
    474       diag_dyn_gr,                           & ! dynamical ice growth  
    475       diag_bot_me,                           & ! vertical bottom melt  
    476       diag_sur_me                              ! vertical surface melt 
     408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   v_newice   !: volume of ice formed in the leads 
     409   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd  !: thermodynamic growth rates  
     410   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   izero, fstroc, fhbricat 
     411   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_sni_gr   ! snow ice growth  
     412   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_lat_gr   ! lateral ice growth  
     413   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_bot_gr   ! bottom ice growth  
     414   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_dyn_gr   ! dynamical ice growth  
     415   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_bot_me   ! vertical bottom melt  
     416   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_sur_me   ! vertical surface melt 
    477417   INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point 
     418 
     419   !!---------------------------------------------------------------------- 
     420   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     421   !! $Id$ 
     422   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     423   !!---------------------------------------------------------------------- 
     424CONTAINS 
     425 
     426   FUNCTION ice_alloc() 
     427      !!----------------------------------------------------------------- 
     428      !!               *** Routine ice_alloc_2 *** 
     429      !!----------------------------------------------------------------- 
     430      INTEGER :: ice_alloc 
     431      ! 
     432      INTEGER :: ierr(20), ii 
     433      !!----------------------------------------------------------------- 
     434 
     435      ierr(:) = 0 
     436 
     437      ! What could be one huge allocate statement is broken-up to try to 
     438      ! stay within Fortran's max-line length limit. 
     439      ii = 1 
     440      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) ,                           & 
     441         &      ahiu     (jpi,jpj) , ahiv     (jpi,jpj) ,                           & 
     442         &      pahu     (jpi,jpj) , pahv     (jpi,jpj) ,                           & 
     443         &      ust2s    (jpi,jpj) , hicol    (jpi,jpj) ,                           & 
     444         &      strp1    (jpi,jpj) , strp2    (jpi,jpj) , strength  (jpi,jpj) ,     & 
     445         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,     & 
     446         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) ) 
     447 
     448      ii = ii + 1 
     449      ALLOCATE( firic    (jpi,jpj) , fcsic  (jpi,jpj) , fleic    (jpi,jpj) , qlatic   (jpi,jpj) ,     & 
     450         &      rdvosif  (jpi,jpj) , rdvobif(jpi,jpj) , fdvolif  (jpi,jpj) , rdvonif  (jpi,jpj) ,     & 
     451         &      sist     (jpi,jpj) , icethi (jpi,jpj) , t_bo     (jpi,jpj) , hicifp   (jpi,jpj) ,     & 
     452         &      frld     (jpi,jpj) , pfrld  (jpi,jpj) , phicif   (jpi,jpj) , fbif     (jpi,jpj) ,     & 
     453         &      rdmsnif  (jpi,jpj) , rdmicif(jpi,jpj) , qldif    (jpi,jpj) , qcmif    (jpi,jpj) ,     & 
     454         &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric   (jpi,jpj) , fscmbq   (jpi,jpj) ,     & 
     455         &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq    (jpi,jpj) , dmgwi    (jpi,jpj) ,     & 
     456         &      fsalt_res(jpi,jpj) , fsbri  (jpi,jpj) , fsalt_rpo(jpi,jpj) , fheat_rpo(jpi,jpj) ,     & 
     457         &      fhbri    (jpi,jpj) , fmmec  (jpi,jpj) , fseqv    (jpi,jpj) , fhmec    (jpi,jpj) ,     & 
     458         &      fheat_res(jpi,jpj)                                                              , STAT=ierr(ii) ) 
     459 
     460      ii = ii + 1 
     461      ALLOCATE( dh_i_surf2D(jpi,jpj) , dh_i_bott2D(jpi,jpj) , fstbif(jpi,jpj) ,     & 
     462         &      fsup2D     (jpi,jpj) , focea2D    (jpi,jpj) , q_s   (jpi,jpj) , STAT=ierr(ii) ) 
     463 
     464      ii = ii + 1 
     465      ALLOCATE( patho_case(jpi, jpj, jpl)  , STAT=ierr(ii) ) 
     466 
     467      ! * Ice global state variables 
     468      ii = ii + 1 
     469      ALLOCATE( ht_i (jpi,jpj,jpl) , a_i  (jpi,jpj,jpl) , v_i  (jpi,jpj,jpl) ,     & 
     470         &      v_s  (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) ,     & 
     471         &      sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i  (jpi,jpj,jpl) ,     & 
     472         &      ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl)                      , STAT=ierr(ii) ) 
     473      ii = ii + 1 
     474      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     & 
     475         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     & 
     476         &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     & 
     477         &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) ) 
     478      ii = ii + 1 
     479      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , at_i_typ(jpi,jpj,jpm) ,                            & 
     480         &      e_s(jpi,jpj,nlay_s,jpl) , vt_i_typ(jpi,jpj,jpm) , e_i_cat(jpi,jpj,jpl) , STAT=ierr(ii) ) 
     481      ii = ii + 1 
     482      ALLOCATE( t_i(jpi,jpj,jkmax,jpl) , e_i(jpi,jpj,jkmax,jpl) , s_i(jpi,jpj,jkmax,jpl) , STAT=ierr(ii) ) 
     483 
     484      ! * Moments for advection 
     485      ii = ii + 1 
     486      ALLOCATE( sxopw(jpi,jpj) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj) , STAT=ierr(ii) ) 
     487      ii = ii + 1 
     488      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
     489         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   & 
     490         &      STAT=ierr(ii) ) 
     491      ii = ii + 1 
     492      ALLOCATE( sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   & 
     493         &      sxc0 (jpi,jpj,jpl) , syc0 (jpi,jpj,jpl) , sxxc0 (jpi,jpj,jpl) , syyc0 (jpi,jpj,jpl) , sxyc0 (jpi,jpj,jpl) ,   & 
     494         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   & 
     495         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   & 
     496         &      STAT=ierr(ii) ) 
     497      ii = ii + 1 
     498      ALLOCATE( sxe (jpi,jpj,jkmax,jpl) , sye (jpi,jpj,jkmax,jpl) , sxxe(jpi,jpj,jkmax,jpl) ,     & 
     499         &      syye(jpi,jpj,jkmax,jpl) , sxye(jpi,jpj,jkmax,jpl)                           , STAT=ierr(ii) ) 
     500 
     501      ! * Old values of global variables 
     502      ii = ii + 1 
     503      ALLOCATE( old_v_s  (jpi,jpj,jpl) , old_v_i  (jpi,jpj,jpl) , old_e_s(jpi,jpj,nlay_s,jpl) ,     & 
     504         &      old_a_i  (jpi,jpj,jpl) , old_smv_i(jpi,jpj,jpl) , old_e_i(jpi,jpj,jkmax ,jpl) ,     & 
     505         &      old_oa_i (jpi,jpj,jpl)                                                        ,     & 
     506         &      old_u_ice(jpi,jpj)     , old_v_ice(jpi,jpj)                                   , STAT=ierr(ii) ) 
     507 
     508      ! * Increment of global variables 
     509      ii = ii + 1 
     510      ALLOCATE( d_a_i_thd(jpi,jpj,jpl) , d_a_i_trp (jpi,jpj,jpl) , d_v_s_thd  (jpi,jpj,jpl) , d_v_s_trp  (jpi,jpj,jpl) ,   & 
     511         &      d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) ,   &      
     512         &      d_sm_i_fl(jpi,jpj,jpl) , d_sm_i_gd (jpi,jpj,jpl) , d_sm_i_se  (jpi,jpj,jpl) , d_sm_i_si  (jpi,jpj,jpl) ,   & 
     513         &      d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) , s_i_newice (jpi,jpj,jpl) ,   & 
     514         &     STAT=ierr(ii) ) 
     515      ii = ii + 1 
     516      ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,jkmax,jpl) , d_u_ice_dyn(jpi,jpj) ,     & 
     517         &      d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,jkmax,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) ) 
     518       
     519      ! * Ice thickness distribution variables 
     520      ii = ii + 1 
     521      ALLOCATE( ice_types(jpl) , ice_cat_bounds(jpm,2) , ice_ncat_types  (jpm) ,     & 
     522         &      hi_max (0:jpl) , hi_mean(jpl)          , hi_max_typ(0:jpl,jpm) , STAT=ierr(ii) ) 
     523 
     524      ! * Ice diagnostics 
     525      ii = ii + 1 
     526      ALLOCATE( dv_dt_thd(jpi,jpj,jpl) , diag_sni_gr(jpi,jpj) , diag_lat_gr(jpi,jpj) ,     & 
     527         &      izero    (jpi,jpj,jpl) , diag_bot_gr(jpi,jpj) , diag_dyn_gr(jpi,jpj) ,     & 
     528         &      fstroc   (jpi,jpj,jpl) , diag_bot_me(jpi,jpj) , diag_sur_me(jpi,jpj) ,     & 
     529         &      fhbricat (jpi,jpj,jpl) , v_newice   (jpi,jpj)                        , STAT=ierr(ii) ) 
     530 
     531      ice_alloc = MAXVAL( ierr(:) ) 
     532      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays.') 
     533      ! 
     534   END FUNCTION ice_alloc 
    478535 
    479536#else 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r2528 r2601  
    1717   USE sbc_oce        ! Surface boundary condition: ocean fields 
    1818   USE sbc_ice        ! Surface boundary condition: ice fields 
     19   USE par_ice        ! LIM: sea-ice parameters 
    1920   USE ice            ! LIM: sea-ice variables 
    2021   USE limmsh         ! LIM: mesh 
     
    2324   USE limthd         ! LIM: ice thermodynamics 
    2425   USE limthd_sal     ! LIM: ice thermodynamics: salinity 
    25    USE par_ice        ! LIM: sea-ice parameters 
    2626   USE limvar         ! LIM: variables 
    2727   USE in_out_manager ! I/O manager 
     
    3131   PRIVATE 
    3232 
    33    PUBLIC   ice_init   ! called by opa.F90 
    34  
    35    !!---------------------------------------------------------------------- 
    36    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     33   PUBLIC   ice_init   ! called by sbcice_lim.F90 
     34 
     35   !!---------------------------------------------------------------------- 
     36   !! NEMO/LIM3 4/0 , UCL - NEMO Consortium (2010) 
    3737   !! $Id$ 
    3838   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    4646      !! ** purpose :    
    4747      !!---------------------------------------------------------------------- 
     48      INTEGER :: ierr 
     49      !!---------------------------------------------------------------------- 
     50 
     51      !                                ! Allocate the ice arrays 
     52      ierr = ice_alloc()                     ! NB: Calls to the _alloc() routines should be in  
     53      !                                      !     the same order as the ice modules are USE'd above 
     54       
     55!     ierr = ierr + ice_alloc_2() 
     56!     ierr = ierr + lim_dia_alloc_2() 
     57!     ierr = ierr + lim_hdf_alloc_2() 
     58!     ierr = ierr + lim_sbc_alloc_2() 
     59!     ierr = ierr + lim_wri_alloc_2() 
     60!     ierr = ierr + thd_ice_alloc_2() 
     61 
     62!     ierr = ierr + lim_rhg_alloc() 
     63!     ierr = ierr + dom_ice_alloc() 
     64!     ierr = ierr + lim_idt_me_alloc() 
     65!     ierr = ierr + thd_ice_alloc() 
     66 
     67      IF( lk_mpp )   CALL mpp_sum( ierr ) 
     68 
     69      IF( ierr > 0 ) THEN 
     70         WRITE(numout,*)  
     71         WRITE(numout,*) 'ERROR: Allocation of memory failed in nemo_alloc' 
     72         IF( lk_mpp ) CALL mppstop() 
     73         STOP 
     74      END IF 
    4875      ! 
    4976      !                                ! Open the namelist file  
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90

    r2528 r2601  
    5252 
    5353   REAL(wp), DIMENSION(jpinfmx) ::   vinfom     ! temporary working space 
    54    REAL(wp), DIMENSION(jpi,jpj) ::   aire       ! masked grid cell area 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aire       ! masked grid cell area 
    5555 
    5656   !! * Substitutions 
     
    6767      !!                  ***  ROUTINE lim_dia  *** 
    6868      !!    
    69       !! ** Purpose : Computation and outputs on file ice.evolu  
    70       !!      the temporal evolution of some key variables 
     69      !! ** Purpose :   Computation and outputs on file ice.evolu  
     70      !!              the temporal evolution of some key variables 
    7171      !!------------------------------------------------------------------- 
    7272      INTEGER  ::   jv, ji, jj, jl   ! dummy loop indices 
     
    410410      !!------------------------------------------------------------------- 
    411411      INTEGER  ::   jv    ! dummy loop indice 
    412       INTEGER  ::   ntot , ndeb , irecl   ! local integers 
     412      INTEGER  ::   ierr, ntot , ndeb , irecl   ! local integers 
    413413      REAL(wp) ::   zxx0, zxx1    ! local scalars 
    414414      CHARACTER(len=jpchinf) ::   titinf 
     
    431431      ENDIF 
    432432 
     433      ALLOCATE( aire(jpi,jpj) , STAT=ierr ) 
     434      IF( ierr /= 0 ) THEN 
     435         CALL ctl_stop( 'lim_dia_init : unable to allocate standard arrays' )   ;   RETURN 
     436      ENDIF 
     437       
    433438      aire(:,:) = area(:,:) * tms(:,:) * tmask_i(:,:)      ! masked grid cell area (interior domain only) 
    434439 
     
    436441      titvar(1) = 'NoIt'  ! iteration number 
    437442      titvar(2) = 'T yr'  ! time step in years 
    438       nbvt = 2            ! number of time variables 
     443      nbvt      = 2       ! number of time variables 
    439444 
    440445      titvar(3) = 'AI_N'  ! sea ice area in the northern Hemisp.(10^12 km2) 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r2528 r2601  
    44   !! LIM ice model : horizontal diffusion of sea-ice quantities 
    55   !!====================================================================== 
     6   !! History :  LIM  !  2000-01 (LIM) Original code 
     7   !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm 
     8   !!            1.0  !  2002-08 (C. Ethe)  F90, free form 
     9   !!---------------------------------------------------------------------- 
    610#if defined key_lim3 
    711   !!---------------------------------------------------------------------- 
     
    1014   !!   lim_hdf  : diffusion trend on sea-ice variable 
    1115   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    13    USE dom_oce 
    14    USE in_out_manager 
    15    USE ice 
    16    USE lbclnk 
    17    USE lib_mpp 
    18    USE prtctl          ! Print control 
     16   USE dom_oce          ! ocean domain 
     17   USE ice              ! LIM-3: ice variables 
     18   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     19   USE lib_mpp          ! MPP library 
     20   USE prtctl           ! Print control 
     21   USE in_out_manager   ! I/O manager 
    1922 
    2023   IMPLICIT NONE 
    2124   PRIVATE 
    2225 
    23    !! * Routine accessibility 
    24    PUBLIC lim_hdf    ! called by lim_tra 
     26   PUBLIC   lim_hdf     ! called by lim_tra 
    2527 
    26    !! * Module variables 
    27    LOGICAL  ::   linit = .TRUE.              ! ??? 
     28   LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call) 
    2829   REAL(wp) ::   epsi04 = 1e-04              ! constant 
    29    REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! ??? 
     30   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient 
    3031 
    3132   !! * Substitution  
    3233#  include "vectopt_loop_substitute.h90" 
    3334   !!---------------------------------------------------------------------- 
    34    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     35   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    3536   !! $Id$ 
    36    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3738   !!---------------------------------------------------------------------- 
    38  
    3939CONTAINS 
    4040 
     
    4343      !!                  ***  ROUTINE lim_hdf  *** 
    4444      !! 
    45       !! ** purpose :   Compute and add the diffusive trend on sea-ice 
    46       !!      variables 
     45      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables 
    4746      !! 
    4847      !! ** method  :   Second order diffusive operator evaluated using a 
    49       !!      Cranck-Nicholson time Scheme. 
     48      !!              Cranck-Nicholson time Scheme. 
    5049      !! 
    5150      !! ** Action  :    update ptab with the diffusive contribution 
    52       !! 
    53       !! History : 
    54       !!        !  00-01 (LIM) Original code 
    55       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    56       !!        !  02-08 (C. Ethe)  F90, free form 
    5751      !!------------------------------------------------------------------- 
    58       ! * Arguments 
    59       REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   & 
    60          ptab                 ! Field on which the diffusion is applied   
    61       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    62          ptab0                ! ??? 
     52      USE wrk_nemo, ONLY:   wrk_use, wrk_release 
     53      USE wrk_nemo, ONLY:   zflu => wrk_2d_11, zdiv  => wrk_2d_13, zrlx  => wrk_2d_15  
     54      USE wrk_nemo, ONLY:   zflv => wrk_2d_12, zdiv0 => wrk_2d_14, ztab0 => wrk_2d_16 
     55      ! 
     56      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied 
     57      ! 
     58      INTEGER  ::  ji, jj            ! dummy loop indices 
     59      INTEGER  ::  its, iter, ierr   ! local integers 
     60      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars 
     61      CHARACTER(lc) ::   charout   ! local character 
     62      !!------------------------------------------------------------------- 
     63       
     64      IF( .NOT. wrk_use(2, 11,12,13,14,15,16) ) THEN 
     65         CALL ctl_stop( 'lim_hdf: ERROR: requested workspace arrays unavailable' )   ;   RETURN 
     66      END IF 
    6367 
    64       ! * Local variables 
    65       INTEGER ::  ji, jj      ! dummy loop indices 
    66       INTEGER ::  & 
    67          its, iter            ! temporary integers 
    68       CHARACTER (len=55) :: charout 
    69       REAL(wp) ::  & 
    70          zalfa, zrlxint, zconv, zeps   ! temporary scalars 
    71       REAL(wp), DIMENSION(jpi,jpj) ::  &  
    72          zrlx, zflu, zflv, &  ! temporary workspaces 
    73          zdiv0, zdiv          !    "           " 
    74       !!------------------------------------------------------------------- 
    75  
    76       ! Initialisation 
    77       ! ---------------    
    78       ! Time integration parameters 
    79       zalfa = 0.5       ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
    80       its   = 100       ! Maximum number of iteration 
    81       zeps  =  2. * epsi04 
    82  
    83       ! Arrays initialization 
    84       ptab0 (:, : ) = ptab(:,:) 
    85       !bug  zflu (:,jpj) = 0.e0 
    86       !bug  zflv (:,jpj) = 0.e0 
    87       zdiv0(:, 1 ) = 0.e0 
    88       zdiv0(:,jpj) = 0.e0 
    89       IF( .NOT.lk_vopt_loop ) THEN 
    90          zflu (jpi,:) = 0.e0    
    91          zflv (jpi,:) = 0.e0 
    92          zdiv0(1,  :) = 0.e0 
    93          zdiv0(jpi,:) = 0.e0 
    94       ENDIF 
    95  
    96       ! Metric coefficient (compute at the first call and saved in 
    97       IF( linit ) THEN 
     68      !                       !==  Initialisation  ==! 
     69      ! 
     70      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact) 
     71         ALLOCATE( efact(jpi,jpj) , STAT=ierr ) 
     72         IF( ierr /= 0 ) THEN 
     73            CALL ctl_stop( 'lim_hdf : unable to allocate standard arrays' )   ;   RETURN 
     74         ENDIF 
    9875         DO jj = 2, jpjm1   
    9976            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    100                zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) & 
     77               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) & 
    10178                  &          / ( e1t(ji,jj) * e2t(ji,jj) ) 
    10279            END DO 
     
    10481         linit = .FALSE. 
    10582      ENDIF 
     83      !                             ! Time integration parameters 
     84      zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
     85      its   = 100                         ! Maximum number of iteration 
     86      zeps  =  2._wp * epsi04 
     87      ! 
     88      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization 
     89      zdiv0(:, 1 ) = 0._wp 
     90      zdiv0(:,jpj) = 0._wp 
     91      IF( .NOT.lk_vopt_loop ) THEN 
     92         zflu (jpi,:) = 0._wp    
     93         zflv (jpi,:) = 0._wp 
     94         zdiv0(1,  :) = 0._wp 
     95         zdiv0(jpi,:) = 0._wp 
     96      ENDIF 
    10697 
    107  
    108       ! Sub-time step loop 
    109       zconv = 1.e0 
     98      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==! 
    11099      iter  = 0 
    111  
    112       !                                                   !=================== 
    113       DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   ! Sub-time step loop 
    114          !                                                !=================== 
    115          ! incrementation of the sub-time step number 
    116          iter = iter + 1 
    117  
    118          ! diffusive fluxes in U- and V- direction 
    119          DO jj = 1, jpjm1 
     100      ! 
     101      DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop 
     102         ! 
     103         iter = iter + 1                                 ! incrementation of the sub-time step number 
     104         ! 
     105         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
    120106            DO ji = 1 , fs_jpim1   ! vector opt. 
    121107               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     
    123109            END DO 
    124110         END DO 
    125  
    126          ! diffusive trend : divergence of the fluxes 
    127          DO jj= 2, jpjm1 
     111         ! 
     112         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
    128113            DO ji = fs_2 , fs_jpim1   ! vector opt.  
    129114               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
     
    131116            END DO 
    132117         END DO 
    133  
    134          ! save the first evaluation of the diffusive trend in zdiv0 
    135          IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        
    136  
    137          ! XXXX iterative evaluation????? 
    138          DO jj = 2, jpjm1 
     118         ! 
     119         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0 
     120         ! 
     121         DO jj = 2, jpjm1                                ! iterative evaluation 
    139122            DO ji = fs_2 , fs_jpim1   ! vector opt. 
    140                zrlxint = (   ptab0(ji,jj)    & 
    141                   &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) )   & 
     123               zrlxint = (   ztab0(ji,jj)    & 
     124                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   & 
    142125                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             &  
    143                   &    / ( 1.0 + zalfa * rdt_ice * zfact(ji,jj) ) 
     126                  &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) ) 
    144127               zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) ) 
    145128            END DO 
    146129         END DO 
    147  
    148          ! lateral boundary condition on zrlx 
    149          CALL lbc_lnk( zrlx, 'T', 1. ) 
    150  
    151          ! convergence test 
    152          zconv = 0.e0 
     130         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition 
     131         ! 
     132         zconv = 0._wp                                   ! convergence test 
    153133         DO jj = 2, jpjm1 
    154134            DO ji = fs_2, fs_jpim1 
     
    156136            END DO 
    157137         END DO 
    158          IF( lk_mpp )   CALL mpp_max( zconv )   ! max over the global domain 
    159  
    160          DO jj = 1, jpj 
    161             DO ji = 1 , jpi 
    162                ptab(ji,jj) = zrlx(ji,jj) 
    163             END DO 
    164          END DO 
    165  
    166          !                                         !========================== 
     138         IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain 
     139         ! 
     140         ptab(:,:) = zrlx(:,:) 
     141         ! 
    167142      END DO                                       ! end of sub-time step loop 
    168       !                                            !========================== 
    169143 
    170144      IF(ln_ctl)   THEN 
    171          zrlx(:,:) = ptab(:,:) - ptab0(:,:) 
     145         zrlx(:,:) = ptab(:,:) - ztab0(:,:) 
    172146         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter 
    173          CALL prt_ctl(tab2d_1=zrlx, clinfo1=charout) 
     147         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout ) 
    174148      ENDIF 
    175  
    176  
     149      ! 
     150      IF( .NOT. wrk_release(2, 11,12,13,14,15,16) ) THEN 
     151         CALL ctl_stop( 'lim_hdf : failed to release workspace arrays.' )   ;   RETURN 
     152      END IF 
     153      ! 
    177154   END SUBROUTINE lim_hdf 
    178155 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r2590 r2601  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limitd_me *** 
    4    !!            Mechanical impact on ice thickness distribution 
    5    !!                     computation of changes in g(h)       
     4   !! LIM-3 : Mechanical impact on ice thickness distribution       
    65   !!====================================================================== 
    76   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code  
     
    1211   !!   'key_lim3' :                                    LIM3 sea-ice model 
    1312   !!---------------------------------------------------------------------- 
    14    USE dom_ice 
    1513   USE par_oce          ! ocean parameters 
    16    USE dom_oce 
    17    USE lbclnk 
     14   USE dom_oce          ! ocean domain 
    1815   USE phycst           ! physical constants (ocean directory)  
    19    USE sbc_oce          ! Surface boundary condition: ocean fields 
    20    USE thd_ice 
    21    USE in_out_manager 
    22    USE ice 
    23    USE par_ice 
    24    USE limthd_lac 
    25    USE limvar 
    26    USE limcons 
     16   USE sbc_oce          ! surface boundary condition: ocean fields 
     17   USE thd_ice          ! LIM-3 thermodynamics 
     18   USE ice              ! LIM-3 variables 
     19   USE par_ice          ! LIM-3 parameters 
     20   USE dom_ice          ! LIM-3 domain 
     21   USE limthd_lac       ! LIM-3 
     22   USE limvar           ! LIM-3 
     23   USE limcons          ! LIM-3 
     24   USE in_out_manager   ! I/O manager 
    2725   USE prtctl           ! Print control 
    28    USE lib_mpp 
     26   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     27   USE lib_mpp          ! MPP library 
    2928   USE wrk_nemo, ONLY: wrk_use, wrk_release 
    3029 
     
    3231   PRIVATE 
    3332 
    34    !! * Routine accessibility 
    35    PUBLIC lim_itd_me        ! called by ice_stp 
    36    PUBLIC lim_itd_me_icestrength 
    37    PUBLIC lim_itd_me_ridgeprep 
    38    PUBLIC lim_itd_me_ridgeshift 
    39    PUBLIC lim_itd_me_asumr 
    40    PUBLIC lim_itd_me_init 
    41    PUBLIC lim_itd_me_zapsmall 
    42    PUBLIC lim_idt_me_alloc  ! called by nemogcm.F90 
    43  
    44    !! * Module variables 
    45    REAL(wp)  ::           &  ! constant values 
    46       epsi20 = 1e-20   ,  & 
    47       epsi13 = 1e-13   ,  & 
    48       epsi11 = 1e-11   ,  & 
    49       zzero  = 0.e0    ,  & 
    50       zone   = 1.e0 
     33   PUBLIC   lim_itd_me               ! called by ice_stp 
     34   PUBLIC   lim_itd_me_icestrength 
     35   PUBLIC   lim_itd_me_init 
     36   PUBLIC   lim_itd_me_zapsmall 
     37   PUBLIC   lim_itd_me_alloc        ! called by nemogcm.F90 
     38 
     39   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
     40   REAL(wp)  ::   epsi11 = 1.e-11_wp   ! constant values 
     41   REAL(wp)  ::   epsi10 = 1.e-10_wp   ! constant values 
    5142 
    5243   !----------------------------------------------------------------------- 
    5344   ! Variables shared among ridging subroutines 
    5445   !----------------------------------------------------------------------- 
    55    REAL(wp), ALLOCATABLE, SAVE, DIMENSION (:,:) ::    & 
    56       asum         , & ! sum of total ice and open water area 
    57       aksum            ! ratio of area removed to area ridged 
    58  
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: &      
    60       athorn           ! participation function; fraction of ridging/ 
    61    !  closing associated w/ category n 
    62  
    63    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  & 
    64       hrmin      , &   ! minimum ridge thickness 
    65       hrmax      , &   ! maximum ridge thickness 
    66       hraft      , &   ! thickness of rafted ice 
    67       krdg       , &   ! mean ridge thickness/thickness of ridging ice  
    68       aridge     , &   ! participating ice ridging 
    69       araft            ! participating ice rafting 
    70  
    71    REAL(wp), PARAMETER :: & 
    72       krdgmin = 1.1, &    ! min ridge thickness multiplier 
    73       kraft   = 2.0       ! rafting multipliyer 
    74  
    75    REAL(wp) :: &                                
    76       Cp  
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
     48 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
     50   !                                                           !  closing associated w/ category n 
     51 
     52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
     53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
     55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
     56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
     57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
     58 
     59   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
     60   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
     61 
     62   REAL(wp) ::   Cp  
    7763   ! 
    7864   !----------------------------------------------------------------------- 
    7965   ! Ridging diagnostic arrays for history files 
    8066   !----------------------------------------------------------------------- 
    81    ! 
    82    REAL (wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: & 
    83       dardg1dt     , & ! rate of fractional area loss by ridging ice (1/s) 
    84       dardg2dt     , & ! rate of fractional area gain by new ridges (1/s) 
    85       dvirdgdt     , & ! rate of ice volume ridged (m/s) 
    86       opening          ! rate of opening due to divergence/shear (1/s) 
    87  
     67   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg1dt   ! rate of fractional area loss by ridging ice (1/s) 
     68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg2dt   ! rate of fractional area gain by new ridges (1/s) 
     69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s) 
     70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s) 
    8871 
    8972   !!---------------------------------------------------------------------- 
    9073   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    9174   !! $Id$ 
    92    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     75   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    9376   !!---------------------------------------------------------------------- 
    94  
    95  
    9677CONTAINS 
    9778 
    98    !!-----------------------------------------------------------------------------! 
    99    !!-----------------------------------------------------------------------------! 
    100  
    101    FUNCTION lim_idt_me_alloc() 
     79   FUNCTION lim_itd_me_alloc() 
    10280      !!---------------------------------------------------------------------! 
    10381      !!                ***  ROUTINE lim_itd_me_alloc *** 
    10482      !!---------------------------------------------------------------------! 
    105       INTEGER :: lim_idt_me_alloc 
     83      INTEGER :: lim_itd_me_alloc 
    10684      !!---------------------------------------------------------------------! 
    107  
    108       ALLOCATE(asum(jpi,jpj), aksum(jpi,jpj), athorn(jpi,jpj,0:jpl), & 
    109                ! 
    110                hrmin(jpi,jpj,jpl),  hrmax(jpi,jpj,jpl)      , & 
    111                hraft(jpi,jpj,jpl),  krdg(jpi,jpj,jpl)       , & 
    112                aridge(jpi,jpj,jpl), araft(jpi,jpj,jpl)      , & 
    113                ! 
    114                dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)       , &  
    115                dvirdgdt(jpi,jpj)  , opening(jpi,jpj)        , & 
    116                !  
    117                Stat=lim_idt_me_alloc) 
    118  
    119       IF(lim_idt_me_alloc /= 0)THEN 
    120          CALL ctl_warn('lim_idt_me_alloc: failed to allocate arrays.') 
     85      ! 
     86      ALLOCATE(                                                                     & 
     87         !* Variables shared among ridging subroutines 
     88         &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl)                    ,     & 
     89         &      aksum(jpi,jpj)                                                ,     & 
     90         ! 
     91         &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) ,     & 
     92         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) ,     & 
     93         ! 
     94         !* Ridging diagnostic arrays for history files 
     95         &      dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)                        ,     &  
     96         &      dvirdgdt(jpi,jpj)  , opening(jpi,jpj)                         , STAT=lim_itd_me_alloc ) 
     97      ! 
     98      IF( lim_itd_me_alloc /= 0 ) THEN 
     99         CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays.' ) 
    121100      END IF 
    122  
    123    END FUNCTION lim_idt_me_alloc 
    124  
    125  
    126    SUBROUTINE lim_itd_me ! (subroutine 1/6) 
     101      ! 
     102   END FUNCTION lim_itd_me_alloc 
     103 
     104 
     105   SUBROUTINE lim_itd_me 
    127106      !!---------------------------------------------------------------------! 
    128107      !!                ***  ROUTINE lim_itd_me *** 
    129       !! ** Purpose : 
    130       !!        This routine computes the mechanical redistribution 
    131       !!                      of ice thickness 
    132       !! 
    133       !! ** Method  : a very simple method :-) 
    134       !! 
    135       !! ** Arguments : 
    136       !!           kideb , kiut : Starting and ending points on which the  
    137       !!                         the computation is applied 
    138       !! 
    139       !! ** Inputs / Ouputs : (global commons) 
    140       !! 
    141       !! ** External :  
    142       !! 
    143       !! ** Steps : 
    144       !!  1) Thickness categories boundaries, ice / o.w. concentrations 
    145       !!     Ridge preparation 
    146       !!  2) Dynamical inputs (closing rate, divu_adv, opning) 
    147       !!  3) Ridging iteration 
    148       !!  4) Ridging diagnostics 
    149       !!  5) Heat, salt and freshwater fluxes 
    150       !!  6) Compute increments of tate variables and come back to old values 
    151       !! 
    152       !! ** References : There are a lot of references and can be difficult /  
    153       !!                 boring to read 
    154       !! 
    155       !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength 
    156       !!  in modeling the thickness distribution of Arctic sea ice, 
    157       !!  J. Geophys. Res., 100, 18,611-18,626. 
    158       !! 
    159       !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice 
    160       !!  cover, Mon. Wea. Rev., 108, 1943-1973, 1980. 
    161       !! 
    162       !! Rothrock, D. A., 1975: The energetics of the plastic deformation of 
    163       !!  pack ice by ridging, J. Geophys. Res., 80, 4514-4519. 
    164       !! 
    165       !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,  
    166       !!  1975: The thickness distribution of sea ice, J. Geophys. Res.,  
    167       !!  80, 4501-4513.  
    168       !! 
    169       !! Bitz et al., JGR 2001 
    170       !! 
    171       !! Amundrud and Melling, JGR 2005 
    172       !! 
    173       !! Babko et al., JGR 2002  
     108      !! 
     109      !! ** Purpose :   computes the mechanical redistribution of ice thickness 
     110      !! 
     111      !! ** Method  :   Steps : 
     112      !!       1) Thickness categories boundaries, ice / o.w. concentrations 
     113      !!          Ridge preparation 
     114      !!       2) Dynamical inputs (closing rate, divu_adv, opning) 
     115      !!       3) Ridging iteration 
     116      !!       4) Ridging diagnostics 
     117      !!       5) Heat, salt and freshwater fluxes 
     118      !!       6) Compute increments of tate variables and come back to old values 
     119      !! 
     120      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626. 
     121      !!                Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980. 
     122      !!                Rothrock, D. A., 1975: JGR, 80, 4514-4519. 
     123      !!                Thorndike et al., 1975, JGR, 80, 4501-4513.  
     124      !!                Bitz et al., JGR, 2001 
     125      !!                Amundrud and Melling, JGR 2005 
     126      !!                Babko et al., JGR 2002  
    174127      !! 
    175128      !!     This routine is based on CICE code and authors William H. Lipscomb, 
    176129      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged 
    177130      !!--------------------------------------------------------------------! 
    178       USE wrk_nemo, ONLY: & 
    179           closing_net   => wrk_2d_1, &  ! net rate at which area is removed    (1/s) 
    180                                         ! (ridging ice area - area of new ridges) / dt 
    181           divu_adv      => wrk_2d_2, &  ! divu as implied by transport scheme  (1/s) 
    182           opning        => wrk_2d_3, &  ! rate of opening due to divergence/shear 
    183           closing_gross => wrk_2d_4, &  ! rate at which area removed, not counting 
    184                                         ! area of new ridges 
    185           msnow_mlt     => wrk_2d_5, &  ! mass of snow added to ocean (kg m-2) 
    186           esnow_mlt     => wrk_2d_6       ! energy needed to melt snow in ocean (J m-2) 
    187       USE wrk_nemo, ONLY: vt_i_init  => wrk_2d_7, &  !  ice volume summed over  
    188                           vt_i_final => wrk_2d_8     !  categories 
    189  
    190       !! * Arguments 
    191  
    192       !! * Local variables 
    193       INTEGER ::   ji,       &   ! spatial dummy loop index 
    194          jj,       &   ! spatial dummy loop index 
    195          jk,       &   ! vertical layering dummy loop index 
    196          jl,       &   ! ice category dummy loop index 
    197          niter,    &   ! iteration counter 
    198          nitermax = 20 ! max number of ridging iterations  
    199  
    200       REAL(wp)  ::             &  ! constant values 
    201          zeps      =  1.0e-10, & 
    202          epsi10    =  1.0e-10, & 
    203          epsi06    =  1.0e-6 
    204  
    205       REAL(wp) ::            & 
    206          w1,                 &  ! temporary variable 
    207          tmpfac,             &  ! factor by which opening/closing rates are cut 
    208          dti                    ! 1 / dt 
    209  
    210       LOGICAL   ::           & 
    211          asum_error              ! flag for asum .ne. 1 
    212  
    213       INTEGER :: iterate_ridging ! if true, repeat the ridging 
    214  
    215       REAL(wp) ::  &          
    216          big = 1.0e8 
    217  
    218       CHARACTER (len = 15) :: fieldid 
    219  
    220       !!-- End of declarations 
    221       !-----------------------------------------------------------------------------! 
    222  
    223       IF(.NOT. wrk_use(2, 1,2,3,4,5,6,7,8))THEN 
    224          CALL ctl_stop(' : requested workspace arrays unavailable.') 
    225          RETURN 
     131      USE wrk_nemo, ONLY:   closing_net   => wrk_2d_1   ! net rate at which area is removed    (1/s) 
     132      !                                                 ! (ridging ice area - area of new ridges) / dt 
     133      USE wrk_nemo, ONLY:   divu_adv      => wrk_2d_2   ! divu as implied by transport scheme  (1/s) 
     134      USE wrk_nemo, ONLY:   opning        => wrk_2d_3   ! rate of opening due to divergence/shear 
     135      USE wrk_nemo, ONLY:   closing_gross => wrk_2d_4   ! rate at which area removed, not counting area of new ridges 
     136      USE wrk_nemo, ONLY:   msnow_mlt     => wrk_2d_5   ! mass of snow added to ocean (kg m-2) 
     137      USE wrk_nemo, ONLY:   esnow_mlt     => wrk_2d_6   ! energy needed to melt snow in ocean (J m-2) 
     138      USE wrk_nemo, ONLY:   vt_i_init     => wrk_2d_7   !  ice volume summed over  
     139      USE wrk_nemo, ONLY:   vt_i_final    => wrk_2d_8   !  categories 
     140      ! 
     141      INTEGER ::   ji, jj, jk, jl   ! dummy loop index 
     142      INTEGER ::   niter, nitermax = 20   ! local integer  
     143 
     144      LOGICAL  ::   asum_error              ! flag for asum .ne. 1 
     145      INTEGER  ::   iterate_ridging ! if true, repeat the ridging 
     146      REAL(wp) ::   w1, tmpfac, dti   ! local scalar 
     147      REAL(wp) ::   big = 1.0e8 
     148      CHARACTER (len = 15) ::   fieldid 
     149      !!----------------------------------------------------------------------------- 
     150 
     151      IF( .NOT. wrk_use(2, 1,2,3,4,5,6,7,8) ) THEN 
     152         CALL ctl_stop(' : requested workspace arrays unavailable.')   ;   RETURN 
    226153      END IF 
    227154 
    228       IF( numit == nstart  ) CALL lim_itd_me_init ! Initialization (first time-step only) 
     155      IF( numit == nstart  )   CALL lim_itd_me_init  ! Initialization (first time-step only) 
    229156 
    230157      IF(ln_ctl) THEN 
     
    241168      hi_max(jpl) = 999.99 
    242169 
    243       Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0    ! proport const for PE 
     170      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0    ! proport const for PE 
    244171      CALL lim_itd_me_ridgeprep ! prepare ridging 
    245172 
    246173      ! conservation check 
    247       IF ( con_i) CALL lim_column_sum (jpl,   v_i, vt_i_init) 
     174      IF ( con_i)   CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    248175 
    249176      ! Initialize arrays. 
    250177      DO jj = 1, jpj 
    251178         DO ji = 1, jpi 
    252  
    253             msnow_mlt(ji,jj) = 0.0 
    254             esnow_mlt(ji,jj) = 0.0 
    255             dardg1dt(ji,jj)  = 0.0 
    256             dardg2dt(ji,jj)  = 0.0 
    257             dvirdgdt(ji,jj)  = 0.0 
    258             opening (ji,jj)  = 0.0 
     179            msnow_mlt(ji,jj) = 0._wp 
     180            esnow_mlt(ji,jj) = 0._wp 
     181            dardg1dt(ji,jj)  = 0._wp 
     182            dardg2dt(ji,jj)  = 0._wp 
     183            dvirdgdt(ji,jj)  = 0._wp 
     184            opening (ji,jj)  = 0._wp 
    259185 
    260186            !-----------------------------------------------------------------------------! 
     
    277203            !  (thick, newly ridged ice). 
    278204 
    279             closing_net(ji,jj) = & 
    280                Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) 
     205            closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    281206 
    282207            ! 2.2 divu_adv 
     
    289214            ! to give asum = 1.0 after ridging. 
    290215 
    291             divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice  ! asum found in ridgeprep 
    292  
    293             IF (divu_adv(ji,jj) .LT. 0.0) & 
    294                closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj)) 
     216            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) / rdt_ice  ! asum found in ridgeprep 
     217 
     218            IF(divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
    295219 
    296220            ! 2.3 opning 
     
    306230      ! 3) Ridging iteration 
    307231      !-----------------------------------------------------------------------------! 
    308       niter = 1                 ! iteration counter 
     232      niter           = 1                 ! iteration counter 
    309233      iterate_ridging = 1 
    310  
    311234 
    312235      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 
     
    361284         !-----------------------------------------------------------------------------! 
    362285 
    363          CALL lim_itd_me_ridgeshift (opning,    closing_gross, & 
    364             msnow_mlt, esnow_mlt) 
     286         CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 
    365287 
    366288         ! 3.4 Compute total area of ice plus open water after ridging. 
     
    415337      ! Update fresh water and heat fluxes due to snow melt. 
    416338 
    417       dti = 1.0/rdt_ice 
     339      dti = 1._wp / rdt_ice 
    418340 
    419341      asum_error = .false.  
     
    432354            ! 5) Heat, salt and freshwater fluxes 
    433355            !-----------------------------------------------------------------------------! 
    434             ! fresh water source for ocean 
    435             fmmec(ji,jj)      = fmmec(ji,jj)      + msnow_mlt(ji,jj)*dti   
    436  
    437             ! heat sink for ocean 
    438             fhmec(ji,jj)      = fhmec(ji,jj)      + esnow_mlt(ji,jj)*dti 
     356            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * dti     ! fresh water source for ocean 
     357            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * dti     ! heat sink for ocean 
    439358 
    440359         END DO 
     
    477396      !----------------- 
    478397 
    479       d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 
    480       d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 
    481       d_a_i_trp(:,:,:)    = a_i(:,:,:)   - old_a_i(:,:,:) 
    482       d_v_s_trp(:,:,:)    = v_s(:,:,:)   - old_v_s(:,:,:)   
    483       d_v_i_trp(:,:,:)    = v_i(:,:,:)   - old_v_i(:,:,:)    
    484       d_e_s_trp(:,:,:,:)  = e_s(:,:,:,:) - old_e_s(:,:,:,:)   
    485       d_e_i_trp(:,:,:,:)  = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    486       d_oa_i_trp(:,:,:)   = oa_i(:,:,:)  - old_oa_i(:,:,:) 
    487       d_smv_i_trp(:,:,:)   = 0.0 
     398      d_u_ice_dyn(:,:)     = u_ice(:,:)    - old_u_ice(:,:) 
     399      d_v_ice_dyn(:,:)     = v_ice(:,:)    - old_v_ice(:,:) 
     400      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
     401      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
     402      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
     403      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
     404      d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
     405      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
     406      d_smv_i_trp(:,:,:)   = 0._wp 
    488407      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    489408         d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
     
    534453      e_i(:,:,:,:)  = old_e_i(:,:,:,:) 
    535454      oa_i(:,:,:)   = old_oa_i(:,:,:) 
    536       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    537          smv_i(:,:,:)  = old_smv_i(:,:,:) 
     455      IF ( ( num_sal == 2 ) .OR. ( num_sal == 4 ) )   smv_i(:,:,:)  = old_smv_i(:,:,:) 
    538456 
    539457      !----------------------------------------------------! 
     
    582500      END DO 
    583501 
    584       IF(.NOT. wrk_release(2, 1,2,3,4,5,6,7,8))THEN 
    585          CALL ctl_stop('lim_itd_me : failed to release workspace arrays.') 
     502      IF( .NOT. wrk_release(2, 1,2,3,4,5,6,7,8) ) THEN 
     503         CALL ctl_stop( 'lim_itd_me : failed to release workspace arrays.' ) 
    586504      END IF 
    587  
     505      ! 
    588506   END SUBROUTINE lim_itd_me 
    589507 
    590    !=============================================================================== 
    591  
    592    SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 
    593  
     508 
     509   SUBROUTINE lim_itd_me_icestrength (kstrngth) 
    594510      !!---------------------------------------------------------------------- 
    595511      !!                ***  ROUTINE lim_itd_me_icestrength *** 
     
    614530      USE wrk_nemo, ONLY: zworka => wrk_2d_1 !: temporary array used here 
    615531      ! 
    616       !! * Arguments 
    617  
    618532      INTEGER, INTENT(in) :: & 
    619533         kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979) 
     
    630544         zp,       &         !: temporary ice strength  
    631545         zdummy 
    632  
    633       IF(.NOT. wrk_use(2, 1))THEN 
    634          CALL ctl_stop('lim_itd_me_icestrength : requested workspace array unavailable.') 
    635          RETURN 
     546      !!---------------------------------------------------------------------- 
     547 
     548      IF( .NOT. wrk_use(2, 1) ) THEN 
     549         CALL ctl_stop('lim_itd_me_icestrength : requested workspace array unavailable.')   ;   RETURN 
    636550      END IF 
    637551 
     
    639553      ! 1) Initialize 
    640554      !------------------------------------------------------------------------------! 
    641       strength(:,:) = 0.0 
     555      strength(:,:) = 0._wp 
    642556 
    643557      !------------------------------------------------------------------------------! 
     
    655569               DO ji = 1, jpi 
    656570 
    657                   IF(     ( a_i(ji,jj,jl)    .GT. epsi11 )                     & 
    658                      .AND. ( athorn(ji,jj,jl) .GT. 0.0    ) ) THEN 
     571                  IF(  a_i(ji,jj,jl)    .GT. epsi11  .AND.                    & 
     572                       athorn(ji,jj,jl) .GT. 0._wp  ) THEN 
    659573                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    660574                     !---------------------------- 
    661575                     ! PE loss from deforming ice 
    662576                     !---------------------------- 
    663                      strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) *    & 
    664                         hi * hi 
     577                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi 
    665578 
    666579                     !-------------------------- 
    667580                     ! PE gain from rafting ice 
    668581                     !-------------------------- 
    669                      strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) & 
    670                         * hi * hi 
     582                     strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) * hi * hi 
    671583 
    672584                     !---------------------------- 
     
    752664                     + strength(ji,jj+1) * tms(ji,jj+1)     
    753665 
    754                   zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj)            & 
    755                      + tms(ji,jj-1) + tms(ji,jj+1) 
     666                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 
    756667                  zworka(ji,jj) = zworka(ji,jj) / zw1 
    757668               ELSE 
     
    811722   END SUBROUTINE lim_itd_me_icestrength 
    812723 
    813    !=============================================================================== 
    814  
    815    SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) 
    816  
     724 
     725   SUBROUTINE lim_itd_me_ridgeprep 
    817726      !!---------------------------------------------------------------------! 
    818727      !!                ***  ROUTINE lim_itd_me_ridgeprep *** 
     
    824733      !! participating in ridging and of the resulting ridges. 
    825734      !! 
    826       !! ** Arguments : 
    827       !! 
    828       !! ** External :  
    829       !! 
    830735      !!---------------------------------------------------------------------! 
    831       !! * Arguments 
    832  
    833736      INTEGER :: & 
    834737         ji,jj,  &          ! horizontal indices 
     
    850753         zworka            ! temporary array used here 
    851754 
    852       REAL(wp)            ::     & 
    853          zdummy,                 & 
    854          epsi06 = 1.0e-6 
     755      REAL(wp)            ::   zdummy 
    855756 
    856757      !------------------------------------------------------------------------------! 
     
    11041005   END SUBROUTINE lim_itd_me_ridgeprep 
    11051006 
    1106    !=============================================================================== 
    1107  
    1108    SUBROUTINE lim_itd_me_ridgeshift(opning,    closing_gross,       & 
    1109       msnow_mlt, esnow_mlt) ! (subroutine 4/6) 
    1110  
     1007 
     1008   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 
    11111009      !!----------------------------------------------------------------------------- 
    11121010      !!                ***  ROUTINE lim_itd_me_icestrength *** 
     
    11191017      !! and add to thicker ice categories. 
    11201018      !! 
    1121       !! ** Arguments : 
    1122       !! 
    1123       !! ** Inputs / Ouputs :  
    1124       !! 
    1125       !! ** External :  
    1126       !! 
    1127  
     1019      !!----------------------------------------------------------------------------- 
    11281020      REAL (wp), DIMENSION(jpi,jpj), INTENT(IN)   :: & 
    11291021         opning,         & ! rate of opening due to divergence/shear 
     
    12191111      REAL(wp) ::        & 
    12201112         zeps          , & 
    1221          epsi10        , & 
    12221113         zindb              ! switch for the presence of ridge poros or not 
    12231114 
     
    12351126 
    12361127      zeps   = 1.0d-20 
    1237       epsi10 = 1.0d-10 
    12381128 
    12391129      !------------------------------------------------------------------------------- 
     
    16561546   END SUBROUTINE lim_itd_me_ridgeshift 
    16571547 
    1658    !============================================================================== 
    1659  
    1660    SUBROUTINE lim_itd_me_asumr !(subroutine 5/6) 
    1661  
     1548 
     1549   SUBROUTINE lim_itd_me_asumr 
    16621550      !!----------------------------------------------------------------------------- 
    16631551      !!                ***  ROUTINE lim_itd_me_asumr *** 
     
    16721560      !! included in the sum instead of being computed as a residual.  
    16731561      !! 
    1674       !! ** Arguments : 
    1675  
    1676       INTEGER :: ji, jj, jl 
    1677  
    1678       !----------------------------------------------------------------- 
    1679       ! open water 
    1680       !----------------------------------------------------------------- 
    1681  
    1682       DO jj = 1, jpj 
    1683          DO ji = 1, jpi 
    1684             asum(ji,jj) = ato_i(ji,jj) 
    1685          END DO 
     1562      !!----------------------------------------------------------------------------- 
     1563      INTEGER ::   jl   ! dummy loop index 
     1564      !!----------------------------------------------------------------------------- 
     1565      ! 
     1566      asum(:,:) = ato_i(:,:)                    ! open water 
     1567      ! 
     1568      DO jl = 1, jpl                            ! ice categories 
     1569         asum(:,:) = asum(:,:) + a_i(:,:,jl) 
    16861570      END DO 
    1687  
    1688       !----------------------------------------------------------------- 
    1689       ! ice categories 
    1690       !----------------------------------------------------------------- 
    1691  
    1692       DO jl = 1, jpl 
    1693          DO jj= 1, jpj 
    1694             DO ji = 1, jpi 
    1695                asum(ji,jj) = asum(ji,jj) + a_i(ji,jj,jl) 
    1696             END DO !ji 
    1697          END DO !jj 
    1698       END DO ! jl 
    1699  
     1571      ! 
    17001572   END SUBROUTINE lim_itd_me_asumr 
    17011573 
    1702    !============================================================================== 
    1703  
    1704    SUBROUTINE lim_itd_me_init ! (subroutine 6/6) 
     1574 
     1575   SUBROUTINE lim_itd_me_init 
    17051576      !!------------------------------------------------------------------- 
    17061577      !!                   ***  ROUTINE lim_itd_me_init *** 
     
    17241595         brinstren_swi 
    17251596      !!------------------------------------------------------------------- 
    1726  
    1727       ! Define the initial parameters 
    1728       ! ------------------------- 
    1729       REWIND( numnam_ice ) 
     1597      ! 
     1598      REWIND( numnam_ice )                   ! read namiceitdme namelist 
    17301599      READ  ( numnam_ice , namiceitdme) 
    1731       IF (lwp) THEN 
     1600      ! 
     1601      IF (lwp) THEN                          ! control print 
    17321602         WRITE(numout,*) 
    17331603         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 
     
    17501620         WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi 
    17511621      ENDIF 
    1752  
     1622      ! 
    17531623   END SUBROUTINE lim_itd_me_init 
    17541624 
    1755    !============================================================================== 
    17561625 
    17571626   SUBROUTINE lim_itd_me_zapsmall 
     
    17691638      !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code 
    17701639      !!------------------------------------------------------------------- 
    1771  
    17721640      INTEGER ::   & 
    17731641         ji,jj,  & ! horizontal indices 
     
    17841652 
    17851653 
    1786       REAL(wp) :: & 
    1787          xtmp      ! temporary variable 
     1654      REAL(wp) ::   xtmp      ! temporary variable 
     1655      !!------------------------------------------------------------------- 
    17881656 
    17891657      DO jl = 1, jpl 
     
    18881756 
    18891757      END DO                 ! jl  
    1890  
     1758      ! 
    18911759   END SUBROUTINE lim_itd_me_zapsmall 
    18921760 
    18931761#else 
    1894    !!====================================================================== 
    1895    !!                       ***  MODULE limitd_me    *** 
    1896    !!                            no sea ice model 
    1897    !!====================================================================== 
    1898  
     1762   !!---------------------------------------------------------------------- 
     1763   !!   Default option         Empty module          NO LIM-3 sea-ice model 
     1764   !!---------------------------------------------------------------------- 
    18991765CONTAINS 
    1900  
    19011766   SUBROUTINE lim_itd_me           ! Empty routines 
    19021767   END SUBROUTINE lim_itd_me 
    19031768   SUBROUTINE lim_itd_me_icestrength 
    19041769   END SUBROUTINE lim_itd_me_icestrength 
    1905    SUBROUTINE lim_itd_me_ridgeprep 
    1906    END SUBROUTINE lim_itd_me_ridgeprep 
    1907    SUBROUTINE lim_itd_me_ridgeshift 
    1908    END SUBROUTINE lim_itd_me_ridgeshift 
    1909    SUBROUTINE lim_itd_me_asumr 
    1910    END SUBROUTINE lim_itd_me_asumr 
    19111770   SUBROUTINE lim_itd_me_sort 
    19121771   END SUBROUTINE lim_itd_me_sort 
     
    19161775   END SUBROUTINE lim_itd_me_zapsmall 
    19171776#endif 
     1777   !!====================================================================== 
    19181778END MODULE limitd_me 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r2528 r2601  
    2828   USE in_out_manager   ! I/O manager 
    2929   USE prtctl           ! Print control 
     30   USE cpl_oasis3, ONLY : lk_cpl 
    3031 
    3132   IMPLICIT NONE 
     
    4041   REAL(wp)  ::   rone   = 1._wp 
    4142 
    42    REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2] 
    43    REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s] 
    44  
    45    REAL(wp), DIMENSION(jpi,jpj) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2] 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s] 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice)  
    4646 
    4747   !! * Substitutions 
     
    7676      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
    7777      !!--------------------------------------------------------------------- 
     78      USE wrk_nemo, ONLY:   wrk_release, wrk_use 
     79      USE wrk_nemo, ONLY:   zfcm1 => wrk_2d_1 , zfcm2 => wrk_2d_2   ! 2D workspace 
     80      USE wrk_nemo, ONLY:   wrk_3d_4, wrk_3d_5                      ! 3D workspace 
     81      ! 
    7882      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    79       !! 
     83      ! 
    8084      INTEGER  ::   ji, jj           ! dummy loop indices 
     85      INTEGER  ::   ierr             ! local integer 
    8186      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    8287      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    83       REAL(wp) ::   zinda            ! switch for testing the values of ice concentration 
    84       REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    85       REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface 
    86       REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
    87 #if defined key_coupled     
    88       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky 
    89       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky 
    90 #endif 
     88      REAL(wp) ::   zinda, zfons, zpme              ! local scalars 
     89      ! 
     90      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace 
    9191      !!--------------------------------------------------------------------- 
     92 
     93      IF( .NOT.wrk_use(2, 1,2) .OR. .NOT.wrk_use(3, 4,5) ) THEN 
     94         CALL ctl_stop( 'lim_sbc_flx_2 : requested workspace arrays unavailable.' )   ;   RETURN 
     95      ENDIF 
     96      ! Set-up pointers to sub-arrays of 3d workspaces 
     97      zalb  => wrk_3d_4(:,:,1:jpl) 
     98      zalbp => wrk_3d_5(:,:,1:jpl) 
    9299 
    93100      IF( kt == nit000 ) THEN 
     
    97104         ! 
    98105         r1_rdtice = 1. / rdt_ice 
     106         ! 
     107         ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                        & 
     108            &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj) , STAT=ierr ) 
     109         ! 
     110         IF( ierr /= 0 ) THEN 
     111            CALL ctl_stop( 'lim_sbc_flx: failed to allocate arrays.' )   ;   RETURN 
     112         END IF 
    99113         ! 
    100114         soce_0(:,:) = soce 
     
    168182            ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead 
    169183 
    170             IF ( num_sal .EQ. 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + & 
     184            IF ( num_sal == 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + & 
    171185               fhbri(ji,jj) ! new contribution due to brine drainage  
    172186 
     
    181195 
    182196            !!gm   this IF prevents the vertorisation of the whole loop 
    183             IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN 
     197            IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    184198               WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    185199               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
     
    274288      !   Storing the transmitted variables           ! 
    275289      !-----------------------------------------------! 
    276  
    277290      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction             
    278291      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                       
    279292 
    280 #if defined key_coupled             
    281293      !------------------------------------------------! 
    282294      !    Computation of snow/ice and ocean albedo    ! 
    283295      !------------------------------------------------! 
    284       zalb  (:,:,:) = 0.e0 
    285       zalbp (:,:,:) = 0.e0 
    286  
    287       CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb ) 
    288  
    289       alb_ice(:,:,:) =  0.5 * zalbp(:,:,:) + 0.5 * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys) 
    290 #endif 
     296      IF( lk_cpl ) THEN          ! coupled case 
     297         CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo 
     298         ! 
     299         alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys) 
     300      ENDIF 
    291301 
    292302      IF(ln_ctl) THEN 
     
    296306         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 
    297307      ENDIF 
     308      ! 
     309      IF(  .NOT. wrk_release(2, 1)  .OR.  .NOT. wrk_release(3, 4,5)  ) THEN 
     310         CALL ctl_stop( 'lim_sbc_flx_2 : failed to release workspace arrays.' ) 
     311      END IF 
    298312      !  
    299313   END SUBROUTINE lim_sbc_flx 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r2528 r2601  
    44   !! LIM transport ice model : sea-ice advection/diffusion 
    55   !!====================================================================== 
     6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code 
     7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations 
     8   !!---------------------------------------------------------------------- 
    69#if defined key_lim3 
    710   !!---------------------------------------------------------------------- 
     
    912   !!---------------------------------------------------------------------- 
    1013   !!   lim_trp      : advection/diffusion process of sea ice 
    11    !!   lim_trp_init : initialization and namelist read 
    12    !!---------------------------------------------------------------------- 
    13    USE phycst 
    14    USE dom_oce 
     14   !!---------------------------------------------------------------------- 
     15   USE phycst          ! physical constant 
     16   USE dom_oce         ! ocean domain 
     17   USE sbc_oce         ! ocean surface boundary condition 
     18   USE par_ice         ! LIM-3 parameter 
     19   USE dom_ice         ! LIM-3 domain 
     20   USE ice             ! LIM-3 variables 
     21   USE limadv          ! LIM-3 advection 
     22   USE limhdf          ! LIM-3 horizontal diffusion 
    1523   USE in_out_manager  ! I/O manager 
    16    USE sbc_oce         ! Surface boundary condition: ocean fields 
    17    USE dom_ice 
    18    USE ice 
    19    USE limadv 
    20    USE limhdf 
    21    USE lbclnk 
    22    USE lib_mpp 
    23    USE par_ice 
     24   USE lbclnk          ! lateral boundary conditions -- MPP exchanges 
     25   USE lib_mpp         ! MPP library 
    2426   USE prtctl          ! Print control 
    2527 
     
    2729   PRIVATE 
    2830 
    29    !! * Routine accessibility 
    30    PUBLIC lim_trp       ! called by ice_step 
    31  
    32    !! * Shared module variables 
    33    REAL(wp), PUBLIC  ::   &  !: 
    34       bound  = 0.e0           !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
    35  
    36    !! * Module variables 
    37    REAL(wp)  ::           &  ! constant values 
    38       epsi06 = 1.e-06  ,  & 
    39       epsi03 = 1.e-03  ,  & 
    40       epsi16 = 1.e-16  ,  & 
    41       rzero  = 0.e0    ,  & 
    42       rone   = 1.e0    ,  & 
    43       zeps10 = 1.e-10 
     31   PUBLIC   lim_trp    ! called by ice_step 
     32 
     33   REAL(wp), PUBLIC ::   bound  = 0._wp   !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
     34 
     35   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
     36   REAL(wp)  ::   epsi03 = 1.e-03_wp   
     37   REAL(wp)  ::   zeps10 = 1.e-10_wp   
     38   REAL(wp)  ::   epsi16 = 1.e-16_wp 
     39   REAL(wp)  ::   rzero  = 0._wp    
     40   REAL(wp)  ::   rone   = 1._wp 
    4441 
    4542   !! * Substitution 
    4643#  include "vectopt_loop_substitute.h90" 
    4744   !!---------------------------------------------------------------------- 
    48    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     45   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    4946   !! $Id$ 
    50    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    51    !!---------------------------------------------------------------------- 
    52  
     47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     48   !!---------------------------------------------------------------------- 
    5349CONTAINS 
    5450 
     
    6460      !! 
    6561      !! ** action : 
    66       !! 
    67       !! History : 
    68       !!   1.0  !  00-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code 
    69       !!        !  01-05 (G. Madec, R. Hordoir) opa norm 
    70       !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp 
    71       !!   3.0  !  05-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations 
    7262      !!--------------------------------------------------------------------- 
    73       INTEGER, INTENT(in) ::   kt     ! number of iteration 
    74       !! * Local Variables 
    75       INTEGER  ::   ji, jj, jk, jl, layer, &  ! dummy loop indices 
    76          initad           ! number of sub-timestep for the advection 
    77       INTEGER  ::   ji_maxu, ji_maxv, jj_maxu, jj_maxv 
    78  
    79       REAL(wp) ::  &                               
    80          zindb  ,  & 
    81          zindsn ,  & 
    82          zindic ,  & 
    83          zusvosn,  & 
    84          zusvoic,  & 
    85          zvbord ,  & 
    86          zcfl   ,  & 
    87          zusnit ,  & 
    88          zrtt, zsal, zage, & 
    89          zbigval, ze, & 
    90          zmaxu, zmaxv 
    91  
    92       REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace 
    93          zui_u , zvi_v , zsm   ,         & 
    94          zs0at, zs0ow 
    95  
    96       REAL(wp), DIMENSION(jpi,jpj,jpl):: &  ! temporary workspace 
    97          zs0ice, zs0sn, zs0a   ,         & 
    98          zs0c0 ,                         & 
    99          zs0sm , zs0oi 
    100  
    101       ! MHE Multilayer heat content 
    102       REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl)  ::   &  ! temporary workspace 
    103          zs0e 
    104  
    105       !--------------------------------------------------------------------- 
    106  
    107       IF( numit == nstart  )   CALL lim_trp_init      ! Initialization (first time-step only) 
    108  
     63      USE wrk_nemo, ONLY:   wrk_use, wrk_release 
     64      USE wrk_nemo, ONLY:   zs0at => wrk_2d_1 , zsm => wrk_2d_2 , zs0ow  => wrk_2d_3      ! 2D workspace 
     65      USE wrk_nemo, ONLY:   wrk_3d_1, wrk_3d_2, wrk_3d_3, wrk_3d_4, wrk_3d_5, wrk_3d_6    ! 3D workspace 
     66      USE wrk_nemo, ONLY:   wrk_4d_1                                                      ! 4D workspace 
     67      ! 
     68      INTEGER, INTENT(in) ::   kt   ! number of iteration 
     69      ! 
     70      INTEGER  ::   ji, jj, jk, jl, layer   ! dummy loop indices 
     71      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
     72      REAL(wp) ::   zindb  , zindsn , zindic      ! local scalar 
     73      REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      - 
     74      REAL(wp) ::   zcfl , zusnit , zrtt          !   -      - 
     75      REAL(wp) ::   ze   , zsal   , zage          !   -      - 
     76      ! 
     77      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi   ! 3D pointer 
     78      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zs0e                                         ! 4D pointer 
     79      !!--------------------------------------------------------------------- 
     80 
     81      IF( .NOT.wrk_use( 2, 1,2,3,4,5 ) ) THEN 
     82         CALL ctl_stop( 'lim_trp : requested workspace arrays unavailable.' )   ;   RETURN 
     83      END IF 
     84 
     85      zs0ice => wrk_3d_1(:,:,1:jpl)   ;   zs0a  => wrk_3d_3   ;   zs0sm => wrk_3d_3 
     86      zs0sn  => wrk_3d_2(:,:,1:jpl)   ;   zs0c0 => wrk_3d_3   ;   zs0oi => wrk_3d_3 
     87      zs0e   => wrk_4d_1(:,:,1:jkmax,1:jpl) 
     88 
     89 
     90      IF( numit == nstart .AND. lwp ) THEN 
     91         WRITE(numout,*) 
     92         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport ' 
     93         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn 
     94         ENDIF 
     95         WRITE(numout,*) '~~~~~~~~~~~~' 
     96      ENDIF 
     97       
    10998      zsm(:,:) = area(:,:) 
    11099 
    111       IF( ln_limdyn ) THEN 
    112          IF( kt == nit000 .AND. lwp ) THEN 
    113             WRITE(numout,*) ' lim_trp : Ice Advection' 
    114             WRITE(numout,*) ' ~~~~~~~' 
    115          ENDIF 
    116  
    117          !-----------------------------------------------------------------------------! 
    118          ! 1) CFL Test                                                              
    119          !-----------------------------------------------------------------------------! 
    120  
    121          !------------------------------------------ 
    122          ! ice velocities at ocean U- and V-points  
    123          !------------------------------------------ 
    124  
    125          ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.         
    126          zvbord = 1.0 + ( 1.0 - bound ) 
    127          DO jj = 1, jpjm1 
    128             DO ji = 1, fs_jpim1 
    129                zui_u(ji,jj) = u_ice(ji,jj) 
    130                zvi_v(ji,jj) = v_ice(ji,jj) 
    131             END DO 
    132          END DO 
    133  
    134          ! Lateral boundary conditions 
    135          CALL lbc_lnk( zui_u, 'U', -1. ) 
    136          CALL lbc_lnk( zvi_v, 'V', -1. ) 
     100      !                             !-------------------------------------! 
     101      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   ! 
     102         !                          !-------------------------------------! 
     103         ! 
    137104 
    138105         !------------------------- 
     106         ! transported fields                                         
     107         !------------------------- 
     108         ! Snow vol, ice vol, salt and age contents, area 
     109         zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area  
     110         DO jl = 1, jpl 
     111            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume 
     112            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume 
     113            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area 
     114            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content 
     115            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content 
     116            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content 
     117            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content 
     118         END DO 
     119 
     120         !-------------------------- 
     121         ! Advection of Ice fields  (Prather scheme)                                             
     122         !-------------------------- 
     123         ! If ice drift field is too fast, use an appropriate time step for advection.          
    139124         ! CFL test for stability 
    140          !------------------------- 
    141  
    142          zcfl  = 0.e0 
    143          zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) ) 
    144          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) ) 
    145  
    146          zmaxu = 0.0 
    147          zmaxv = 0.0 
    148          DO ji = 1, jpim1 
    149             DO jj = 1, jpjm1 
    150                IF ( (ABS(zui_u(ji,jj)) .GT. zmaxu) ) THEN  
    151                   zmaxu = MAX(zui_u(ji,jj), zmaxu ) 
    152                   ji_maxu = ji 
    153                   jj_maxu = jj 
    154                ENDIF 
    155                IF ( (ABS(zvi_v(ji,jj)) .GT. zmaxv) ) THEN  
    156                   zmaxv = MAX(zvi_v(ji,jj), zmaxv ) 
    157                   ji_maxv = ji 
    158                   jj_maxv = jj 
    159                ENDIF 
    160             END DO 
    161          END DO 
    162  
    163          IF (lk_mpp ) CALL mpp_max(zcfl) 
    164  
    165          IF ( zcfl > 0.5 .AND. lwp ) & 
    166             WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl 
    167  
    168          !-----------------------------------------------------------------------------! 
    169          ! 2) Computation of transported fields                                         
    170          !-----------------------------------------------------------------------------! 
    171  
    172          !------------------------------------------------------ 
    173          ! 1.1) Snow vol, ice vol, salt and age contents, area 
    174          !------------------------------------------------------ 
    175  
    176          zs0ow (:,:) =  ato_i(:,:)    * area(:,:)           ! Open water area  
    177          DO jl = 1, jpl  !sum over thickness categories 
    178             ! area -> is the unmasked and masked area of T-grid cell 
    179             zs0sn (:,:,jl) =  v_s(:,:,jl)    * area(:,:)    ! Snow volume. 
    180             zs0ice(:,:,jl) =  v_i(:,:,jl)    * area(:,:)    ! Ice volume. 
    181             zs0a  (:,:,jl) =  a_i(:,:,jl)    * area(:,:)    ! Ice area 
    182             zs0sm (:,:,jl) =  smv_i(:,:,jl)  * area(:,:)    ! Salt content 
    183             zs0oi (:,:,jl) =  oa_i (:,:,jl)  * area(:,:)    ! Age content 
    184  
    185             !---------------------------------- 
    186             ! 1.2) Ice and snow heat contents 
    187             !---------------------------------- 
    188  
    189             zs0c0 (:,:,jl)     = e_s(:,:,1,jl)              ! Snow heat cont. 
    190             DO jk = 1, nlay_i 
    191                zs0e(:,:,jk,jl) = e_i(:,:,jk,jl)             ! Ice heat content 
    192             END DO 
    193          END DO 
    194  
    195          !-----------------------------------------------------------------------------! 
    196          ! 3) Advection of Ice fields                                               
    197          !-----------------------------------------------------------------------------! 
    198  
    199          ! If ice drift field is too fast, use an appropriate time step for advection.          
     125         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 
     126         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 
     127         IF(lk_mpp )   CALL mpp_max( zcfl ) 
     128!!gm more readability: 
     129!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp 
     130!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp 
     131!         ENDIF 
     132!!gm end 
    200133         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 
    201134         zusnit = 1.0 / REAL( initad )  
    202  
    203          IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==! 
     135         IF( zcfl > 0.5 .AND. lwp )   & 
     136            WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
     137               &                        ': the ice time stepping is split in two' 
     138 
     139         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    204140            DO jk = 1,initad 
    205                !--- ice open water area 
    206                CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ow(:,:) , sxopw(:,:) , &  
    207                   sxxopw(:,:), syopw(:,:) , &  
    208                   syyopw(:,:), sxyopw(:,:) ) 
    209                CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ow(:,:) , sxopw (:,:) , & 
    210                   sxxopw(:,:), syopw (:,:) , &  
    211                   syyopw(:,:), sxyopw(:,:) ) 
     141               CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
     142                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     143               CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   & 
     144                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    212145               DO jl = 1, jpl 
    213                   !--- ice volume  --- 
    214                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &  
    215                      sxxice(:,:,jl) , syice (:,:,jl) , &  
    216                      syyice(:,:,jl) , sxyice(:,:,jl) ) 
    217                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
    218                      sxxice(:,:,jl) , syice (:,:,jl) , &  
    219                      syyice(:,:,jl) , sxyice(:,:,jl) ) 
    220                   !--- snow volume  --- 
    221                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
    222                      sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
    223                      syysn (:,:,jl) , sxysn (:,:,jl) ) 
    224                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
    225                      sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
    226                      syysn (:,:,jl) , sxysn (:,:,jl) ) 
    227                   !--- ice salinity --- 
    228                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 
    229                      sxxsal(:,:,jl) , sysal (:,:,jl) , & 
    230                      syysal(:,:,jl) , sxysal(:,:,jl)  ) 
    231                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 
    232                      sxxsal(:,:,jl) , sysal (:,:,jl) , & 
    233                      syysal(:,:,jl) , sxysal(:,:,jl)  ) 
    234                   !--- ice age      ---      
    235                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 
    236                      sxxage(:,:,jl) , syage (:,:,jl) , & 
    237                      syyage(:,:,jl) , sxyage(:,:,jl)  ) 
    238                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 
    239                      sxxage(:,:,jl) , syage (:,:,jl) , & 
    240                      syyage(:,:,jl) , sxyage(:,:,jl)  ) 
    241                   !--- ice concentrations --- 
    242                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
    243                      sxxa  (:,:,jl) , sya   (:,:,jl) , &  
    244                      syya  (:,:,jl) , sxya  (:,:,jl)  ) 
    245                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , &  
    246                      sxxa  (:,:,jl) , sya   (:,:,jl) , &  
    247                      syya  (:,:,jl) , sxya  (:,:,jl)  ) 
    248                   !--- ice / snow thermal energetic contents --- 
    249                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &  
    250                      sxxc0 (:,:,jl) , syc0  (:,:,jl) , & 
    251                      syyc0 (:,:,jl) , sxyc0 (:,:,jl)  ) 
    252                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
    253                      sxxc0 (:,:,jl) , syc0  (:,:,jl) , & 
    254                      syyc0 (:,:,jl) , sxyc0 (:,:,jl)  ) 
    255                   DO layer = 1, nlay_i 
    256                      CALL lim_adv_x( zusnit, zui_u, rone , zsm, & 
    257                         zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , &  
    258                         sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , & 
    259                         syye(:,:,layer,jl) , sxye(:,:,layer,jl) ) 
    260                      CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, &  
    261                         zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , &  
    262                         sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , & 
    263                         syye(:,:,layer,jl) , sxye(:,:,layer,jl) ) 
     146                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     147                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     148                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     149                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     150                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     151                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     152                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     153                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     154                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     155                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     156                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     157                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     158                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---      
     159                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     160                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     161                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     162                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     163                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     164                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &  
     165                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     166                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     167                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     168                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     169                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     170                  DO layer = 1, nlay_i                                                           !--- ice heat contents --- 
     171                     CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     172                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
     173                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
     174                     CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     175                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
     176                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
    264177                  END DO 
    265178               END DO 
     
    267180         ELSE 
    268181            DO jk = 1, initad 
    269                !--- ice volume  --- 
    270                CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ow (:,:) , sxopw (:,:) , & 
    271                   sxxopw(:,:) , syopw (:,:) , &  
    272                   syyopw(:,:) , sxyopw(:,:) ) 
    273                CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ow (:,:) , sxopw (:,:) , &  
    274                   sxxopw(:,:) , syopw (:,:) , & 
    275                   syyopw(:,:) , sxyopw(:,:) ) 
     182               CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
     183                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     184               CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   & 
     185                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    276186               DO jl = 1, jpl 
    277                   !--- ice volume  --- 
    278                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
    279                      sxxice(:,:,jl) , syice (:,:,jl) , &  
    280                      syyice(:,:,jl) , sxyice(:,:,jl) ) 
    281                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &  
    282                      sxxice(:,:,jl) , syice (:,:,jl) , & 
    283                      syyice(:,:,jl) , sxyice(:,:,jl) ) 
    284                   !--- snow volume  --- 
    285                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &  
    286                      sxxsn (:,:,jl) , sysn  (:,:,jl) , &  
    287                      syysn (:,:,jl) , sxysn (:,:,jl) ) 
    288                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &  
    289                      sxxsn (:,:,jl) , sysn  (:,:,jl) , &  
    290                      syysn (:,:,jl) , sxysn (:,:,jl) ) 
    291                   !--- ice salinity --- 
    292                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 
    293                      sxxsal(:,:,jl) , sysal (:,:,jl) , & 
    294                      syysal(:,:,jl) , sxysal(:,:,jl) ) 
    295                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , & 
    296                      sxxsal(:,:,jl) , sysal (:,:,jl) , & 
    297                      syysal(:,:,jl) , sxysal(:,:,jl) ) 
    298                   !--- ice age      --- 
    299                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 
    300                      sxxage(:,:,jl) , syage (:,:,jl) , &  
    301                      syyage(:,:,jl) , sxyage(:,:,jl)  ) 
    302                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , & 
    303                      sxxage(:,:,jl) , syage (:,:,jl) , & 
    304                      syyage(:,:,jl) , sxyage(:,:,jl)   ) 
    305                   !--- ice concentration --- 
    306                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , &  
    307                      sxxa  (:,:,jl) , sya   (:,:,jl) , &  
    308                      syya  (:,:,jl) , sxya  (:,:,jl)  ) 
    309                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , &  
    310                      sxxa  (:,:,jl) , sya   (:,:,jl) , &  
    311                      syya  (:,:,jl) , sxya  (:,:,jl)  ) 
    312                   !--- ice / snow thermal energetic contents --- 
    313                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &  
    314                      sxxc0 (:,:,jl) , syc0  (:,:,jl) , & 
    315                      syyc0 (:,:,jl) , sxyc0 (:,:,jl)  ) 
    316                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
    317                      sxxc0 (:,:,jl) , syc0  (:,:,jl) , & 
    318                      syyc0 (:,:,jl) , sxyc0 (:,:,jl)  ) 
    319                   DO layer = 1, nlay_i 
    320                      CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0e(:,:,layer,jl) , & 
    321                         sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , & 
    322                         syye (:,:,layer,jl), sxye (:,:,layer,jl) ) 
    323                      CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0e(:,:,layer,jl) , & 
    324                         sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , & 
    325                         syye (:,:,layer,jl), sxye (:,:,layer,jl)  ) 
     187                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     188                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     189                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     190                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
     191                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     192                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     193                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     194                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
     195                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     196                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     197                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     198                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
     199 
     200                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     201                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     202                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     203                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
     204                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     205                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     206                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
     207                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
     208                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     209                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     210                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     211                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
     212                  DO layer = 1, nlay_i                                                           !--- ice heat contents --- 
     213                     CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     214                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
     215                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
     216                     CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     217                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
     218                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
    326219                  END DO 
    327  
    328220               END DO 
    329221            END DO 
     
    333225         ! Recover the properties from their contents 
    334226         !------------------------------------------- 
    335  
    336          zs0ow (:,:)       = zs0ow(:,:) / area(:,:) 
     227         zs0ow(:,:) = zs0ow(:,:) / area(:,:) 
    337228         DO jl = 1, jpl 
    338229            zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) 
     
    351242         !------------------------------------------------------------------------------! 
    352243 
     244         !-------------------------------- 
     245         !  diffusion of open water area 
     246         !-------------------------------- 
     247         zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction 
     248         DO jl = 2, jpl 
     249            zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl) 
     250         END DO 
     251         ! 
     252         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
     253         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
     254            DO ji = 1 , fs_jpim1   ! vector opt. 
     255               pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji  ,jj) ) ) )   & 
     256                  &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 
     257               pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj  ) ) ) )   & 
     258                  &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 
     259            END DO 
     260         END DO 
     261         ! 
     262         CALL lim_hdf( zs0ow (:,:) )   ! Diffusion 
     263 
    353264         !------------------------------------ 
    354          ! 4.1) diffusion of open water area 
     265         !  Diffusion of other ice variables 
    355266         !------------------------------------ 
    356  
    357          ! Compute total ice fraction 
    358          zs0at(:,:) = 0.0 
    359          DO jl = 1, jpl 
    360             DO jj = 1, jpj 
    361                DO ji = 1, jpi 
    362                   zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) ! 
    363                END DO 
    364             END DO 
    365          END DO 
    366  
    367          ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    368          DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row 
    369             DO ji = 1 , fs_jpim1   ! vector opt. 
    370                pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji  ,jj) ) ) )   & 
    371                   &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 
    372                pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji,jj  ) ) ) )   & 
    373                   &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 
    374             END DO !jj 
    375          END DO ! ji 
    376  
    377          ! Diffusion 
    378          CALL lim_hdf( zs0ow (:,:) ) 
    379  
    380          !---------------------------------------- 
    381          ! 4.2) Diffusion of other ice variables 
    382          !---------------------------------------- 
    383          DO jl = 1, jpl 
    384  
    385             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    386             DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row 
     267         DO jl = 1, jpl 
     268         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
     269            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
    387270               DO ji = 1 , fs_jpim1   ! vector opt. 
    388                   pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj,jl) ) ) )   & 
    389                      &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
    390                   pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj,jl  ) ) ) )   & 
    391                      &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
    392                END DO !jj 
    393             END DO ! ji 
     271                  pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj,jl) ) ) )   & 
     272                     &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
     273                  pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ,jl) ) ) )   & 
     274                     &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
     275               END DO 
     276            END DO 
    394277 
    395278            CALL lim_hdf( zs0ice (:,:,jl) ) 
     
    401284            DO jk = 1, nlay_i 
    402285               CALL lim_hdf( zs0e (:,:,jk,jl) ) 
    403             END DO ! jk 
    404          END DO !jl 
     286            END DO 
     287         END DO 
    405288 
    406289         !----------------------------------------- 
    407          ! 4.3) Remultiply everything by ice area 
     290         ! Remultiply everything by ice area 
    408291         !----------------------------------------- 
    409          zs0ow(:,:) = MAX(rzero, zs0ow(:,:) * area(:,:) ) 
     292         zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) ) 
    410293         DO jl = 1, jpl 
    411294            zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) )    !!bug:  est-ce utile 
     
    432315               DO jj = 1, jpj 
    433316                  DO ji = 1, jpi 
    434                      zs0e (ji,jj,jk,jl) =  & 
    435                         MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) ) 
     317                     zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) ) 
    436318                  END DO 
    437319               END DO 
     
    441323         DO jj = 1, jpj 
    442324            DO ji = 1, jpi 
    443                zs0ow (ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) ) 
    444             END DO 
    445          END DO 
    446  
    447          zs0at(:,:) = 0.0 
     325               zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) ) 
     326            END DO 
     327         END DO 
     328 
     329         zs0at(:,:) = 0._wp 
    448330         DO jl = 1, jpl 
    449331            DO jj = 1, jpj 
     
    463345         ! 5.2) Snow thickness, Ice thickness, Ice concentrations 
    464346         !--------------------------------------------------------- 
    465  
    466347         DO jj = 1, jpj 
    467348            DO ji = 1, jpi 
    468                zindb         = MAX( 0.0 , SIGN( 1.0, zs0at(ji,jj) - zeps10) ) 
    469                zs0ow(ji,jj)  = (1.0 - zindb) + zindb*MAX( zs0ow(ji,jj), 0.00 ) 
    470                ato_i(ji,jj)  = zs0ow(ji,jj) 
    471             END DO 
    472          END DO 
    473  
    474          ! Remove very small areas  
    475          DO jl = 1, jpl 
     349               zindb        = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) ) 
     350               zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 
     351               ato_i(ji,jj) = zs0ow(ji,jj) 
     352            END DO 
     353         END DO 
     354 
     355         DO jl = 1, jpl         ! Remove very small areas  
    476356            DO jj = 1, jpj 
    477357               DO ji = 1, jpi 
    478358                  zindb         = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) ) 
    479  
    480                   zs0a(ji,jj,jl)  = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 
    481                   v_s(ji,jj,jl)   = zindb * zs0sn (ji,jj,jl)  
    482                   v_i(ji,jj,jl)   = zindb * zs0ice(ji,jj,jl) 
    483  
    484                   zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) ) 
    485                   zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) ) 
    486                   zindb           = MAX( zindsn, zindic ) 
    487                   zs0a (ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration 
    488                   a_i  (ji,jj,jl) = zs0a(ji,jj,jl) 
    489                   v_s  (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 
    490                   v_i  (ji,jj,jl) = zindic * v_i(ji,jj,jl) 
     359                  ! 
     360                  zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 
     361                  v_s(ji,jj,jl)  = zindb * zs0sn (ji,jj,jl)  
     362                  v_i(ji,jj,jl)  = zindb * zs0ice(ji,jj,jl) 
     363                  ! 
     364                  zindsn         = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) ) 
     365                  zindic         = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) ) 
     366                  zindb          = MAX( zindsn, zindic ) 
     367                  zs0a(ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration 
     368                  a_i (ji,jj,jl) = zs0a(ji,jj,jl) 
     369                  v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 
     370                  v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 
    491371               END DO 
    492372            END DO 
     
    495375         DO jj = 1, jpj 
    496376            DO ji = 1, jpi 
    497                zs0at(ji,jj)       = SUM( zs0a(ji,jj,1:jpl) ) 
     377               zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) 
    498378            END DO 
    499379         END DO 
     
    503383         !---------------------- 
    504384 
    505          zbigval         =  1.0d+13 
     385         zbigval = 1.d+13 
    506386 
    507387         DO jl = 1, jpl 
     
    518398 
    519399                  ! Ice salinity and age 
    520                   zsal            = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , & 
    521                      zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * & 
    522                      v_i(ji,jj,jl) 
     400                  zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , & 
     401                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 
    523402                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    524403                     smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0 
    525404 
    526                   zage            = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &  
    527                      MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * & 
    528                      a_i(ji,jj,jl) 
     405                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &  
     406                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl) 
    529407                  oa_i (ji,jj,jl)  = zindic*zage  
    530408 
     
    583461         END DO 
    584462      ENDIF 
    585  
     463      ! 
     464      IF( .NOT.wrk_release( 2, 1,2,3,4,5 ) )   CALL ctl_stop('lim_trp_2 : failed to release workspace arrays.') 
     465      ! 
    586466   END SUBROUTINE lim_trp 
    587  
    588  
    589    SUBROUTINE lim_trp_init 
    590       !!------------------------------------------------------------------- 
    591       !!                  ***  ROUTINE lim_trp_init  *** 
    592       !! 
    593       !! ** Purpose :   initialization of ice advection parameters 
    594       !! 
    595       !! ** Method  : Read the namicetrp namelist and check the parameter  
    596       !!       values called at the first timestep (nit000) 
    597       !! 
    598       !! ** input   :   Namelist namicetrp 
    599       !! 
    600       !! history : 
    601       !!   2.0  !  03-08 (C. Ethe)  Original code 
    602       !!------------------------------------------------------------------- 
    603       NAMELIST/namicetrp/ bound 
    604       !!------------------------------------------------------------------- 
    605  
    606       ! Read Namelist namicetrp 
    607       REWIND ( numnam_ice ) 
    608       READ   ( numnam_ice  , namicetrp ) 
    609       IF(lwp) THEN 
    610          WRITE(numout,*) 
    611          WRITE(numout,*) 'lim_trp_init : Ice parameters for advection ' 
    612          WRITE(numout,*) '~~~~~~~~~~~~' 
    613          WRITE(numout,*) ' boundary conditions (0.0 no-slip, 1.0 free-slip) bound  = ', bound 
    614          WRITE(numout,*)  
    615       ENDIF 
    616  
    617    END SUBROUTINE lim_trp_init 
    618467 
    619468#else 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r2528 r2601  
    2727   PRIVATE 
    2828 
    29    !! * Accessibility 
    3029   PUBLIC lim_wri        ! routine called by lim_step.F90 
    3130 
    32    !! * Module variables 
    33    INTEGER, PARAMETER ::   &  !: 
    34       jpnoumax = 40             !: maximum number of variable for ice output 
    35    INTEGER  ::                                & 
    36       noumef          ,                       &  ! number of fields 
    37       noumefa         ,                       &  ! number of additional fields 
    38       add_diag_swi    ,                       &  ! additional diagnostics 
    39       nz                                         ! dimension for the itd field 
    40  
    41    REAL(wp)           , DIMENSION(jpnoumax) ::  & 
    42       cmulti          ,                       &  ! multiplicative constant 
    43       cadd            ,                       &  ! additive constant 
    44       cmultia         ,                       &  ! multiplicative constant 
    45       cadda                                      ! additive constant 
    46    CHARACTER(len = 35), DIMENSION(jpnoumax) ::  & 
    47       titn, titna                                ! title of the field 
    48    CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::  & 
    49       nam, nama                                  ! name of the field 
    50    CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::  & 
    51       uni, unia                                  ! unit of the field 
    52    INTEGER            , DIMENSION(jpnoumax) ::  & 
    53       nc, nca                                    ! switch for saving field ( = 1 ) or not ( = 0 ) 
    54  
    55    REAL(wp)  ::            &  ! constant values 
    56       epsi16 = 1e-16   ,  & 
    57       zzero  = 0.e0     ,  & 
    58       zone   = 1.e0 
     31   INTEGER, PARAMETER ::   jpnoumax = 40   !: maximum number of variable for ice output 
     32    
     33   INTEGER  ::   noumef             ! number of fields 
     34   INTEGER  ::   noumefa            ! number of additional fields 
     35   INTEGER  ::   add_diag_swi       ! additional diagnostics 
     36   INTEGER  ::   nz                                         ! dimension for the itd field 
     37 
     38   REAL(wp) , DIMENSION(jpnoumax) ::   cmulti         ! multiplicative constant 
     39   REAL(wp) , DIMENSION(jpnoumax) ::   cadd           ! additive constant 
     40   REAL(wp) , DIMENSION(jpnoumax) ::   cmultia        ! multiplicative constant 
     41   REAL(wp) , DIMENSION(jpnoumax) ::   cadda          ! additive constant 
     42   CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn, titna   ! title of the field 
     43   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam , nama    ! name of the field 
     44   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   uni , unia    ! unit of the field 
     45   INTEGER            , DIMENSION(jpnoumax) ::   nc  , nca     ! switch for saving field ( = 1 ) or not ( = 0 ) 
     46 
     47   REAL(wp)  ::   epsi16 = 1e-16_wp 
     48   REAL(wp)  ::   zzero  = 0._wp 
     49   REAL(wp)  ::   zone   = 1._wp 
    5950       
    6051   !!---------------------------------------------------------------------- 
    61    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     52   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    6253   !! $Id$ 
    63    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6455   !!---------------------------------------------------------------------- 
    6556CONTAINS 
     
    7970      !!  modif : 03/06/98 
    8071      !!------------------------------------------------------------------- 
    81       INTEGER, INTENT(in) :: & 
    82          kindic                 ! if kindic < 0 there has been an error somewhere 
    83  
    84       !! * Local variables 
     72      USE wrk_nemo, ONLY: wrk_release, wrk_use 
     73      USE wrk_nemo, ONLY: zfield => wrk_2d_1             ! 2D workspace 
     74      USE wrk_nemo, ONLY: wrk_3d_1, wrk_3D_2, wrk_3d_3   ! 3D workspace 
     75      ! 
     76      INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere 
     77      ! 
     78      INTEGER ::  ji, jj, jk, jl, jf, ipl ! dummy loop indices 
     79      INTEGER ::  ierr 
    8580      REAL(wp),DIMENSION(1) ::   zdept 
    86  
    87       REAL(wp) :: & 
    88          zsto, zjulian,zout, & 
    89          zindh,zinda,zindb 
    90       REAL(wp), DIMENSION(jpi,jpj,jpnoumax) :: & 
    91          zcmo,               & 
    92          zcmoa                   ! additional fields 
    93  
    94       REAL(wp), DIMENSION(jpi,jpj) ::  & 
    95          zfield 
    96  
    97       REAL(wp), DIMENSION(jpi,jpj,jpl) ::  & 
    98          zmaskitd, zoi, zei 
    99  
    100       INTEGER ::  ji, jj, jk, jl, jf, ipl ! dummy loop indices 
    101  
    102       CHARACTER(len = 40)  :: & 
    103          clhstnam, clop, & 
    104          clhstnama 
    105  
    106       INTEGER , SAVE ::      & 
    107          nice, nhorid, ndim, niter, ndepid 
    108       INTEGER , SAVE ::      & 
    109          nicea, nhorida, ndimitd 
    110       INTEGER , DIMENSION( jpij ) , SAVE ::  & 
    111          ndex51 
    112       INTEGER , DIMENSION( jpij*jpl ) , SAVE ::  & 
    113          ndexitd 
     81      REAL(wp) ::  zsto, zjulian, zout, zindh, zinda, zindb 
     82      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zcmo, zcmoa   ! additional fields 
     83      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zmaskitd, zoi, zei 
     84 
     85      CHARACTER(len = 40) ::   clhstnam, clop, clhstnama 
     86 
     87      INTEGER , SAVE ::   nice, nhorid, ndim, niter, ndepid 
     88      INTEGER , SAVE ::   nicea, nhorida, ndimitd 
     89      INTEGER , ALLOCATABLE, DIMENSION(:), SAVE ::   ndex51 
     90      INTEGER , ALLOCATABLE, DIMENSION(:), SAVE ::   ndexitd 
    11491      !!------------------------------------------------------------------- 
    11592 
    11693      ipl = jpl 
    11794 
    118       IF ( numit == nstart ) THEN  
     95      zcmo     => wrk_3d_1(:,:,1:jpnoumax) 
     96      zcmoa    => wrk_3d_2(:,:,1:jpnoumax) 
     97      zmaskitd => wrk_3d_2(:,:,1:jpl) 
     98      zoi      => wrk_3d_2(:,:,1:jpl) 
     99      zei      => wrk_3d_2(:,:,1:jpl) 
     100 
     101 
     102      IF( numit == nstart ) THEN  
     103 
     104         ALLOCATE( ndex51(jpij) , ndexitd(jpij*jpl) , STAT=ierr ) 
     105         IF( ierr /= 0 ) THEN 
     106            CALL ctl_stop( 'lim_wri : unable to allocate standard arrays' )   ;   RETURN 
     107         ENDIF 
    119108 
    120109         CALL lim_wri_init  
     
    209198 
    210199      !-- calculs des valeurs instantanees 
    211       zcmo( 1:jpi, 1:jpj, 1:jpnoumax ) = 0.0  
     200      zcmo ( 1:jpi, 1:jpj, 1:jpnoumax ) = 0.0  
    212201      zcmoa( 1:jpi, 1:jpj, 1:jpnoumax ) = 0.0  
    213202 
     
    391380      !! ** input   :   Namelist namicewri 
    392381      !! 
    393       !! history : 
    394       !!  8.5  ! 03-08 (C. Ethe) original code 
    395       !!------------------------------------------------------------------- 
    396       !! * Local declarations 
     382      !! history :  8.5  ! 03-08 (C. Ethe) original code 
     383      !!------------------------------------------------------------------- 
    397384      INTEGER ::   nf      ! ??? 
    398385 
     
    416403 
    417404      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield 
    418  
     405      ! 
    419406      NAMELIST/namiceout/ noumef, & 
    420407         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   & 
  • branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r2590 r2601  
    44   !! LIM sea-ice :   Ice thermodynamics in 1D 
    55   !!===================================================================== 
    6    !! History : 
    7    !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module 
     6   !! History :  3.0  !  2002-11  (C. Ethe)  F90: Free form and module 
    87   !!---------------------------------------------------------------------- 
    9    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    10    !! $Id$ 
    11    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    12    !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    14    USE par_ice 
     8   USE par_ice        ! LIM-3 parameters 
     9   USE in_out_manager ! I/O manager 
    1510 
    1611   IMPLICIT NONE 
     
    2217   !! * Share Module variables 
    2318   !!--------------------------- 
    24    REAL(wp) , PUBLIC ::   & !!! ** ice-thermo namelist (namicethd) ** 
    25       hmelt   = -0.15  ,  &  !: maximum melting at the bottom; active only for 
    26                                 !: one category 
    27       hicmin  = 0.2    ,  &  !: (REMOVE) 
    28       hiclim  = 0.05   ,  &  !: minimum ice thickness 
    29       amax    = 0.999  ,  &  !: maximum lead fraction 
    30       sbeta   = 1.0    ,  &  !: numerical scheme for diffusion in ice  (REMOVE) 
    31       parlat  = 0.0    ,  &  !: (REMOVE) 
    32       hakspl  = 0.5    ,  &  !: (REMOVE) 
    33       hibspl  = 0.5    ,  &  !: (REMOVE) 
    34       exld    = 2.0    ,  &  !: (REMOVE) 
    35       hakdif  = 1.0    ,  &  !: (REMOVE) 
    36       thth    = 0.2    ,  &  !: (REMOVE) 
    37       hnzst   = 0.1    ,  &  !: thick. of the surf. layer in temp. comp. 
    38       parsub  = 1.0    ,  &  !: switch for snow sublimation or not 
    39       alphs   = 1.0    ,  &  !: coef. for snow density when snow-ice formation 
    40       fraz_swi= 1.0    ,  &  !: use of frazil ice collection in function of wind (1.0) or not (0.0) 
    41       maxfrazb= 0.7    ,  &  !: maximum portion of frazil ice collecting at the ice bottom 
    42       vfrazb  = 0.41667,  &  !: threshold drift speed for collection of bottom frazil ice 
    43       Cfrazb  = 5.0          !: squeezing coefficient for collection of bottom frazil ice 
    44  
    45    REAL(wp), PUBLIC, DIMENSION(2)  ::  &  !:    
    46       hiccrit = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m) 
     19   !                                         !!! ** ice-thermo namelist (namicethd) ** 
     20   REAL(wp), PUBLIC ::   hmelt   = -0.15     !: maximum melting at the bottom; active only for one category 
     21   REAL(wp), PUBLIC ::   hicmin  = 0.2       !: (REMOVE) 
     22   REAL(wp), PUBLIC ::   hiclim  = 0.05      !: minimum ice thickness 
     23   REAL(wp), PUBLIC ::   amax    = 0.999     !: maximum lead fraction 
     24   REAL(wp), PUBLIC ::   sbeta   = 1.0       !: numerical scheme for diffusion in ice  (REMOVE) 
     25   REAL(wp), PUBLIC ::   parlat  = 0.0       !: (REMOVE) 
     26   REAL(wp), PUBLIC ::   hakspl  = 0.5       !: (REMOVE) 
     27   REAL(wp), PUBLIC ::   hibspl  = 0.5       !: (REMOVE) 
     28   REAL(wp), PUBLIC ::   exld    = 2.0       !: (REMOVE) 
     29   REAL(wp), PUBLIC ::   hakdif  = 1.0       !: (REMOVE) 
     30   REAL(wp), PUBLIC ::   thth    = 0.2       !: (REMOVE) 
     31   REAL(wp), PUBLIC ::   hnzst   = 0.1       !: thick. of the surf. layer in temp. comp. 
     32   REAL(wp), PUBLIC ::   parsub  = 1.0       !: switch for snow sublimation or not 
     33   REAL(wp), PUBLIC ::   alphs   = 1.0       !: coef. for snow density when snow-ice formation 
     34   REAL(wp), PUBLIC ::   fraz_swi= 1.0       !: use of frazil ice collection in function of wind (1.0) or not (0.0) 
     35   REAL(wp), PUBLIC ::   maxfrazb= 0.7       !: maximum portion of frazil ice collecting at the ice bottom 
     36   REAL(wp), PUBLIC ::   vfrazb  = 0.41667   !: threshold drift speed for collection of bottom frazil ice 
     37   REAL(wp), PUBLIC ::   Cfrazb  = 5.0       !: squeezing coefficient for collection of bottom frazil ice 
     38 
     39   REAL(wp), PUBLIC, DIMENSION(2) ::   hiccrit = (/0.3,0.3/)  !: ice th. for lateral accretion in the NH (SH) (m) 
    4740 
    4841   !!----------------------------- 
     
    5346   !: are the variables corresponding to 2d vectors 
    5447 
    55    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   &  !: 
    56       npb     ,   &   !: number of points where computations has to be done 
    57       npac            !: correspondance between the points (lateral accretion) 
    58  
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   &  !:  
    60       qldif_1d    ,     &  !: corresponding to the 2D var  qldif 
    61       qcmif_1d    ,     &  !: corresponding to the 2D var  qcmif 
    62       fstbif_1d   ,     &  !:    "                  "      fstric 
    63       fltbif_1d   ,     &  !:    "                  "      ffltbif 
    64       fscbq_1d    ,     &  !:    "                  "      fscmcbq 
    65       qsr_ice_1d  ,     &  !:    "                  "      qsr_ice 
    66       fr1_i0_1d   ,     &  !:    "                  "      fr1_i0 
    67       fr2_i0_1d   ,     &  !:    "                  "      fr2_i0 
    68       qnsr_ice_1d ,     &  !:    "                  "      qns_ice 
    69       qfvbq_1d    ,     &  !:    "                  "      qfvbq 
    70       t_bo_b               !:    "                  "      t_bo 
    71  
    72    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   &  !:  
    73       sprecip_1d  ,     &  !:    "                  "      sprecip 
    74       frld_1d     ,     &  !:    "                  "      frld 
    75       at_i_b      ,     &  !:    "                  "      frld 
    76       fbif_1d     ,     &  !:    "                  "      fbif 
    77       rdmicif_1d  ,     &  !:    "                  "      rdmicif 
    78       rdmsnif_1d  ,     &  !:    "                  "      rdmsnif 
    79       qlbbq_1d    ,     &  !:    "                  "      qlbsbq 
    80       dmgwi_1d    ,     &  !:    "                  "      dmgwi 
    81       dvsbq_1d    ,     &  !:    "                  "      rdvosif 
    82       dvbbq_1d    ,     &  !:    "                  "      rdvobif 
    83       dvlbq_1d    ,     &  !:    "                  "      rdvolif 
    84       dvnbq_1d    ,     &  !:    "                  "      rdvolif 
    85       dqns_ice_1d ,     &  !:    "                  "      dqns_ice 
    86       qla_ice_1d  ,     &  !:    "                  "      qla_ice 
    87       dqla_ice_1d ,     &  !:    "                  "      dqla_ice 
    88                                 ! to reintegrate longwave flux inside the ice thermodynamics 
    89 !!sm: not used      qtur_ice_1d ,     &  !:    "                  "      qtur_ice 
    90 !!sm: not used      dqtu_ice_1d ,     &  !:    "                  "      dqtu_ice 
    91 !!sm: not used      catm_ice_1d ,     &  !:    "                  "      catm_ice 
    92       tatm_ice_1d ,     &  !:    "                  "      tatm_ice 
    93 !!sm: not used      evsq_ice_1d ,     &  !:    "                  "      evsq_ice 
    94 !!sm: not used      sbud_ice_1d ,     &  !:    "                  "      sbud_ice 
    95       fsup        ,     &  !:    Energy flux sent from bottom to lateral ablation if |dhb|> 0.15 m 
    96       focea       ,     &  !:    Remaining energy in case of total ablation 
    97       i0          ,     &  !:    fraction of radiation transmitted to the ice interior 
    98       old_ht_i_b  ,     &  !:    Ice thickness at the beginnning of the time step [m] 
    99       old_ht_s_b  ,     &  !:    Snow thickness at the beginning of the time step [m] 
    100       fsbri_1d    ,     &  !:    Salt flux due to brine drainage 
    101       fhbri_1d    ,     &  !:    Heat flux due to brine drainage 
    102       fseqv_1d    ,     &  !:    Equivalent Salt flux due to ice growth/decay 
    103       dsm_i_fl_1d ,     &  !:    Ice salinity variations due to flushing 
    104       dsm_i_gd_1d ,     &  !:    Ice salinity variations due to gravity drainage 
    105       dsm_i_se_1d ,     &  !:    Ice salinity variations due to basal salt entrapment 
    106 !!sm: not used      dsm_i_la_1d ,     &  !:    Ice salinity variations due to lateral accretion     
    107       dsm_i_si_1d ,     &  !:    Ice salinity variations due to lateral accretion     
    108       hicol_b              !:    Ice collection thickness accumulated in fleads 
    109  
    110    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   &  !: 
    111       t_su_b      ,     &  !:    "                  "      t_su 
    112       a_i_b       ,     &  !:                              a_i 
    113       ht_i_b      ,     &  !:    "                  "      ht_s 
    114       ht_s_b      ,     &  !:    "                  "      ht_i 
    115       fc_su       ,     &  !:    Surface Conduction flux  
    116       fc_bo_i     ,     &  !:    Bottom  Conduction flux  
    117       dh_s_tot    ,     &  !:    Snow accretion/ablation        [m] 
    118       dh_i_surf   ,     &  !:    Ice surface accretion/ablation [m] 
    119       dh_i_bott   ,     &  !:    Ice bottom accretion/ablation  [m] 
    120       dh_snowice  ,     &  !:    Snow ice formation             [m of ice] 
    121       sm_i_b      ,     &  !:    Ice bulk salinity [ppt] 
    122       s_i_new     ,     &  !:    Salinity of new ice at the bottom 
    123       s_snowice   ,     &  !:    Salinity of new snow ice on top of the ice 
    124       o_i_b                !:    Ice age                        [days] 
    125  
    126    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &  !: 
    127       t_s_b              !: corresponding to the 2D var  t_s 
    128    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &  !: 
    129       t_i_b,            &  !: corresponding to the 2D var  t_i 
    130       s_i_b,            &  !: profiled ice salinity 
    131       q_i_b,            &  !:    Ice  enthalpy per unit volume 
    132       q_s_b                !:    Snow enthalpy per unit volume 
     48   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: number of points where computations has to be done 
     49   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: correspondance between points (lateral accretion) 
     50 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qldif_1d      !: <==> the 2D  qldif 
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qcmif_1d      !: <==> the 2D  qcmif 
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fstbif_1d     !: <==> the 2D  fstric 
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fltbif_1d     !: <==> the 2D  ffltbif 
     55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fscbq_1d      !: <==> the 2D  fscmcbq 
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d    !: <==> the 2D  qsr_ice 
     57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d     !: <==> the 2D  fr1_i0 
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d     !: <==> the 2D  fr2_i0 
     59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qnsr_ice_1d   !: <==> the 2D  qns_ice 
     60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qfvbq_1d      !: <==> the 2D  qfvbq 
     61   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_b        !: <==> the 2D  t_bo 
     62 
     63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip 
     64   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld 
     65   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld 
     66   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif 
     67   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmicif_1d    !: <==> the 2D  rdmicif 
     68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmsnif_1d    !: <==> the 2D  rdmsnif 
     69   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq 
     70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi 
     71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dvsbq_1d      !: <==> the 2D  rdvosif 
     72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dvbbq_1d      !: <==> the 2D  rdvobif 
     73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dvlbq_1d      !: <==> the 2D  rdvolif 
     74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dvnbq_1d      !: <==> the 2D  rdvolif 
     75   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqns_ice_1d   !: <==> the 2D  dqns_ice 
     76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qla_ice_1d    !: <==> the 2D  qla_ice 
     77   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqla_ice_1d   !: <==> the 2D  dqla_ice 
     78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tatm_ice_1d   !: <==> the 2D  tatm_ice 
     79   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
     80   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fsup          !: Energy flux sent from bottom to lateral ablation if |dhb|> 0.15 m 
     81   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   focea         !: Remaining energy in case of total ablation 
     82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
     83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   old_ht_i_b    !: Ice thickness at the beginnning of the time step [m] 
     84    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::  old_ht_s_b    !: Snow thickness at the beginning of the time step [m] 
     85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fsbri_1d      !: Salt flux due to brine drainage 
     86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhbri_1d      !: Heat flux due to brine drainage 
     87   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fseqv_1d      !: Equivalent Salt flux due to ice growth/decay 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_fl_1d   !: Ice salinity variations due to flushing 
     89   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_gd_1d   !: Ice salinity variations due to gravity drainage 
     90   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_se_1d   !: Ice salinity variations due to basal salt entrapment 
     91   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_si_1d   !: Ice salinity variations due to lateral accretion     
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hicol_b       !: Ice collection thickness accumulated in fleads 
     93 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_b      !: <==> the 2D  t_su 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_b       !: <==> the 2D  a_i 
     96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_b      !: <==> the 2D  ht_s 
     97   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_b      !: <==> the 2D  ht_i 
     98   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su       !: Surface Conduction flux  
     99   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i     !: Bottom  Conduction flux  
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot    !: Snow accretion/ablation        [m] 
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf   !: Ice surface accretion/ablation [m] 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott   !: Ice bottom accretion/ablation  [m] 
     103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice  !: Snow ice formation             [m of ice] 
     104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_b      !: Ice bulk salinity [ppt] 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new     !: Salinity of new ice at the bottom 
     106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_snowice   !: Salinity of new snow ice on top of the ice 
     107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   o_i_b       !: Ice age                        [days] 
     108 
     109   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_b   !: corresponding to the 2D var  t_s 
     110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_b   !: corresponding to the 2D var  t_i 
     111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   s_i_b   !: profiled ice salinity 
     112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_i_b   !:    Ice  enthalpy per unit volume 
     113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_s_b   !:    Snow enthalpy per unit volume 
    133114 
    134115   ! Clean the following ... 
    135116   ! These variables are coded for conservation checks 
    136    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   &  ! 
    137       qt_i_in   ,           &  !: ice energy summed over categories (initial) 
    138       qt_i_fin  ,           &  !: ice energy summed over categories (final) 
    139       qt_s_in, qt_s_fin  ,  &  !: snow energy summed over categories 
    140       dq_i, sum_fluxq    ,  &  !: increment of energy, sum of fluxes 
    141       fatm, foce,           &  !: atmospheric, oceanic, heat flux 
    142       cons_error, surf_error   !: conservation, surface error 
    143  
    144    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)::   &  !:  goes to trash 
    145       q_i_layer_in,         & 
    146       q_i_layer_fin,        & 
    147       dq_i_layer, radab 
    148  
    149    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   &  !: 
    150       ftotal_in  ,          &  !: initial total heat flux 
    151       ftotal_fin               !: final total heat flux 
    152  
    153    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &  !: 
    154       fc_s 
    155    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   &  !: 
    156       fc_i 
    157    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &  !: 
    158       de_s_lay 
    159    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   &  !: 
    160       de_i_lay 
    161    INTEGER , PUBLIC ::                           & 
    162       jiindex_1d   ! 1D index of debugging point 
    163  
    164    !!====================================================================== 
     117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_i_in                  !: ice energy summed over categories (initial) 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_i_fin                 !: ice energy summed over categories (final) 
     119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_s_in, qt_s_fin        !: snow energy summed over categories 
     120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dq_i, sum_fluxq          !: increment of energy, sum of fluxes 
     121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fatm, foce               !: atmospheric, oceanic, heat flux 
     122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cons_error, surf_error   !: conservation, surface error 
     123 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_i_layer_in        !: goes to trash 
     125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_i_layer_fin       !: goes to trash 
     126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dq_i_layer, radab   !: goes to trash 
     127 
     128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftotal_in    !: initial total heat flux 
     129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftotal_fin   !: final total heat flux 
     130 
     131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fc_s 
     132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fc_i 
     133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   de_s_lay 
     134   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   de_i_lay 
     135    
     136   INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point 
     137 
     138   !!---------------------------------------------------------------------- 
     139   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     140   !! $Id$ 
     141   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     142   !!---------------------------------------------------------------------- 
    165143CONTAINS 
    166144 
     
    169147      !!                ***  ROUTINE thd_ice_alloc *** 
    170148      !!---------------------------------------------------------------------! 
    171       INTEGER :: thd_ice_alloc 
    172       INTEGER :: ierr(4) 
     149      INTEGER ::   thd_ice_alloc   ! return value 
     150      INTEGER ::   ierr(4) 
    173151      !!---------------------------------------------------------------------! 
    174152 
    175       ALLOCATE(npb(jpij)     , npac(jpij),                           & 
    176                ! 
    177                qldif_1d(jpij) , qcmif_1d(jpij) , fstbif_1d(jpij)   , &    
    178                fltbif_1d(jpij), fscbq_1d(jpij) , qsr_ice_1d(jpij)  , &    
    179                fr1_i0_1d(jpij), fr2_i0_1d(jpij), qnsr_ice_1d(jpij) , &     
    180                qfvbq_1d(jpij) , t_bo_b(jpij)   ,                     & 
    181                Stat=ierr(1)) 
    182                ! 
    183       ALLOCATE(sprecip_1d(jpij), frld_1d(jpij)   , at_i_b(jpij)    , &     
    184                fbif_1d(jpij)   , rdmicif_1d(jpij), rdmsnif_1d(jpij), & 
    185                qlbbq_1d(jpij)  , dmgwi_1d(jpij)  , dvsbq_1d(jpij)  , &    
    186                dvbbq_1d(jpij)  , dvlbq_1d(jpij)  , dvnbq_1d(jpij)  , &    
    187                dqns_ice_1d(jpij),qla_ice_1d(jpij), dqla_ice_1d(jpij),& 
    188                tatm_ice_1d(jpij),fsup(jpij)      , focea(jpij)     , &    
    189                i0(jpij)        , old_ht_i_b(jpij), old_ht_s_b(jpij), &   
    190                fsbri_1d(jpij)  , fhbri_1d(jpij)  , fseqv_1d(jpij)  , & 
    191                dsm_i_fl_1d(jpij),dsm_i_gd_1d(jpij),dsm_i_se_1d(jpij),&      
    192                dsm_i_si_1d(jpij),hicol_b(jpij)                     , & 
    193                Stat=ierr(2)) 
    194                ! 
    195       ALLOCATE(t_su_b(jpij)     , a_i_b(jpij)    , ht_i_b(jpij)    , &    
    196                ht_s_b(jpij)     , fc_su(jpij)    , fc_bo_i(jpij)   , &     
    197                dh_s_tot(jpij)   , dh_i_surf(jpij), dh_i_bott(jpij) , &     
    198                dh_snowice(jpij) , sm_i_b(jpij)   , s_i_new(jpij)   , &     
    199                s_snowice(jpij)  , o_i_b(jpij)                      , & 
    200                ! 
    201                t_s_b(jpij,nlay_s),                                   & 
    202                ! 
    203                t_i_b(jpij,jkmax), s_i_b(jpij,jkmax)                , &             
    204                q_i_b(jpij,jkmax), q_s_b(jpij,jkmax)                , & 
    205                Stat=ierr(3)) 
    206                ! 
    207       ALLOCATE(qt_i_in(jpij,jpl) , qt_i_fin(jpij,jpl), qt_s_in(jpij,jpl),   & 
    208                qt_s_fin(jpij,jpl), dq_i(jpij,jpl)    , sum_fluxq(jpij,jpl), & 
    209                fatm(jpij,jpl),     foce(jpij,jpl)    , cons_error(jpij,jpl),& 
    210                surf_error(jpij,jpl),                                        & 
    211                ! 
    212                q_i_layer_in(jpij,jkmax), q_i_layer_fin(jpij,jkmax),        & 
    213                dq_i_layer(jpij,jkmax)  , radab(jpij,jkmax),                & 
    214                ! 
    215                ftotal_in(jpij), ftotal_fin(jpij),                          & 
    216                ! 
    217                fc_s(jpij,0:nlay_s),   fc_i(jpij,0:jkmax)                 , & 
    218                de_s_lay(jpij,nlay_s), de_i_lay(jpij,jkmax)               , & 
    219                ! 
    220                Stat=ierr(4)) 
    221  
    222       thd_ice_alloc = MAXVAL(ierr) 
    223  
    224       IF(thd_ice_alloc /= 0)THEN 
    225          CALL ctl_warn('thd_ice_alloc: failed to allocate arrays.') 
     153      ALLOCATE( npb      (jpij) , npac     (jpij),                          & 
     154         !                                                                  ! 
     155         &      qldif_1d (jpij) , qcmif_1d (jpij) , fstbif_1d  (jpij) ,     & 
     156         &      fltbif_1d(jpij) , fscbq_1d (jpij) , qsr_ice_1d (jpij) ,     & 
     157         &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qnsr_ice_1d(jpij) ,     & 
     158         &      qfvbq_1d (jpij) , t_bo_b   (jpij)                     , STAT=ierr(1) ) 
     159      ! 
     160      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     & 
     161         &      fbif_1d    (jpij) , rdmicif_1d (jpij) , rdmsnif_1d (jpij) ,     & 
     162         &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     & 
     163         &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     & 
     164         &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,     & 
     165         &      tatm_ice_1d(jpij) , fsup       (jpij) , focea      (jpij) ,     &    
     166         &      i0         (jpij) , old_ht_i_b (jpij) , old_ht_s_b (jpij) ,     &   
     167         &      fsbri_1d   (jpij) , fhbri_1d   (jpij) , fseqv_1d   (jpij) ,     & 
     168         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &      
     169         &      dsm_i_si_1d(jpij) , hicol_b    (jpij)                     , STAT=ierr(2) ) 
     170      ! 
     171      ALLOCATE( t_su_b    (jpij) , a_i_b    (jpij) , ht_i_b   (jpij) ,    &    
     172         &      ht_s_b    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
     173         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
     174         &      dh_snowice(jpij) , sm_i_b   (jpij) , s_i_new  (jpij) ,    &     
     175         &      s_snowice (jpij) , o_i_b    (jpij)                   ,    & 
     176         !                                                                ! 
     177         &      t_s_b(jpij,nlay_s),                                       & 
     178         !                                                                ! 
     179         &      t_i_b(jpij,jkmax), s_i_b(jpij,jkmax)                ,     &             
     180         &      q_i_b(jpij,jkmax), q_s_b(jpij,jkmax)                , STAT=ierr(3)) 
     181      ! 
     182      ALLOCATE( qt_i_in   (jpij,jpl) , qt_i_fin(jpij,jpl) , qt_s_in   (jpij,jpl) ,     & 
     183         &      qt_s_fin  (jpij,jpl) , dq_i    (jpij,jpl) , sum_fluxq (jpij,jpl) ,     & 
     184         &      fatm      (jpij,jpl) , foce    (jpij,jpl) , cons_error(jpij,jpl) ,     & 
     185         &      surf_error(jpij,jpl)                                             ,     & 
     186         !                                                                             ! 
     187         &      q_i_layer_in(jpij,jkmax) , q_i_layer_fin(jpij,jkmax)             ,     & 
     188         &      dq_i_layer  (jpij,jkmax) , radab        (jpij,jkmax)             ,     & 
     189         !                                                                             ! 
     190         &      ftotal_in(jpij), ftotal_fin(jpij)                                ,     & 
     191         !                                                                             ! 
     192         &      fc_s(jpij,0:nlay_s) , de_s_lay(jpij,nlay_s)                      ,     & 
     193         &      fc_i(jpij,0:jkmax)  , de_i_lay(jpij,jkmax)                       , STAT=ierr(4) ) 
     194 
     195      thd_ice_alloc = MAXVAL( ierr ) 
     196 
     197      IF( thd_ice_alloc /= 0 ) THEN 
     198         CALL ctl_warn( 'thd_ice_alloc: failed to allocate arrays.' ) 
    226199      END IF 
    227  
     200      ! 
    228201   END FUNCTION thd_ice_alloc 
    229  
     202    
     203   !!====================================================================== 
    230204END MODULE thd_ice 
Note: See TracChangeset for help on using the changeset viewer.