source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icedyn_rhg_evp.F90 @ 12808

Last change on this file since 12808 was 12489, checked in by davestorkey, 9 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

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