Ticket #78: limrhg.F90

File limrhg.F90, 38.0 KB (added by nemo_user, 13 years ago)

modified limrhg.F90

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