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 @ 9596

Last change on this file since 9596 was 9596, checked in by nicolasmartin, 6 years ago

Reorganisation plan for NEMO repository: preliminary implementation of the new trunk structure (ref #2015)

  • Deletion of useless 1st level, add new explicit names for folders
  • Extract dev tools (SETTE and Trusting) to /utils/CI
  • Extract figures folder from documentation to /utils/figures
  • Extrac TEST_CASES to /utils/test_cases
  • Move ./TOOLS/COMPILE to ./mk
  • Add few wiki files at the root: README, INSTALL and LICENSE with full text of CeCILL v2.0
  • Style: uppercase for wiki files, lowercase for folder names

Additionally, inclusion of routines renaming in NST: agrif_si3 -> agrif_ice

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