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

Context Navigation

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

Last change on this file since 12 was 12, checked in by opalod, 21 years ago
CT : UPDATE001 : First major NEMO update
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 35.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: