source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 4333

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

remove remaining bugs in LIM3, so that it can run in both regional and global config

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