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

Context Navigation

source: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7278

Last change on this file since 7278 was 7278, checked in by flavoni, 8 years ago
update branch CNRS-2016 to trunk 6720
Property svn:keywords set to `Id`
File size: 48.8 KB

Line
1	MODULE limitd_me
2	!!======================================================================
3	!! * MODULE limitd_me *
4	!! LIM-3 : Mechanical impact on ice thickness distribution
5	!!======================================================================
6	!! History : LIM ! 2006-02 (M. Vancoppenolle) Original code
7	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
9	!!----------------------------------------------------------------------
10	#if defined key_lim3
11	!!----------------------------------------------------------------------
12	!! 'key_lim3' LIM-3 sea-ice model
13	!!----------------------------------------------------------------------
14	USE par_oce ! ocean parameters
15	USE dom_oce ! ocean domain
16	USE phycst ! physical constants (ocean directory)
17	USE sbc_oce ! surface boundary condition: ocean fields
18	USE thd_ice ! LIM thermodynamics
19	USE ice ! LIM variables
20	USE dom_ice ! LIM domain
21	USE limvar ! LIM
22	USE lbclnk ! lateral boundary condition - MPP exchanges
23	USE lib_mpp ! MPP library
24	USE wrk_nemo ! work arrays
25	USE prtctl ! Print control
26
27	USE in_out_manager ! I/O manager
28	USE iom ! I/O manager
29	USE lib_fortran ! glob_sum
30	USE limdiahsb
31	USE timing ! Timing
32	USE limcons ! conservation tests
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC lim_itd_me ! called by ice_stp
38	PUBLIC lim_itd_me_icestrength
39	PUBLIC lim_itd_me_init
40	PUBLIC lim_itd_me_alloc ! called by sbc_lim_init
41
42	!-----------------------------------------------------------------------
43	! Variables shared among ridging subroutines
44	!-----------------------------------------------------------------------
45	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
46	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
47	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n
48	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
49	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
50	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
51	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness
52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
54
55	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
56	REAL(wp), PARAMETER :: kraft = 0.5_wp ! rafting multipliyer
57
58	REAL(wp) :: Cp !
59	!
60	!
61	!!----------------------------------------------------------------------
62	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
63	!! $Id$
64	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
65	!!----------------------------------------------------------------------
66	CONTAINS
67
68	INTEGER FUNCTION lim_itd_me_alloc()
69	!!---------------------------------------------------------------------!
70	!! * ROUTINE lim_itd_me_alloc *
71	!!---------------------------------------------------------------------!
72	ALLOCATE( &
73	!* Variables shared among ridging subroutines
74	& asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , &
75	& aksum(jpi,jpj) , &
76	& hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
77	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc )
78	!
79	IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
80	!
81	END FUNCTION lim_itd_me_alloc
82
83
84	SUBROUTINE lim_itd_me
85	!!---------------------------------------------------------------------!
86	!! * ROUTINE lim_itd_me *
87	!!
88	!! ** Purpose : computes the mechanical redistribution of ice thickness
89	!!
90	!! ** Method : Steps :
91	!! 1) Thickness categories boundaries, ice / o.w. concentrations
92	!! Ridge preparation
93	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
94	!! 3) Ridging iteration
95	!! 4) Ridging diagnostics
96	!! 5) Heat, salt and freshwater fluxes
97	!! 6) Compute increments of tate variables and come back to old values
98	!!
99	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
100	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
101	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
102	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
103	!! Bitz et al., JGR, 2001
104	!! Amundrud and Melling, JGR 2005
105	!! Babko et al., JGR 2002
106	!!
107	!! This routine is based on CICE code and authors William H. Lipscomb,
108	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
109	!!--------------------------------------------------------------------!
110	INTEGER :: ji, jj, jk, jl ! dummy loop index
111	INTEGER :: niter ! local integer
112	INTEGER :: iterate_ridging ! if true, repeat the ridging
113	REAL(wp) :: za, zfac ! local scalar
114	CHARACTER (len = 15) :: fieldid
115	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
116	! (ridging ice area - area of new ridges) / dt
117	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
118	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
119	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
120	!
121	INTEGER, PARAMETER :: nitermax = 20
122	!
123	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
124	!!-----------------------------------------------------------------------------
125	IF( nn_timing == 1 ) CALL timing_start('limitd_me')
126
127	CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross )
128
129	IF(ln_ctl) THEN
130	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
131	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
132	ENDIF
133
134	IF( ln_limdyn ) THEN ! Start ridging and rafting !
135
136	! conservation test
137	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
138
139	!-----------------------------------------------------------------------------!
140	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
141	!-----------------------------------------------------------------------------!
142	Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
143	!
144	CALL lim_itd_me_ridgeprep ! prepare ridging
145	!
146
147	DO jj = 1, jpj ! Initialize arrays.
148	DO ji = 1, jpi
149
150	!-----------------------------------------------------------------------------!
151	! 2) Dynamical inputs (closing rate, divu_adv, opning)
152	!-----------------------------------------------------------------------------!
153	!
154	! 2.1 closing_net
155	!-----------------
156	! Compute the net rate of closing due to convergence
157	! and shear, based on Flato and Hibler (1995).
158	!
159	! The energy dissipation rate is equal to the net closing rate
160	! times the ice strength.
161	!
162	! NOTE: The NET closing rate is equal to the rate that open water
163	! area is removed, plus the rate at which ice area is removed by
164	! ridging, minus the rate at which area is added in new ridges.
165	! The GROSS closing rate is equal to the first two terms (open
166	! water closing and thin ice ridging) without the third term
167	! (thick, newly ridged ice).
168
169	closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
170
171	! 2.2 divu_adv
172	!--------------
173	! Compute divu_adv, the divergence rate given by the transport/
174	! advection scheme, which may not be equal to divu as computed
175	! from the velocity field.
176	!
177	! If divu_adv < 0, make sure the closing rate is large enough
178	! to give asum = 1.0 after ridging.
179
180	divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
181
182	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
183
184	! 2.3 opning
185	!------------
186	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
187	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
188	END DO
189	END DO
190
191	!-----------------------------------------------------------------------------!
192	! 3) Ridging iteration
193	!-----------------------------------------------------------------------------!
194	niter = 1 ! iteration counter
195	iterate_ridging = 1
196
197	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
198
199	! 3.2 closing_gross
200	!-----------------------------------------------------------------------------!
201	! Based on the ITD of ridging and ridged ice, convert the net
202	! closing rate to a gross closing rate.
203	! NOTE: 0 < aksum <= 1
204	closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
205
206	! correction to closing rate and opening if closing rate is excessive
207	!---------------------------------------------------------------------
208	! Reduce the closing rate if more than 100% of the open water
209	! would be removed. Reduce the opening rate proportionately.
210	DO jj = 1, jpj
211	DO ji = 1, jpi
212	za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
213	IF( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i
214	zfac = - ato_i(ji,jj) / za
215	opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice
216	ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum
217	zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za
218	opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice
219	ENDIF
220	END DO
221	END DO
222
223	! correction to closing rate / opening if excessive ice removal
224	!---------------------------------------------------------------
225	! Reduce the closing rate if more than 100% of any ice category
226	! would be removed. Reduce the opening rate proportionately.
227	DO jl = 1, jpl
228	DO jj = 1, jpj
229	DO ji = 1, jpi
230	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
231	IF( za > a_i(ji,jj,jl) ) THEN
232	zfac = a_i(ji,jj,jl) / za
233	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
234	ENDIF
235	END DO
236	END DO
237	END DO
238
239	! 3.3 Redistribute area, volume, and energy.
240	!-----------------------------------------------------------------------------!
241
242	CALL lim_itd_me_ridgeshift( opning, closing_gross )
243
244
245	! 3.4 Compute total area of ice plus open water after ridging.
246	!-----------------------------------------------------------------------------!
247	! This is in general not equal to one because of divergence during transport
248	asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
249
250	! 3.5 Do we keep on iterating ???
251	!-----------------------------------------------------------------------------!
252	! Check whether asum = 1. If not (because the closing and opening
253	! rates were reduced above), ridge again with new rates.
254
255	iterate_ridging = 0
256	DO jj = 1, jpj
257	DO ji = 1, jpi
258	IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
259	closing_net(ji,jj) = 0._wp
260	opning (ji,jj) = 0._wp
261	ELSE
262	iterate_ridging = 1
263	divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
264	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
265	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
266	ENDIF
267	END DO
268	END DO
269
270	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
271
272	! Repeat if necessary.
273	! NOTE: If strength smoothing is turned on, the ridging must be
274	! iterated globally because of the boundary update in the
275	! smoothing.
276
277	niter = niter + 1
278
279	IF( iterate_ridging == 1 ) THEN
280	CALL lim_itd_me_ridgeprep
281	IF( niter > nitermax ) THEN
282	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
283	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
284	ENDIF
285	ENDIF
286
287	END DO !! on the do while over iter
288
289	CALL lim_var_agg( 1 )
290
291	!-----------------------------------------------------------------------------!
292	! control prints
293	!-----------------------------------------------------------------------------!
294	IF(ln_ctl) THEN
295	CALL lim_var_glo2eqv
296
297	CALL prt_ctl_info(' ')
298	CALL prt_ctl_info(' - Cell values : ')
299	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
300	CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' lim_itd_me : cell area :')
301	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
302	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
303	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
304	DO jl = 1, jpl
305	CALL prt_ctl_info(' ')
306	CALL prt_ctl_info(' - Category : ', ivar1=jl)
307	CALL prt_ctl_info(' ~~~~~~~~~~')
308	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
309	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
310	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
311	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
312	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
313	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
314	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
315	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
316	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
317	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
318	DO jk = 1, nlay_i
319	CALL prt_ctl_info(' ')
320	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
321	CALL prt_ctl_info(' ~~~~~~~')
322	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
323	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
324	END DO
325	END DO
326	ENDIF
327
328	! conservation test
329	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
330
331	ENDIF ! ln_limdyn=.true.
332	!
333	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross )
334	!
335	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
336	END SUBROUTINE lim_itd_me
337
338	SUBROUTINE lim_itd_me_ridgeprep
339	!!---------------------------------------------------------------------!
340	!! * ROUTINE lim_itd_me_ridgeprep *
341	!!
342	!! ** Purpose : preparation for ridging and strength calculations
343	!!
344	!! ** Method : Compute the thickness distribution of the ice and open water
345	!! participating in ridging and of the resulting ridges.
346	!!---------------------------------------------------------------------!
347	INTEGER :: ji,jj, jl ! dummy loop indices
348	REAL(wp) :: Gstari, astari, hrmean, zdummy ! local scalar
349	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
350	!------------------------------------------------------------------------------!
351
352	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
353
354	Gstari = 1.0/rn_gstar
355	astari = 1.0/rn_astar
356	aksum(:,:) = 0.0
357	athorn(:,:,:) = 0.0
358	aridge(:,:,:) = 0.0
359	araft (:,:,:) = 0.0
360
361	! Zero out categories with very small areas
362	CALL lim_var_zapsmall
363
364	! Ice thickness needed for rafting
365	DO jl = 1, jpl
366	DO jj = 1, jpj
367	DO ji = 1, jpi
368	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
369	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
370	END DO
371	END DO
372	END DO
373
374	!------------------------------------------------------------------------------!
375	! 1) Participation function
376	!------------------------------------------------------------------------------!
377
378	! Compute total area of ice plus open water.
379	! This is in general not equal to one because of divergence during transport
380	asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
381
382	! Compute cumulative thickness distribution function
383	! Compute the cumulative thickness distribution function Gsum,
384	! where Gsum(n) is the fractional area in categories 0 to n.
385	! initial value (in h = 0) equals open water area
386	Gsum(:,:,-1) = 0._wp
387	Gsum(:,:,0 ) = ato_i(:,:)
388	! for each value of h, you have to add ice concentration then
389	DO jl = 1, jpl
390	Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
391	END DO
392
393	! Normalize the cumulative distribution to 1
394	DO jl = 0, jpl
395	Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
396	END DO
397
398	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
399	!--------------------------------------------------------------------------------------------------
400	! Compute the participation function athorn; this is analogous to
401	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
402	! area lost from category n due to ridging/closing
403	! athorn(n) = total area lost due to ridging/closing
404	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
405	!
406	! The expressions for athorn are found by integrating b(h)g(h) between
407	! the category boundaries.
408	! athorn is always >= 0 and SUM(athorn(0:jpl))=1
409	!-----------------------------------------------------------------
410
411	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
412	DO jl = 0, jpl
413	DO jj = 1, jpj
414	DO ji = 1, jpi
415	IF ( Gsum(ji,jj,jl) < rn_gstar ) THEN
416	athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
417	& ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
418	ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
419	athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &
420	& ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
421	ELSE
422	athorn(ji,jj,jl) = 0._wp
423	ENDIF
424	END DO
425	END DO
426	END DO
427
428	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
429	!
430	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
431	DO jl = -1, jpl
432	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
433	END DO
434	DO jl = 0, jpl
435	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
436	END DO
437	!
438	ENDIF
439
440	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
441	!
442	DO jl = 1, jpl
443	DO jj = 1, jpj
444	DO ji = 1, jpi
445	zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
446	aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
447	araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl)
448	END DO
449	END DO
450	END DO
451
452	ELSE
453	!
454	DO jl = 1, jpl
455	aridge(:,:,jl) = athorn(:,:,jl)
456	END DO
457	!
458	ENDIF
459
460	!-----------------------------------------------------------------
461	! 2) Transfer function
462	!-----------------------------------------------------------------
463	! Compute max and min ridged ice thickness for each ridging category.
464	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
465	!
466	! This parameterization is a modified version of Hibler (1980).
467	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
468	! and for very thick ridging ice must be >= krdgmin*hi
469	!
470	! The minimum ridging thickness, hrmin, is equal to 2*hi
471	! (i.e., rafting) and for very thick ridging ice is
472	! constrained by hrmin <= (hrmean + hi)/2.
473	!
474	! The maximum ridging thickness, hrmax, is determined by
475	! hrmean and hrmin.
476	!
477	! These modifications have the effect of reducing the ice strength
478	! (relative to the Hibler formulation) when very thick ice is
479	! ridging.
480	!
481	! aksum = net area removed/ total area removed
482	! where total area removed = area of ice that ridges
483	! net area removed = total area removed - area of new ridges
484	!-----------------------------------------------------------------
485
486	aksum(:,:) = athorn(:,:,0)
487	! Transfer function
488	DO jl = 1, jpl !all categories have a specific transfer function
489	DO jj = 1, jpj
490	DO ji = 1, jpi
491
492	IF( athorn(ji,jj,jl) > 0._wp ) THEN
493	hrmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
494	hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
495	hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
496	hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
497	krdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
498
499	! Normalization factor : aksum, ensures mass conservation
500	aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) ) &
501	& + araft (ji,jj,jl) * ( 1._wp - kraft )
502
503	ELSE
504	hrmin(ji,jj,jl) = 0._wp
505	hrmax(ji,jj,jl) = 0._wp
506	hraft(ji,jj,jl) = 0._wp
507	krdg (ji,jj,jl) = 1._wp
508	ENDIF
509
510	END DO
511	END DO
512	END DO
513	!
514	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
515	!
516	END SUBROUTINE lim_itd_me_ridgeprep
517
518
519	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross )
520	!!----------------------------------------------------------------------
521	!! * ROUTINE lim_itd_me_icestrength *
522	!!
523	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
524	!!
525	!! ** Method : Remove area, volume, and energy from each ridging category
526	!! and add to thicker ice categories.
527	!!----------------------------------------------------------------------
528	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
529	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
530	!
531	CHARACTER (len=80) :: fieldid ! field identifier
532	!
533	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
534	INTEGER :: ij ! horizontal index, combines i and j loops
535	INTEGER :: icells ! number of cells with a_i > puny
536	REAL(wp) :: hL, hR, farea ! left and right limits of integration
537
538	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
539	REAL(wp), POINTER, DIMENSION(:) :: zswitch, fvol ! new ridge volume going to n2
540
541	REAL(wp), POINTER, DIMENSION(:) :: afrac ! fraction of category area ridged
542	REAL(wp), POINTER, DIMENSION(:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
543	REAL(wp), POINTER, DIMENSION(:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
544	REAL(wp), POINTER, DIMENSION(:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
545
546	REAL(wp), POINTER, DIMENSION(:) :: vrdg1 ! volume of ice ridged
547	REAL(wp), POINTER, DIMENSION(:) :: vrdg2 ! volume of new ridges
548	REAL(wp), POINTER, DIMENSION(:) :: vsw ! volume of seawater trapped into ridges
549	REAL(wp), POINTER, DIMENSION(:) :: srdg1 ! sal*volume of ice ridged
550	REAL(wp), POINTER, DIMENSION(:) :: srdg2 ! sal*volume of new ridges
551	REAL(wp), POINTER, DIMENSION(:) :: smsw ! sal*volume of water trapped into ridges
552	REAL(wp), POINTER, DIMENSION(:) :: oirdg1, oirdg2 ! ice age of ice ridged
553
554	REAL(wp), POINTER, DIMENSION(:) :: afrft ! fraction of category area rafted
555	REAL(wp), POINTER, DIMENSION(:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
556	REAL(wp), POINTER, DIMENSION(:) :: virft , vsrft ! ice & snow volume of rafting ice
557	REAL(wp), POINTER, DIMENSION(:) :: esrft , smrft ! snow energy & salinity of rafting ice
558	REAL(wp), POINTER, DIMENSION(:) :: oirft1, oirft2 ! ice age of ice rafted
559
560	REAL(wp), POINTER, DIMENSION(:,:) :: eirft ! ice energy of rafting ice
561	REAL(wp), POINTER, DIMENSION(:,:) :: erdg1 ! enth*volume of ice ridged
562	REAL(wp), POINTER, DIMENSION(:,:) :: erdg2 ! enth*volume of new ridges
563	REAL(wp), POINTER, DIMENSION(:,:) :: ersw ! enth of water trapped into ridges
564	!!----------------------------------------------------------------------
565
566	CALL wrk_alloc( jpij, indxi, indxj )
567	CALL wrk_alloc( jpij, zswitch, fvol )
568	CALL wrk_alloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
569	CALL wrk_alloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
570	CALL wrk_alloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
571	CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
572
573	!-------------------------------------------------------------------------------
574	! 1) Compute change in open water area due to closing and opening.
575	!-------------------------------------------------------------------------------
576	DO jj = 1, jpj
577	DO ji = 1, jpi
578	ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + &
579	& ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
580	END DO
581	END DO
582
583	!-----------------------------------------------------------------
584	! 3) Pump everything from ice which is being ridged / rafted
585	!-----------------------------------------------------------------
586	! Compute the area, volume, and energy of ice ridging in each
587	! category, along with the area of the resulting ridge.
588
589	DO jl1 = 1, jpl !jl1 describes the ridging category
590
591	!------------------------------------------------
592	! 3.1) Identify grid cells with nonzero ridging
593	!------------------------------------------------
594	icells = 0
595	DO jj = 1, jpj
596	DO ji = 1, jpi
597	IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
598	icells = icells + 1
599	indxi(icells) = ji
600	indxj(icells) = jj
601	ENDIF
602	END DO
603	END DO
604
605	DO ij = 1, icells
606	ji = indxi(ij) ; jj = indxj(ij)
607
608	!--------------------------------------------------------------------
609	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
610	!--------------------------------------------------------------------
611	ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
612	arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
613
614	!---------------------------------------------------------------
615	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
616	!---------------------------------------------------------------
617	afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
618	afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
619	ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
620	arft2(ij) = arft1(ij) * kraft
621
622	!--------------------------------------------------------------------------
623	! 3.4) Subtract area, volume, and energy from ridging
624	! / rafting category n1.
625	!--------------------------------------------------------------------------
626	vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
627	vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
628	vsw (ij) = vrdg1(ij) * rn_por_rdg
629
630	vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij)
631	esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij)
632	srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij)
633	oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij)
634	oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
635
636	! rafting volumes, heat contents ...
637	virft (ij) = v_i (ji,jj, jl1) * afrft(ij)
638	vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij)
639	esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij)
640	smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij)
641	oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij)
642	oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * kraft
643
644	!-----------------------------------------------------------------
645	! 3.5) Compute properties of new ridges
646	!-----------------------------------------------------------------
647	smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids
648	srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge
649
650	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
651	wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids
652
653	! virtual salt flux to keep salinity constant
654	IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
655	srdg2(ij) = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) ) ! ridge salinity = sm_i
656	sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj) * vsw(ij) * rhoic * r1_rdtice & ! put back sss_m into the ocean
657	& - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean
658	ENDIF
659
660	!------------------------------------------
661	! 3.7 Put the snow somewhere in the ocean
662	!------------------------------------------
663	! Place part of the snow lost by ridging into the ocean.
664	! Note that esrdg > 0; the ocean must cool to melt snow.
665	! If the ocean temp = Tf already, new ice must grow.
666	! During the next time step, thermo_rates will determine whether
667	! the ocean cools or new ice grows.
668	wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg ) &
669	& + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! fresh water source for ocean
670
671	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg ) &
672	& - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
673
674	!-----------------------------------------------------------------
675	! 3.8 Compute quantities used to apportion ice among categories
676	! in the n2 loop below
677	!-----------------------------------------------------------------
678	dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) )
679	dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
680
681
682	! update jl1 (removing ridged/rafted area)
683	a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij)
684	v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij)
685	v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij)
686	e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
687	smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij)
688	oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij)
689
690	END DO
691
692	!--------------------------------------------------------------------
693	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
694	! compute ridged ice enthalpy
695	!--------------------------------------------------------------------
696	DO jk = 1, nlay_i
697	DO ij = 1, icells
698	ji = indxi(ij) ; jj = indxj(ij)
699	! heat content of ridged ice
700	erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij)
701	eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)
702
703	! enthalpy of the trapped seawater (J/m2, >0)
704	! clem: if sst>0, then ersw <0 (is that possible?)
705	ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
706
707	! heat flux to the ocean
708	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
709
710	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
711	erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
712
713	! update jl1
714	e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
715
716	END DO
717	END DO
718
719	!-------------------------------------------------------------------------------
720	! 4) Add area, volume, and energy of new ridge to each category jl2
721	!-------------------------------------------------------------------------------
722	DO jl2 = 1, jpl
723	! over categories to which ridged/rafted ice is transferred
724	DO ij = 1, icells
725	ji = indxi(ij) ; jj = indxj(ij)
726
727	! Compute the fraction of ridged ice area and volume going to thickness category jl2.
728	IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
729	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
730	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
731	farea = ( hR - hL ) * dhr(ij)
732	fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
733	ELSE
734	farea = 0._wp
735	fvol(ij) = 0._wp
736	ENDIF
737
738	! Compute the fraction of rafted ice area and volume going to thickness category jl2
739	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
740	zswitch(ij) = 1._wp
741	ELSE
742	zswitch(ij) = 0._wp
743	ENDIF
744
745	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) )
746	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) )
747	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
748	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
749	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij) + &
750	& vsrft (ij) * rn_fsnowrft * zswitch(ij) )
751	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij) + &
752	& esrft (ij) * rn_fsnowrft * zswitch(ij) )
753
754	END DO
755
756	! Transfer ice energy to category jl2 by ridging
757	DO jk = 1, nlay_i
758	DO ij = 1, icells
759	ji = indxi(ij) ; jj = indxj(ij)
760	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)
761	END DO
762	END DO
763	!
764	END DO ! jl2
765
766	END DO ! jl1 (deforming categories)
767
768	!
769	CALL wrk_dealloc( jpij, indxi, indxj )
770	CALL wrk_dealloc( jpij, zswitch, fvol )
771	CALL wrk_dealloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
772	CALL wrk_dealloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
773	CALL wrk_dealloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
774	CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
775	!
776	END SUBROUTINE lim_itd_me_ridgeshift
777
778	SUBROUTINE lim_itd_me_icestrength( kstrngth )
779	!!----------------------------------------------------------------------
780	!! * ROUTINE lim_itd_me_icestrength *
781	!!
782	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
783	!!
784	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
785	!! dissipated per unit area removed from the ice pack under compression,
786	!! and assumed proportional to the change in potential energy caused
787	!! by ridging. Note that only Hibler's formulation is stable and that
788	!! ice strength has to be smoothed
789	!!
790	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
791	!!----------------------------------------------------------------------
792	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
793	INTEGER :: ji,jj, jl ! dummy loop indices
794	INTEGER :: ksmooth ! smoothing the resistance to deformation
795	INTEGER :: numts_rm ! number of time steps for the P smoothing
796	REAL(wp) :: zp, z1_3 ! local scalars
797	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
798	!!----------------------------------------------------------------------
799
800	CALL wrk_alloc( jpi, jpj, zworka )
801
802	!------------------------------------------------------------------------------!
803	! 1) Initialize
804	!------------------------------------------------------------------------------!
805	strength(:,:) = 0._wp
806
807	!------------------------------------------------------------------------------!
808	! 2) Compute thickness distribution of ridging and ridged ice
809	!------------------------------------------------------------------------------!
810	CALL lim_itd_me_ridgeprep
811
812	!------------------------------------------------------------------------------!
813	! 3) Rothrock(1975)'s method
814	!------------------------------------------------------------------------------!
815	IF( kstrngth == 1 ) THEN
816	z1_3 = 1._wp / 3._wp
817	DO jl = 1, jpl
818	DO jj= 1, jpj
819	DO ji = 1, jpi
820	!
821	IF( athorn(ji,jj,jl) > 0._wp ) THEN
822	!----------------------------
823	! PE loss from deforming ice
824	!----------------------------
825	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
826
827	!--------------------------
828	! PE gain from rafting ice
829	!--------------------------
830	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
831
832	!----------------------------
833	! PE gain from ridging ice
834	!----------------------------
835	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 * &
836	& ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) + &
837	& hrmin(ji,jj,jl) * hrmin(ji,jj,jl) + &
838	& hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
839	!!(a3-b3)/(a-b) = aa+ab+bb
840	ENDIF
841	!
842	END DO
843	END DO
844	END DO
845
846	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
847	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
848	ksmooth = 1
849
850	!------------------------------------------------------------------------------!
851	! 4) Hibler (1979)' method
852	!------------------------------------------------------------------------------!
853	ELSE ! kstrngth ne 1: Hibler (1979) form
854	!
855	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
856	!
857	ksmooth = 1
858	!
859	ENDIF ! kstrngth
860	!
861	!------------------------------------------------------------------------------!
862	! 5) Impact of brine volume
863	!------------------------------------------------------------------------------!
864	! CAN BE REMOVED
865	IF( ln_icestr_bvf ) THEN
866	DO jj = 1, jpj
867	DO ji = 1, jpi
868	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
869	END DO
870	END DO
871	ENDIF
872	!
873	!------------------------------------------------------------------------------!
874	! 6) Smoothing ice strength
875	!------------------------------------------------------------------------------!
876	!
877	!-------------------
878	! Spatial smoothing
879	!-------------------
880	IF ( ksmooth == 1 ) THEN
881
882	CALL lbc_lnk( strength, 'T', 1. )
883
884	DO jj = 2, jpjm1
885	DO ji = 2, jpim1
886	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
887	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
888	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
889	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
890	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
891	ELSE
892	zworka(ji,jj) = 0._wp
893	ENDIF
894	END DO
895	END DO
896
897	DO jj = 2, jpjm1
898	DO ji = 2, jpim1
899	strength(ji,jj) = zworka(ji,jj)
900	END DO
901	END DO
902	CALL lbc_lnk( strength, 'T', 1. )
903
904	ENDIF
905
906	!--------------------
907	! Temporal smoothing
908	!--------------------
909	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
910	strp1(:,:) = 0.0
911	strp2(:,:) = 0.0
912	ENDIF
913
914	IF ( ksmooth == 2 ) THEN
915
916	CALL lbc_lnk( strength, 'T', 1. )
917
918	DO jj = 1, jpj - 1
919	DO ji = 1, jpi - 1
920	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
921	numts_rm = 1 ! number of time steps for the running mean
922	IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
923	IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
924	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
925	strp2(ji,jj) = strp1(ji,jj)
926	strp1(ji,jj) = strength(ji,jj)
927	strength(ji,jj) = zp
928
929	ENDIF
930	END DO
931	END DO
932
933	ENDIF ! ksmooth
934
935	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
936
937	CALL wrk_dealloc( jpi, jpj, zworka )
938	!
939	END SUBROUTINE lim_itd_me_icestrength
940
941	SUBROUTINE lim_itd_me_init
942	!!-------------------------------------------------------------------
943	!! * ROUTINE lim_itd_me_init *
944	!!
945	!! ** Purpose : Physical constants and parameters linked
946	!! to the mechanical ice redistribution
947	!!
948	!! ** Method : Read the namiceitdme namelist
949	!! and check the parameters values
950	!! called at the first timestep (nit000)
951	!!
952	!! ** input : Namelist namiceitdme
953	!!-------------------------------------------------------------------
954	INTEGER :: ios ! Local integer output status for namelist read
955	NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, &
956	& rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
957	& nn_partfun
958	!!-------------------------------------------------------------------
959	!
960	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
961	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
962	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
963
964	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
965	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
966	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
967	IF(lwm) WRITE ( numoni, namiceitdme )
968	!
969	IF (lwp) THEN ! control print
970	WRITE(numout,*)
971	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
972	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
973	WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs
974	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg
975	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft
976	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar
977	WRITE(numout,)' Equivalent to G for an exponential part function rn_astar = ', rn_astar
978	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar
979	WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting
980	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
981	WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
982	WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg
983	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun
984	ENDIF
985	!
986	END SUBROUTINE lim_itd_me_init
987
988	#else
989	!!----------------------------------------------------------------------
990	!! Default option Empty module NO LIM-3 sea-ice model
991	!!----------------------------------------------------------------------
992	CONTAINS
993	SUBROUTINE lim_itd_me ! Empty routines
994	END SUBROUTINE lim_itd_me
995	SUBROUTINE lim_itd_me_icestrength
996	END SUBROUTINE lim_itd_me_icestrength
997	SUBROUTINE lim_itd_me_init
998	END SUBROUTINE lim_itd_me_init
999	#endif
1000	!!======================================================================
1001	END MODULE limitd_me

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

Download in other formats: