source: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4269

Last change on this file since 4269 was 4269, checked in by clem, 7 years ago

update lim3 for bdy purpose

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