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

Last change on this file since 833 was 825, checked in by ctlod, 16 years ago
dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 90.4 KB

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

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

Download in other formats: