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

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

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

Context Navigation

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

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

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

Download in other formats: