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 branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90 @ 8752

Last change on this file since 8752 was 8752, checked in by dancopsey, 6 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

File size: 55.3 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) 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_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                       ESIM sea-ice model
17   !!----------------------------------------------------------------------
18   !!   ice_dyn_rhg_evp       : computes ice velocities from EVP rheology
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE dom_oce        ! Ocean domain
22   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
23   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
24   USE ice            ! sea-ice: ice variables
25   USE icedyn_rdgrft  ! sea-ice: ice strength
26   !
27   USE in_out_manager ! I/O manager
28   USE iom            ! I/O manager library
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
31   USE lbclnk         ! lateral boundary conditions (or mpp links)
32   USE prtctl         ! Print control
33   !
34#if defined key_agrif
35   USE agrif_lim3_interp
36#endif
37   USE bdy_oce   , ONLY: ln_bdy 
38   USE bdyice 
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#if defined key_agrif 
180      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
181      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
182#endif
183      !
184      !------------------------------------------------------------------------------!
185      ! 0) mask at F points for the ice
186      !------------------------------------------------------------------------------!
187      ! ocean/land mask
188      DO jj = 1, jpjm1
189         DO ji = 1, jpim1      ! NO vector opt.
190            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
191         END DO
192      END DO
193      CALL lbc_lnk( zfmask, 'F', 1._wp )
194
195      ! Lateral boundary conditions on velocity (modify zfmask)
196      zwf(:,:) = zfmask(:,:)
197      DO jj = 2, jpjm1
198         DO ji = fs_2, fs_jpim1   ! vector opt.
199            IF( zfmask(ji,jj) == 0._wp ) THEN
200               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) ) )
201            ENDIF
202         END DO
203      END DO
204      DO jj = 2, jpjm1
205         IF( zfmask(1,jj) == 0._wp ) THEN
206            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
207         ENDIF
208         IF( zfmask(jpi,jj) == 0._wp ) THEN
209            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
210         ENDIF
211      END DO
212      DO ji = 2, jpim1
213         IF( zfmask(ji,1) == 0._wp ) THEN
214            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
215         ENDIF
216         IF( zfmask(ji,jpj) == 0._wp ) THEN
217            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
218         ENDIF
219      END DO
220      CALL lbc_lnk( zfmask, 'F', 1._wp )
221
222      !------------------------------------------------------------------------------!
223      ! 1) define some variables and initialize arrays
224      !------------------------------------------------------------------------------!
225      zrhoco = rau0 * rn_cio 
226
227      ! ecc2: square of yield ellipse eccenticrity
228      ecc2    = rn_ecc * rn_ecc
229      z1_ecc2 = 1._wp / ecc2
230
231      ! Time step for subcycling
232      zdtevp   = rdt_ice / REAL( nn_nevp )
233      z1_dtevp = 1._wp / zdtevp
234
235      ! alpha parameters (Bouillon 2009)
236      IF( .NOT. ln_aEVP ) THEN
237         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
238         zalph2 = zalph1 * z1_ecc2
239
240         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
241         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
242      ENDIF
243         
244      ! Initialise stress tensor
245      zs1 (:,:) = pstress1_i (:,:) 
246      zs2 (:,:) = pstress2_i (:,:)
247      zs12(:,:) = pstress12_i(:,:)
248
249      ! Ice strength
250      CALL ice_strength
251
252      ! scale factors
253      DO jj = 2, jpjm1
254         DO ji = fs_2, fs_jpim1
255            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
256            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
257         END DO
258      END DO
259           
260      !
261      !------------------------------------------------------------------------------!
262      ! 2) Wind / ocean stress, mass terms, coriolis terms
263      !------------------------------------------------------------------------------!
264
265      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
266         !                                           
267         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
268         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
269         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
270         !
271         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
272         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
273         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
274         !
275         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
276         !
277      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
278         zpice(:,:) = ssh_m(:,:)
279      ENDIF
280
281      DO jj = 2, jpjm1
282         DO ji = fs_2, fs_jpim1
283
284            ! ice fraction at U-V points
285            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)
286            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)
287
288            ! Ice/snow mass at U-V points
289            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
290            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
291            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
292            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
293            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
294
295            ! Ocean currents at U-V points
296            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
297               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
298           
299            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
300               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
301
302            ! Coriolis at T points (m*f)
303            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
304
305            ! dt/m at T points (for alpha and beta coefficients)
306            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
307           
308            ! m/dt
309            zmU_t(ji,jj)    = zmassU * z1_dtevp
310            zmV_t(ji,jj)    = zmassV * z1_dtevp
311           
312            ! Drag ice-atm.
313            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
314            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
315
316            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
317            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
318            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
319
320            ! masks
321            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
322            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
323
324            ! switches
325            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
326            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
327
328         END DO
329      END DO
330      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
331      !
332      !------------------------------------------------------------------------------!
333      ! 3) Solution of the momentum equation, iterative procedure
334      !------------------------------------------------------------------------------!
335      !
336      !                                               !----------------------!
337      DO jter = 1 , nn_nevp                           !    loop over jter    !
338         !                                            !----------------------!       
339         IF(ln_ctl) THEN   ! Convergence test
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
345
346         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
347         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
348            DO ji = 1, jpim1
349
350               ! shear at F points
351               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)   &
352                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
353                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
354
355            END DO
356         END DO
357         CALL lbc_lnk( zds, 'F', 1. )
358
359         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
360            DO ji = 2, jpi ! no vector loop
361
362               ! shear**2 at T points (doc eq. A16)
363               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
364                  &   + 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)  &
365                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
366             
367               ! divergence at T points
368               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
369                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
370                  &    ) * r1_e1e2t(ji,jj)
371               zdiv2 = zdiv * zdiv
372               
373               ! tension at T points
374               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)   &
375                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
376                  &   ) * r1_e1e2t(ji,jj)
377               zdt2 = zdt * zdt
378               
379               ! delta at T points
380               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
381
382               ! P/delta at T points
383               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
384
385               ! alpha & beta for aEVP
386               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
387               !   alpha = beta = sqrt(4*gamma)
388               IF( ln_aEVP ) THEN
389                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
390                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
391                  zalph2   = zalph1
392                  z1_alph2 = z1_alph1
393               ENDIF
394               
395               ! stress at T points
396               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
397               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
398             
399            END DO
400         END DO
401         CALL lbc_lnk( zp_delt, 'T', 1. )
402
403         DO jj = 1, jpjm1
404            DO ji = 1, jpim1
405
406               ! alpha & beta for aEVP
407               IF( ln_aEVP ) THEN
408                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
409                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
410                  zbeta(ji,jj) = zalph2
411               ENDIF
412               
413               ! P/delta at F points
414               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) )
415               
416               ! stress at F points
417               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
418
419            END DO
420         END DO
421
422         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
423         DO jj = 2, jpjm1
424            DO ji = fs_2, fs_jpim1               
425               !                   !--- U points
426               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
427                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
428                  &                    ) * r1_e2u(ji,jj)                                                                      &
429                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
430                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
431                  &                  ) * r1_e1e2u(ji,jj)
432               !
433               !                !--- V points
434               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
435                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
436                  &                    ) * r1_e1v(ji,jj)                                                                      &
437                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
438                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
439                  &                  ) * r1_e1e2v(ji,jj)
440               !
441               !                !--- u_ice at V point
442               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
443                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
444               !
445               !                !--- v_ice at U point
446               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
447                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
448            END DO
449         END DO
450         !
451         ! --- Computation of ice velocity --- !
452         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
453         !  Bouillon et al. 2009 (eq 34-35) => stable
454         IF( MOD(jter,2) == 0 ) THEN ! even iterations
455            !
456            DO jj = 2, jpjm1
457               DO ji = fs_2, fs_jpim1
458                  !                 !--- tau_io/(v_oce - v_ice)
459                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
460                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
461                  !                 !--- Ocean-to-Ice stress
462                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
463                  !
464                  !                 !--- tau_bottom/v_ice
465                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
466                  zTauB = - tau_icebfr(ji,jj) / zvel
467                  !
468                  !                 !--- Coriolis at V-points (energy conserving formulation)
469                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
470                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
471                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
472                  !
473                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
474                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
475                  !
476                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
477                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
478                  !
479                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
480                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
481                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
482                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
483                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
484                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
485                     &           ) * zmaskV(ji,jj)
486                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
487                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
488                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
489                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
490                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
491                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
492                     &            ) * zmaskV(ji,jj)
493                  ENDIF
494               END DO
495            END DO
496            CALL lbc_lnk( v_ice, 'V', -1. )
497            !
498#if defined key_agrif
499!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
500            CALL agrif_interp_lim3( 'V' )
501#endif
502            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
503            !
504            DO jj = 2, jpjm1
505               DO ji = fs_2, fs_jpim1         
506                  !                 !--- tau_io/(u_oce - u_ice)
507                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
508                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
509                  !                 !--- Ocean-to-Ice stress
510                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
511                  !
512                  !                 !--- tau_bottom/u_ice
513                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
514                  zTauB = - tau_icebfr(ji,jj) / zvel
515                  !
516                  !                 !--- Coriolis at U-points (energy conserving formulation)
517                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
518                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
519                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
520                  !
521                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
522                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
523                  !
524                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
525                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
526                  !
527                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
528                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
529                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
530                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
531                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
532                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
533                     &            ) * zmaskU(ji,jj)
534                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
535                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
536                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
537                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
538                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
539                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
540                     &            ) * zmaskU(ji,jj)
541                  ENDIF
542               END DO
543            END DO
544            CALL lbc_lnk( u_ice, 'U', -1. )
545            !
546#if defined key_agrif
547!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
548            CALL agrif_interp_lim3( 'U' )
549#endif
550            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
551            !
552         ELSE ! odd iterations
553            !
554            DO jj = 2, jpjm1
555               DO ji = fs_2, fs_jpim1
556                  !                 !--- tau_io/(u_oce - u_ice)
557                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
558                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
559                  !                 !--- Ocean-to-Ice stress
560                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
561                  !
562                  !                 !--- tau_bottom/u_ice
563                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
564                  zTauB = - tau_icebfr(ji,jj) / zvel
565                  !
566                  !                 !--- Coriolis at U-points (energy conserving formulation)
567                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
568                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
569                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
570                  !
571                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
572                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
573                  !
574                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
575                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
576                  !
577                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
578                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
579                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
580                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + 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                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
583                     &            ) * zmaskU(ji,jj)
584                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
585                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
586                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
587                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
588                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
589                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
590                     &            ) * zmaskU(ji,jj)
591                  ENDIF
592               END DO
593            END DO
594            CALL lbc_lnk( u_ice, 'U', -1. )
595            !
596#if defined key_agrif
597!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
598            CALL agrif_interp_lim3( 'U' )
599#endif
600            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
601            !
602            DO jj = 2, jpjm1
603               DO ji = fs_2, fs_jpim1
604                  !                 !--- tau_io/(v_oce - v_ice)
605                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
606                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
607                  !                 !--- Ocean-to-Ice stress
608                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
609                  !
610                  !                 !--- tau_bottom/v_ice
611                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
612                  ztauB = - tau_icebfr(ji,jj) / zvel
613                  !
614                  !                 !--- Coriolis at v-points (energy conserving formulation)
615                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
616                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
617                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
618                  !
619                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
620                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
621                  !
622                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
623                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
624                  !
625                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
626                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
627                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
628                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
629                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
630                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
631                     &           ) * zmaskV(ji,jj)
632                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
633                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
634                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
635                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
636                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
637                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
638                     &            ) * zmaskV(ji,jj)
639                  ENDIF
640               END DO
641            END DO
642            CALL lbc_lnk( v_ice, 'V', -1. )
643            !
644#if defined key_agrif
645!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
646            CALL agrif_interp_lim3( 'V' )
647#endif
648            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
649            !
650         ENDIF
651         
652         IF(ln_ctl) THEN   ! Convergence test
653            DO jj = 2 , jpjm1
654               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
655            END DO
656            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
657            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
658         ENDIF
659         !
660         !                                                ! ==================== !
661      END DO                                              !  end loop over jter  !
662      !                                                   ! ==================== !
663      !
664      !------------------------------------------------------------------------------!
665      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
666      !------------------------------------------------------------------------------!
667      DO jj = 1, jpjm1
668         DO ji = 1, jpim1
669
670            ! shear at F points
671            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)   &
672               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
673               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
674
675         END DO
676      END DO           
677     
678      DO jj = 2, jpjm1
679         DO ji = 2, jpim1 ! no vector loop
680           
681            ! tension**2 at T points
682            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)   &
683               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
684               &   ) * r1_e1e2t(ji,jj)
685            zdt2 = zdt * zdt
686           
687            ! shear**2 at T points (doc eq. A16)
688            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
689               &   + 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)  &
690               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
691           
692            ! shear at T points
693            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
694
695            ! divergence at T points
696            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
697               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
698               &             ) * r1_e1e2t(ji,jj)
699           
700            ! delta at T points
701            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
702            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
703            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
704
705         END DO
706      END DO
707      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
708     
709      ! --- Store the stress tensor for the next time step --- !
710      CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
711      pstress1_i (:,:) = zs1 (:,:)
712      pstress2_i (:,:) = zs2 (:,:)
713      pstress12_i(:,:) = zs12(:,:)
714      !
715
716      !------------------------------------------------------------------------------!
717      ! 5) diagnostics
718      !------------------------------------------------------------------------------!
719      DO jj = 1, jpj
720         DO ji = 1, jpi
721            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
722         END DO
723      END DO
724
725      ! --- divergence, shear and strength --- !
726      IF( iom_use('idive')  )   CALL iom_put( "idive"  , pdivu_i(:,:)  * zswi(:,:) )   ! divergence
727      IF( iom_use('ishear') )   CALL iom_put( "ishear" , pshear_i(:,:) * zswi(:,:) )   ! shear
728      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
729
730      ! --- charge ellipse --- !
731      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
732         !
733         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
734         !         
735         DO jj = 2, jpjm1
736            DO ji = 2, jpim1
737               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
738                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
739                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
740
741               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
742
743               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
744
745!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
746!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
747!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
748!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
749               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
750               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
751               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
752            END DO
753         END DO
754         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
755         !
756         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
757         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
758         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
759         !
760         DEALLOCATE( zsig1 , zsig2 , zsig3 )
761      ENDIF
762     
763      ! --- SIMIP --- !
764      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
765         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
766         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
767         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
768
769         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
770            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
771            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
772            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
773         
774         DO jj = 2, jpjm1
775            DO ji = 2, jpim1
776               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
777               
778               ! Stress tensor invariants (normal and shear stress N/m)
779               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
780               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
781               
782               ! Stress terms of the momentum equation (N/m2)
783               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
784               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
785               
786               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
787               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
788               
789               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
790               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
791               
792               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
793               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
794               
795               ! 2D ice mass, snow mass, area transport arrays (X, Y)
796               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
797               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
798               
799               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
800               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
801               
802               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
803               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
804               
805               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
806               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
807               
808            END DO
809         END DO
810         
811         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
812            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
813            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
814            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
815         
816         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    )
817         
818         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
819            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
820            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
821         
822         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
823         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
824         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
825         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
826         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
827         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
828         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
829         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
830         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
831         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
832         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
833         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
834         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
835         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
836         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
837         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
838
839         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
840            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
841            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
842            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
843
844      ENDIF
845      !
846   END SUBROUTINE ice_dyn_rhg_evp
847
848   SUBROUTINE rhg_evp_rst( cdrw, kt )
849      !!---------------------------------------------------------------------
850      !!                   ***  ROUTINE rhg_evp_rst  ***
851      !!                     
852      !! ** Purpose :   Read or write RHG file in restart file
853      !!
854      !! ** Method  :   use of IOM library
855      !!----------------------------------------------------------------------
856      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
857      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
858      !
859      INTEGER  ::   iter            ! local integer
860      INTEGER  ::   id1, id2, id3   ! local integers
861      !!----------------------------------------------------------------------
862      !
863      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
864         !                                   ! ---------------
865         IF( ln_rstart ) THEN                   !* Read the restart file
866            !
867            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
868            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
869            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
870            !
871            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
872               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
873               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
874               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
875            ELSE                                     ! start rheology from rest
876               IF(lwp) WRITE(numout,*) '   ==>>   previous run without rheology, set stresses to 0'
877               stress1_i (:,:) = 0._wp
878               stress2_i (:,:) = 0._wp
879               stress12_i(:,:) = 0._wp
880            ENDIF
881         ELSE                                   !* Start from rest
882            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set stresses to 0'
883            stress1_i (:,:) = 0._wp
884            stress2_i (:,:) = 0._wp
885            stress12_i(:,:) = 0._wp
886         ENDIF
887         !
888      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
889         !                                   ! -------------------
890         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
891         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
892         !
893         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
894         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
895         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
896         !
897      ENDIF
898      !
899   END SUBROUTINE rhg_evp_rst
900
901#else
902   !!----------------------------------------------------------------------
903   !!   Default option         Empty module           NO ESIM sea-ice model
904   !!----------------------------------------------------------------------
905#endif
906
907   !!==============================================================================
908END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.