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.

limitd_me.F90 in branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5064

Last change on this file since 5064 was 5064, checked in by clem, 9 years ago
LIM3 now includes specific physics to run with only 1 sea ice category (i.e. LIM2 type)
Property svn:keywords set to `Id`
File size: 63.8 KB

Line
1	MODULE limitd_me
2	!!======================================================================
3	!! * MODULE limitd_me *
4	!! LIM-3 : Mechanical impact on ice thickness distribution
5	!!======================================================================
6	!! History : LIM ! 2006-02 (M. Vancoppenolle) Original code
7	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
9	!!----------------------------------------------------------------------
10	#if defined key_lim3
11	!!----------------------------------------------------------------------
12	!! 'key_lim3' LIM-3 sea-ice model
13	!!----------------------------------------------------------------------
14	USE par_oce ! ocean parameters
15	USE dom_oce ! ocean domain
16	USE phycst ! physical constants (ocean directory)
17	USE sbc_oce ! surface boundary condition: ocean fields
18	USE thd_ice ! LIM thermodynamics
19	USE ice ! LIM variables
20	USE dom_ice ! LIM domain
21	USE limthd_lac ! LIM
22	USE limvar ! LIM
23	USE in_out_manager ! I/O manager
24	USE lbclnk ! lateral boundary condition - MPP exchanges
25	USE lib_mpp ! MPP library
26	USE wrk_nemo ! work arrays
27	USE prtctl ! Print control
28
29	USE iom ! I/O manager
30	USE lib_fortran ! glob_sum
31	USE limdiahsb
32	USE timing ! Timing
33	USE limcons ! conservation tests
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC lim_itd_me ! called by ice_stp
39	PUBLIC lim_itd_me_icestrength
40	PUBLIC lim_itd_me_init
41	PUBLIC lim_itd_me_alloc ! called by sbc_lim_init
42
43	!-----------------------------------------------------------------------
44	! Variables shared among ridging subroutines
45	!-----------------------------------------------------------------------
46	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
47	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
48	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/
49	! ! closing associated w/ category n
50	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
51	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice
54	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
55	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
56
57	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
58	REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer
59	REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
60
61	REAL(wp) :: Cp !
62	!
63	!-----------------------------------------------------------------------
64	! Ridging diagnostic arrays for history files
65	!-----------------------------------------------------------------------
66	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s)
67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s)
68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s)
69	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s)
70	!
71	!!----------------------------------------------------------------------
72	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
73	!! $Id$
74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
75	!!----------------------------------------------------------------------
76	CONTAINS
77
78	INTEGER FUNCTION lim_itd_me_alloc()
79	!!---------------------------------------------------------------------!
80	!! * ROUTINE lim_itd_me_alloc *
81	!!---------------------------------------------------------------------!
82	ALLOCATE( &
83	!* Variables shared among ridging subroutines
84	& asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , &
85	& aksum(jpi,jpj) , &
86	!
87	& hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
88	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , &
89	!
90	!* Ridging diagnostic arrays for history files
91	& dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , &
92	& dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc )
93	!
94	IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
95	!
96	END FUNCTION lim_itd_me_alloc
97
98
99	SUBROUTINE lim_itd_me
100	!!---------------------------------------------------------------------!
101	!! * ROUTINE lim_itd_me *
102	!!
103	!! ** Purpose : computes the mechanical redistribution of ice thickness
104	!!
105	!! ** Method : Steps :
106	!! 1) Thickness categories boundaries, ice / o.w. concentrations
107	!! Ridge preparation
108	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
109	!! 3) Ridging iteration
110	!! 4) Ridging diagnostics
111	!! 5) Heat, salt and freshwater fluxes
112	!! 6) Compute increments of tate variables and come back to old values
113	!!
114	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
115	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
116	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
117	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
118	!! Bitz et al., JGR, 2001
119	!! Amundrud and Melling, JGR 2005
120	!! Babko et al., JGR 2002
121	!!
122	!! This routine is based on CICE code and authors William H. Lipscomb,
123	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
124	!!--------------------------------------------------------------------!
125	INTEGER :: ji, jj, jk, jl ! dummy loop index
126	INTEGER :: niter, nitermax = 20 ! local integer
127	LOGICAL :: asum_error ! flag for asum .ne. 1
128	INTEGER :: iterate_ridging ! if true, repeat the ridging
129	REAL(wp) :: w1, tmpfac ! local scalar
130	CHARACTER (len = 15) :: fieldid
131	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
132	! (ridging ice area - area of new ridges) / dt
133	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
134	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
135	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
136	REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
137	REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
138	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
139	!
140	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
141	!!-----------------------------------------------------------------------------
142	IF( nn_timing == 1 ) CALL timing_start('limitd_me')
143
144	CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
145
146	IF(ln_ctl) THEN
147	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
148	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
149	ENDIF
150
151	IF( ln_limdyn ) THEN ! Start ridging and rafting !
152
153	! conservation test
154	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
155
156	!-----------------------------------------------------------------------------!
157	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
158	!-----------------------------------------------------------------------------!
159	Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0 ! proport const for PE
160	!
161	CALL lim_itd_me_ridgeprep ! prepare ridging
162	!
163	IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check
164
165	DO jj = 1, jpj ! Initialize arrays.
166	DO ji = 1, jpi
167	msnow_mlt(ji,jj) = 0._wp
168	esnow_mlt(ji,jj) = 0._wp
169	dardg1dt (ji,jj) = 0._wp
170	dardg2dt (ji,jj) = 0._wp
171	dvirdgdt (ji,jj) = 0._wp
172	opening (ji,jj) = 0._wp
173
174	!-----------------------------------------------------------------------------!
175	! 2) Dynamical inputs (closing rate, divu_adv, opning)
176	!-----------------------------------------------------------------------------!
177	!
178	! 2.1 closing_net
179	!-----------------
180	! Compute the net rate of closing due to convergence
181	! and shear, based on Flato and Hibler (1995).
182	!
183	! The energy dissipation rate is equal to the net closing rate
184	! times the ice strength.
185	!
186	! NOTE: The NET closing rate is equal to the rate that open water
187	! area is removed, plus the rate at which ice area is removed by
188	! ridging, minus the rate at which area is added in new ridges.
189	! The GROSS closing rate is equal to the first two terms (open
190	! water closing and thin ice ridging) without the third term
191	! (thick, newly ridged ice).
192
193	closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
194
195	! 2.2 divu_adv
196	!--------------
197	! Compute divu_adv, the divergence rate given by the transport/
198	! advection scheme, which may not be equal to divu as computed
199	! from the velocity field.
200	!
201	! If divu_adv < 0, make sure the closing rate is large enough
202	! to give asum = 1.0 after ridging.
203
204	divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
205
206	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
207
208	! 2.3 opning
209	!------------
210	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
211	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
212	END DO
213	END DO
214
215	!-----------------------------------------------------------------------------!
216	! 3) Ridging iteration
217	!-----------------------------------------------------------------------------!
218	niter = 1 ! iteration counter
219	iterate_ridging = 1
220
221	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
222
223	DO jj = 1, jpj
224	DO ji = 1, jpi
225
226	! 3.2 closing_gross
227	!-----------------------------------------------------------------------------!
228	! Based on the ITD of ridging and ridged ice, convert the net
229	! closing rate to a gross closing rate.
230	! NOTE: 0 < aksum <= 1
231	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
232
233	! correction to closing rate and opening if closing rate is excessive
234	!---------------------------------------------------------------------
235	! Reduce the closing rate if more than 100% of the open water
236	! would be removed. Reduce the opening rate proportionately.
237	IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
238	w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
239	IF ( w1 .GT. ato_i(ji,jj)) THEN
240	tmpfac = ato_i(ji,jj) / w1
241	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
242	opning(ji,jj) = opning(ji,jj) * tmpfac
243	ENDIF !w1
244	ENDIF !at0i and athorn
245
246	END DO ! ji
247	END DO ! jj
248
249	! correction to closing rate / opening if excessive ice removal
250	!---------------------------------------------------------------
251	! Reduce the closing rate if more than 100% of any ice category
252	! would be removed. Reduce the opening rate proportionately.
253
254	DO jl = 1, jpl
255	DO jj = 1, jpj
256	DO ji = 1, jpi
257	IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN
258	w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
259	IF ( w1 > a_i(ji,jj,jl) ) THEN
260	tmpfac = a_i(ji,jj,jl) / w1
261	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
262	opning (ji,jj) = opning (ji,jj) * tmpfac
263	ENDIF
264	ENDIF
265	END DO !ji
266	END DO ! jj
267	END DO !jl
268
269	! 3.3 Redistribute area, volume, and energy.
270	!-----------------------------------------------------------------------------!
271
272	CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
273
274	! 3.4 Compute total area of ice plus open water after ridging.
275	!-----------------------------------------------------------------------------!
276
277	CALL lim_itd_me_asumr
278
279	! 3.5 Do we keep on iterating ???
280	!-----------------------------------------------------------------------------!
281	! Check whether asum = 1. If not (because the closing and opening
282	! rates were reduced above), ridge again with new rates.
283
284	iterate_ridging = 0
285
286	DO jj = 1, jpj
287	DO ji = 1, jpi
288	IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN
289	closing_net(ji,jj) = 0._wp
290	opning (ji,jj) = 0._wp
291	ELSE
292	iterate_ridging = 1
293	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
294	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
295	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
296	ENDIF
297	END DO
298	END DO
299
300	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
301
302	! Repeat if necessary.
303	! NOTE: If strength smoothing is turned on, the ridging must be
304	! iterated globally because of the boundary update in the
305	! smoothing.
306
307	niter = niter + 1
308
309	IF( iterate_ridging == 1 ) THEN
310	IF( niter > nitermax ) THEN
311	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
312	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
313	ENDIF
314	CALL lim_itd_me_ridgeprep
315	ENDIF
316
317	END DO !! on the do while over iter
318
319	!-----------------------------------------------------------------------------!
320	! 4) Ridging diagnostics
321	!-----------------------------------------------------------------------------!
322	! Convert ridging rate diagnostics to correct units.
323	! Update fresh water and heat fluxes due to snow melt.
324
325	asum_error = .false.
326
327	DO jj = 1, jpj
328	DO ji = 1, jpi
329
330	IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.
331
332	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
333	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
334	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
335	opening (ji,jj) = opening (ji,jj) * r1_rdtice
336
337	!-----------------------------------------------------------------------------!
338	! 5) Heat, salt and freshwater fluxes
339	!-----------------------------------------------------------------------------!
340	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
341	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
342
343	END DO
344	END DO
345
346	! Check if there is a ridging error
347	DO jj = 1, jpj
348	DO ji = 1, jpi
349	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
350	WRITE(numout,*) ' '
351	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
352	WRITE(numout,*) ' limitd_me '
353	WRITE(numout,*) ' POINT : ', ji, jj
354	WRITE(numout,*) ' jpl, a_i, athorn '
355	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
356	DO jl = 1, jpl
357	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
358	END DO
359	ENDIF ! asum
360
361	END DO !ji
362	END DO !jj
363
364	! Conservation check
365	IF ( con_i ) THEN
366	CALL lim_column_sum (jpl, v_i, vt_i_final)
367	fieldid = ' v_i : limitd_me '
368	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
369	ENDIF
370
371	!-----------------------------------------------------------------------------!
372	! 6) Updating state variables and trend terms (done in limupdate)
373	!-----------------------------------------------------------------------------!
374	CALL lim_var_glo2eqv
375	CALL lim_var_zapsmall
376	CALL lim_var_agg( 1 )
377
378
379	IF(ln_ctl) THEN ! Control print
380	CALL prt_ctl_info(' ')
381	CALL prt_ctl_info(' - Cell values : ')
382	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
383	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :')
384	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
385	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
386	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
387	DO jl = 1, jpl
388	CALL prt_ctl_info(' ')
389	CALL prt_ctl_info(' - Category : ', ivar1=jl)
390	CALL prt_ctl_info(' ~~~~~~~~~~')
391	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
392	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
393	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
394	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
395	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
396	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
397	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
398	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
399	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
400	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
401	DO jk = 1, nlay_i
402	CALL prt_ctl_info(' ')
403	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
404	CALL prt_ctl_info(' ~~~~~~~')
405	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
406	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
407	END DO
408	END DO
409	ENDIF
410
411	! conservation test
412	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
413
414	ENDIF ! ln_limdyn=.true.
415	!
416	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
417	!
418	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
419	END SUBROUTINE lim_itd_me
420
421
422	SUBROUTINE lim_itd_me_icestrength( kstrngth )
423	!!----------------------------------------------------------------------
424	!! * ROUTINE lim_itd_me_icestrength *
425	!!
426	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
427	!!
428	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
429	!! dissipated per unit area removed from the ice pack under compression,
430	!! and assumed proportional to the change in potential energy caused
431	!! by ridging. Note that only Hibler's formulation is stable and that
432	!! ice strength has to be smoothed
433	!!
434	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
435	!!----------------------------------------------------------------------
436	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
437
438	INTEGER :: ji,jj, jl ! dummy loop indices
439	INTEGER :: ksmooth ! smoothing the resistance to deformation
440	INTEGER :: numts_rm ! number of time steps for the P smoothing
441	REAL(wp) :: hi, zw1, zp, zdummy, zzc, z1_3 ! local scalars
442	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
443	!!----------------------------------------------------------------------
444
445	CALL wrk_alloc( jpi, jpj, zworka )
446
447	!------------------------------------------------------------------------------!
448	! 1) Initialize
449	!------------------------------------------------------------------------------!
450	strength(:,:) = 0._wp
451
452	!------------------------------------------------------------------------------!
453	! 2) Compute thickness distribution of ridging and ridged ice
454	!------------------------------------------------------------------------------!
455	CALL lim_itd_me_ridgeprep
456
457	!------------------------------------------------------------------------------!
458	! 3) Rothrock(1975)'s method
459	!------------------------------------------------------------------------------!
460	IF( kstrngth == 1 ) THEN
461	z1_3 = 1._wp / 3._wp
462	DO jl = 1, jpl
463	DO jj= 1, jpj
464	DO ji = 1, jpi
465	!
466	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
467	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
468	!----------------------------
469	! PE loss from deforming ice
470	!----------------------------
471	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi
472
473	!--------------------------
474	! PE gain from rafting ice
475	!--------------------------
476	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi
477
478	!----------------------------
479	! PE gain from ridging ice
480	!----------------------------
481	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl) &
482	* z1_3 * (hrmax(ji,jj,jl)3 - hrmin(ji,jj,jl)3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )
483	!!gm Optimization: (a3-b3)/(a-b) = aa+ab+bb ==> less costly operations even if a*3 is replaced by aa*a...
484	ENDIF ! aicen > epsi10
485	!
486	END DO ! ji
487	END DO !jj
488	END DO !jl
489
490	zzc = Cf * Cp ! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and Cf accounts for frictional dissipation
491	strength(:,:) = zzc * strength(:,:) / aksum(:,:)
492
493	ksmooth = 1
494
495	!------------------------------------------------------------------------------!
496	! 4) Hibler (1979)' method
497	!------------------------------------------------------------------------------!
498	ELSE ! kstrngth ne 1: Hibler (1979) form
499	!
500	strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) ) )
501	!
502	ksmooth = 1
503	!
504	ENDIF ! kstrngth
505
506	!
507	!------------------------------------------------------------------------------!
508	! 5) Impact of brine volume
509	!------------------------------------------------------------------------------!
510	! CAN BE REMOVED
511	!
512	IF ( brinstren_swi == 1 ) THEN
513
514	DO jj = 1, jpj
515	DO ji = 1, jpi
516	IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
517	zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
518	ELSE
519	zdummy = 0.0
520	ENDIF
521	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
522	END DO ! j
523	END DO ! i
524
525	ENDIF
526
527	!
528	!------------------------------------------------------------------------------!
529	! 6) Smoothing ice strength
530	!------------------------------------------------------------------------------!
531	!
532	!-------------------
533	! Spatial smoothing
534	!-------------------
535	IF ( ksmooth == 1 ) THEN
536
537	CALL lbc_lnk( strength, 'T', 1. )
538
539	DO jj = 2, jpj - 1
540	DO ji = 2, jpi - 1
541	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
542	! present
543	zworka(ji,jj) = 4.0 * strength(ji,jj) &
544	& + strength(ji-1,jj) * tms(ji-1,jj) &
545	& + strength(ji+1,jj) * tms(ji+1,jj) &
546	& + strength(ji,jj-1) * tms(ji,jj-1) &
547	& + strength(ji,jj+1) * tms(ji,jj+1)
548
549	zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)
550	zworka(ji,jj) = zworka(ji,jj) / zw1
551	ELSE
552	zworka(ji,jj) = 0._wp
553	ENDIF
554	END DO
555	END DO
556
557	DO jj = 2, jpj - 1
558	DO ji = 2, jpi - 1
559	strength(ji,jj) = zworka(ji,jj)
560	END DO
561	END DO
562	CALL lbc_lnk( strength, 'T', 1. )
563
564	ENDIF ! ksmooth
565
566	!--------------------
567	! Temporal smoothing
568	!--------------------
569	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
570	strp1(:,:) = 0.0
571	strp2(:,:) = 0.0
572	ENDIF
573
574	IF ( ksmooth == 2 ) THEN
575
576
577	CALL lbc_lnk( strength, 'T', 1. )
578
579	DO jj = 1, jpj - 1
580	DO ji = 1, jpi - 1
581	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is present
582	numts_rm = 1 ! number of time steps for the running mean
583	IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
584	IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
585	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
586	strp2(ji,jj) = strp1(ji,jj)
587	strp1(ji,jj) = strength(ji,jj)
588	strength(ji,jj) = zp
589
590	ENDIF
591	END DO
592	END DO
593
594	ENDIF ! ksmooth
595
596	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
597
598	CALL wrk_dealloc( jpi, jpj, zworka )
599	!
600	END SUBROUTINE lim_itd_me_icestrength
601
602
603	SUBROUTINE lim_itd_me_ridgeprep
604	!!---------------------------------------------------------------------!
605	!! * ROUTINE lim_itd_me_ridgeprep *
606	!!
607	!! ** Purpose : preparation for ridging and strength calculations
608	!!
609	!! ** Method : Compute the thickness distribution of the ice and open water
610	!! participating in ridging and of the resulting ridges.
611	!!---------------------------------------------------------------------!
612	INTEGER :: ji,jj, jl ! dummy loop indices
613	REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar
614	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
615	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
616	!------------------------------------------------------------------------------!
617
618	CALL wrk_alloc( jpi,jpj, zworka )
619	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
620
621	Gstari = 1.0/Gstar
622	astari = 1.0/astar
623	aksum(:,:) = 0.0
624	athorn(:,:,:) = 0.0
625	aridge(:,:,:) = 0.0
626	araft (:,:,:) = 0.0
627	hrmin(:,:,:) = 0.0
628	hrmax(:,:,:) = 0.0
629	hraft(:,:,:) = 0.0
630	krdg (:,:,:) = 1.0
631
632	! ! Zero out categories with very small areas
633	CALL lim_var_zapsmall
634
635	!------------------------------------------------------------------------------!
636	! 1) Participation function
637	!------------------------------------------------------------------------------!
638
639	! Compute total area of ice plus open water.
640	CALL lim_itd_me_asumr
641	! This is in general not equal to one
642	! because of divergence during transport
643
644	! Compute cumulative thickness distribution function
645	! Compute the cumulative thickness distribution function Gsum,
646	! where Gsum(n) is the fractional area in categories 0 to n.
647	! initial value (in h = 0) equals open water area
648
649	Gsum(:,:,-1) = 0._wp
650
651	DO jj = 1, jpj
652	DO ji = 1, jpi
653	IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj)
654	ELSE ; Gsum(ji,jj,0) = 0._wp
655	ENDIF
656	END DO
657	END DO
658
659	! for each value of h, you have to add ice concentration then
660	DO jl = 1, jpl
661	DO jj = 1, jpj
662	DO ji = 1, jpi
663	IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
664	ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
665	ENDIF
666	END DO
667	END DO
668	END DO
669
670	! Normalize the cumulative distribution to 1
671	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
672	DO jl = 0, jpl
673	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
674	END DO
675
676	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
677	!--------------------------------------------------------------------------------------------------
678	! Compute the participation function athorn; this is analogous to
679	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
680	! area lost from category n due to ridging/closing
681	! athorn(n) = total area lost due to ridging/closing
682	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
683	!
684	! The expressions for athorn are found by integrating b(h)g(h) between
685	! the category boundaries.
686	!-----------------------------------------------------------------
687
688	IF( partfun_swi == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
689	DO jl = 0, jpl
690	DO jj = 1, jpj
691	DO ji = 1, jpi
692	IF( Gsum(ji,jj,jl) < Gstar) THEN
693	athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
694	(2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
695	ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
696	athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * &
697	(2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
698	ELSE
699	athorn(ji,jj,jl) = 0.0
700	ENDIF
701	END DO ! ji
702	END DO ! jj
703	END DO ! jl
704
705	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
706	!
707	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
708
709	DO jl = -1, jpl
710	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
711	END DO !jl
712	DO jl = 0, jpl
713	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
714	END DO
715	!
716	ENDIF ! partfun_swi
717
718	IF( raft_swi == 1 ) THEN ! Ridging and rafting ice participation functions
719	!
720	DO jl = 1, jpl
721	DO jj = 1, jpj
722	DO ji = 1, jpi
723	IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN
724	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
725	aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
726	araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
727	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
728	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
729	ENDIF ! athorn
730	END DO ! ji
731	END DO ! jj
732	END DO ! jl
733
734	ELSE ! raft_swi = 0
735	!
736	DO jl = 1, jpl
737	aridge(:,:,jl) = athorn(:,:,jl)
738	END DO
739	!
740	ENDIF
741
742	IF ( raft_swi == 1 ) THEN
743
744	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
745	DO jl = 1, jpl
746	DO jj = 1, jpj
747	DO ji = 1, jpi
748	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
749	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
750	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
751	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
752	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
753	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
754	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
755	ENDIF
756	END DO
757	END DO
758	END DO
759	ENDIF
760
761	ENDIF
762
763	!-----------------------------------------------------------------
764	! 2) Transfer function
765	!-----------------------------------------------------------------
766	! Compute max and min ridged ice thickness for each ridging category.
767	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
768	!
769	! This parameterization is a modified version of Hibler (1980).
770	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
771	! and for very thick ridging ice must be >= krdgmin*hi
772	!
773	! The minimum ridging thickness, hrmin, is equal to 2*hi
774	! (i.e., rafting) and for very thick ridging ice is
775	! constrained by hrmin <= (hrmean + hi)/2.
776	!
777	! The maximum ridging thickness, hrmax, is determined by
778	! hrmean and hrmin.
779	!
780	! These modifications have the effect of reducing the ice strength
781	! (relative to the Hibler formulation) when very thick ice is
782	! ridging.
783	!
784	! aksum = net area removed/ total area removed
785	! where total area removed = area of ice that ridges
786	! net area removed = total area removed - area of new ridges
787	!-----------------------------------------------------------------
788
789	! Transfer function
790	DO jl = 1, jpl !all categories have a specific transfer function
791	DO jj = 1, jpj
792	DO ji = 1, jpi
793
794	IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
795	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
796	hrmean = MAX(SQRT(Hstarhi), hikrdgmin)
797	hrmin(ji,jj,jl) = MIN(2.0hi, 0.5(hrmean + hi))
798	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
799	hraft(ji,jj,jl) = kraft*hi
800	krdg(ji,jj,jl) = hrmean / hi
801	ELSE
802	hraft(ji,jj,jl) = 0.0
803	hrmin(ji,jj,jl) = 0.0
804	hrmax(ji,jj,jl) = 0.0
805	krdg (ji,jj,jl) = 1.0
806	ENDIF
807
808	END DO ! ji
809	END DO ! jj
810	END DO ! jl
811
812	! Normalization factor : aksum, ensures mass conservation
813	aksum(:,:) = athorn(:,:,0)
814	DO jl = 1, jpl
815	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
816	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
817	END DO
818	!
819	CALL wrk_dealloc( jpi,jpj, zworka )
820	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
821	!
822	END SUBROUTINE lim_itd_me_ridgeprep
823
824
825	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
826	!!----------------------------------------------------------------------
827	!! * ROUTINE lim_itd_me_icestrength *
828	!!
829	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
830	!!
831	!! ** Method : Remove area, volume, and energy from each ridging category
832	!! and add to thicker ice categories.
833	!!----------------------------------------------------------------------
834	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
835	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
836	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
837	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
838	!
839	CHARACTER (len=80) :: fieldid ! field identifier
840	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
841	!
842	LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny
843	LOGICAL :: large_afrac ! flag for afrac > 1
844	LOGICAL :: large_afrft ! flag for afrac > 1
845	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
846	INTEGER :: ij ! horizontal index, combines i and j loops
847	INTEGER :: icells ! number of cells with aicen > puny
848	REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration
849	REAL(wp) :: zsstK ! SST in Kelvin
850
851	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
852
853	REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories
854	REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers
855
856	REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging
857	REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging
858	REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging
859
860	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging
861
862	REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2
863	REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
864	REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
865	REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice
866	REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
867
868	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg1 ! volume of ice ridged
869	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg2 ! volume of new ridges
870	REAL(wp), POINTER, DIMENSION(:,:) :: vsw ! volume of seawater trapped into ridges
871	REAL(wp), POINTER, DIMENSION(:,:) :: srdg1 ! sal*volume of ice ridged
872	REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges
873	REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges
874
875	REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted
876	REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
877	REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice
878	REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice
879	REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice
880
881	REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice
882	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged
883	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges
884	REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges
885	!!----------------------------------------------------------------------
886
887	CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )
888	CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
889	CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
890	CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
891	CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
892	CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
893	CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
894	CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
895
896	! Conservation check
897	eice_init(:,:) = 0._wp
898
899	IF( con_i ) THEN
900	CALL lim_column_sum (jpl, v_i, vice_init )
901	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
902	DO ji = mi0(jiindx), mi1(jiindx)
903	DO jj = mj0(jjindx), mj1(jjindx)
904	WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)
905	WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)
906	END DO
907	END DO
908	ENDIF
909
910	!-------------------------------------------------------------------------------
911	! 1) Compute change in open water area due to closing and opening.
912	!-------------------------------------------------------------------------------
913
914	neg_ato_i = .false.
915
916	DO jj = 1, jpj
917	DO ji = 1, jpi
918	ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &
919	& + opning(ji,jj) * rdt_ice
920	IF( ato_i(ji,jj) < -epsi10 ) THEN
921	neg_ato_i = .TRUE.
922	ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error
923	ato_i(ji,jj) = 0._wp
924	ENDIF
925	END DO !jj
926	END DO !ji
927
928	! if negative open water area alert it
929	IF( neg_ato_i ) THEN ! there is a bug
930	DO jj = 1, jpj
931	DO ji = 1, jpi
932	IF( ato_i(ji,jj) < -epsi10 ) THEN
933	WRITE(numout,*) ''
934	WRITE(numout,*) 'Ridging error: ato_i < 0'
935	WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
936	ENDIF ! ato_i < -epsi10
937	END DO
938	END DO
939	ENDIF
940
941	!-----------------------------------------------------------------
942	! 2) Save initial state variables
943	!-----------------------------------------------------------------
944
945	DO jl = 1, jpl
946	aicen_init(:,:,jl) = a_i(:,:,jl)
947	vicen_init(:,:,jl) = v_i(:,:,jl)
948	vsnwn_init(:,:,jl) = v_s(:,:,jl)
949	!
950	smv_i_init(:,:,jl) = smv_i(:,:,jl)
951	oa_i_init (:,:,jl) = oa_i (:,:,jl)
952	END DO !jl
953
954	esnwn_init(:,:,:) = e_s(:,:,1,:)
955
956	DO jl = 1, jpl
957	DO jk = 1, nlay_i
958	eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
959	END DO
960	END DO
961
962	!
963	!-----------------------------------------------------------------
964	! 3) Pump everything from ice which is being ridged / rafted
965	!-----------------------------------------------------------------
966	! Compute the area, volume, and energy of ice ridging in each
967	! category, along with the area of the resulting ridge.
968
969	DO jl1 = 1, jpl !jl1 describes the ridging category
970
971	!------------------------------------------------
972	! 3.1) Identify grid cells with nonzero ridging
973	!------------------------------------------------
974
975	icells = 0
976	DO jj = 1, jpj
977	DO ji = 1, jpi
978	IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &
979	& .AND. closing_gross(ji,jj) > 0._wp ) THEN
980	icells = icells + 1
981	indxi(icells) = ji
982	indxj(icells) = jj
983	ENDIF ! test on a_icen_init
984	END DO ! ji
985	END DO ! jj
986
987	large_afrac = .false.
988	large_afrft = .false.
989
990	!CDIR NODEP
991	DO ij = 1, icells
992	ji = indxi(ij)
993	jj = indxj(ij)
994
995	!--------------------------------------------------------------------
996	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
997	!--------------------------------------------------------------------
998
999	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1000	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1001	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1002	arft2(ji,jj) = arft1(ji,jj) / kraft
1003
1004	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1005	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1006	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1007	oirft2(ji,jj)= oirft1(ji,jj) / kraft
1008
1009	!---------------------------------------------------------------
1010	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1011	!---------------------------------------------------------------
1012
1013	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1014	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1015
1016	IF (afrac(ji,jj) > kamax + epsi10) THEN !riging
1017	large_afrac = .true.
1018	ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error
1019	afrac(ji,jj) = kamax
1020	ENDIF
1021	IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
1022	large_afrft = .true.
1023	ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error
1024	afrft(ji,jj) = kamax
1025	ENDIF
1026
1027	!--------------------------------------------------------------------------
1028	! 3.4) Subtract area, volume, and energy from ridging
1029	! / rafting category n1.
1030	!--------------------------------------------------------------------------
1031	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
1032	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1033	vsw (ji,jj) = vrdg1(ji,jj) * ridge_por
1034
1035	vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
1036	esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
1037	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
1038	srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
1039
1040	! rafting volumes, heat contents ...
1041	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1042	vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
1043	esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
1044	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
1045
1046	! substract everything
1047	a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj)
1048	v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj)
1049	v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj)
1050	e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj)
1051	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
1052	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
1053
1054	!-----------------------------------------------------------------
1055	! 3.5) Compute properties of new ridges
1056	!-----------------------------------------------------------------
1057	!-------------
1058	! Salinity
1059	!-------------
1060	smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014
1061	srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge
1062
1063	!srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity
1064
1065	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
1066	wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! gurvan: increase in ice volume du to seawater frozen in voids
1067
1068	!------------------------------------
1069	! 3.6 Increment ridging diagnostics
1070	!------------------------------------
1071
1072	! jl1 looping 1-jpl
1073	! ij looping 1-icells
1074
1075	dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1076	dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1077	opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
1078
1079	IF( con_i ) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1080
1081	!------------------------------------------
1082	! 3.7 Put the snow somewhere in the ocean
1083	!------------------------------------------
1084	! Place part of the snow lost by ridging into the ocean.
1085	! Note that esnow_mlt < 0; the ocean must cool to melt snow.
1086	! If the ocean temp = Tf already, new ice must grow.
1087	! During the next time step, thermo_rates will determine whether
1088	! the ocean cools or new ice grows.
1089	! jl1 looping 1-jpl
1090	! ij looping 1-icells
1091
1092	msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosnvsrdg(ji,jj)(1.0-fsnowrdg) & ! rafting included
1093	& + rhosnvsrft(ji,jj)(1.0-fsnowrft)
1094
1095	! in J/m2 (same as e_s)
1096	esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included
1097	& - esrft(ji,jj)*(1.0-fsnowrft)
1098
1099	!-----------------------------------------------------------------
1100	! 3.8 Compute quantities used to apportion ice among categories
1101	! in the n2 loop below
1102	!-----------------------------------------------------------------
1103
1104	! jl1 looping 1-jpl
1105	! ij looping 1-icells
1106
1107	dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1108	dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1109
1110	END DO ! ij
1111
1112	!--------------------------------------------------------------------
1113	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1114	! compute ridged ice enthalpy
1115	!--------------------------------------------------------------------
1116	DO jk = 1, nlay_i
1117	!CDIR NODEP
1118	DO ij = 1, icells
1119	ji = indxi(ij)
1120	jj = indxj(ij)
1121	! heat content of ridged ice
1122	erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)
1123	eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1124	e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
1125
1126
1127	! enthalpy of the trapped seawater (J/m2, >0)
1128	! clem: if sst>0, then ersw <0 (is that possible?)
1129	zsstK = sst_m(ji,jj) + rt0
1130	ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i )
1131
1132	! heat flux to the ocean
1133	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
1134
1135	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
1136	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1137
1138	END DO ! ij
1139	END DO !jk
1140
1141
1142	IF( con_i ) THEN
1143	DO jk = 1, nlay_i
1144	!CDIR NODEP
1145	DO ij = 1, icells
1146	ji = indxi(ij)
1147	jj = indxj(ij)
1148	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
1149	END DO ! ij
1150	END DO !jk
1151	ENDIF
1152
1153	IF( large_afrac ) THEN ! there is a bug
1154	!CDIR NODEP
1155	DO ij = 1, icells
1156	ji = indxi(ij)
1157	jj = indxj(ij)
1158	IF( afrac(ji,jj) > kamax + epsi10 ) THEN
1159	WRITE(numout,*) ''
1160	WRITE(numout,*) ' ardg > a_i'
1161	WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
1162	ENDIF
1163	END DO
1164	ENDIF
1165	IF( large_afrft ) THEN ! there is a bug
1166	!CDIR NODEP
1167	DO ij = 1, icells
1168	ji = indxi(ij)
1169	jj = indxj(ij)
1170	IF( afrft(ji,jj) > kamax + epsi10 ) THEN
1171	WRITE(numout,*) ''
1172	WRITE(numout,*) ' arft > a_i'
1173	WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
1174	ENDIF
1175	END DO
1176	ENDIF
1177
1178	!-------------------------------------------------------------------------------
1179	! 4) Add area, volume, and energy of new ridge to each category jl2
1180	!-------------------------------------------------------------------------------
1181	! jl1 looping 1-jpl
1182	DO jl2 = 1, jpl
1183	! over categories to which ridged ice is transferred
1184	!CDIR NODEP
1185	DO ij = 1, icells
1186	ji = indxi(ij)
1187	jj = indxj(ij)
1188
1189	! Compute the fraction of ridged ice area and volume going to
1190	! thickness category jl2.
1191	! Transfer area, volume, and energy accordingly.
1192
1193	IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. &
1194	hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1195	hL = 0._wp
1196	hR = 0._wp
1197	ELSE
1198	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1199	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
1200	ENDIF
1201
1202	! fraction of ridged ice area and volume going to n2
1203	farea = ( hR - hL ) / dhr(ji,jj)
1204	fvol(ji,jj) = ( hRhR - hLhL ) / dhr2(ji,jj)
1205
1206	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea
1207	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1208	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1209	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
1210	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1211	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea
1212
1213	END DO ! ij
1214
1215	! Transfer ice energy to category jl2 by ridging
1216	DO jk = 1, nlay_i
1217	!CDIR NODEP
1218	DO ij = 1, icells
1219	ji = indxi(ij)
1220	jj = indxj(ij)
1221	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
1222	END DO
1223	END DO
1224	!
1225	END DO ! jl2 (new ridges)
1226
1227	DO jl2 = 1, jpl
1228
1229	!CDIR NODEP
1230	DO ij = 1, icells
1231	ji = indxi(ij)
1232	jj = indxj(ij)
1233	! Compute the fraction of rafted ice area and volume going to
1234	! thickness category jl2, transfer area, volume, and energy accordingly.
1235	!
1236	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
1237	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
1238	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj)
1239	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj)
1240	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * fsnowrft
1241	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft
1242	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj)
1243	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)
1244	ENDIF ! hraft
1245	!
1246	END DO ! ij
1247
1248	! Transfer rafted ice energy to category jl2
1249	DO jk = 1, nlay_i
1250	!CDIR NODEP
1251	DO ij = 1, icells
1252	ji = indxi(ij)
1253	jj = indxj(ij)
1254	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
1255	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
1256	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
1257	ENDIF
1258	END DO ! ij
1259	END DO !jk
1260
1261	END DO ! jl2
1262
1263	END DO ! jl1 (deforming categories)
1264
1265	! Conservation check
1266	IF ( con_i ) THEN
1267	CALL lim_column_sum (jpl, v_i, vice_final)
1268	fieldid = ' v_i : limitd_me '
1269	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
1270
1271	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
1272	fieldid = ' e_i : limitd_me '
1273	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
1274
1275	DO ji = mi0(jiindx), mi1(jiindx)
1276	DO jj = mj0(jjindx), mj1(jjindx)
1277	WRITE(numout,*) ' vice_init : ', vice_init (ji,jj)
1278	WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1279	WRITE(numout,*) ' eice_init : ', eice_init (ji,jj)
1280	WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1281	END DO
1282	END DO
1283	ENDIF
1284	!
1285	CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )
1286	CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
1287	CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1288	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
1289	CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1290	CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
1291	CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
1292	CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
1293	!
1294	END SUBROUTINE lim_itd_me_ridgeshift
1295
1296
1297	SUBROUTINE lim_itd_me_asumr
1298	!!-----------------------------------------------------------------------------
1299	!! * ROUTINE lim_itd_me_asumr *
1300	!!
1301	!! ** Purpose : finds total fractional area
1302	!!
1303	!! ** Method : Find the total area of ice plus open water in each grid cell.
1304	!! This is similar to the aggregate_area subroutine except that the
1305	!! total area can be greater than 1, so the open water area is
1306	!! included in the sum instead of being computed as a residual.
1307	!!-----------------------------------------------------------------------------
1308	INTEGER :: jl ! dummy loop index
1309	!!-----------------------------------------------------------------------------
1310	!
1311	asum(:,:) = ato_i(:,:) ! open water
1312	DO jl = 1, jpl ! ice categories
1313	asum(:,:) = asum(:,:) + a_i(:,:,jl)
1314	END DO
1315	!
1316	END SUBROUTINE lim_itd_me_asumr
1317
1318
1319	SUBROUTINE lim_itd_me_init
1320	!!-------------------------------------------------------------------
1321	!! * ROUTINE lim_itd_me_init *
1322	!!
1323	!! ** Purpose : Physical constants and parameters linked
1324	!! to the mechanical ice redistribution
1325	!!
1326	!! ** Method : Read the namiceitdme namelist
1327	!! and check the parameters values
1328	!! called at the first timestep (nit000)
1329	!!
1330	!! ** input : Namelist namiceitdme
1331	!!-------------------------------------------------------------------
1332	INTEGER :: ios ! Local integer output status for namelist read
1333	NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft, &
1334	& Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, &
1335	& partfun_swi, brinstren_swi
1336	!!-------------------------------------------------------------------
1337	!
1338	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
1339	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
1340	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
1341
1342	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
1343	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
1344	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
1345	IF(lwm) WRITE ( numoni, namiceitdme )
1346	!
1347	IF (lwp) THEN ! control print
1348	WRITE(numout,*)
1349	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1350	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1351	WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi
1352	WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs
1353	WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf
1354	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg
1355	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft
1356	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar
1357	WRITE(numout,)' Equivalent to G for an exponential part function astar ', astar
1358	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar
1359	WRITE(numout,*)' Rafting of ice sheets or not raft_swi ', raft_swi
1360	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter
1361	WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft
1362	WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por
1363	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi
1364	WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi
1365	ENDIF
1366	!
1367	END SUBROUTINE lim_itd_me_init
1368
1369	#else
1370	!!----------------------------------------------------------------------
1371	!! Default option Empty module NO LIM-3 sea-ice model
1372	!!----------------------------------------------------------------------
1373	CONTAINS
1374	SUBROUTINE lim_itd_me ! Empty routines
1375	END SUBROUTINE lim_itd_me
1376	SUBROUTINE lim_itd_me_icestrength
1377	END SUBROUTINE lim_itd_me_icestrength
1378	SUBROUTINE lim_itd_me_sort
1379	END SUBROUTINE lim_itd_me_sort
1380	SUBROUTINE lim_itd_me_init
1381	END SUBROUTINE lim_itd_me_init
1382	#endif
1383	!!======================================================================
1384	END MODULE limitd_me

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

Download in other formats: