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

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

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

Download in other formats: