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

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 3625

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

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

Download in other formats: