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

Last change on this file since 834 was 834, checked in by ctlod, 16 years ago
Clean comments and useless lines, see ticket:#72
File size: 71.1 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	CALL prt_ctl(tab3d_1=a_i , clinfo1=' lim_itd_me: a_i : ', kdim=jpl)
209	CALL prt_ctl(tab3d_1=v_i , clinfo1=' lim_itd_me: v_i : ', kdim=jpl)
210	ENDIF
211
212	!-----------------------------------------------------------------------------!
213	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
214	!-----------------------------------------------------------------------------!
215	! Set hi_max(ncat) to a big value to ensure that all ridged ice
216	! is thinner than hi_max(ncat).
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	!-------------------------!
465	! Back to initial values
466	!-------------------------!
467
468	! update of fields will be made later in lim update
469	u_ice(:,:) = old_u_ice(:,:)
470	v_ice(:,:) = old_v_ice(:,:)
471	a_i(:,:,:) = old_a_i(:,:,:)
472	v_s(:,:,:) = old_v_s(:,:,:)
473	v_i(:,:,:) = old_v_i(:,:,:)
474	e_s(:,:,:,:) = old_e_s(:,:,:,:)
475	e_i(:,:,:,:) = old_e_i(:,:,:,:)
476	oa_i(:,:,:) = old_oa_i(:,:,:)
477	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
478	smv_i(:,:,:) = old_smv_i(:,:,:)
479
480
481	!----------------------------------------------------!
482	! Advection of ice in a free cell, newly ridged ice
483	!----------------------------------------------------!
484
485	! to allow for thermodynamics to melt new ice
486	! we immediately advect ice in free cells
487
488	! heat content has to be corrected before ice volume
489	DO jl = 1, jpl
490	DO jk = 1, nlay_i
491	DO jj = 1, jpj
492	DO ji = 1, jpi
493	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
494	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
495	old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl)
496	d_e_i_trp(ji,jj,jk,jl) = 0.0
497	ENDIF
498	END DO
499	END DO
500	END DO
501	END DO
502
503	DO jl = 1, jpl
504	DO jj = 1, jpj
505	DO ji = 1, jpi
506	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
507	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
508	old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl)
509	d_v_i_trp(ji,jj,jl) = 0.0
510	old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl)
511	d_a_i_trp(ji,jj,jl) = 0.0
512	old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl)
513	d_v_s_trp(ji,jj,jl) = 0.0
514	old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)
515	d_e_s_trp(ji,jj,1,jl) = 0.0
516	old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl)
517	d_oa_i_trp(ji,jj,jl) = 0.0
518	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
519	old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)
520	d_smv_i_trp(ji,jj,jl) = 0.0
521	ENDIF
522	END DO
523	END DO
524	END DO
525
526	END SUBROUTINE lim_itd_me
527
528	!===============================================================================
529
530	SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6)
531
532	!!----------------------------------------------------------------------
533	!! * ROUTINE lim_itd_me_icestrength *
534	!! ** Purpose :
535	!! This routine computes ice strength used in dynamics routines
536	!! of ice thickness
537	!!
538	!! ** Method :
539	!! Compute the strength of the ice pack, defined as the energy (J m-2)
540	!! dissipated per unit area removed from the ice pack under compression,
541	!! and assumed proportional to the change in potential energy caused
542	!! by ridging. Note that only Hibler's formulation is stable and that
543	!! ice strength has to be smoothed
544	!!
545	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
546	!!
547	!! ** External :
548	!!
549	!! ** References :
550	!!
551	!!----------------------------------------------------------------------
552	!! * Arguments
553
554	INTEGER, INTENT(in) :: &
555	kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
556
557	INTEGER :: &
558	ji,jj, & !: horizontal indices
559	jl, & !: thickness category index
560	ksmooth, & !: smoothing the resistance to deformation
561	numts_rm !: number of time steps for the P smoothing
562
563	REAL(wp) :: &
564	hi, & !: ice thickness (m)
565	zw1, & !: temporary variable
566	zp, & !: temporary ice strength
567	zdummy
568
569	REAL(wp), DIMENSION(jpi,jpj) :: &
570	zworka !: temporary array used here
571
572	!------------------------------------------------------------------------------!
573	! 1) Initialize
574	!------------------------------------------------------------------------------!
575	strength(:,:) = 0.0
576
577	!------------------------------------------------------------------------------!
578	! 2) Compute thickness distribution of ridging and ridged ice
579	!------------------------------------------------------------------------------!
580	CALL lim_itd_me_ridgeprep
581
582	!------------------------------------------------------------------------------!
583	! 3) Rothrock(1975)'s method
584	!------------------------------------------------------------------------------!
585	IF (kstrngth == 1) then
586
587	DO jl = 1, jpl
588	DO jj= 1, jpj
589	DO ji = 1, jpi
590
591	IF &
592	( ( a_i(ji,jj,jl) .GT. epsi11 ) &
593	.AND. ( athorn(ji,jj,jl) .GT. 0.0 ) ) THEN
594	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
595	!----------------------------
596	! PE loss from deforming ice
597	!----------------------------
598	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * &
599	hi * hi
600
601	!--------------------------
602	! PE gain from rafting ice
603	!--------------------------
604	strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) &
605	* hi * hi
606
607	!----------------------------
608	! PE gain from ridging ice
609	!----------------------------
610	strength(ji,jj) = strength(ji,jj) &
611	+ aridge(ji,jj,jl)/krdg(ji,jj,jl) &
612	* 1.0/3.0 * (hrmax(ji,jj,jl)3 - hrmin(ji,jj,jl)3) &
613	/ (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))
614	ENDIF ! aicen > epsi11
615
616	END DO ! ji
617	END DO !jj
618	END DO !jl
619
620	DO jj = 1, jpj
621	DO ji = 1, jpi
622	strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj)
623	! Cp = (g/2)(rhow-rhoi)(rhoi/rhow)
624	! Cf accounts for frictional dissipation
625
626	END DO ! j
627	END DO ! i
628
629	ksmooth = 1
630
631	!------------------------------------------------------------------------------!
632	! 4) Hibler (1979)' method
633	!------------------------------------------------------------------------------!
634	ELSE ! kstrngth ne 1: Hibler (1979) form
635
636	DO jj = 1, jpj
637	DO ji = 1, jpi
638	strength(ji,jj) = Pstarvt_i(ji,jj)exp(-C_rhg*(1.0-at_i(ji,jj)))
639	END DO ! j
640	END DO ! i
641
642	ksmooth = 1
643
644	ENDIF ! kstrngth
645
646	!
647	!------------------------------------------------------------------------------!
648	! 5) Impact of brine volume
649	!------------------------------------------------------------------------------!
650	! CAN BE REMOVED
651	!
652	IF ( brinstren_swi .EQ. 1 ) THEN
653
654	DO jj = 1, jpj
655	DO ji = 1, jpi
656	IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
657	zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
658	ELSE
659	zdummy = 0.0
660	ENDIF
661	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
662	END DO ! j
663	END DO ! i
664
665	ENDIF
666
667	!
668	!------------------------------------------------------------------------------!
669	! 6) Smoothing ice strength
670	!------------------------------------------------------------------------------!
671	!
672	!-------------------
673	! Spatial smoothing
674	!-------------------
675	IF ( ksmooth .EQ. 1 ) THEN
676
677	CALL lbc_lnk( strength, 'T', 1. )
678
679	DO jj = 1, jpj - 1
680	DO ji = 1, jpi - 1
681	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is
682	! present
683	zworka(ji,jj) = 4.0 * strength(ji,jj) &
684	+ strength(ji-1,jj) * tms(ji-1,jj) &
685	+ strength(ji+1,jj) * tms(ji+1,jj) &
686	+ strength(ji,jj-1) * tms(ji,jj-1) &
687	+ strength(ji,jj+1) * tms(ji,jj+1)
688
689	zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) &
690	+ tms(ji,jj-1) + tms(ji,jj+1)
691	zworka(ji,jj) = zworka(ji,jj) / zw1
692	ELSE
693	zworka(ji,jj) = 0.0
694	ENDIF
695	END DO
696	END DO
697
698	DO jj = 1, jpj - 1
699	DO ji = 1, jpi - 1
700	strength(ji,jj) = zworka(ji,jj)
701	END DO
702	END DO
703
704	ENDIF ! ksmooth
705
706	!--------------------
707	! Temporal smoothing
708	!--------------------
709	IF ( numit .EQ. nit000 + nfice - 1 ) THEN
710	strp1(:,:) = 0.0
711	strp2(:,:) = 0.0
712	ENDIF
713
714	IF ( ksmooth .EQ. 2 ) THEN
715
716
717	CALL lbc_lnk( strength, 'T', 1. )
718
719	DO jj = 1, jpj - 1
720	DO ji = 1, jpi - 1
721	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is
722	! present
723	numts_rm = 1 ! number of time steps for the running mean
724	IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
725	IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
726	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / &
727	numts_rm
728	strp2(ji,jj) = strp1(ji,jj)
729	strp1(ji,jj) = strength(ji,jj)
730	strength(ji,jj) = zp
731
732	ENDIF
733	END DO
734	END DO
735
736	ENDIF ! ksmooth
737
738	! Boundary conditions
739	CALL lbc_lnk( strength, 'T', 1. )
740
741	END SUBROUTINE lim_itd_me_icestrength
742
743	!===============================================================================
744
745	SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6)
746
747	!!---------------------------------------------------------------------!
748	!! * ROUTINE lim_itd_me_ridgeprep *
749	!! ** Purpose :
750	!! preparation for ridging and strength calculations
751	!!
752	!! ** Method :
753	!! Compute the thickness distribution of the ice and open water
754	!! participating in ridging and of the resulting ridges.
755	!!
756	!! ** Arguments :
757	!!
758	!! ** External :
759	!!
760	!!---------------------------------------------------------------------!
761	!! * Arguments
762
763	INTEGER :: &
764	ji,jj, & ! horizontal indices
765	jl, & ! thickness category index
766	krdg_index ! which participation function using
767
768	REAL(wp) :: &
769	Gstari, & ! = 1.0/Gstar
770	astari ! = 1.0/astar
771
772	REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: &
773	Gsum ! Gsum(n) = sum of areas in categories 0 to n
774
775	REAL(wp) :: &
776	hi, & ! ice thickness for each cat (m)
777	hrmean ! mean ridge thickness (m)
778
779	REAL(wp), DIMENSION(jpi,jpj) :: &
780	zworka ! temporary array used here
781
782	REAL(wp) :: &
783	zdummy, &
784	epsi06 = 1.0e-6
785
786	!------------------------------------------------------------------------------!
787
788	Gstari = 1.0/Gstar
789	astari = 1.0/astar
790	aksum(:,:) = 0.0
791	athorn(:,:,:) = 0.0
792	aridge(:,:,:) = 0.0
793	araft (:,:,:) = 0.0
794	hrmin(:,:,:) = 0.0
795	hrmax(:,:,:) = 0.0
796	hraft(:,:,:) = 0.0
797	krdg (:,:,:) = 1.0
798
799	! ! Zero out categories with very small areas
800	CALL lim_itd_me_zapsmall
801
802	!------------------------------------------------------------------------------!
803	! 1) Participation function
804	!------------------------------------------------------------------------------!
805
806	! Compute total area of ice plus open water.
807	CALL lim_itd_me_asumr
808	! This is in general not equal to one
809	! because of divergence during transport
810
811	! Compute cumulative thickness distribution function
812	! Compute the cumulative thickness distribution function Gsum,
813	! where Gsum(n) is the fractional area in categories 0 to n.
814	! initial value (in h = 0) equals open water area
815
816	Gsum(:,:,-1) = 0.0
817
818	DO jj = 1, jpj
819	DO ji = 1, jpi
820	IF (ato_i(ji,jj) .GT. epsi11) THEN
821	Gsum(ji,jj,0) = ato_i(ji,jj)
822	ELSE
823	Gsum(ji,jj,0) = 0.0
824	ENDIF
825	END DO
826	END DO
827
828	! for each value of h, you have to add ice concentration then
829	DO jl = 1, jpl
830	DO jj = 1, jpj
831	DO ji = 1, jpi
832	IF ( a_i(ji,jj,jl) .GT. epsi11 ) THEN
833	Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
834	ELSE
835	Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
836	ENDIF
837	END DO
838	END DO
839	END DO
840
841	! Normalize the cumulative distribution to 1
842	DO jj = 1, jpj
843	DO ji = 1, jpi
844	zworka(ji,jj) = 1.0 / Gsum(ji,jj,jpl)
845	END DO
846	END DO
847
848	DO jl = 0, jpl
849	DO jj = 1, jpj
850	DO ji = 1, jpi
851	Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj)
852	END DO
853	END DO
854	END DO
855
856	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
857	!--------------------------------------------------------------------------------------------------
858	! Compute the participation function athorn; this is analogous to
859	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
860	! area lost from category n due to ridging/closing
861	! athorn(n) = total area lost due to ridging/closing
862	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
863	!
864	! The expressions for athorn are found by integrating b(h)g(h) between
865	! the category boundaries.
866	!-----------------------------------------------------------------
867
868	krdg_index = 1
869
870	IF ( krdg_index .EQ. 0 ) THEN
871
872	!--- Linear formulation (Thorndike et al., 1975)
873	DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates
874	DO jj = 1, jpj
875	DO ji = 1, jpi
876	IF (Gsum(ji,jj,jl) < Gstar) THEN
877	athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
878	(2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
879	ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN
880	athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * &
881	(2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)
882	ELSE
883	athorn(ji,jj,jl) = 0.0
884	ENDIF
885	END DO ! ji
886	END DO ! jj
887	END DO ! jl
888
889	ELSE ! krdg_index = 1
890
891	!--- Exponential, more stable formulation (Lipscomb et al, 2007)
892	! precompute exponential terms using Gsum as a work array
893	zdummy = 1.0 / (1.0-EXP(-astari))
894
895	DO jl = -1, jpl
896	DO jj = 1, jpj
897	DO ji = 1, jpi
898	Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)astari)zdummy
899	END DO !ji
900	END DO !jj
901	END DO !jl
902
903	! compute athorn
904	DO jl = 0, ice_cat_bounds(1,2)
905	DO jj = 1, jpj
906	DO ji = 1, jpi
907	athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl)
908	END DO !ji
909	END DO ! jj
910	END DO !jl
911
912	ENDIF ! krdg_index
913
914	! Ridging and rafting ice participation functions
915	IF ( raftswi .EQ. 1 ) THEN
916
917	DO jl = 1, jpl
918	DO jj = 1, jpj
919	DO ji = 1, jpi
920	IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN
921	aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - &
922	hparmeter ) ) + 1.0 ) / 2.0 * &
923	athorn(ji,jj,jl)
924	araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - &
925	hparmeter ) ) + 1.0 ) / 2.0 * &
926	athorn(ji,jj,jl)
927	IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl) = 0.0
928	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0)
929	ENDIF ! athorn
930	END DO ! ji
931	END DO ! jj
932	END DO ! jl
933
934	ELSE ! raftswi = 0
935
936	DO jl = 1, jpl
937	DO jj = 1, jpj
938	DO ji = 1, jpi
939	aridge(ji,jj,jl) = 1.0*athorn(ji,jj,jl)
940	END DO
941	END DO
942	END DO
943
944	ENDIF
945
946	IF ( raftswi .EQ. 1 ) THEN
947
948	DO jl = 1, jpl
949	DO jj = 1, jpj
950	DO ji = 1, jpi
951	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. &
952	epsi11 ) THEN
953	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
954	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
955	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
956	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
957	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
958	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
959	ENDIF
960	END DO
961	END DO
962	END DO
963
964	ENDIF
965
966	!-----------------------------------------------------------------
967	! 2) Transfer function
968	!-----------------------------------------------------------------
969	! Compute max and min ridged ice thickness for each ridging category.
970	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
971	!
972	! This parameterization is a modified version of Hibler (1980).
973	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
974	! and for very thick ridging ice must be >= krdgmin*hi
975	!
976	! The minimum ridging thickness, hrmin, is equal to 2*hi
977	! (i.e., rafting) and for very thick ridging ice is
978	! constrained by hrmin <= (hrmean + hi)/2.
979	!
980	! The maximum ridging thickness, hrmax, is determined by
981	! hrmean and hrmin.
982	!
983	! These modifications have the effect of reducing the ice strength
984	! (relative to the Hibler formulation) when very thick ice is
985	! ridging.
986	!
987	! aksum = net area removed/ total area removed
988	! where total area removed = area of ice that ridges
989	! net area removed = total area removed - area of new ridges
990	!-----------------------------------------------------------------
991
992	! Transfer function
993	DO jl = 1, jpl !all categories have a specific transfer function
994	DO jj = 1, jpj
995	DO ji = 1, jpi
996
997	IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
998	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
999	hrmean = MAX(SQRT(Hstarhi), hikrdgmin)
1000	hrmin(ji,jj,jl) = MIN(2.0hi, 0.5(hrmean + hi))
1001	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
1002	hraft(ji,jj,jl) = kraft*hi
1003	krdg(ji,jj,jl) = hrmean / hi
1004	ELSE
1005	hraft(ji,jj,jl) = 0.0
1006	hrmin(ji,jj,jl) = 0.0
1007	hrmax(ji,jj,jl) = 0.0
1008	krdg (ji,jj,jl) = 1.0
1009	ENDIF
1010
1011	END DO ! ji
1012	END DO ! jj
1013	END DO ! jl
1014
1015	! Normalization factor : aksum, ensures mass conservation
1016	DO jj = 1, jpj
1017	DO ji = 1, jpi
1018	aksum(ji,jj) = athorn(ji,jj,0)
1019	END DO
1020	END DO
1021
1022	DO jl = 1, jpl
1023	DO jj = 1, jpj
1024	DO ji = 1, jpi
1025	aksum(ji,jj) = aksum(ji,jj) &
1026	+ aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) &
1027	+ araft (ji,jj,jl) * (1.0 - 1.0/kraft)
1028	END DO
1029	END DO
1030	END DO
1031
1032	END SUBROUTINE lim_itd_me_ridgeprep
1033
1034	!===============================================================================
1035
1036	SUBROUTINE lim_itd_me_ridgeshift(opning, closing_gross, &
1037	msnow_mlt, esnow_mlt) ! (subroutine 4/6)
1038
1039	!!-----------------------------------------------------------------------------
1040	!! * ROUTINE lim_itd_me_icestrength *
1041	!! ** Purpose :
1042	!! This routine shift ridging ice among thickness categories
1043	!! of ice thickness
1044	!!
1045	!! ** Method :
1046	!! Remove area, volume, and energy from each ridging category
1047	!! and add to thicker ice categories.
1048	!!
1049	!! ** Arguments :
1050	!!
1051	!! ** Inputs / Ouputs :
1052	!!
1053	!! ** External :
1054	!!
1055
1056	REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: &
1057	opning, & ! rate of opening due to divergence/shear
1058	closing_gross ! rate at which area removed, not counting
1059	! area of new ridges
1060
1061	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: &
1062	msnow_mlt, & ! mass of snow added to ocean (kg m-2)
1063	esnow_mlt ! energy needed to melt snow in ocean (J m-2)
1064
1065	INTEGER :: &
1066	ji, jj, & ! horizontal indices
1067	jl, jl1, jl2, & ! thickness category indices
1068	jk, & ! ice layer index
1069	ij, & ! horizontal index, combines i and j loops
1070	icells ! number of cells with aicen > puny
1071
1072	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
1073	indxi, indxj ! compressed indices
1074
1075	REAL(wp), DIMENSION(jpi,jpj) :: &
1076	vice_init, vice_final, & ! ice volume summed over categories
1077	eice_init, eice_final ! ice energy summed over layers
1078
1079	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
1080	aicen_init, & ! ice area before ridging
1081	vicen_init, & ! ice volume before ridging
1082	vsnon_init, & ! snow volume before ridging
1083	esnon_init, & ! snow energy before ridging
1084	smv_i_init, & ! ice salinity before ridging
1085	oa_i_init ! ice age before ridging
1086
1087	REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: &
1088	eicen_init ! ice energy before ridging
1089
1090	REAL(wp), DIMENSION(jpi,jpj) :: &
1091	afrac , & ! fraction of category area ridged
1092	ardg1 , & ! area of ice ridged
1093	ardg2 , & ! area of new ridges
1094	vsrdg , & ! snow volume of ridging ice
1095	esrdg , & ! snow energy of ridging ice
1096	oirdg1 , & ! areal age content of ridged ice
1097	oirdg2 , & ! areal age content of ridging ice
1098	dhr , & ! hrmax - hrmin
1099	dhr2 , & ! hrmax^2 - hrmin^2
1100	fvol ! fraction of new ridge volume going to n2
1101
1102	REAL(wp), DIMENSION(jpi,jpj) :: &
1103	vrdg1 , & ! volume of ice ridged
1104	vrdg2 , & ! volume of new ridges
1105	vsw , & ! volume of seawater trapped into ridges
1106	srdg1 , & ! sal*volume of ice ridged
1107	srdg2 , & ! sal*volume of new ridges
1108	smsw ! sal*volume of water trapped into ridges
1109
1110	REAL(wp), DIMENSION(jpi,jpj) :: &
1111	afrft , & ! fraction of category area rafted
1112	arft1 , & ! area of ice rafted
1113	arft2 , & ! area of new rafted zone
1114	virft , & ! ice volume of rafting ice
1115	vsrft , & ! snow volume of rafting ice
1116	esrft , & ! snow energy of rafting ice
1117	smrft , & ! salinity of rafting ice
1118	oirft1 , & ! areal age content of rafted ice
1119	oirft2 ! areal age content of rafting ice
1120
1121	REAL(wp), DIMENSION(jpi,jpj,jkmax) :: &
1122	eirft , & ! ice energy of rafting ice
1123	erdg1 , & ! enth*volume of ice ridged
1124	erdg2 , & ! enth*volume of new ridges
1125	ersw ! enth of water trapped into ridges
1126
1127	REAL(wp) :: &
1128	hL, hR , & ! left and right limits of integration
1129	farea , & ! fraction of new ridge area going to n2
1130	zdummy , & ! dummy argument
1131	zdummy0 , & ! dummy argument
1132	ztmelts , & ! ice melting point
1133	sm_newridge ! new ridged ice salinity
1134
1135	CHARACTER (len=80) :: &
1136	fieldid ! field identifier
1137
1138	LOGICAL, PARAMETER :: &
1139	l_conservation_check = .true. ! if true, check conservation
1140	! (useful for debugging)
1141	LOGICAL :: &
1142	neg_ato_i , & ! flag for ato_i(i,j) < -puny
1143	large_afrac , & ! flag for afrac > 1
1144	large_afrft ! flag for afrac > 1
1145
1146	REAL(wp) :: &
1147	zeps , &
1148	epsi10 , &
1149	zindb ! switch for the presence of ridge poros or not
1150
1151	!----------------------------------------------------------------------------
1152
1153	! Conservation check
1154	eice_init(:,:) = 0.0
1155
1156	IF ( con_i ) THEN
1157	CALL lim_column_sum (jpl, v_i, vice_init )
1158	WRITE(numout,*) ' vice_init : ', vice_init(jiindex,jjindex)
1159	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
1160	WRITE(numout,*) ' eice_init : ', eice_init(jiindex,jjindex)
1161	ENDIF
1162
1163	zeps = 1.0d-20
1164	epsi10 = 1.0d-10
1165
1166	!-------------------------------------------------------------------------------
1167	! 1) Compute change in open water area due to closing and opening.
1168	!-------------------------------------------------------------------------------
1169
1170	neg_ato_i = .false.
1171
1172	DO jj = 1, jpj
1173	DO ji = 1, jpi
1174	ato_i(ji,jj) = ato_i(ji,jj) &
1175	- athorn(ji,jj,0)closing_gross(ji,jj)rdt_ice &
1176	+ opning(ji,jj)*rdt_ice
1177	IF (ato_i(ji,jj) .LT. -epsi11) THEN
1178	neg_ato_i = .true.
1179	ELSEIF (ato_i(ji,jj) .LT. 0.0) THEN ! roundoff error
1180	ato_i(ji,jj) = 0.0
1181	ENDIF
1182	END DO !jj
1183	END DO !ji
1184
1185	! if negative open water area alert it
1186	IF (neg_ato_i) THEN ! there is a bug
1187	DO jj = 1, jpj
1188	DO ji = 1, jpi
1189	IF (ato_i(ji,jj) .LT. -epsi11) THEN
1190	WRITE(numout,*) ''
1191	WRITE(numout,*) 'Ridging error: ato_i < 0'
1192	WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
1193	ENDIF ! ato_i < -epsi11
1194	END DO ! ji
1195	END DO ! jj
1196	ENDIF ! neg_ato_i
1197
1198	!-----------------------------------------------------------------
1199	! 2) Save initial state variables
1200	!-----------------------------------------------------------------
1201
1202	DO jl = 1, jpl
1203	DO jj = 1, jpj
1204	DO ji = 1, jpi
1205	aicen_init(ji,jj,jl) = a_i(ji,jj,jl)
1206	vicen_init(ji,jj,jl) = v_i(ji,jj,jl)
1207	vsnon_init(ji,jj,jl) = v_s(ji,jj,jl)
1208
1209	esnon_init(ji,jj,jl) = e_s(ji,jj,1,jl)
1210	smv_i_init(ji,jj,jl) = smv_i(ji,jj,jl)
1211	oa_i_init (ji,jj,jl) = oa_i(ji,jj,jl)
1212	END DO !ji
1213	END DO ! jj
1214	END DO !jl
1215
1216	DO jl = 1, jpl
1217	DO jk = 1, nlay_i
1218	DO jj = 1, jpj
1219	DO ji = 1, jpi
1220	eicen_init(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)
1221	END DO !ji
1222	END DO !jj
1223	END DO !jk
1224	END DO !jl
1225
1226	!
1227	!-----------------------------------------------------------------
1228	! 3) Pump everything from ice which is being ridged / rafted
1229	!-----------------------------------------------------------------
1230	! Compute the area, volume, and energy of ice ridging in each
1231	! category, along with the area of the resulting ridge.
1232
1233	DO jl1 = 1, jpl !jl1 describes the ridging category
1234
1235	!------------------------------------------------
1236	! 3.1) Identify grid cells with nonzero ridging
1237	!------------------------------------------------
1238
1239	icells = 0
1240	DO jj = 1, jpj
1241	DO ji = 1, jpi
1242	IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0 &
1243	.AND. closing_gross(ji,jj) > 0.0) THEN
1244	icells = icells + 1
1245	indxi(icells) = ji
1246	indxj(icells) = jj
1247	ENDIF ! test on a_icen_init
1248	END DO ! ji
1249	END DO ! jj
1250
1251	large_afrac = .false.
1252	large_afrft = .false.
1253
1254	DO ij = 1, icells
1255	ji = indxi(ij)
1256	jj = indxj(ij)
1257
1258	!--------------------------------------------------------------------
1259	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
1260	!--------------------------------------------------------------------
1261
1262	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1263	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1264	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
1265	arft2(ji,jj) = arft1(ji,jj) / kraft
1266
1267	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1268	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
1269	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
1270	oirft2(ji,jj)= oirft1(ji,jj) / kraft
1271
1272	!---------------------------------------------------------------
1273	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
1274	!---------------------------------------------------------------
1275
1276	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
1277	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
1278
1279	IF (afrac(ji,jj) > 1.0 + epsi11) THEN !riging
1280	large_afrac = .true.
1281	ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error
1282	afrac(ji,jj) = 1.0
1283	ENDIF
1284	IF (afrft(ji,jj) > 1.0 + epsi11) THEN !rafting
1285	large_afrft = .true.
1286	ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error
1287	afrft(ji,jj) = 1.0
1288	ENDIF
1289
1290	!--------------------------------------------------------------------------
1291	! 3.4) Subtract area, volume, and energy from ridging
1292	! / rafting category n1.
1293	!--------------------------------------------------------------------------
1294	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / &
1295	( 1.0 + ridge_por )
1296	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
1297	vsw (ji,jj) = vrdg1(ji,jj) * ridge_por
1298
1299	vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
1300	esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
1301	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / &
1302	( 1. + ridge_por )
1303	srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
1304
1305	! rafting volumes, heat contents ...
1306	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
1307	vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
1308	esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
1309	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
1310
1311	! substract everything
1312	a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj)
1313	v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj)
1314	v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj)
1315	e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj)
1316	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
1317	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
1318
1319	!-----------------------------------------------------------------
1320	! 3.5) Compute properties of new ridges
1321	!-----------------------------------------------------------------
1322	!-------------
1323	! Salinity
1324	!-------------
1325	smsw(ji,jj) = sss_io(ji,jj) * vsw(ji,jj) * ridge_por
1326
1327	! salinity of new ridge
1328	sm_newridge = ( srdg1(ji,jj) + smsw(ji,jj) ) / vrdg2(ji,jj)
1329	zdummy = sm_newridge * vrdg2(ji,jj)
1330	! has to be smaller than s_i_max
1331	sm_newridge = MIN( s_i_max, sm_newridge )
1332
1333	! salt flux due to ridge creation
1334	fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + &
1335	MAX ( zdummy - srdg2(ji,jj) , 0.0 ) &
1336	* rhoic / rdt_ice
1337
1338	! sal times volume for new ridges
1339	srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj)
1340
1341	!------------------------------------
1342	! 3.6 Increment ridging diagnostics
1343	!------------------------------------
1344
1345	! jl1 looping 1-jpl
1346	! ij looping 1-icells
1347
1348	dardg1dt(ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1349	dardg2dt(ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1350	diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice
1351	opening(ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice
1352
1353	IF (con_i) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1354
1355	!------------------------------------------
1356	! 3.7 Put the snow somewhere in the ocean
1357	!------------------------------------------
1358
1359	! Place part of the snow lost by ridging into the ocean.
1360	! Note that esnow_mlt < 0; the ocean must cool to melt snow.
1361	! If the ocean temp = Tf already, new ice must grow.
1362	! During the next time step, thermo_rates will determine whether
1363	! the ocean cools or new ice grows.
1364	! jl1 looping 1-jpl
1365	! ij looping 1-icells
1366
1367	msnow_mlt(ji,jj) = msnow_mlt(ji,jj) &
1368	+ rhosnvsrdg(ji,jj)(1.0-fsnowrdg) &
1369	!rafting included
1370	+ rhosnvsrft(ji,jj)(1.0-fsnowrft)
1371
1372	esnow_mlt(ji,jj) = esnow_mlt(ji,jj) &
1373	+ esrdg(ji,jj)*(1.0-fsnowrdg) &
1374	!rafting included
1375	+ esrft(ji,jj)*(1.0-fsnowrft)
1376
1377	!-----------------------------------------------------------------
1378	! 3.8 Compute quantities used to apportion ice among categories
1379	! in the n2 loop below
1380	!-----------------------------------------------------------------
1381
1382	! jl1 looping 1-jpl
1383	! ij looping 1-icells
1384
1385	dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1386	dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) &
1387	- hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1388
1389
1390	END DO ! ij
1391
1392	!--------------------------------------------------------------------
1393	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1394	! compute ridged ice enthalpy
1395	!--------------------------------------------------------------------
1396	DO jk = 1, nlay_i
1397	DO ij = 1, icells
1398	ji = indxi(ij)
1399	jj = indxj(ij)
1400	! heat content of ridged ice
1401	erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / &
1402	( 1.0 + ridge_por )
1403	eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1404	e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) &
1405	- erdg1(ji,jj,jk) &
1406	- eirft(ji,jj,jk)
1407	! sea water heat content
1408	ztmelts = - tmut * sss_io(ji,jj) + rtt
1409	! heat content per unit volume
1410	zdummy0 = - rcp * ( sst_io(ji,jj) - rtt ) * vsw(ji,jj)
1411
1412	! corrected sea water salinity
1413	zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) )
1414	zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / &
1415	MAX( ridge_por * vsw(ji,jj), zeps )
1416
1417	ztmelts = - tmut * zdummy + rtt
1418	ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)
1419
1420	! heat flux
1421	fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / &
1422	rdt_ice
1423
1424	! Correct dimensions to avoid big values
1425	ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+09
1426
1427	! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
1428	ersw(ji,jj,jk) = ersw(ji,jj,jk) * &
1429	area(ji,jj) * vsw(ji,jj) / &
1430	nlay_i
1431
1432	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1433	END DO ! ij
1434	END DO !jk
1435
1436
1437	IF ( con_i ) THEN
1438	DO jk = 1, nlay_i
1439	DO ij = 1, icells
1440	ji = indxi(ij)
1441	jj = indxj(ij)
1442	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - &
1443	erdg1(ji,jj,jk)
1444	END DO ! ij
1445	END DO !jk
1446	ENDIF
1447
1448	IF (large_afrac) THEN ! there is a bug
1449	DO ij = 1, icells
1450	ji = indxi(ij)
1451	jj = indxj(ij)
1452	IF ( afrac(ji,jj) > 1.0 + epsi11 ) THEN
1453	WRITE(numout,*) ''
1454	WRITE(numout,*) ' ardg > a_i'
1455	WRITE(numout,*) ' ardg, aicen_init : ', &
1456	ardg1(ji,jj), aicen_init(ji,jj,jl1)
1457	ENDIF ! afrac > 1 + puny
1458	ENDDO ! if
1459	ENDIF ! large_afrac
1460	IF (large_afrft) THEN ! there is a bug
1461	DO ij = 1, icells
1462	ji = indxi(ij)
1463	jj = indxj(ij)
1464	IF ( afrft(ji,jj) > 1.0 + epsi11 ) THEN
1465	WRITE(numout,*) ''
1466	WRITE(numout,*) ' arft > a_i'
1467	WRITE(numout,*) ' arft, aicen_init : ', &
1468	arft1(ji,jj), aicen_init(ji,jj,jl1)
1469	ENDIF ! afrft > 1 + puny
1470	ENDDO ! if
1471	ENDIF ! large_afrft
1472
1473	!-------------------------------------------------------------------------------
1474	! 4) Add area, volume, and energy of new ridge to each category jl2
1475	!-------------------------------------------------------------------------------
1476	! jl1 looping 1-jpl
1477	DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
1478	! over categories to which ridged ice is transferred
1479	DO ij = 1, icells
1480	ji = indxi(ij)
1481	jj = indxj(ij)
1482
1483	! Compute the fraction of ridged ice area and volume going to
1484	! thickness category jl2.
1485	! Transfer area, volume, and energy accordingly.
1486
1487	IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR. &
1488	hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN
1489	hL = 0.0
1490	hR = 0.0
1491	ELSE
1492	hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1))
1493	hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2))
1494	ENDIF
1495
1496	! fraction of ridged ice area and volume going to n2
1497	farea = (hR-hL) / dhr(ji,jj)
1498	fvol(ji,jj) = (hRhR - hLhL) / dhr2(ji,jj)
1499
1500	a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + farea * ardg2(ji,jj)
1501	v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj)
1502	v_s(ji,jj,jl2) = v_s(ji,jj,jl2) &
1503	+ fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg
1504	e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) &
1505	+ fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg
1506	smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj)
1507	oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + farea * oirdg2(ji,jj)
1508
1509	END DO ! ij
1510
1511	! Transfer ice energy to category jl2 by ridging
1512	DO jk = 1, nlay_i
1513	DO ij = 1, icells
1514	ji = indxi(ij)
1515	jj = indxj(ij)
1516	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) &
1517	+ fvol(ji,jj)*erdg2(ji,jj,jk)
1518	END DO ! ij
1519	END DO !jk
1520
1521
1522	END DO ! jl2 (new ridges)
1523
1524	DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
1525
1526	DO ij = 1, icells
1527	ji = indxi(ij)
1528	jj = indxj(ij)
1529	! Compute the fraction of rafted ice area and volume going to
1530	! thickness category jl2, transfer area, volume, and energy accordingly.
1531
1532	IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. &
1533	hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
1534	a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj)
1535	v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj)
1536	v_s(ji,jj,jl2) = v_s(ji,jj,jl2) &
1537	+ vsrft(ji,jj)*fsnowrft
1538	e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) &
1539	+ esrft(ji,jj)*fsnowrft
1540	smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) &
1541	+ smrft(ji,jj)
1542	oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) &
1543	+ oirft2(ji,jj)
1544	ENDIF ! hraft
1545
1546	END DO ! ij
1547
1548	! Transfer rafted ice energy to category jl2
1549	DO jk = 1, nlay_i
1550	DO ij = 1, icells
1551	ji = indxi(ij)
1552	jj = indxj(ij)
1553	IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. &
1554	hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
1555	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) &
1556	+ eirft(ji,jj,jk)
1557	ENDIF
1558	END DO ! ij
1559	END DO !jk
1560
1561	END DO ! jl2
1562
1563	END DO ! jl1 (deforming categories)
1564
1565	! Conservation check
1566	IF ( con_i ) THEN
1567	CALL lim_column_sum (jpl, v_i, vice_final)
1568	fieldid = ' v_i : limitd_me '
1569	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
1570	WRITE(numout,*) ' vice_init : ', vice_init(jiindex,jjindex)
1571	WRITE(numout,*) ' vice_final : ', vice_final(jiindex,jjindex)
1572
1573	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
1574	fieldid = ' e_i : limitd_me '
1575	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
1576	WRITE(numout,*) ' eice_init : ', eice_init(jiindex,jjindex)
1577	WRITE(numout,*) ' eice_final : ', eice_final(jiindex,jjindex)
1578	ENDIF
1579
1580	END SUBROUTINE lim_itd_me_ridgeshift
1581
1582	!==============================================================================
1583
1584	SUBROUTINE lim_itd_me_asumr !(subroutine 5/6)
1585
1586	!!-----------------------------------------------------------------------------
1587	!! * ROUTINE lim_itd_me_asumr *
1588	!! ** Purpose :
1589	!! This routine finds total fractional area
1590	!!
1591	!! ** Method :
1592	!! Find the total area of ice plus open water in each grid cell.
1593	!!
1594	!! This is similar to the aggregate_area subroutine except that the
1595	!! total area can be greater than 1, so the open water area is
1596	!! included in the sum instead of being computed as a residual.
1597	!!
1598	!! ** Arguments :
1599
1600	INTEGER :: ji, jj, jl
1601
1602	!-----------------------------------------------------------------
1603	! open water
1604	!-----------------------------------------------------------------
1605
1606	DO jj = 1, jpj
1607	DO ji = 1, jpi
1608	asum(ji,jj) = ato_i(ji,jj)
1609	END DO
1610	END DO
1611
1612	!-----------------------------------------------------------------
1613	! ice categories
1614	!-----------------------------------------------------------------
1615
1616	DO jl = 1, jpl
1617	DO jj= 1, jpj
1618	DO ji = 1, jpi
1619	asum(ji,jj) = asum(ji,jj) + a_i(ji,jj,jl)
1620	END DO !ji
1621	END DO !jj
1622	END DO ! jl
1623
1624	END SUBROUTINE lim_itd_me_asumr
1625
1626	!==============================================================================
1627
1628	SUBROUTINE lim_itd_me_init ! (subroutine 6/6)
1629	!!-------------------------------------------------------------------
1630	!! * ROUTINE lim_itd_me_init *
1631	!!
1632	!! ** Purpose : Physical constants and parameters linked
1633	!! to the mechanical ice redistribution
1634	!!
1635	!! ** Method : Read the namiceitdme namelist
1636	!! and check the parameters values
1637	!! called at the first timestep (nit000)
1638	!!
1639	!! ** input : Namelist namiceitdme
1640	!!
1641	!! history :
1642	!! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1643	!!-------------------------------------------------------------------
1644	NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,&
1645	Gstar, astar, &
1646	Hstar, raftswi, hparmeter, Craft, ridge_por, &
1647	sal_max_ridge, partfun_swi, transfun_swi, &
1648	brinstren_swi
1649	!!-------------------------------------------------------------------
1650
1651	! Define the initial parameters
1652	! -------------------------
1653	REWIND( numnam_ice )
1654	READ ( numnam_ice , namiceitdme)
1655	IF (lwp) THEN
1656	WRITE(numout,*)
1657	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1658	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1659	WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi
1660	WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs
1661	WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf
1662	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg
1663	WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft
1664	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar
1665	WRITE(numout,)' Equivalent to G for an exponential part function astar ', astar
1666	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar
1667	WRITE(numout,*)' Rafting of ice sheets or not raftswi ', raftswi
1668	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter
1669	WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft
1670	WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por
1671	WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge
1672	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi
1673	WRITE(numout,*)' Switch for tran. function (0) linear (1) exponential transfun_swi ', transfun_swi
1674	WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi
1675	ENDIF
1676
1677	END SUBROUTINE lim_itd_me_init
1678
1679	!==============================================================================
1680
1681	SUBROUTINE lim_itd_me_zapsmall
1682	!!-------------------------------------------------------------------
1683	!! * ROUTINE lim_itd_me_zapsmall *
1684	!!
1685	!! ** Purpose : Remove too small sea ice areas and correct salt fluxes
1686	!!
1687	!!
1688	!! history :
1689	!! author: William H. Lipscomb, LANL
1690	!! Nov 2003: Modified by Julie Schramm to conserve volume and energy
1691	!! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
1692	!! additions to local freshwater, salt, and heat fluxes
1693	!! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code
1694	!!-------------------------------------------------------------------
1695
1696	INTEGER :: &
1697	ji,jj, & ! horizontal indices
1698	jl, & ! ice category index
1699	jk, & ! ice layer index
1700	icells, & ! number of cells with ice to zap
1701	ij ! combined i/j horizontal index
1702
1703	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
1704	indxi, & ! compressed indices for i/j directions
1705	indxj
1706
1707	REAL(wp) :: &
1708	xtmp ! temporary variable
1709
1710	DO jl = 1, jpl
1711
1712	!-----------------------------------------------------------------
1713	! Count categories to be zapped.
1714	! Abort model in case of negative area.
1715	!-----------------------------------------------------------------
1716
1717	icells = 0
1718	DO jj = 1, jpj
1719	DO ji = 1, jpi
1720	IF ( a_i(ji,jj,jl) .LT. -1.0e-11 ) THEN
1721	WRITE (numout,*) ' ALERTE 98 '
1722	WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl
1723	WRITE (numout,*) ' a_i ', a_i(ji,jj,jl)
1724	ELSEIF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) &
1725	.OR. &
1726	( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) &
1727	.OR. &
1728	!new line
1729	( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) &
1730	.OR. &
1731	( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN
1732	icells = icells + 1
1733	indxi(icells) = ji
1734	indxj(icells) = jj
1735	ENDIF
1736	END DO
1737	END DO
1738	WRITE(numout,*) icells, ' cells of ice zapped in the ocean '
1739
1740	!-----------------------------------------------------------------
1741	! Zap ice energy and use ocean heat to melt ice
1742	!-----------------------------------------------------------------
1743
1744	DO jk = 1, nlay_i
1745	DO ij = 1, icells
1746	ji = indxi(ij)
1747	jj = indxj(ij)
1748
1749	xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice
1750	xtmp = xtmp * unit_fac
1751	! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1752	e_i(ji,jj,jk,jl) = 0.0
1753
1754	END DO ! ij
1755	END DO ! jk
1756
1757	DO ij = 1, icells
1758	ji = indxi(ij)
1759	jj = indxj(ij)
1760
1761	!-----------------------------------------------------------------
1762	! Zap snow energy and use ocean heat to melt snow
1763	!-----------------------------------------------------------------
1764
1765	! xtmp = esnon(i,j,n) / dt ! < 0
1766	! fhnet(i,j) = fhnet(i,j) + xtmp
1767	! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
1768	! xtmp is greater than 0
1769	! fluxes are positive to the ocean
1770	! here the flux has to be negative for the ocean
1771	xtmp = ( rhosncpic( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice
1772	! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
1773
1774	t_s(ji,jj,1,jl) = rtt
1775
1776	!-----------------------------------------------------------------
1777	! zap ice and snow volume, add water and salt to ocean
1778	!-----------------------------------------------------------------
1779
1780	! xtmp = (rhoivicen(i,j,n) + rhosvsnon(i,j,n)) / dt
1781	! fresh(i,j) = fresh(i,j) + xtmp
1782	! fresh_hist(i,j) = fresh_hist(i,j) + xtmp
1783
1784	! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_io(ji,jj) ) * &
1785	! rhosn * v_s(ji,jj,jl) / rdt_ice
1786
1787	! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_io(ji,jj) - sm_i(ji,jj,jl) ) * &
1788	! rhoic * v_i(ji,jj,jl) / rdt_ice
1789
1790	! fsalt(i,j) = fsalt(i,j) + xtmp
1791	! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp
1792
1793	ato_i(ji,jj) = ato_i(ji,jj) + a_i(ji,jj,jl)
1794	a_i(ji,jj,jl) = 0.0
1795	v_i(ji,jj,jl) = 0.0
1796	v_s(ji,jj,jl) = 0.0
1797	t_su(ji,jj,jl) = t_bo(ji,jj)
1798	oa_i(ji,jj,jl) = 0.0
1799	smv_i(ji,jj,jl) = 0.0
1800
1801	END DO ! ij
1802
1803	END DO ! jl
1804
1805	END SUBROUTINE lim_itd_me_zapsmall
1806
1807	#else
1808	!!======================================================================
1809	!! * MODULE limitd_me *
1810	!! no sea ice model
1811	!!======================================================================
1812
1813	CONTAINS
1814
1815	SUBROUTINE lim_itd_me ! Empty routines
1816	END SUBROUTINE lim_itd_me
1817	SUBROUTINE lim_itd_me_icestrength
1818	END SUBROUTINE lim_itd_me_icestrength
1819	SUBROUTINE lim_itd_me_ridgeprep
1820	END SUBROUTINE lim_itd_me_ridgeprep
1821	SUBROUTINE lim_itd_th_me_ridgeshift
1822	END SUBROUTINE lim_itd_me_ridgeshift
1823	SUBROUTINE lim_itd_me_asumr
1824	END SUBROUTINE lim_itd_me_asumr
1825	SUBROUTINE lim_itd_me_sort
1826	END SUBROUTINE lim_itd_me_sort
1827	SUBROUTINE lim_itd_me_init
1828	END SUBROUTINE lim_itd_me_init
1829	SUBROUTINE lim_itd_me_zapsmall
1830	END SUBROUTINE lim_itd_me_zapsmall
1831	#endif
1832	END MODULE limitd_me

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

Download in other formats: