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/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limrhg.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.5 KB
Line 
1MODULE limrhg
2   !!======================================================================
3   !!                     ***  MODULE  limrhg  ***
4   !!   Ice rheology :  performs sea ice rheology
5   !!======================================================================
6#if defined key_ice_lim
7   !!----------------------------------------------------------------------
8   !!   'key_ice_lim'                                     LIM sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_rhg   : computes ice velocities
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE phycst
14   USE par_oce
15   USE ice_oce         ! ice variables
16   USE dom_ice
17   USE ice
18   USE lbclnk
19   USE lib_mpp
20   USE in_out_manager  ! I/O manager
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC lim_rhg  ! routine called by lim_dyn
27
28   !! * Module variables
29   REAL(wp)  ::           &  ! constant values
30      rzero   = 0.e0   ,  &
31      rone    = 1.e0
32   !!----------------------------------------------------------------------
33   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
34   !! $Header$
35   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE lim_rhg( k_j1, k_jpj )
41      !!-------------------------------------------------------------------
42      !!                 ***  SUBROUTINR lim_rhg  ***
43      !!
44      !! ** purpose :   determines the velocity field of sea ice by using
45      !!  atmospheric (wind stress) and oceanic (water stress and surface
46      !!  tilt) forcings. Ice-ice interaction is described by a non-linear
47      !!  viscous-plastic law including shear strength and a bulk rheology.   
48      !!
49      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity
50      !!
51      !! History :
52      !!   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
53      !!   1.0  !  94-12  (H. Goosse)
54      !!   2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
55      !!-------------------------------------------------------------------
56      ! * Arguments
57      INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation
58         &                     k_jpj      ! northern j-index for ice computation
59
60      ! * Local variables
61      INTEGER ::   ji, jj              ! dummy loop indices
62
63      INTEGER  :: &
64         iim1, ijm1, iip1 , ijp1   , & ! temporary integers
65         iter, jter                    !    "          "
66
67      REAL(wp) :: &
68         ze11  , ze12  , ze22  , ze21  ,   &  ! temporary scalars
69         zt11  , zt12  , zt21  , zt22  ,   &  !    "         "
70         zvis11, zvis21, zvis12, zvis22,   &  !    "         "
71         zgphsx, ztagnx, zusw  ,           &  !    "         "
72         zgphsy, ztagny                       !    "         "
73      REAL(wp) :: &
74         zresm, zunw, zvnw, zur, zvr, zmod, za, zac, &
75         zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  &
76         ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag
77      REAL(wp),DIMENSION(jpi,jpj) ::   &
78         zpresh, zfrld, zmass, zcorl,     &
79         zu0, zv0, zviszeta, zviseta,     &
80         zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2,  &
81         zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr
82      REAL(wp),DIMENSION(jpi,jpj,2,2) :: &
83         zs11, zs12, zs22, zs21
84      !!-------------------------------------------------------------------
85     
86      !  Store initial velocities
87      !  ------------------------
88      zu0(:,:) = u_ice(:,:)
89      zv0(:,:) = v_ice(:,:)
90
91      ! Ice mass, ice strength, and wind stress at the center            |
92      ! of the grid squares.                                             |
93      !-------------------------------------------------------------------
94
95      DO jj = k_j1 , k_jpj-1
96         DO ji = 1 , jpi
97            za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
98            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
99#if defined key_lim_cp1 && defined key_coupled
100            zb1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) )
101            zb2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) )
102#else
103            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
104            zb2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
105#endif
106         END DO
107      END DO
108
109
110      !---------------------------------------------------------------------------
111      !  Wind stress, coriolis and mass terms at the corners of the grid squares |
112      !  Gradient of ice strenght.                                               |
113      !---------------------------------------------------------------------------
114         
115      DO jj = k_j1+1, k_jpj-1
116         DO ji = 2, jpi
117            zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   &
118               &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1)
119            zusw  = 1.0 / MAX( zstms, epsd )
120
121            zt11 = tms(ji  ,jj  ) * frld(ji  ,jj  ) 
122            zt12 = tms(ji-1,jj  ) * frld(ji-1,jj  ) 
123            zt21 = tms(ji  ,jj-1) * frld(ji  ,jj-1) 
124            zt22 = tms(ji-1,jj-1) * frld(ji-1,jj-1)
125
126            ! Leads area.
127            zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   &
128               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
129
130            ! Mass and coriolis coeff.
131            zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   &
132               &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
133            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj)
134
135            ! Wind stress.
136#if defined key_lim_cp1 && defined key_coupled
137            ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
138               &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
139            ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
140               &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
141#else
142            ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
143               &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj)
144            ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
145               &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj)
146#endif
147
148            ! Gradient of ice strength
149            zgphsx =   ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji  ,jj-1)   &
150               &     + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji  ,jj  )   &
151               &     - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1)   &
152               &     - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj  )
153
154            zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji  ,jj-1)   &
155               &     - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1)   &
156               &     + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji  ,jj  )   &
157               &     + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj  )
158
159            ! Computation of the velocity field taking into account the ice-ice interaction.                                 
160            ! Terms that are independent of the velocity field.
161            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_oce(ji,jj) - zgphsx
162            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(ji,jj) - zgphsy
163         END DO
164      END DO
165
166!! inutile!!
167!!??    CALL lbc_lnk( za1ct, 'I', -1. )
168!!??    CALL lbc_lnk( za2ct, 'I', -1. )
169
170
171      ! SOLUTION OF THE MOMENTUM EQUATION.
172      !------------------------------------------
173      !                                                   ! ==================== !
174      DO iter = 1 , 2 * nbiter                            !    loop over iter    !
175         !                                                ! ==================== !       
176         zindu = MOD( iter , 2 )
177         zusdtp = ( zindu * 2.0 + ( 1.0 - zindu ) * 1.0 )  * REAL( nbiter ) / rdt_ice
178
179         ! Computation of free drift field for free slip boundary conditions.
180
181           DO jj = k_j1, k_jpj-1
182              DO ji = 1, jpim1
183                 !- Rate of strain tensor.
184                 zt11 =   akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj  ) - u_ice(ji  ,jj+1) )  &
185                    &   + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj  ) + v_ice(ji  ,jj+1) )
186                 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) )  &
187                    &   - akappa(ji,jj,2,1) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) )
188                 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) )  &
189                    &   + akappa(ji,jj,2,1) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) )
190                 zt21 =   akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj  ) - v_ice(ji  ,jj+1) )  &
191                    &   - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj  ) + u_ice(ji  ,jj+1) )
192
193                 !- Rate of strain tensor.
194                 zdgp = zt11 + zt22
195                 zdgi = zt12 + zt21
196                 ztrace2 = zdgp * zdgp 
197                 zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
198
199                 !  Creep limit depends on the size of the grid.
200                 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl)
201
202                 !-  Computation of viscosities.
203                 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
204                 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
205              END DO
206           END DO
207!!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile
208!!??       CALL lbc_lnk( zviseta , 'I', -1. )
209
210
211           !-  Determination of zc1u, zc2u, zc1v and zc2v.
212           DO jj = k_j1+1, k_jpj-1
213              DO ji = 2, jpim1
214                 ze11   =  akappa(ji-1,jj-1,1,1)
215                 ze12   = +akappa(ji-1,jj-1,2,2)
216                 ze22   =  akappa(ji-1,jj-1,2,1)
217                 ze21   = -akappa(ji-1,jj-1,1,2)
218                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
219                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
220                 zvis12 =       zviseta (ji-1,jj-1) + dm
221                 zvis21 =       zviseta (ji-1,jj-1)
222
223                 zdiag = zvis22 * ( ze11 + ze22 )
224                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
225                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
226                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
227                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
228
229                 ze11   = -akappa(ji,jj-1,1,1)
230                 ze12   = +akappa(ji,jj-1,2,2)
231                 ze22   =  akappa(ji,jj-1,2,1)
232                 ze21   = -akappa(ji,jj-1,1,2)
233                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
234                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
235                 zvis12 =       zviseta (ji,jj-1) + dm
236                 zvis21 =       zviseta (ji,jj-1)
237
238                 zdiag = zvis22 * ( ze11 + ze22 )
239                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
240                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
241                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
242                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
243
244                 ze11   =  akappa(ji-1,jj,1,1)
245                 ze12   = -akappa(ji-1,jj,2,2)
246                 ze22   =  akappa(ji-1,jj,2,1)
247                 ze21   = -akappa(ji-1,jj,1,2)
248                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
249                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
250                 zvis12 =       zviseta (ji-1,jj) + dm
251                 zvis21 =       zviseta (ji-1,jj)
252
253                 zdiag = zvis22 * ( ze11 + ze22 ) 
254                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
255                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
256                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
257                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
258
259                 ze11   = -akappa(ji,jj,1,1)
260                 ze12   = -akappa(ji,jj,2,2)
261                 ze22   =  akappa(ji,jj,2,1)
262                 ze21   = -akappa(ji,jj,1,2)
263                 zvis11 = 2.0 * zviseta (ji,jj) + dm
264                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
265                 zvis12 =       zviseta (ji,jj) + dm
266                 zvis21 =       zviseta (ji,jj)
267
268                 zdiag = zvis22 * ( ze11 + ze22 )
269                 zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
270                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
271                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
272                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
273              END DO
274           END DO
275
276           DO jj = k_j1+1, k_jpj-1
277              DO ji = 2, jpim1
278                 zc1u(ji,jj) =   &
279                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
280                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
281                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
282                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
283                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
284                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
285                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
286                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
287                 
288                 zc2u(ji,jj) =   &
289                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
290                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
291                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
292                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
293                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
294                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
295                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
296                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
297             END DO
298           END DO
299
300           DO jj = k_j1+1, k_jpj-1
301              DO ji = 2, jpim1
302                 !  zc1v , zc2v.
303                 ze11   =  akappa(ji-1,jj-1,1,2)
304                 ze12   = -akappa(ji-1,jj-1,2,1)
305                 ze22   = +akappa(ji-1,jj-1,2,2)
306                 ze21   =  akappa(ji-1,jj-1,1,1)
307                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
308                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
309                 zvis12 =       zviseta (ji-1,jj-1) + dm
310                 zvis21 =       zviseta (ji-1,jj-1)
311
312                 zdiag = zvis22 * ( ze11 + ze22 )
313                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
314                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
315                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
316                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
317 
318                 ze11   =  akappa(ji,jj-1,1,2)
319                 ze12   = -akappa(ji,jj-1,2,1)
320                 ze22   = +akappa(ji,jj-1,2,2)
321                 ze21   = -akappa(ji,jj-1,1,1)
322                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
323                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
324                 zvis12 =       zviseta (ji,jj-1) + dm
325                 zvis21 =       zviseta (ji,jj-1)
326
327                 zdiag = zvis22 * ( ze11 + ze22 )
328                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
329                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
330                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
331                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
332
333                 ze11   =  akappa(ji-1,jj,1,2)
334                 ze12   = -akappa(ji-1,jj,2,1)
335                 ze22   = -akappa(ji-1,jj,2,2)
336                 ze21   =  akappa(ji-1,jj,1,1)
337                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
338                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
339                 zvis12 =       zviseta (ji-1,jj) + dm
340                 zvis21 =       zviseta (ji-1,jj)
341
342                 zdiag = zvis22 * ( ze11 + ze22 )
343                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
344                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
345                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
346                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
347
348                 ze11   =  akappa(ji,jj,1,2)
349                 ze12   = -akappa(ji,jj,2,1)
350                 ze22   = -akappa(ji,jj,2,2)
351                 ze21   = -akappa(ji,jj,1,1)
352                 zvis11 = 2.0 * zviseta (ji,jj) + dm
353                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
354                 zvis12 =       zviseta (ji,jj) + dm
355                 zvis21 =       zviseta (ji,jj)
356
357                 zdiag = zvis22 * ( ze11 + ze22 )
358                 zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
359                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
360                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
361                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
362
363              END DO
364           END DO
365
366           DO jj = k_j1+1, k_jpj-1
367              DO ji = 2, jpim1
368                 zc1v(ji,jj) =   &
369                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
370                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
371                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
372                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
373                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
374                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
375                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
376                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
377                 zc2v(ji,jj) =   &
378                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
379                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
380                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
381                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
382                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
383                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
384                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
385                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
386              END DO
387           END DO
388
389         ! Relaxation.
390           
391iflag:   DO jter = 1 , nbitdr
392
393            !  Store previous drift field.   
394            DO jj = k_j1, k_jpj-1
395               zu_ice(:,jj) = u_ice(:,jj)
396               zv_ice(:,jj) = v_ice(:,jj)
397            END DO
398
399            DO jj = k_j1+1, k_jpj-1
400               zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
401               DO ji = 2, jpim1
402                 zur     = u_ice(ji,jj) - u_oce(ji,jj)
403                 zvr     = v_ice(ji,jj) - v_oce(ji,jj)
404                 zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) )
405                 za      = rhoco * zmod
406                 zac     = za * cangvg
407                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang
408                  zmassdt = zusdtp * zmass(ji,jj)
409                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj)
410
411                  za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
412                     &        + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) )
413
414                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
415                     &        + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) )
416
417                  zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj)
418                  zb2(ji,jj)  = zmpzas  - zc2u(ji,jj)
419                  zc1(ji,jj)  = zmpzas  + zc1v(ji,jj)
420                  zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj) 
421                  zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj)
422                  zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
423               END DO
424            END DO
425
426            ! The computation of ice interaction term is splitted into two parts
427            !-------------------------------------------------------------------------
428
429            ! Terms that do not involve already up-dated velocities.
430         
431            DO jj = k_j1+1, k_jpj-1
432               DO ji = 2, jpim1
433                  iim1 = ji
434                  ijm1 = jj - 1
435                  iip1 = ji + 1
436                  ijp1 = jj
437                  ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1)
438                  ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1)
439                  ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1)
440                  ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1)
441                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
442                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
443                  zvis12 =       zviseta (iim1,ijm1) + dm
444                  zvis21 =       zviseta (iim1,ijm1)
445                  zdiag = zvis22 * ( ze11 + ze22 )
446                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
447                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
448                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
449                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
450
451
452                  iim1 = ji - 1
453                  ijm1 = jj
454                  iip1 = ji
455                  ijp1 = jj + 1                   
456                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
457                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
458                  ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   &
459                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
460                  ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   &
461                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
462                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
463                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )
464                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
465                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
466                  zvis12 =       zviseta (iim1,ijm1) + dm
467                  zvis21 =       zviseta (iim1,ijm1)
468                  zdiag = zvis22 * ( ze11 + ze22 )
469                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
470                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
471                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
472                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
473
474                  iim1 = ji
475                  ijm1 = jj
476                  iip1 = ji + 1
477                  ijp1 = jj + 1
478                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
479                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
480                  ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   &
481                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
482                  ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   &
483                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
484                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
485                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
486                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
487                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
488                  zvis12 =       zviseta (iim1,ijm1) + dm
489                  zvis21 =       zviseta (iim1,ijm1)
490
491                  zdiag = zvis22 * ( ze11 + ze22 )
492                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
493                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
494                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
495                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
496
497               END DO
498            END DO
499
500            ! Terms involving already up-dated velocities.     
501            !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method;
502            ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method.   
503             
504            DO jj = k_j1+1, k_jpj-1
505               DO ji = 2, jpim1
506                  iim1 = ji - 1
507                  ijm1 = jj - 1
508                  iip1 = ji
509                  ijp1 = jj
510                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   &
511                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) )
512                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   &
513                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) )
514                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   &
515                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) )
516                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   &
517                     &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) )
518                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
519                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
520                  zvis12 =       zviseta (iim1,ijm1) + dm
521                  zvis21 =       zviseta (iim1,ijm1)
522
523                  zdiag = zvis22 * ( ze11 + ze22 )
524                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
525                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
526                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
527                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
528
529
530                  iim1 = ji
531                  ijm1 = jj - 1
532                  iip1 = ji + 1
533                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   &
534                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) )
535                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   &
536                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )
537                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   &
538                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )
539                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   &
540                     &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) )
541                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
542                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
543                  zvis12 =       zviseta (iim1,ijm1) + dm
544                  zvis21 =       zviseta (iim1,ijm1)
545
546                  zdiag = zvis22 * ( ze11 + ze22 )
547                  zs11(ji,jj,2,1) =  zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag
548                  zs12(ji,jj,2,1) =  zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21
549                  zs22(ji,jj,2,1) =  zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag
550                  zs21(ji,jj,2,1) =  zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12
551
552
553                  iim1 = ji - 1
554                  ijm1 = jj 
555                  ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1)
556                  ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1)
557                  ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1)
558                  ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1)
559                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
560                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
561                  zvis12 =       zviseta (iim1,ijm1) + dm
562                  zvis21 =       zviseta (iim1,ijm1)
563
564                  zdiag = zvis22 * ( ze11 + ze22 )
565                  zs11(ji,jj,1,2) =  zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 
566                  zs12(ji,jj,1,2) =  zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21
567                  zs22(ji,jj,1,2) =  zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag
568                  zs21(ji,jj,1,2) =  zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12
569
570!i             END DO
571!i          END DO
572
573!i          DO jj = k_j1+1, k_jpj-1
574!i             DO ji = 2, jpim1
575                  zd1(ji,jj) =   &
576                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)  &
577                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)  &
578                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)  &
579                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)  &
580                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)  &
581                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)  &
582                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)  &
583                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
584                  zd2(ji,jj) =   &
585                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)  &
586                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)  &
587                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)  &
588                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)  &
589                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)  &
590                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)  &
591                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)  &
592                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
593               END DO
594            END DO
595
596            DO jj = k_j1+1, k_jpj-1
597               DO ji = 2, jpim1
598                  zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        &
599                     &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj)
600
601                  zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        &
602                     &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj)
603
604                  zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
605
606                  u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
607                  v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
608               END DO
609            END DO
610
611            CALL lbc_lnk( u_ice, 'I', -1. )
612            CALL lbc_lnk( v_ice, 'I', -1. )
613
614            !---  5.2.5.4. Convergence test.
615            DO jj = k_j1+1 , k_jpj-1
616               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
617            END DO
618            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
619            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
620
621            IF ( zresm <= resl) EXIT iflag
622
623         END DO iflag
624
625         zindu1 = 1.0 - zindu
626         DO jj = k_j1 , k_jpj-1
627            zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj)
628            zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj)
629         END DO
630      !                                                   ! ==================== !
631      END DO                                              !  end loop over iter  !
632      !                                                   ! ==================== !
633
634      IF(l_ctl) THEN
635         WRITE(numout,*) ' lim_rhg  : res= ', zresm, 'iter= ', jter,' u_ice ', SUM( u_ice(2:nictl,2:njctl) ) , &
636         &                                                          ' v_ice ', SUM( v_ice(2:nictl,2:njctl) )
637      ENDIF
638
639   END SUBROUTINE lim_rhg
640
641#else
642   !!----------------------------------------------------------------------
643   !!   Default option          Dummy module           NO LIM sea-ice model
644   !!----------------------------------------------------------------------
645CONTAINS
646   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
647      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
648   END SUBROUTINE lim_rhg
649#endif
650
651   !!==============================================================================
652END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.