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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago
Update NEMOGCM from branch nemo_v3_3_beta
Property svn:keywords set to `Id`
File size: 74.6 KB

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

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

Download in other formats: