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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_evp.F90 in NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 @ 12263

Last change on this file since 12263 was 12263, checked in by stefryn, 4 years ago

add test case for rheology, will change lateral bounday conditions later

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