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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 2760

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

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

Download in other formats: