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

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

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

Download in other formats: