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

Last change on this file since 869 was 869, checked in by rblod, 16 years ago
Parallelisation of LIM3. This commit seems to ensure the reproducibility mono/mpp. See ticket #77.
File size: 73.8 KB

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

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

Download in other formats: