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/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 3938

Last change on this file since 3938 was 3938, checked in by flavoni, 11 years ago

dev_r3406_CNRS_LIM3: update LIM3, see ticket #1116
  • Property svn:keywords set to Id
File size: 38.9 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   !!            4.0  !  2011-01  (A Porter)  dynamical allocation
11   !!----------------------------------------------------------------------
12#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
13   !!----------------------------------------------------------------------
14   !!   'key_lim3'               OR                     LIM-3 sea-ice model
15   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
16   !!----------------------------------------------------------------------
17   !!   lim_rhg   : computes ice velocities
18   !!----------------------------------------------------------------------
19   USE phycst           ! Physical constant
20   USE par_oce          ! Ocean parameters
21   USE dom_oce          ! Ocean domain
22   USE sbc_oce          ! Surface boundary condition: ocean fields
23   USE sbc_ice          ! Surface boundary condition: ice fields
24   USE lbclnk           ! Lateral Boundary Condition / MPP link
25   USE lib_mpp          ! MPP library
26   USE wrk_nemo         ! work arrays
27   USE in_out_manager   ! I/O manager
28   USE prtctl           ! Print control
29#if defined key_lim3
30   USE ice              ! LIM-3: ice variables
31   USE dom_ice          ! LIM-3: ice domain
32   USE limitd_me        ! LIM-3:
33#else
34   USE ice_2            ! LIM2: ice variables
35   USE dom_ice_2        ! LIM2: ice domain
36#endif
37   USE lib_fortran      ! to use key_nosignedzero
38
39#if defined key_bdy
40   USE bdyice_lim
41#endif
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
47
48   REAL(wp) ::   rzero   = 0._wp   ! constant values
49   REAL(wp) ::   rone    = 1._wp   ! constant values
50   REAL(wp) ::   epsi20  = 1.e-20_wp   ! constant values
51     
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE lim_rhg( k_j1, k_jpj )
62      !!-------------------------------------------------------------------
63      !!                 ***  SUBROUTINE lim_rhg  ***
64      !!                          EVP-C-grid
65      !!
66      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
67      !!  stress and sea-surface slope. Ice-ice interaction is described by
68      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
69      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
70      !!
71      !!  The points in the C-grid look like this, dear reader
72      !!
73      !!                              (ji,jj)
74      !!                                 |
75      !!                                 |
76      !!                      (ji-1,jj)  |  (ji,jj)
77      !!                             ---------   
78      !!                            |         |
79      !!                            | (ji,jj) |------(ji,jj)
80      !!                            |         |
81      !!                             ---------   
82      !!                     (ji-1,jj-1)     (ji,jj-1)
83      !!
84      !! ** Inputs  : - wind forcing (stress), oceanic currents
85      !!                ice total volume (vt_i) per unit area
86      !!                snow total volume (vt_s) per unit area
87      !!
88      !! ** Action  : - compute u_ice, v_ice : the components of the
89      !!                sea-ice velocity vector
90      !!              - compute delta_i, shear_i, divu_i, which are inputs
91      !!                of the ice thickness distribution
92      !!
93      !! ** Steps   : 1) Compute ice snow mass, ice strength
94      !!              2) Compute wind, oceanic stresses, mass terms and
95      !!                 coriolis terms of the momentum equation
96      !!              3) Solve the momentum equation (iterative procedure)
97      !!              4) Prevent high velocities if the ice is thin
98      !!              5) Recompute invariants of the strain rate tensor
99      !!                 which are inputs of the ITD, store stress
100      !!                 for the next time step
101      !!              6) Control prints of residual (convergence)
102      !!                 and charge ellipse.
103      !!                 The user should make sure that the parameters
104      !!                 nevp, telast and creepl maintain stress state
105      !!                 on the charge ellipse for plastic flow
106      !!                 e.g. in the Canadian Archipelago
107      !!
108      !! References : Hunke and Dukowicz, JPO97
109      !!              Bouillon et al., Ocean Modelling 2009
110      !!              Vancoppenolle et al., Ocean Modelling 2008
111      !!-------------------------------------------------------------------
112      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
113      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
114      !!
115      INTEGER ::   ji, jj   ! dummy loop indices
116      INTEGER ::   jter     ! local integers
117      CHARACTER (len=50) ::   charout
118      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
119      REAL(wp) ::   za, zstms, zsang, zmask   ! local scalars
120
121      REAL(wp) ::   dtevp              ! time step for subcycling
122      REAL(wp) ::   dtotel, ecc2       ! square of yield ellipse eccenticity
123      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
124      REAL(wp) ::   zu_ice2, zv_ice1   !
125      REAL(wp) ::   zddc, zdtc, zdst   ! delta on corners and on centre
126      REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
127      REAL(wp) ::   sigma1, sigma2     ! internal ice stress
128
129      REAL(wp) ::   zresm         ! Maximal error on ice velocity
130      REAL(wp) ::   zindb         ! ice (1) or not (0)     
131      REAL(wp) ::   zdummy        ! dummy argument
132
133      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
134      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
135      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
138      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
139      REAL(wp), POINTER, DIMENSION(:,:) ::   zc1              ! ice mass
140      REAL(wp), POINTER, DIMENSION(:,:) ::   zusw             ! temporary weight for ice strength computation
141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                           
142      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points
143      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
144      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
145     
146      REAL(wp), POINTER, DIMENSION(:,:) ::   zdd   , zdt      ! Divergence and tension at centre of grid cells
147      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
148      REAL(wp), POINTER, DIMENSION(:,:) ::   deltat, deltac   ! Delta at centre and corners of grid cells
149      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2
150      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
151      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
152     
153      !!-------------------------------------------------------------------
154
155      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
156      CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
157      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
158      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         )
159
160#if  defined key_lim2 && ! defined key_lim2_vp
161# if defined key_agrif
162     USE ice_2, vt_s => hsnm
163     USE ice_2, vt_i => hicm
164# else
165     vt_s => hsnm
166     vt_i => hicm
167# endif
168     at_i(:,:) = 1. - frld(:,:)
169#endif
170      !
171      !------------------------------------------------------------------------------!
172      ! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                !
173      !------------------------------------------------------------------------------!
174      !
175      ! Put every vector to 0
176      zpresh (:,:) = 0._wp   ;   zc1   (:,:) = 0._wp
177      zpreshc(:,:) = 0._wp
178      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
179      zdd    (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
180
181#if defined key_lim3
182      CALL lim_itd_me_icestrength( ridge_scheme_swi )      ! LIM-3: Ice strength on T-points
183#endif
184
185!CDIR NOVERRCHK
186      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
187!CDIR NOVERRCHK
188         DO ji = 1 , jpi
189            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
190#if defined key_lim3
191            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp
192#endif
193#if defined key_lim2
194            zpresh(ji,jj) = tms(ji,jj) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
195#endif
196            ! tmi = 1 where there is ice or on land
197            tmi(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) ) ) ) * tms(ji,jj)
198         END DO
199      END DO
200
201      ! Ice strength on grid cell corners (zpreshc)
202      ! needed for calculation of shear stress
203!CDIR NOVERRCHK
204      DO jj = k_j1+1, k_jpj-1
205!CDIR NOVERRCHK
206         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
207            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
208               &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
209               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &
210               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1)
211            zusw(ji,jj)    = 1.0 / MAX( zstms, epsd )
212            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + &
213               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + &
214               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
215               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   &
216               &             ) * zusw(ji,jj)
217         END DO
218      END DO
219
220      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
221      !
222      !------------------------------------------------------------------------------!
223      ! 2) Wind / ocean stress, mass terms, coriolis terms
224      !------------------------------------------------------------------------------!
225      !
226      !  Wind stress, coriolis and mass terms on the sides of the squares       
227      !  zfrld1: lead fraction on U-points                                     
228      !  zfrld2: lead fraction on V-points                                     
229      !  zmass1: ice/snow mass on U-points                                   
230      !  zmass2: ice/snow mass on V-points                                   
231      !  zcorl1: Coriolis parameter on U-points                             
232      !  zcorl2: Coriolis parameter on V-points                           
233      !  (ztagnx,ztagny): wind stress on U/V points                       
234      !  u_oce1: ocean u component on u points                           
235      !  v_oce1: ocean v component on u points                         
236      !  u_oce2: ocean u component on v points                         
237      !  v_oce2: ocean v component on v points                       
238
239      DO jj = k_j1+1, k_jpj-1
240         DO ji = fs_2, fs_jpim1
241
242            zt11 = tms(ji  ,jj) * e1t(ji  ,jj)
243            zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
244            zt21 = tms(ji,jj  ) * e2t(ji,jj  )
245            zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
246
247            ! Leads area.
248            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
249            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
250
251            ! Mass, coriolis coeff. and currents
252            zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)
253            zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)
254            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
255               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
256            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   &
257               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
258            !
259            u_oce1(ji,jj)  = u_oce(ji,jj)
260            v_oce2(ji,jj)  = v_oce(ji,jj)
261
262            ! Ocean has no slip boundary condition
263            v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    &
264               &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
265               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
266
267            u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     &
268               &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
269               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
270
271            ! Wind stress at U,V-point
272            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
273            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
274
275            ! Computation of the velocity field taking into account the ice internal interaction.
276            ! Terms that are independent of the velocity field.
277
278            ! SB On utilise maintenant le gradient de la pente de l'ocean
279            ! include it later
280
281            zdsshx =  ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj)
282            zdsshy =  ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj)
283
284            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
285            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
286
287         END DO
288      END DO
289
290      !
291      !------------------------------------------------------------------------------!
292      ! 3) Solution of the momentum equation, iterative procedure
293      !------------------------------------------------------------------------------!
294      !
295      ! Time step for subcycling
296      dtevp  = rdt_ice / nevp
297      dtotel = dtevp / ( 2._wp * telast )
298
299      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
300      ecc2 = ecc * ecc
301
302      !-Initialise stress tensor
303      zs1 (:,:) = stress1_i (:,:) 
304      zs2 (:,:) = stress2_i (:,:)
305      zs12(:,:) = stress12_i(:,:)
306
307      !                                               !----------------------!
308      DO jter = 1 , nevp                              !    loop over jter    !
309         !                                            !----------------------!       
310         DO jj = k_j1, k_jpj-1
311            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
312            zv_ice(:,jj) = v_ice(:,jj)
313         END DO
314
315         DO jj = k_j1+1, k_jpj-1
316            DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
317
318               ! 
319               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
320               !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
321               !- zds(:,:): shear on northeast corner of grid cells
322               !
323               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
324               !                      there are many repeated calculations.
325               !                      Speed could be improved by regrouping terms. For
326               !                      the moment, however, the stress is on clarity of coding to avoid
327               !                      bugs (Martin, for Miguel).
328               !
329               !- ALSO: arrays zdd, zdt, zds and delta could
330               !  be removed in the future to minimise memory demand.
331               !
332               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
333               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
334               !              the corners is the same as in the B grid.
335               !
336               !
337               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
338                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
339                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
340                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
341                  &          )                                             &
342                  &         / area(ji,jj)
343
344               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
345                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
346                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
347                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
348                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
349                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
350                  &         )                                              &
351                  &        / area(ji,jj)
352
353               !
354               zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
355                  &            -u_ice(ji,jj)/e1u(ji,jj)                    &
356                  &           )*e1f(ji,jj)*e1f(ji,jj)                      &
357                  &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                &
358                  &            -v_ice(ji,jj)/e2v(ji,jj)                    &
359                  &           )*e2f(ji,jj)*e2f(ji,jj)                      &
360                  &         )                                              &
361                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
362                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
363                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
364
365
366               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
367                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
368                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
369
370               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
371                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
372                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
373
374            END DO
375         END DO
376         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
377
378!CDIR NOVERRCHK
379         DO jj = k_j1+1, k_jpj-1
380!CDIR NOVERRCHK
381            DO ji = fs_2, fs_jpim1
382
383               !- Calculate Delta at centre of grid cells
384               zdst       = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
385                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
386                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
387                  &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          &
388                  &          )                                          &
389                  &         / area(ji,jj)
390
391               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 
392               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
393
394               !-Calculate stress tensor components zs1 and zs2
395               !-at centre of grid cells (see section 3.5 of CICE user's guide).
396               zs1(ji,jj) = ( zs1(ji,jj) &
397                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
398                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
399                  * zpresh(ji,jj) ) )                          &       
400                  &        / ( 1.0 + alphaevp * dtotel )
401
402               zs2(ji,jj) = ( zs2(ji,jj)   &
403                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
404                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
405                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
406
407            END DO
408         END DO
409
410         CALL lbc_lnk( zs1(:,:), 'T', 1. )
411         CALL lbc_lnk( zs2(:,:), 'T', 1. )
412
413!CDIR NOVERRCHK
414         DO jj = k_j1+1, k_jpj-1
415!CDIR NOVERRCHK
416            DO ji = fs_2, fs_jpim1
417               !- Calculate Delta on corners
418               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
419                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
420                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
421                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
422                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
423                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
424                  &         )                                               &
425                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
426
427               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
428                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
429                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
430                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
431                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
432                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
433                  &         )                                               &
434                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
435
436               deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
437
438               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
439               zs12(ji,jj) = ( zs12(ji,jj)      &
440                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
441                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
442                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
443
444            END DO ! ji
445         END DO ! jj
446
447         CALL lbc_lnk( zs12(:,:), 'F', 1. )
448
449         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
450         DO jj = k_j1+1, k_jpj-1
451            DO ji = fs_2, fs_jpim1
452               !- contribution of zs1, zs2 and zs12 to zf1
453               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
454                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
455                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
456                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
457               ! contribution of zs1, zs2 and zs12 to zf2
458               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
459                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
460                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
461                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
462                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
463            END DO
464         END DO
465         !
466         ! Computation of ice velocity
467         !
468         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
469         !
470         IF (MOD(jter,2).eq.0) THEN 
471
472!CDIR NOVERRCHK
473            DO jj = k_j1+1, k_jpj-1
474!CDIR NOVERRCHK
475               DO ji = fs_2, fs_jpim1
476                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
477                  zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
478                  z0           = zmass1(ji,jj)/dtevp
479
480                  ! SB modif because ocean has no slip boundary condition
481                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
482                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
483                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
484                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
485                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
486                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
487                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
488                  zcca         = z0+za*cangvg
489                  zccb         = zcorl1(ji,jj)+za*zsang
490                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
491
492               END DO
493            END DO
494
495            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
496
497!CDIR NOVERRCHK
498            DO jj = k_j1+1, k_jpj-1
499!CDIR NOVERRCHK
500               DO ji = fs_2, fs_jpim1
501
502                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
503                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
504                  z0           = zmass2(ji,jj)/dtevp
505                  ! SB modif because ocean has no slip boundary condition
506                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
507                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
508                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
509                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
510                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
511                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
512                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
513                  zcca         = z0+za*cangvg
514                  zccb         = zcorl2(ji,jj)+za*zsang
515                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
516
517               END DO
518            END DO
519
520            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
521
522         ELSE 
523!CDIR NOVERRCHK
524            DO jj = k_j1+1, k_jpj-1
525!CDIR NOVERRCHK
526               DO ji = fs_2, fs_jpim1
527                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
528                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
529                  z0           = zmass2(ji,jj)/dtevp
530                  ! SB modif because ocean has no slip boundary condition
531                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      &
532                     &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
533                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
534
535                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
536                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
537                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
538                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
539                  zcca         = z0+za*cangvg
540                  zccb         = zcorl2(ji,jj)+za*zsang
541                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
542
543               END DO
544            END DO
545
546            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
547
548!CDIR NOVERRCHK
549            DO jj = k_j1+1, k_jpj-1
550!CDIR NOVERRCHK
551               DO ji = fs_2, fs_jpim1
552                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
553                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
554                  z0           = zmass1(ji,jj)/dtevp
555                  ! SB modif because ocean has no slip boundary condition
556                  ! GG Bug
557                  !                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      &
558                  !                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   &
559                  !                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
560                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
561                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
562                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
563
564                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
565                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
566                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
567                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
568                  zcca         = z0+za*cangvg
569                  zccb         = zcorl1(ji,jj)+za*zsang
570                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
571               END DO ! ji
572            END DO ! jj
573
574            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
575
576         ENDIF
577         
578!#if defined key_bdy
579!         ! clem: change u_ice and v_ice at the boundary for each iteration
580!         CALL bdy_ice_lim_dyn()
581!#endif         
582
583         IF(ln_ctl) THEN
584            !---  Convergence test.
585            DO jj = k_j1+1 , k_jpj-1
586               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
587                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
588            END DO
589            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
590            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
591         ENDIF
592
593         !                                                ! ==================== !
594      END DO                                              !  end loop over jter  !
595      !                                                   ! ==================== !
596      !
597      !------------------------------------------------------------------------------!
598      ! 4) Prevent ice velocities when the ice is thin
599      !------------------------------------------------------------------------------!
600      !clem : add hminrhg in the namelist
601      !
602      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
603      ! ocean velocity,
604      ! This prevents high velocity when ice is thin
605!CDIR NOVERRCHK
606      DO jj = k_j1+1, k_jpj-1
607!CDIR NOVERRCHK
608         DO ji = fs_2, fs_jpim1
609            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
610            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
611            IF ( zdummy .LE. hminrhg ) THEN
612               u_ice(ji,jj) = u_oce(ji,jj)
613               v_ice(ji,jj) = v_oce(ji,jj)
614            ENDIF ! zdummy
615         END DO
616      END DO
617
618      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
619      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
620
621      ! clem: change u_ice and v_ice at the boundary
622#if defined key_bdy
623      CALL bdy_ice_lim_dyn()
624#endif         
625
626      DO jj = k_j1+1, k_jpj-1 
627         DO ji = fs_2, fs_jpim1
628            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
629            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
630            IF ( zdummy .LE. hminrhg ) THEN
631               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
632                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
633                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
634
635               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
636                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
637                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
638            ENDIF ! zdummy
639         END DO
640      END DO
641
642      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
643      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
644
645      ! Recompute delta, shear and div, inputs for mechanical redistribution
646!CDIR NOVERRCHK
647      DO jj = k_j1+1, k_jpj-1
648!CDIR NOVERRCHK
649         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
650            !- zdd(:,:), zdt(:,:): divergence and tension at centre
651            !- zds(:,:): shear on northeast corner of grid cells
652            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
653            zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
654
655            IF ( zdummy .LE. hminrhg ) THEN
656
657               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
658                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
659                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
660                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
661                  &         )                                              &
662                  &         / area(ji,jj)
663
664               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
665                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
666                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
667                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
668                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
669                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
670                  &         )                                              &
671                  &        / area(ji,jj)
672               !
673               ! SB modif because ocean has no slip boundary condition
674               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              &
675                  &           - u_ice(ji,jj)   / e1u(ji,jj) )              &
676                  &           * e1f(ji,jj) * e1f(ji,jj)                    &
677                  &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             &
678                  &            - v_ice(ji,jj)  / e2v(ji,jj) )              &
679                  &           * e2f(ji,jj) * e2f(ji,jj) )                  &
680                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
681                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
682                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
683
684               zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)          &
685                  &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &
686                  &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &
687                  &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)         & 
688                  &          )                                             &
689                  &         / area(ji,jj)
690
691               deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 
692                  &                          ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 & 
693                  &                          ) + creepl
694
695            ENDIF ! zdummy
696
697         END DO !jj
698      END DO !ji
699      !
700      !------------------------------------------------------------------------------!
701      ! 5) Store stress tensor and its invariants
702      !------------------------------------------------------------------------------!
703      !
704      ! * Invariants of the stress tensor are required for limitd_me
705      ! accelerates convergence and improves stability
706      DO jj = k_j1+1, k_jpj-1
707         DO ji = fs_2, fs_jpim1
708            divu_i (ji,jj) = zdd   (ji,jj)
709            delta_i(ji,jj) = deltat(ji,jj)
710            ! begin TECLIM change
711            ! shear_i(ji,jj) = zds   (ji,jj)
712            zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
713               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
714               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
715               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)         &   
716               &          )                                             &   
717               &         / area(ji,jj) 
718            shear_i(ji,jj) = SQRT( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst )
719            ! end TECLIM change
720         END DO
721      END DO
722
723      ! Lateral boundary condition
724      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
725      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
726      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
727      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
728
729      ! * Store the stress tensor for the next time step
730      stress1_i (:,:) = zs1 (:,:)
731      stress2_i (:,:) = zs2 (:,:)
732      stress12_i(:,:) = zs12(:,:)
733
734      !
735      !------------------------------------------------------------------------------!
736      ! 6) Control prints of residual and charge ellipse
737      !------------------------------------------------------------------------------!
738      !
739      ! print the residual for convergence
740      IF(ln_ctl) THEN
741         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
742         CALL prt_ctl_info(charout)
743         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
744      ENDIF
745
746      ! print charge ellipse
747      ! This can be desactivated once the user is sure that the stress state
748      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
749      IF(ln_ctl) THEN
750         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
751         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
752         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
753         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
754            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
755            CALL prt_ctl_info(charout)
756            DO jj = k_j1+1, k_jpj-1
757               DO ji = 2, jpim1
758                  IF (zpresh(ji,jj) .GT. 1.0) THEN
759                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
760                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
761                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
762                     CALL prt_ctl_info(charout)
763                  ENDIF
764               END DO
765            END DO
766            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
767            CALL prt_ctl_info(charout)
768         ENDIF
769      ENDIF
770      !
771      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
772      CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
773      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
774      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         )
775
776   END SUBROUTINE lim_rhg
777
778#else
779   !!----------------------------------------------------------------------
780   !!   Default option          Dummy module           NO LIM sea-ice model
781   !!----------------------------------------------------------------------
782CONTAINS
783   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
784      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
785   END SUBROUTINE lim_rhg
786#endif
787
788   !!==============================================================================
789END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.