Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/README
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/README	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/README	(revision 12263)
@@ -0,0 +1,43 @@
+
+-------
+purpose
+-------
+This demonstration case allows for examining the sea ice rheology under idealised forcing.
+The set-up follows the paper "Stress and deformation characteristics of sea ice in a
+high-resolution, anisotropic sea ice model" by Heorton et al, Phil. Trans, R. Soc. 2018.
+
+-----------
+Experiments
+-----------
+A range of experiments can be configured:
+
+- The suggested experiment is to examine the different shear patterns formed by different rheologies. The default set-up reproduces the leftmost panels of figure 3a (ln_rhg_EAP true in namelist_ice_cfg) and 3b (ln_rhg_EVP true in namelist_ice_cfg). Increase the last time step in namelist_cfg appropriately to generate the other panels. 
+
+- The resolution of the box can be changed analogous to figure 2 of the paper (change rn_dx and rn_dy in namelist_cfg, default is 2000m)
+
+- The ratio of the windstress component can be changed to cause different modes of failure (Rwind in usrdef_sbc.F90, default is -0.8) 
+
+It takes about 10 minutes to run a simulation on 16 processors
+
+----------
+How to run
+----------
+
+a) Compile and run the model once to get mesh_mask_xxxx.nc files  
+../../../makenemo -a ICE_RHEO -m X64_MOBILIS -j 4
+mpirun ./nemo -np 1
+
+b) Combine the mesh mask files into one
+../../../tools/REBUILD_NEMO/rebuild_nemo mesh_mask 16
+
+c) Create the initial condition file for sea-ice (initice.nc) by running this python script:
+python ./make_INITICE.py
+
+d) Run the model a second time
+mpirun ./nemo -np 16
+
+---------------
+What to look at
+---------------
+EVP rheology has less defined shear strain maxima (linear kinematic features). If the model is run longer one can see that the LKFs appearing with EVP rheology intersect at a different angle than with EAP rheology
+
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/context_nemo.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/context_nemo.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/context_nemo.xml	(revision 12263)
@@ -0,0 +1,47 @@
+<!--
+ ============================================================================================== 
+    NEMO context
+============================================================================================== 
+-->
+<context id="nemo">
+<!-- $id$ -->
+    <variable_definition>
+    <!-- Year of time origin for NetCDF files; defaults to 1800 -->
+       <variable id="ref_year" type="int"   > 1800 </variable>
+       <variable id="rau0"     type="float" > 1026.0 </variable>
+       <variable id="cpocean"  type="float" > 3991.86795711963 </variable>
+       <variable id="convSpsu" type="float" > 0.99530670233846  </variable>
+       <variable id="rhoic"    type="float" > 917.0 </variable>
+       <variable id="rhosn"    type="float" > 330.0 </variable>
+       <variable id="missval"  type="float" > 1.e20 </variable>
+    </variable_definition>
+<!-- Fields definition -->
+    <field_definition src="./field_def_nemo-oce.xml"/>    <!--  NEMO ocean dynamics     -->
+    <field_definition src="./field_def_nemo-ice.xml"/>    <!--  NEMO sea-ice model      -->
+
+<!-- Files definition -->
+    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics      -->
+    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO sea-ice model       -->
+    <!-- 
+============================================================================================================
+= grid definition = = DO NOT CHANGE =
+============================================================================================================
+    -->
+    
+    <axis_definition>
+      <axis id="deptht"  long_name="Vertical T levels"  unit="m" positive="down" />
+      <axis id="depthu"  long_name="Vertical U levels"  unit="m" positive="down" />
+      <axis id="depthv"  long_name="Vertical V levels"  unit="m" positive="down" />
+      <axis id="depthw"  long_name="Vertical W levels"  unit="m" positive="down" />
+      <axis id="nfloat"  long_name="Float number"       unit="-"                 />
+      <axis id="icbcla"  long_name="Iceberg class"      unit="1"                 />
+      <axis id="ncatice" long_name="Ice category"       unit="1"                 />
+      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"              />
+      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"              />
+    </axis_definition>
+ 
+    <domain_definition src="./domain_def_nemo.xml"/>
+  
+    <grid_definition src="./grid_def_nemo.xml"/>
+
+</context>
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/domain_def_nemo.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/domain_def_nemo.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/domain_def_nemo.xml	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/domain_def_nemo.xml
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-ice.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-ice.xml	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-ice.xml
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-oce.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/field_def_nemo-oce.xml	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-oce.xml
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-ice.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-ice.xml	(revision 12263)
@@ -0,0 +1,60 @@
+<?xml version="1.0"?>
+    <!-- $id$ -->
+    
+ <!--
+============================================================================================================
+=                                           output files definition                                        =
+============================================================================================================
+  -->
+    
+    <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1d" min_digits="4">
+    
+    
+      <file_group id="5ts" output_freq="5ts"  output_level="10" enabled=".TRUE."/>  <!-- 1 day files -->   	
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
+
+      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."> <!-- 1h files -->
+	<file id="file1" name_suffix="_icemod" description="ice variables" enabled=".true." >
+
+        <!-- ice mask -->
+        <field field_ref="icemask"          name="simsk"   />
+
+        <!-- general -->
+        <field field_ref="snwvolu"          name="snvolu" />
+        <field field_ref="icethic"          name="sithic" />
+        <field field_ref="iceconc"          name="siconc" />
+
+        <!-- momentum -->
+        <field field_ref="uice"             name="sivelu"  />
+        <field field_ref="vice"             name="sivelv"  />
+        <field field_ref="icevel"           name="sivelo"  />
+        <field field_ref="utau_ai"          name="utau_ai" />
+        <field field_ref="vtau_ai"          name="vtau_ai" />
+
+        <!-- rheology -->
+        <field field_ref="icediv"           name="sidive"  />
+        <field field_ref="iceshe"           name="sishea"  />
+        <field field_ref="icestr"           name="sistre"  />
+        <field field_ref="normstr"          name="normstr" />
+        <field field_ref="sheastr"          name="sheastr" />
+        <field field_ref="isig1"            name="isig1"   />
+        <field field_ref="isig2"            name="isig2"   />
+        <field field_ref="isig3"            name="isig3"   />
+
+        </file>
+      </file_group>
+
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->        
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->   	
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml	(revision 12263)
@@ -0,0 +1,44 @@
+<?xml version="1.0"?>
+
+<!--
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own filesfor ocean dynamics context                                         =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+   
+    <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1d" min_digits="4">
+    
+
+      <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
+      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."> <!-- 1h files -->
+        <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="uoce"         name="vozocrtx"  />
+          <field field_ref="utau"         name="sozotaux"  />
+        </file>
+
+        <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="voce"         name="vomecrty"  />
+          <field field_ref="vtau"         name="sometauy"  />
+        </file>
+      </file_group>
+
+      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->   
+      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m" output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m" output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+      <file_group id="2y"  output_freq="2y" output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y" output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+
+
+
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/grid_def_nemo.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/grid_def_nemo.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/grid_def_nemo.xml	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/grid_def_nemo.xml
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/iodef.xml
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/iodef.xml	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/iodef.xml	(revision 12263)
@@ -0,0 +1,30 @@
+<?xml version="1.0"?>
+<simulation> 
+
+<!-- ============================================================================================ -->
+<!-- XIOS context                                                                                 -->
+<!-- ============================================================================================ -->
+
+  <context id="xios" >
+
+      <variable_definition>
+
+	  <!--
+          <variable id="optimal_buffer_size"       type="string">memory</variable>
+          <variable id="buffer_size_factor"        type="double">1.0</variable>
+	  -->
+	  <variable id="info_level"                type="int">10</variable>
+	  <variable id="using_server"              type="bool">false</variable>
+	  <variable id="using_oasis"               type="bool">false</variable>
+	  <variable id="oasis_codes_id"            type="string" >oceanx</variable>
+	
+      </variable_definition>
+  </context>
+
+<!-- ============================================================================================ -->
+<!-- NEMO  CONTEXT add and suppress the components you need                                       -->
+<!-- ============================================================================================ -->
+
+  <context id="nemo" src="./context_nemo.xml"/>       <!--  NEMO       -->
+
+</simulation>
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/layout.dat
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/layout.dat	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/layout.dat	(revision 12263)
@@ -0,0 +1,19 @@
+   jpnij   jpimax  jpjmax    jpk  jpiglo  jpjglo ( local:    narea     jpi     jpj )
+      16     252     252       1    1000    1000 ( local:        1     252     252 )
+nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj 
+    0  252  252    1    1  251  251    1    1    4   -1   -1    1     -1     -1
+    1  252  252    2    1  251  251  251    1    5   -1    0    2      0     -1
+    2  251  252    2    1  250  251  501    1    6   -1    1    3      0     -1
+    3  251  252    2    1  251  251  750    1    7   -1    2    4      1     -1
+    4  252  252    1    2  251  251    1  251    8    0    3    5     -1      0
+    5  252  252    2    2  251  251  251  251    9    1    4    6      0      0
+    6  251  252    2    2  250  251  501  251   10    2    5    7      0      0
+    7  251  252    2    2  251  251  750  251   11    3    6    8      1      0
+    8  252  251    1    2  251  250    1  501   12    4    7    9     -1      0
+    9  252  251    2    2  251  250  251  501   13    5    8   10      0      0
+   10  251  251    2    2  250  250  501  501   14    6    9   11      0      0
+   11  251  251    2    2  251  250  750  501   15    7   10   12      1      0
+   12  252  251    1    2  251  251    1  750   -1    8   11   13     -1      1
+   13  252  251    2    2  251  251  251  750   -1    9   12   14      0      1
+   14  251  251    2    2  250  251  501  750   -1   10   13   15      0      1
+   15  251  251    2    2  251  251  750  750   -1   11   14   -1      1      1
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/make_INITICE.py
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/make_INITICE.py	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/make_INITICE.py	(revision 12263)
@@ -0,0 +1,93 @@
+#!/usr/bin/python
+
+import os,sys
+from netCDF4 import Dataset as netcdf
+import numpy as np
+import matplotlib.pyplot as plt
+from math import exp
+from math import ceil
+
+resname=''
+
+# input file
+fcoord='mesh_mask.nc'
+
+# output file
+fflx='initice.nc'
+
+print '   creating init ice file  ' +fflx
+
+# Reading coordinates file
+nccoord=netcdf(fcoord,'r')
+nav_lon=nccoord.variables['nav_lon']
+nav_lat=nccoord.variables['nav_lat']
+time_counter=1
+LON1= nav_lon.shape[1]
+LAT1= nav_lon.shape[0]
+print 'nav_lon.shape[1]' ,nav_lon.shape[1]
+print 'LON1 ', LON1
+print 'LAT1 ', LAT1
+
+# Creating INITICE netcdf file
+nc=netcdf(fflx,'w')
+nc.createDimension('y',LAT1)
+nc.createDimension('x',LON1)
+nc.createDimension('time_counter',None)    # Setting dimension size to 0 or None makes it unlimited. 
+
+cdflon=nc.createVariable('nav_lon','f',('y','x'))
+cdflat=nc.createVariable('nav_lat','f',('y','x'))
+cdftimecounter=nc.createVariable('time_counter','f',('time_counter'))
+
+# ati : Fraction of open waters in sea ice - units %
+# hti : Sea ice thickness - units m
+# hts : Snow thickness - units m
+# smi : Sea ice salinity:
+# tmi : Sea ice internal temperature - units K
+# tsu : Sea ice surface temperature - units K
+#
+# Take constant values from namelist &namiceini of NEMO
+rn_hti_ini=2.0
+rn_hts_ini=0.1            #  initial real snow thickness (m)
+rn_ati_ini=0.99           #  initial ice concentration   (-)
+rn_smi_ini=6.3            #  initial ice salinity     (g/kg)
+rn_tmi_ini=270.           #  initial ice/snw temperature (K)
+rn_tsu_ini=270.           #  initial sea ice temperature (K)
+#
+cdfati=nc.createVariable('ati','f',('time_counter','y','x'))
+cdfati.units='Percentage'
+cdfati.long_name='Sea ice concentration'
+cdfhti=nc.createVariable('hti','f',('time_counter','y','x'))
+cdfhti.long_name='Sea ice thickness'
+cdfhti.units='m'
+cdfhts=nc.createVariable('hts','f',('time_counter','y','x'))
+cdfhts.long_name='Snow thickness'
+cdfhts.units='m'
+cdfsmi=nc.createVariable('smi','f',('time_counter','y','x'))
+cdfsmi.long_name='Sea ice salinity'
+cdfsmi.units='pss'
+cdftmi=nc.createVariable('tmi','f',('time_counter','y','x'))
+cdftmi.long_name='Sea ice internal temperature'
+cdftmi.units='Kelvin'
+cdftsu=nc.createVariable('tsu','f',('time_counter','y','x'))
+cdftsu.long_name='Sea ice surface temperature'
+cdftsu.units='Kelvin'
+
+cdflon[:,:]=nav_lon[:,:]
+cdflat[:,:]=nav_lat[:,:]
+cdftimecounter[0]=1
+
+cdfati[:,:,:]=rn_ati_ini
+cdfhts[:,:,:]=rn_hts_ini
+cdfsmi[:,:,:]=rn_smi_ini
+cdftmi[:,:,:]=rn_tmi_ini
+cdftsu[:,:,:]=rn_tsu_ini
+
+# --- add noise in initial thickness to help nucleation of deformation ---
+for y in np.arange(1,1000,1) :
+    for x in np.arange(1,1000,1) :
+        cdfhti[:,y,x] = rn_hti_ini+0.02*rn_hti_ini*np.random.uniform(-1,1)
+
+nc.close()
+nccoord.close()
+
+#sys.exit()
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_cfg	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_cfg	(revision 12263)
@@ -0,0 +1,273 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!               ICE_RHEO configuration                              !!
+!!======================================================================
+!!              ***  Domain & Run management namelists  ***           !!
+!!                                                                    !!
+!!   namrun       parameters of the run
+!!   namdom       space and time domain
+!!   namcfg       parameters of the configuration                       (default: user defined GYRE)
+!!   namwad       Wetting and drying                                    (default: OFF)
+!!   namtsd       data: temperature & salinity                          (default: OFF)
+!!   namcrs       coarsened grid (for outputs and/or TOP)               (ln_crs =T)
+!!   namc1d       1D configuration options                              ("key_c1d")
+!!   namc1d_dyndmp 1D newtonian damping applied on currents             ("key_c1d")
+!!   namc1d_uvd   1D data (currents)                                    ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namusr_def    !   ICE_RHEO user defined namelist
+!-----------------------------------------------------------------------
+   rn_dx       =    2000.   ! horizontal resolution in meters
+   rn_dy       =    2000.   ! horizontal resolution in meters
+   ln_corio    =    .false. ! set coriolis=0 (false) or not (true)
+      rn_ppgphi0  =  90.    ! Reference latitude [degrees]
+/
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   nn_no       =       0      !  job number (no more used...)
+   cn_exp      =  "ICE_RHEO"  !  experience name
+   nn_it000    =       1      !  first time step
+   nn_itend    =     720     !  last  time step (std 5475)
+   nn_istate   =       1      !  output the initial state (1) or not (0)
+/
+!-----------------------------------------------------------------------
+&namdom        !   time and space domain
+!-----------------------------------------------------------------------
+   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time
+   !
+   rn_rdt      = 30.       !  time step for the dynamics (and tracer if nn_acc=0)
+   !
+   ln_meshmask = .true.    !  =T  create a mesh file
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+   ln_read_cfg = .false.    !  (=T) read the domain configuration file
+      !                     !  (=F) user defined configuration  ==>>>  see usrdef(_...) modules
+      cn_domcfg = "ICE_RHEO_domcfg"    ! domain configuration filename
+/
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***          !!
+!!                                                                    !!
+!!   namsbc          surface boundary condition manager                 (default: NO selection)
+!!   namsbc_flx      flux               formulation                     (ln_flx     =T)
+!!   namsbc_blk      Bulk formulae formulation                          (ln_blk     =T)
+!!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
+!!   namsbc_sas      Stand-Alone Surface module                         (SAS_SRC  only)
+!!   namsbc_iif      Ice-IF: use observed ice cover                     (nn_ice = 1   )
+!!   namtra_qsr      penetrative solar radiation                        (ln_traqsr  =T)
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
+!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
+!!   namsbc_isf      ice shelf melting/freezing                         (ln_isfcav  =T : read (ln_read_cfg=T) or set or usr_def_zgr )
+!!   namsbc_iscpl    coupling option between land ice model and ocean   (ln_isfcav  =T)
+!!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
+!!   namberg         iceberg floats                                     (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_usr      = .true.   !  user defined formulation                  (T => check usrdef_sbc)
+   ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx )
+   ln_blk      = .false.    !  Bulk formulation                          (T => fill namsbc_blk )
+   nn_fsbc     =  1        !  frequency of surface boundary condition computation
+   nn_ice      =  2        !  sea-ice model
+/
+!-----------------------------------------------------------------------
+&namsbc_blk    !   namsbc_blk  generic Bulk formula                     (ln_blk =T)
+!-----------------------------------------------------------------------
+   !                    !  bulk algorithm :
+   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008)
+   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
+   ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
+      !
+      rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
+      rn_zu       = 10.       !  Wind vector reference height (m)
+      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012)
+      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
+      ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
+      rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+      rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
+      rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to
+      !                       !  calculate the wind stress (0.=absolute or 1.=relative winds)
+
+   cn_dir      = './'      !  root directory for the bulk data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !       weights filename               ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
+   sn_wndi     = 'u_10.15JUNE2009_fill'       ,    3.        , 'U_10_MOD',   .true.    , .true. , 'yearly'  , ''                                   , 'Uwnd'   , ''
+   sn_wndj     = 'v_10.15JUNE2009_fill'       ,    3.        , 'V_10_MOD',   .true.    , .true. , 'yearly'  , ''                                   , 'Vwnd'   , ''
+   sn_qsr      = ''                           ,    0.        , 'SWDN_MOD',   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_qlw      = ''                           ,    0.        , 'LWDN_MOD',   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_tair     = ''                           ,    0.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_humi     = ''                           ,    0.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_prec     = ''                           ,    0.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_snow     = ''                           ,    0.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_slp      = ''                           ,    0.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+   sn_tdif     = ''                           ,    0.        , 'taudif'  ,   .false.   , .true. , 'yearly'  , ''                                   , ''       , ''
+
+/
+!-----------------------------------------------------------------------
+&namsbc_sas    !   Stand-Alone Surface module: ocean data               (SAS_SRC  only)
+!-----------------------------------------------------------------------
+   l_sasread   = .false.  ! =T Read in file ;  =F set all to 0. (see sbcssm.F90)
+/
+!!======================================================================
+!!               ***  Lateral boundary condition  ***                 !!
+!!                                                                    !!
+!!   namlbc        lateral momentum boundary condition                  (default: NO selection)
+!!   namagrif      agrif nested grid   (read by child model only)       ("key_agrif")
+!!   nam_tide      Tidal forcing                                        (default: OFF)
+!!   nambdy        Unstructured open boundaries                         (default: OFF)
+!!   nambdy_dta    Unstructured open boundaries - external data         (see  nambdy)
+!!   nambdy_tide   tidal forcing at open boundaries                     (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+   rn_shlat    =    0.     !  free slip
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+/
+!!======================================================================
+!!                ***  Top/Bottom boundary condition  ***             !!
+!!                                                                    !!
+!!   namdrg        top/bottom drag coefficient                          (default: NO selection)
+!!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      (ln_OFF=F)
+!!   nambbc        bottom temperature boundary condition                (default: OFF)
+!!   nambbl        bottom boundary layer scheme                         (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdrg        !   top/bottom drag coefficient                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_OFF     = .true.     !  free-slip       : Cd = 0                  
+/
+!!======================================================================
+!!                        Tracer (T-S) namelists                      !!
+!!                                                                    !!
+!!   nameos        equation of state                                    (default: NO selection)
+!!   namtra_adv    advection scheme                                     (default: NO selection)
+!!   namtra_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namtra_mle    mixed layer eddy param. (Fox-Kemper param.)          (default: OFF)
+!!   namtra_eiv    eddy induced velocity param.                         (default: OFF)
+!!   namtra_dmp    T & S newtonian damping                              (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos        !   ocean Equation Of Seawater                           (default: OFF)
+!-----------------------------------------------------------------------
+   ln_teos10    = .true.         !  = Use TEOS-10 equation of state
+/
+!-----------------------------------------------------------------------
+&namtra_adv    !   advection scheme for tracer                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_traadv_fct =  .true.    !  FCT scheme
+      nn_fct_h   =  4               !  =2/4, horizontal 2nd / 4th order 
+      nn_fct_v   =  2               !  =2/4, vertical   2nd / COMPACT 4th order 
+/
+!----------------------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection)
+!----------------------------------------------------------------------------------
+   ln_traldf_OFF    =  .true.   !    laplacian operator
+/
+!!======================================================================
+!!                      ***  Dynamics namelists  ***                  !!
+!!                                                                    !!
+!!   nam_vvl       vertical coordinate options                          (default: z-star)
+!!   namdyn_adv    formulation of the momentum advection                (default: NO selection)
+!!   namdyn_vor    advection scheme                                     (default: NO selection)
+!!   namdyn_hpg    hydrostatic pressure gradient                        (default: NO selection)
+!!   namdyn_spg    surface pressure gradient                            (default: NO selection)
+!!   namdyn_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namdta_dyn    offline TOP: dynamics read in files                  (OFF_SRC only)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynadv_OFF = .true. !  linear dynamics (no momentum advection)
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg    !   Hydrostatic pressure gradient option                 (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_hpg_zco  = .true.   !  z-coordinate - full steps
+/
+!-----------------------------------------------------------------------
+&namdyn_spg    !   surface pressure gradient                            (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynspg_ts  = .true.  !  split-explicit free surface
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf    !   lateral diffusion on momentum                        (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynldf_OFF = .true.     !  No operator (i.e. no explicit diffusion)
+/
+!!======================================================================
+!!                     vertical physics namelists                     !!
+!!                                                                    !!
+!!    namzdf        vertical physics manager                            (default: NO selection)
+!!    namzdf_ric    richardson number vertical mixing                   (ln_zdfric=T)
+!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
+!!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
+!!    namzdf_osm    OSM vertical diffusion                              (ln_zdfosm=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namzdf        !   vertical physics manager                             (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_zdfcst   = .true.       !  constant mixing
+/
+!!======================================================================
+!!                  ***  Diagnostics namelists  ***                   !!
+!!                                                                    !!
+!!   namtrd       dynamics and/or tracer trends                         (default: OFF)
+!!   namptr       Poleward Transport Diagnostics                        (default: OFF)
+!!   namhsb       Heat and salt budgets                                 (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namflo       float parameters                                      (default: OFF)
+!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF)
+!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF)
+!!   nam_dia25h   25h Mean Output                                       (default: OFF)
+!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4")
+!!======================================================================
+!
+!!======================================================================
+!!               ***  Observation & Assimilation  ***                 !!
+!!                                                                    !!
+!!   namobs       observation and model comparison                      (default: OFF)
+!!   nam_asminc   assimilation increments                               ('key_asminc')
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Miscellaneous namelists  ***                 !!
+!!                                                                    !!
+!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
+!!   namctl            Control prints                                   (default: OFF)
+!!   namsto            Stochastic parametrization of EOS                (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints                                       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsto        ! Stochastic parametrization of EOS                      (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_cfg	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_cfg	(revision 12263)
@@ -0,0 +1,124 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! SI3 configuration namelist: Overwrites SHARED/namelist_ice_ref
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface boundary conditions    (namsbc)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (namthd_pnd)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
+!------------------------------------------------------------------------------
+&nampar         !   Generic parameters
+!------------------------------------------------------------------------------
+   jpl              =   1             !  number of ice  categories
+   ln_icedyn        = .true.          !  ice dynamics (T) or not (F)
+   ln_icethd        = .false.         !  ice thermo   (T) or not (F)
+/
+!------------------------------------------------------------------------------
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn         !   Ice dynamics
+!------------------------------------------------------------------------------
+   ln_dynALL        = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
+   ln_dynADV1D      = .false.         !  dyn.: only advection 1D                  (Schar & Smolarkiewicz 1996 test case)
+   ln_dynADV2D      = .false.         !  dyn.: only advection 2D w prescribed vel.(rn_uvice + advection)
+/
+!------------------------------------------------------------------------------
+&namdyn_rdgrft  !   Ice ridging/rafting
+!------------------------------------------------------------------------------
+   ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)
+      rn_pstar      =   2700          !     ice strength thickness parameter [N/m2]
+      rn_crhg       =   20.0          !     ice strength conc. parameter (-)
+/
+!------------------------------------------------------------------------------
+&namdyn_rhg     !   Ice rheology
+!------------------------------------------------------------------------------
+   ln_rhg_EVP       = .false.          !  EVP rheology
+      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
+      rn_creepl     =   2.0e-9        !     creep limit [1/s]
+      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
+      nn_nevp       = 120             !     number of EVP subcycles                             
+      rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
+                                      !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_EAP       = .true.          !  EVP rheology
+/
+!------------------------------------------------------------------------------
+&namdyn_adv     !   Ice advection
+!------------------------------------------------------------------------------
+   ln_adv_Pra       = .true.         !  Advection scheme (Prather)
+   ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho)
+      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
+/
+!------------------------------------------------------------------------------
+&namsbc         !   Ice surface boundary conditions
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd         !   Ice thermodynamics
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_zdf     !   Ice heat diffusion
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_da      !   Ice lateral melting
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_do      !   Ice growth in open water
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+   ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
+   ln_iceini_file   = .true.         !  netcdf file provided for initialization (T) or not (F)
+   rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze)
+   rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North
+   rn_hts_ini_n     =   0.1           !  initial snow thickness      (m), North
+   rn_ati_ini_n     =   0.99          !  initial ice concentration   (-), North
+   rn_smi_ini_n     =   6.3           !  initial ice salinity     (g/kg), North
+   rn_tmi_ini_n     = 270.            !  initial ice temperature     (K), North
+   rn_tsu_ini_n     = 270.            !  initial surface temperature (K), North
+   rn_tms_ini_n     = 270.            !  initial snw temperature     (K), North
+   rn_apd_ini_n     =   0.0           !  initial pond fraction       (-), North
+   rn_hpd_ini_n     =   0.0           !  initial pond depth          (m), North
+   sn_hti = 'initice'                 , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hts = 'initice'                 , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_ati = 'initice'                 , -12 ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tsu = 'initice'                 , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tmi = 'initice'                 , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_smi = 'initice'                 , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   cn_dir='./'
+/
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+/
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_ref
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_ref	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_ref	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ice_ref
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ref
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ref	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ref	(revision 12263)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ref
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/nemo
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/nemo	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/nemo	(revision 12263)
@@ -0,0 +1,1 @@
+link /hpcdata/working/stefryn/my_branches/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/BLD/bin/nemo.exe
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 12263)
@@ -0,0 +1,1699 @@
+MODULE icedyn_rhg_eap
+   !!======================================================================
+   !!                     ***  MODULE  icedyn_rhg_eap  ***
+   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
+   !!======================================================================
+   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
+   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
+   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
+   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
+   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila) AGRIF
+   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
+   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
+   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
+   !!                 !  2019     (S. Rynders)  change into eap rheology from
+   !!                                           CICE code (Tsamados, Heorton) 
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3'                                       SI3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_dyn_rhg_eap : computes ice velocities from EVP rheology
+   !!   rhg_eap_rst     : read/write EVP fields in ice restart
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constant
+   USE dom_oce        ! Ocean domain
+   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
+   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
+   USE ice            ! sea-ice: ice variables
+   USE icevar         ! ice_var_sshdyn
+   USE icedyn_rdgrft  ! sea-ice: ice strength
+   USE bdy_oce , ONLY : ln_bdy 
+   USE bdyice 
+#if defined key_agrif
+   USE agrif_ice_interp
+#endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+   USE prtctl         ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
+   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
+
+   REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg)
+
+   ! look-up table for calculating structure tensor
+   INTEGER, PARAMETER ::   nx_yield = 41
+   INTEGER, PARAMETER ::   ny_yield = 41
+   INTEGER, PARAMETER ::   na_yield = 21
+
+   REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
+   !! $Id: icedyn_rhg_eap.F90 11536 2019-09-11 13:54:18Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12 )
+      !!-------------------------------------------------------------------
+      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  ***
+      !!                             EAP-C-grid
+      !!
+      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
+      !!  stress and sea-surface slope. Ice-ice interaction is described by 
+      !!  a non-linear elasto-anisotropic-plastic (EAP) law including shear 
+      !!  strength and a bulk rheology .	
+      !!
+      !!  The points in the C-grid look like this, dear reader
+      !!
+      !!                              (ji,jj)
+      !!                                 |
+      !!                                 |
+      !!                      (ji-1,jj)  |  (ji,jj)
+      !!                             ---------   
+      !!                            |         |
+      !!                            | (ji,jj) |------(ji,jj)
+      !!                            |         |
+      !!                             ---------   
+      !!                     (ji-1,jj-1)     (ji,jj-1)
+      !!
+      !! ** Inputs  : - wind forcing (stress), oceanic currents
+      !!                ice total volume (vt_i) per unit area
+      !!                snow total volume (vt_s) per unit area
+      !!
+      !! ** Action  : - compute u_ice, v_ice : the components of the 
+      !!                sea-ice velocity vector
+      !!              - compute delta_i, shear_i, divu_i, which are inputs
+      !!                of the ice thickness distribution
+      !!
+      !! ** Steps   : 0) compute mask at F point
+      !!              1) Compute ice snow mass, ice strength 
+      !!              2) Compute wind, oceanic stresses, mass terms and
+      !!                 coriolis terms of the momentum equation
+      !!              3) Solve the momentum equation (iterative procedure)
+      !!              4) Recompute delta, shear and divergence
+      !!                 (which are inputs of the ITD) & store stress
+      !!                 for the next time step
+      !!              5) Diagnostics including charge ellipse
+      !!
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
+      !!
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
+      !!-------------------------------------------------------------------
+      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
+      !!
+      INTEGER ::   ji, jj, ij0, ij1, ii0, ii1     ! dummy loop indices
+      INTEGER ::   jter         ! local integers
+      !
+      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
+      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
+      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
+      REAL(wp) ::   zalph1, z1_alph1                                    ! alpha coef from Bouillon 2009 or Kimmritz 2017
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2, zdsT ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
+      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
+      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
+      !
+      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                                        ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
+      REAL(wp) ::   zshear, zdum1, zdum2
+      REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components
+      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution
+      REAL(wp) ::   invw                                                ! for test case
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
+!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
+      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
+      !                                                                 !    ice bottom surface if ice is embedded   
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
+
+      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
+      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
+      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
+      !! --- diags
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      !! --- SIMIP diags
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
+      !!-------------------------------------------------------------------
+
+      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
+      !
+!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      !------------------------------------------------------------------------------!
+      ! 0) mask at F points for the ice
+      !------------------------------------------------------------------------------!
+      ! ocean/land mask
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1      ! NO vector opt.
+            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      zwf(:,:) = zfmask(:,:)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         IF( zfmask(1,jj) == 0._wp ) THEN
+            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+         ENDIF
+      END DO
+      CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) define some variables and initialize arrays
+      !------------------------------------------------------------------------------!
+      zrhoco = rau0 * rn_cio 
+!extra code for test case boundary conditions
+      invw=1._wp/(zrhoco*0.5_wp)
+
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
+
+      ! alpha parameters (Bouillon 2009)
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+      ENDIF
+         
+      ! Initialise stress tensor 
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
+
+      ! Ice strength
+      CALL ice_strength
+
+      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      ELSE                         ;   zkt = 0._wp
+      ENDIF
+      !
+      !------------------------------------------------------------------------------!
+      ! 2) Wind / ocean stress, mass terms, coriolis terms
+      !------------------------------------------------------------------------------!
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+
+            ! ice fraction at U-V points
+            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ice/snow mass at U-V points
+            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
+            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+
+            ! Coriolis at T points (m*f)
+            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+            
+            ! Drag ice-atm.
+            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zmf, 'T', 1., zdt_m, 'T', 1. )
+      !
+      !                                  !== Landfast ice parameterization ==!
+      !
+      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ! ice thickness at U-V points
+               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+               ! ice-bottom stress at U points
+               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ! ice_bottom stress at T points
+               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', tau_icebfr(:,:), 'T', 1. )
+         !
+      ELSE                               !-- no landfast
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ztaux_base(ji,jj) = 0._wp
+               ztauy_base(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+
+      !------------------------------------------------------------------------------!
+      ! 3) Solution of the momentum equation, iterative procedure
+      !------------------------------------------------------------------------------!
+      !
+      !                                               ! ==================== !
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
+         !                                            ! ==================== !        
+         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
+         !
+!!$         IF(ln_ctl) THEN   ! Convergence test
+!!$            DO jj = 1, jpjm1
+!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+!!$               zv_ice(:,jj) = v_ice(:,jj)
+!!$            END DO
+!!$         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpim1
+
+               ! shear at F points
+               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)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1. )
+
+         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+            DO ji = 2, jpi ! no vector loop
+
+               ! shear at T points 
+               zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+               !WRITE(numout,*) 'shear output', ji, jj, zdsT
+               
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               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)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e1e2t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! --- anisotropic stress calculation --- !
+               CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, &
+                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), &
+                                       stressptmp, stressmtmp, &
+                                       stress12tmp(ji,jj), strength(ji,jj), &
+                                       alphar, alphas)
+               IF (jter == nn_nevp) THEN
+                  rdg_conv(ji,jj) = -min( alphar, 0._wp)    
+               ENDIF
+
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+               ENDIF
+               
+               ! stress at T points (zkt/=0 if landfast)
+               !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
+               !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+               zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1
+               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1
+               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1
+               !WRITE(numout,*) 'stress output', ji, jj, zs1(ji,jj)
+            END DO
+         END DO
+         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. )
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! stress12tmp at F points 
+               stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
+                  &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
+                  &   ) * 0.25_wp * r1_e1e2f(ji,jj)
+
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+               ENDIF
+               
+               ! stress at F points (zkt/=0 if landfast)
+               !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1
+               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1
+
+            END DO
+         END DO
+
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+               !                   !--- U points
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e1e2u(ji,jj)
+               !
+               !                !--- V points
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e1e2v(ji,jj)
+               !
+               !                !--- ice currents at U-V point
+               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
+               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
+               !
+            END DO
+         END DO
+         !
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) == 0 ) THEN ! even iterations
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/v_ice
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj))
+                  END IF
+
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( 'V' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1          
+                  !                 !--- tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/u_ice
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj))
+                  END IF
+
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( 'U' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            !
+         ELSE ! odd iterations
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/u_ice
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj))
+                  END IF
+
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( 'U' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/v_ice
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj))
+                  END IF
+
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( 'V' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            !
+         ENDIF
+
+!!$         IF(ln_ctl) THEN   ! Convergence test
+!!$            DO jj = 2 , jpjm1
+!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+!!$            END DO
+!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
+!!$            CALL mpp_max( 'icedyn_rhg_eap', zresm )   ! max over the global domain
+!!$         ENDIF
+         !
+         !                                                ! ==================== !
+      END DO                                              !  end loop over jter  !
+      !                                                   ! ==================== !
+      !
+      !------------------------------------------------------------------------------!
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+
+            ! shear at F points
+            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)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            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)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e1e2t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+            
+            ! shear at T points
+            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &             ) * r1_e1e2t(ji,jj)
+            
+            ! delta at T points
+            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      !CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 5) diagnostics
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+         END DO
+      END DO
+
+      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
+         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
+            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         !
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+      ENDIF
+       
+      ! --- divergence, shear and strength --- !
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+
+      ! --- stress tensor --- !
+      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+         !         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
+
+               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+
+               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
+               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         !
+         CALL iom_put( 'isig1' , zsig1 )
+         CALL iom_put( 'isig2' , zsig2 )
+         CALL iom_put( 'isig3' , zsig3 )
+         !
+         ! Stress tensor invariants (normal and shear stress N/m)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+
+         DEALLOCATE( zsig1 , zsig2 , zsig3 )
+      ENDIF
+
+      ! --- anisotropy tensor --- !
+      IF( iom_use('aniso') ) THEN
+         CALL iom_put( 'aniso' , aniso_11 * zmsk00 )
+      ENDIF
+      
+      ! --- SIMIP --- !
+      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
+         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zspgU, 'U', -1., zspgV, 'V', -1., &
+            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+      ENDIF
+
+      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
+         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
+         !
+         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
+            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
+         !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ! 2D ice mass, snow mass, area transport arrays (X, Y)
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
+
+               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+
+         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
+         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
+         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
+         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
+         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
+         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
+
+      ENDIF
+      !
+   END SUBROUTINE ice_dyn_rhg_eap
+
+   SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, &
+                                 shear, a11, a12, &
+                                 stressp,  stressm, &
+                                 stress12, strength, &
+                                 alphar, alphas)
+      !!---------------------------------------------------------------------
+      !!                   ***  SUBROUTINE rhg_eap_rst  ***
+      !! Updates the stress depending on values of strain rate and structure
+      !! tensor and for ksub=ndte it computes closing and sliding rate
+      !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   ) ::   ksub, ndte
+      REAL(wp), INTENT(in   ) ::   phi, strength
+      REAL(wp), INTENT(in   ) ::   divu, tension, shear
+      REAL(wp), INTENT(in   ) ::   a11, a12
+      REAL(wp), INTENT(  out) ::   stressp, stressm, stress12
+      REAL(wp), INTENT(  out) ::   alphar, alphas
+
+      INTEGER  ::   kx ,ky, ka
+
+      REAL(wp) ::   stemp11r, stemp12r, stemp22r   
+      REAL(wp) ::   stemp11s, stemp12s, stemp22s 
+      REAL(wp) ::   a22, Qd11Qd11, Qd11Qd12, Qd12Qd12 
+      REAL(wp) ::   Q11Q11, Q11Q12, Q12Q12 
+      REAL(wp) ::   dtemp11, dtemp12, dtemp22 
+      REAL(wp) ::   rotstemp11r, rotstemp12r, rotstemp22r   
+      REAL(wp) ::   rotstemp11s, rotstemp12s, rotstemp22s
+      REAL(wp) ::   sig11, sig12, sig22 
+      REAL(wp) ::   sgprm11, sgprm12, sgprm22 
+      REAL(wp) ::   invstressconviso 
+      REAL(wp) ::   Angle_denom_gamma,  Angle_denom_alpha 
+      REAL(wp) ::   Tany_1, Tany_2 
+      REAL(wp) ::   x, y, dx, dy, da, kxw, kyw, kaw
+      REAL(wp) ::   invdx, invdy, invda   
+      REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin
+
+      REAL(wp), PARAMETER ::   kfriction = 0.45_wp
+      !!---------------------------------------------------------------------
+      ! Factor to maintain the same stress as in EVP (see Section 3)
+      ! Can be set to 1 otherwise
+      invstressconviso = 1._wp/(1._wp+kfriction*kfriction)
+ 
+      invsin = 1._wp/sin(2._wp*phi) * invstressconviso 
+      !now uses phi as set in higher code
+       
+      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
+      ! rates
+
+      ! 1) structure tensor
+      a22 = 1._wp-a11
+      Q11Q11 = 1._wp
+      Q12Q12 = rsmall
+      Q11Q12 = rsmall
+ 
+      ! gamma: angle between general coordiantes and principal axis of A 
+      ! here Tan2gamma = 2 a12 / (a11 - a22) 
+      if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then      
+         Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + &
+                             4._wp*a12*a12 )
+   
+         Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Cos^2 
+         Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Sin^2
+         Q11Q12 = a12*Angle_denom_gamma                     !CosSin 
+      endif
+ 
+      ! rotation Q*atemp*Q^T
+      atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22 
+      
+      ! make first principal value the largest
+      atempprime = max(atempprime, 1.0_wp - atempprime)
+ 
+      ! 2) strain rate
+      dtemp11 = 0.5_wp*(divu + tension)
+      dtemp12 = shear*0.5_wp
+      dtemp22 = 0.5_wp*(divu - tension)
+
+      !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear
+
+      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
+
+      Qd11Qd11 = 1.0_wp
+      Qd12Qd12 = rsmall
+      Qd11Qd12 = rsmall
+
+      if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then
+
+         Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* &
+                             ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12)
+
+         Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2 
+         Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2
+         Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin
+      endif
+
+      dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22
+      dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22
+      ! In cos and sin values
+      x = 0._wp
+      if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then
+         x = atan2(dtemp2,dtemp1)
+      endif      
+               
+      ! to ensure the angle lies between pi/4 and 9 pi/4 
+      if (x < rpi*0.25_wp) x = x + rpi*2.0_wp
+      
+      ! y: angle between major principal axis of strain rate and structure
+      ! tensor
+      ! y = gamma - alpha 
+      ! Expressesed componently with
+      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
+      
+      Tany_1 = Q11Q12 - Qd11Qd12
+      Tany_2 = Q11Q11 - Qd12Qd12
+      
+      y = 0._wp
+      
+      if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then
+         y = atan2(Tany_1,Tany_2)
+      endif
+      
+      ! to make sure y is between 0 and pi
+      if (y > rpi) y = y - rpi
+      if (y < 0)  y = y + rpi
+
+      ! 3) update anisotropic stress tensor 
+      dx   = rpi/real(nx_yield-1,kind=wp)
+      dy   = rpi/real(ny_yield-1,kind=wp)
+      da   = 0.5_wp/real(na_yield-1,kind=wp)
+      invdx = 1._wp/dx
+      invdy = 1._wp/dy
+      invda = 1._wp/da
+
+      ! % need 8 coords and 8 weights
+      ! % range in kx
+      kx  = int((x-rpi*0.25_wp-rpi)*invdx) + 1
+      kxw = kx - (x-rpi*0.25_wp-rpi)*invdx 
+      
+      ky  = int(y*invdy) + 1
+      kyw = ky - y*invdy 
+      
+      ka  = int((atempprime-0.5_wp)*invda) + 1
+      kaw = ka - (atempprime-0.5_wp)*invda
+      
+      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013)
+      stemp11r =  kxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)
+      stemp12r =  kxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)
+      stemp22r = kxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)
+      
+      stemp11s =  kxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)
+      stemp12s =  kxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)
+      stemp22s =  kxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) &
+        & + (1._wp-kxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) &
+        & + kxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) &
+        & + kxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) &
+        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) &
+        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) &
+        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)
+                   
+      ! Calculate mean ice stress over a collection of floes (Equation 3 in
+      ! Tsamados 2013)
+
+      sig11  = strength*(stemp11r + kfriction*stemp11s) * invsin
+      sig12  = strength*(stemp12r + kfriction*stemp12s) * invsin
+      sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin
+
+      !WRITE(numout,*) 'strength inside loop', strength
+      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22
+
+      ! Back - rotation of the stress from principal axes into general coordinates
+
+      ! Update stress
+      sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 -        2._wp*Q11Q12 *sig12
+      sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12
+      sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 +        2._wp*Q11Q12 *sig12
+
+      stressp  = sgprm11 + sgprm22
+      stress12 = sgprm12
+      stressm  = sgprm11 - sgprm22
+
+      !WRITE(numout,*) 'stress output inside loop', ksub, stressp
+
+      ! Compute ridging and sliding functions in general coordinates 
+      ! (Equation 11 in Tsamados 2013)
+      if (ksub == ndte) then
+         rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r &
+                     + Q12Q12*stemp22r
+         rotstemp12r = Q11Q11*stemp12r +    Q11Q12*(stemp11r-stemp22r) &
+                     - Q12Q12*stemp12r
+         rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r & 
+                     + Q11Q11*stemp22r
+
+         rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s &
+                     + Q12Q12*stemp22s
+         rotstemp12s = Q11Q11*stemp12s +    Q11Q12*(stemp11s-stemp22s) &
+                     - Q12Q12*stemp12s
+         rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s & 
+                     + Q11Q11*stemp22s
+
+         alphar =  rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 &
+                 + rotstemp22r*dtemp22
+         alphas =  rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 &
+                 + rotstemp22s*dtemp22
+      endif
+   END SUBROUTINE update_stress_rdg 
+
+   SUBROUTINE rhg_eap_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_eap_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!      
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter                      ! local integer
+      INTEGER  ::   id1, id2, id3, id4, id5   ! local integers
+      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers
+    
+      INTEGER, PARAMETER            ::    nz = 100
+      
+      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit
+      REAL(wp) ::   da, dx, dy, dp, dz, a1
+
+      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
+            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
+            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
+            id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. )
+            id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+               CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11 )
+               CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12 )
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
+               stress1_i (:,:) = 0._wp
+               stress2_i (:,:) = 0._wp
+               stress12_i(:,:) = 0._wp
+               aniso_11 (:,:) = 0.5_wp
+               aniso_12 (:,:) = 0._wp
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
+            stress1_i (:,:) = 0._wp
+            stress2_i (:,:) = 0._wp
+            stress12_i(:,:) = 0._wp
+            aniso_11 (:,:) = 0.5_wp
+            aniso_12 (:,:) = 0._wp
+         ENDIF
+         !
+
+      da = 0.5_wp/real(na_yield-1,kind=wp)
+      ainit = 0.5_wp - da
+      dx = rpi/real(nx_yield-1,kind=wp)
+      xinit = rpi + 0.25_wp*rpi - dx
+      dz = rpi/real(nz,kind=wp)
+      zinit = -rpi*0.5_wp
+      dy = rpi/real(ny_yield-1,kind=wp)
+      yinit = -dy
+
+      DO ia=1,na_yield
+       DO ix=1,nx_yield
+        DO iy=1,ny_yield
+         s11r(ix,iy,ia) = 0._wp
+         s12r(ix,iy,ia) = 0._wp
+         s22r(ix,iy,ia) = 0._wp
+         s11s(ix,iy,ia) = 0._wp
+         s12s(ix,iy,ia) = 0._wp
+         s22s(ix,iy,ia) = 0._wp
+         IF (ia <= na_yield-1) then
+          DO iz=1,nz
+          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
+          s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          ENDDO
+          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+         ELSE
+          s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+         ENDIF
+        ENDDO
+       ENDDO
+      ENDDO
+
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+         CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11 )
+         CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12 )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_eap_rst
+
+
+   FUNCTION w1(a)
+      !!-------------------------------------------------------------------
+      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp) ::   w1
+      !!-------------------------------------------------------------------
+   
+      w1 = - 223.87569446_wp &
+           + 2361.2198663_wp*a &
+           - 10606.56079975_wp*a*a &
+           + 26315.50025642_wp*a*a*a &
+           - 38948.30444297_wp*a*a*a*a &
+           + 34397.72407466_wp*a*a*a*a*a &
+           - 16789.98003081_wp*a*a*a*a*a*a &
+           + 3495.82839237_wp*a*a*a*a*a*a*a
+   
+   END FUNCTION w1
+
+
+   FUNCTION w2(a)
+      !!-------------------------------------------------------------------
+      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp) ::   w2
+      !!-------------------------------------------------------------------
+   
+      w2 = - 6670.68911883_wp &
+           + 70222.33061536_wp*a &
+           - 314871.71525448_wp*a*a &
+           + 779570.02793492_wp*a*a*a &
+           - 1151098.82436864_wp*a*a*a*a &
+           + 1013896.59464498_wp*a*a*a*a*a &
+           - 493379.44906738_wp*a*a*a*a*a*a &
+           + 102356.551518_wp*a*a*a*a*a*a*a
+   
+   END FUNCTION w2
+
+   FUNCTION s11kr(x,y,z,phi) 
+      !!-------------------------------------------------------------------
+      !! Function : s11kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+   
+      REAL(wp) ::   s11kr, p, pih
+   
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+   
+      pih = 0.5_wp*rpi
+      p = phi
+   
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
+   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else
+   ! x=x-pi/2
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+   
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+   
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+   
+      s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11)
+
+   END FUNCTION s11kr
+
+   FUNCTION s12kr(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s12kr, s12r0, s21r0, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+   
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+   
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+   
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+   
+      s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12)
+      s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21)
+      s12kr=0.5_wp*(s12r0+s21r0)
+   
+   END FUNCTION s12kr
+
+   FUNCTION s22kr(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s22kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s22kr, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22)
+
+   END FUNCTION s22kr
+
+   FUNCTION s11ks(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s11ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s11ks, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11)
+
+   END FUNCTION s11ks
+
+   FUNCTION s12ks(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s12ks, s12s0, s21s0, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p =phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12)
+      s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21)
+      s12ks=0.5_wp*(s12s0+s21s0)
+
+   END FUNCTION s12ks
+
+   FUNCTION s22ks(x,y,z,phi) 
+      !!-------------------------------------------------------------------
+      !! Function : s22ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s22ks, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p =phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) THEN
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) THEN
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22)
+
+   END FUNCTION s22ks
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module           NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+
+#endif
+   !!==============================================================================
+END MODULE icedyn_rhg_eap
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90	(revision 12263)
@@ -0,0 +1,938 @@
+MODULE icedyn_rhg_evp
+   !!======================================================================
+   !!                     ***  MODULE  icedyn_rhg_evp  ***
+   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
+   !!======================================================================
+   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
+   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
+   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
+   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
+   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila) AGRIF
+   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
+   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
+   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3'                                       SI3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
+   !!   rhg_evp_rst     : read/write EVP fields in ice restart
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constant
+   USE dom_oce        ! Ocean domain
+   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
+   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
+   USE ice            ! sea-ice: ice variables
+   USE icevar         ! ice_var_sshdyn
+   USE icedyn_rdgrft  ! sea-ice: ice strength
+   USE bdy_oce , ONLY : ln_bdy 
+   USE bdyice 
+#if defined key_agrif
+   USE agrif_ice_interp
+#endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+   USE prtctl         ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
+   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
+   !! $Id: icedyn_rhg_evp.F90 11536 2019-09-11 13:54:18Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
+      !!-------------------------------------------------------------------
+      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
+      !!                             EVP-C-grid
+      !!
+      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
+      !!  stress and sea-surface slope. Ice-ice interaction is described by 
+      !!  a non-linear elasto-viscous-plastic (EVP) law including shear 
+      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).	
+      !!
+      !!  The points in the C-grid look like this, dear reader
+      !!
+      !!                              (ji,jj)
+      !!                                 |
+      !!                                 |
+      !!                      (ji-1,jj)  |  (ji,jj)
+      !!                             ---------   
+      !!                            |         |
+      !!                            | (ji,jj) |------(ji,jj)
+      !!                            |         |
+      !!                             ---------   
+      !!                     (ji-1,jj-1)     (ji,jj-1)
+      !!
+      !! ** Inputs  : - wind forcing (stress), oceanic currents
+      !!                ice total volume (vt_i) per unit area
+      !!                snow total volume (vt_s) per unit area
+      !!
+      !! ** Action  : - compute u_ice, v_ice : the components of the 
+      !!                sea-ice velocity vector
+      !!              - compute delta_i, shear_i, divu_i, which are inputs
+      !!                of the ice thickness distribution
+      !!
+      !! ** Steps   : 0) compute mask at F point
+      !!              1) Compute ice snow mass, ice strength 
+      !!              2) Compute wind, oceanic stresses, mass terms and
+      !!                 coriolis terms of the momentum equation
+      !!              3) Solve the momentum equation (iterative procedure)
+      !!              4) Recompute delta, shear and divergence
+      !!                 (which are inputs of the ITD) & store stress
+      !!                 for the next time step
+      !!              5) Diagnostics including charge ellipse
+      !!
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
+      !!
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
+      !!-------------------------------------------------------------------
+      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
+      !!
+      INTEGER ::   ji, jj, ij0, ij1, ii0, ii1        ! dummy loop indices
+      INTEGER ::   jter         ! local integers
+      !
+      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
+      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
+      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
+      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
+      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
+      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
+      !
+      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                                        ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
+      REAL(wp) ::   zshear, zdum1, zdum2
+      REAL(wp) ::   invw                                                ! for test case
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
+!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
+      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
+      !                                                                 !    ice bottom surface if ice is embedded   
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
+
+      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
+      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
+      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
+      !! --- diags
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      !! --- SIMIP diags
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
+      !!-------------------------------------------------------------------
+
+      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
+      !
+!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      !------------------------------------------------------------------------------!
+      ! 0) mask at F points for the ice
+      !------------------------------------------------------------------------------!
+      ! ocean/land mask
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1      ! NO vector opt.
+            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+         END DO
+      END DO
+      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      zwf(:,:) = zfmask(:,:)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         IF( zfmask(1,jj) == 0._wp ) THEN
+            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+         ENDIF
+      END DO
+      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) define some variables and initialize arrays
+      !------------------------------------------------------------------------------!
+      zrhoco = rau0 * rn_cio 
+!extra code for test case boundary conditions
+      invw=1._wp/(zrhoco*0.5_wp)
+
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
+
+      ! alpha parameters (Bouillon 2009)
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ENDIF
+         
+      ! Initialise stress tensor 
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
+
+      ! Ice strength
+      CALL ice_strength
+
+      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      ELSE                         ;   zkt = 0._wp
+      ENDIF
+      !
+      !------------------------------------------------------------------------------!
+      ! 2) Wind / ocean stress, mass terms, coriolis terms
+      !------------------------------------------------------------------------------!
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+
+            ! ice fraction at U-V points
+            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ice/snow mass at U-V points
+            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
+            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+
+            ! Coriolis at T points (m*f)
+            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+            
+            ! Drag ice-atm.
+            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
+      !
+      !                                  !== Landfast ice parameterization ==!
+      !
+      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ! ice thickness at U-V points
+               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+               ! ice-bottom stress at U points
+               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ! ice_bottom stress at T points
+               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
+         !
+      ELSE                               !-- no landfast
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ztaux_base(ji,jj) = 0._wp
+               ztauy_base(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+
+      !------------------------------------------------------------------------------!
+      ! 3) Solution of the momentum equation, iterative procedure
+      !------------------------------------------------------------------------------!
+      !
+      !                                               ! ==================== !
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
+         !                                            ! ==================== !        
+         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
+         !
+!!$         IF(ln_ctl) THEN   ! Convergence test
+!!$            DO jj = 1, jpjm1
+!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+!!$               zv_ice(:,jj) = v_ice(:,jj)
+!!$            END DO
+!!$         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpim1
+
+               ! shear at F points
+               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)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
+
+         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+            DO ji = 2, jpi ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               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)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e1e2t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+                  zalph2   = zalph1
+                  z1_alph2 = z1_alph1
+               ENDIF
+               
+               ! stress at T points (zkt/=0 if landfast)
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  zbeta(ji,jj) = zalph2
+               ENDIF
+               
+               ! P/delta at F points
+               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               
+               ! stress at F points (zkt/=0 if landfast)
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+
+            END DO
+         END DO
+
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+               !                   !--- U points
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e1e2u(ji,jj)
+               !
+               !                !--- V points
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e1e2v(ji,jj)
+               !
+               !                !--- ice currents at U-V point
+               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
+               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
+               !
+            END DO
+         END DO
+         !
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) == 0 ) THEN ! even iterations
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/v_ice
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj))
+                  END IF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( 'V' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1          
+                  !                 !--- tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/u_ice
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj))
+                  END IF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( 'U' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            !
+         ELSE ! odd iterations
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/u_ice
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj))
+                  END IF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( 'U' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  !                 !--- tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  !                 !--- Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  !
+                  !                 !--- tau_bottom/v_ice
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
+                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+                  !
+                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
+                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+!extra code for test case boundary conditions
+                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN
+                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj))
+                  END IF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( 'V' )
+#endif
+            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            !
+         ENDIF
+
+!!$         IF(ln_ctl) THEN   ! Convergence test
+!!$            DO jj = 2 , jpjm1
+!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+!!$            END DO
+!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
+!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
+!!$         ENDIF
+         !
+         !                                                ! ==================== !
+      END DO                                              !  end loop over jter  !
+      !                                                   ! ==================== !
+      !
+      !------------------------------------------------------------------------------!
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+
+            ! shear at F points
+            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)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            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)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e1e2t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+            
+            ! shear at T points
+            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &             ) * r1_e1e2t(ji,jj)
+            
+            ! delta at T points
+            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 5) diagnostics
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+         END DO
+      END DO
+
+      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
+         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
+            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         !
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+      ENDIF
+       
+      ! --- divergence, shear and strength --- !
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+
+      ! --- stress tensor --- !
+      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+         !         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
+
+               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+
+               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
+               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         !
+         CALL iom_put( 'isig1' , zsig1 )
+         CALL iom_put( 'isig2' , zsig2 )
+         CALL iom_put( 'isig3' , zsig3 )
+         !
+         ! Stress tensor invariants (normal and shear stress N/m)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+
+         DEALLOCATE( zsig1 , zsig2 , zsig3 )
+      ENDIF
+      
+      ! --- SIMIP --- !
+      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
+         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., &
+            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+      ENDIF
+
+      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
+         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
+         !
+         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
+            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
+         !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ! 2D ice mass, snow mass, area transport arrays (X, Y)
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
+
+               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+
+         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
+         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
+         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
+         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
+         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
+         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
+
+      ENDIF
+      !
+   END SUBROUTINE ice_dyn_rhg_evp
+
+
+   SUBROUTINE rhg_evp_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_evp_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter            ! local integer
+      INTEGER  ::   id1, id2, id3   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
+            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
+            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
+               stress1_i (:,:) = 0._wp
+               stress2_i (:,:) = 0._wp
+               stress12_i(:,:) = 0._wp
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
+            stress1_i (:,:) = 0._wp
+            stress2_i (:,:) = 0._wp
+            stress12_i(:,:) = 0._wp
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_evp_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module           NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!==============================================================================
+END MODULE icedyn_rhg_evp
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90	(revision 12263)
@@ -0,0 +1,151 @@
+MODULE usrdef_hgr
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_hgr  ***
+   !!
+   !!                      ===  ICE_RHEO configuration  ===
+   !!
+   !! User defined :   mesh and Coriolis parameter of a user configuration
+   !!======================================================================
+   !! History :  NEMO  ! 2016-08  (S. Flavoni, G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr    : initialize the horizontal mesh for ICE_RHEO configuration
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE par_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE usrdef_nam, ONLY: rn_dx, rn_dy, ln_corio, rn_ppgphi0   ! horizontal resolution in meters
+   !                                                            coriolis and reference latitude
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_hgr   ! called by domhgr.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_hgr.F90 10074 2018-08-28 16:15:49Z nicolasmartin $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
+      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
+      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain not on the sphere)
+      &                    pe1t  , pe1u  , pe1v   , pe1f   ,   &   ! scale factors       (required)
+      &                    pe2t  , pe2u  , pe2v   , pe2f   ,   &   !
+      &                    ke1e2u_v      , pe1e2u , pe1e2v     )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE usr_def_hgr  ***
+      !!
+      !! ** Purpose :   user defined mesh and Coriolis parameter
+      !!
+      !! ** Method  :   set all intent(out) argument to a proper value
+      !!                ICE_RHEO configuration : uniform grid spacing (rn_dx)
+      !!                without Coriolis force (f=0)
+      !!
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
+      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
+      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors                             [m]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors                             [m]
+      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces computed here, =0 otherwise 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
+      !
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zphi0, zlam0, zbeta, zf0
+      REAL(wp) ::   zti, zui, ztj, zvj   ! local scalars
+      !!-------------------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_hgr : ICE_RHEO configuration bassin'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          f-plane with regular gridcell size'
+
+      !                          ==========
+      zlam0 = -(jpiglo-1)/2 * 1.e-3 * rn_dx
+      zphi0 = -(jpjglo-1)/2 * 1.e-3 * rn_dy
+
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 )
+            zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp
+
+            plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti
+            plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui
+            plamv(ji,jj) = plamt(ji,jj) 
+            plamf(ji,jj) = plamu(ji,jj) 
+   
+            pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj
+            pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj
+            pphiu(ji,jj) = pphit(ji,jj) 
+            pphif(ji,jj) = pphiv(ji,jj) 
+         END DO
+      END DO
+         
+         ! Horizontal scale factors (in meters)
+         !                              ======
+!! ==> EITHER 1) variable scale factors
+!! clem: This can be used with a 1proc simulation but I think it breaks repro when >1procs are used      
+!!         DO jj = 1, jpj
+!!            DO ji = 1, jpi
+!!               !!pe1t(ji,jj) = rn_dx * EXP( -0.8/REAL(jpiglo**2) * (mi0(ji)-REAL(jpiglo+1)*0.5)**2 )  ! gaussian shape
+!!               !!pe2t(ji,jj) = rn_dy * EXP( -0.8/REAL(jpjglo**2) * (mj0(jj)-REAL(jpjglo+1)*0.5)**2 )  ! gaussian shape
+!!               pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5) ) ! linear shape
+!!               pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5) ) ! linear shape
+!!            END DO
+!!         END DO
+!!#if defined key_agrif 
+!!         IF( .NOT. Agrif_Root() ) THEN ! only works if the zoom is positioned at the center of the parent grid
+!!            DO jj = 1, jpj
+!!               DO ji = 1, jpi
+!!                  pe1t(ji,jj) = rn_dx * ( 1. -0.1 * ABS(REAL(mi0(ji))-REAL(jpiglo+1)*0.5) / (1.-REAL(jpiglo+1)*0.5)  &
+!!                     &                            * REAL(jpiglo) / REAL(Agrif_Parent(jpiglo) * Agrif_Rhox()) )       ! factor to match parent grid
+!!                  pe2t(ji,jj) = rn_dy * ( 1. -0.1 * ABS(REAL(mj0(jj))-REAL(jpjglo+1)*0.5) / (1.-REAL(jpjglo+1)*0.5)  &
+!!                     &                            * REAL(jpjglo) / REAL(Agrif_Parent(jpjglo) * Agrif_Rhoy()) )       ! factor to match parent grid
+!!               END DO
+!!            END DO
+!!         ENDIF
+!!#endif
+!! ==> OR 2) constant scale factors
+         pe1t(:,:) = rn_dx
+         pe2t(:,:) = rn_dy
+!! ==> END
+         
+      pe1u(:,:) = pe1t(:,:)      ;      pe2u(:,:) = pe2t(:,:)
+      pe1v(:,:) = pe1t(:,:)      ;      pe2v(:,:) = pe2t(:,:)
+      pe1f(:,:) = pe1t(:,:)      ;      pe2f(:,:) = pe2t(:,:)
+
+      !                             ! NO reduction of grid size in some straits 
+      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
+      pe1e2u(:,:) = 0._wp           !    CAUTION: set to zero to avoid error with some compilers that
+      pe1e2v(:,:) = 0._wp           !             require an initialization of INTENT(out) arguments
+      !
+      !
+      !                       !==  Coriolis parameter  ==!
+      kff = 1                       !  indicate not to compute Coriolis parameter afterward
+      !
+      IF( ln_corio ) THEN
+         zbeta = 2._wp * omega * COS( rad * rn_ppgphi0 ) / ra
+         zf0   = 2._wp * omega * SIN( rad * rn_ppgphi0 )
+         pff_f(:,:) = zf0 + zbeta * pphif(:,:) * 1.e+3
+         pff_t(:,:) = zf0 + zbeta * pphit(:,:) * 1.e+3
+      ELSE
+         pff_f(:,:) = 0.
+         pff_t(:,:) = 0.
+      ENDIF
+      !
+   END SUBROUTINE usr_def_hgr
+
+   !!======================================================================
+END MODULE usrdef_hgr
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_nam.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_nam.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_nam.F90	(revision 12263)
@@ -0,0 +1,108 @@
+MODULE usrdef_nam
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_nam  ***
+   !!
+   !!                      ===  ICE_RHEO configuration  ===
+   !!
+   !! User defined : set the domain characteristics of a user configuration
+   !!======================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni, G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_nam   : read user defined namelist and set global domain size
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce  , ONLY: nimpp , njmpp , Agrif_Root           ! i- & j-indices of the local domain
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_nam   ! called by nemogcm.F90
+
+   !                               !!* namusr_def namelist *!!
+   REAL(wp), PUBLIC ::   rn_dx      ! resolution in meters defining the horizontal domain size
+   REAL(wp), PUBLIC ::   rn_dy      ! resolution in meters defining the horizontal domain size
+   REAL(wp), PUBLIC ::   rn_ppgphi0 ! reference latitude for beta-plane 
+   LOGICAL , PUBLIC ::   ln_corio   ! set coriolis at 0 (ln_corio=F) or not 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_nam.F90 10161 2018-10-01 09:48:55Z clem $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_nam( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read user defined namelist and define the domain size
+      !!
+      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
+      !!
+      !!                Here ICE_RHEO configuration
+      !!
+      !! ** input   : - namusr_def namelist found in namelist_cfg
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name
+      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution
+      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
+      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
+      !
+      INTEGER ::   ios       ! Local integer
+      REAL(wp)::   zlx, zly  ! Local scalars
+      !!
+      NAMELIST/namusr_def/ rn_dx, rn_dy, ln_corio, rn_ppgphi0
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
+      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' )
+      !
+      IF(lwm)   WRITE( numond, namusr_def )
+      !
+      cd_cfg = 'ICE_RHEO'           ! name & resolution (not used)
+      kk_cfg = NINT( rn_dx )
+      !
+      ! Global Domain size:  domain is  2000 km x 2000 Km x 10 m
+         kpi = NINT( 2000.e3 / rn_dx ) 
+         kpj = NINT( 2000.e3 / rn_dy ) 
+      kpk = 1
+      !
+!!      zlx = (kpi-2)*rn_dx*1.e-3
+!!      zly = (kpj-2)*rn_dy*1.e-3
+      zlx = kpi*rn_dx*1.e-3
+      zly = kpj*rn_dy*1.e-3
+      !
+      kperio = 0     ! closed periodicity for the zoom
+      !                             ! control print
+      IF(lwp) THEN
+         WRITE(numout,*) '   '
+         WRITE(numout,*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '   Namelist namusr_def : ICE_RHEO test case'
+         WRITE(numout,*) '      horizontal resolution                    rn_dx  = ', rn_dx, ' meters'
+         WRITE(numout,*) '      horizontal resolution                    rn_dy  = ', rn_dy, ' meters'
+         WRITE(numout,*) '      ICE_RHEO domain = 2000 km x 2000 km'
+         WRITE(numout,*) '         LX [km]: ', zlx
+         WRITE(numout,*) '         LY [km]: ', zly
+         WRITE(numout,*) '         resulting global domain size :        jpiglo = ', kpi
+         WRITE(numout,*) '                                               jpjglo = ', kpj
+         WRITE(numout,*) '                                               jpkglo = ', kpk
+         WRITE(numout,*) '         Coriolis:', ln_corio
+         WRITE(numout,*) '   '
+         WRITE(numout,*) '   Lateral boundary condition of the global domain'
+         WRITE(numout,*) '      ICE_RHEO closed basin                    jperio = ', kperio
+      ENDIF
+      !
+   END SUBROUTINE usr_def_nam
+
+   !!======================================================================
+END MODULE usrdef_nam
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90	(revision 12263)
@@ -0,0 +1,225 @@
+MODULE usrdef_sbc
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_sbc  ***
+   !! 
+   !!                      ===  ICE_RHEO configuration  ===
+   !!
+   !! User defined :   surface forcing of a user configuration
+   !!======================================================================
+   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_sbc    : user defined surface bounday conditions in ICE_RHEO case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbc_ice         ! Surface boundary condition: ice fields
+   USE phycst          ! physical constants
+   USE ice, ONLY       : at_i_b, a_i_b, at_i, u_ice, v_ice
+   USE icethd_dh       ! for CALL ice_thd_snwblow
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usrdef_sbc_oce      ! routine called by sbcmod.F90 for sbc ocean
+   PUBLIC   usrdef_sbc_ice_tau  ! routine called by icestp.F90 for ice dynamics
+   PUBLIC   usrdef_sbc_ice_flx  ! routine called by icestp.F90 for ice thermo
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_sbc.F90 10074 2018-08-28 16:15:49Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usrdef_sbc_oce( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE usr_def_sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the surface boundary
+      !!              condition, i.e. the momentum, heat and freshwater fluxes.
+      !!
+      !! ** Method  :   all 0 fields, for ICE_RHEO case
+      !!                CAUTION : never mask the surface stress field !
+      !!
+      !! ** Action  : - set to ZERO all the ocean surface boundary condition, i.e.   
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER ::    ij0, ij1, ii0, ii1, jj, ji   ! loop indices
+      REAL(wp) ::   zrhoco          ! ocean density and drag coefficient product
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+         !
+         !IF(lwp)   WRITE(numout,*)' usrdef_sbc_oce : ICE_RHEO case: ocean boudary conditions'
+
+         utau(:,:) = 0._wp
+         utau(:,:) = 0._wp
+
+         !ij0 = 1   ;   ij1 = 25                       ! set boundary condition
+         !ii0 = 975   ;   ii1 = 1000
+         !DO jj = mj0(ij0), mj1(ij1)
+         !   DO ji = mi0(ii0), mi1(ii1)
+         !      utau(ji,jj) = -utau_ice(ji,jj) 
+         !      vtau(ji,jj) = -vtau_ice(ji,jj)
+         !   END DO
+         !END DO
+
+         taum(:,:) = 0._wp   ! assume these are not used
+         wndm(:,:) = 0._wp
+         !
+         emp (:,:) = 0._wp
+         sfx (:,:) = 0._wp
+         qns (:,:) = 0._wp
+         qsr (:,:) = 0._wp
+         !
+         utau_b(:,:) = 0._wp 
+         vtau_b(:,:) = 0._wp
+         emp_b (:,:) = 0._wp
+         sfx_b (:,:) = 0._wp
+         qns_b (:,:) = 0._wp
+         !
+      ENDIF
+      !
+   END SUBROUTINE usrdef_sbc_oce
+
+   SUBROUTINE usrdef_sbc_ice_tau( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE usrdef_sbc_ice_tau  ***
+      !!
+      !! ** Purpose :   provide the surface boundary (momentum) condition over
+      !sea-ice
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER ::   jj, ji   ! loop indices
+
+      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f       ! relative wind module and components at F-point
+      REAL(wp) ::   zwndi_t , zwndj_t                 ! relative wind components at T-point
+      REAL(wp), DIMENSION(jpi,jpj) ::   windu, windv  ! wind components (idealised forcing)
+      REAL(wp), PARAMETER ::   r_vfac = 1._wp         ! relative velocity (make 0 for absolute velocity) 
+      REAL(wp), PARAMETER ::   Rwind = -0.8_wp        ! ratio of wind components
+      REAL(wp), PARAMETER ::   Umax = 15._wp          ! maximum wind speed (m/s)
+      REAL(wp), PARAMETER ::   d = 2000._wp           ! size of the domain (km)
+      REAL(wp), PARAMETER ::   res = 2._wp            ! gridcell size
+      REAL(wp), PARAMETER ::   zrhoa  = 1.22          ! Air density kg/m3
+      REAL(wp), PARAMETER ::   Cd_atm =  1.4e-3       ! transfer coefficient over ice
+      !!---------------------------------------------------------------------
+      ! extra code for test case
+      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: analytical stress forcing'
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1   
+            ! wind spins up over 6 hours, factor 1000 to balance the units
+            windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
+            windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
+
+      wndm_ice(:,:) = 0._wp      !!gm brutal....
+
+      ! ------------------------------------------------------------ !
+      !    Wind module relative to the moving ice ( U10m - U_ice )   !
+      ! ------------------------------------------------------------ !
+      ! C-grid ice dynamics :   U & V-points (same as ocean)
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1   
+            zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+            zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
+            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1. )
+
+      !!gm brutal....
+      utau_ice  (:,:) = 0._wp
+      vtau_ice  (:,:) = 0._wp
+      !!gm end
+
+      ! ------------------------------------------------------------ !
+      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
+      ! ------------------------------------------------------------ !
+      ! C-grid ice dynamics :   U & V-points (same as ocean)
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1   
+            utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
+               &          * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
+            vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
+               &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+      !
+   END SUBROUTINE usrdef_sbc_ice_tau
+
+   SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE usrdef_sbc_ice_flx  ***
+      !!
+      !! ** Purpose :   provide the surface boundary (flux) condition over sea-ice
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi    ! ice thickness
+      !!
+      REAL(wp) ::   zfr1, zfr2                 ! local variables
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsnw   ! snw distribution after wind blowing
+      !!---------------------------------------------------------------------
+      !
+      IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: NO flux forcing'
+      !
+      ! ocean variables (renaming)
+      emp_oce (:,:)   = 0._wp   ! uniform value for freshwater budget (E-P)
+      qsr_oce (:,:)   = 0._wp   ! uniform value for     solar radiation
+      qns_oce (:,:)   = 0._wp   ! uniform value for non-solar radiation
+
+      ! ice variables
+      alb_ice (:,:,:) = 0.7_wp  ! useless
+      qsr_ice (:,:,:) = 0._wp   ! uniform value for     solar radiation
+      qns_ice (:,:,:) = 0._wp   ! uniform value for non-solar radiation
+      sprecip (:,:)   = 0._wp   ! uniform value for snow precip
+      evap_ice(:,:,:) = 0._wp   ! uniform value for sublimation
+
+      ! ice fields deduced from above
+      zsnw(:,:) = 1._wp
+      !!CALL lim_thd_snwblow( at_i_b, zsnw )  ! snow distribution over ice after wind blowing 
+      emp_ice  (:,:)   = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw(:,:)
+      emp_oce  (:,:)   = emp_oce(:,:) - sprecip(:,:) * (1._wp - zsnw(:,:) )
+      qevap_ice(:,:,:) =   0._wp
+      qprec_ice(:,:)   =   rhos * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) !  in J/m3
+      qemp_oce (:,:)   = - emp_oce(:,:) * sst_m(:,:) * rcp
+      qemp_ice (:,:)   =   sprecip(:,:) * zsnw * ( sst_m(:,:) * rcpi - rLfus ) * tmask(:,:,1) ! solid precip (only)
+
+      ! total fluxes
+      emp_tot (:,:) = emp_ice  + emp_oce
+      qns_tot (:,:) = at_i_b(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
+      qsr_tot (:,:) = at_i_b(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
+
+      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
+      !
+      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
+      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
+         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
+      ELSEWHERE                                                         ! zero when hs>0
+         qtr_ice_top(:,:,:) = 0._wp 
+      END WHERE
+          
+   END SUBROUTINE usrdef_sbc_ice_flx
+
+
+   !!======================================================================
+END MODULE usrdef_sbc
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_zgr.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_zgr.F90	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_zgr.F90	(revision 12263)
@@ -0,0 +1,111 @@
+MODULE usrdef_zgr
+   !!======================================================================
+   !!                   ***  MODULE  usrdef_zgr  ***
+   !!
+   !!                   ===      ICE_RHEO case     ===
+   !!
+   !! Ocean domain : user defined vertical coordinate system 
+   !!======================================================================
+   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_zgr   : user defined vertical coordinate system (required)
+   !!---------------------------------------------------------------------
+   USE oce            ! ocean variables
+   USE usrdef_nam     ! User defined : namelist variables
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! distributed memory computing library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_zgr   ! called by domzgr.F90
+
+  !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_zgr.F90 10074 2018-08-28 16:15:49Z nicolasmartin $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS             
+
+   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
+      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
+      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
+      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
+      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
+      &                    k_top  , k_bot    )                             ! top & bottom ocean level
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE usr_def_zgr  ***
+      !!
+      !! ** Purpose :   User defined the vertical coordinates
+      !!
+      !!----------------------------------------------------------------------
+      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
+      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors 
+      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
+      !
+      INTEGER  ::   jk, k_dz  ! dummy indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : ICE_RHEO configuration '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (slab ocean - ice cover in square box domain)'
+      !
+      !
+      ! type of vertical coordinate  ==>>>   here ICE_RHEO : slab ocean always
+      ! ---------------------------
+      ld_zco    = .TRUE.       ! z-full-step coordinate
+      ld_zps    = .FALSE.      ! z-partial-step coordinate
+      ld_sco    = .FALSE.      ! s-coordinate
+      ld_isfcav = .FALSE.      ! ISF Ice Shelves Flag
+      !
+      !
+      ! Build the vertical coordinate system
+      ! ------------------------------------
+      !
+      !                       !==  UNmasked meter bathymetry  ==!
+      !
+      ! 
+      k_dz = 1
+      DO jk = 1, jpk
+         pdepw_1d(jk) =    k_dz
+         pdept_1d(jk) =    k_dz
+         pe3w_1d (jk) =    k_dz
+         pe3t_1d (jk) =    k_dz
+      END DO
+      !                       !==  top masked level bathymetry  ==!  (all coordinates)
+      !
+      ! no ocean cavities : top ocean level is ONE, except over land
+      k_top(:,:) = 1
+      !
+      !                       !==  z-coordinate  ==!   (step-like topography)
+      !                                !* bottom ocean compute from the depth of grid-points
+      jpkm1 = jpk
+      k_bot(:,:) = 1    ! here use k_top as a land mask
+      !                                !* horizontally uniform coordinate (reference z-co everywhere)
+      DO jk = 1, jpk
+         pdept(:,:,jk) = pdept_1d(jk)
+         pdepw(:,:,jk) = pdepw_1d(jk)
+         pe3t (:,:,jk) = pe3t_1d (jk)
+         pe3u (:,:,jk) = pe3t_1d (jk)
+         pe3v (:,:,jk) = pe3t_1d (jk)
+         pe3f (:,:,jk) = pe3t_1d (jk)
+         pe3w (:,:,jk) = pe3w_1d (jk)
+         pe3uw(:,:,jk) = pe3w_1d (jk)
+         pe3vw(:,:,jk) = pe3w_1d (jk)
+      END DO
+      !
+   END SUBROUTINE usr_def_zgr
+
+   !!======================================================================
+END MODULE usrdef_zgr
Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/cpp_ICE_RHEO.fcm
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/cpp_ICE_RHEO.fcm	(revision 12263)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/cpp_ICE_RHEO.fcm	(revision 12263)
@@ -0,0 +1,1 @@
+bld::tool::fppkeys key_si3 key_mpp_mpi key_nosignedzero key_iomput
