source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 7256

Last change on this file since 7256 was 7256, checked in by cbricaud, 4 years ago

phaze NEMO routines in CRS branch with nemo_v3_6_STABLE branch at rev 7213 (09-09-2016) (merge -r 5519:7213 )

  • Property svn:keywords set to Id
File size: 35.3 KB
Line 
1MODULE limrhg
2   !!======================================================================
3   !!                     ***  MODULE  limrhg  ***
4   !!   Ice rheology : sea ice rheology
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
10   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
12   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy)
13   !!----------------------------------------------------------------------
14#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'               OR                     LIM-3 sea-ice model
17   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
18   !!----------------------------------------------------------------------
19   !!   lim_rhg       : computes ice velocities
20   !!----------------------------------------------------------------------
21   USE phycst         ! Physical constant
22   USE oce     , ONLY :  snwice_mass, snwice_mass_b
23   USE par_oce        ! Ocean parameters
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce        ! Surface boundary condition: ocean fields
26   USE sbc_ice        ! Surface boundary condition: ice fields
27#if defined key_lim3
28   USE ice            ! LIM-3: ice variables
29   USE dom_ice        ! LIM-3: ice domain
30   USE limitd_me      ! LIM-3:
31#else
32   USE ice_2          ! LIM-2: ice variables
33   USE dom_ice_2      ! LIM-2: ice domain
34#endif
35   USE lbclnk         ! Lateral Boundary Condition / MPP link
36   USE lib_mpp        ! MPP library
37   USE wrk_nemo       ! work arrays
38   USE in_out_manager ! I/O manager
39   USE prtctl         ! Print control
40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
41#if defined key_agrif && defined key_lim2
42   USE agrif_lim2_interp
43#endif
44#if defined key_bdy
45   USE bdyice_lim
46#endif
47
48   IMPLICIT NONE
49   PRIVATE
50
51   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
52
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE lim_rhg( k_j1, k_jpj )
63      !!-------------------------------------------------------------------
64      !!                 ***  SUBROUTINE lim_rhg  ***
65      !!                          EVP-C-grid
66      !!
67      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
68      !!  stress and sea-surface slope. Ice-ice interaction is described by
69      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
70      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
71      !!
72      !!  The points in the C-grid look like this, dear reader
73      !!
74      !!                              (ji,jj)
75      !!                                 |
76      !!                                 |
77      !!                      (ji-1,jj)  |  (ji,jj)
78      !!                             ---------   
79      !!                            |         |
80      !!                            | (ji,jj) |------(ji,jj)
81      !!                            |         |
82      !!                             ---------   
83      !!                     (ji-1,jj-1)     (ji,jj-1)
84      !!
85      !! ** Inputs  : - wind forcing (stress), oceanic currents
86      !!                ice total volume (vt_i) per unit area
87      !!                snow total volume (vt_s) per unit area
88      !!
89      !! ** Action  : - compute u_ice, v_ice : the components of the
90      !!                sea-ice velocity vector
91      !!              - compute delta_i, shear_i, divu_i, which are inputs
92      !!                of the ice thickness distribution
93      !!
94      !! ** Steps   : 1) Compute ice snow mass, ice strength
95      !!              2) Compute wind, oceanic stresses, mass terms and
96      !!                 coriolis terms of the momentum equation
97      !!              3) Solve the momentum equation (iterative procedure)
98      !!              4) Recompute invariants of the strain rate tensor
99      !!                 which are inputs of the ITD, store stress
100      !!                 for the next time step
101      !!              5) Control prints of residual (convergence)
102      !!                 and charge ellipse.
103      !!                 The user should make sure that the parameters
104      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
105      !!                 on the charge ellipse for plastic flow
106      !!                 e.g. in the Canadian Archipelago
107      !!
108      !! ** Notes   : Boundary condition for ice is chosen no-slip
109      !!              but can be adjusted with param rn_shlat
110      !!
111      !! References : Hunke and Dukowicz, JPO97
112      !!              Bouillon et al., Ocean Modelling 2009
113      !!-------------------------------------------------------------------
114      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
115      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
116      !!
117      INTEGER ::   ji, jj   ! dummy loop indices
118      INTEGER ::   jter     ! local integers
119      CHARACTER (len=50) ::   charout
120
121      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
122      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
123      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
124      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
125      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
126      REAL(wp) ::   zTauO, zTauE, zCor                                       ! temporary scalars
127
128      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
129      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
130      REAL(wp) ::   zintb, zintn                                             ! dummy argument
131     
132      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength
133      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors
134      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points
135      !
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
138      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points
139      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
140      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
141      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
142      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses
143     
144      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear
145      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components
146      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence
147      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope:
148                                                                             !   ocean surface (ssh_m) if ice is not embedded
149                                                                             !   ice top surface if ice is embedded   
150      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
151      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
152      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
153
154      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
155      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
156      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
157      !!-------------------------------------------------------------------
158
159      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
160      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
161      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
162      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
163      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
164
165#if  defined key_lim2 && ! defined key_lim2_vp
166# if defined key_agrif
167      USE ice_2, vt_s => hsnm
168      USE ice_2, vt_i => hicm
169# else
170      vt_s => hsnm
171      vt_i => hicm
172# endif
173      at_i(:,:) = 1. - frld(:,:)
174#endif
175#if defined key_agrif && defined key_lim2 
176      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
177#endif
178      !
179      !------------------------------------------------------------------------------!
180      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj)
181      !------------------------------------------------------------------------------!
182      ! ocean/land mask
183      DO jj = 1, jpjm1
184         DO ji = 1, jpim1      ! NO vector opt.
185            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
186         END DO
187      END DO
188      CALL lbc_lnk( zfmask, 'F', 1._wp )
189
190      ! Lateral boundary conditions on velocity (modify zfmask)
191      zwf(:,:) = zfmask(:,:)
192      DO jj = 2, jpjm1
193         DO ji = fs_2, fs_jpim1   ! vector opt.
194            IF( zfmask(ji,jj) == 0._wp ) THEN
195               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
196            ENDIF
197         END DO
198      END DO
199      DO jj = 2, jpjm1
200         IF( zfmask(1,jj) == 0._wp ) THEN
201            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
202         ENDIF
203         IF( zfmask(jpi,jj) == 0._wp ) THEN
204            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
205         ENDIF
206      END DO
207      DO ji = 2, jpim1
208         IF( zfmask(ji,1) == 0._wp ) THEN
209            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
210         ENDIF
211         IF( zfmask(ji,jpj) == 0._wp ) THEN
212            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
213         ENDIF
214      END DO
215      CALL lbc_lnk( zfmask, 'F', 1._wp )
216
217      !------------------------------------------------------------------------------!
218      ! 1) define some variables and initialize arrays
219      !------------------------------------------------------------------------------!
220      ! ecc2: square of yield ellipse eccenticrity
221      ecc2    = rn_ecc * rn_ecc
222      z1_ecc2 = 1._wp / ecc2
223
224      ! Time step for subcycling
225      zdtevp   = rdt_ice / REAL( nn_nevp )
226      z1_dtevp = 1._wp / zdtevp
227
228      ! alpha parameters (Bouillon 2009)
229#if defined key_lim3
230      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
231#else
232      zalph1 = ( 2._wp * telast ) * z1_dtevp
233#endif
234      zalph2 = zalph1 * z1_ecc2
235
236      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
237      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
238
239      ! Initialise stress tensor
240      zs1 (:,:) = stress1_i (:,:) 
241      zs2 (:,:) = stress2_i (:,:)
242      zs12(:,:) = stress12_i(:,:)
243
244      ! Ice strength
245#if defined key_lim3
246      CALL lim_itd_me_icestrength( nn_icestr )
247      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
248#else
249      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) )
250#endif
251
252      ! scale factors
253      DO jj = k_j1+1, k_jpj-1
254         DO ji = fs_2, fs_jpim1
255            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
256            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
257         END DO
258      END DO
259           
260      !
261      !------------------------------------------------------------------------------!
262      ! 2) Wind / ocean stress, mass terms, coriolis terms
263      !------------------------------------------------------------------------------!
264
265      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
266         !                                           
267         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
268         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
269         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
270         !
271         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
272         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
273         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
274         !
275         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
276         !
277      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
278         zpice(:,:) = ssh_m(:,:)
279      ENDIF
280
281      DO jj = k_j1+1, k_jpj-1
282         DO ji = fs_2, fs_jpim1
283
284            ! ice fraction at U-V points
285            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
286            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
287
288            ! Ice/snow mass at U-V points
289            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
290            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
291            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
292            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
293            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
294
295            ! Ocean currents at U-V points
296            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
297               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
298           
299            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
300               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
301
302            ! Coriolis at T points (m*f)
303            zmf(ji,jj)      = zm1 * fcor(ji,jj)
304
305            ! m/dt
306            zmU_t(ji,jj)    = zmassU * z1_dtevp
307            zmV_t(ji,jj)    = zmassV * z1_dtevp
308
309            ! Drag ice-atm.
310            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
311            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
312
313            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
314            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
315            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
316
317            ! masks
318            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
319            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
320
321            ! switches
322            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
323            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
324
325         END DO
326      END DO
327      CALL lbc_lnk( zmf, 'T', 1. )
328      !
329      !------------------------------------------------------------------------------!
330      ! 3) Solution of the momentum equation, iterative procedure
331      !------------------------------------------------------------------------------!
332      !
333      !                                               !----------------------!
334      DO jter = 1 , nn_nevp                           !    loop over jter    !
335         !                                            !----------------------!       
336         IF(ln_ctl) THEN   ! Convergence test
337            DO jj = k_j1, k_jpj-1
338               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
339               zv_ice(:,jj) = v_ice(:,jj)
340            END DO
341         ENDIF
342
343         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
344         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
345            DO ji = 1, jpim1
346
347               ! shear at F points
348               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)   &
349                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
350                  &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)
351
352            END DO
353         END DO
354         CALL lbc_lnk( zds, 'F', 1. )
355
356         DO jj = k_j1+1, k_jpj-1
357            DO ji = 2, jpim1 ! no vector loop
358
359               ! shear**2 at T points (doc eq. A16)
360               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
361                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
362                  &   ) * 0.25_wp * r1_e12t(ji,jj)
363             
364               ! divergence at T points
365               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
366                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
367                  &    ) * r1_e12t(ji,jj)
368               zdiv2 = zdiv * zdiv
369               
370               ! tension at T points
371               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)   &
372                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
373                  &   ) * r1_e12t(ji,jj)
374               zdt2 = zdt * zdt
375               
376               ! delta at T points
377               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 ) 
378
379               ! P/delta at T points
380               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
381               
382               ! stress at T points
383               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
384               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
385             
386            END DO
387         END DO
388         CALL lbc_lnk( zp_delt, 'T', 1. )
389
390         DO jj = k_j1, k_jpj-1
391            DO ji = 1, jpim1
392
393               ! P/delta at F points
394               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) )
395               
396               ! stress at F points
397               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
398
399            END DO
400         END DO
401         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
402 
403         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
404         DO jj = k_j1+1, k_jpj-1
405            DO ji = fs_2, fs_jpim1               
406
407               ! U points
408               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
409                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
410                  &                    ) * r1_e2u(ji,jj)                                                                      &
411                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
412                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
413                  &                  ) * r1_e12u(ji,jj)
414
415               ! V points
416               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
417                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
418                  &                    ) * r1_e1v(ji,jj)                                                                      &
419                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
420                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
421                  &                  ) * r1_e12v(ji,jj)
422
423               ! u_ice at V point
424               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
425                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
426               
427               ! v_ice at U point
428               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
429                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
430
431            END DO
432         END DO
433         !
434         ! --- Computation of ice velocity --- !
435         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
436         !  Bouillon et al. 2009 (eq 34-35) => stable
437         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
438           
439            DO jj = k_j1+1, k_jpj-1
440               DO ji = fs_2, fs_jpim1
441
442                  ! tau_io/(v_oce - v_ice)
443                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
444                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
445
446                  ! Coriolis at V-points (energy conserving formulation)
447                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
448                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
449                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
450
451                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
452                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
453                 
454                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
455                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
456                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
457                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
458                     &           ) * zmaskV(ji,jj)
459               END DO
460            END DO
461            CALL lbc_lnk( v_ice, 'V', -1. )
462           
463#if defined key_agrif && defined key_lim2
464            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
465#endif
466#if defined key_bdy
467            CALL bdy_ice_lim_dyn( 'V' )
468#endif         
469
470            DO jj = k_j1+1, k_jpj-1
471               DO ji = fs_2, fs_jpim1
472                               
473                  ! tau_io/(u_oce - u_ice)
474                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
475                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
476
477                  ! Coriolis at U-points (energy conserving formulation)
478                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
479                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
480                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
481                 
482                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
483                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
484
485                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
486                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
487                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
488                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
489                     &           ) * zmaskU(ji,jj)
490               END DO
491            END DO
492            CALL lbc_lnk( u_ice, 'U', -1. )
493           
494#if defined key_agrif && defined key_lim2
495            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
496#endif
497#if defined key_bdy
498            CALL bdy_ice_lim_dyn( 'U' )
499#endif         
500
501         ELSE ! odd iterations
502
503            DO jj = k_j1+1, k_jpj-1
504               DO ji = fs_2, fs_jpim1
505                               
506                  ! tau_io/(u_oce - u_ice)
507                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
508                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
509
510                  ! Coriolis at U-points (energy conserving formulation)
511                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
512                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
513                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
514                 
515                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
516                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
517
518                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
519                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
520                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
521                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
522                     &           ) * zmaskU(ji,jj)
523               END DO
524            END DO
525            CALL lbc_lnk( u_ice, 'U', -1. )
526           
527#if defined key_agrif && defined key_lim2
528            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
529#endif
530#if defined key_bdy
531            CALL bdy_ice_lim_dyn( 'U' )
532#endif         
533
534           DO jj = k_j1+1, k_jpj-1
535               DO ji = fs_2, fs_jpim1
536
537                  ! tau_io/(v_oce - v_ice)
538                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
539                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
540
541                  ! Coriolis at V-points (energy conserving formulation)
542                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
543                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
544                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
545
546                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
547                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
548                 
549                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
550                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
551                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
552                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
553                     &           ) * zmaskV(ji,jj)
554               END DO
555            END DO
556            CALL lbc_lnk( v_ice, 'V', -1. )
557           
558#if defined key_agrif && defined key_lim2
559            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
560#endif
561#if defined key_bdy
562            CALL bdy_ice_lim_dyn( 'V' )
563#endif         
564
565         ENDIF
566         
567         IF(ln_ctl) THEN   ! Convergence test
568            DO jj = k_j1+1, k_jpj-1
569               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
570            END DO
571            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
572            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
573         ENDIF
574         !
575         !                                                ! ==================== !
576      END DO                                              !  end loop over jter  !
577      !                                                   ! ==================== !
578      !
579      !------------------------------------------------------------------------------!
580      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
581      !------------------------------------------------------------------------------!
582      DO jj = k_j1, k_jpj-1 
583         DO ji = 1, jpim1
584
585            ! shear at F points
586            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)   &
587               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
588               &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)
589
590         END DO
591      END DO           
592      CALL lbc_lnk( zds, 'F', 1. )
593     
594      DO jj = k_j1+1, k_jpj-1 
595         DO ji = 2, jpim1 ! no vector loop
596           
597            ! tension**2 at T points
598            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)   &
599               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
600               &   ) * r1_e12t(ji,jj)
601            zdt2 = zdt * zdt
602           
603            ! shear**2 at T points (doc eq. A16)
604            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
605               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
606               &   ) * 0.25_wp * r1_e12t(ji,jj)
607           
608            ! shear at T points
609            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
610
611            ! divergence at T points
612            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
613               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
614               &            ) * r1_e12t(ji,jj)
615           
616            ! delta at T points
617            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 ) 
618            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
619            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
620
621         END DO
622      END DO
623      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
624     
625      ! --- Store the stress tensor for the next time step --- !
626      stress1_i (:,:) = zs1 (:,:)
627      stress2_i (:,:) = zs2 (:,:)
628      stress12_i(:,:) = zs12(:,:)
629
630      !
631      !------------------------------------------------------------------------------!
632      ! 5) Control prints of residual and charge ellipse
633      !------------------------------------------------------------------------------!
634      !
635      ! print the residual for convergence
636      IF(ln_ctl) THEN
637         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
638         CALL prt_ctl_info(charout)
639         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
640      ENDIF
641
642      ! print charge ellipse
643      ! This can be desactivated once the user is sure that the stress state
644      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
645      IF(ln_ctl) THEN
646         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
647         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
648         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
649         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
650            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
651            CALL prt_ctl_info(charout)
652            DO jj = k_j1+1, k_jpj-1
653               DO ji = 2, jpim1
654                  IF (zpresh(ji,jj) > 1.0) THEN
655                     zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
656                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
657                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
658                     CALL prt_ctl_info(charout)
659                  ENDIF
660               END DO
661            END DO
662            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
663            CALL prt_ctl_info(charout)
664         ENDIF
665      ENDIF
666      !
667      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
668      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
669      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
670      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
671      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
672
673   END SUBROUTINE lim_rhg
674
675#else
676   !!----------------------------------------------------------------------
677   !!   Default option          Dummy module           NO LIM sea-ice model
678   !!----------------------------------------------------------------------
679CONTAINS
680   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
681      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
682   END SUBROUTINE lim_rhg
683#endif
684
685   !!==============================================================================
686END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.