source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/ICE/icedyn_rhg_evp.F90 @ 11869

Last change on this file since 11869 was 11869, checked in by acc, 11 months ago

Branch 2019/fix_sn_cfctl_ticket2328. Changes to enable correct functionality for the sn_cfctl%l_mppout and sn_cfctl%l_mpptop options. These changes also introduce a sn_cfctl%l_oasout option to toggle the OASIS setup information (sbccpl.F90, only) which was yet another misuse of ln_ctl. The next step may be to remove most references to ln_ctl altogether and provide a single control mechanism. TBD. See ticket #2328

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