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

Last change on this file since 1465 was 1465, checked in by smasson, 15 years ago
supress ice_oce module, see ticket:448
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 sbc_oce ! Surface boundary condition: ocean fields
22	USE thd_ice
23	USE limistate
24	USE in_out_manager
25	USE ice
26	USE par_ice
27	USE limthd_lac
28	USE limvar
29	USE iceini
30	USE limcons
31	USE prtctl ! Print control
32	USE lib_mpp
33
34	IMPLICIT NONE
35	PRIVATE
36
37	!! * Routine accessibility
38	PUBLIC lim_itd_me ! called by ice_stp
39	PUBLIC lim_itd_me_icestrength
40	PUBLIC lim_itd_me_ridgeprep
41	PUBLIC lim_itd_me_ridgeshift
42	PUBLIC lim_itd_me_asumr
43	PUBLIC lim_itd_me_init
44	PUBLIC lim_itd_me_zapsmall
45
46	!! * Module variables
47	REAL(wp) :: & ! constant values
48	epsi20 = 1e-20 , &
49	epsi13 = 1e-13 , &
50	epsi11 = 1e-11 , &
51	zzero = 0.e0 , &
52	zone = 1.e0
53
54	!-----------------------------------------------------------------------
55	! Variables shared among ridging subroutines
56	!-----------------------------------------------------------------------
57	REAL(wp), DIMENSION (jpi,jpj) :: &
58	asum , & ! sum of total ice and open water area
59	aksum ! ratio of area removed to area ridged
60
61	REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &
62	athorn ! participation function; fraction of ridging/
63	! closing associated w/ category n
64
65	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
66	hrmin , & ! minimum ridge thickness
67	hrmax , & ! maximum ridge thickness
68	hraft , & ! thickness of rafted ice
69	krdg , & ! mean ridge thickness/thickness of ridging ice
70	aridge , & ! participating ice ridging
71	araft ! participating ice rafting
72
73	REAL(wp), PARAMETER :: &
74	krdgmin = 1.1, & ! min ridge thickness multiplier
75	kraft = 2.0 ! rafting multipliyer
76
77	REAL(wp) :: &
78	Cp
79	!
80	!-----------------------------------------------------------------------
81	! Ridging diagnostic arrays for history files
82	!-----------------------------------------------------------------------
83	!
84	REAL (wp), DIMENSION(jpi,jpj) :: &
85	dardg1dt , & ! rate of fractional area loss by ridging ice (1/s)
86	dardg2dt , & ! rate of fractional area gain by new ridges (1/s)
87	dvirdgdt , & ! rate of ice volume ridged (m/s)
88	opening ! rate of opening due to divergence/shear (1/s)
89
90
91	!!----------------------------------------------------------------------
92	!! LIM 3.0, UCL-LOCEAN-IPSL (2008)
93	!! $Id$
94	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
95	!!----------------------------------------------------------------------
96
97
98	CONTAINS
99
100	!!-----------------------------------------------------------------------------!
101	!!-----------------------------------------------------------------------------!
102
103	SUBROUTINE lim_itd_me ! (subroutine 1/6)
104	!!---------------------------------------------------------------------!
105	!! * ROUTINE lim_itd_me *
106	!! ** Purpose :
107	!! This routine computes the mechanical redistribution
108	!! of ice thickness
109	!!
110	!! ** Method : a very simple method :-)
111	!!
112	!! ** Arguments :
113	!! kideb , kiut : Starting and ending points on which the
114	!! the computation is applied
115	!!
116	!! ** Inputs / Ouputs : (global commons)
117	!!
118	!! ** External :
119	!!
120	!! ** Steps :
121	!! 1) Thickness categories boundaries, ice / o.w. concentrations
122	!! Ridge preparation
123	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
124	!! 3) Ridging iteration
125	!! 4) Ridging diagnostics
126	!! 5) Heat, salt and freshwater fluxes
127	!! 6) Compute increments of tate variables and come back to old values
128	!!
129	!! ** References : There are a lot of references and can be difficult /
130	!! boring to read
131	!!
132	!! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength
133	!! in modeling the thickness distribution of Arctic sea ice,
134	!! J. Geophys. Res., 100, 18,611-18,626.
135	!!
136	!! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice
137	!! cover, Mon. Wea. Rev., 108, 1943-1973, 1980.
138	!!
139	!! Rothrock, D. A., 1975: The energetics of the plastic deformation of
140	!! pack ice by ridging, J. Geophys. Res., 80, 4514-4519.
141	!!
142	!! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,
143	!! 1975: The thickness distribution of sea ice, J. Geophys. Res.,
144	!! 80, 4501-4513.
145	!!
146	!! Bitz et al., JGR 2001
147	!!
148	!! Amundrud and Melling, JGR 2005
149	!!
150	!! Babko et al., JGR 2002
151	!!
152	!! ** History :
153	!! This routine is based on CICE code
154	!! and authors William H. Lipscomb,
155	!! and Elizabeth C. Hunke, LANL
156	!! are gratefully acknowledged
157	!!
158	!! (02-2006) Martin Vancoppenolle, UCL-ASTR
159	!!
160	!!--------------------------------------------------------------------!
161	!! * Arguments
162
163	!! * Local variables
164	INTEGER :: ji, & ! spatial dummy loop index
165	jj, & ! spatial dummy loop index
166	jk, & ! vertical layering dummy loop index
167	jl, & ! ice category dummy loop index
168	niter, & ! iteration counter
169	nitermax = 20 ! max number of ridging iterations
170
171	REAL(wp) :: & ! constant values
172	zeps = 1.0e-10, &
173	epsi10 = 1.0e-10, &
174	epsi06 = 1.0e-6
175
176	REAL(wp), DIMENSION(jpi,jpj) :: &
177	closing_net, & ! net rate at which area is removed (1/s)
178	! (ridging ice area - area of new ridges) / dt
179	divu_adv , & ! divu as implied by transport scheme (1/s)
180	opning , & ! rate of opening due to divergence/shear
181	closing_gross, & ! rate at which area removed, not counting
182	! area of new ridges
183	msnow_mlt , & ! mass of snow added to ocean (kg m-2)
184	esnow_mlt ! energy needed to melt snow in ocean (J m-2)
185
186	REAL(wp) :: &
187	w1, & ! temporary variable
188	tmpfac, & ! factor by which opening/closing rates are cut
189	dti ! 1 / dt
190
191	LOGICAL :: &
192	asum_error ! flag for asum .ne. 1
193
194	INTEGER :: iterate_ridging ! if true, repeat the ridging
195
196	REAL(wp) :: &
197	big = 1.0e8
198
199	REAL (wp), DIMENSION(jpi,jpj) :: & !
200	vt_i_init, vt_i_final ! ice volume summed over categories
201
202	CHARACTER (len = 15) :: fieldid
203
204	!!-- End of declarations
205	!-----------------------------------------------------------------------------!
206
207	IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only)
208
209	IF(ln_ctl) THEN
210	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
211	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
212	ENDIF
213
214	!-----------------------------------------------------------------------------!
215	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
216	!-----------------------------------------------------------------------------!
217	! Set hi_max(ncat) to a big value to ensure that all ridged ice
218	! is thinner than hi_max(ncat).
219
220	hi_max(jpl) = 999.99
221
222	Cp = 0.5* grav * (rau0-rhoic)*rhoic/rau0 ! proport const for PE
223	CALL lim_itd_me_ridgeprep ! prepare ridging
224
225	! conservation check
226	IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init)
227
228	! Initialize arrays.
229	DO jj = 1, jpj
230	DO ji = 1, jpi
231
232	msnow_mlt(ji,jj) = 0.0
233	esnow_mlt(ji,jj) = 0.0
234	dardg1dt(ji,jj) = 0.0
235	dardg2dt(ji,jj) = 0.0
236	dvirdgdt(ji,jj) = 0.0
237	opening (ji,jj) = 0.0
238
239	!-----------------------------------------------------------------------------!
240	! 2) Dynamical inputs (closing rate, divu_adv, opning)
241	!-----------------------------------------------------------------------------!
242	!
243	! 2.1 closing_net
244	!-----------------
245	! Compute the net rate of closing due to convergence
246	! and shear, based on Flato and Hibler (1995).
247	!
248	! The energy dissipation rate is equal to the net closing rate
249	! times the ice strength.
250	!
251	! NOTE: The NET closing rate is equal to the rate that open water
252	! area is removed, plus the rate at which ice area is removed by
253	! ridging, minus the rate at which area is added in new ridges.
254	! The GROSS closing rate is equal to the first two terms (open
255	! water closing and thin ice ridging) without the third term
256	! (thick, newly ridged ice).
257
258	closing_net(ji,jj) = &
259	Cs0.5(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0)
260
261	! 2.2 divu_adv
262	!--------------
263	! Compute divu_adv, the divergence rate given by the transport/
264	! advection scheme, which may not be equal to divu as computed
265	! from the velocity field.
266	!
267	! If divu_adv < 0, make sure the closing rate is large enough
268	! to give asum = 1.0 after ridging.
269
270	divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep
271
272	IF (divu_adv(ji,jj) .LT. 0.0) &
273	closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj))
274
275	! 2.3 opning
276	!------------
277	! Compute the (non-negative) opening rate that will give
278	! asum = 1.0 after ridging.
279	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
280
281	END DO
282	END DO
283
284	!-----------------------------------------------------------------------------!
285	! 3) Ridging iteration
286	!-----------------------------------------------------------------------------!
287	niter = 1 ! iteration counter
288	iterate_ridging = 1
289
290
291	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
292
293	DO jj = 1, jpj
294	DO ji = 1, jpi
295
296	! 3.2 closing_gross
297	!-----------------------------------------------------------------------------!
298	! Based on the ITD of ridging and ridged ice, convert the net
299	! closing rate to a gross closing rate.
300	! NOTE: 0 < aksum <= 1
301	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
302
303	! correction to closing rate and opening if closing rate is excessive
304	!---------------------------------------------------------------------
305	! Reduce the closing rate if more than 100% of the open water
306	! would be removed. Reduce the opening rate proportionately.
307	IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
308	w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
309	IF ( w1 .GT. ato_i(ji,jj)) THEN
310	tmpfac = ato_i(ji,jj) / w1
311	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
312	opning(ji,jj) = opning(ji,jj) * tmpfac
313	ENDIF !w1
314	ENDIF !at0i and athorn
315
316	END DO ! ji
317	END DO ! jj
318
319	! correction to closing rate / opening if excessive ice removal
320	!---------------------------------------------------------------
321	! Reduce the closing rate if more than 100% of any ice category
322	! would be removed. Reduce the opening rate proportionately.
323
324	DO jl = 1, jpl
325	DO jj = 1, jpj
326	DO ji = 1, jpi
327	IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
328	w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
329	IF ( w1 .GT. a_i(ji,jj,jl) ) THEN
330	tmpfac = a_i(ji,jj,jl) / w1
331	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
332	opning(ji,jj) = opning(ji,jj) * tmpfac
333	ENDIF
334	ENDIF
335	END DO !ji
336	END DO ! jj
337	END DO !jl
338
339	! 3.3 Redistribute area, volume, and energy.
340	!-----------------------------------------------------------------------------!
341
342	CALL lim_itd_me_ridgeshift (opning, closing_gross, &
343	msnow_mlt, esnow_mlt)
344
345	! 3.4 Compute total area of ice plus open water after ridging.
346	!-----------------------------------------------------------------------------!
347
348	CALL lim_itd_me_asumr
349
350	! 3.5 Do we keep on iterating ???
351	!-----------------------------------------------------------------------------!
352	! Check whether asum = 1. If not (because the closing and opening
353	! rates were reduced above), ridge again with new rates.
354
355	iterate_ridging = 0
356
357	DO jj = 1, jpj
358	DO ji = 1, jpi
359	IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN
360	closing_net(ji,jj) = 0.0
361	opning(ji,jj) = 0.0
362	ELSE
363	iterate_ridging = 1
364	divu_adv(ji,jj) = (1.0 - asum(ji,jj)) / rdt_ice
365	closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj))
366	opning(ji,jj) = MAX(0.0, divu_adv(ji,jj))
367	ENDIF
368	END DO
369	END DO
370
371	IF( lk_mpp ) CALL mpp_max(iterate_ridging)
372
373	! Repeat if necessary.
374	! NOTE: If strength smoothing is turned on, the ridging must be
375	! iterated globally because of the boundary update in the
376	! smoothing.
377
378	niter = niter + 1
379
380	IF (iterate_ridging == 1) THEN
381	IF (niter .GT. nitermax) THEN
382	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
383	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
384	ENDIF
385	CALL lim_itd_me_ridgeprep
386	ENDIF
387
388	END DO !! on the do while over iter
389
390	!-----------------------------------------------------------------------------!
391	! 4) Ridging diagnostics
392	!-----------------------------------------------------------------------------!
393	! Convert ridging rate diagnostics to correct units.
394	! Update fresh water and heat fluxes due to snow melt.
395
396	dti = 1.0/rdt_ice
397
398	asum_error = .false.
399
400	DO jj = 1, jpj
401	DO ji = 1, jpi
402
403	IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true.
404
405	dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti
406	dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti
407	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti
408	opening (ji,jj) = opening (ji,jj) * dti
409
410	!-----------------------------------------------------------------------------!
411	! 5) Heat, salt and freshwater fluxes
412	!-----------------------------------------------------------------------------!
413	! fresh water source for ocean
414	fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti
415
416	! heat sink for ocean
417	fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti
418
419	END DO
420	END DO
421
422	! Check if there is a ridging error
423	DO jj = 1, jpj
424	DO ji = 1, jpi
425	IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug
426	WRITE(numout,*) ' '
427	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
428	WRITE(numout,*) ' limitd_me '
429	WRITE(numout,*) ' POINT : ', ji, jj
430	WRITE(numout,*) ' jpl, a_i, athorn '
431	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
432	DO jl = 1, jpl
433	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
434	END DO
435	ENDIF ! asum
436
437	END DO !ji
438	END DO !jj
439
440	! Conservation check
441	IF ( con_i ) THEN
442	CALL lim_column_sum (jpl, v_i, vt_i_final)
443	fieldid = ' v_i : limitd_me '
444	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
445	ENDIF
446
447	!-----------------------------------------------------------------------------!
448	! 6) Updating state variables and trend terms
449	!-----------------------------------------------------------------------------!
450
451	CALL lim_var_glo2eqv
452	CALL lim_itd_me_zapsmall
453
454	!-----------------
455	! Trend terms
456	!-----------------
457
458	d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:)
459	d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:)
460	d_a_i_trp(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)
461	d_v_s_trp(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)
462	d_v_i_trp(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)
463	d_e_s_trp(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)
464	d_e_i_trp(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
465	d_oa_i_trp(:,:,:) = oa_i(:,:,:) - old_oa_i(:,:,:)
466	d_smv_i_trp(:,:,:) = 0.0
467	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
468	d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
469
470	IF(ln_ctl) THEN ! Control print
471	CALL prt_ctl_info(' ')
472	CALL prt_ctl_info(' - Cell values : ')
473	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
474	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me : cell area :')
475	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
476	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
477	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
478	DO jl = 1, jpl
479	CALL prt_ctl_info(' ')
480	CALL prt_ctl_info(' - Category : ', ivar1=jl)
481	CALL prt_ctl_info(' ~~~~~~~~~~')
482	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
483	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
484	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
485	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
486	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
487	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
488	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
489	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
490	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
491	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
492	DO jk = 1, nlay_i
493	CALL prt_ctl_info(' ')
494	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
495	CALL prt_ctl_info(' ~~~~~~~')
496	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
497	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
498	END DO
499	END DO
500	ENDIF
501
502	!-------------------------!
503	! Back to initial values
504	!-------------------------!
505
506	! update of fields will be made later in lim update
507	u_ice(:,:) = old_u_ice(:,:)
508	v_ice(:,:) = old_v_ice(:,:)
509	a_i(:,:,:) = old_a_i(:,:,:)
510	v_s(:,:,:) = old_v_s(:,:,:)
511	v_i(:,:,:) = old_v_i(:,:,:)
512	e_s(:,:,:,:) = old_e_s(:,:,:,:)
513	e_i(:,:,:,:) = old_e_i(:,:,:,:)
514	oa_i(:,:,:) = old_oa_i(:,:,:)
515	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
516	smv_i(:,:,:) = old_smv_i(:,:,:)
517
518	!----------------------------------------------------!
519	! Advection of ice in a free cell, newly ridged ice
520	!----------------------------------------------------!
521
522	! to allow for thermodynamics to melt new ice
523	! we immediately advect ice in free cells
524
525	! heat content has to be corrected before ice volume
526	DO jl = 1, jpl
527	DO jk = 1, nlay_i
528	DO jj = 1, jpj
529	DO ji = 1, jpi
530	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
531	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
532	old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl)
533	d_e_i_trp(ji,jj,jk,jl) = 0.0
534	ENDIF
535	END DO
536	END DO
537	END DO
538	END DO
539
540	DO jl = 1, jpl
541	DO jj = 1, jpj
542	DO ji = 1, jpi
543	IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. &
544	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
545	old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl)
546	d_v_i_trp(ji,jj,jl) = 0.0
547	old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl)
548	d_a_i_trp(ji,jj,jl) = 0.0
549	old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl)
550	d_v_s_trp(ji,jj,jl) = 0.0
551	old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)
552	d_e_s_trp(ji,jj,1,jl) = 0.0
553	old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl)
554	d_oa_i_trp(ji,jj,jl) = 0.0
555	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
556	old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)
557	d_smv_i_trp(ji,jj,jl) = 0.0
558	ENDIF
559	END DO
560	END DO
561	END DO
562
563	END SUBROUTINE lim_itd_me
564
565	!===============================================================================
566
567	SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6)
568
569	!!----------------------------------------------------------------------
570	!! * ROUTINE lim_itd_me_icestrength *
571	!! ** Purpose :
572	!! This routine computes ice strength used in dynamics routines
573	!! of ice thickness
574	!!
575	!! ** Method :
576	!! Compute the strength of the ice pack, defined as the energy (J m-2)
577	!! dissipated per unit area removed from the ice pack under compression,
578	!! and assumed proportional to the change in potential energy caused
579	!! by ridging. Note that only Hibler's formulation is stable and that
580	!! ice strength has to be smoothed
581	!!
582	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
583	!!
584	!! ** External :
585	!!
586	!! ** References :
587	!!
588	!!----------------------------------------------------------------------
589	!! * Arguments
590
591	INTEGER, INTENT(in) :: &
592	kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
593
594	INTEGER :: &
595	ji,jj, & !: horizontal indices
596	jl, & !: thickness category index
597	ksmooth, & !: smoothing the resistance to deformation
598	numts_rm !: number of time steps for the P smoothing
599
600	REAL(wp) :: &
601	hi, & !: ice thickness (m)
602	zw1, & !: temporary variable
603	zp, & !: temporary ice strength
604	zdummy
605
606	REAL(wp), DIMENSION(jpi,jpj) :: &
607	zworka !: temporary array used here
608
609	!------------------------------------------------------------------------------!
610	! 1) Initialize
611	!------------------------------------------------------------------------------!
612	strength(:,:) = 0.0
613
614	!------------------------------------------------------------------------------!
615	! 2) Compute thickness distribution of ridging and ridged ice
616	!------------------------------------------------------------------------------!
617	CALL lim_itd_me_ridgeprep
618
619	!------------------------------------------------------------------------------!
620	! 3) Rothrock(1975)'s method
621	!------------------------------------------------------------------------------!
622	IF (kstrngth == 1) then
623
624	DO jl = 1, jpl
625	DO jj= 1, jpj
626	DO ji = 1, jpi
627
628	IF( ( 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 .AND. ln_nicep ) 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	IF( ln_nicep ) 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_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: