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

Context Navigation

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago
[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]
Property svn:keywords set to `Id`
File size: 72.8 KB

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

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

Download in other formats: