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

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

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

Download in other formats: