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

Last change on this file since 1264 was 1156, checked in by rblod, 16 years ago
Update Id and licence information, see ticket #210
Property svn:keywords set to `Id`
File size: 74.8 KB

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

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

Download in other formats: