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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_evp.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

source: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90 @ 15499

Last change on this file since 15499 was 15499, checked in by vancop, 8 months ago

northfold and stress diagnostic issues in VP rheology

  • Property svn:keywords set to Id
File size: 64.3 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
10   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila) AGRIF
12   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icevar         ! ice_var_sshdyn
29   USE icedyn_rdgrft  ! sea-ice: ice strength
30   USE bdy_oce , ONLY : ln_bdy 
31   USE bdyice 
32#if defined key_agrif
33   USE agrif_ice_interp
34#endif
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE lib_mpp        ! MPP library
39   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
40   USE lbclnk         ! lateral boundary conditions (or mpp links)
41   USE prtctl         ! Print control
42
43   USE netcdf         ! NetCDF library for convergence test
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
48   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
49
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52
53   !! for convergence tests
54   INTEGER ::   ncvgid   ! netcdf file id
55   INTEGER ::   nvarid   ! netcdf variable id
56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
57   !!----------------------------------------------------------------------
58   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
59   !! $Id$
60   !! Software governed by the CeCILL license (see ./LICENSE)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
65      !!-------------------------------------------------------------------
66      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
67      !!                             EVP-C-grid
68      !!
69      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
70      !!  stress and sea-surface slope. Ice-ice interaction is described by
71      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
72      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
73      !!
74      !!  The points in the C-grid look like this, dear reader
75      !!
76      !!                              (ji,jj)
77      !!                                 |
78      !!                                 |
79      !!                      (ji-1,jj)  |  (ji,jj)
80      !!                             ---------   
81      !!                            |         |
82      !!                            | (ji,jj) |------(ji,jj)
83      !!                            |         |
84      !!                             ---------   
85      !!                     (ji-1,jj-1)     (ji,jj-1)
86      !!
87      !! ** Inputs  : - wind forcing (stress), oceanic currents
88      !!                ice total volume (vt_i) per unit area
89      !!                snow total volume (vt_s) per unit area
90      !!
91      !! ** Action  : - compute u_ice, v_ice : the components of the
92      !!                sea-ice velocity vector
93      !!              - compute delta_i, shear_i, divu_i, which are inputs
94      !!                of the ice thickness distribution
95      !!
96      !! ** Steps   : 0) compute mask at F point
97      !!              1) Compute ice snow mass, ice strength
98      !!              2) Compute wind, oceanic stresses, mass terms and
99      !!                 coriolis terms of the momentum equation
100      !!              3) Solve the momentum equation (iterative procedure)
101      !!              4) Recompute delta, shear and divergence
102      !!                 (which are inputs of the ITD) & store stress
103      !!                 for the next time step
104      !!              5) Diagnostics including charge ellipse
105      !!
106      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
107      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
108      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
109      !!              (i.e. more stable and better convergence)
110      !!
111      !! References : Hunke and Dukowicz, JPO97
112      !!              Bouillon et al., Ocean Modelling 2009
113      !!              Bouillon et al., Ocean Modelling 2013
114      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
115      !!-------------------------------------------------------------------
116      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
117      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
118      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
119      !!
120      INTEGER ::   ji, jj       ! dummy loop indices
121      INTEGER ::   jter         ! local integers
122      !
123      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
124      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
125      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
126      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
127      REAl(wp) ::   zbetau, zbetav
128      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
129      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
130      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
131      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
132      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
133      !
134      REAL(wp) ::   zintb, zintn                                        ! dummy argument
135      REAL(wp) ::   zfac_x, zfac_y
136      REAL(wp) ::   zshear, zdum1, zdum2
137      !
138      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points
140      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension
141      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
142      !
143      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
144      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
145      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
146      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
147      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
148      !
149      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
150      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
151      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
152      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
153      !                                                                 !    ice bottom surface if ice is embedded   
154      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
155      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
156      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
157      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
158      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
159      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
160      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
161      !
162      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
164      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
165
166      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
167      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
168      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
169      !! --- check convergence
170      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
171      !! --- diags
172      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
174      !! --- SIMIP diags
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
181      !!-------------------------------------------------------------------
182
183      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
184      !
185      ! for diagnostics and convergence tests
186      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
187      DO jj = 1, jpj
188         DO ji = 1, jpi
189            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
190            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
191         END DO
192      END DO
193      !
194      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
195      !------------------------------------------------------------------------------!
196      ! 0) mask at F points for the ice
197      !------------------------------------------------------------------------------!
198      ! ocean/land mask
199      DO jj = 1, jpjm1
200         DO ji = 1, jpim1      ! NO vector opt.
201            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
202         END DO
203      END DO
204      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
205
206      ! Lateral boundary conditions on velocity (modify zfmask)
207      DO jj = 2, jpjm1
208         DO ji = fs_2, fs_jpim1   ! vector opt.
209            IF( zfmask(ji,jj) == 0._wp ) THEN
210               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
211                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
212            ENDIF
213         END DO
214      END DO
215      DO jj = 2, jpjm1
216         IF( zfmask(1,jj) == 0._wp ) THEN
217            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
218         ENDIF
219         IF( zfmask(jpi,jj) == 0._wp ) THEN
220            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
221        ENDIF
222      END DO
223      DO ji = 2, jpim1
224         IF( zfmask(ji,1) == 0._wp ) THEN
225            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
226         ENDIF
227         IF( zfmask(ji,jpj) == 0._wp ) THEN
228            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
229         ENDIF
230      END DO
231      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
232
233      !------------------------------------------------------------------------------!
234      ! 1) define some variables and initialize arrays
235      !------------------------------------------------------------------------------!
236      zrhoco = rau0 * rn_cio 
237
238      ! ecc2: square of yield ellipse eccenticrity
239      ecc2    = rn_ecc * rn_ecc
240      z1_ecc2 = 1._wp / ecc2
241
242      ! alpha parameters (Bouillon 2009)
243      IF( .NOT. ln_aEVP ) THEN
244         zdtevp   = rdt_ice / REAL( nn_nevp )
245         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
246         zalph2 = zalph1 * z1_ecc2
247
248         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
249         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
250      ELSE
251         zdtevp   = rdt_ice
252         ! zalpha parameters set later on adaptatively
253      ENDIF
254      z1_dtevp = 1._wp / zdtevp
255         
256      ! Initialise stress tensor
257      zs1 (:,:) = pstress1_i (:,:) 
258      zs2 (:,:) = pstress2_i (:,:)
259      zs12(:,:) = pstress12_i(:,:)
260
261      ! Ice strength
262      CALL ice_strength
263
264      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
265      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
266      ELSE                         ;   zkt = 0._wp
267      ENDIF
268      !
269      !------------------------------------------------------------------------------!
270      ! 2) Wind / ocean stress, mass terms, coriolis terms
271      !------------------------------------------------------------------------------!
272      ! sea surface height
273      !    embedded sea ice: compute representative ice top surface
274      !    non-embedded sea ice: use ocean surface for slope calculation
275      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
276
277      DO jj = 2, jpjm1
278         DO ji = fs_2, fs_jpim1
279
280            ! ice fraction at U-V points
281            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)
282            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)
283
284            ! Ice/snow mass at U-V points
285            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
286            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
287            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
288            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
289            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
290
291            ! Ocean currents at U-V points
292            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)
293            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)
294
295            ! Coriolis at T points (m*f)
296            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
297
298            ! dt/m at T points (for alpha and beta coefficients)
299            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
300           
301            ! m/dt
302            zmU_t(ji,jj)    = zmassU * z1_dtevp
303            zmV_t(ji,jj)    = zmassV * z1_dtevp
304           
305            ! Drag ice-atm.
306            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
307            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
308
309            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
310            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
311            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
312
313            ! masks
314            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
315            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
316
317            ! switches
318            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
319            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
320            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
321            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
322
323         END DO
324      END DO
325      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
326      !
327      !                                  !== Landfast ice parameterization ==!
328      !
329      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
330         DO jj = 2, jpjm1
331            DO ji = fs_2, fs_jpim1
332               ! ice thickness at U-V points
333               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)
334               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)
335               ! ice-bottom stress at U points
336               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj)
337               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
338               ! ice-bottom stress at V points
339               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj)
340               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
341               ! ice_bottom stress at T points
342               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj)
343               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
344            END DO
345         END DO
346         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
347         !
348      ELSE                               !-- no landfast
349         DO jj = 2, jpjm1
350            DO ji = fs_2, fs_jpim1
351               ztaux_base(ji,jj) = 0._wp
352               ztauy_base(ji,jj) = 0._wp
353            END DO
354         END DO
355      ENDIF
356
357      !------------------------------------------------------------------------------!
358      ! 3) Solution of the momentum equation, iterative procedure
359      !------------------------------------------------------------------------------!
360      !
361      !                                               ! ==================== !
362      DO jter = 1 , nn_nevp                           !    loop over jter    !
363         !                                            ! ==================== !       
364         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
365         !
366         ! convergence test
367         IF( ln_rhg_chkcvg ) THEN
368            DO jj = 1, jpj
369               DO ji = 1, jpi
370                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
371                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
372               END DO
373            END DO
374         ENDIF
375
376         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
377         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
378            DO ji = 1, jpim1
379
380               ! shear at F points
381               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)   &
382                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
383                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
384
385            END DO
386         END DO
387         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
388
389         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
390            DO ji = 2, jpi ! no vector loop
391
392               ! shear**2 at T points (doc eq. A16)
393               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
394                  &   + 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)  &
395                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
396             
397               ! divergence at T points
398               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
399                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
400                  &    ) * r1_e1e2t(ji,jj)
401               zdiv2 = zdiv * zdiv
402               
403               ! tension at T points
404               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)   &
405                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
406                  &   ) * r1_e1e2t(ji,jj)
407               zdt2 = zdt * zdt
408               
409               ! delta at T points
410               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
411               
412               zdeltastar_t(ji,jj) = zdelta + rn_creepl                     ! store delta* at previous iterate (for ellipse computation)
413
414               ! P/delta at T points
415               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
416
417               ! alpha for aEVP
418               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
419               !   alpha = beta = sqrt(4*gamma)
420               IF( ln_aEVP ) THEN
421                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
422                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
423                  zalph2   = zalph1
424                  z1_alph2 = z1_alph1
425                  ! explicit:
426                  ! z1_alph1 = 1._wp / zalph1
427                  ! z1_alph2 = 1._wp / zalph1
428                  ! zalph1 = zalph1 - 1._wp
429                  ! zalph2 = zalph1
430               ENDIF
431               
432               ! stress at T points (zkt/=0 if landfast)
433               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
434               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
435             
436            END DO
437         END DO
438         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
439
440         ! Save beta at T-points for further computations
441         IF( ln_aEVP ) THEN
442            DO jj = 1, jpj
443               DO ji = 1, jpi
444                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
445               END DO
446            END DO
447         ENDIF
448         
449         DO jj = 1, jpjm1
450            DO ji = 1, jpim1
451
452               ! alpha for aEVP
453               IF( ln_aEVP ) THEN
454                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
455                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
456                  ! explicit:
457                  ! z1_alph2 = 1._wp / zalph2
458                  ! zalph2 = zalph2 - 1._wp
459               ENDIF
460               
461               ! P/delta at F points
462               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) )
463               
464               ! stress at F points (zkt/=0 if landfast)
465               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
466
467            END DO
468         END DO
469
470         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
471         DO jj = 2, jpjm1
472            DO ji = fs_2, fs_jpim1               
473               !                   !--- U points
474               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
475                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
476                  &                    ) * r1_e2u(ji,jj)                                                                      &
477                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
478                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
479                  &                  ) * r1_e1e2u(ji,jj)
480               !
481               !                !--- V points
482               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
483                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
484                  &                    ) * r1_e1v(ji,jj)                                                                      &
485                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
486                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
487                  &                  ) * r1_e1e2v(ji,jj)
488               !
489               !                !--- ice currents at U-V point
490               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)
491               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)
492               !
493            END DO
494         END DO
495         !
496         ! --- Computation of ice velocity --- !
497         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
498         !  Bouillon et al. 2009 (eq 34-35) => stable
499         IF( MOD(jter,2) == 0 ) THEN ! even iterations
500            !
501            DO jj = 2, jpjm1
502               DO ji = fs_2, fs_jpim1
503                  !                 !--- tau_io/(v_oce - v_ice)
504                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
505                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
506                  !                 !--- Ocean-to-Ice stress
507                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
508                  !
509                  !                 !--- tau_bottom/v_ice
510                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
511                  zTauB = ztauy_base(ji,jj) / zvel
512                  !                 !--- OceanBottom-to-Ice stress
513                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
514                  !
515                  !                 !--- Coriolis at V-points (energy conserving formulation)
516                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
517                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
518                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
519                  !
520                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
521                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
522                  !
523                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
524                  !                                         1 = sliding friction : TauB < RHS
525                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
526                  !
527                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
528                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
529                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
530                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
531                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
532                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
533                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
534                        &                                    ) / ( zbetav + 1._wp )                                              &
535                        &             ) * 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
536                        &           )   * zmsk00y(ji,jj)
537                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
538                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
539                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
540                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
541                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
542                        &             ) * 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
543                        &            )  * zmsk00y(ji,jj)
544                  ENDIF
545               END DO
546            END DO
547            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
548            !
549#if defined key_agrif
550!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
551            CALL agrif_interp_ice( 'V' )
552#endif
553            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
554            !
555            DO jj = 2, jpjm1
556               DO ji = fs_2, fs_jpim1         
557                  !                 !--- tau_io/(u_oce - u_ice)
558                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
559                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
560                  !                 !--- Ocean-to-Ice stress
561                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
562                  !
563                  !                 !--- tau_bottom/u_ice
564                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
565                  zTauB = ztaux_base(ji,jj) / zvel
566                  !                 !--- OceanBottom-to-Ice stress
567                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
568                  !
569                  !                 !--- Coriolis at U-points (energy conserving formulation)
570                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
571                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
572                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
573                  !
574                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
575                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
576                  !
577                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
578                  !                                         1 = sliding friction : TauB < RHS
579                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
580                  !
581                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
582                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
583                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
584                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
585                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
586                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
587                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
588                        &                                    ) / ( zbetau + 1._wp )                                              &
589                        &             ) * 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
590                        &           )   * zmsk00x(ji,jj)
591                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
592                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
593                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
594                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
595                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
596                        &             ) * 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
597                        &           )   * zmsk00x(ji,jj)
598                  ENDIF
599               END DO
600            END DO
601            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
602            !
603#if defined key_agrif
604!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
605            CALL agrif_interp_ice( 'U' )
606#endif
607            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
608            !
609         ELSE ! odd iterations
610            !
611            DO jj = 2, jpjm1
612               DO ji = fs_2, fs_jpim1
613                  !                 !--- tau_io/(u_oce - u_ice)
614                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
615                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
616                  !                 !--- Ocean-to-Ice stress
617                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
618                  !
619                  !                 !--- tau_bottom/u_ice
620                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
621                  zTauB = ztaux_base(ji,jj) / zvel
622                  !                 !--- OceanBottom-to-Ice stress
623                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
624                  !
625                  !                 !--- Coriolis at U-points (energy conserving formulation)
626                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
627                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
628                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
629                  !
630                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
631                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
632                  !
633                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
634                  !                                         1 = sliding friction : TauB < RHS
635                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
636                  !
637                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
638                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
639                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
640                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
641                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
642                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
643                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
644                        &                                    ) / ( zbetau + 1._wp )                                              &
645                        &             ) * 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
646                        &           )   * zmsk00x(ji,jj)
647                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
648                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
649                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
650                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
651                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
652                        &             ) * 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
653                        &           )   * zmsk00x(ji,jj)
654                  ENDIF
655               END DO
656            END DO
657            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
658            !
659#if defined key_agrif
660!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
661            CALL agrif_interp_ice( 'U' )
662#endif
663            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
664            !
665            DO jj = 2, jpjm1
666               DO ji = fs_2, fs_jpim1
667                  !                 !--- tau_io/(v_oce - v_ice)
668                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
669                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
670                  !                 !--- Ocean-to-Ice stress
671                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
672                  !
673                  !                 !--- tau_bottom/v_ice
674                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
675                  zTauB = ztauy_base(ji,jj) / zvel
676                  !                 !--- OceanBottom-to-Ice stress
677                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
678                  !
679                  !                 !--- Coriolis at v-points (energy conserving formulation)
680                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
681                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
682                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
683                  !
684                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
685                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
686                  !
687                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
688                  !                                         1 = sliding friction : TauB < RHS
689                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
690                  !
691                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
692                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
693                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
694                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
695                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
696                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
697                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
698                        &                                    ) / ( zbetav + 1._wp )                                              & 
699                        &             ) * 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
700                        &           )   * zmsk00y(ji,jj)
701                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
702                                     !---> MV : comment above is misleading as EVP is a purely explicit method
703                                     !---> MV : I found code below quite unreadable :-)
704                                     !---> MV : two maskings + landfast / no landfast options, I think this is a bit too much
705                                     !---> MV : I would suggest to split at least the landfast / masking steps if possible
706                                     !---> MV : Should also comment why 1% of ocean velocity is assumed - not sure it makes sense btw
707                                     !---> MV : I would say 100% of ocean current makes much more sense for free-drift
708                                     !---> MV : of very low concentration sea ice..
709                                     !---> MV : but maybe it is a numerical trick... in brief it's n
710                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
711                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
712                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
713                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
714                        &             ) * 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
715                        &           )   * zmsk00y(ji,jj)
716                  ENDIF
717               END DO
718            END DO
719            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
720            !
721#if defined key_agrif
722!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
723            CALL agrif_interp_ice( 'V' )
724#endif
725            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
726            !
727         ENDIF
728
729         ! convergence test
730         IF( ln_rhg_chkcvg )   CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
731         !
732         !                                                ! ==================== !
733      END DO                                              !  end loop over jter  !
734      !                                                   ! ==================== !
735      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
736      !
737      !------------------------------------------------------------------------------!
738      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
739      !------------------------------------------------------------------------------!
740      DO jj = 1, jpjm1
741         DO ji = 1, jpim1
742
743            ! shear at F points
744            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)   &
745               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
746               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
747
748         END DO
749      END DO           
750     
751      DO jj = 2, jpjm1
752         DO ji = 2, jpim1 ! no vector loop
753           
754            ! tension**2 at T points
755            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)   &
756               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
757               &   ) * r1_e1e2t(ji,jj)
758            zdt2 = zdt * zdt
759           
760            zten_i(ji,jj) = zdt
761           
762            ! shear**2 at T points (doc eq. A16)
763            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
764               &   + 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)  &
765               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
766           
767            ! shear at T points
768            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
769
770            ! divergence at T points
771            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
772               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
773               &             ) * r1_e1e2t(ji,jj)
774           
775            ! delta at T points
776            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
777            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
778            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
779
780         END DO
781      END DO
782      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
783     
784      ! --- Store the stress tensor for the next time step --- !
785      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
786      pstress1_i (:,:) = zs1 (:,:)
787      pstress2_i (:,:) = zs2 (:,:)
788      pstress12_i(:,:) = zs12(:,:)
789      !
790
791      !------------------------------------------------------------------------------!
792      ! 5) diagnostics
793      !------------------------------------------------------------------------------!
794      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
795      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
796         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
797         !
798         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
799            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
800         !
801         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
802         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
803         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
804         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
805         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
806         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
807      ENDIF
808       
809      ! --- divergence, shear and strength --- !
810      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
811      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
812      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
813      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
814
815      ! --- Stress tensor invariants (SIMIP diags) --- !
816      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
817         !
818         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
819         !         
820         DO jj = 2, jpj - 1
821            DO ji = 2, jpi - 1
822           
823               ! Ice stresses
824               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
825               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
826               ! I know, this can be confusing...
827               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
828               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
829               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
830               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
831               
832               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
833               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st stress invariant, aka average normal stress, aka negative pressure
834               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )   ! 2nd  ''       '', aka maximum shear stress
835               
836            END DO
837         END DO
838
839         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
840         
841         !
842         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
843         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
844         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
845         
846         DEALLOCATE ( zsig_I, zsig_II )
847         
848      ENDIF
849
850      ! --- Normalized stress tensor principal components --- !
851      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
852      ! Recommendation 1 : we use ice strength, not replacement pressure
853      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
854      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
855         !
856         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
857         !         
858         DO jj = 2, jpj - 1
859            DO ji = 2, jpi - 1
860           
861               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
862               !                        and **deformations** at current iterates
863               !                        following Lemieux & Dupont (2020)
864               zfac             =   zp_delt(ji,jj)
865               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
866               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
867               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
868               
869               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
870               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure
871               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )     ! 2nd  ''       '', aka maximum shear stress
872
873               ! Normalized  principal stresses (used to display the ellipse)
874               z1_strength      =   1._wp / MAX ( 1.0, strength(ji,jj) )
875               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
876               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
877            END DO
878         END DO
879         
880         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
881               
882         !
883         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
884         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
885
886         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
887         
888      ENDIF
889      ! --- SIMIP --- !
890      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
891         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
892         !
893         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., &
894            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
895
896         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
897         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
898         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
899         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
900         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
901         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
902      ENDIF
903
904      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
905         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
906         !
907         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
908            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
909         !
910         DO jj = 2, jpjm1
911            DO ji = 2, jpim1
912               ! 2D ice mass, snow mass, area transport arrays (X, Y)
913               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
914               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
915
916               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
917               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
918
919               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
920               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
921
922               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
923               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
924
925            END DO
926         END DO
927
928         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
929            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
930            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
931
932         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
933         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
934         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
935         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
936         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
937         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
938
939         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
940            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
941
942      ENDIF
943      !
944      ! --- convergence tests --- !
945      IF( ln_rhg_chkcvg ) THEN
946         IF( iom_use('uice_cvg') ) THEN
947            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
948               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
949                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
950            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
951               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
952                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
953            ENDIF
954         ENDIF
955      ENDIF     
956      !
957      DEALLOCATE( zmsk00, zmsk15 )
958      !
959   END SUBROUTINE ice_dyn_rhg_evp
960
961
962   SUBROUTINE rhg_cvg_evp( kt, kiter, kitermax, pu, pv, pub, pvb )
963      !!----------------------------------------------------------------------
964      !!                    ***  ROUTINE rhg_cvg_evp  ***
965      !!                     
966      !! ** Purpose :   check convergence of evp ice rheology
967      !!
968      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
969      !!                during the sub timestepping of rheology so as:
970      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
971      !!                This routine is called every sub-iteration, so it is cpu expensive
972      !!
973      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
974      !!----------------------------------------------------------------------
975      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
976      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
977      !!
978      INTEGER           ::   it, idtime, istatus
979      INTEGER           ::   ji, jj          ! dummy loop indices
980      REAL(wp)          ::   zresm           ! local real
981      CHARACTER(len=20) ::   clname
982      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
983      !!----------------------------------------------------------------------
984
985      ! create file
986      IF( kt == nit000 .AND. kiter == 1 ) THEN
987         !
988         IF( lwp ) THEN
989            WRITE(numout,*)
990            WRITE(numout,*) 'rhg_cvg_evp : ice rheology convergence control'
991            WRITE(numout,*) '~~~~~~~~~~~'
992         ENDIF
993         !
994         IF( lwm ) THEN
995            clname = 'ice_cvg.nc'
996            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
997            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
998            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
999            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
1000            istatus = NF90_ENDDEF(ncvgid)
1001         ENDIF
1002         !
1003      ENDIF
1004
1005      ! time
1006      it = ( kt - 1 ) * kitermax + kiter
1007     
1008      ! convergence
1009      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
1010         zresm = 0._wp
1011      ELSE
1012         DO jj = 1, jpj
1013            DO ji = 1, jpi
1014               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
1015                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
1016            END DO
1017         END DO
1018         zresm = MAXVAL( zres )
1019         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
1020      ENDIF
1021
1022      IF( lwm ) THEN
1023         ! write variables
1024         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
1025         ! close file
1026         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
1027      ENDIF
1028     
1029   END SUBROUTINE rhg_cvg_evp
1030
1031
1032   SUBROUTINE rhg_evp_rst( cdrw, kt )
1033      !!---------------------------------------------------------------------
1034      !!                   ***  ROUTINE rhg_evp_rst  ***
1035      !!                     
1036      !! ** Purpose :   Read or write RHG file in restart file
1037      !!
1038      !! ** Method  :   use of IOM library
1039      !!----------------------------------------------------------------------
1040      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1041      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
1042      !
1043      INTEGER  ::   iter            ! local integer
1044      INTEGER  ::   id1, id2, id3   ! local integers
1045      !!----------------------------------------------------------------------
1046      !
1047      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
1048         !                                   ! ---------------
1049         IF( ln_rstart ) THEN                   !* Read the restart file
1050            !
1051            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
1052            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
1053            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
1054            !
1055            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
1056               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
1057               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
1058               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
1059            ELSE                                     ! start rheology from rest
1060               IF(lwp) WRITE(numout,*)
1061               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
1062               stress1_i (:,:) = 0._wp
1063               stress2_i (:,:) = 0._wp
1064               stress12_i(:,:) = 0._wp
1065            ENDIF
1066         ELSE                                   !* Start from rest
1067            IF(lwp) WRITE(numout,*)
1068            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
1069            stress1_i (:,:) = 0._wp
1070            stress2_i (:,:) = 0._wp
1071            stress12_i(:,:) = 0._wp
1072         ENDIF
1073         !
1074      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1075         !                                   ! -------------------
1076         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
1077         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1078         !
1079         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
1080         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
1081         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
1082         !
1083      ENDIF
1084      !
1085   END SUBROUTINE rhg_evp_rst
1086
1087   
1088#else
1089   !!----------------------------------------------------------------------
1090   !!   Default option         Empty module           NO SI3 sea-ice model
1091   !!----------------------------------------------------------------------
1092#endif
1093
1094   !!==============================================================================
1095END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.