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

Last change on this file since 10413 was 10413, checked in by clem, 2 years ago

merge dev_r9947_SI3_advection with the trunk. All sette tests passed. There is probably a conservation issue with the new advection scheme that I should solve but the routine icedyn_adv_umx.F90 is going to change anyway in the next couple of weeks. At worst, the old routine can be plugged without harm

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