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

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

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

Last change on this file since 13286 was 13286, checked in by smasson, 4 years ago

trunk: merge extra halos branch in trunk, see #2366

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