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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 5429

Last change on this file since 5429 was 5429, checked in by smasson, 9 years ago

merge dev_r5302_CNRS18_HPC_scalability into the trunk, see #1523

  • Property svn:keywords set to Id
File size: 36.4 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
[3680]10   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
[1244]12   !!----------------------------------------------------------------------
[2528]13#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
[825]14   !!----------------------------------------------------------------------
[2528]15   !!   'key_lim3'               OR                     LIM-3 sea-ice model
[2717]16   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
[825]17   !!----------------------------------------------------------------------
[3625]18   !!   lim_rhg       : computes ice velocities
[825]19   !!----------------------------------------------------------------------
[3625]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
[2528]26#if defined key_lim3
[3625]27   USE ice            ! LIM-3: ice variables
28   USE dom_ice        ! LIM-3: ice domain
29   USE limitd_me      ! LIM-3:
[2528]30#else
[3625]31   USE ice_2          ! LIM-2: ice variables
32   USE dom_ice_2      ! LIM-2: ice domain
[2528]33#endif
[3625]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) 
[3680]40#if defined key_agrif && defined key_lim2
41   USE agrif_lim2_interp
42#endif
[4161]43#if defined key_bdy
44   USE bdyice_lim
45#endif
[825]46
47   IMPLICIT NONE
48   PRIVATE
49
[2715]50   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
[825]51
[868]52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
[825]54   !!----------------------------------------------------------------------
[4161]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
[5123]104      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
[834]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      !!-------------------------------------------------------------------
111      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
112      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
[834]113      !!
[2528]114      INTEGER ::   ji, jj   ! dummy loop indices
115      INTEGER ::   jter     ! local integers
[825]116      CHARACTER (len=50) ::   charout
[2528]117      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
[5123]118      REAL(wp) ::   za, zstms          ! local scalars
119      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
[825]120
[5123]121      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
122      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! 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                    ! delta on corners and on centre
126      REAL(wp) ::   zdst                          ! shear at the center of the grid point
127      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
128      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
[825]129
[2715]130      REAL(wp) ::   zresm         ! Maximal error on ice velocity
[3625]131      REAL(wp) ::   zintb, zintn  ! 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
[5123]139      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
140      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
[3294]141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
142      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
[5123]143      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
[3294]144     
[4990]145      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
[3294]146      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
147      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2
148      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
149      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
[3625]150      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
151                                                              !   ocean surface (ssh_m) if ice is not embedded
[3791]152                                                              !   ice top surface if ice is embedded   
[5123]153
154      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
155      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
[2528]156      !!-------------------------------------------------------------------
[3294]157
158      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
[5123]159      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
[4990]160      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
161      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
[3294]162
[2528]163#if  defined key_lim2 && ! defined key_lim2_vp
164# if defined key_agrif
[5123]165      USE ice_2, vt_s => hsnm
166      USE ice_2, vt_i => hicm
[2528]167# else
[5123]168      vt_s => hsnm
169      vt_i => hicm
[2528]170# endif
[5123]171      at_i(:,:) = 1. - frld(:,:)
[2528]172#endif
[3680]173#if defined key_agrif && defined key_lim2 
[5123]174      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
[3680]175#endif
[921]176      !
177      !------------------------------------------------------------------------------!
[4990]178      ! 1) Ice strength (zpresh)                                !
[921]179      !------------------------------------------------------------------------------!
180      !
[825]181      ! Put every vector to 0
[4990]182      delta_i(:,:) = 0._wp   ;
183      zpresh (:,:) = 0._wp   ; 
[2528]184      zpreshc(:,:) = 0._wp
185      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
[4990]186      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
187      shear_i(:,:) = 0._wp
[825]188
[2528]189#if defined key_lim3
[5123]190      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
[2528]191#endif
[825]192
[2528]193      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
[825]194         DO ji = 1 , jpi
[2528]195#if defined key_lim3
[5123]196            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
[2528]197#endif
[2580]198#if defined key_lim2
[5123]199            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
[2580]200#endif
[5123]201            ! zmask = 1 where there is ice or on land
202            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
[825]203         END DO
204      END DO
205
[834]206      ! Ice strength on grid cell corners (zpreshc)
207      ! needed for calculation of shear stress
[825]208      DO jj = k_j1+1, k_jpj-1
[868]209         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
[5123]210            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
211               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
212            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
213               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
214               &             ) / MAX( zstms, zepsi )
[825]215         END DO
216      END DO
217      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
[921]218      !
219      !------------------------------------------------------------------------------!
220      ! 2) Wind / ocean stress, mass terms, coriolis terms
221      !------------------------------------------------------------------------------!
222      !
[825]223      !  Wind stress, coriolis and mass terms on the sides of the squares       
224      !  zfrld1: lead fraction on U-points                                     
225      !  zfrld2: lead fraction on V-points                                     
226      !  zmass1: ice/snow mass on U-points                                   
227      !  zmass2: ice/snow mass on V-points                                   
228      !  zcorl1: Coriolis parameter on U-points                             
229      !  zcorl2: Coriolis parameter on V-points                           
230      !  (ztagnx,ztagny): wind stress on U/V points                       
231      !  v_oce1: ocean v component on u points                         
232      !  u_oce2: ocean u component on v points                         
[921]233
[3625]234      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
[5123]235         !                                           
236         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
237         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
[3625]238         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
[5123]239         !
240         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
241         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
[3625]242         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]243         !
[3625]244         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
[5123]245         !
[3625]246      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
247         zpice(:,:) = ssh_m(:,:)
248      ENDIF
249
[825]250      DO jj = k_j1+1, k_jpj-1
[868]251         DO ji = fs_2, fs_jpim1
[825]252
[5123]253            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
254            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
255            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
[4990]256
[5123]257            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
258            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
259            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
260            zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
[825]261
262            ! Leads area.
[5123]263            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
264            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
[825]265
266            ! Mass, coriolis coeff. and currents
[5123]267            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
268            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
269            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
270               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
271            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
272               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
[825]273            !
[834]274            ! Ocean has no slip boundary condition
[5123]275            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
276               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
277               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
[825]278
[5123]279            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
280               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
281               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
[825]282
[1469]283            ! Wind stress at U,V-point
284            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
285            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
[825]286
[834]287            ! Computation of the velocity field taking into account the ice internal interaction.
[825]288            ! Terms that are independent of the velocity field.
289
290            ! SB On utilise maintenant le gradient de la pente de l'ocean
291            ! include it later
[834]292
[5123]293            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
294            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
[825]295
296            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
297            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
298
299         END DO
300      END DO
301
[921]302      !
303      !------------------------------------------------------------------------------!
304      ! 3) Solution of the momentum equation, iterative procedure
305      !------------------------------------------------------------------------------!
306      !
[825]307      ! Time step for subcycling
[5123]308      dtevp  = rdt_ice / nn_nevp
309#if defined key_lim3
310      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
311#else
[2528]312      dtotel = dtevp / ( 2._wp * telast )
[5123]313#endif
314      z1_dtotel = 1._wp / ( 1._wp + dtotel )
315      z1_dtevp  = 1._wp / dtevp
[825]316      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
[5123]317      ecc2 = rn_ecc * rn_ecc
[4205]318      ecci = 1. / ecc2
[825]319
320      !-Initialise stress tensor
[2528]321      zs1 (:,:) = stress1_i (:,:) 
322      zs2 (:,:) = stress2_i (:,:)
[866]323      zs12(:,:) = stress12_i(:,:)
[825]324
[2528]325      !                                               !----------------------!
[5123]326      DO jter = 1 , nn_nevp                           !    loop over jter    !
[2528]327         !                                            !----------------------!       
[825]328         DO jj = k_j1, k_jpj-1
[2528]329            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
[825]330            zv_ice(:,jj) = v_ice(:,jj)
[921]331         END DO
[825]332
[834]333         DO jj = k_j1+1, k_jpj-1
[5123]334            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
[825]335
[921]336               
337               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
[4990]338               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
[921]339               !- zds(:,:): shear on northeast corner of grid cells
340               !
341               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,
342               !                      there are many repeated calculations.
343               !                      Speed could be improved by regrouping terms. For
344               !                      the moment, however, the stress is on clarity of coding to avoid
345               !                      bugs (Martin, for Miguel).
346               !
[4990]347               !- ALSO: arrays zdt, zds and delta could
[921]348               !  be removed in the future to minimise memory demand.
349               !
350               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
351               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
352               !              the corners is the same as in the B grid.
353               !
354               !
[5123]355               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
356                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
357                  &            ) * r1_e12t(ji,jj)
[825]358
[5123]359               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
360                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
361                  &         ) * r1_e12t(ji,jj)
[825]362
[921]363               !
[5123]364               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
365                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
366                  &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
367                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
[825]368
369
[5123]370               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
371                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
372                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
[825]373
[5123]374               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
375                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
376                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
[921]377            END DO
378         END DO
379
[5429]380         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond.
381         
[921]382         DO jj = k_j1+1, k_jpj-1
383            DO ji = fs_2, fs_jpim1
[825]384
[921]385               !- Calculate Delta at centre of grid cells
[5123]386               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
387                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
388                  &            ) * r1_e12t(ji,jj)
[825]389
[5123]390               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 
391               delta_i(ji,jj) = delta + rn_creepl
[4205]392
[921]393               !- Calculate Delta on corners
[5123]394               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
395                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
396                  &    ) * r1_e12f(ji,jj)
[825]397
[5123]398               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
399                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
400                  &    ) * r1_e12f(ji,jj)
[825]401
[5123]402               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
[825]403
[5123]404               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
405               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
406                  &          ) * z1_dtotel
407               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
408                  &          ) * z1_dtotel
409               !-Calculate stress tensor component zs12 at corners
410               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
411                  &          ) * z1_dtotel 
[4205]412
[5123]413            END DO
414         END DO
[825]415
[5429]416         CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
417 
[921]418         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
419         DO jj = k_j1+1, k_jpj-1
420            DO ji = fs_2, fs_jpim1
421               !- contribution of zs1, zs2 and zs12 to zf1
[5123]422               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
423                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
424                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
425                  &                ) * r1_e12u(ji,jj)
[921]426               ! contribution of zs1, zs2 and zs12 to zf2
[5123]427               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
428                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
429                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
430                  &               )  * r1_e12v(ji,jj)
[921]431            END DO
432         END DO
[825]433         !
434         ! Computation of ice velocity
435         !
436         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
437         !
[921]438         IF (MOD(jter,2).eq.0) THEN
[825]439
[921]440            DO jj = k_j1+1, k_jpj-1
441               DO ji = fs_2, fs_jpim1
[5123]442                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
443                  z0           = zmass1(ji,jj) * z1_dtevp
[825]444
[921]445                  ! SB modif because ocean has no slip boundary condition
[5123]446                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
447                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
448                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
449                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
450                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
451                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
452                  zcca         = z0 + za
[4990]453                  zccb         = zcorl1(ji,jj)
[5123]454                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
[921]455               END DO
456            END DO
[825]457
[921]458            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]459#if defined key_agrif && defined key_lim2
[5123]460            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]461#endif
[4333]462#if defined key_bdy
463         CALL bdy_ice_lim_dyn( 'U' )
464#endif         
[825]465
[921]466            DO jj = k_j1+1, k_jpj-1
467               DO ji = fs_2, fs_jpim1
[834]468
[5123]469                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
470                  z0           = zmass2(ji,jj) * z1_dtevp
[921]471                  ! SB modif because ocean has no slip boundary condition
[5123]472                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
473                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
474                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
475                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
476                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
477                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
478                  zcca         = z0 + za
[4990]479                  zccb         = zcorl2(ji,jj)
[5123]480                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
[921]481               END DO
482            END DO
[825]483
[921]484            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]485#if defined key_agrif && defined key_lim2
[5123]486            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]487#endif
[4333]488#if defined key_bdy
489         CALL bdy_ice_lim_dyn( 'V' )
490#endif         
[825]491
[834]492         ELSE
[921]493            DO jj = k_j1+1, k_jpj-1
494               DO ji = fs_2, fs_jpim1
[5123]495                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
496                  z0           = zmass2(ji,jj) * z1_dtevp
[921]497                  ! SB modif because ocean has no slip boundary condition
[5123]498                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
499                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
500                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
[825]501
[5123]502                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
503                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
504                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
505                  zcca         = z0 + za
[4990]506                  zccb         = zcorl2(ji,jj)
[5123]507                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
[921]508               END DO
509            END DO
[825]510
[921]511            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
[4161]512#if defined key_agrif && defined key_lim2
[5123]513            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]514#endif
[4333]515#if defined key_bdy
516         CALL bdy_ice_lim_dyn( 'V' )
517#endif         
[825]518
[921]519            DO jj = k_j1+1, k_jpj-1
520               DO ji = fs_2, fs_jpim1
[5123]521                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
522                  z0           = zmass1(ji,jj) * z1_dtevp
523                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
524                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
525                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[825]526
[5123]527                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
528                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
529                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
530                  zcca         = z0 + za
[4990]531                  zccb         = zcorl1(ji,jj)
[5123]532                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
533               END DO
534            END DO
[825]535
[921]536            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
[4161]537#if defined key_agrif && defined key_lim2
[5123]538            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]539#endif
[4333]540#if defined key_bdy
541         CALL bdy_ice_lim_dyn( 'U' )
542#endif         
[825]543
[921]544         ENDIF
[4161]545         
[921]546         IF(ln_ctl) THEN
547            !---  Convergence test.
548            DO jj = k_j1+1 , k_jpj-1
[5123]549               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
[921]550            END DO
[5123]551            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
[921]552            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
553         ENDIF
554
[4161]555         !                                                ! ==================== !
[868]556      END DO                                              !  end loop over jter  !
[825]557      !                                                   ! ==================== !
[921]558      !
559      !------------------------------------------------------------------------------!
560      ! 4) Prevent ice velocities when the ice is thin
561      !------------------------------------------------------------------------------!
[5123]562      ! If the ice volume is below zvmin then ice velocity should equal the
563      ! ocean velocity. This prevents high velocity when ice is thin
[825]564      DO jj = k_j1+1, k_jpj-1
[868]565         DO ji = fs_2, fs_jpim1
[5123]566            IF ( vt_i(ji,jj) <= zvmin ) THEN
[888]567               u_ice(ji,jj) = u_oce(ji,jj)
568               v_ice(ji,jj) = v_oce(ji,jj)
[5123]569            ENDIF
[825]570         END DO
571      END DO
[866]572
[5429]573      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
574
[4161]575#if defined key_agrif && defined key_lim2
[5123]576      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
577      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
[3680]578#endif
[4161]579#if defined key_bdy
[4333]580      CALL bdy_ice_lim_dyn( 'U' )
581      CALL bdy_ice_lim_dyn( 'V' )
[4161]582#endif         
[869]583
[868]584      DO jj = k_j1+1, k_jpj-1 
585         DO ji = fs_2, fs_jpim1
[5123]586            IF ( vt_i(ji,jj) <= zvmin ) THEN
587               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     &
588                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
589                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
[868]590
[5123]591               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
592                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
593                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
594            ENDIF
[868]595         END DO
596      END DO
597
[5429]598      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
[868]599
[866]600      ! Recompute delta, shear and div, inputs for mechanical redistribution
[825]601      DO jj = k_j1+1, k_jpj-1
[5123]602         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
[4990]603            !- divu_i(:,:), zdt(:,:): divergence and tension at centre
[825]604            !- zds(:,:): shear on northeast corner of grid cells
[5123]605            IF ( vt_i(ji,jj) <= zvmin ) THEN
[825]606
[5123]607               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
608                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
609                  &            ) * r1_e12t(ji,jj)
[825]610
[5123]611               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
612                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
613                  &         ) * r1_e12t(ji,jj)
[921]614               !
615               ! SB modif because ocean has no slip boundary condition
[5123]616               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
617                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
618                  &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
619                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
[825]620
[5123]621               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
622                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
[825]623
[5123]624               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 
625               delta_i(ji,jj) = delta + rn_creepl
[4161]626           
[5123]627            ENDIF
628         END DO
629      END DO
[921]630      !
631      !------------------------------------------------------------------------------!
632      ! 5) Store stress tensor and its invariants
633      !------------------------------------------------------------------------------!
[866]634      ! * Invariants of the stress tensor are required for limitd_me
[3791]635      !   (accelerates convergence and improves stability)
[866]636      DO jj = k_j1+1, k_jpj-1
[868]637         DO ji = fs_2, fs_jpim1
[5123]638            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
639               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
[4990]640            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
[825]641         END DO
[866]642      END DO
[4161]643
644      ! Lateral boundary condition
[5429]645      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. )
[866]646
[868]647      ! * Store the stress tensor for the next time step
648      stress1_i (:,:) = zs1 (:,:)
649      stress2_i (:,:) = zs2 (:,:)
650      stress12_i(:,:) = zs12(:,:)
651
[921]652      !
653      !------------------------------------------------------------------------------!
654      ! 6) Control prints of residual and charge ellipse
655      !------------------------------------------------------------------------------!
656      !
[834]657      ! print the residual for convergence
658      IF(ln_ctl) THEN
[868]659         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
[834]660         CALL prt_ctl_info(charout)
661         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
662      ENDIF
[825]663
[834]664      ! print charge ellipse
665      ! This can be desactivated once the user is sure that the stress state
666      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]667      IF(ln_ctl) THEN
668         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
669         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
670         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
671         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
672            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
673            CALL prt_ctl_info(charout)
674            DO jj = k_j1+1, k_jpj-1
675               DO ji = 2, jpim1
[5123]676                  IF (zpresh(ji,jj) > 1.0) THEN
[825]677                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
678                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
679                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
680                     CALL prt_ctl_info(charout)
681                  ENDIF
682               END DO
683            END DO
684            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
685            CALL prt_ctl_info(charout)
686         ENDIF
687      ENDIF
[2715]688      !
[3294]689      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
[5123]690      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
[4990]691      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
692      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
[3294]693
[825]694   END SUBROUTINE lim_rhg
695
696#else
697   !!----------------------------------------------------------------------
698   !!   Default option          Dummy module           NO LIM sea-ice model
699   !!----------------------------------------------------------------------
700CONTAINS
701   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
702      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
703   END SUBROUTINE lim_rhg
704#endif
705
706   !!==============================================================================
707END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.