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/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4298

Last change on this file since 4298 was 4298, checked in by clem, 10 years ago

Property svn:keywords set to `Id`
File size: 76.3 KB

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

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

Download in other formats: