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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 11 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

  • Property svn:keywords set to Id
File size: 41.8 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       ! 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) * 0.5_wp
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
328      !-Initialise stress tensor
329      zs1 (:,:) = stress1_i (:,:) 
330      zs2 (:,:) = stress2_i (:,:)
331      zs12(:,:) = stress12_i(:,:)
332
333      !                                               !----------------------!
334      DO jter = 1 , nevp                              !    loop over jter    !
335         !                                            !----------------------!       
336         DO jj = k_j1, k_jpj-1
337            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
338            zv_ice(:,jj) = v_ice(:,jj)
339         END DO
340
341         DO jj = k_j1+1, k_jpj-1
342            DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
343
344               
345               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
346               !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells
347               !- zds(:,:): shear on northeast corner of grid cells
348               !
349               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
350               !                      there are many repeated calculations.
351               !                      Speed could be improved by regrouping terms. For
352               !                      the moment, however, the stress is on clarity of coding to avoid
353               !                      bugs (Martin, for Miguel).
354               !
355               !- ALSO: arrays zdd, zdt, zds and delta could
356               !  be removed in the future to minimise memory demand.
357               !
358               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
359               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
360               !              the corners is the same as in the B grid.
361               !
362               !
363               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
364                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
365                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
366                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
367                  &          )                                             &
368                  &         / area(ji,jj)
369
370               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
371                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
372                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
373                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
374                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
375                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
376                  &         )                                              &
377                  &        / area(ji,jj)
378
379               !
380               zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                &
381                  &            -u_ice(ji,jj)/e1u(ji,jj)                    &
382                  &           )*e1f(ji,jj)*e1f(ji,jj)                      &
383                  &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                &
384                  &            -v_ice(ji,jj)/e2v(ji,jj)                    &
385                  &           )*e2f(ji,jj)*e2f(ji,jj)                      &
386                  &         )                                              &
387                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
388                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
389                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
390
391
392               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
393                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
394                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
395
396               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
397                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
398                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
399
400            END DO
401         END DO
402         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
403
404!CDIR NOVERRCHK
405         DO jj = k_j1+1, k_jpj-1
406!CDIR NOVERRCHK
407            DO ji = fs_2, fs_jpim1
408
409               !- Calculate Delta at centre of grid cells
410               zzdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
411                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
412                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
413                  &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          &
414                  &          )                                          &
415                  &         / area(ji,jj)
416
417               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 
418               ! MV rewriting
419               ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl )
420               !!gm faster to replace the line above with simply:
421               !!                deltat(ji,jj) = MAX( delta, creepl )
422               !!gm end 
423               deltat(ji,jj) = delta + creepl
424               ! END MV
425               !-Calculate stress tensor components zs1 and zs2
426               !-at centre of grid cells (see section 3.5 of CICE user's guide).
427               zs1(ji,jj) = ( zs1(ji,jj) &
428                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
429                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
430                  * zpresh(ji,jj) ) )                          &       
431                  &        / ( 1.0 + alphaevp * dtotel )
432
433               zs2(ji,jj) = ( zs2(ji,jj)   &
434                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
435                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
436                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
437
438            END DO
439         END DO
440
441         CALL lbc_lnk( zs1(:,:), 'T', 1. )
442         CALL lbc_lnk( zs2(:,:), 'T', 1. )
443
444!CDIR NOVERRCHK
445         DO jj = k_j1+1, k_jpj-1
446!CDIR NOVERRCHK
447            DO ji = fs_2, fs_jpim1
448               !- Calculate Delta on corners
449               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
450                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
451                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
452                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
453                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
454                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
455                  &         )                                               &
456                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
457
458               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
459                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
460                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
461                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
462                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
463                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
464                  &         )                                               &
465                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
466
467               deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
468
469               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
470               zs12(ji,jj) = ( zs12(ji,jj)      &
471                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
472                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
473                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
474
475            END DO ! ji
476         END DO ! jj
477
478         CALL lbc_lnk( zs12(:,:), 'F', 1. )
479
480!#if defined key_bdy
481!         ! clem: change zs1, zs2, zs12 at the boundary for each iteration
482!         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 )
483!         CALL lbc_lnk( zs1 (:,:), 'T', 1. )
484!         CALL lbc_lnk( zs2 (:,:), 'T', 1. )
485!         CALL lbc_lnk( zs12(:,:), 'F', 1. )
486!#endif         
487
488         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
489         DO jj = k_j1+1, k_jpj-1
490            DO ji = fs_2, fs_jpim1
491               !- contribution of zs1, zs2 and zs12 to zf1
492               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
493                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
494                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
495                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
496               ! contribution of zs1, zs2 and zs12 to zf2
497               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
498                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
499                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
500                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
501                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
502            END DO
503         END DO
504         !
505         ! Computation of ice velocity
506         !
507         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
508         !
509         IF (MOD(jter,2).eq.0) THEN 
510
511!CDIR NOVERRCHK
512            DO jj = k_j1+1, k_jpj-1
513!CDIR NOVERRCHK
514               DO ji = fs_2, fs_jpim1
515                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
516                  zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg
517                  z0           = zmass1(ji,jj)/dtevp
518
519                  ! SB modif because ocean has no slip boundary condition
520                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
521                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
522                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
523                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
524                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
525                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
526                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
527                  zcca         = z0+za*cangvg
528                  zccb         = zcorl1(ji,jj)+za*zsang
529                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
530
531               END DO
532            END DO
533
534            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
535#if defined key_agrif && defined key_lim2
536            CALL agrif_rhg_lim2( jter, nevp, 'U' )
537#endif
538
539!CDIR NOVERRCHK
540            DO jj = k_j1+1, k_jpj-1
541!CDIR NOVERRCHK
542               DO ji = fs_2, fs_jpim1
543
544                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
545                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
546                  z0           = zmass2(ji,jj)/dtevp
547                  ! SB modif because ocean has no slip boundary condition
548                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
549                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
550                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
551                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
552                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
553                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
554                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
555                  zcca         = z0+za*cangvg
556                  zccb         = zcorl2(ji,jj)+za*zsang
557                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
558
559               END DO
560            END DO
561
562            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
563#if defined key_agrif && defined key_lim2
564            CALL agrif_rhg_lim2( jter, nevp, 'V' )
565#endif
566
567         ELSE 
568!CDIR NOVERRCHK
569            DO jj = k_j1+1, k_jpj-1
570!CDIR NOVERRCHK
571               DO ji = fs_2, fs_jpim1
572                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass2(ji,jj))))*tmv(ji,jj)
573                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
574                  z0           = zmass2(ji,jj)/dtevp
575                  ! SB modif because ocean has no slip boundary condition
576                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      &
577                     &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
578                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
579
580                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
581                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
582                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
583                     za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj))
584                  zcca         = z0+za*cangvg
585                  zccb         = zcorl2(ji,jj)+za*zsang
586                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask
587
588               END DO
589            END DO
590
591            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
592#if defined key_agrif && defined key_lim2
593            CALL agrif_rhg_lim2( jter, nevp , 'V' )
594#endif
595
596!CDIR NOVERRCHK
597            DO jj = k_j1+1, k_jpj-1
598!CDIR NOVERRCHK
599               DO ji = fs_2, fs_jpim1
600                  zmask        = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj)
601                  zsang        = SIGN(1.0,fcor(ji,jj))*sangvg
602                  z0           = zmass1(ji,jj)/dtevp
603                  ! SB modif because ocean has no slip boundary condition
604                  ! GG Bug
605                  !                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      &
606                  !                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   &
607                  !                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
608                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
609                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
610                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
611
612                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
613                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
614                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
615                     za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj))
616                  zcca         = z0+za*cangvg
617                  zccb         = zcorl1(ji,jj)+za*zsang
618                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 
619               END DO ! ji
620            END DO ! jj
621
622            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
623#if defined key_agrif && defined key_lim2
624            CALL agrif_rhg_lim2( jter, nevp, 'U' )
625#endif
626
627         ENDIF
628         
629!#if defined key_bdy
630!         ! clem: change u_ice and v_ice at the boundary for each iteration
631!         CALL bdy_ice_lim_dyn( 1 )
632!#endif         
633
634         IF(ln_ctl) THEN
635            !---  Convergence test.
636            DO jj = k_j1+1 , k_jpj-1
637               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
638                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
639            END DO
640            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
641            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
642         ENDIF
643
644         !                                                ! ==================== !
645      END DO                                              !  end loop over jter  !
646      !                                                   ! ==================== !
647      !
648      !------------------------------------------------------------------------------!
649      ! 4) Prevent ice velocities when the ice is thin
650      !------------------------------------------------------------------------------!
651      !clem : add hminrhg in the namelist
652      !
653      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the
654      ! ocean velocity,
655      ! This prevents high velocity when ice is thin
656!CDIR NOVERRCHK
657      DO jj = k_j1+1, k_jpj-1
658!CDIR NOVERRCHK
659         DO ji = fs_2, fs_jpim1
660            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
661            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
662            zdummy = vt_i(ji,jj)
663            IF ( zdummy .LE. hminrhg ) THEN
664               u_ice(ji,jj) = u_oce(ji,jj)
665               v_ice(ji,jj) = v_oce(ji,jj)
666            ENDIF ! zdummy
667         END DO
668      END DO
669
670      CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
671      CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
672#if defined key_agrif && defined key_lim2
673      CALL agrif_rhg_lim2( nevp , nevp, 'U' )
674      CALL agrif_rhg_lim2( nevp , nevp, 'V' )
675#endif
676#if defined key_bdy
677      ! clem: change u_ice and v_ice at the boundary
678      CALL bdy_ice_lim_dyn( 1 )
679#endif         
680
681      DO jj = k_j1+1, k_jpj-1 
682         DO ji = fs_2, fs_jpim1
683            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
684            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
685            zdummy = vt_i(ji,jj)
686            IF ( zdummy .LE. hminrhg ) THEN
687               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   &
688                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) &
689                  &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
690
691               u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   &
692                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) &
693                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
694            ENDIF ! zdummy
695         END DO
696      END DO
697
698      CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 
699      CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
700
701      ! Recompute delta, shear and div, inputs for mechanical redistribution
702!CDIR NOVERRCHK
703      DO jj = k_j1+1, k_jpj-1
704!CDIR NOVERRCHK
705         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
706            !- zdd(:,:), zdt(:,:): divergence and tension at centre
707            !- zds(:,:): shear on northeast corner of grid cells
708            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 
709            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 )
710            zdummy = vt_i(ji,jj)
711            IF ( zdummy .LE. hminrhg ) THEN
712
713               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
714                  &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
715                  &          +e1v(ji,jj)*v_ice(ji,jj)                      &
716                  &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
717                  &         )                                              &
718                  &         / area(ji,jj)
719
720               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
721                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
722                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
723                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
724                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
725                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
726                  &         )                                              &
727                  &        / area(ji,jj)
728               !
729               ! SB modif because ocean has no slip boundary condition
730               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              &
731                  &           - u_ice(ji,jj)   / e1u(ji,jj) )              &
732                  &           * e1f(ji,jj) * e1f(ji,jj)                    &
733                  &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             &
734                  &            - v_ice(ji,jj)  / e2v(ji,jj) )              &
735                  &           * e2f(ji,jj) * e2f(ji,jj) )                  &
736                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
737                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
738                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
739
740               zdst(ji,jj) = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
741                  &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    &
742                  &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    &
743                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
744
745!              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &
746!                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &
747!                  &                          ) + creepl
748               ! MV rewriting
749               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 
750               deltat(ji,jj) = delta + creepl
751               ! END MV
752           
753            ENDIF ! zdummy
754
755         END DO !jj
756      END DO !ji
757      !
758      !------------------------------------------------------------------------------!
759      ! 5) Store stress tensor and its invariants
760      !------------------------------------------------------------------------------!
761      !
762      ! * Invariants of the stress tensor are required for limitd_me
763      !   (accelerates convergence and improves stability)
764      DO jj = k_j1+1, k_jpj-1
765         DO ji = fs_2, fs_jpim1
766            divu_i (ji,jj) = zdd   (ji,jj)
767            delta_i(ji,jj) = deltat(ji,jj)
768            ! begin TECLIM change
769            zdst(ji,jj)= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
770               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
771               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
772               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
773            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) )
774            ! end TECLIM change
775         END DO
776      END DO
777
778      ! Lateral boundary condition
779      CALL lbc_lnk( divu_i (:,:), 'T', 1. )
780      CALL lbc_lnk( delta_i(:,:), 'T', 1. )
781      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. )
782      CALL lbc_lnk( shear_i(:,:), 'T', 1. )
783
784      ! * Store the stress tensor for the next time step
785      stress1_i (:,:) = zs1 (:,:)
786      stress2_i (:,:) = zs2 (:,:)
787      stress12_i(:,:) = zs12(:,:)
788
789      !
790      !------------------------------------------------------------------------------!
791      ! 6) Control prints of residual and charge ellipse
792      !------------------------------------------------------------------------------!
793      !
794      ! print the residual for convergence
795      IF(ln_ctl) THEN
796         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
797         CALL prt_ctl_info(charout)
798         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
799      ENDIF
800
801      ! print charge ellipse
802      ! This can be desactivated once the user is sure that the stress state
803      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
804      IF(ln_ctl) THEN
805         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
806         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
807         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
808         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
809            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
810            CALL prt_ctl_info(charout)
811            DO jj = k_j1+1, k_jpj-1
812               DO ji = 2, jpim1
813                  IF (zpresh(ji,jj) .GT. 1.0) THEN
814                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
815                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
816                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
817                     CALL prt_ctl_info(charout)
818                  ENDIF
819               END DO
820            END DO
821            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
822            CALL prt_ctl_info(charout)
823         ENDIF
824      ENDIF
825      !
826      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
827      CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
828      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  )
829      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
830
831   END SUBROUTINE lim_rhg
832
833#else
834   !!----------------------------------------------------------------------
835   !!   Default option          Dummy module           NO LIM sea-ice model
836   !!----------------------------------------------------------------------
837CONTAINS
838   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
839      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
840   END SUBROUTINE lim_rhg
841#endif
842
843   !!==============================================================================
844END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.