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

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4072

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

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

Download in other formats: