New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limitd_me.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

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

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

Download in other formats: