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 trunk/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limitd_me.F90 @ 867

Last change on this file since 867 was 867, checked in by ctlod, 16 years ago
Add control prints for sea-ice
File size: 72.8 KB

Line
1	MODULE limitd_me
2
3	#if defined key_lim3
4	!!----------------------------------------------------------------------
5	!! 'key_lim3' : LIM3 sea-ice model
6	!!----------------------------------------------------------------------
7	!!
8	!!======================================================================
9	!! * MODULE limitd_me *
10	!! Mechanical impact on ice thickness distribution
11	!! computation of changes in g(h)
12	!!======================================================================
13	!!
14	!!----------------------------------------------------------------------
15	!! * Modules used
16	USE dom_ice
17	USE par_oce ! ocean parameters
18	USE dom_oce
19	USE lbclnk
20	USE phycst ! physical constants (ocean directory)
21	USE ice_oce ! ice variables
22	USE thd_ice
23	USE limistate
24	USE in_out_manager
25	USE ice
26	USE par_ice
27	USE limthd_lac
28	USE limvar
29	USE iceini
30	USE limcons
31	USE prtctl ! Print control
32
33	IMPLICIT NONE
34	PRIVATE
35
36	!! * Routine accessibility
37	PUBLIC lim_itd_me ! called by ice_stp
38	PUBLIC lim_itd_me_icestrength
39	PUBLIC lim_itd_me_ridgeprep
40	PUBLIC lim_itd_me_ridgeshift
41	PUBLIC lim_itd_me_asumr
42	PUBLIC lim_itd_me_init
43	PUBLIC lim_itd_me_zapsmall
44
45	!! * Module variables
46	REAL(wp) :: & ! constant values
47	epsi20 = 1e-20 , &
48	epsi13 = 1e-13 , &
49	epsi11 = 1e-11 , &
50	zzero = 0.e0 , &
51	zone = 1.e0
52
53	!-----------------------------------------------------------------------
54	! Variables shared among ridging subroutines
55	!-----------------------------------------------------------------------
56	REAL(wp), DIMENSION (jpi,jpj) :: &
57	asum , & ! sum of total ice and open water area
58	aksum ! ratio of area removed to area ridged
59
60	REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &
61	athorn ! participation function; fraction of ridging/
62	! closing associated w/ category n
63
64	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
65	hrmin , & ! minimum ridge thickness
66	hrmax , & ! maximum ridge thickness
67	hraft , & ! thickness of rafted ice
68	krdg , & ! mean ridge thickness/thickness of ridging ice
69	aridge , & ! participating ice ridging
70	araft ! participating ice rafting
71
72	REAL(wp), PARAMETER :: &
73	krdgmin = 1.1, & ! min ridge thickness multiplier
74	kraft = 2.0 ! rafting multipliyer
75
76	REAL(wp) :: &
77	Cp
78	!
79	!-----------------------------------------------------------------------
80	! Ridging diagnostic arrays for history files
81	!-----------------------------------------------------------------------
82	!
83	REAL (wp), DIMENSION(jpi,jpj) :: &
84	dardg1dt , & ! rate of fractional area loss by ridging ice (1/s)
85	dardg2dt , & ! rate of fractional area gain by new ridges (1/s)
86	dvirdgdt , & ! rate of ice volume ridged (m/s)
87	opening ! rate of opening due to divergence/shear (1/s)
88
89
90	!!----------------------------------------------------------------------
91	!! LIM-@ 3.0, UCL-ASTR (2008)
92	!! (c) UCL-ASTR and Martin Vancoppenolle
93	!!----------------------------------------------------------------------
94
95	CONTAINS
96
97	!!-----------------------------------------------------------------------------!
98	!!-----------------------------------------------------------------------------!
99
100	SUBROUTINE lim_itd_me ! (subroutine 1/6)
101	!!---------------------------------------------------------------------!
102	!! * ROUTINE lim_itd_me *
103	!! ** Purpose :
104	!! This routine computes the mechanical redistribution
105	!! of ice thickness
106	!!
107	!! ** Method : a very simple method :-)
108	!!
109	!! ** Arguments :
110	!! kideb , kiut : Starting and ending points on which the
111	!! the computation is applied
112	!!
113	!! ** Inputs / Ouputs : (global commons)
114	!!
115	!! ** External :
116	!!
117	!! ** Steps :
118	!! 1) Thickness categories boundaries, ice / o.w. concentrations
119	!! Ridge preparation
120	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
121	!! 3) Ridging iteration
122	!! 4) Ridging diagnostics
123	!! 5) Heat, salt and freshwater fluxes
124	!! 6) Compute increments of tate variables and come back to old values
125	!!
126	!! ** References : There are a lot of references and can be difficult /
127	!! boring to read
128	!!
129	!! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength
130	!! in modeling the thickness distribution of Arctic sea ice,
131	!! J. Geophys. Res., 100, 18,611-18,626.
132	!!
133	!! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice
134	!! cover, Mon. Wea. Rev., 108, 1943-1973, 1980.
135	!!
136	!! Rothrock, D. A., 1975: The energetics of the plastic deformation of
137	!! pack ice by ridging, J. Geophys. Res., 80, 4514-4519.
138	!!
139	!! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,
140	!! 1975: The thickness distribution of sea ice, J. Geophys. Res.,
141	!! 80, 4501-4513.
142	!!
143	!! Bitz et al., JGR 2001
144	!!
145	!! Amundrud and Melling, JGR 2005
146	!!
147	!! Babko et al., JGR 2002
148	!!
149	!! ** History :
150	!! This routine is based on CICE code
151	!! and authors William H. Lipscomb,
152	!! and Elizabeth C. Hunke, LANL
153	!! are gratefully acknowledged
154	!!
155	!! (02-2006) Martin Vancoppenolle, UCL-ASTR
156	!!
157	!!--------------------------------------------------------------------!
158	!! * Arguments
159
160	!! * Local variables
161	INTEGER :: ji, & ! spatial dummy loop index
162	jj, & ! spatial dummy loop index
163	jk, & ! vertical layering dummy loop index
164	jl, & ! ice category dummy loop index
165	niter, & ! iteration counter
166	nitermax = 20 ! max number of ridging iterations
167
168	REAL(wp) :: & ! constant values
169	zeps = 1.0e-10, &
170	epsi10 = 1.0e-10, &
171	epsi06 = 1.0e-6
172
173	REAL(wp), DIMENSION(jpi,jpj) :: &
174	closing_net, & ! net rate at which area is removed (1/s)
175	! (ridging ice area - area of new ridges) / dt
176	divu_adv , & ! divu as implied by transport scheme (1/s)
177	opning , & ! rate of opening due to divergence/shear
178	closing_gross, & ! rate at which area removed, not counting
179	! area of new ridges
180	msnow_mlt , & ! mass of snow added to ocean (kg m-2)
181	esnow_mlt ! energy needed to melt snow in ocean (J m-2)
182
183	REAL(wp) :: &
184	w1, & ! temporary variable
185	tmpfac, & ! factor by which opening/closing rates are cut
186	dti ! 1 / dt
187
188	LOGICAL :: &
189	iterate_ridging, & ! if true, repeat the ridging
190	asum_error ! flag for asum .ne. 1
191
192	REAL(wp) :: &
193	big = 1.0e8
194
195	REAL (wp), DIMENSION(jpi,jpj) :: & !
196	vt_i_init, vt_i_final ! ice volume summed over categories
197
198	CHARACTER (len = 15) :: fieldid
199
200	!!-- End of declarations
201	!-----------------------------------------------------------------------------!
202
203	IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only)
204
205	IF(ln_ctl) THEN
206	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
207	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
208	ENDIF
209
210	!-----------------------------------------------------------------------------!
211	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
212	!-----------------------------------------------------------------------------!
213	! Set hi_max(ncat) to a big value to ensure that all ridged ice
214	! is thinner than hi_max(ncat).
215
216	hi_max(jpl) = 999.99
217
218	Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0 ! proport const for PE
219	CALL lim_itd_me_ridgeprep ! prepare ridging
220
221	! conservation check
222	IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init)
223
224	! Initialize arrays.
225	DO jj = 1, jpj
226	DO ji = 1, jpi
227
228	msnow_mlt(ji,jj) = 0.0
229	esnow_mlt(ji,jj) = 0.0
230	dardg1dt(ji,jj) = 0.0
231	dardg2dt(ji,jj) = 0.0
232	dvirdgdt(ji,jj) = 0.0
233	opening (ji,jj) = 0.0
234
235	!-----------------------------------------------------------------------------!
236	! 2) Dynamical inputs (closing rate, divu_adv, opning)
237	!-----------------------------------------------------------------------------!
238	!
239	! 2.1 closing_net
240	!-----------------
241	! Compute the net rate of closing due to convergence
242	! and shear, based on Flato and Hibler (1995).
243	!
244	! The energy dissipation rate is equal to the net closing rate
245	! times the ice strength.
246	!
247	! NOTE: The NET closing rate is equal to the rate that open water
248	! area is removed, plus the rate at which ice area is removed by
249	! ridging, minus the rate at which area is added in new ridges.
250	! The GROSS closing rate is equal to the first two terms (open
251	! water closing and thin ice ridging) without the third term
252	! (thick, newly ridged ice).
253
254	closing_net(ji,jj) = &
255	Cs0.5(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0)
256
257	! 2.2 divu_adv
258	!--------------
259	! Compute divu_adv, the divergence rate given by the transport/
260	! advection scheme, which may not be equal to divu as computed
261	! from the velocity field.
262	!
263	! If divu_adv < 0, make sure the closing rate is large enough
264	! to give asum = 1.0 after ridging.
265
266	divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep
267
268	IF (divu_adv(ji,jj) .LT. 0.0) &
269	closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj))
270
271	! 2.3 opning
272	!------------
273	! Compute the (non-negative) opening rate that will give
274	! asum = 1.0 after ridging.
275	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
276
277	END DO
278	END DO
279
280	!-----------------------------------------------------------------------------!
281	! 3) Ridging iteration
282	!-----------------------------------------------------------------------------!
283	niter = 1 ! iteration counter
284	iterate_ridging = .true.
285
286
287	DO WHILE ( iterate_ridging .AND. niter < nitermax )
288
289	DO jj = 1, jpj
290	DO ji = 1, jpi
291
292	! 3.2 closing_gross
293	!-----------------------------------------------------------------------------!
294	! Based on the ITD of ridging and ridged ice, convert the net
295	! closing rate to a gross closing rate.
296	! NOTE: 0 < aksum <= 1
297	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
298
299	! correction to closing rate and opening if closing rate is excessive
300	!---------------------------------------------------------------------
301	! Reduce the closing rate if more than 100% of the open water
302	! would be removed. Reduce the opening rate proportionately.
303	IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
304	w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
305	IF ( w1 .GT. ato_i(ji,jj)) THEN
306	tmpfac = ato_i(ji,jj) / w1
307	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
308	opning(ji,jj) = opning(ji,jj) * tmpfac
309	ENDIF !w1
310	ENDIF !at0i and athorn
311
312	END DO ! ji
313	END DO ! jj
314
315	! correction to closing rate / opening if excessive ice removal
316	!---------------------------------------------------------------
317	! Reduce the closing rate if more than 100% of any ice category
318	! would be removed. Reduce the opening rate proportionately.
319
320	DO jl = 1, jpl
321	DO jj = 1, jpj
322	DO ji = 1, jpi
323	IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
324	w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
325	IF ( w1 .GT. a_i(ji,jj,jl) ) THEN
326	tmpfac = a_i(ji,jj,jl) / w1
327	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
328	opning(ji,jj) = opning(ji,jj) * tmpfac
329	ENDIF
330	ENDIF
331	END DO !ji
332	END DO ! jj
333	END DO !jl
334
335	! 3.3 Redistribute area, volume, and energy.
336	!-----------------------------------------------------------------------------!
337
338	CALL lim_itd_me_ridgeshift (opning, closing_gross, &
339	msnow_mlt, esnow_mlt)
340
341	! 3.4 Compute total area of ice plus open water after ridging.
342	!-----------------------------------------------------------------------------!
343
344	CALL lim_itd_me_asumr
345
346	! 3.5 Do we keep on iterating ???
347	!-----------------------------------------------------------------------------!
348	! Check whether asum = 1. If not (because the closing and opening
349	! rates were reduced above), ridge again with new rates.
350
351	iterate_ridging = .false.
352
353	DO jj = 1, jpj
354	DO ji = 1, jpi
355	IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN
356	closing_net(ji,jj) = 0.0
357	opning(ji,jj) = 0.0
358	ELSE
359	iterate_ridging = .true.
360	divu_adv(ji,jj) = (1.0 - asum(ji,jj)) / rdt_ice
361	closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj))
362	opning(ji,jj) = MAX(0.0, divu_adv(ji,jj))
363	ENDIF
364	END DO
365	END DO
366
367	! Repeat if necessary.
368	! NOTE: If strength smoothing is turned on, the ridging must be
369	! iterated globally because of the boundary update in the
370	! smoothing.
371
372	niter = niter + 1
373
374	IF (iterate_ridging) THEN
375	IF (niter .GT. nitermax) THEN
376	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
377	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
378	ENDIF
379	CALL lim_itd_me_ridgeprep
380	ENDIF
381
382	END DO !! on the do while over iter
383
384	!-----------------------------------------------------------------------------!
385	! 4) Ridging diagnostics
386	!-----------------------------------------------------------------------------!
387	! Convert ridging rate diagnostics to correct units.
388	! Update fresh water and heat fluxes due to snow melt.
389
390	dti = 1.0/rdt_ice
391
392	asum_error = .false.
393
394	DO jj = 1, jpj
395	DO ji = 1, jpi
396
397	IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true.
398
399	dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti
400	dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti
401	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti
402	opening (ji,jj) = opening (ji,jj) * dti
403
404	!-----------------------------------------------------------------------------!
405	! 5) Heat, salt and freshwater fluxes
406	!-----------------------------------------------------------------------------!
407	! fresh water source for ocean
408	fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti
409
410	! heat sink for ocean
411	fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti
412
413	END DO
414	END DO
415
416	! Check if there is a ridging error
417	DO jj = 1, jpj
418	DO ji = 1, jpi
419	IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug
420	WRITE(numout,*) ' '
421	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
422	WRITE(numout,*) ' limitd_me '
423	WRITE(numout,*) ' POINT : ', ji, jj
424	WRITE(numout,*) ' jpl, a_i, athorn '
425	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
426	DO jl = 1, jpl
427	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
428	END DO
429	ENDIF ! asum
430
431	END DO !ji
432	END DO !jj
433
434	! Conservation check
435	IF ( con_i ) THEN
436	CALL lim_column_sum (jpl, v_i, vt_i_final)
437	fieldid = ' v_i : limitd_me '
438	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
439	ENDIF
440
441	!-----------------------------------------------------------------------------!
442	! 6) Updating state variables and trend terms
443	!-----------------------------------------------------------------------------!
444
445	CALL lim_var_glo2eqv
446	CALL lim_itd_me_zapsmall
447
448	!-----------------
449	! Trend terms
450	!-----------------
451
452	d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:)
453	d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:)
454	d_a_i_trp(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)
455	d_v_s_trp(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)
456	d_v_i_trp(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)
457	d_e_s_trp(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)
458	d_e_i_trp(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
459	d_oa_i_trp(:,:,:) = oa_i(:,:,:) - old_oa_i(:,:,:)
460	d_smv_i_trp(:,:,:) = 0.0
461	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
462	d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
463
464	IF(ln_ctl) THEN ! Control print
465	CALL prt_ctl_info(' ')
466	CALL prt_ctl_info(' - Cell values : ')
467	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
468	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me : cell area :')
469	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
470	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
471	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
472	DO jl = 1, jpl
473	CALL prt_ctl_info(' ')
474	CALL prt_ctl_info(' - Category : ', ivar1=jl)
475	CALL prt_ctl_info(' ~~~~~~~~~~')
476	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
477	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
478	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
479	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
480	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
481	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
482	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
483	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
484	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
485	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
486	DO jk = 1, nlay_i
487	CALL prt_ctl_info(' ')
488	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
489	CALL prt_ctl_info(' ~~~~~~~')
490	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
491	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
492	END DO
493	END DO
494	ENDIF
495
496	!-------------------------!
497	! Back to initial values
498	!-------------------------!
499
500	! update of fields will be made later in lim update
501	u_ice(:,:) = old_u_ice(:,:)
502	v_ice(:,:) = old_v_ice(:,:)
503	a_i(:,:,:) = old_a_i(:,:,:)
504	v_s(:,:,:) = old_v_s(:,:,:)
505	v_i(:,:,:) = old_v_i(:,:,:)
506	e_s(:,:,:,:) = old_e_s(:,:,:,:)
507	e_i(:,:,:,:) = old_e_i(:,:,:,:)
508	oa_i(:,:,:) = old_oa_i(:,:,:)
509	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
510	smv_i(:,:,:) = old_smv_i(:,:,:)
511
512	!----------------------------------------------------!
513	! Advection of ice in a free cell, newly ridged ice
514	!----------------------------------------------------!
515
516	! to allow for thermodynamics to melt new ice
517	! we immediately advect ice in free cells
518
519	! heat content has to be corrected before ice volume
520	DO jl = 1, jpl
521	DO jk = 1, nlay_i
522	DO jj = 1, jpj
523	DO ji = 1, jpi
524	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
525	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
526	old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl)
527	d_e_i_trp(ji,jj,jk,jl) = 0.0
528	ENDIF
529	END DO
530	END DO
531	END DO
532	END DO
533
534	DO jl = 1, jpl
535	DO jj = 1, jpj
536	DO ji = 1, jpi
537	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
538	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
539	old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl)
540	d_v_i_trp(ji,jj,jl) = 0.0
541	old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl)
542	d_a_i_trp(ji,jj,jl) = 0.0
543	old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl)
544	d_v_s_trp(ji,jj,jl) = 0.0
545	old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)
546	d_e_s_trp(ji,jj,1,jl) = 0.0
547	old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl)
548	d_oa_i_trp(ji,jj,jl) = 0.0
549	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
550	old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)
551	d_smv_i_trp(ji,jj,jl) = 0.0
552	ENDIF
553	END DO
554	END DO
555	END DO
556
557	END SUBROUTINE lim_itd_me
558
559	!===============================================================================
560
561	SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6)
562
563	!!----------------------------------------------------------------------
564	!! * ROUTINE lim_itd_me_icestrength *
565	!! ** Purpose :
566	!! This routine computes ice strength used in dynamics routines
567	!! of ice thickness
568	!!
569	!! ** Method :
570	!! Compute the strength of the ice pack, defined as the energy (J m-2)
571	!! dissipated per unit area removed from the ice pack under compression,
572	!! and assumed proportional to the change in potential energy caused
573	!! by ridging. Note that only Hibler's formulation is stable and that
574	!! ice strength has to be smoothed
575	!!
576	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
577	!!
578	!! ** External :
579	!!
580	!! ** References :
581	!!
582	!!----------------------------------------------------------------------
583	!! * Arguments
584
585	INTEGER, INTENT(in) :: &
586	kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
587
588	INTEGER :: &
589	ji,jj, & !: horizontal indices
590	jl, & !: thickness category index
591	ksmooth, & !: smoothing the resistance to deformation
592	numts_rm !: number of time steps for the P smoothing
593
594	REAL(wp) :: &
595	hi, & !: ice thickness (m)
596	zw1, & !: temporary variable
597	zp, & !: temporary ice strength
598	zdummy
599
600	REAL(wp), DIMENSION(jpi,jpj) :: &
601	zworka !: temporary array used here
602
603	!------------------------------------------------------------------------------!
604	! 1) Initialize
605	!------------------------------------------------------------------------------!
606	strength(:,:) = 0.0
607
608	!------------------------------------------------------------------------------!
609	! 2) Compute thickness distribution of ridging and ridged ice
610	!------------------------------------------------------------------------------!
611	CALL lim_itd_me_ridgeprep
612
613	!------------------------------------------------------------------------------!
614	! 3) Rothrock(1975)'s method
615	!------------------------------------------------------------------------------!
616	IF (kstrngth == 1) then
617
618	DO jl = 1, jpl
619	DO jj= 1, jpj
620	DO ji = 1, jpi
621
622	IF &
623	( ( a_i(ji,jj,jl) .GT. epsi11 ) &
624	.AND. ( athorn(ji,jj,jl) .GT. 0.0 ) ) THEN
625	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
626	!----------------------------
627	! PE loss from deforming ice
628	!----------------------------
629	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * &
630	hi * hi
631
632	!--------------------------
633	! PE gain from rafting ice
634	!--------------------------
635	strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) &
636	* hi * hi
637
638	!----------------------------
639	! PE gain from ridging ice
640	!----------------------------
641	strength(ji,jj) = strength(ji,jj) &
642	+ aridge(ji,jj,jl)/krdg(ji,jj,jl) &
643	* 1.0/3.0 * (hrmax(ji,jj,jl)3 - hrmin(ji,jj,jl)3) &
644	/ (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))
645	ENDIF ! aicen > epsi11
646
647	END DO ! ji
648	END DO !jj
649	END DO !jl
650
651	DO jj = 1, jpj
652	DO ji = 1, jpi
653	strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj)
654	! Cp = (g/2)(rhow-rhoi)(rhoi/rhow)
655	! Cf accounts for frictional dissipation
656
657	END DO ! j
658	END DO ! i
659
660	ksmooth = 1
661
662	!------------------------------------------------------------------------------!
663	! 4) Hibler (1979)' method
664	!------------------------------------------------------------------------------!
665	ELSE ! kstrngth ne 1: Hibler (1979) form
666
667	DO jj = 1, jpj
668	DO ji = 1, jpi
669	strength(ji,jj) = Pstarvt_i(ji,jj)exp(-C_rhg*(1.0-at_i(ji,jj)))
670	END DO ! j
671	END DO ! i
672
673	ksmooth = 1
674
675	ENDIF ! kstrngth
676
677	!
678	!------------------------------------------------------------------------------!
679	! 5) Impact of brine volume
680	!------------------------------------------------------------------------------!
681	! CAN BE REMOVED
682	!
683	IF ( brinstren_swi .EQ. 1 ) THEN
684
685	DO jj = 1, jpj
686	DO ji = 1, jpi
687	IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
688	zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
689	ELSE
690	zdummy = 0.0
691	ENDIF
692	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
693	END DO ! j
694	END DO ! i
695
696	ENDIF
697
698	!
699	!------------------------------------------------------------------------------!
700	! 6) Smoothing ice strength
701	!------------------------------------------------------------------------------!
702	!
703	!-------------------
704	! Spatial smoothing
705	!-------------------
706	IF ( ksmooth .EQ. 1 ) THEN
707
708	CALL lbc_lnk( strength, 'T', 1. )
709
710	DO jj = 1, jpj - 1
711	DO ji = 1, jpi - 1
712	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is
713	! present
714	zworka(ji,jj) = 4.0 * strength(ji,jj) &
715	+ strength(ji-1,jj) * tms(ji-1,jj) &
716	+ strength(ji+1,jj) * tms(ji+1,jj) &
717	+ strength(ji,jj-1) * tms(ji,jj-1) &
718	+ strength(ji,jj+1) * tms(ji,jj+1)
719
720	zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) &
721	+ tms(ji,jj-1) + tms(ji,jj+1)
722	zworka(ji,jj) = zworka(ji,jj) / zw1
723	ELSE
724	zworka(ji,jj) = 0.0
725	ENDIF
726	END DO
727	END DO
728
729	DO jj = 1, jpj - 1
730	DO ji = 1, jpi - 1
731	strength(ji,jj) = zworka(ji,jj)
732	END DO
733	END DO
734
735	ENDIF ! ksmooth
736
737	!--------------------
738	! Temporal smoothing
739	!--------------------
740	IF ( numit .EQ. nit000 + nfice - 1 ) THEN
741	strp1(:,:) = 0.0
742	strp2(:,:) = 0.0
743	ENDIF
744
745	IF ( ksmooth .EQ. 2 ) THEN
746
747
748	CALL lbc_lnk( strength, 'T', 1. )
749
750	DO jj = 1, jpj - 1
751	DO ji = 1, jpi - 1
752	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is
753	! present
754	numts_rm = 1 ! number of time steps for the running mean
755	IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
756	IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
757	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / &
758	numts_rm
759	strp2(ji,jj) = strp1(ji,jj)
760	strp1(ji,jj) = strength(ji,jj)
761	strength(ji,jj) = zp
762
763	ENDIF
764	END DO
765	END DO
766
767	ENDIF ! ksmooth
768
769	! Boundary conditions
770	CALL lbc_lnk( strength, 'T', 1. )
771
772	END SUBROUTINE lim_itd_me_icestrength
773
774	!===============================================================================
775
776	SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6)
777
778	!!---------------------------------------------------------------------!
779	!! * ROUTINE lim_itd_me_ridgeprep *
780	!! ** Purpose :
781	!! preparation for ridging and strength calculations
782	!!
783	!! ** Method :
784	!! Compute the thickness distribution of the ice and open water
785	!! participating in ridging and of the resulting ridges.
786	!!
787	!! ** Arguments :
788	!!
789	!! ** External :
790	!!
791	!!---------------------------------------------------------------------!
792	!! * Arguments
793
794	INTEGER :: &
795	ji,jj, & ! horizontal indices
796	jl, & ! thickness category index
797	krdg_index ! which participation function using
798
799	REAL(wp) :: &
800	Gstari, & ! = 1.0/Gstar
801	astari ! = 1.0/astar
802
803	REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: &
804	Gsum ! Gsum(n) = sum of areas in categories 0 to n
805
806	REAL(wp) :: &
807	hi, & ! ice thickness for each cat (m)
808	hrmean ! mean ridge thickness (m)
809
810	REAL(wp), DIMENSION(jpi,jpj) :: &
811	zworka ! temporary array used here
812
813	REAL(wp) :: &
814	zdummy, &
815	epsi06 = 1.0e-6
816
817	!------------------------------------------------------------------------------!
818
819	Gstari = 1.0/Gstar
820	astari = 1.0/astar
821	aksum(:,:) = 0.0
822	athorn(:,:,:) = 0.0
823	aridge(:,:,:) = 0.0
824	araft (:,:,:) = 0.0
825	hrmin(:,:,:) = 0.0
826	hrmax(:,:,:) = 0.0
827	hraft(:,:,:) = 0.0
828	krdg (:,:,:) = 1.0
829
830	! ! Zero out categories with very small areas
831	CALL lim_itd_me_zapsmall
832
833	!------------------------------------------------------------------------------!
834	! 1) Participation function
835	!------------------------------------------------------------------------------!
836
837	! Compute total area of ice plus open water.
838	CALL lim_itd_me_asumr
839	! This is in general not equal to one
840	! because of divergence during transport
841
842	! Compute cumulative thickness distribution function
843	! Compute the cumulative thickness distribution function Gsum,
844	! where Gsum(n) is the fractional area in categories 0 to n.
845	! initial value (in h = 0) equals open water area
846
847	Gsum(:,:,-1) = 0.0
848
849	DO jj = 1, jpj
850	DO ji = 1, jpi
851	IF (ato_i(ji,jj) .GT. epsi11) THEN
852	Gsum(ji,jj,0) = ato_i(ji,jj)
853	ELSE
854	Gsum(ji,jj,0) = 0.0
855	ENDIF
856	END DO
857	END DO
858
859	! for each value of h, you have to add ice concentration then
860	DO jl = 1, jpl
861	DO jj = 1, jpj
862	DO ji = 1, jpi
863	IF ( a_i(ji,jj,jl) .GT. epsi11 ) THEN
864	Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
865	ELSE
866	Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
867	ENDIF
868	END DO
869	END DO
870	END DO
871
872	! Normalize the cumulative distribution to 1
873	DO jj = 1, jpj
874	DO ji = 1, jpi
875	zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl)
876	END DO
877	END DO
878
879	DO jl = 0, jpl
880	DO jj = 1, jpj
881	DO ji = 1, jpi
882	Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj)
883	END DO
884	END DO
885	END DO
886
887	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
888	!--------------------------------------------------------------------------------------------------
889	! Compute the participation function athorn; this is analogous to
890	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
891	! area lost from category n due to ridging/closing
892	! athorn(n) = total area lost due to ridging/closing
893	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
894	!
895	! The expressions for athorn are found by integrating b(h)g(h) between
896	! the category boundaries.
897	!-----------------------------------------------------------------
898
899	krdg_index = 1
900
901	IF ( krdg_index .EQ. 0 ) THEN
902
903	!--- Linear formulation (Thorndike et al., 1975)
904	DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates
905	DO jj = 1, jpj
906	DO ji = 1, jpi
907	IF (Gsum(ji,jj,jl) < Gstar) THEN
908	athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
909	(2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
910	ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
911	athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * &
912	(2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
913	ELSE
914	athorn(ji,jj,jl) = 0.0
915	ENDIF
916	END DO ! ji
917	END DO ! jj
918	END DO ! jl
919
920	ELSE ! krdg_index = 1
921
922	!--- Exponential, more stable formulation (Lipscomb et al, 2007)
923	! precompute exponential terms using Gsum as a work array
924	zdummy = 1.0 / (1.0-EXP(-astari))
925
926	DO jl = -1, jpl
927	DO jj = 1, jpj
928	DO ji = 1, jpi
929	Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)astari)zdummy
930	END DO !ji
931	END DO !jj
932	END DO !jl
933
934	! compute athorn
935	DO jl = 0, ice_cat_bounds(1,2)
936	DO jj = 1, jpj
937	DO ji = 1, jpi
938	athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl)
939	END DO !ji
940	END DO ! jj
941	END DO !jl
942
943	ENDIF ! krdg_index
944
945	! Ridging and rafting ice participation functions
946	IF ( raftswi .EQ. 1 ) THEN
947
948	DO jl = 1, jpl
949	DO jj = 1, jpj
950	DO ji = 1, jpi
951	IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN
952	aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - &
953	hparmeter ) ) + 1.0 ) / 2.0 * &
954	athorn(ji,jj,jl)
955	araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - &
956	hparmeter ) ) + 1.0 ) / 2.0 * &
957	athorn(ji,jj,jl)
958	IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl) = 0.0
959	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0)
960	ENDIF ! athorn
961	END DO ! ji
962	END DO ! jj
963	END DO ! jl
964
965	ELSE ! raftswi = 0
966
967	DO jl = 1, jpl
968	DO jj = 1, jpj
969	DO ji = 1, jpi
970	aridge(ji,jj,jl) = 1.0*athorn(ji,jj,jl)
971	END DO
972	END DO
973	END DO
974
975	ENDIF
976
977	IF ( raftswi .EQ. 1 ) THEN
978
979	DO jl = 1, jpl
980	DO jj = 1, jpj
981	DO ji = 1, jpi
982	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. &
983	epsi11 ) THEN
984	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
985	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
986	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
987	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
988	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
989	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
990	ENDIF
991	END DO
992	END DO
993	END DO
994
995	ENDIF
996
997	!-----------------------------------------------------------------
998	! 2) Transfer function
999	!-----------------------------------------------------------------
1000	! Compute max and min ridged ice thickness for each ridging category.
1001	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
1002	!
1003	! This parameterization is a modified version of Hibler (1980).
1004	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
1005	! and for very thick ridging ice must be >= krdgmin*hi
1006	!
1007	! The minimum ridging thickness, hrmin, is equal to 2*hi
1008	! (i.e., rafting) and for very thick ridging ice is
1009	! constrained by hrmin <= (hrmean + hi)/2.
1010	!
1011	! The maximum ridging thickness, hrmax, is determined by
1012	! hrmean and hrmin.
1013	!
1014	! These modifications have the effect of reducing the ice strength
1015	! (relative to the Hibler formulation) when very thick ice is
1016	! ridging.
1017	!
1018	! aksum = net area removed/ total area removed
1019	! where total area removed = area of ice that ridges
1020	! net area removed = total area removed - area of new ridges
1021	!-----------------------------------------------------------------
1022
1023	! Transfer function
1024	DO jl = 1, jpl !all categories have a specific transfer function
1025	DO jj = 1, jpj
1026	DO ji = 1, jpi
1027
1028	IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
1029	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
1030	hrmean = MAX(SQRT(Hstarhi), hikrdgmin)
1031	hrmin(ji,jj,jl) = MIN(2.0hi, 0.5(hrmean + hi))
1032	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
1033	hraft(ji,jj,jl) = kraft*hi
1034	krdg(ji,jj,jl) = hrmean / hi
1035	ELSE
1036	hraft(ji,jj,jl) = 0.0
1037	hrmin(ji,jj,jl) = 0.0
1038	hrmax(ji,jj,jl) = 0.0
1039	krdg (ji,jj,jl) = 1.0
1040	ENDIF
1041
1042	END DO ! ji
1043	END DO ! jj
1044	END DO ! jl
1045
1046	! Normalization factor : aksum, ensures mass conservation
1047	DO jj = 1, jpj
1048	DO ji = 1, jpi
1049	aksum(ji,jj) = athorn(ji,jj,0)
1050	END DO
1051	END DO
1052
1053	DO jl = 1, jpl
1054	DO jj = 1, jpj
1055	DO ji = 1, jpi
1056	aksum(ji,jj) = aksum(ji,jj) &
1057	+ aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) &
1058	+ araft (ji,jj,jl) * (1.0 - 1.0/kraft)
1059	END DO
1060	END DO
1061	END DO
1062
1063	END SUBROUTINE lim_itd_me_ridgeprep
1064
1065	!===============================================================================
1066
1067	SUBROUTINE lim_itd_me_ridgeshift(opning, closing_gross, &
1068	msnow_mlt, esnow_mlt) ! (subroutine 4/6)
1069
1070	!!-----------------------------------------------------------------------------
1071	!! * ROUTINE lim_itd_me_icestrength *
1072	!! ** Purpose :
1073	!! This routine shift ridging ice among thickness categories
1074	!! of ice thickness
1075	!!
1076	!! ** Method :
1077	!! Remove area, volume, and energy from each ridging category
1078	!! and add to thicker ice categories.
1079	!!
1080	!! ** Arguments :
1081	!!
1082	!! ** Inputs / Ouputs :
1083	!!
1084	!! ** External :
1085	!!
1086
1087	REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: &
1088	opning, & ! rate of opening due to divergence/shear
1089	closing_gross ! rate at which area removed, not counting
1090	! area of new ridges
1091
1092	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: &
1093	msnow_mlt, & ! mass of snow added to ocean (kg m-2)
1094	esnow_mlt ! energy needed to melt snow in ocean (J m-2)
1095
1096	INTEGER :: &
1097	ji, jj, & ! horizontal indices
1098	jl, jl1, jl2, & ! thickness category indices
1099	jk, & ! ice layer index
1100	ij, & ! horizontal index, combines i and j loops
1101	icells ! number of cells with aicen > puny
1102
1103	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
1104	indxi, indxj ! compressed indices
1105
1106	REAL(wp), DIMENSION(jpi,jpj) :: &
1107	vice_init, vice_final, & ! ice volume summed over categories
1108	eice_init, eice_final ! ice energy summed over layers
1109
1110	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
1111	aicen_init, & ! ice area before ridging
1112	vicen_init, & ! ice volume before ridging
1113	vsnon_init, & ! snow volume before ridging
1114	esnon_init, & ! snow energy before ridging
1115	smv_i_init, & ! ice salinity before ridging
1116	oa_i_init ! ice age before ridging
1117
1118	REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: &
1119	eicen_init ! ice energy before ridging
1120
1121	REAL(wp), DIMENSION(jpi,jpj) :: &
1122	afrac , & ! fraction of category area ridged
1123	ardg1 , & ! area of ice ridged
1124	ardg2 , & ! area of new ridges
1125	vsrdg , & ! snow volume of ridging ice
1126	esrdg , & ! snow energy of ridging ice
1127	oirdg1 , & ! areal age content of ridged ice
1128	oirdg2 , & ! areal age content of ridging ice
1129	dhr , & ! hrmax - hrmin
1130	dhr2 , & ! hrmax^2 - hrmin^2
1131	fvol ! fraction of new ridge volume going to n2
1132
1133	REAL(wp), DIMENSION(jpi,jpj) :: &
1134	vrdg1 , & ! volume of ice ridged
1135	vrdg2 , & ! volume of new ridges
1136	vsw , & ! volume of seawater trapped into ridges
1137	srdg1 , & ! sal*volume of ice ridged
1138	srdg2 , & ! sal*volume of new ridges
1139	smsw ! sal*volume of water trapped into ridges
1140
1141	REAL(wp), DIMENSION(jpi,jpj) :: &
1142	afrft , & ! fraction of category area rafted
1143	arft1 , & ! area of ice rafted
1144	arft2 , & ! area of new rafted zone
1145	virft , & ! ice volume of rafting ice
1146	vsrft , & ! snow volume of rafting ice
1147	esrft , & ! snow energy of rafting ice
1148	smrft , & ! salinity of rafting ice
1149	oirft1 , & ! areal age content of rafted ice
1150	oirft2 ! areal age content of rafting ice
1151
1152	REAL(wp), DIMENSION(jpi,jpj,jkmax) :: &
1153	eirft , & ! ice energy of rafting ice
1154	erdg1 , & ! enth*volume of ice ridged
1155	erdg2 , & ! enth*volume of new ridges
1156	ersw ! enth of water trapped into ridges
1157
1158	REAL(wp) :: &
1159	hL, hR , & ! left and right limits of integration
1160	farea , & ! fraction of new ridge area going to n2
1161	zdummy , & ! dummy argument
1162	zdummy0 , & ! dummy argument
1163	ztmelts , & ! ice melting point
1164	sm_newridge ! new ridged ice salinity
1165
1166	CHARACTER (len=80) :: &
1167	fieldid ! field identifier
1168
1169	LOGICAL, PARAMETER :: &
1170	l_conservation_check = .true. ! if true, check conservation
1171	! (useful for debugging)
1172	LOGICAL :: &
1173	neg_ato_i , & ! flag for ato_i(i,j) < -puny
1174	large_afrac , & ! flag for afrac > 1
1175	large_afrft ! flag for afrac > 1
1176
1177	REAL(wp) :: &
1178	zeps , &
1179	epsi10 , &
1180	zindb ! switch for the presence of ridge poros or not
1181
1182	!----------------------------------------------------------------------------
1183
1184	! Conservation check
1185	eice_init(:,:) = 0.0
1186
1187	IF ( con_i ) THEN
1188	CALL lim_column_sum (jpl, v_i, vice_init )
1189	WRITE(numout,*) ' vice_init : ', vice_init(jiindex,jjindex)
1190	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
1191	WRITE(numout,*) ' eice_init : ', eice_init(jiindex,jjindex)
1192	ENDIF
1193
1194	zeps = 1.0d-20
1195	epsi10 = 1.0d-10
1196
1197	!-------------------------------------------------------------------------------
1198	! 1) Compute change in open water area due to closing and opening.
1199	!-------------------------------------------------------------------------------
1200
1201	neg_ato_i = .false.
1202
1203	DO jj = 1, jpj
1204	DO ji = 1, jpi
1205	ato_i(ji,jj) = ato_i(ji,jj) &
1206	- athorn(ji,jj,0)closing_gross(ji,jj)rdt_ice &
1207	+ opning(ji,jj)*rdt_ice
1208	IF (ato_i(ji,jj) .LT. -epsi11) THEN
1209	neg_ato_i = .true.
1210	ELSEIF (ato_i(ji,jj) .LT. 0.0) THEN ! roundoff error
1211	ato_i(ji,jj) = 0.0
1212	ENDIF
1213	END DO !jj
1214	END DO !ji
1215
1216	! if negative open water area alert it
1217	IF (neg_ato_i) THEN ! there is a bug
1218	DO jj = 1, jpj
1219	DO ji = 1, jpi
1220	IF (ato_i(ji,jj) .LT. -epsi11) THEN
1221	WRITE(numout,*) ''
1222	WRITE(numout,*) 'Ridging error: ato_i < 0'
1223	WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
1224	ENDIF ! ato_i < -epsi11
1225	END DO ! ji
1226	END DO ! jj
1227	ENDIF ! neg_ato_i
1228
1229	!-----------------------------------------------------------------
1230	! 2) Save initial state variables
1231	!-----------------------------------------------------------------
1232
1233	DO jl = 1, jpl
1234	DO jj = 1, jpj
1235	DO ji = 1, jpi
1236	aicen_init(ji,jj,jl) = a_i(ji,jj,jl)
1237	vicen_init(ji,jj,jl) = v_i(ji,jj,jl)
1238	vsnon_init(ji,jj,jl) = v_s(ji,jj,jl)
1239
1240	esnon_init(ji,jj,jl) = e_s(ji,jj,1,jl)
1241	smv_i_init(ji,jj,jl) = smv_i(ji,jj,jl)
1242	oa_i_init (ji,jj,jl) = oa_i(ji,jj,jl)
1243	END DO !ji
1244	END DO ! jj
1245	END DO !jl
1246
1247	DO jl = 1, jpl
1248	DO jk = 1, nlay_i
1249	DO jj = 1, jpj
1250	DO ji = 1, jpi
1251	eicen_init(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)
1252	END DO !ji
1253	END DO !jj
1254	END DO !jk
1255	END DO !jl
1256
1257	!
1258	!-----------------------------------------------------------------
1259	! 3) Pump everything from ice which is being ridged / rafted
1260	!-----------------------------------------------------------------
1261	! Compute the area, volume, and energy of ice ridging in each
1262	! category, along with the area of the resulting ridge.
1263
1264	DO jl1 = 1, jpl !jl1 describes the ridging category
1265
1266	!------------------------------------------------
1267	! 3.1) Identify grid cells with nonzero ridging
1268	!------------------------------------------------
1269
1270	icells = 0
1271	DO jj = 1, jpj
1272	DO ji = 1, jpi
1273	IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0 &
1274	.AND. closing_gross(ji,jj) > 0.0) THEN
1275	icells = icells + 1
1276	indxi(icells) = ji
1277	indxj(icells) = jj
1278	ENDIF ! test on a_icen_init
1279	END DO ! ji
1280	END DO ! jj
1281
1282	large_afrac = .false.
1283	large_afrft = .false.
1284
1285	DO ij = 1, icells
1286	ji = indxi(ij)
1287	jj = indxj(ij)
1288
1289	!--------------------------------------------------------------------
1290	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
1291	!--------------------------------------------------------------------
1292
1293	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1294	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1295	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1296	arft2(ji,jj) = arft1(ji,jj) / kraft
1297
1298	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1299	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1300	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1301	oirft2(ji,jj)= oirft1(ji,jj) / kraft
1302
1303	!---------------------------------------------------------------
1304	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1305	!---------------------------------------------------------------
1306
1307	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1308	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1309
1310	IF (afrac(ji,jj) > 1.0 + epsi11) THEN !riging
1311	large_afrac = .true.
1312	ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error
1313	afrac(ji,jj) = 1.0
1314	ENDIF
1315	IF (afrft(ji,jj) > 1.0 + epsi11) THEN !rafting
1316	large_afrft = .true.
1317	ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error
1318	afrft(ji,jj) = 1.0
1319	ENDIF
1320
1321	!--------------------------------------------------------------------------
1322	! 3.4) Subtract area, volume, and energy from ridging
1323	! / rafting category n1.
1324	!--------------------------------------------------------------------------
1325	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / &
1326	( 1.0 + ridge_por )
1327	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1328	vsw (ji,jj) = vrdg1(ji,jj) * ridge_por
1329
1330	vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
1331	esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
1332	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / &
1333	( 1. + ridge_por )
1334	srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
1335
1336	! rafting volumes, heat contents ...
1337	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1338	vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
1339	esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
1340	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
1341
1342	! substract everything
1343	a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj)
1344	v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj)
1345	v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj)
1346	e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj)
1347	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
1348	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
1349
1350	!-----------------------------------------------------------------
1351	! 3.5) Compute properties of new ridges
1352	!-----------------------------------------------------------------
1353	!-------------
1354	! Salinity
1355	!-------------
1356	smsw(ji,jj) = sss_io(ji,jj) * vsw(ji,jj) * ridge_por
1357
1358	! salinity of new ridge
1359	sm_newridge = ( srdg1(ji,jj) + smsw(ji,jj) ) / vrdg2(ji,jj)
1360	zdummy = sm_newridge * vrdg2(ji,jj)
1361	! has to be smaller than s_i_max
1362	sm_newridge = MIN( s_i_max, sm_newridge )
1363
1364	! salt flux due to ridge creation
1365	fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + &
1366	MAX ( zdummy - srdg2(ji,jj) , 0.0 ) &
1367	* rhoic / rdt_ice
1368
1369	! sal times volume for new ridges
1370	srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj)
1371
1372	!------------------------------------
1373	! 3.6 Increment ridging diagnostics
1374	!------------------------------------
1375
1376	! jl1 looping 1-jpl
1377	! ij looping 1-icells
1378
1379	dardg1dt(ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1380	dardg2dt(ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1381	diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice
1382	opening(ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice
1383
1384	IF (con_i) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1385
1386	!------------------------------------------
1387	! 3.7 Put the snow somewhere in the ocean
1388	!------------------------------------------
1389
1390	! Place part of the snow lost by ridging into the ocean.
1391	! Note that esnow_mlt < 0; the ocean must cool to melt snow.
1392	! If the ocean temp = Tf already, new ice must grow.
1393	! During the next time step, thermo_rates will determine whether
1394	! the ocean cools or new ice grows.
1395	! jl1 looping 1-jpl
1396	! ij looping 1-icells
1397
1398	msnow_mlt(ji,jj) = msnow_mlt(ji,jj) &
1399	+ rhosnvsrdg(ji,jj)(1.0-fsnowrdg) &
1400	!rafting included
1401	+ rhosnvsrft(ji,jj)(1.0-fsnowrft)
1402
1403	esnow_mlt(ji,jj) = esnow_mlt(ji,jj) &
1404	+ esrdg(ji,jj)*(1.0-fsnowrdg) &
1405	!rafting included
1406	+ esrft(ji,jj)*(1.0-fsnowrft)
1407
1408	!-----------------------------------------------------------------
1409	! 3.8 Compute quantities used to apportion ice among categories
1410	! in the n2 loop below
1411	!-----------------------------------------------------------------
1412
1413	! jl1 looping 1-jpl
1414	! ij looping 1-icells
1415
1416	dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1417	dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) &
1418	- hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1419
1420
1421	END DO ! ij
1422
1423	!--------------------------------------------------------------------
1424	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1425	! compute ridged ice enthalpy
1426	!--------------------------------------------------------------------
1427	DO jk = 1, nlay_i
1428	DO ij = 1, icells
1429	ji = indxi(ij)
1430	jj = indxj(ij)
1431	! heat content of ridged ice
1432	erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / &
1433	( 1.0 + ridge_por )
1434	eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1435	e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) &
1436	- erdg1(ji,jj,jk) &
1437	- eirft(ji,jj,jk)
1438	! sea water heat content
1439	ztmelts = - tmut * sss_io(ji,jj) + rtt
1440	! heat content per unit volume
1441	zdummy0 = - rcp * ( sst_io(ji,jj) - rtt ) * vsw(ji,jj)
1442
1443	! corrected sea water salinity
1444	zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) )
1445	zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / &
1446	MAX( ridge_por * vsw(ji,jj), zeps )
1447
1448	ztmelts = - tmut * zdummy + rtt
1449	ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)
1450
1451	! heat flux
1452	fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / &
1453	rdt_ice
1454
1455	! Correct dimensions to avoid big values
1456	ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+09
1457
1458	! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
1459	ersw(ji,jj,jk) = ersw(ji,jj,jk) * &
1460	area(ji,jj) * vsw(ji,jj) / &
1461	nlay_i
1462
1463	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1464	END DO ! ij
1465	END DO !jk
1466
1467
1468	IF ( con_i ) THEN
1469	DO jk = 1, nlay_i
1470	DO ij = 1, icells
1471	ji = indxi(ij)
1472	jj = indxj(ij)
1473	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - &
1474	erdg1(ji,jj,jk)
1475	END DO ! ij
1476	END DO !jk
1477	ENDIF
1478
1479	IF (large_afrac) THEN ! there is a bug
1480	DO ij = 1, icells
1481	ji = indxi(ij)
1482	jj = indxj(ij)
1483	IF ( afrac(ji,jj) > 1.0 + epsi11 ) THEN
1484	WRITE(numout,*) ''
1485	WRITE(numout,*) ' ardg > a_i'
1486	WRITE(numout,*) ' ardg, aicen_init : ', &
1487	ardg1(ji,jj), aicen_init(ji,jj,jl1)
1488	ENDIF ! afrac > 1 + puny
1489	ENDDO ! if
1490	ENDIF ! large_afrac
1491	IF (large_afrft) THEN ! there is a bug
1492	DO ij = 1, icells
1493	ji = indxi(ij)
1494	jj = indxj(ij)
1495	IF ( afrft(ji,jj) > 1.0 + epsi11 ) THEN
1496	WRITE(numout,*) ''
1497	WRITE(numout,*) ' arft > a_i'
1498	WRITE(numout,*) ' arft, aicen_init : ', &
1499	arft1(ji,jj), aicen_init(ji,jj,jl1)
1500	ENDIF ! afrft > 1 + puny
1501	ENDDO ! if
1502	ENDIF ! large_afrft
1503
1504	!-------------------------------------------------------------------------------
1505	! 4) Add area, volume, and energy of new ridge to each category jl2
1506	!-------------------------------------------------------------------------------
1507	! jl1 looping 1-jpl
1508	DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
1509	! over categories to which ridged ice is transferred
1510	DO ij = 1, icells
1511	ji = indxi(ij)
1512	jj = indxj(ij)
1513
1514	! Compute the fraction of ridged ice area and volume going to
1515	! thickness category jl2.
1516	! Transfer area, volume, and energy accordingly.
1517
1518	IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR. &
1519	hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN
1520	hL = 0.0
1521	hR = 0.0
1522	ELSE
1523	hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1))
1524	hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2))
1525	ENDIF
1526
1527	! fraction of ridged ice area and volume going to n2
1528	farea = (hR-hL) / dhr(ji,jj)
1529	fvol(ji,jj) = (hRhR - hLhL) / dhr2(ji,jj)
1530
1531	a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + farea * ardg2(ji,jj)
1532	v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj)
1533	v_s(ji,jj,jl2) = v_s(ji,jj,jl2) &
1534	+ fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg
1535	e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) &
1536	+ fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg
1537	smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj)
1538	oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + farea * oirdg2(ji,jj)
1539
1540	END DO ! ij
1541
1542	! Transfer ice energy to category jl2 by ridging
1543	DO jk = 1, nlay_i
1544	DO ij = 1, icells
1545	ji = indxi(ij)
1546	jj = indxj(ij)
1547	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) &
1548	+ fvol(ji,jj)*erdg2(ji,jj,jk)
1549	END DO ! ij
1550	END DO !jk
1551
1552
1553	END DO ! jl2 (new ridges)
1554
1555	DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
1556
1557	DO ij = 1, icells
1558	ji = indxi(ij)
1559	jj = indxj(ij)
1560	! Compute the fraction of rafted ice area and volume going to
1561	! thickness category jl2, transfer area, volume, and energy accordingly.
1562
1563	IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. &
1564	hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
1565	a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj)
1566	v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj)
1567	v_s(ji,jj,jl2) = v_s(ji,jj,jl2) &
1568	+ vsrft(ji,jj)*fsnowrft
1569	e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) &
1570	+ esrft(ji,jj)*fsnowrft
1571	smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) &
1572	+ smrft(ji,jj)
1573	oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) &
1574	+ oirft2(ji,jj)
1575	ENDIF ! hraft
1576
1577	END DO ! ij
1578
1579	! Transfer rafted ice energy to category jl2
1580	DO jk = 1, nlay_i
1581	DO ij = 1, icells
1582	ji = indxi(ij)
1583	jj = indxj(ij)
1584	IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. &
1585	hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
1586	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) &
1587	+ eirft(ji,jj,jk)
1588	ENDIF
1589	END DO ! ij
1590	END DO !jk
1591
1592	END DO ! jl2
1593
1594	END DO ! jl1 (deforming categories)
1595
1596	! Conservation check
1597	IF ( con_i ) THEN
1598	CALL lim_column_sum (jpl, v_i, vice_final)
1599	fieldid = ' v_i : limitd_me '
1600	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
1601	WRITE(numout,*) ' vice_init : ', vice_init(jiindex,jjindex)
1602	WRITE(numout,*) ' vice_final : ', vice_final(jiindex,jjindex)
1603
1604	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
1605	fieldid = ' e_i : limitd_me '
1606	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
1607	WRITE(numout,*) ' eice_init : ', eice_init(jiindex,jjindex)
1608	WRITE(numout,*) ' eice_final : ', eice_final(jiindex,jjindex)
1609	ENDIF
1610
1611	END SUBROUTINE lim_itd_me_ridgeshift
1612
1613	!==============================================================================
1614
1615	SUBROUTINE lim_itd_me_asumr !(subroutine 5/6)
1616
1617	!!-----------------------------------------------------------------------------
1618	!! * ROUTINE lim_itd_me_asumr *
1619	!! ** Purpose :
1620	!! This routine finds total fractional area
1621	!!
1622	!! ** Method :
1623	!! Find the total area of ice plus open water in each grid cell.
1624	!!
1625	!! This is similar to the aggregate_area subroutine except that the
1626	!! total area can be greater than 1, so the open water area is
1627	!! included in the sum instead of being computed as a residual.
1628	!!
1629	!! ** Arguments :
1630
1631	INTEGER :: ji, jj, jl
1632
1633	!-----------------------------------------------------------------
1634	! open water
1635	!-----------------------------------------------------------------
1636
1637	DO jj = 1, jpj
1638	DO ji = 1, jpi
1639	asum(ji,jj) = ato_i(ji,jj)
1640	END DO
1641	END DO
1642
1643	!-----------------------------------------------------------------
1644	! ice categories
1645	!-----------------------------------------------------------------
1646
1647	DO jl = 1, jpl
1648	DO jj= 1, jpj
1649	DO ji = 1, jpi
1650	asum(ji,jj) = asum(ji,jj) + a_i(ji,jj,jl)
1651	END DO !ji
1652	END DO !jj
1653	END DO ! jl
1654
1655	END SUBROUTINE lim_itd_me_asumr
1656
1657	!==============================================================================
1658
1659	SUBROUTINE lim_itd_me_init ! (subroutine 6/6)
1660	!!-------------------------------------------------------------------
1661	!! * ROUTINE lim_itd_me_init *
1662	!!
1663	!! ** Purpose : Physical constants and parameters linked
1664	!! to the mechanical ice redistribution
1665	!!
1666	!! ** Method : Read the namiceitdme namelist
1667	!! and check the parameters values
1668	!! called at the first timestep (nit000)
1669	!!
1670	!! ** input : Namelist namiceitdme
1671	!!
1672	!! history :
1673	!! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1674	!!-------------------------------------------------------------------
1675	NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,&
1676	Gstar, astar, &
1677	Hstar, raftswi, hparmeter, Craft, ridge_por, &
1678	sal_max_ridge, partfun_swi, transfun_swi, &
1679	brinstren_swi
1680	!!-------------------------------------------------------------------
1681
1682	! Define the initial parameters
1683	! -------------------------
1684	REWIND( numnam_ice )
1685	READ ( numnam_ice , namiceitdme)
1686	IF (lwp) THEN
1687	WRITE(numout,*)
1688	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1689	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1690	WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi
1691	WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs
1692	WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf
1693	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg
1694	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft
1695	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar
1696	WRITE(numout,)' Equivalent to G for an exponential part function astar ', astar
1697	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar
1698	WRITE(numout,*)' Rafting of ice sheets or not raftswi ', raftswi
1699	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter
1700	WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft
1701	WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por
1702	WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge
1703	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi
1704	WRITE(numout,*)' Switch for tran. function (0) linear (1) exponential transfun_swi ', transfun_swi
1705	WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi
1706	ENDIF
1707
1708	END SUBROUTINE lim_itd_me_init
1709
1710	!==============================================================================
1711
1712	SUBROUTINE lim_itd_me_zapsmall
1713	!!-------------------------------------------------------------------
1714	!! * ROUTINE lim_itd_me_zapsmall *
1715	!!
1716	!! ** Purpose : Remove too small sea ice areas and correct salt fluxes
1717	!!
1718	!!
1719	!! history :
1720	!! author: William H. Lipscomb, LANL
1721	!! Nov 2003: Modified by Julie Schramm to conserve volume and energy
1722	!! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
1723	!! additions to local freshwater, salt, and heat fluxes
1724	!! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1725	!!-------------------------------------------------------------------
1726
1727	INTEGER :: &
1728	ji,jj, & ! horizontal indices
1729	jl, & ! ice category index
1730	jk, & ! ice layer index
1731	icells, & ! number of cells with ice to zap
1732	ij ! combined i/j horizontal index
1733
1734	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
1735	indxi, & ! compressed indices for i/j directions
1736	indxj
1737
1738	REAL(wp) :: &
1739	xtmp ! temporary variable
1740
1741	DO jl = 1, jpl
1742
1743	!-----------------------------------------------------------------
1744	! Count categories to be zapped.
1745	! Abort model in case of negative area.
1746	!-----------------------------------------------------------------
1747
1748	icells = 0
1749	DO jj = 1, jpj
1750	DO ji = 1, jpi
1751	IF ( a_i(ji,jj,jl) .LT. -1.0e-11 ) THEN
1752	WRITE (numout,*) ' ALERTE 98 '
1753	WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl
1754	WRITE (numout,*) ' a_i ', a_i(ji,jj,jl)
1755	ELSEIF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) &
1756	.OR. &
1757	( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) &
1758	.OR. &
1759	!new line
1760	( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) &
1761	.OR. &
1762	( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN
1763	icells = icells + 1
1764	indxi(icells) = ji
1765	indxj(icells) = jj
1766	ENDIF
1767	END DO
1768	END DO
1769	WRITE(numout,*) icells, ' cells of ice zapped in the ocean '
1770
1771	!-----------------------------------------------------------------
1772	! Zap ice energy and use ocean heat to melt ice
1773	!-----------------------------------------------------------------
1774
1775	DO jk = 1, nlay_i
1776	DO ij = 1, icells
1777	ji = indxi(ij)
1778	jj = indxj(ij)
1779
1780	xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice
1781	xtmp = xtmp * unit_fac
1782	! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1783	e_i(ji,jj,jk,jl) = 0.0
1784
1785	END DO ! ij
1786	END DO ! jk
1787
1788	DO ij = 1, icells
1789	ji = indxi(ij)
1790	jj = indxj(ij)
1791
1792	!-----------------------------------------------------------------
1793	! Zap snow energy and use ocean heat to melt snow
1794	!-----------------------------------------------------------------
1795
1796	! xtmp = esnon(i,j,n) / dt ! < 0
1797	! fhnet(i,j) = fhnet(i,j) + xtmp
1798	! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
1799	! xtmp is greater than 0
1800	! fluxes are positive to the ocean
1801	! here the flux has to be negative for the ocean
1802	xtmp = ( rhosncpic( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice
1803	! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1804
1805	t_s(ji,jj,1,jl) = rtt
1806
1807	!-----------------------------------------------------------------
1808	! zap ice and snow volume, add water and salt to ocean
1809	!-----------------------------------------------------------------
1810
1811	! xtmp = (rhoivicen(i,j,n) + rhosvsnon(i,j,n)) / dt
1812	! fresh(i,j) = fresh(i,j) + xtmp
1813	! fresh_hist(i,j) = fresh_hist(i,j) + xtmp
1814
1815	! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_io(ji,jj) ) * &
1816	! rhosn * v_s(ji,jj,jl) / rdt_ice
1817
1818	! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_io(ji,jj) - sm_i(ji,jj,jl) ) * &
1819	! rhoic * v_i(ji,jj,jl) / rdt_ice
1820
1821	! fsalt(i,j) = fsalt(i,j) + xtmp
1822	! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp
1823
1824	ato_i(ji,jj) = ato_i(ji,jj) + a_i(ji,jj,jl)
1825	a_i(ji,jj,jl) = 0.0
1826	v_i(ji,jj,jl) = 0.0
1827	v_s(ji,jj,jl) = 0.0
1828	t_su(ji,jj,jl) = t_bo(ji,jj)
1829	oa_i(ji,jj,jl) = 0.0
1830	smv_i(ji,jj,jl) = 0.0
1831
1832	END DO ! ij
1833
1834	END DO ! jl
1835
1836	END SUBROUTINE lim_itd_me_zapsmall
1837
1838	#else
1839	!!======================================================================
1840	!! * MODULE limitd_me *
1841	!! no sea ice model
1842	!!======================================================================
1843
1844	CONTAINS
1845
1846	SUBROUTINE lim_itd_me ! Empty routines
1847	END SUBROUTINE lim_itd_me
1848	SUBROUTINE lim_itd_me_icestrength
1849	END SUBROUTINE lim_itd_me_icestrength
1850	SUBROUTINE lim_itd_me_ridgeprep
1851	END SUBROUTINE lim_itd_me_ridgeprep
1852	SUBROUTINE lim_itd_th_me_ridgeshift
1853	END SUBROUTINE lim_itd_me_ridgeshift
1854	SUBROUTINE lim_itd_me_asumr
1855	END SUBROUTINE lim_itd_me_asumr
1856	SUBROUTINE lim_itd_me_sort
1857	END SUBROUTINE lim_itd_me_sort
1858	SUBROUTINE lim_itd_me_init
1859	END SUBROUTINE lim_itd_me_init
1860	SUBROUTINE lim_itd_me_zapsmall
1861	END SUBROUTINE lim_itd_me_zapsmall
1862	#endif
1863	END MODULE limitd_me

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

Download in other formats: