Changeset 8540


Ignore:
Timestamp:
2017-09-19T10:01:20+02:00 (3 years ago)
Author:
timgraham
Message:

Changes required to get GO6 running with LIM3 (updated LIM3 and associated files in SBC to head of v3.6 stable, and changes to field_def and namelist_ref)

Location:
branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM
Files:
32 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/CONFIG/SHARED/field_def.xml

    r8308 r8540  
    373373         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   /> 
    374374         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
     375         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   /> 
    375376         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
    376377         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(<20cm) + open water"      unit="m/day"   /> 
     
    599600      <!-- LIM3 scalar variables --> 
    600601 
    601       <field_group id="SBC_scalar"  grid_ref="scalar" > 
     602      <field_group id="SBC_scalar"  domain_ref="1point" > 
    602603         <!-- available with ln_limdiaout --> 
    603          <field id="ibgvoltot"    long_name="global mean ice volume"                                 unit="km3"        /> 
    604          <field id="sbgvoltot"    long_name="global mean snow volume"                                unit="km3"        /> 
    605          <field id="ibgarea"      long_name="global mean ice area"                                   unit="km2"        /> 
    606          <field id="ibgsaline"    long_name="global mean ice salinity"                               unit="0.001"       /> 
    607          <field id="ibgtemper"    long_name="global mean ice temperature"                            unit="degree_C"       /> 
    608          <field id="ibgheatco"    long_name="global mean ice heat content"                           unit="10^20J"     /> 
    609          <field id="sbgheatco"    long_name="global mean snow heat content"                          unit="10^20J"     /> 
    610          <field id="ibgsaltco"    long_name="global mean ice salt content"                           unit="0.001*km3"   /> 
    611  
    612          <field id="ibgvfx"       long_name="global mean volume flux (emp)"                          unit="m/day"      /> 
    613          <field id="ibgvfxbog"    long_name="global mean volume flux (bottom growth)"                unit="m/day"      /> 
    614          <field id="ibgvfxopw"    long_name="global mean volume flux (open water growth)"            unit="m/day"      /> 
    615          <field id="ibgvfxsni"    long_name="global mean volume flux (snow-ice growth)"              unit="m/day"      /> 
    616          <field id="ibgvfxdyn"    long_name="global mean volume flux (dynamic growth)"               unit="m/day"      /> 
    617          <field id="ibgvfxbom"    long_name="global mean volume flux (bottom melt)"                  unit="m/day"      /> 
    618          <field id="ibgvfxsum"    long_name="global mean volume flux (surface melt)"                 unit="m/day"      /> 
    619          <field id="ibgvfxres"    long_name="global mean volume flux (resultant)"                    unit="m/day"      /> 
    620          <field id="ibgvfxspr"    long_name="global mean volume flux (snow precip)"                  unit="m/day"      /> 
    621          <field id="ibgvfxsnw"    long_name="global mean volume flux (snow melt)"                    unit="m/day"      /> 
    622          <field id="ibgvfxsub"    long_name="global mean volume flux (snow sublimation)"             unit="m/day"      /> 
    623  
    624          <field id="ibgsfx"       long_name="global mean salt flux (total)"                          unit="0.001*m/day" /> 
    625          <field id="ibgsfxbri"    long_name="global mean salt flux (brines)"                         unit="0.001*m/day" /> 
    626          <field id="ibgsfxdyn"    long_name="global mean salt flux (dynamic)"                        unit="0.001*m/day" /> 
    627          <field id="ibgsfxres"    long_name="global mean salt flux (resultant)"                      unit="0.001*m/day" /> 
    628          <field id="ibgsfxbog"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    629          <field id="ibgsfxopw"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    630          <field id="ibgsfxsni"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    631          <field id="ibgsfxbom"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    632          <field id="ibgsfxsum"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    633          <field id="ibgsfxsub"    long_name="global mean salt flux (thermo)"                         unit="0.001*m/day" /> 
    634  
    635          <field id="ibghfxdhc"    long_name="Heat content variation in snow and ice"                 unit="W"          /> 
    636          <field id="ibghfxspr"    long_name="Heat content of snow precip"                            unit="W"          /> 
    637  
    638          <field id="ibghfxthd"    long_name="heat fluxes from ice-ocean exchange during thermo"      unit="W"          /> 
    639          <field id="ibghfxsum"    long_name="heat fluxes causing surface ice melt"                   unit="W"          /> 
    640          <field id="ibghfxbom"    long_name="heat fluxes causing bottom ice melt"                    unit="W"          /> 
    641          <field id="ibghfxbog"    long_name="heat fluxes causing bottom ice growth"                  unit="W"          /> 
    642          <field id="ibghfxdif"    long_name="heat fluxes causing ice temperature change"             unit="W"          /> 
    643          <field id="ibghfxopw"    long_name="heat fluxes causing open water ice formation"           unit="W"          /> 
    644          <field id="ibghfxdyn"    long_name="heat fluxes from ice-ocean exchange during dynamic"     unit="W"          /> 
    645          <field id="ibghfxres"    long_name="heat fluxes from ice-ocean exchange during resultant"   unit="W"          /> 
    646          <field id="ibghfxsub"    long_name="heat fluxes from sublimation"                           unit="W"          /> 
    647          <field id="ibghfxsnw"    long_name="heat fluxes from snow-ocean exchange"                   unit="W"          /> 
    648          <field id="ibghfxout"    long_name="non solar heat fluxes received by the ocean"            unit="W"          /> 
    649          <field id="ibghfxin"     long_name="total heat fluxes at the ice surface"                   unit="W"          /> 
    650  
    651          <field id="ibgfrcvol"    long_name="global mean forcing volume (emp)"                       unit="km3"        /> 
    652          <field id="ibgfrcsfx"    long_name="global mean forcing salt   (sfx)"                       unit="0.001*km3"   /> 
    653          <field id="ibgvolgrm"    long_name="global mean ice growth+melt volume"                     unit="km3"        /> 
     604         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"       /> 
     605         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"       /> 
     606         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"     /> 
     607         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"     /> 
     608         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"   /> 
     609         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"      /> 
     610         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"      /> 
     611  
     612         <field id="ibgvolume"    long_name="drift in ice/snow volume (equivalent ocean volume)"            unit="km3"        /> 
     613         <field id="ibgsaltco"    long_name="drift in ice salt content (equivalent ocean volume)"           unit="pss*km3"    /> 
     614         <field id="ibgheatco"    long_name="drift in ice/snow heat content"                                unit="1e20J"      /> 
     615         <field id="ibgheatfx"    long_name="drift in ice/snow heat flux"                                   unit="W/m2"       /> 
     616 
     617         <field id="ibgvol_tot"    long_name="global mean ice volume"                                       unit="km3"        /> 
     618         <field id="sbgvol_tot"    long_name="global mean snow volume"                                      unit="km3"        /> 
     619         <field id="ibgarea_tot"   long_name="global mean ice area"                                         unit="km2"        /> 
     620         <field id="ibgsalt_tot"   long_name="global mean ice salt content"                                 unit="1e-3*km3"   /> 
     621         <field id="ibgheat_tot"   long_name="global mean ice heat content"                                 unit="1e20J"      /> 
     622         <field id="sbgheat_tot"   long_name="global mean snow heat content"                                unit="1e20J"      /> 
    654623      </field_group> 
    655624   
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r6498 r8540  
    6969   rn_relast      =    0.333       !  ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    7070                                   !     advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
    71    nn_ahi0        =    2           !  horizontal diffusivity computation 
    72                                    !     0: use rn_ahi0_ref 
    73                                    !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
    74                                    !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
    75    rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s)  
    76                                    !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    7771/ 
    7872!------------------------------------------------------------------------------ 
    7973&namicehdf     !   Ice horizontal diffusion 
    8074!------------------------------------------------------------------------------ 
     75   nn_ahi0        =    -1          !  horizontal diffusivity computation 
     76                                   !    -1: no diffusion (bypass limhdf) 
     77                                   !     0: use rn_ahi0_ref 
     78                                   !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
     79                                   !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
     80   rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s) 
     81                                   !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    8182   nn_convfrq     = 5              !  convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 
    8283/ 
     
    9899                                   !     0: k = k0 + beta.S/T (Untersteiner, 1964) 
    99100                                   !     1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 
     101   rn_cdsn     = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 
     102                                   !  Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    100103   nn_monocat  = 0                 !  virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 
    101104                                   !     2: simple piling instead of ridging --- temporary option 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90

    r6486 r8540  
    2525   !!---------------------------------------------------------------------- 
    2626   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    27    !! $Id$ 
     27   !! $Id: dom_ice.F90 5123 2015-03-04 16:06:03Z clem $ 
    2828   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    2929   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r7993 r8540  
    212212   REAL(wp), PUBLIC ::   rn_betas         !: coef. for partitioning of snowfall between leads and sea ice 
    213213   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
     214   REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K] 
    214215   REAL(wp), PUBLIC ::   nn_conv_dif      !: maximal number of iterations for heat diffusion 
    215216   REAL(wp), PUBLIC ::   rn_terr_dif      !: maximal tolerated error (C) for heat diffusion 
     
    243244   ! 
    244245   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin] 
    245    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icethi      !: total ice thickness (for all categories) (diag only) 
    246246   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
    247247   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction 
     
    320320   !                                                                  !  this is an extensive variable that has to be transported 
    321321   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days) 
    322    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i    !: Sea-Ice Age times volume per area (days.m) 
    323322   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days) 
     323 
     324   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume 
    324325 
    325326   !! Variables summed over all categories, or associated to all the ice in a single grid cell 
    326327   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice   !: components of the ice velocity (m/s) 
    327    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v   !: components of the ice-ocean stress (N/m2) 
    328328   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s    !: ice and snow total volume per unit area (m) 
    329329   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i           !: ice total fractional area (ice concentration) 
    330330   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i          !: =1-at_i ; total open water fractional area 
    331331   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s    !: ice and snow total heat content 
    332    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ot_i           !: mean age over all categories 
    333    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i           !: mean ice temperature over all categories 
    334    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bv_i           !: brine volume averaged over all categories 
    335    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories [PSU] 
     332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories 
     333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories 
     334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i        !: mean sea ice salinity averaged over all categories [PSU] 
     335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_i        !: mean ice  thickness over all categories 
     337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories 
     338   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories 
    336339 
    337340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s        !: Snow temperatures [K] 
     
    405408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s]  
    406409   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s]  
     410 
    407411   ! 
    408412   !!---------------------------------------------------------------------- 
    409413   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    410    !! $Id$ 
     414   !! $Id: ice.F90 7814 2017-03-20 16:21:42Z clem $ 
    411415   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    412416   !!---------------------------------------------------------------------- 
     
    415419   FUNCTION ice_alloc() 
    416420      !!----------------------------------------------------------------- 
    417       !!               *** Routine ice_alloc_2 *** 
     421      !!               *** Routine ice_alloc *** 
    418422      !!----------------------------------------------------------------- 
    419423      INTEGER :: ice_alloc 
     
    435439 
    436440      ii = ii + 1 
    437       ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,                        & 
     441      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                        & 
    438442         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        & 
    439443         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        & 
     
    456460         &      v_s  (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) ,     & 
    457461         &      sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i  (jpi,jpj,jpl) ,     & 
    458          &      ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl)                      , STAT=ierr(ii) ) 
    459       ii = ii + 1 
    460       ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     & 
     462         &      oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) ) 
     463      ii = ii + 1 
     464      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,      & 
    461465         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     & 
    462          &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     & 
    463          &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) ) 
     466         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     & 
     467         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     & 
     468         &      om_i (jpi,jpj) , STAT=ierr(ii) ) 
    464469      ii = ii + 1 
    465470      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 
     
    501506 
    502507      ice_alloc = MAXVAL( ierr(:) ) 
    503       IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays.') 
     508      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc: failed to allocate arrays.') 
    504509      ! 
    505510   END FUNCTION ice_alloc 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r6486 r8540  
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    37    !! $Id$ 
     37   !! $Id: limadv.F90 5429 2015-06-16 09:57:07Z smasson $ 
    3838   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3939   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r6498 r8540  
    3737   !!---------------------------------------------------------------------- 
    3838   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    39    !! $Id$ 
     39   !! $Id: limcons.F90 6963 2016-09-30 12:40:04Z clem $ 
    4040   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4141   !!---------------------------------------------------------------------- 
     
    288288#if ! defined key_bdy 
    289289      ! heat flux 
    290       zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  & 
    291          &              * e12t * tmask(:,:,1) * zconv )  
     290      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es  & 
     291      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
     292         &              ) * e12t * tmask(:,:,1) * zconv )  
    292293      ! salt flux 
    293294      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

    r6486 r8540  
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    42    !! $Id$ 
     42   !! $Id: limctl.F90 5043 2015-01-28 16:44:18Z clem $ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r6498 r8540  
    3131 
    3232   PUBLIC   lim_diahsb        ! routine called by ice_step.F90 
    33  
    34    real(wp) ::   frc_sal, frc_vol   ! global forcing trends 
    35    real(wp) ::   bg_grme            ! global ice growth+melt trends 
    36  
     33   PUBLIC   lim_diahsb_init   ! routine called in sbcice_lim.F90 
     34 
     35   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents 
     36   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends 
     37    
    3738   !! * Substitutions 
    3839#  include "vectopt_loop_substitute.h90" 
     
    4041   !!---------------------------------------------------------------------- 
    4142   !! NEMO/OPA 3.4 , NEMO Consortium (2012) 
    42    !! $Id$ 
     43   !! $Id: limdiahsb.F90 6963 2016-09-30 12:40:04Z clem $ 
    4344   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4445   !!---------------------------------------------------------------------- 
     
    4647CONTAINS 
    4748 
    48    SUBROUTINE lim_diahsb 
     49   SUBROUTINE lim_diahsb( kt ) 
    4950      !!--------------------------------------------------------------------------- 
    5051      !!                  ***  ROUTINE lim_diahsb  *** 
     
    5354      !!  
    5455      !!--------------------------------------------------------------------------- 
     56      INTEGER, INTENT(in) :: kt    ! number of iteration 
    5557      !! 
    56       real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 
    57       real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   & 
    58       &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub  
    59       real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 
    60       real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub   
    61       real(wp)   ::   zbg_hfx_dhc, zbg_hfx_spr 
    62       real(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in    
    63       real(wp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw 
    64       real(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme  
    65       real(wp)   ::   z1_area                     !    -     - 
    66       REAL(wp)   ::   ztmp 
     58      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem 
     59      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot   
     60      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem   
    6761      !!--------------------------------------------------------------------------- 
    6862      IF( nn_timing == 1 )   CALL timing_start('lim_diahsb') 
    6963 
    70       IF( numit == nstart ) CALL lim_diahsb_init  
    71  
    72       ! 1/area 
    73       z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 ) 
    74  
    75       rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 
    76       ! ----------------------- ! 
    77       ! 1 -  Content variations ! 
    78       ! ----------------------- ! 
    79       zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice  
    80       zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow 
    81       zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area 
    82       zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )       ! mean salt content 
    83       zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content 
    84  
    85       !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 
    86       !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 
    87  
    88       ! Volume 
    89       ztmp = rswitch * z1_area * r1_rau0 * rday 
    90       zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    91       zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    92       zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    93       zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    94       zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    95       zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    96       zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    97       zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    98       zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    99       zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    100       zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    101  
    102       ! Salt 
    103       zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    104       zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    105       zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    106       zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    107  
    108       zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    109       zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    110       zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    111       zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    112       zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    113       zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    114  
    115       ! Heat budget 
    116       zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J] 
    117       zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 
    118       zbg_hfx_dhc  = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    119       zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    120  
    121       zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    122       zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    123       zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    124       zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    125       zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    126       zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    127       zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    128       zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    129       zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    130       zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    131       zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    132       zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    133      
    134       ! --------------------------------------------- ! 
    135       ! 2 - Trends due to forcing and ice growth/melt ! 
    136       ! --------------------------------------------- ! 
    137       z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 
    138       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 
    139       z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
    140                           &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 
    141                           &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 
    142       ! 
    143       frc_vol  = frc_vol  + z_frc_vol  * rdt_ice 
    144       frc_sal  = frc_sal  + z_frc_sal  * rdt_ice 
    145       bg_grme  = bg_grme  + z_bg_grme  * rdt_ice 
     64      ! ----------------------- ! 
     65      ! 1 -  Contents ! 
     66      ! ----------------------- ! 
     67      zbg_ivol = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3) 
     68      zbg_svol = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3) 
     69      zbg_area = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2) 
     70      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3) 
     71      zbg_item = glob_sum( et_i * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
     72      zbg_stem = glob_sum( et_s * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
    14673       
    147       ! difference 
    148       !frc_vol = zbg_ivo - frc_vol 
    149       !frc_sal = zbg_sal - frc_sal 
    150        
    151       ! ----------------------- ! 
    152       ! 3 - Diagnostics writing ! 
    153       ! ----------------------- ! 
    154       rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) ) 
    155       ! 
    156       IF( iom_use('ibgvoltot') )   & 
    157       CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)          
    158       IF( iom_use('sbgvoltot') )   & 
    159       CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9        )   ! snw volume (km3 equivalent liquid)        
    160       IF( iom_use('ibgarea') )   & 
    161       CALL iom_put( 'ibgarea'   , zbg_are * 1.e-6                          )   ! ice area   (km2) 
    162       IF( iom_use('ibgsaline') )   & 
    163       CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) )   ! ice saline (psu) 
    164       IF( iom_use('ibgtemper') )   & 
    165       CALL iom_put( 'ibgtemper' , rswitch * zbg_tem / MAX( zbg_ivo, epsi06 ) )   ! ice temper (C) 
    166       CALL iom_put( 'ibgheatco' , zbg_ihc                                  )   ! ice heat content (1.e20 J)         
    167       CALL iom_put( 'sbgheatco' , zbg_shc                                  )   ! snw heat content (1.e20 J) 
    168       IF( iom_use('ibgsaltco') )   & 
    169       CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9        )   ! ice salt content (psu*km3 equivalent liquid)         
    170  
    171       CALL iom_put( 'ibgvfx'    , zbg_vfx                                  )   ! volume flux emp (m/day liquid) 
    172       CALL iom_put( 'ibgvfxbog' , zbg_vfx_bog                              )   ! volume flux bottom growth     -(m/day equivalent liquid) 
    173       CALL iom_put( 'ibgvfxopw' , zbg_vfx_opw                              )   ! volume flux open water growth - 
    174       CALL iom_put( 'ibgvfxsni' , zbg_vfx_sni                              )   ! volume flux snow ice growth   - 
    175       CALL iom_put( 'ibgvfxdyn' , zbg_vfx_dyn                              )   ! volume flux dynamic growth    - 
    176       CALL iom_put( 'ibgvfxbom' , zbg_vfx_bom                              )   ! volume flux bottom melt       - 
    177       CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum                              )   ! volume flux surface melt      - 
    178       CALL iom_put( 'ibgvfxres' , zbg_vfx_res                              )   ! volume flux resultant         - 
    179       CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr                              )   ! volume flux from snow precip         - 
    180       CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw                              )   ! volume flux from snow melt         - 
    181       CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub                              )   ! volume flux from sublimation         - 
    182            
    183       CALL iom_put( 'ibgsfx'    , zbg_sfx                                  )   ! salt flux         -(psu*m/day equivalent liquid)        
    184       CALL iom_put( 'ibgsfxbri' , zbg_sfx_bri                              )   ! salt flux brines  -       
    185       CALL iom_put( 'ibgsfxdyn' , zbg_sfx_dyn                              )   ! salt flux dynamic -     
    186       CALL iom_put( 'ibgsfxres' , zbg_sfx_res                              )   ! salt flux result  -     
    187       CALL iom_put( 'ibgsfxbog' , zbg_sfx_bog                              )   ! salt flux bottom growth    
    188       CALL iom_put( 'ibgsfxopw' , zbg_sfx_opw                              )   ! salt flux open water growth - 
    189       CALL iom_put( 'ibgsfxsni' , zbg_sfx_sni                              )   ! salt flux snow ice growth   - 
    190       CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       - 
    191       CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      - 
    192       CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      - 
    193  
    194       CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W] 
    195       CALL iom_put( 'ibghfxspr' , zbg_hfx_spr                              )   ! Heat content of snow precip [W] 
    196  
    197       CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   !  
    198       CALL iom_put( 'ibghfxsub' , zbg_hfx_sub                              )   !  
    199       CALL iom_put( 'ibghfxdyn' , zbg_hfx_dyn                              )   !  
    200       CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   !  
    201       CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   !  
    202       CALL iom_put( 'ibghfxsum' , zbg_hfx_sum                              )   !  
    203       CALL iom_put( 'ibghfxbom' , zbg_hfx_bom                              )   !  
    204       CALL iom_put( 'ibghfxbog' , zbg_hfx_bog                              )   !  
    205       CALL iom_put( 'ibghfxdif' , zbg_hfx_dif                              )   !  
    206       CALL iom_put( 'ibghfxopw' , zbg_hfx_opw                              )   !  
    207       CALL iom_put( 'ibghfxout' , zbg_hfx_out                              )   !  
    208       CALL iom_put( 'ibghfxin'  , zbg_hfx_in                               )   !  
    209  
    210       CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9                          )   ! vol - forcing     (km3 equivalent liquid)  
    211       CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9                          )   ! sal - forcing     (psu*km3 equivalent liquid)    
    212       IF( iom_use('ibgvolgrm') )   & 
    213       CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)          
    214  
     74      ! ---------------------------! 
     75      ! 2 - Trends due to forcing  ! 
     76      ! ---------------------------! 
     77      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean  
     78      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm 
     79      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean 
     80      z_frc_tembot =           glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice) 
     81      z_frc_temtop =           glob_sum( hfx_in (:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean 
     82      ! 
     83      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3 
     84      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3 
     85      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss 
     86      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J 
     87      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J 
     88             
     89      ! ----------------------- ! 
     90      ! 3 -  Content variations ! 
     91      ! ----------------------- ! 
     92      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3)  
     93         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  
     94      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss) 
     95         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 
     96      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J) 
     97      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check) 
     98         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) 
     99 
     100      ! ----------------------- ! 
     101      ! 4 -  Drifts             ! 
     102      ! ----------------------- ! 
     103      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 
     104      zdiff_sal = zdiff_sal - frc_sal 
     105      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 
     106 
     107      ! ----------------------- ! 
     108      ! 5 - Diagnostics writing ! 
     109      ! ----------------------- ! 
     110      ! 
     111      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
     112      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
     113      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J) 
     114      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2) 
     115         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     116 
     117      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
     118      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
     119      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
     120      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
     121      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
     122      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2)  
     123         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     124      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2)  
     125         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     126 
     127      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3) 
     128      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3) 
     129      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2) 
     130      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3) 
     131      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J) 
     132      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J) 
    215133      ! 
    216134      IF( lrst_ice )   CALL lim_diahsb_rst( numit, 'WRITE' ) 
    217135      ! 
    218136      IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb') 
    219 ! 
     137      ! 
    220138   END SUBROUTINE lim_diahsb 
    221139 
     
    233151      !!             - Compute coefficients for conversion 
    234152      !!--------------------------------------------------------------------------- 
    235       INTEGER            ::   jk       ! dummy loop indice 
    236153      INTEGER            ::   ierror   ! local integer 
    237154      !! 
     
    247164         WRITE(numout,*) '~~~~~~~~~~~~' 
    248165      ENDIF 
    249       ! 
     166      !       
     167      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 
     168      IF( ierror > 0 )  THEN 
     169         CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' ) 
     170         RETURN 
     171      ENDIF 
     172 
    250173      CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files 
    251174      ! 
     
    263186     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    264187     ! 
    265      INTEGER ::   id1, id2, id3   ! local integers 
    266188     !!---------------------------------------------------------------------- 
    267189     ! 
    268190     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    269191        IF( ln_rstart ) THEN                   !* Read the restart file 
    270            !id1 = iom_varid( numrir, 'frc_vol'  , ldstop = .TRUE. ) 
    271192           ! 
    272193           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    273            IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp 
    274            IF(lwp) WRITE(numout,*) '~~~~~~~' 
    275            CALL iom_get( numrir, 'frc_vol', frc_vol ) 
    276            CALL iom_get( numrir, 'frc_sal', frc_sal ) 
    277            CALL iom_get( numrir, 'bg_grme', bg_grme ) 
     194           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp 
     195           IF(lwp) WRITE(numout,*) '~~~~~~~' 
     196           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  ) 
     197           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  ) 
     198           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  ) 
     199           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  ) 
     200           CALL iom_get( numrir, 'frc_sal'    , frc_sal     ) 
     201           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 
     202           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 
     203           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 
    278204        ELSE 
    279205           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    280206           IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 
    281207           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    282            frc_vol  = 0._wp                                           
    283            frc_sal  = 0._wp                                                  
    284            bg_grme  = 0._wp                                        
     208           ! set trends to 0 
     209           frc_voltop  = 0._wp                                           
     210           frc_volbot  = 0._wp                                           
     211           frc_temtop  = 0._wp                                                  
     212           frc_tembot  = 0._wp                                                  
     213           frc_sal     = 0._wp                                                  
     214           ! record initial ice volume, salt and temp 
     215           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
     216           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
     217           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
     218            
    285219       ENDIF 
    286220 
     
    288222        !                                   ! ------------------- 
    289223        IF(lwp) WRITE(numout,*) '~~~~~~~' 
    290         IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp 
     224        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp 
    291225        IF(lwp) WRITE(numout,*) '~~~~~~~' 
    292         CALL iom_rstput( kt, nitrst, numriw, 'frc_vol'   , frc_vol     ) 
    293         CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'   , frc_sal     ) 
    294         CALL iom_rstput( kt, nitrst, numriw, 'bg_grme'   , bg_grme     ) 
     226        CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop  ) 
     227        CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot  ) 
     228        CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop  ) 
     229        CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot  ) 
     230        CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'    , frc_sal     ) 
     231        CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 
     232        CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 
     233        CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 
    295234        ! 
    296235     ENDIF 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r6486 r8540  
    4141   !!---------------------------------------------------------------------- 
    4242   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    43    !! $Id$ 
     43   !! $Id: limdyn.F90 7607 2017-01-25 15:37:31Z cetlod $ 
    4444   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4545   !!---------------------------------------------------------------------- 
     
    244244      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    245245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
    246          &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 
    247       INTEGER  ::   ji, jj 
    248       REAL(wp) ::   za00, zd_max 
     246         &                nn_nevp, rn_relast 
    249247      !!------------------------------------------------------------------- 
    250248 
     
    272270         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
    273271         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast 
    274          WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0 
    275          WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref 
    276272      ENDIF 
    277273      ! 
     
    279275      rhoco  = rau0  * rn_cio 
    280276      ! 
    281       !  Diffusion coefficients 
    282       SELECT CASE( nn_ahi0 ) 
    283  
    284       CASE( 0 ) 
    285          ahiu(:,:) = rn_ahi0_ref 
    286          ahiv(:,:) = rn_ahi0_ref 
    287  
    288          IF(lwp) WRITE(numout,*) '' 
    289          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
    290  
    291       CASE( 1 )  
    292  
    293          zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    294          IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
    295           
    296          ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    297                                                         !                    (60° = min latitude for ice cover)   
    298          ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
    299  
    300          IF(lwp) WRITE(numout,*) '' 
    301          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
    302          IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
    303           
    304       CASE( 2 )  
    305  
    306          zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    307          IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    308           
    309          za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    310                                                  !                    (60° = min latitude for ice cover)   
    311          DO jj = 1, jpj 
    312             DO ji = 1, jpi 
    313                ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
    314                ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
    315             END DO 
    316          END DO 
    317          ! 
    318          IF(lwp) WRITE(numout,*) '' 
    319          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
    320          IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
    321           
    322       END SELECT 
    323  
    324277   END SUBROUTINE lim_dyn_init 
    325278 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r7993 r8540  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    41    !! $Id$ 
     41   !! $Id: limhdf.F90 7621 2017-01-31 08:43:47Z cetlod $ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
     
    202202      ! ----------------------- 
    203203 
    204       IF(ln_ctl)   THEN 
    205          DO jk = 1 , isize 
    206             zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk) 
    207             WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter 
    208             CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout ) 
    209          END DO 
    210       ENDIF 
     204 !     IF(ln_ctl)   THEN 
     205 !        DO jk = 1 , isize 
     206 !           zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk) 
     207 !           WRITE(charout,FMT="('lim_hdf  : zconv =',D23.16, ' iter =',I4)") zconv, iter 
     208 !           CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout ) 
     209 !        END DO 
     210  !    ENDIF 
    211211      ! 
    212212      CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 ) 
     
    233233      !!------------------------------------------------------------------- 
    234234      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    235       NAMELIST/namicehdf/ nn_convfrq  
    236       !!------------------------------------------------------------------- 
    237       ! 
    238       IF(lwp) THEN 
    239          WRITE(numout,*) 
    240          WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 
    241          WRITE(numout,*) '~~~~~~~' 
    242       ENDIF 
     235      NAMELIST/namicehdf/  nn_ahi0, rn_ahi0_ref, nn_convfrq  
     236      INTEGER  ::   ji, jj 
     237      REAL(wp) ::   za00, zd_max 
     238      !!------------------------------------------------------------------- 
    243239      ! 
    244240      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 
     
    253249      IF(lwp) THEN                          ! control print 
    254250         WRITE(numout,*) 
    255          WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation ' 
    256          WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
     251         WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 
     252         WRITE(numout,*) '~~~~~~~~~~~' 
     253         WRITE(numout,*) '   horizontal diffusivity calculation                          nn_ahi0      = ', nn_ahi0 
     254         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)                  rn_ahi0_ref  = ', rn_ahi0_ref 
     255         WRITE(numout,*) '   convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
    257256      ENDIF 
     257      ! 
     258      !  Diffusion coefficients 
     259      SELECT CASE( nn_ahi0 ) 
     260 
     261      CASE( -1 ) 
     262         ahiu(:,:) = 0._wp 
     263         ahiv(:,:) = 0._wp 
     264 
     265         IF(lwp) WRITE(numout,*) '' 
     266         IF(lwp) WRITE(numout,*) '   No sea-ice diffusion applied' 
     267 
     268      CASE( 0 ) 
     269         ahiu(:,:) = rn_ahi0_ref 
     270         ahiv(:,:) = rn_ahi0_ref 
     271 
     272         IF(lwp) WRITE(numout,*) '' 
     273         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
     274 
     275      CASE( 1 )  
     276 
     277         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     278         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
     279          
     280         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     281                                                        !                    (60deg = min latitude for ice cover)   
     282         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
     283 
     284         IF(lwp) WRITE(numout,*) '' 
     285         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
     286         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
     287          
     288      CASE( 2 )  
     289 
     290         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     291         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     292          
     293         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     294                                                 !                    (60deg = min latitude for ice cover)   
     295         DO jj = 1, jpj 
     296            DO ji = 1, jpi 
     297               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
     298               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     299            END DO 
     300         END DO 
     301         ! 
     302         IF(lwp) WRITE(numout,*) '' 
     303         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
     304         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
     305          
     306      END SELECT 
    258307      ! 
    259308   END SUBROUTINE lim_hdf_init 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r7993 r8540  
    5151   !!---------------------------------------------------------------------- 
    5252   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
    53    !! $Id$ 
     53   !! $Id: limistate.F90 6696 2016-06-13 14:55:42Z clem $ 
    5454   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5555   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r7993 r8540  
    6161   !!---------------------------------------------------------------------- 
    6262   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    63    !! $Id$ 
     63   !! $Id: limitd_me.F90 7814 2017-03-20 16:21:42Z clem $ 
    6464   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6565   !!---------------------------------------------------------------------- 
     
    259259                  closing_net(ji,jj) = 0._wp 
    260260                  opning     (ji,jj) = 0._wp 
     261                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 
    261262               ELSE 
    262263                  iterate_ridging    = 1 
     
    844845         END DO 
    845846    
    846          strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 
     847         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
    847848                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
    848849         ksmooth = 1 
     
    853854      ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    854855         ! 
    855          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) 
     856         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1) 
    856857         ! 
    857858         ksmooth = 1 
     
    866867         DO jj = 1, jpj 
    867868            DO ji = 1, jpi 
    868                strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 
     869               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0))) 
    869870            END DO 
    870871         END DO 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r6486 r8540  
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    42    !! $Id$ 
     42   !! $Id: limitd_th.F90 5407 2015-06-11 19:13:22Z smasson $ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90

    r6486 r8540  
    2727   !!---------------------------------------------------------------------- 
    2828   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    29    !! $Id$ 
     29   !! $Id: limmsh.F90 5123 2015-03-04 16:06:03Z clem $ 
    3030   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3131   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6487 r8540  
    1010   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    1111   !!            3.5  !  2012-08  (R. Benshila)  AGRIF  
     12   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy) 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
     
    5455   !!---------------------------------------------------------------------- 
    5556   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    56    !! $Id$ 
     57   !! $Id: limrhg.F90 6964 2016-09-30 12:41:39Z clem $ 
    5758   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5859   !!---------------------------------------------------------------------- 
     
    9596      !!                 coriolis terms of the momentum equation 
    9697      !!              3) Solve the momentum equation (iterative procedure) 
    97       !!              4) Prevent high velocities if the ice is thin 
    98       !!              5) Recompute invariants of the strain rate tensor 
     98      !!              4) Recompute invariants of the strain rate tensor 
    9999      !!                 which are inputs of the ITD, store stress 
    100100      !!                 for the next time step 
    101       !!              6) Control prints of residual (convergence) 
     101      !!              5) Control prints of residual (convergence) 
    102102      !!                 and charge ellipse. 
    103103      !!                 The user should make sure that the parameters 
     
    106106      !!                 e.g. in the Canadian Archipelago 
    107107      !! 
     108      !! ** Notes   : Boundary condition for ice is chosen no-slip  
     109      !!              but can be adjusted with param rn_shlat 
     110      !! 
    108111      !! References : Hunke and Dukowicz, JPO97 
    109112      !!              Bouillon et al., Ocean Modelling 2009 
     
    115118      INTEGER ::   jter     ! local integers 
    116119      CHARACTER (len=50) ::   charout 
    117       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    118       REAL(wp) ::   za, zstms          ! local scalars 
    119       REAL(wp) ::   zc1, zc2, zc3      ! ice mass 
    120  
    121       REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling 
    122       REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    123       REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars 
    124       REAL(wp) ::   zu_ice2, zv_ice1              ! 
    125       REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre 
    126       REAL(wp) ::   zdst                          ! shear at the center of the grid point 
    127       REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface 
    128       REAL(wp) ::   sigma1, sigma2                ! internal ice stress 
    129  
    130       REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    131       REAL(wp) ::   zintb, zintn  ! dummy argument 
    132  
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc) 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                            
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points 
     120 
     121      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
     122      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
     123      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
     124      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
     125      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
     126      REAL(wp) ::   zTauO, zTauE, zCor                                       ! temporary scalars 
     127 
     128      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
     129      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
     130      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    144131       
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
    146       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells 
    147       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2  
    148       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    149       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
    150       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
    151                                                               !   ocean surface (ssh_m) if ice is not embedded 
    152                                                               !   ice top surface if ice is embedded    
    153  
    154       REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter 
    155       REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors 
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points 
     135      ! 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     140      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
     141      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses 
     143       
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear 
     145      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components 
     146      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence 
     147      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope: 
     148                                                                             !   ocean surface (ssh_m) if ice is not embedded 
     149                                                                             !   ice top surface if ice is embedded    
     150      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays 
     151      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence 
     152      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice 
     153 
     154      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
     155      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     156      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 
    156157      !!------------------------------------------------------------------- 
    157158 
    158       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    159       CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    160       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    161       CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     159      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     160      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     161      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     162      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     163      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    162164 
    163165#if  defined key_lim2 && ! defined key_lim2_vp 
     
    176178      ! 
    177179      !------------------------------------------------------------------------------! 
    178       ! 1) Ice strength (zpresh)                                ! 
    179       !------------------------------------------------------------------------------! 
    180       ! 
    181       ! Put every vector to 0 
    182       delta_i(:,:) = 0._wp   ; 
    183       zpresh (:,:) = 0._wp   ;   
    184       zpreshc(:,:) = 0._wp 
    185       u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp 
    186       divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
    187       shear_i(:,:) = 0._wp 
    188  
     180      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj)  
     181      !------------------------------------------------------------------------------! 
     182      ! ocean/land mask 
     183      DO jj = 1, jpjm1 
     184         DO ji = 1, jpim1      ! NO vector opt. 
     185            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
     186         END DO 
     187      END DO 
     188      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     189 
     190      ! Lateral boundary conditions on velocity (modify zfmask) 
     191      zwf(:,:) = zfmask(:,:) 
     192      DO jj = 2, jpjm1 
     193         DO ji = fs_2, fs_jpim1   ! vector opt. 
     194            IF( zfmask(ji,jj) == 0._wp ) THEN 
     195               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     196            ENDIF 
     197         END DO 
     198      END DO 
     199      DO jj = 2, jpjm1 
     200         IF( zfmask(1,jj) == 0._wp ) THEN 
     201            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     202         ENDIF 
     203         IF( zfmask(jpi,jj) == 0._wp ) THEN 
     204            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
     205         ENDIF 
     206      END DO 
     207      DO ji = 2, jpim1 
     208         IF( zfmask(ji,1) == 0._wp ) THEN 
     209            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     210         ENDIF 
     211         IF( zfmask(ji,jpj) == 0._wp ) THEN 
     212            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     213         ENDIF 
     214      END DO 
     215      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     216 
     217      !------------------------------------------------------------------------------! 
     218      ! 1) define some variables and initialize arrays 
     219      !------------------------------------------------------------------------------! 
     220      ! ecc2: square of yield ellipse eccenticrity 
     221      ecc2    = rn_ecc * rn_ecc 
     222      z1_ecc2 = 1._wp / ecc2 
     223 
     224      ! Time step for subcycling 
     225      zdtevp   = rdt_ice / REAL( nn_nevp ) 
     226      z1_dtevp = 1._wp / zdtevp 
     227 
     228      ! alpha parameters (Bouillon 2009) 
    189229#if defined key_lim3 
    190       CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points 
    191 #endif 
    192  
    193       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    194          DO ji = 1 , jpi 
     230      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     231#else 
     232      zalph1 = ( 2._wp * telast ) * z1_dtevp 
     233#endif 
     234      zalph2 = zalph1 * z1_ecc2 
     235 
     236      z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     237      z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     238 
     239      ! Initialise stress tensor  
     240      zs1 (:,:) = stress1_i (:,:)  
     241      zs2 (:,:) = stress2_i (:,:) 
     242      zs12(:,:) = stress12_i(:,:) 
     243 
     244      ! Ice strength 
    195245#if defined key_lim3 
    196             zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj) 
    197 #endif 
    198 #if defined key_lim2 
    199             zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 
    200 #endif 
    201             ! zmask = 1 where there is ice or on land 
    202             zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 
     246      CALL lim_itd_me_icestrength( nn_icestr ) 
     247      zpresh(:,:) = tmask(:,:,1) *  strength(:,:) 
     248#else 
     249      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) ) 
     250#endif 
     251 
     252      ! scale factors 
     253      DO jj = k_j1+1, k_jpj-1 
     254         DO ji = fs_2, fs_jpim1 
     255            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
     256            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    203257         END DO 
    204258      END DO 
    205  
    206       ! Ice strength on grid cell corners (zpreshc) 
    207       ! needed for calculation of shear stress  
    208       DO jj = k_j1+1, k_jpj-1 
    209          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    210             zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   & 
    211                &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1) 
    212             zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
    213                &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     & 
    214                &             ) / MAX( zstms, zepsi ) 
    215          END DO 
    216       END DO 
    217       CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
     259             
    218260      ! 
    219261      !------------------------------------------------------------------------------! 
    220262      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    221263      !------------------------------------------------------------------------------! 
    222       ! 
    223       !  Wind stress, coriolis and mass terms on the sides of the squares         
    224       !  zfrld1: lead fraction on U-points                                       
    225       !  zfrld2: lead fraction on V-points                                      
    226       !  zmass1: ice/snow mass on U-points                                     
    227       !  zmass2: ice/snow mass on V-points                                    
    228       !  zcorl1: Coriolis parameter on U-points                              
    229       !  zcorl2: Coriolis parameter on V-points                             
    230       !  (ztagnx,ztagny): wind stress on U/V points                        
    231       !  v_oce1: ocean v component on u points                           
    232       !  u_oce2: ocean u component on v points                          
    233264 
    234265      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
     
    242273         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    243274         ! 
    244          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:) ) * r1_rau0 
     275         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 
    245276         ! 
    246277      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
     
    251282         DO ji = fs_2, fs_jpim1 
    252283 
    253             zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    254             zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    255             zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    256  
    257             zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj) 
    258             zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 
    259             zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  ) 
    260             zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 
    261  
    262             ! Leads area. 
    263             zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 
    264             zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 
    265  
    266             ! Mass, coriolis coeff. and currents 
    267             zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 
    268             zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 
    269             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   & 
    270                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 
    271             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   & 
    272                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 
    273             ! 
    274             ! Ocean has no slip boundary condition 
    275             v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      & 
    276                &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  & 
    277                &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)   
    278  
    279             u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      & 
    280                &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  & 
    281                &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    282  
    283             ! Wind stress at U,V-point 
    284             ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 
    285             ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 
    286  
    287             ! Computation of the velocity field taking into account the ice internal interaction. 
    288             ! Terms that are independent of the velocity field. 
    289  
    290             ! SB On utilise maintenant le gradient de la pente de l'ocean 
    291             ! include it later 
    292  
    293             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
    294             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    295  
    296             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    297             za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 
     284            ! ice fraction at U-V points 
     285            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 
     286            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 
     287 
     288            ! Ice/snow mass at U-V points 
     289            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
     290            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
     291            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     292            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 
     293            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 
     294 
     295            ! Ocean currents at U-V points 
     296            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
     297               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     298             
     299            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
     300               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     301 
     302            ! Coriolis at T points (m*f) 
     303            zmf(ji,jj)      = zm1 * fcor(ji,jj) 
     304 
     305            ! m/dt 
     306            zmU_t(ji,jj)    = zmassU * z1_dtevp 
     307            zmV_t(ji,jj)    = zmassV * z1_dtevp 
     308 
     309            ! Drag ice-atm. 
     310            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     311            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     312 
     313            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     314            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     315            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
     316 
     317            ! masks 
     318            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     319            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     320 
     321            ! switches 
     322            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
     323            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
    298324 
    299325         END DO 
    300326      END DO 
    301  
     327      CALL lbc_lnk( zmf, 'T', 1. ) 
    302328      ! 
    303329      !------------------------------------------------------------------------------! 
     
    305331      !------------------------------------------------------------------------------! 
    306332      ! 
    307       ! Time step for subcycling 
    308       dtevp  = rdt_ice / nn_nevp 
    309 #if defined key_lim3 
    310       dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 
    311 #else 
    312       dtotel = dtevp / ( 2._wp * telast ) 
    313 #endif 
    314       z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
    315       z1_dtevp  = 1._wp / dtevp 
    316       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    317       ecc2 = rn_ecc * rn_ecc 
    318       ecci = 1. / ecc2 
    319  
    320       !-Initialise stress tensor  
    321       zs1 (:,:) = stress1_i (:,:)  
    322       zs2 (:,:) = stress2_i (:,:) 
    323       zs12(:,:) = stress12_i(:,:) 
    324  
    325333      !                                               !----------------------! 
    326334      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    327335         !                                            !----------------------!         
    328          DO jj = k_j1, k_jpj-1 
    329             zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step 
    330             zv_ice(:,jj) = v_ice(:,jj) 
    331          END DO 
    332  
    333          DO jj = k_j1+1, k_jpj-1 
    334             DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask 
    335  
    336                !   
    337                !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 
    338                !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 
    339                !- zds(:,:): shear on northeast corner of grid cells 
    340                ! 
    341                !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,  
    342                !                      there are many repeated calculations.  
    343                !                      Speed could be improved by regrouping terms. For 
    344                !                      the moment, however, the stress is on clarity of coding to avoid 
    345                !                      bugs (Martin, for Miguel). 
    346                ! 
    347                !- ALSO: arrays zdt, zds and delta could  
    348                !  be removed in the future to minimise memory demand. 
    349                ! 
    350                !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 
    351                !              grid cells, exactly as in the B grid case. For simplicity, the indexation on 
    352                !              the corners is the same as in the B grid. 
    353                ! 
    354                ! 
    355                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    356                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    357                   &            ) * r1_e12t(ji,jj) 
    358  
    359                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    360                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    361                   &         ) * r1_e12t(ji,jj) 
    362  
    363                ! 
     336         IF(ln_ctl) THEN   ! Convergence test 
     337            DO jj = k_j1, k_jpj-1 
     338               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     339               zv_ice(:,jj) = v_ice(:,jj) 
     340            END DO 
     341         ENDIF 
     342 
     343         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     344         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     345            DO ji = 1, jpim1 
     346 
     347               ! shear at F points 
    364348               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    365349                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    366                   &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   & 
    367                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    368  
    369  
    370                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    371                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    372                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
    373  
    374                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    375                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    376                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    377             END DO 
    378          END DO 
    379  
    380          CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond. 
    381           
     350                  &         ) * r1_e12f(ji,jj) * zfmask(ji,jj) 
     351 
     352            END DO 
     353         END DO 
     354         CALL lbc_lnk( zds, 'F', 1. ) 
     355 
    382356         DO jj = k_j1+1, k_jpj-1 
    383             DO ji = fs_2, fs_jpim1 
    384  
    385                !- Calculate Delta at centre of grid cells 
    386                zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
    387                   &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
    388                   &            ) * r1_e12t(ji,jj) 
    389  
    390                delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    391                delta_i(ji,jj) = delta + rn_creepl 
    392  
    393                !- Calculate Delta on corners 
    394                zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    395                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    396                   &    ) * r1_e12f(ji,jj) 
    397  
    398                zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    399                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    400                   &    ) * r1_e12f(ji,jj) 
    401  
    402                zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 
    403  
    404                !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 
    405                zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    406                   &          ) * z1_dtotel 
    407                zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    408                   &          ) * z1_dtotel 
    409                !-Calculate stress tensor component zs12 at corners 
    410                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  & 
    411                   &          ) * z1_dtotel  
    412  
    413             END DO 
    414          END DO 
    415  
    416          CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     357            DO ji = 2, jpim1 ! no vector loop 
     358 
     359               ! shear**2 at T points (doc eq. A16) 
     360               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     361                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     362                  &   ) * 0.25_wp * r1_e12t(ji,jj) 
     363               
     364               ! divergence at T points 
     365               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     366                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     367                  &    ) * r1_e12t(ji,jj) 
     368               zdiv2 = zdiv * zdiv 
     369                
     370               ! tension at T points 
     371               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     372                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     373                  &   ) * r1_e12t(ji,jj) 
     374               zdt2 = zdt * zdt 
     375                
     376               ! delta at T points 
     377               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )   
     378 
     379               ! P/delta at T points 
     380               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 
     381                
     382               ! stress at T points 
     383               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
     384               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     385              
     386            END DO 
     387         END DO 
     388         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     389 
     390         DO jj = k_j1, k_jpj-1 
     391            DO ji = 1, jpim1 
     392 
     393               ! P/delta at F points 
     394               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     395                
     396               ! stress at F points 
     397               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     398 
     399            END DO 
     400         END DO 
     401         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    417402  
    418          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
     403         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    419404         DO jj = k_j1+1, k_jpj-1 
    420             DO ji = fs_2, fs_jpim1 
    421                !- contribution of zs1, zs2 and zs12 to zf1 
    422                zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  & 
    423                   &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          & 
    424                   &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  & 
    425                   &                ) * r1_e12u(ji,jj) 
    426                ! contribution of zs1, zs2 and zs12 to zf2 
    427                zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
    428                   &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          & 
    429                   &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  & 
    430                   &               )  * r1_e12v(ji,jj) 
     405            DO ji = fs_2, fs_jpim1                
     406 
     407               ! U points 
     408               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     409                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     410                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     411                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     412                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     413                  &                  ) * r1_e12u(ji,jj) 
     414 
     415               ! V points 
     416               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     417                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     418                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     419                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     420                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     421                  &                  ) * r1_e12v(ji,jj) 
     422 
     423               ! u_ice at V point 
     424               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     425                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     426                
     427               ! v_ice at U point 
     428               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     429                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     430 
    431431            END DO 
    432432         END DO 
    433433         ! 
    434          ! Computation of ice velocity 
    435          ! 
    436          ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 
    437          ! 
    438          IF (MOD(jter,2).eq.0) THEN  
    439  
     434         ! --- Computation of ice velocity --- ! 
     435         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
     436         !  Bouillon et al. 2009 (eq 34-35) => stable 
     437         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
     438             
    440439            DO jj = k_j1+1, k_jpj-1 
    441440               DO ji = fs_2, fs_jpim1 
    442                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    443                   z0           = zmass1(ji,jj) * z1_dtevp 
    444  
    445                   ! SB modif because ocean has no slip boundary condition 
    446                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     & 
    447                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    448                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    449                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    450                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    451                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    452                   zcca         = z0 + za 
    453                   zccb         = zcorl1(ji,jj) 
    454                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     441 
     442                  ! tau_io/(v_oce - v_ice) 
     443                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     444                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     445 
     446                  ! Coriolis at V-points (energy conserving formulation) 
     447                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     448                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     449                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     450 
     451                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     452                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     453                   
     454                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     455                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     456                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     457                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     458                     &           ) * zmaskV(ji,jj) 
    455459               END DO 
    456460            END DO 
    457  
    458             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
     461            CALL lbc_lnk( v_ice, 'V', -1. ) 
     462             
     463#if defined key_agrif && defined key_lim2 
     464            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
     465#endif 
     466#if defined key_bdy 
     467            CALL bdy_ice_lim_dyn( 'V' ) 
     468#endif          
     469 
     470            DO jj = k_j1+1, k_jpj-1 
     471               DO ji = fs_2, fs_jpim1 
     472                                
     473                  ! tau_io/(u_oce - u_ice) 
     474                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     475                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     476 
     477                  ! Coriolis at U-points (energy conserving formulation) 
     478                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     479                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     480                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     481                   
     482                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     483                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     484 
     485                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     486                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     487                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     488                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     489                     &           ) * zmaskU(ji,jj) 
     490               END DO 
     491            END DO 
     492            CALL lbc_lnk( u_ice, 'U', -1. ) 
     493             
    459494#if defined key_agrif && defined key_lim2 
    460495            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    461496#endif 
    462497#if defined key_bdy 
    463          CALL bdy_ice_lim_dyn( 'U' ) 
     498            CALL bdy_ice_lim_dyn( 'U' ) 
    464499#endif          
     500 
     501         ELSE ! odd iterations 
    465502 
    466503            DO jj = k_j1+1, k_jpj-1 
    467504               DO ji = fs_2, fs_jpim1 
    468  
    469                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    470                   z0           = zmass2(ji,jj) * z1_dtevp 
    471                   ! SB modif because ocean has no slip boundary condition 
    472                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    473                      &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    474                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    475                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &  
    476                      &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    477                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    478                   zcca         = z0 + za 
    479                   zccb         = zcorl2(ji,jj) 
    480                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
     505                                
     506                  ! tau_io/(u_oce - u_ice) 
     507                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     508                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     509 
     510                  ! Coriolis at U-points (energy conserving formulation) 
     511                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     512                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     513                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     514                   
     515                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     516                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     517 
     518                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     519                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     520                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     521                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     522                     &           ) * zmaskU(ji,jj) 
    481523               END DO 
    482524            END DO 
    483  
    484             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
     525            CALL lbc_lnk( u_ice, 'U', -1. ) 
     526             
     527#if defined key_agrif && defined key_lim2 
     528            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
     529#endif 
     530#if defined key_bdy 
     531            CALL bdy_ice_lim_dyn( 'U' ) 
     532#endif          
     533 
     534           DO jj = k_j1+1, k_jpj-1 
     535               DO ji = fs_2, fs_jpim1 
     536 
     537                  ! tau_io/(v_oce - v_ice) 
     538                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     539                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     540 
     541                  ! Coriolis at V-points (energy conserving formulation) 
     542                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     543                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     544                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     545 
     546                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     547                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     548                   
     549                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     550                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     551                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     552                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     553                     &           ) * zmaskV(ji,jj) 
     554               END DO 
     555            END DO 
     556            CALL lbc_lnk( v_ice, 'V', -1. ) 
     557             
    485558#if defined key_agrif && defined key_lim2 
    486559            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    487560#endif 
    488561#if defined key_bdy 
    489          CALL bdy_ice_lim_dyn( 'V' ) 
     562            CALL bdy_ice_lim_dyn( 'V' ) 
    490563#endif          
    491564 
    492          ELSE  
     565         ENDIF 
     566          
     567         IF(ln_ctl) THEN   ! Convergence test 
    493568            DO jj = k_j1+1, k_jpj-1 
    494                DO ji = fs_2, fs_jpim1 
    495                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    496                   z0           = zmass2(ji,jj) * z1_dtevp 
    497                   ! SB modif because ocean has no slip boundary condition 
    498                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    499                      &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    500                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)    
    501  
    502                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
    503                      &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    504                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    505                   zcca         = z0 + za 
    506                   zccb         = zcorl2(ji,jj) 
    507                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    508                END DO 
    509             END DO 
    510  
    511             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    512 #if defined key_agrif && defined key_lim2 
    513             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    514 #endif 
    515 #if defined key_bdy 
    516          CALL bdy_ice_lim_dyn( 'V' ) 
    517 #endif          
    518  
    519             DO jj = k_j1+1, k_jpj-1 
    520                DO ji = fs_2, fs_jpim1 
    521                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    522                   z0           = zmass1(ji,jj) * z1_dtevp 
    523                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       & 
    524                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    525                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    526  
    527                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    528                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    529                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    530                   zcca         = z0 + za 
    531                   zccb         = zcorl1(ji,jj) 
    532                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    533                END DO 
    534             END DO 
    535  
    536             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    537 #if defined key_agrif && defined key_lim2 
    538             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    539 #endif 
    540 #if defined key_bdy 
    541          CALL bdy_ice_lim_dyn( 'U' ) 
    542 #endif          
    543  
    544          ENDIF 
    545           
    546          IF(ln_ctl) THEN 
    547             !---  Convergence test. 
    548             DO jj = k_j1+1 , k_jpj-1 
    549569               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    550570            END DO 
     
    552572            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    553573         ENDIF 
    554  
     574         ! 
    555575         !                                                ! ==================== ! 
    556576      END DO                                              !  end loop over jter  ! 
     
    558578      ! 
    559579      !------------------------------------------------------------------------------! 
    560       ! 4) Prevent ice velocities when the ice is thin 
    561       !------------------------------------------------------------------------------! 
    562       ! If the ice volume is below zvmin then ice velocity should equal the 
    563       ! ocean velocity. This prevents high velocity when ice is thin 
    564       DO jj = k_j1+1, k_jpj-1 
    565          DO ji = fs_2, fs_jpim1 
    566             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    567                u_ice(ji,jj) = u_oce(ji,jj) 
    568                v_ice(ji,jj) = v_oce(ji,jj) 
    569             ENDIF 
     580      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
     581      !------------------------------------------------------------------------------! 
     582      DO jj = k_j1, k_jpj-1  
     583         DO ji = 1, jpim1 
     584 
     585            ! shear at F points 
     586            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     587               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     588               &         ) * r1_e12f(ji,jj) * zfmask(ji,jj) 
     589 
     590         END DO 
     591      END DO            
     592      CALL lbc_lnk( zds, 'F', 1. ) 
     593       
     594      DO jj = k_j1+1, k_jpj-1  
     595         DO ji = 2, jpim1 ! no vector loop 
     596             
     597            ! tension**2 at T points 
     598            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     599               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     600               &   ) * r1_e12t(ji,jj) 
     601            zdt2 = zdt * zdt 
     602             
     603            ! shear**2 at T points (doc eq. A16) 
     604            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     605               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     606               &   ) * 0.25_wp * r1_e12t(ji,jj) 
     607             
     608            ! shear at T points 
     609            shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     610 
     611            ! divergence at T points 
     612            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     613               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     614               &            ) * r1_e12t(ji,jj) 
     615             
     616            ! delta at T points 
     617            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
     618            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
     619            delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     620 
    570621         END DO 
    571622      END DO 
    572  
    573       CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
    574  
    575 #if defined key_agrif && defined key_lim2 
    576       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
    577       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 
    578 #endif 
    579 #if defined key_bdy 
    580       CALL bdy_ice_lim_dyn( 'U' ) 
    581       CALL bdy_ice_lim_dyn( 'V' ) 
    582 #endif          
    583  
    584       DO jj = k_j1+1, k_jpj-1  
    585          DO ji = fs_2, fs_jpim1 
    586             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    587                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     & 
    588                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    589                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    590  
    591                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    592                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    593                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    594             ENDIF  
    595          END DO 
    596       END DO 
    597  
    598       CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    599  
    600       ! Recompute delta, shear and div, inputs for mechanical redistribution  
    601       DO jj = k_j1+1, k_jpj-1 
    602          DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask 
    603             !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    604             !- zds(:,:): shear on northeast corner of grid cells 
    605             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    606  
    607                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   & 
    608                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   & 
    609                   &            ) * r1_e12t(ji,jj) 
    610  
    611                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
    612                   &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
    613                   &         ) * r1_e12t(ji,jj) 
    614                ! 
    615                ! SB modif because ocean has no slip boundary condition  
    616                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    617                   &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    618                   &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     & 
    619                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    620  
    621                zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    & 
    622                   &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj) 
    623  
    624                delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    625                delta_i(ji,jj) = delta + rn_creepl 
    626              
    627             ENDIF 
    628          END DO 
    629       END DO 
    630       ! 
    631       !------------------------------------------------------------------------------! 
    632       ! 5) Store stress tensor and its invariants 
    633       !------------------------------------------------------------------------------! 
    634       ! * Invariants of the stress tensor are required for limitd_me 
    635       !   (accelerates convergence and improves stability) 
    636       DO jj = k_j1+1, k_jpj-1 
    637          DO ji = fs_2, fs_jpim1 
    638             zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &    
    639                &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj)  
    640             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    641          END DO 
    642       END DO 
    643  
    644       ! Lateral boundary condition 
    645       CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    646  
    647       ! * Store the stress tensor for the next time step 
     623      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
     624       
     625      ! --- Store the stress tensor for the next time step --- ! 
    648626      stress1_i (:,:) = zs1 (:,:) 
    649627      stress2_i (:,:) = zs2 (:,:) 
     
    652630      ! 
    653631      !------------------------------------------------------------------------------! 
    654       ! 6) Control prints of residual and charge ellipse 
     632      ! 5) Control prints of residual and charge ellipse 
    655633      !------------------------------------------------------------------------------! 
    656634      ! 
     
    675653               DO ji = 2, jpim1 
    676654                  IF (zpresh(ji,jj) > 1.0) THEN 
    677                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    678                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
     655                     zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
     656                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
    679657                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    680658                     CALL prt_ctl_info(charout) 
     
    687665      ENDIF 
    688666      ! 
    689       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    690       CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    691       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    692       CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     667      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     668      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     669      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     670      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     671      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    693672 
    694673   END SUBROUTINE lim_rhg 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r6486 r8540  
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    42    !! $Id$ 
     42   !! $Id: limrst.F90 7814 2017-03-20 16:21:42Z clem $ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
     
    108108      INTEGER ::   iter 
    109109      CHARACTER(len=15) ::   znam 
    110       CHARACTER(len=1)  ::   zchar, zchar1 
     110      CHARACTER(len=2)  ::   zchar, zchar1 
    111111      REAL(wp), POINTER, DIMENSION(:,:) :: z2d 
    112112      !!---------------------------------------------------------------------- 
     
    130130      ! Prognostic variables  
    131131      DO jl = 1, jpl  
    132          WRITE(zchar,'(I1)') jl 
    133          znam = 'v_i'//'_htc'//zchar 
     132         WRITE(zchar,'(I2)') jl 
     133         znam = 'v_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    134134         z2d(:,:) = v_i(:,:,jl) 
    135135         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    136          znam = 'v_s'//'_htc'//zchar 
     136         znam = 'v_s'//'_htc'//TRIM(ADJUSTL(zchar)) 
    137137         z2d(:,:) = v_s(:,:,jl) 
    138138         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    139          znam = 'smv_i'//'_htc'//zchar 
     139         znam = 'smv_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    140140         z2d(:,:) = smv_i(:,:,jl) 
    141141         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    142          znam = 'oa_i'//'_htc'//zchar 
     142         znam = 'oa_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    143143         z2d(:,:) = oa_i(:,:,jl) 
    144144         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    145          znam = 'a_i'//'_htc'//zchar 
     145         znam = 'a_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    146146         z2d(:,:) = a_i(:,:,jl) 
    147147         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    148          znam = 't_su'//'_htc'//zchar 
     148         znam = 't_su'//'_htc'//TRIM(ADJUSTL(zchar)) 
    149149         z2d(:,:) = t_su(:,:,jl) 
    150150         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    151       END DO 
    152  
    153       DO jl = 1, jpl  
    154          WRITE(zchar,'(I1)') jl 
    155          znam = 'tempt_sl1'//'_htc'//zchar 
     151         znam = 'tempt_sl1'//'_htc'//TRIM(ADJUSTL(zchar)) 
    156152         z2d(:,:) = e_s(:,:,1,jl) 
    157153         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    158       END DO 
    159  
    160       DO jl = 1, jpl  
    161          WRITE(zchar,'(I1)') jl 
    162154         DO jk = 1, nlay_i  
    163             WRITE(zchar1,'(I1)') jk 
    164             znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 
     155            WRITE(zchar1,'(I2)') jk 
     156            znam = 'tempt'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    165157            z2d(:,:) = e_i(:,:,jk,jl) 
    166158            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     
    177169 
    178170      DO jl = 1, jpl  
    179          WRITE(zchar,'(I1)') jl 
    180          znam = 'sxice'//'_htc'//zchar 
     171         WRITE(zchar,'(I2)') jl 
     172         znam = 'sxice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    181173         z2d(:,:) = sxice(:,:,jl) 
    182174         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    183          znam = 'syice'//'_htc'//zchar 
     175         znam = 'syice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    184176         z2d(:,:) = syice(:,:,jl) 
    185177         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    186          znam = 'sxxice'//'_htc'//zchar 
     178         znam = 'sxxice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    187179         z2d(:,:) = sxxice(:,:,jl) 
    188180         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    189          znam = 'syyice'//'_htc'//zchar 
     181         znam = 'syyice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    190182         z2d(:,:) = syyice(:,:,jl) 
    191183         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    192          znam = 'sxyice'//'_htc'//zchar 
     184         znam = 'sxyice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    193185         z2d(:,:) = sxyice(:,:,jl) 
    194186         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    195          znam = 'sxsn'//'_htc'//zchar 
     187         znam = 'sxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    196188         z2d(:,:) = sxsn(:,:,jl) 
    197189         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    198          znam = 'sysn'//'_htc'//zchar 
     190         znam = 'sysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    199191         z2d(:,:) = sysn(:,:,jl) 
    200192         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    201          znam = 'sxxsn'//'_htc'//zchar 
     193         znam = 'sxxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    202194         z2d(:,:) = sxxsn(:,:,jl) 
    203195         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    204          znam = 'syysn'//'_htc'//zchar 
     196         znam = 'syysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    205197         z2d(:,:) = syysn(:,:,jl) 
    206198         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    207          znam = 'sxysn'//'_htc'//zchar 
     199         znam = 'sxysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    208200         z2d(:,:) = sxysn(:,:,jl) 
    209201         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    210          znam = 'sxa'//'_htc'//zchar 
     202         znam = 'sxa'//'_htc'//TRIM(ADJUSTL(zchar)) 
    211203         z2d(:,:) = sxa(:,:,jl) 
    212204         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    213          znam = 'sya'//'_htc'//zchar 
     205         znam = 'sya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    214206         z2d(:,:) = sya(:,:,jl) 
    215207         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    216          znam = 'sxxa'//'_htc'//zchar 
     208         znam = 'sxxa'//'_htc'//TRIM(ADJUSTL(zchar)) 
    217209         z2d(:,:) = sxxa(:,:,jl) 
    218210         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    219          znam = 'syya'//'_htc'//zchar 
     211         znam = 'syya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    220212         z2d(:,:) = syya(:,:,jl) 
    221213         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    222          znam = 'sxya'//'_htc'//zchar 
     214         znam = 'sxya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    223215         z2d(:,:) = sxya(:,:,jl) 
    224216         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    225          znam = 'sxc0'//'_htc'//zchar 
     217         znam = 'sxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    226218         z2d(:,:) = sxc0(:,:,jl) 
    227219         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    228          znam = 'syc0'//'_htc'//zchar 
     220         znam = 'syc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    229221         z2d(:,:) = syc0(:,:,jl) 
    230222         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    231          znam = 'sxxc0'//'_htc'//zchar 
     223         znam = 'sxxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    232224         z2d(:,:) = sxxc0(:,:,jl) 
    233225         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    234          znam = 'syyc0'//'_htc'//zchar 
     226         znam = 'syyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    235227         z2d(:,:) = syyc0(:,:,jl) 
    236228         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    237          znam = 'sxyc0'//'_htc'//zchar 
     229         znam = 'sxyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    238230         z2d(:,:) = sxyc0(:,:,jl) 
    239231         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    240          znam = 'sxsal'//'_htc'//zchar 
     232         znam = 'sxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    241233         z2d(:,:) = sxsal(:,:,jl) 
    242234         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    243          znam = 'sysal'//'_htc'//zchar 
     235         znam = 'sysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    244236         z2d(:,:) = sysal(:,:,jl) 
    245237         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    246          znam = 'sxxsal'//'_htc'//zchar 
     238         znam = 'sxxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    247239         z2d(:,:) = sxxsal(:,:,jl) 
    248240         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    249          znam = 'syysal'//'_htc'//zchar 
     241         znam = 'syysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    250242         z2d(:,:) = syysal(:,:,jl) 
    251243         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    252          znam = 'sxysal'//'_htc'//zchar 
     244         znam = 'sxysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    253245         z2d(:,:) = sxysal(:,:,jl) 
    254246         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    255          znam = 'sxage'//'_htc'//zchar 
     247         znam = 'sxage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    256248         z2d(:,:) = sxage(:,:,jl) 
    257249         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    258          znam = 'syage'//'_htc'//zchar 
     250         znam = 'syage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    259251         z2d(:,:) = syage(:,:,jl) 
    260252         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    261          znam = 'sxxage'//'_htc'//zchar 
     253         znam = 'sxxage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    262254         z2d(:,:) = sxxage(:,:,jl) 
    263255         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    264          znam = 'syyage'//'_htc'//zchar 
     256         znam = 'syyage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    265257         z2d(:,:) = syyage(:,:,jl) 
    266258         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    267          znam = 'sxyage'//'_htc'//zchar 
     259         znam = 'sxyage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    268260         z2d(:,:) = sxyage(:,:,jl) 
    269261         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     
    277269 
    278270      DO jl = 1, jpl  
    279          WRITE(zchar,'(I1)') jl 
     271         WRITE(zchar,'(I2)') jl 
    280272         DO jk = 1, nlay_i  
    281             WRITE(zchar1,'(I1)') jk 
    282             znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     273            WRITE(zchar1,'(I2)') jk 
     274            znam = 'sxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    283275            z2d(:,:) = sxe(:,:,jk,jl) 
    284276            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    285             znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     277            znam = 'sye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    286278            z2d(:,:) = sye(:,:,jk,jl) 
    287279            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    288             znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     280            znam = 'sxxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    289281            z2d(:,:) = sxxe(:,:,jk,jl) 
    290282            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    291             znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     283            znam = 'syye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    292284            z2d(:,:) = syye(:,:,jk,jl) 
    293285            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    294             znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     286            znam = 'sxye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    295287            z2d(:,:) = sxye(:,:,jk,jl) 
    296288            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
     
    318310      REAL(wp), POINTER, DIMENSION(:,:) ::   z2d 
    319311      CHARACTER(len=15) ::   znam 
    320       CHARACTER(len=1)  ::   zchar, zchar1 
     312      CHARACTER(len=2)  ::   zchar, zchar1 
    321313      INTEGER           ::   jlibalt = jprstlib 
    322314      LOGICAL           ::   llok 
     
    356348         &                   '   control of time parameter  nrstdt' ) 
    357349 
     350      ! Prognostic variables  
    358351      DO jl = 1, jpl  
    359          WRITE(zchar,'(I1)') jl 
    360          znam = 'v_i'//'_htc'//zchar 
     352         WRITE(zchar,'(I2)') jl 
     353         znam = 'v_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    361354         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    362355         v_i(:,:,jl) = z2d(:,:) 
    363          znam = 'v_s'//'_htc'//zchar 
     356         znam = 'v_s'//'_htc'//TRIM(ADJUSTL(zchar)) 
    364357         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    365358         v_s(:,:,jl) = z2d(:,:)  
    366          znam = 'smv_i'//'_htc'//zchar 
     359         znam = 'smv_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    367360         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    368361         smv_i(:,:,jl) = z2d(:,:) 
    369          znam = 'oa_i'//'_htc'//zchar 
     362         znam = 'oa_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    370363         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    371364         oa_i(:,:,jl) = z2d(:,:) 
    372          znam = 'a_i'//'_htc'//zchar 
     365         znam = 'a_i'//'_htc'//TRIM(ADJUSTL(zchar)) 
    373366         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    374367         a_i(:,:,jl) = z2d(:,:) 
    375          znam = 't_su'//'_htc'//zchar 
     368         znam = 't_su'//'_htc'//TRIM(ADJUSTL(zchar)) 
    376369         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    377370         t_su(:,:,jl) = z2d(:,:) 
    378       END DO 
    379  
    380       DO jl = 1, jpl  
    381          WRITE(zchar,'(I1)') jl 
    382          znam = 'tempt_sl1'//'_htc'//zchar 
     371         znam = 'tempt_sl1'//'_htc'//TRIM(ADJUSTL(zchar)) 
    383372         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    384373         e_s(:,:,1,jl) = z2d(:,:) 
    385       END DO 
    386  
    387       DO jl = 1, jpl  
    388          WRITE(zchar,'(I1)') jl 
    389374         DO jk = 1, nlay_i  
    390             WRITE(zchar1,'(I1)') jk 
    391             znam = 'tempt'//'_il'//zchar1//'_htc'//zchar 
     375            WRITE(zchar1,'(I2)') jk 
     376            znam = 'tempt'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    392377            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    393378            e_i(:,:,jk,jl) = z2d(:,:) 
     
    404389 
    405390      DO jl = 1, jpl  
    406          WRITE(zchar,'(I1)') jl 
    407          znam = 'sxice'//'_htc'//zchar 
     391         WRITE(zchar,'(I2)') jl 
     392         znam = 'sxice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    408393         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    409394         sxice(:,:,jl) = z2d(:,:) 
    410          znam = 'syice'//'_htc'//zchar 
     395         znam = 'syice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    411396         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    412397         syice(:,:,jl) = z2d(:,:) 
    413          znam = 'sxxice'//'_htc'//zchar 
     398         znam = 'sxxice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    414399         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    415400         sxxice(:,:,jl) = z2d(:,:) 
    416          znam = 'syyice'//'_htc'//zchar 
     401         znam = 'syyice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    417402         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    418403         syyice(:,:,jl) = z2d(:,:) 
    419          znam = 'sxyice'//'_htc'//zchar 
     404         znam = 'sxyice'//'_htc'//TRIM(ADJUSTL(zchar)) 
    420405         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    421406         sxyice(:,:,jl) = z2d(:,:) 
    422          znam = 'sxsn'//'_htc'//zchar 
     407         znam = 'sxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    423408         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    424409         sxsn(:,:,jl) = z2d(:,:) 
    425          znam = 'sysn'//'_htc'//zchar 
     410         znam = 'sysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    426411         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    427412         sysn(:,:,jl) = z2d(:,:) 
    428          znam = 'sxxsn'//'_htc'//zchar 
     413         znam = 'sxxsn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    429414         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    430415         sxxsn(:,:,jl) = z2d(:,:) 
    431          znam = 'syysn'//'_htc'//zchar 
     416         znam = 'syysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    432417         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    433418         syysn(:,:,jl) = z2d(:,:) 
    434          znam = 'sxysn'//'_htc'//zchar 
     419         znam = 'sxysn'//'_htc'//TRIM(ADJUSTL(zchar)) 
    435420         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    436421         sxysn(:,:,jl) = z2d(:,:) 
    437          znam = 'sxa'//'_htc'//zchar 
     422         znam = 'sxa'//'_htc'//TRIM(ADJUSTL(zchar)) 
    438423         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    439424         sxa(:,:,jl) = z2d(:,:) 
    440          znam = 'sya'//'_htc'//zchar 
     425         znam = 'sya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    441426         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    442427         sya(:,:,jl) = z2d(:,:) 
    443          znam = 'sxxa'//'_htc'//zchar 
     428         znam = 'sxxa'//'_htc'//TRIM(ADJUSTL(zchar)) 
    444429         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    445430         sxxa(:,:,jl) = z2d(:,:) 
    446          znam = 'syya'//'_htc'//zchar 
     431         znam = 'syya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    447432         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    448433         syya(:,:,jl) = z2d(:,:) 
    449          znam = 'sxya'//'_htc'//zchar 
     434         znam = 'sxya'//'_htc'//TRIM(ADJUSTL(zchar)) 
    450435         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    451436         sxya(:,:,jl) = z2d(:,:) 
    452          znam = 'sxc0'//'_htc'//zchar 
     437         znam = 'sxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    453438         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    454439         sxc0(:,:,jl) = z2d(:,:) 
    455          znam = 'syc0'//'_htc'//zchar 
     440         znam = 'syc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    456441         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    457442         syc0(:,:,jl) = z2d(:,:) 
    458          znam = 'sxxc0'//'_htc'//zchar 
     443         znam = 'sxxc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    459444         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    460445         sxxc0(:,:,jl) = z2d(:,:) 
    461          znam = 'syyc0'//'_htc'//zchar 
     446         znam = 'syyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    462447         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    463448         syyc0(:,:,jl) = z2d(:,:) 
    464          znam = 'sxyc0'//'_htc'//zchar 
     449         znam = 'sxyc0'//'_htc'//TRIM(ADJUSTL(zchar)) 
    465450         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    466451         sxyc0(:,:,jl) = z2d(:,:) 
    467          znam = 'sxsal'//'_htc'//zchar 
     452         znam = 'sxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    468453         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    469454         sxsal(:,:,jl) = z2d(:,:) 
    470          znam = 'sysal'//'_htc'//zchar 
     455         znam = 'sysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    471456         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    472457         sysal(:,:,jl) = z2d(:,:) 
    473          znam = 'sxxsal'//'_htc'//zchar 
     458         znam = 'sxxsal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    474459         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    475460         sxxsal(:,:,jl) = z2d(:,:) 
    476          znam = 'syysal'//'_htc'//zchar 
     461         znam = 'syysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    477462         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    478463         syysal(:,:,jl) = z2d(:,:) 
    479          znam = 'sxysal'//'_htc'//zchar 
     464         znam = 'sxysal'//'_htc'//TRIM(ADJUSTL(zchar)) 
    480465         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    481466         sxysal(:,:,jl) = z2d(:,:) 
    482          znam = 'sxage'//'_htc'//zchar 
     467         znam = 'sxage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    483468         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    484469         sxage(:,:,jl) = z2d(:,:) 
    485          znam = 'syage'//'_htc'//zchar 
     470         znam = 'syage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    486471         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    487472         syage(:,:,jl) = z2d(:,:) 
    488          znam = 'sxxage'//'_htc'//zchar 
     473         znam = 'sxxage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    489474         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    490475         sxxage(:,:,jl) = z2d(:,:) 
    491          znam = 'syyage'//'_htc'//zchar 
     476         znam = 'syyage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    492477         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    493478         syyage(:,:,jl) = z2d(:,:) 
    494          znam = 'sxyage'//'_htc'//zchar 
     479         znam = 'sxyage'//'_htc'//TRIM(ADJUSTL(zchar)) 
    495480         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    496481         sxyage(:,:,jl)= z2d(:,:) 
     
    504489 
    505490      DO jl = 1, jpl  
    506          WRITE(zchar,'(I1)') jl 
     491         WRITE(zchar,'(I2)') jl 
    507492         DO jk = 1, nlay_i  
    508             WRITE(zchar1,'(I1)') jk 
    509             znam = 'sxe'//'_il'//zchar1//'_htc'//zchar 
     493            WRITE(zchar1,'(I2)') jk 
     494            znam = 'sxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    510495            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    511496            sxe(:,:,jk,jl) = z2d(:,:) 
    512             znam = 'sye'//'_il'//zchar1//'_htc'//zchar 
     497            znam = 'sye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    513498            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    514499            sye(:,:,jk,jl) = z2d(:,:) 
    515             znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar 
     500            znam = 'sxxe'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    516501            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    517502            sxxe(:,:,jk,jl) = z2d(:,:) 
    518             znam = 'syye'//'_il'//zchar1//'_htc'//zchar 
     503            znam = 'syye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    519504            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    520505            syye(:,:,jk,jl) = z2d(:,:) 
    521             znam = 'sxye'//'_il'//zchar1//'_htc'//zchar 
     506            znam = 'sxye'//'_il'//TRIM(ADJUSTL(zchar1))//'_htc'//TRIM(ADJUSTL(zchar)) 
    522507            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    523508            sxye(:,:,jk,jl) = z2d(:,:) 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r6498 r8540  
    6060   !!---------------------------------------------------------------------- 
    6161   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    62    !! $Id$ 
     62   !! $Id: limsbc.F90 6963 2016-09-30 12:40:04Z clem $ 
    6363   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6464   !!---------------------------------------------------------------------- 
     
    110110      !!--------------------------------------------------------------------- 
    111111 
    112       ! make calls for heat fluxes before it is modified 
    113       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
    114       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
    115       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
    116       IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
    117       IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
    118       IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
    119       IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
    120       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
    121          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
    122       IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
    123       IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
    124       IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice) 
    125       IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
    126  
    127       ! albedo output 
     112      ! make call for albedo output before it is modified 
    128113      CALL wrk_alloc( jpi,jpj, zalb )     
    129114 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90

    r6486 r8540  
    2121   !!---------------------------------------------------------------------- 
    2222   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    23    !! $Id$ 
     23   !! $Id: limtab.F90 4161 2013-11-07 10:01:27Z cetlod $ 
    2424   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    2525   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r6498 r8540  
    5656   !!---------------------------------------------------------------------- 
    5757   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    58    !! $Id$ 
     58   !! $Id: limthd.F90 7607 2017-01-25 15:37:31Z cetlod $ 
    5959   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6060   !!---------------------------------------------------------------------- 
     
    613613         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    614614         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
    615          !          
    616615      END SELECT 
    617616 
     
    633632      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    634633      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
    635          &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
    636          &                nn_monocat, ln_it_qnsice 
     634         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon,         & 
     635         &                rn_cdsn, nn_monocat, ln_it_qnsice 
    637636      !!------------------------------------------------------------------- 
    638637      ! 
     
    673672         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif 
    674673         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon 
     674         WRITE(numout,*)'      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn 
    675675         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i 
    676676         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r7993 r8540  
    3838   !!---------------------------------------------------------------------- 
    3939   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    40    !! $Id$ 
     40   !! $Id: limthd_dh.F90 6469 2016-04-13 15:15:13Z clem $ 
    4141   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4242   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r6486 r8540  
    3232   !!---------------------------------------------------------------------- 
    3333   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    34    !! $Id$ 
     34   !! $Id: limthd_dif.F90 7607 2017-01-25 15:37:31Z cetlod $ 
    3535   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3636   !!---------------------------------------------------------------------- 
     
    376376 
    377377               ! Effective thickness he (zhe) 
    378                zfac     = 1._wp / ( rcdsn + zkimean ) 
    379                zratio_s = rcdsn   * zfac 
     378               zfac     = 1._wp / ( rn_cdsn + zkimean ) 
     379               zratio_s = rn_cdsn   * zfac 
    380380               zratio_i = zkimean * zfac 
    381381               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     
    400400         DO ji = kideb, kiut 
    401401            zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
    402             zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac 
    403             zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac 
     402            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac 
     403            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac 
    404404         END DO 
    405405 
    406406         DO jk = 1, nlay_s-1 
    407407            DO ji = kideb , kiut 
    408                zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
     408               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
    409409            END DO 
    410410         END DO 
     
    422422            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
    423423            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
    424             zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
    425            &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
     424            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / &  
     425           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) ) 
    426426            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    427427         END DO 
     
    697697               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    698698         END DO 
     699 
    699700         ! 
    700701         !-------------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r6486 r8540  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    41    !! $Id$ 
     41   !! $Id: limthd_ent.F90 5134 2015-03-09 17:27:34Z clem $ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r6498 r8540  
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    42    !! $Id$ 
     42   !! $Id: limthd_lac.F90 6316 2016-02-15 13:35:37Z cetlod $ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r7993 r8540  
    3333   !!---------------------------------------------------------------------- 
    3434   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    35    !! $Id$ 
     35   !! $Id: limthd_sal.F90 6469 2016-04-13 15:15:13Z clem $ 
    3636   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3737   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r7993 r8540  
    4444   !!---------------------------------------------------------------------- 
    4545   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    46    !! $Id$ 
     46   !! $Id: limtrp.F90 6476 2016-04-15 12:50:29Z mcastril $ 
    4747   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4848   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r6498 r8540  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    41    !! $Id$ 
     41   !! $Id: limupdate1.F90 7814 2017-03-20 16:21:42Z clem $ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
     
    6262 
    6363      IF( kt == nit000 .AND. lwp ) THEN 
    64          WRITE(numout,*) ' lim_update1 '  
    65          WRITE(numout,*) ' ~~~~~~~~~~~ ' 
     64         WRITE(numout,*)''  
     65         WRITE(numout,*)' lim_update1 '  
     66         WRITE(numout,*)' ~~~~~~~~~~~ ' 
    6667      ENDIF 
    6768 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r6498 r8540  
    4141   !!---------------------------------------------------------------------- 
    4242   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    43    !! $Id$ 
     43   !! $Id: limupdate2.F90 7814 2017-03-20 16:21:42Z clem $ 
    4444   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4545   !!---------------------------------------------------------------------- 
     
    6262 
    6363      IF( kt == nit000 .AND. lwp ) THEN 
    64          WRITE(numout,*) ' lim_update2 ' 
    65          WRITE(numout,*) ' ~~~~~~~~~~~ ' 
     64         WRITE(numout,*)'' 
     65         WRITE(numout,*)' lim_update2 ' 
     66         WRITE(numout,*)' ~~~~~~~~~~~ ' 
    6667      ENDIF 
    6768 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r7993 r8540  
    2727   !!                        - et_i(jpi,jpj)  !total ice thermal content  
    2828   !!                        - smt_i(jpi,jpj) !mean ice salinity 
    29    !!                        - ot_i(jpi,jpj)  !average ice age 
     29   !!                        - tm_i (jpi,jpj) !mean ice temperature 
    3030   !!====================================================================== 
    3131   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code 
     
    5454   PUBLIC   lim_var_eqv2glo       
    5555   PUBLIC   lim_var_salprof       
    56    PUBLIC   lim_var_icetm         
    5756   PUBLIC   lim_var_bv            
    5857   PUBLIC   lim_var_salprof1d     
     
    6261   !!---------------------------------------------------------------------- 
    6362   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 
    64    !! $Id$ 
     63   !! $Id: limvar.F90 7814 2017-03-20 16:21:42Z clem $ 
    6564   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6665   !!---------------------------------------------------------------------- 
     
    8988      ! Compute variables 
    9089      !-------------------- 
    91       vt_i (:,:) = 0._wp 
    92       vt_s (:,:) = 0._wp 
    93       at_i (:,:) = 0._wp 
    94       ato_i(:,:) = 1._wp 
    95       ! 
    96       DO jl = 1, jpl 
    97          DO jj = 1, jpj 
    98             DO ji = 1, jpi 
    99                ! 
    100                vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 
    101                vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 
    102                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    103                ! 
    104                rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
    105                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch  ! ice thickness 
    106             END DO 
    107          END DO 
    108       END DO 
    109  
     90      ! integrated values 
     91      vt_i (:,:) = SUM( v_i, dim=3 ) 
     92      vt_s (:,:) = SUM( v_s, dim=3 ) 
     93      at_i (:,:) = SUM( a_i, dim=3 ) 
     94      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
     95      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
     96      ! 
    11097      DO jj = 1, jpj 
    11198         DO ji = 1, jpi 
     
    115102 
    116103      IF( kn > 1 ) THEN 
    117          et_s (:,:) = 0._wp 
    118          ot_i (:,:) = 0._wp 
     104         ! 
     105         ! mean ice/snow thickness 
     106         DO jj = 1, jpj 
     107            DO ji = 1, jpi 
     108               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
     109               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 
     110               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 
     111            ENDDO 
     112         ENDDO 
     113 
     114         ! mean temperature (K), salinity and age 
    119115         smt_i(:,:) = 0._wp 
    120          et_i (:,:) = 0._wp 
    121          ! 
     116         tm_i(:,:)  = 0._wp 
     117         tm_su(:,:) = 0._wp 
     118         om_i (:,:) = 0._wp 
    122119         DO jl = 1, jpl 
     120             
    123121            DO jj = 1, jpj 
    124122               DO ji = 1, jpi 
    125                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
    126                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) )  
    127                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity 
    128                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) )  
    129                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age 
    130                END DO 
    131             END DO 
    132          END DO 
    133          ! 
    134          DO jl = 1, jpl 
     123                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     124                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 
     125                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 ) 
     126               END DO 
     127            END DO 
     128             
    135129            DO jk = 1, nlay_i 
    136                et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content 
    137             END DO 
    138          END DO 
     130               DO jj = 1, jpj 
     131                  DO ji = 1, jpi 
     132                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     133                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     134                        &            / MAX( vt_i(ji,jj) , epsi10 ) 
     135                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  & 
     136                        &            / MAX( vt_i(ji,jj) , epsi10 ) 
     137                  END DO 
     138               END DO 
     139            END DO 
     140         END DO 
     141         tm_i  = tm_i  + rt0 
     142         tm_su = tm_su + rt0 
    139143         ! 
    140144      ENDIF 
     
    246250      ! Mean temperature 
    247251      !------------------- 
    248       vt_i (:,:) = 0._wp 
    249       DO jl = 1, jpl 
    250          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    251       END DO 
     252      ! integrated values 
     253      vt_i (:,:) = SUM( v_i, dim=3 ) 
     254      vt_s (:,:) = SUM( v_s, dim=3 ) 
     255      at_i (:,:) = SUM( a_i, dim=3 ) 
    252256 
    253257      tm_i(:,:) = 0._wp 
     
    398402 
    399403 
    400    SUBROUTINE lim_var_icetm 
    401       !!------------------------------------------------------------------ 
    402       !!                ***  ROUTINE lim_var_icetm *** 
    403       !! 
    404       !! ** Purpose :   computes mean sea ice temperature 
     404   SUBROUTINE lim_var_bv 
     405      !!------------------------------------------------------------------ 
     406      !!                ***  ROUTINE lim_var_bv *** 
     407      !! 
     408      !! ** Purpose :   computes mean brine volume (%) in sea ice 
     409      !! 
     410      !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
     411      !! 
     412      !! References : Vancoppenolle et al., JGR, 2007 
    405413      !!------------------------------------------------------------------ 
    406414      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    407415      !!------------------------------------------------------------------ 
    408  
    409       ! Mean sea ice temperature 
    410       vt_i (:,:) = 0._wp 
    411       DO jl = 1, jpl 
    412          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    413       END DO 
    414  
    415       tm_i(:,:) = 0._wp 
     416      ! 
     417      bvm_i(:,:)   = 0._wp 
     418      bv_i (:,:,:) = 0._wp 
    416419      DO jl = 1, jpl 
    417420         DO jk = 1, nlay_i 
    418421            DO jj = 1, jpj 
    419422               DO ji = 1, jpi 
    420                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    421                   tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
    422                      &            / MAX( vt_i(ji,jj) , epsi10 ) 
    423                END DO 
    424             END DO 
    425          END DO 
    426       END DO 
    427       tm_i = tm_i + rt0 
    428  
    429    END SUBROUTINE lim_var_icetm 
    430  
    431  
    432    SUBROUTINE lim_var_bv 
    433       !!------------------------------------------------------------------ 
    434       !!                ***  ROUTINE lim_var_bv *** 
    435       !! 
    436       !! ** Purpose :   computes mean brine volume (%) in sea ice 
    437       !! 
    438       !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
    439       !! 
    440       !! References : Vancoppenolle et al., JGR, 2007 
    441       !!------------------------------------------------------------------ 
    442       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    443       REAL(wp) ::   zbvi             ! local scalars 
    444       !!------------------------------------------------------------------ 
    445       ! 
    446       vt_i (:,:) = 0._wp 
    447       DO jl = 1, jpl 
    448          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    449       END DO 
    450  
    451       bv_i(:,:) = 0._wp 
    452       DO jl = 1, jpl 
    453          DO jk = 1, nlay_i 
    454             DO jj = 1, jpj 
    455                DO ji = 1, jpi 
    456                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) 
    457                   zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   & 
    458                      &                   * v_i(ji,jj,jl) * r1_nlay_i 
    459                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  ) 
    460                   bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 ) 
    461                END DO 
     423                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) 
     424                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  & 
     425                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) 
     426               END DO 
     427            END DO 
     428         END DO 
     429          
     430         DO jj = 1, jpj 
     431            DO ji = 1, jpi 
     432               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     433               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 ) 
    462434            END DO 
    463435         END DO 
     
    683655      INTEGER  :: ji, jk, jl             ! dummy loop indices 
    684656      INTEGER  :: ijpij, i_fill, jl0   
    685       REAL(wp) :: zarg, zV, zconv, zdh 
     657      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    686658      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
    687659      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
     
    704676         IF( zhti(ji) > 0._wp ) THEN 
    705677 
    706          ! initialisation of tests 
    707          itest(:)  = 0 
     678            ! find which category (jl0) the input ice thickness falls into 
     679            jl0 = jpl 
     680            DO jl = 1, jpl 
     681               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     682                  jl0 = jl 
     683                  CYCLE 
     684               ENDIF 
     685            END DO 
     686 
     687            ! initialisation of tests 
     688            itest(:)  = 0 
    708689          
    709          i_fill = jpl + 1                                             !==================================== 
    710          DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
    711             ! iteration                                               !==================================== 
    712             i_fill = i_fill - 1 
     690            i_fill = jpl + 1                                             !==================================== 
     691            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
     692               ! iteration                                               !==================================== 
     693               i_fill = i_fill - 1 
     694                
     695               ! initialisation of ice variables for each try 
     696               zht_i(ji,1:jpl) = 0._wp 
     697               za_i (ji,1:jpl) = 0._wp 
     698               itest(:)        = 0       
     699                
     700               ! *** case very thin ice: fill only category 1 
     701               IF ( i_fill == 1 ) THEN 
     702                  zht_i(ji,1) = zhti(ji) 
     703                  za_i (ji,1) = zai (ji) 
     704                   
     705               ! *** case ice is thicker: fill categories >1 
     706               ELSE 
     707 
     708                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean  
     709                  DO jl = 1, i_fill - 1 
     710                     zht_i(ji,jl) = hi_mean(jl) 
     711                  END DO 
     712                   
     713                  ! Concentrations in the (i_fill-1) categories  
     714                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
     715                  DO jl = 1, i_fill - 1 
     716                     IF ( jl /= jl0 ) THEN 
     717                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     718                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     719                     ENDIF 
     720                  END DO 
     721                   
     722                  ! Concentration in the last (i_fill) category 
     723                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
     724                   
     725                  ! Ice thickness in the last (i_fill) category 
     726                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
     727                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )  
     728                   
     729                  ! clem: correction if concentration of upper cat is greater than lower cat 
     730                  !       (it should be a gaussian around jl0 but sometimes it is not) 
     731                  IF ( jl0 /= jpl ) THEN 
     732                     DO jl = jpl, jl0+1, -1 
     733                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 
     734                           zdv = zht_i(ji,jl) * za_i(ji,jl) 
     735                           zht_i(ji,jl    ) = 0._wp 
     736                           za_i (ji,jl    ) = 0._wp 
     737                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 
     738                        END IF 
     739                     ENDDO 
     740                  ENDIF 
     741                
     742               ENDIF ! case ice is thick or thin 
    713743             
    714             ! initialisation of ice variables for each try 
    715             zht_i(ji,1:jpl) = 0._wp 
    716             za_i (ji,1:jpl) = 0._wp 
     744               !--------------------- 
     745               ! Compatibility tests 
     746               !---------------------  
     747               ! Test 1: area conservation 
     748               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
     749               IF ( zconv < epsi06 ) itest(1) = 1 
    717750             
    718             ! *** case very thin ice: fill only category 1 
    719             IF ( i_fill == 1 ) THEN 
    720                zht_i(ji,1) = zhti(ji) 
    721                za_i (ji,1) = zai (ji) 
    722  
    723             ! *** case ice is thicker: fill categories >1 
    724             ELSE 
    725  
    726                ! Fill ice thicknesses except the last one (i_fill) by hmean  
    727                DO jl = 1, i_fill - 1 
    728                   zht_i(ji,jl) = hi_mean(jl) 
    729                END DO 
     751               ! Test 2: volume conservation 
     752               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
     753               IF ( zconv < epsi06 ) itest(2) = 1 
    730754                
    731                ! find which category (jl0) the input ice thickness falls into 
    732                jl0 = i_fill 
     755               ! Test 3: thickness of the last category is in-bounds ? 
     756               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     757                
     758               ! Test 4: positivity of ice concentrations 
     759               itest(4) = 1 
    733760               DO jl = 1, i_fill 
    734                   IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
    735                      jl0 = jl 
    736            CYCLE 
    737                   ENDIF 
    738                END DO 
    739                 
    740                ! Concentrations in the (i_fill-1) categories  
    741                za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
    742                DO jl = 1, i_fill - 1 
    743                   IF ( jl == jl0 ) CYCLE 
    744                   zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    745                   za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
    746                END DO 
    747                 
    748                ! Concentration in the last (i_fill) category 
    749                za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
    750                 
    751                ! Ice thickness in the last (i_fill) category 
    752                zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
    753                zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
    754                 
    755             ENDIF ! case ice is thick or thin 
    756              
    757             !--------------------- 
    758             ! Compatibility tests 
    759             !---------------------  
    760             ! Test 1: area conservation 
    761             zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
    762             IF ( zconv < epsi06 ) itest(1) = 1 
    763              
    764             ! Test 2: volume conservation 
    765             zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
    766             IF ( zconv < epsi06 ) itest(2) = 1 
    767              
    768             ! Test 3: thickness of the last category is in-bounds ? 
    769             IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
    770              
    771             ! Test 4: positivity of ice concentrations 
    772             itest(4) = 1 
    773             DO jl = 1, i_fill 
    774                IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
    775             END DO             
    776                                                            !============================ 
    777          END DO                                            ! end iteration on categories 
    778                                                            !============================ 
     761                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     762               END DO 
     763               !                                         !============================ 
     764            END DO                                       ! end iteration on categories 
     765               !                                         !============================ 
    779766         ENDIF ! if zhti > 0 
    780767      END DO ! i loop 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6498 r8540  
    1717   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1818   USE sbc_ice         ! Surface boundary condition: ice fields 
    19    USE dom_ice 
    2019   USE ice 
    2120   USE limvar 
     
    3635   !!---------------------------------------------------------------------- 
    3736   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    38    !! $Id$ 
     37   !! $Id: limwri.F90 6965 2016-09-30 13:34:19Z clem $ 
    3938   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4039   !!---------------------------------------------------------------------- 
    4140CONTAINS 
    42  
    43 #if defined key_dimgout 
    44 # include "limwri_dimg.h90" 
    45 #else 
    4641 
    4742   SUBROUTINE lim_wri( kindic ) 
     
    5954      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices 
    6055      REAL(wp) ::  z1_365 
    61       REAL(wp) ::  ztmp 
    62       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s 
    63       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace 
     56      REAL(wp) ::  z2da, z2db, ztmp 
     57      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2 
     58      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace 
    6459      !!------------------------------------------------------------------- 
    6560 
    6661      IF( nn_timing == 1 )  CALL timing_start('limwri') 
    6762 
    68       CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    69       CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi ) 
     63      CALL wrk_alloc( jpi, jpj, jpl, zswi2 ) 
     64      CALL wrk_alloc( jpi, jpj     , z2d, zswi ) 
    7065 
    7166      !----------------------------- 
     
    7469      z1_365 = 1._wp / 365._wp 
    7570 
    76       CALL lim_var_icetm      ! mean sea ice temperature 
    77  
    78       CALL lim_var_bv         ! brine volume 
    79  
    80       DO jj = 1, jpj          ! presence indicator of ice 
     71      ! brine volume 
     72      CALL lim_var_bv  
     73 
     74      ! tresholds for outputs 
     75      DO jj = 1, jpj 
    8176         DO ji = 1, jpi 
    8277            zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) 
    8378         END DO 
    8479      END DO 
    85       ! 
    86       ! 
    87       !                                              
    88       IF ( iom_use( "icethic_cea" ) ) THEN                       ! mean ice thickness 
    89          DO jj = 1, jpj  
     80      DO jl = 1, jpl 
     81         DO jj = 1, jpj 
    9082            DO ji = 1, jpi 
    91                z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 
     83               zswi2(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    9284            END DO 
    9385         END DO 
    94          CALL iom_put( "icethic_cea"  , z2d              ) 
    95       ENDIF 
    96  
    97       IF ( iom_use( "snowthic_cea" ) ) THEN                      ! snow thickness = mean snow thickness over the cell  
    98          DO jj = 1, jpj                                             
    99             DO ji = 1, jpi 
    100                z2d(ji,jj)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 
    101             END DO 
    102          END DO 
    103          CALL iom_put( "snowthic_cea" , z2d              )        
    104       ENDIF 
     86      END DO 
    10587      ! 
     88      ! fluxes  
     89      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
     90      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
     91      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
     92      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
     93      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
     94      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
     95      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
     96      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
     97         &                                                      * a_i_b(:,:,:),dim=3 ) + qemp_ice(:,:) ) 
     98      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
     99      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
     100      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   !emp over ocean (taking into account the snow blown away from the ice) 
     101      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   !emp over ice   (taking into account the snow blown away from the ice) 
     102 
     103      ! velocity 
    106104      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN  
    107105         DO jj = 2 , jpjm1 
    108106            DO ji = 2 , jpim1 
    109                z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
    110                z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     107               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
     108               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     109               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db ) 
    111110           END DO 
    112111         END DO 
    113          CALL lbc_lnk( z2da, 'T', -1. ) 
    114          CALL lbc_lnk( z2db, 'T', -1. ) 
    115          CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component 
    116          CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component 
    117          DO jj = 1, jpj                                  
    118             DO ji = 1, jpi 
    119                z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) )  
    120             END DO 
    121          END DO 
    122          CALL iom_put( "icevel"       , z2d              )       ! ice velocity module 
     112         CALL lbc_lnk( z2d, 'T', 1. ) 
     113         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component 
     114         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component 
     115         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module 
    123116      ENDIF 
    124117      ! 
    125       IF ( iom_use( "miceage" ) ) THEN  
    126          z2d(:,:) = 0.e0 
    127          DO jl = 1, jpl 
    128             DO jj = 1, jpj 
    129                DO ji = 1, jpi 
    130                   rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
    131                   z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 ) 
    132                END DO 
    133             END DO 
    134          END DO 
    135          CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age 
    136       ENDIF 
    137  
    138       IF ( iom_use( "micet" ) ) THEN  
    139          DO jj = 1, jpj 
    140             DO ji = 1, jpi 
    141                z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj) 
    142             END DO 
    143          END DO 
    144          CALL iom_put( "micet"       , z2d               )        ! mean ice temperature 
    145       ENDIF 
     118      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age 
     119      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean 
     120      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean 
     121      IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature 
     122      IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature 
     123      IF ( iom_use( "icecolf" ) )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness 
    146124      ! 
    147       IF ( iom_use( "icest" ) ) THEN  
    148          z2d(:,:) = 0.e0 
    149          DO jl = 1, jpl 
    150             DO jj = 1, jpj 
    151                DO ji = 1, jpi 
    152                   z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) 
    153                END DO 
    154             END DO 
    155          END DO 
    156          CALL iom_put( "icest"       , z2d              )        ! ice surface temperature 
    157       ENDIF 
    158  
    159       IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness 
    160  
    161125      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature 
    162126      CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity 
    163       CALL iom_put( "iceconc"     , at_i                )        ! ice concentration 
    164       CALL iom_put( "icevolu"     , vt_i                )        ! ice volume = mean ice thickness over the cell 
    165       CALL iom_put( "icehc"       , et_i                )        ! ice total heat content 
    166       CALL iom_put( "isnowhc"     , et_s                )        ! snow total heat content 
    167       CALL iom_put( "ibrinv"      , bv_i * 100._wp      )        ! brine volume 
     127      CALL iom_put( "iceconc"     , at_i  * zswi        )        ! ice concentration 
     128      CALL iom_put( "icevolu"     , vt_i  * zswi        )        ! ice volume = mean ice thickness over the cell 
     129      CALL iom_put( "icehc"       , et_i  * zswi        )        ! ice total heat content 
     130      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content 
     131      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume 
    168132      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point 
    169133      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point 
    170134      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation  
    171       CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity 
    172  
    173       CALL iom_put( "icestr"      , strength * 0.001    )        ! ice strength 
    174       CALL iom_put( "idive"       , divu_i * 1.0e8      )        ! divergence 
    175       CALL iom_put( "ishear"      , shear_i * 1.0e8     )        ! shear 
    176       CALL iom_put( "snowvol"     , vt_s                )        ! snow volume 
     135      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity 
     136 
     137      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength 
     138      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence 
     139      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear 
     140      CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume 
    177141       
    178142      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport 
     
    183147 
    184148      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth 
    185       CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt 
    186       CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt 
     149      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting 
     150      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting 
    187151      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation 
    188152      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation 
    189153      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting 
    190       CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual 
     154      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant) 
    191155      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines 
    192156      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation 
     
    202166      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt  
    203167      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt  
     168 
     169      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
     170         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
     171         ELSEWHERE                                       ; z2d = 0._wp 
     172         END WHERE 
     173         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 
     174      ENDIF 
     175 
     176      ztmp = rday / rhosn 
     177      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    204178      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
    205       CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow)  
    206       CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    207        
     179      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice)  
     180      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean       
     181  
    208182      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total) 
    209183      CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics) 
     
    225199      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    226200      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    227       CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base  
     201      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base  
    228202      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    229203      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    230204 
    231  
    232       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
    233          DO jj = 1, jpj  
    234             DO ji = 1, jpi 
    235                z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness 
    236             END DO 
    237          END DO 
    238          WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog 
    239          ELSEWHERE                                   ; z2da = 0._wp 
    240          END WHERE 
    241          CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp ) 
    242       ENDIF 
    243205       
    244206      !-------------------------------- 
    245207      ! Output values for each category 
    246208      !-------------------------------- 
    247       CALL iom_put( "iceconc_cat"      , a_i         )        ! area for categories 
    248       CALL iom_put( "icethic_cat"      , ht_i        )        ! thickness for categories 
    249       CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories 
    250       CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories 
    251  
     209      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories 
     210      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories 
     211      IF ( iom_use( "snowthic_cat" ) )  CALL iom_put( "snowthic_cat"     , ht_s  * zswi2   )        ! snow depth for categories 
     212      IF ( iom_use( "salinity_cat" ) )  CALL iom_put( "salinity_cat"     , sm_i  * zswi2   )        ! salinity for categories 
    252213      ! ice temperature 
    253       IF ( iom_use( "icetemp_cat" ) ) THEN  
    254          zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i 
    255          CALL iom_put( "icetemp_cat"   , zt_i - rt0  ) 
    256       ENDIF 
    257        
     214      IF ( iom_use( "icetemp_cat"  ) )  CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 ) 
    258215      ! snow temperature 
    259       IF ( iom_use( "snwtemp_cat" ) ) THEN  
    260          zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s 
    261          CALL iom_put( "snwtemp_cat"   , zt_s - rt0  ) 
    262       ENDIF 
    263  
    264       ! Compute ice age 
    265       IF ( iom_use( "iceage_cat" ) ) THEN  
    266          DO jl = 1, jpl  
    267             DO jj = 1, jpj 
    268                DO ji = 1, jpi 
    269                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
    270                   rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) ) 
    271                   zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch 
    272                END DO 
    273             END DO 
    274          END DO 
    275          CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories 
    276       ENDIF 
    277  
    278       ! Compute brine volume 
    279       IF ( iom_use( "brinevol_cat" ) ) THEN  
    280          zei(:,:,:) = 0._wp 
    281          DO jl = 1, jpl  
    282             DO jk = 1, nlay_i 
    283                DO jj = 1, jpj 
    284                   DO ji = 1, jpi 
    285                      rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    286                      zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 *  & 
    287                         ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & 
    288                         rswitch * r1_nlay_i 
    289                   END DO 
    290                END DO 
    291             END DO 
    292          END DO 
    293          CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories 
    294       ENDIF 
     216      IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 ) 
     217      ! ice age 
     218      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 ) 
     219      ! brine volume 
     220      IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 ) 
    295221 
    296222      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s 
     
    298224      !     not yet implemented 
    299225       
    300       CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    301       CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db ) 
     226      CALL wrk_dealloc( jpi, jpj, jpl, zswi2 ) 
     227      CALL wrk_dealloc( jpi, jpj     , z2d, zswi ) 
    302228 
    303229      IF( nn_timing == 1 )  CALL timing_stop('limwri') 
    304230       
    305231   END SUBROUTINE lim_wri 
    306 #endif 
    307232 
    308233  
     
    319244      !!   4.0  !  2013-06  (C. Rousset) 
    320245      !!---------------------------------------------------------------------- 
    321       INTEGER, INTENT( in ) ::   kt               ! ocean time-step index) 
    322       INTEGER, INTENT( in ) ::   kid , kh_i        
     246      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index) 
     247      INTEGER, INTENT( in )   ::   kid , kh_i 
     248      INTEGER                 ::   nz_i, jl 
     249      REAL(wp), DIMENSION(jpl) :: jcat 
    323250      !!---------------------------------------------------------------------- 
    324  
    325       CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   & 
    326       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    327       CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   & 
    328       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    329       CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   & 
    330       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    331       CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   & 
    332       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    333       CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   & 
    334       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    335       CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   & 
    336       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    337       CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   & 
    338       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    339       CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   & 
    340       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    341       CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   & 
     251      DO jl = 1, jpl 
     252         jcat(jl) = REAL(jl) 
     253      ENDDO 
     254       
     255      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 
     256 
     257      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   & 
     258      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     259      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   & 
     260      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     261      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   & 
     262      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     263      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   & 
     264      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     265      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   & 
     266      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     267      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   & 
     268      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     269      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   & 
     270      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     271      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   & 
     272      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     273      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   & 
    342274      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    343275      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   & 
    344276      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    345       CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   & 
    346       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    347       CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   & 
    348       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    349       CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   & 
    350       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    351       CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   & 
    352       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    353       CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   & 
    354       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    355       CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   & 
     277      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   & 
     278      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     279      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   & 
     280      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     281      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   & 
     282      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     283 
     284      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   & 
     285      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     286      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   & 
     287      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     288      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   & 
    356289      &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    357       CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   & 
    358       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    359       CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   & 
    360       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    361       CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   & 
    362       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    363       CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   & 
    364       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    365       CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   & 
    366       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    367       CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   & 
    368       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     290      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   & 
     291      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     292      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   & 
     293      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     294      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   & 
     295      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     296      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   & 
     297      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     298 
     299      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   & 
     300      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     301      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   & 
     302      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     303      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   & 
     304      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     305      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   & 
     306      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     307      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   & 
     308      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     309      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   & 
     310      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    369311 
    370312      CALL histend( kid, snc4set )   ! end of the file definition 
    371313 
    372       CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) )     
    373       CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) ) 
    374       CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
    375       CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) ) 
    376       CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) ) 
    377       CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) ) 
    378       CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
    379       CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) ) 
    380       CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) ) 
     314      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )     
     315      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) ) 
     316      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
     317      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) ) 
     318      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) ) 
     319      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) ) 
     320      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
     321      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) ) 
     322      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) ) 
    381323      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) ) 
    382       CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) ) 
    383       CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) ) 
    384       CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    385  
    386       CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) ) 
    387       CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
    388       CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) ) 
    389       CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) ) 
    390       CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) ) 
    391       CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    392       CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    393       CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) ) 
    394       CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
     324      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) ) 
     325      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) ) 
     326      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
     327 
     328      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) ) 
     329      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
     330      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) ) 
     331      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) ) 
     332      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) ) 
     333      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) ) 
     334      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) ) 
     335 
     336      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )     
     337      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )     
     338      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )     
     339      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )     
     340      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )     
     341      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )     
    395342 
    396343      ! Close the file 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90

    r6486 r8540  
    22   !!---------------------------------------------------------------------- 
    33   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    4    !! $Id$ 
     4   !! $Id: limwri_dimg.h90 5123 2015-03-04 16:06:03Z clem $ 
    55   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    66   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r6498 r8540  
    130130   !!---------------------------------------------------------------------- 
    131131   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    132    !! $Id$ 
     132   !! $Id: thd_ice.F90 6399 2016-03-22 17:17:23Z clem $ 
    133133   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    134134   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_lim3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6498 r8540  
    7878   !!---------------------------------------------------------------------- 
    7979   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011) 
    80    !! $Id$ 
     80   !! $Id: sbcice_lim.F90 7778 2017-03-09 14:19:31Z clem $ 
    8181   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    8282   !!---------------------------------------------------------------------- 
     
    229229         CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes 
    230230         ! 
    231          IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs  
     231         IF(ln_limdiaout) CALL lim_diahsb( kt )     ! Diagnostics and outputs  
    232232         ! 
    233233         CALL lim_wri( 1 )                          ! Ice outputs  
     
    310310         numit = nit000 - 1 
    311311      ENDIF 
    312       CALL lim_var_agg(1) 
     312      CALL lim_var_agg(2) 
    313313      CALL lim_var_glo2eqv 
    314314      ! 
    315315      CALL lim_sbc_init                 ! ice surface boundary condition    
     316      ! 
     317      IF( ln_limdiaout) CALL lim_diahsb_init  ! initialization for diags 
    316318      ! 
    317319      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     
    648650CONTAINS 
    649651   SUBROUTINE sbc_ice_lim ( kt, kblk )     ! Dummy routine 
     652      INTEGER, INTENT(in) ::   kt, kblk 
    650653      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk 
    651654   END SUBROUTINE sbc_ice_lim 
Note: See TracChangeset for help on using the changeset viewer.