source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 13550

Last change on this file since 13550 was 13550, checked in by clem, 5 months ago

trunk: remove one communication in rheology, so about 100 overall.

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