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 @ 106

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

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