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

Context Navigation

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 3977

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