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.
limrhg.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 6746

Last change on this file since 6746 was 6746, checked in by clem, 8 years ago

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