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

Context Navigation

source: branches/CMIP5_IPSL/NEMO/LIM_SRC_3/limitd_me.F90 @ 8506

Last change on this file since 8506 was 1572, checked in by ctlod, 15 years ago
Cosmetic changes only following the bugs correction related to ticket: #505
Property svn:keywords set to `Id`
File size: 74.7 KB

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

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

Download in other formats: