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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limrhg.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 41.1 KB

Line
1	MODULE limrhg
2	!!======================================================================
3	!! * MODULE limrhg *
4	!! Ice rheology : sea ice rheology
5	!!======================================================================
6	!! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code
7	!! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3
8	!! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9	!! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas
10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
12	!! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
13	!!----------------------------------------------------------------------
14	#if defined key_lim3
15	!!----------------------------------------------------------------------
16	!! 'key_lim3' LIM-3 sea-ice model
17	!!----------------------------------------------------------------------
18	!! lim_rhg : computes ice velocities
19	!!----------------------------------------------------------------------
20	USE phycst ! Physical constant
21	USE oce , ONLY : snwice_mass, snwice_mass_b
22	USE par_oce ! Ocean parameters
23	USE dom_oce ! Ocean domain
24	USE sbc_oce ! Surface boundary condition: ocean fields
25	USE sbc_ice ! Surface boundary condition: ice fields
26	USE ice ! ice variables
27	USE limitd_me ! ice strength
28	USE lbclnk ! Lateral Boundary Condition / MPP link
29	USE lib_mpp ! MPP library
30	USE wrk_nemo ! work arrays
31	USE in_out_manager ! I/O manager
32	USE prtctl ! Print control
33	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
34	#if defined key_agrif
35	USE agrif_lim3_interp
36	#endif
37	USE bdy_oce , ONLY: ln_bdy
38	USE bdyice_lim
39
40	IMPLICIT NONE
41	PRIVATE
42
43	PUBLIC lim_rhg ! routine called by lim_dyn
44
45	!! * Substitutions
46	# include "vectopt_loop_substitute.h90"
47	!!----------------------------------------------------------------------
48	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
49	!! $Id$
50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51	!!----------------------------------------------------------------------
52	CONTAINS
53
54	SUBROUTINE lim_rhg
55	!!-------------------------------------------------------------------
56	!! * SUBROUTINE lim_rhg *
57	!! EVP-C-grid
58	!!
59	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
60	!! stress and sea-surface slope. Ice-ice interaction is described by
61	!! a non-linear elasto-viscous-plastic (EVP) law including shear
62	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
63	!!
64	!! The points in the C-grid look like this, dear reader
65	!!
66	!! (ji,jj)
67	!! \|
68	!! \|
69	!! (ji-1,jj) \| (ji,jj)
70	!! ---------
71	!! \| \|
72	!! \| (ji,jj) \|------(ji,jj)
73	!! \| \|
74	!! ---------
75	!! (ji-1,jj-1) (ji,jj-1)
76	!!
77	!! ** Inputs : - wind forcing (stress), oceanic currents
78	!! ice total volume (vt_i) per unit area
79	!! snow total volume (vt_s) per unit area
80	!!
81	!! ** Action : - compute u_ice, v_ice : the components of the
82	!! sea-ice velocity vector
83	!! - compute delta_i, shear_i, divu_i, which are inputs
84	!! of the ice thickness distribution
85	!!
86	!! ** Steps : 1) Compute ice snow mass, ice strength
87	!! 2) Compute wind, oceanic stresses, mass terms and
88	!! coriolis terms of the momentum equation
89	!! 3) Solve the momentum equation (iterative procedure)
90	!! 4) Prevent high velocities if the ice is thin
91	!! 5) Recompute invariants of the strain rate tensor
92	!! which are inputs of the ITD, store stress
93	!! for the next time step
94	!! 6) Control prints of residual (convergence)
95	!! and charge ellipse.
96	!! The user should make sure that the parameters
97	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
98	!! on the charge ellipse for plastic flow
99	!! e.g. in the Canadian Archipelago
100	!!
101	!! ** Notes : There is the possibility to use mEVP from Bouillon 2013
102	!! (by uncommenting some lines in part 3 and changing alpha and beta parameters)
103	!! but this solution appears very unstable (see Kimmritz et al 2016)
104	!!
105	!! References : Hunke and Dukowicz, JPO97
106	!! Bouillon et al., Ocean Modelling 2009
107	!! Bouillon et al., Ocean Modelling 2013
108	!!-------------------------------------------------------------------
109	INTEGER :: ji, jj ! dummy loop indices
110	INTEGER :: jter ! local integers
111	CHARACTER (len=50) :: charout
112
113	REAL(wp) :: zrhoco ! rau0 * rn_cio
114	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
115	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
116	REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013
117	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass
118	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
119	REAL(wp) :: zTauO, zTauB, zTauE, zCor, zvel ! temporary scalars
120
121	REAL(wp) :: zsig1, zsig2 ! internal ice stress
122	REAL(wp) :: zresm ! Maximal error on ice velocity
123	REAL(wp) :: zintb, zintn ! dummy argument
124
125	REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors
126	REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points
127	!
128	REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points
129	REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points
130	REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points
131	REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points
132	REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points
133	REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
134	REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses
135
136	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear
137	REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components
138	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence
139	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
140	! ocean surface (ssh_m) if ice is not embedded
141	! ice top surface if ice is embedded
142	REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays
143	REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence
144	REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice
145
146	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
147	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity
148	REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
149	!!-------------------------------------------------------------------
150
151	CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
152	CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
153	CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
154	CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
155	CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
156
157	#if defined key_agrif
158	CALL agrif_interp_lim3( 'U', 0, nn_nevp ) ! First interpolation of coarse values
159	CALL agrif_interp_lim3( 'V', 0, nn_nevp )
160	#endif
161	!
162	!------------------------------------------------------------------------------!
163	! 0) mask at F points for the ice
164	!------------------------------------------------------------------------------!
165	! ocean/land mask
166	!$OMP PARALLEL DO schedule(static) private(jj, ji)
167	DO jj = 1, jpjm1
168	DO ji = 1, jpim1 ! NO vector opt.
169	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
170	END DO
171	END DO
172	CALL lbc_lnk( zfmask, 'F', 1._wp )
173
174	! Lateral boundary conditions on velocity (modify zfmask)
175	!$OMP PARALLEL
176	!$OMP DO schedule(static) private(jj, ji)
177	DO jj = 1, jpj
178	DO ji = 1, jpi
179	zwf(ji,jj) = zfmask(ji,jj)
180	END DO
181	END DO
182	!$OMP DO schedule(static) private(jj, ji)
183	DO jj = 2, jpjm1
184	DO ji = fs_2, fs_jpim1 ! vector opt.
185	IF( zfmask(ji,jj) == 0._wp ) THEN
186	zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
187	ENDIF
188	END DO
189	END DO
190	!$OMP DO schedule(static) private(jj)
191	DO jj = 2, jpjm1
192	IF( zfmask(1,jj) == 0._wp ) THEN
193	zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
194	ENDIF
195	IF( zfmask(jpi,jj) == 0._wp ) THEN
196	zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
197	ENDIF
198	END DO
199	!$OMP DO schedule(static) private(ji)
200	DO ji = 2, jpim1
201	IF( zfmask(ji,1) == 0._wp ) THEN
202	zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
203	ENDIF
204	IF( zfmask(ji,jpj) == 0._wp ) THEN
205	zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
206	ENDIF
207	END DO
208	!$OMP END PARALLEL
209	CALL lbc_lnk( zfmask, 'F', 1._wp )
210
211	!------------------------------------------------------------------------------!
212	! 1) define some variables and initialize arrays
213	!------------------------------------------------------------------------------!
214	zrhoco = rau0 * rn_cio
215
216	! ecc2: square of yield ellipse eccenticrity
217	ecc2 = rn_ecc * rn_ecc
218	z1_ecc2 = 1._wp / ecc2
219
220	! Time step for subcycling
221	zdtevp = rdt_ice / REAL( nn_nevp )
222	z1_dtevp = 1._wp / zdtevp
223
224	! alpha parameters (Bouillon 2009)
225	zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
226	zalph2 = zalph1 * z1_ecc2
227
228	! alpha and beta parameters (Bouillon 2013)
229	!!zalph1 = 40.
230	!!zalph2 = 40.
231	!!zbeta = 3000.
232	!!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001)
233
234	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
235	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
236
237	! Initialise stress tensor
238	!$OMP PARALLEL DO schedule(static) private(jj, ji)
239	DO jj = 1, jpj
240	DO ji = 1, jpi
241	zs1 (ji,jj) = stress1_i (ji,jj)
242	zs2 (ji,jj) = stress2_i (ji,jj)
243	zs12(ji,jj) = stress12_i(ji,jj)
244	END DO
245	END DO
246
247	! Ice strength
248	CALL lim_itd_me_icestrength( nn_icestr )
249
250	! scale factors
251	!$OMP PARALLEL DO schedule(static) private(jj, ji)
252	DO jj = 2, jpjm1
253	DO ji = fs_2, fs_jpim1
254	z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )
255	z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )
256	END DO
257	END DO
258
259	!
260	!------------------------------------------------------------------------------!
261	! 2) Wind / ocean stress, mass terms, coriolis terms
262	!------------------------------------------------------------------------------!
263
264	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
265	!
266	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
267	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
268	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
269	!
270	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
271	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
272	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
273	!
274	!$OMP PARALLEL DO schedule(static) private(jj, ji)
275	DO jj = 1, jpj
276	DO ji = 1, jpi
277	zpice(ji,jj) = ssh_m(ji,jj) + ( zintn * snwice_mass(ji,jj) + zintb * snwice_mass_b(ji,jj) ) * r1_rau0
278	END DO
279	END DO
280	!
281	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
282	!$OMP PARALLEL DO schedule(static) private(jj, ji)
283	DO jj = 1, jpj
284	DO ji = 1, jpi
285	zpice(ji,jj) = ssh_m(ji,jj)
286	END DO
287	END DO
288	ENDIF
289
290	!$OMP PARALLEL DO schedule(static) private(jj,ji,zm1,zm2,zm3,zmassU,zmassV)
291	DO jj = 2, jpjm1
292	DO ji = fs_2, fs_jpim1
293
294	! ice fraction at U-V points
295	zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
296	zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
297
298	! Ice/snow mass at U-V points
299	zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
300	zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
301	zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
302	zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
303	zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
304
305	! Ocean currents at U-V points
306	v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) &
307	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
308
309	u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) &
310	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
311
312	! Coriolis at T points (m*f)
313	zmf(ji,jj) = zm1 * ff_t(ji,jj)
314
315	! m/dt
316	zmU_t(ji,jj) = zmassU * z1_dtevp
317	zmV_t(ji,jj) = zmassV * z1_dtevp
318
319	! Drag ice-atm.
320	zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
321	zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
322
323	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
324	zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
325	zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
326
327	! masks
328	zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
329	zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
330
331	! switches
332	zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
333	zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
334
335	END DO
336	END DO
337	CALL lbc_lnk( zmf, 'T', 1. )
338	!
339	!------------------------------------------------------------------------------!
340	! 3) Solution of the momentum equation, iterative procedure
341	!------------------------------------------------------------------------------!
342	!
343	! !----------------------!
344	DO jter = 1 , nn_nevp ! loop over jter !
345	! !----------------------!
346	IF(ln_ctl) THEN ! Convergence test
347	!$OMP PARALLEL DO schedule(static) private(jj, ji)
348	DO jj = 1, jpjm1
349	DO ji = 1, jpi
350	zu_ice(ji,jj) = u_ice(ji,jj) ! velocity at previous time step
351	zv_ice(ji,jj) = v_ice(ji,jj)
352	END DO
353	END DO
354	ENDIF
355
356	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
357	!$OMP PARALLEL DO schedule(static) private(jj, ji)
358	DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
359	DO ji = 1, jpim1
360
361	! shear at F points
362	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
363	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
364	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
365
366	END DO
367	END DO
368	CALL lbc_lnk( zds, 'F', 1. )
369
370	!$OMP PARALLEL DO schedule(static) private(jj,ji,zds2,zdiv,zdiv2,zdt,zdt2,zdelta)
371	DO jj = 2, jpjm1
372	DO ji = 2, jpim1 ! no vector loop
373
374	! shear**2 at T points (doc eq. A16)
375	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
376	& + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
377	& ) * 0.25_wp * r1_e1e2t(ji,jj)
378
379	! divergence at T points
380	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
381	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
382	& ) * r1_e1e2t(ji,jj)
383	zdiv2 = zdiv * zdiv
384
385	! tension at T points
386	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
387	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
388	& ) * r1_e1e2t(ji,jj)
389	zdt2 = zdt * zdt
390
391	! delta at T points
392	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
393
394	! P/delta at T points
395	zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
396
397	! stress at T points
398	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
399	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
400
401	END DO
402	END DO
403	CALL lbc_lnk( zp_delt, 'T', 1. )
404
405	!$OMP PARALLEL DO schedule(static) private(jj,ji,zp_delf)
406	DO jj = 1, jpjm1
407	DO ji = 1, jpim1
408
409	! P/delta at F points
410	zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
411
412	! stress at F points
413	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
414
415	END DO
416	END DO
417	CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
418
419
420	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
421	!$OMP PARALLEL DO schedule(static) private(jj,ji)
422	DO jj = 2, jpjm1
423	DO ji = fs_2, fs_jpim1
424
425	! U points
426	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
427	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
428	& ) * r1_e2u(ji,jj) &
429	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
430	& ) * 2._wp * r1_e1u(ji,jj) &
431	& ) * r1_e1e2u(ji,jj)
432
433	! V points
434	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
435	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
436	& ) * r1_e1v(ji,jj) &
437	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
438	& ) * 2._wp * r1_e2v(ji,jj) &
439	& ) * r1_e1e2v(ji,jj)
440
441	! u_ice at V point
442	u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
443	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
444
445	! v_ice at U point
446	v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
447	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
448
449	END DO
450	END DO
451	!
452	! --- Computation of ice velocity --- !
453	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
454	! Bouillon et al. 2009 (eq 34-35) => stable
455	IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
456
457	!$OMP PARALLEL DO schedule(static) private(jj,ji,zTauO,zvel,zTauB,zCor,zTauE,rswitch)
458	DO jj = 2, jpjm1
459	DO ji = fs_2, fs_jpim1
460
461	! tau_io/(v_oce - v_ice)
462	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
463	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
464
465	! tau_bottom/v_ice
466	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
467	zTauB = - tau_icebfr(ji,jj) / zvel
468
469	! Coriolis at V-points (energy conserving formulation)
470	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
471	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
472	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
473
474	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
475	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
476
477	! landfast switch => 0 = static friction ; 1 = sliding friction
478	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
479
480	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
481	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
482	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
483	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
484	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
485	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
486	& ) * zmaskV(ji,jj)
487	! Bouillon 2013
488	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
489	!! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
490	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
491
492	END DO
493	END DO
494	CALL lbc_lnk( v_ice, 'V', -1. )
495
496	#if defined key_agrif
497	!! CALL agrif_interp_lim3( 'V', jter, nn_nevp )
498	CALL agrif_interp_lim3( 'V' )
499	#endif
500	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
501
502	!$OMP PARALLEL DO schedule(static) private(jj,ji,zTauO,zvel,zTauB,zCor,zTauE,rswitch)
503	DO jj = 2, jpjm1
504	DO ji = fs_2, fs_jpim1
505
506	! tau_io/(u_oce - u_ice)
507	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
508	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
509
510	! tau_bottom/u_ice
511	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
512	zTauB = - tau_icebfr(ji,jj) / zvel
513
514	! Coriolis at U-points (energy conserving formulation)
515	zCor = 0.25_wp * r1_e1u(ji,jj) * &
516	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
517	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
518
519	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
520	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
521
522	! landfast switch => 0 = static friction ; 1 = sliding friction
523	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
524
525	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
526	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
527	& + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
528	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
529	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
530	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
531	& ) * zmaskU(ji,jj)
532	! Bouillon 2013
533	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
534	!! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
535	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
536	END DO
537	END DO
538	CALL lbc_lnk( u_ice, 'U', -1. )
539
540	#if defined key_agrif
541	!! CALL agrif_interp_lim3( 'U', jter, nn_nevp )
542	CALL agrif_interp_lim3( 'U' )
543	#endif
544	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
545
546	ELSE ! odd iterations
547
548	!$OMP PARALLEL DO schedule(static) private(jj,ji,zTauO,zvel,zTauB,zCor,zTauE,rswitch)
549	DO jj = 2, jpjm1
550	DO ji = fs_2, fs_jpim1
551
552	! tau_io/(u_oce - u_ice)
553	zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
554	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
555
556	! tau_bottom/u_ice
557	zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
558	zTauB = - tau_icebfr(ji,jj) / zvel
559
560	! Coriolis at U-points (energy conserving formulation)
561	zCor = 0.25_wp * r1_e1u(ji,jj) * &
562	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
563	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
564
565	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
566	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
567
568	! landfast switch => 0 = static friction ; 1 = sliding friction
569	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
570
571	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
572	u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity
573	& + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
574	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
575	& + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
576	& ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
577	& ) * zmaskU(ji,jj)
578	! Bouillon 2013
579	!!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) &
580	!! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) &
581	!! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
582	END DO
583	END DO
584	CALL lbc_lnk( u_ice, 'U', -1. )
585
586	#if defined key_agrif
587	!! CALL agrif_interp_lim3( 'U', jter, nn_nevp )
588	CALL agrif_interp_lim3( 'U' )
589	#endif
590	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
591
592	!$OMP PARALLEL DO schedule(static) private(jj,ji,zTauO,zvel,zTauB,zCor,zTauE,rswitch)
593	DO jj = 2, jpjm1
594	DO ji = fs_2, fs_jpim1
595
596	! tau_io/(v_oce - v_ice)
597	zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
598	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
599
600	! tau_bottom/v_ice
601	zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
602	ztauB = - tau_icebfr(ji,jj) / zvel
603
604	! Coriolis at V-points (energy conserving formulation)
605	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
606	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
607	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
608
609	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
610	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
611
612	! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
613	rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
614
615	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
616	v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity
617	& + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
618	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
619	& + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0
620	& ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
621	& ) * zmaskV(ji,jj)
622	! Bouillon 2013
623	!!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) &
624	!! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) &
625	!! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
626	END DO
627	END DO
628	CALL lbc_lnk( v_ice, 'V', -1. )
629
630	#if defined key_agrif
631	!! CALL agrif_interp_lim3( 'V', jter, nn_nevp )
632	CALL agrif_interp_lim3( 'V' )
633	#endif
634	IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
635
636	ENDIF
637
638	IF(ln_ctl) THEN ! Convergence test
639	!$OMP PARALLEL DO schedule(static) private(jj, ji)
640	DO jj = 2 , jpjm1
641	DO ji = 1, jpi
642	zresr(ji,jj) = MAX( ABS( u_ice(ji,jj) - zu_ice(ji,jj) ), ABS( v_ice(ji,jj) - zv_ice(ji,jj) ) )
643	END DO
644	END DO
645	zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
646	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
647	ENDIF
648	!
649	! ! ==================== !
650	END DO ! end loop over jter !
651	! ! ==================== !
652	!
653	!------------------------------------------------------------------------------!
654	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
655	!------------------------------------------------------------------------------!
656	!$OMP PARALLEL DO schedule(static) private(jj, ji)
657	DO jj = 1, jpjm1
658	DO ji = 1, jpim1
659
660	! shear at F points
661	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
662	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
663	& ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
664
665	END DO
666	END DO
667	CALL lbc_lnk( zds, 'F', 1. )
668
669	!$OMP PARALLEL DO schedule(static) private(jj,ji,zdt,zdt2,zds2,zdelta,rswitch)
670	DO jj = 2, jpjm1
671	DO ji = 2, jpim1 ! no vector loop
672
673	! tension**2 at T points
674	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
675	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
676	& ) * r1_e1e2t(ji,jj)
677	zdt2 = zdt * zdt
678
679	! shear**2 at T points (doc eq. A16)
680	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
681	& + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
682	& ) * 0.25_wp * r1_e1e2t(ji,jj)
683
684	! shear at T points
685	shear_i(ji,jj) = SQRT( zdt2 + zds2 )
686
687	! divergence at T points
688	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
689	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
690	& ) * r1_e1e2t(ji,jj)
691
692	! delta at T points
693	zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )
694	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
695	delta_i(ji,jj) = zdelta + rn_creepl * rswitch
696
697	END DO
698	END DO
699	CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
700
701	! --- Store the stress tensor for the next time step --- !
702	!$OMP PARALLEL DO schedule(static) private(jj, ji)
703	DO jj = 1, jpj
704	DO ji = 1, jpi
705	stress1_i (ji,jj) = zs1 (ji,jj)
706	stress2_i (ji,jj) = zs2 (ji,jj)
707	stress12_i(ji,jj) = zs12(ji,jj)
708	END DO
709	END DO
710	!
711
712	!------------------------------------------------------------------------------!
713	! 5) Control prints of residual and charge ellipse
714	!------------------------------------------------------------------------------!
715	!
716	! print the residual for convergence
717	IF(ln_ctl) THEN
718	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
719	CALL prt_ctl_info(charout)
720	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
721	ENDIF
722
723	! print charge ellipse
724	! This can be desactivated once the user is sure that the stress state
725	! lie on the charge ellipse. See Bouillon et al. 08 for more details
726	IF(ln_ctl) THEN
727	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
728	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
729	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
730	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
731	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
732	CALL prt_ctl_info(charout)
733	DO jj = 2, jpjm1
734	DO ji = 2, jpim1
735	IF (strength(ji,jj) > 1.0) THEN
736	zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
737	zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)*2 + 4zs12(ji,jj)*2 ) ) / ( 2strength(ji,jj) )
738	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
739	CALL prt_ctl_info(charout)
740	ENDIF
741	END DO
742	END DO
743	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
744	CALL prt_ctl_info(charout)
745	ENDIF
746	ENDIF
747	!
748
749	CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
750	CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
751	CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
752	CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
753	CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
754
755	END SUBROUTINE lim_rhg
756
757	#else
758	!!----------------------------------------------------------------------
759	!! Default option Dummy module NO LIM sea-ice model
760	!!----------------------------------------------------------------------
761	CONTAINS
762	SUBROUTINE lim_rhg ! Dummy routine
763	WRITE(,) 'lim_rhg: You should not have seen this print! error?'
764	END SUBROUTINE lim_rhg
765	#endif
766
767	!!==============================================================================
768	END MODULE limrhg

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

Download in other formats: