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

Last change on this file since 868 was 868, checked in by rblod, 16 years ago
First optimisation of LIM3, limrhg optimisation induces computation change
File size: 73.7 KB

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

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

Download in other formats: