New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limrhg.F90 in branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

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

dev_r3406_CNRS_LIM3: update LIM3, see ticket #1116

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