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

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

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

Download in other formats: