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 NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 10379

Last change on this file since 10379 was 7814, checked in by clem, 7 years ago
cosmetic changes for better phasing with trunk
Property svn:keywords set to `Id`
File size: 49.0 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	ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
262	ELSE
263	iterate_ridging = 1
264	divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
265	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
266	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
267	ENDIF
268	END DO
269	END DO
270
271	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
272
273	! Repeat if necessary.
274	! NOTE: If strength smoothing is turned on, the ridging must be
275	! iterated globally because of the boundary update in the
276	! smoothing.
277
278	niter = niter + 1
279
280	IF( iterate_ridging == 1 ) THEN
281	CALL lim_itd_me_ridgeprep
282	IF( niter > nitermax ) THEN
283	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
284	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
285	ENDIF
286	ENDIF
287
288	END DO !! on the do while over iter
289
290	CALL lim_var_agg( 1 )
291
292	!-----------------------------------------------------------------------------!
293	! control prints
294	!-----------------------------------------------------------------------------!
295	IF(ln_ctl) THEN
296	CALL lim_var_glo2eqv
297
298	CALL prt_ctl_info(' ')
299	CALL prt_ctl_info(' - Cell values : ')
300	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
301	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :')
302	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
303	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
304	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
305	DO jl = 1, jpl
306	CALL prt_ctl_info(' ')
307	CALL prt_ctl_info(' - Category : ', ivar1=jl)
308	CALL prt_ctl_info(' ~~~~~~~~~~')
309	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
310	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
311	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
312	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
313	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
314	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
315	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
316	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
317	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
318	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
319	DO jk = 1, nlay_i
320	CALL prt_ctl_info(' ')
321	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
322	CALL prt_ctl_info(' ~~~~~~~')
323	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
324	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
325	END DO
326	END DO
327	ENDIF
328
329	! conservation test
330	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
331
332	ENDIF ! ln_limdyn=.true.
333	!
334	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross )
335	!
336	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
337	END SUBROUTINE lim_itd_me
338
339	SUBROUTINE lim_itd_me_ridgeprep
340	!!---------------------------------------------------------------------!
341	!! * ROUTINE lim_itd_me_ridgeprep *
342	!!
343	!! ** Purpose : preparation for ridging and strength calculations
344	!!
345	!! ** Method : Compute the thickness distribution of the ice and open water
346	!! participating in ridging and of the resulting ridges.
347	!!---------------------------------------------------------------------!
348	INTEGER :: ji,jj, jl ! dummy loop indices
349	REAL(wp) :: Gstari, astari, hrmean, zdummy ! local scalar
350	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
351	!------------------------------------------------------------------------------!
352
353	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
354
355	Gstari = 1.0/rn_gstar
356	astari = 1.0/rn_astar
357	aksum(:,:) = 0.0
358	athorn(:,:,:) = 0.0
359	aridge(:,:,:) = 0.0
360	araft (:,:,:) = 0.0
361
362	! Zero out categories with very small areas
363	CALL lim_var_zapsmall
364
365	! Ice thickness needed for rafting
366	DO jl = 1, jpl
367	DO jj = 1, jpj
368	DO ji = 1, jpi
369	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
370	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
371	END DO
372	END DO
373	END DO
374
375	!------------------------------------------------------------------------------!
376	! 1) Participation function
377	!------------------------------------------------------------------------------!
378
379	! Compute total area of ice plus open water.
380	! This is in general not equal to one because of divergence during transport
381	asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
382
383	! Compute cumulative thickness distribution function
384	! Compute the cumulative thickness distribution function Gsum,
385	! where Gsum(n) is the fractional area in categories 0 to n.
386	! initial value (in h = 0) equals open water area
387	Gsum(:,:,-1) = 0._wp
388	Gsum(:,:,0 ) = ato_i(:,:)
389	! for each value of h, you have to add ice concentration then
390	DO jl = 1, jpl
391	Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
392	END DO
393
394	! Normalize the cumulative distribution to 1
395	DO jl = 0, jpl
396	Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
397	END DO
398
399	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
400	!--------------------------------------------------------------------------------------------------
401	! Compute the participation function athorn; this is analogous to
402	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
403	! area lost from category n due to ridging/closing
404	! athorn(n) = total area lost due to ridging/closing
405	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
406	!
407	! The expressions for athorn are found by integrating b(h)g(h) between
408	! the category boundaries.
409	! athorn is always >= 0 and SUM(athorn(0:jpl))=1
410	!-----------------------------------------------------------------
411
412	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
413	DO jl = 0, jpl
414	DO jj = 1, jpj
415	DO ji = 1, jpi
416	IF ( Gsum(ji,jj,jl) < rn_gstar ) THEN
417	athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
418	& ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
419	ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
420	athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &
421	& ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
422	ELSE
423	athorn(ji,jj,jl) = 0._wp
424	ENDIF
425	END DO
426	END DO
427	END DO
428
429	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
430	!
431	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
432	DO jl = -1, jpl
433	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
434	END DO
435	DO jl = 0, jpl
436	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
437	END DO
438	!
439	ENDIF
440
441	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
442	!
443	DO jl = 1, jpl
444	DO jj = 1, jpj
445	DO ji = 1, jpi
446	zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
447	aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
448	araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl)
449	END DO
450	END DO
451	END DO
452
453	ELSE
454	!
455	DO jl = 1, jpl
456	aridge(:,:,jl) = athorn(:,:,jl)
457	END DO
458	!
459	ENDIF
460
461	!-----------------------------------------------------------------
462	! 2) Transfer function
463	!-----------------------------------------------------------------
464	! Compute max and min ridged ice thickness for each ridging category.
465	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
466	!
467	! This parameterization is a modified version of Hibler (1980).
468	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
469	! and for very thick ridging ice must be >= krdgmin*hi
470	!
471	! The minimum ridging thickness, hrmin, is equal to 2*hi
472	! (i.e., rafting) and for very thick ridging ice is
473	! constrained by hrmin <= (hrmean + hi)/2.
474	!
475	! The maximum ridging thickness, hrmax, is determined by
476	! hrmean and hrmin.
477	!
478	! These modifications have the effect of reducing the ice strength
479	! (relative to the Hibler formulation) when very thick ice is
480	! ridging.
481	!
482	! aksum = net area removed/ total area removed
483	! where total area removed = area of ice that ridges
484	! net area removed = total area removed - area of new ridges
485	!-----------------------------------------------------------------
486
487	aksum(:,:) = athorn(:,:,0)
488	! Transfer function
489	DO jl = 1, jpl !all categories have a specific transfer function
490	DO jj = 1, jpj
491	DO ji = 1, jpi
492
493	IF( athorn(ji,jj,jl) > 0._wp ) THEN
494	hrmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
495	hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
496	hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
497	hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
498	krdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
499
500	! Normalization factor : aksum, ensures mass conservation
501	aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) ) &
502	& + araft (ji,jj,jl) * ( 1._wp - kraft )
503
504	ELSE
505	hrmin(ji,jj,jl) = 0._wp
506	hrmax(ji,jj,jl) = 0._wp
507	hraft(ji,jj,jl) = 0._wp
508	krdg (ji,jj,jl) = 1._wp
509	ENDIF
510
511	END DO
512	END DO
513	END DO
514	!
515	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
516	!
517	END SUBROUTINE lim_itd_me_ridgeprep
518
519
520	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross )
521	!!----------------------------------------------------------------------
522	!! * ROUTINE lim_itd_me_icestrength *
523	!!
524	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
525	!!
526	!! ** Method : Remove area, volume, and energy from each ridging category
527	!! and add to thicker ice categories.
528	!!----------------------------------------------------------------------
529	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
530	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
531	!
532	CHARACTER (len=80) :: fieldid ! field identifier
533	!
534	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
535	INTEGER :: ij ! horizontal index, combines i and j loops
536	INTEGER :: icells ! number of cells with a_i > puny
537	REAL(wp) :: hL, hR, farea ! left and right limits of integration
538
539	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
540	REAL(wp), POINTER, DIMENSION(:) :: zswitch, fvol ! new ridge volume going to n2
541
542	REAL(wp), POINTER, DIMENSION(:) :: afrac ! fraction of category area ridged
543	REAL(wp), POINTER, DIMENSION(:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
544	REAL(wp), POINTER, DIMENSION(:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
545	REAL(wp), POINTER, DIMENSION(:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
546
547	REAL(wp), POINTER, DIMENSION(:) :: vrdg1 ! volume of ice ridged
548	REAL(wp), POINTER, DIMENSION(:) :: vrdg2 ! volume of new ridges
549	REAL(wp), POINTER, DIMENSION(:) :: vsw ! volume of seawater trapped into ridges
550	REAL(wp), POINTER, DIMENSION(:) :: srdg1 ! sal*volume of ice ridged
551	REAL(wp), POINTER, DIMENSION(:) :: srdg2 ! sal*volume of new ridges
552	REAL(wp), POINTER, DIMENSION(:) :: smsw ! sal*volume of water trapped into ridges
553	REAL(wp), POINTER, DIMENSION(:) :: oirdg1, oirdg2 ! ice age of ice ridged
554
555	REAL(wp), POINTER, DIMENSION(:) :: afrft ! fraction of category area rafted
556	REAL(wp), POINTER, DIMENSION(:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
557	REAL(wp), POINTER, DIMENSION(:) :: virft , vsrft ! ice & snow volume of rafting ice
558	REAL(wp), POINTER, DIMENSION(:) :: esrft , smrft ! snow energy & salinity of rafting ice
559	REAL(wp), POINTER, DIMENSION(:) :: oirft1, oirft2 ! ice age of ice rafted
560
561	REAL(wp), POINTER, DIMENSION(:,:) :: eirft ! ice energy of rafting ice
562	REAL(wp), POINTER, DIMENSION(:,:) :: erdg1 ! enth*volume of ice ridged
563	REAL(wp), POINTER, DIMENSION(:,:) :: erdg2 ! enth*volume of new ridges
564	REAL(wp), POINTER, DIMENSION(:,:) :: ersw ! enth of water trapped into ridges
565	!!----------------------------------------------------------------------
566
567	CALL wrk_alloc( jpij, indxi, indxj )
568	CALL wrk_alloc( jpij, zswitch, fvol )
569	CALL wrk_alloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
570	CALL wrk_alloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
571	CALL wrk_alloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
572	CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
573
574	!-------------------------------------------------------------------------------
575	! 1) Compute change in open water area due to closing and opening.
576	!-------------------------------------------------------------------------------
577	DO jj = 1, jpj
578	DO ji = 1, jpi
579	ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + &
580	& ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
581	END DO
582	END DO
583
584	!-----------------------------------------------------------------
585	! 3) Pump everything from ice which is being ridged / rafted
586	!-----------------------------------------------------------------
587	! Compute the area, volume, and energy of ice ridging in each
588	! category, along with the area of the resulting ridge.
589
590	DO jl1 = 1, jpl !jl1 describes the ridging category
591
592	!------------------------------------------------
593	! 3.1) Identify grid cells with nonzero ridging
594	!------------------------------------------------
595	icells = 0
596	DO jj = 1, jpj
597	DO ji = 1, jpi
598	IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
599	icells = icells + 1
600	indxi(icells) = ji
601	indxj(icells) = jj
602	ENDIF
603	END DO
604	END DO
605
606	DO ij = 1, icells
607	ji = indxi(ij) ; jj = indxj(ij)
608
609	!--------------------------------------------------------------------
610	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
611	!--------------------------------------------------------------------
612	ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
613	arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
614
615	!---------------------------------------------------------------
616	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
617	!---------------------------------------------------------------
618	afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
619	afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
620	ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
621	arft2(ij) = arft1(ij) * kraft
622
623	!--------------------------------------------------------------------------
624	! 3.4) Subtract area, volume, and energy from ridging
625	! / rafting category n1.
626	!--------------------------------------------------------------------------
627	vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
628	vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
629	vsw (ij) = vrdg1(ij) * rn_por_rdg
630
631	vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij)
632	esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij)
633	srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij)
634	oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij)
635	oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
636
637	! rafting volumes, heat contents ...
638	virft (ij) = v_i (ji,jj, jl1) * afrft(ij)
639	vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij)
640	esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij)
641	smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij)
642	oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij)
643	oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * kraft
644
645	!-----------------------------------------------------------------
646	! 3.5) Compute properties of new ridges
647	!-----------------------------------------------------------------
648	smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids
649	srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge
650
651	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
652	wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids
653
654	! virtual salt flux to keep salinity constant
655	IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
656	srdg2(ij) = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) ) ! ridge salinity = sm_i
657	sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj) * vsw(ij) * rhoic * r1_rdtice & ! put back sss_m into the ocean
658	& - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean
659	ENDIF
660
661	!------------------------------------------
662	! 3.7 Put the snow somewhere in the ocean
663	!------------------------------------------
664	! Place part of the snow lost by ridging into the ocean.
665	! Note that esrdg > 0; the ocean must cool to melt snow.
666	! If the ocean temp = Tf already, new ice must grow.
667	! During the next time step, thermo_rates will determine whether
668	! the ocean cools or new ice grows.
669	wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg ) &
670	& + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! fresh water source for ocean
671
672	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg ) &
673	& - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
674
675	!-----------------------------------------------------------------
676	! 3.8 Compute quantities used to apportion ice among categories
677	! in the n2 loop below
678	!-----------------------------------------------------------------
679	dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) )
680	dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
681
682
683	! update jl1 (removing ridged/rafted area)
684	a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij)
685	v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij)
686	v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij)
687	e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
688	smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij)
689	oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij)
690
691	END DO
692
693	!--------------------------------------------------------------------
694	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
695	! compute ridged ice enthalpy
696	!--------------------------------------------------------------------
697	DO jk = 1, nlay_i
698	DO ij = 1, icells
699	ji = indxi(ij) ; jj = indxj(ij)
700	! heat content of ridged ice
701	erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij)
702	eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)
703
704	! enthalpy of the trapped seawater (J/m2, >0)
705	! clem: if sst>0, then ersw <0 (is that possible?)
706	ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
707
708	! heat flux to the ocean
709	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
710
711	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
712	erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
713
714	! update jl1
715	e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
716
717	END DO
718	END DO
719
720	!-------------------------------------------------------------------------------
721	! 4) Add area, volume, and energy of new ridge to each category jl2
722	!-------------------------------------------------------------------------------
723	DO jl2 = 1, jpl
724	! over categories to which ridged/rafted ice is transferred
725	DO ij = 1, icells
726	ji = indxi(ij) ; jj = indxj(ij)
727
728	! Compute the fraction of ridged ice area and volume going to thickness category jl2.
729	IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
730	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
731	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
732	farea = ( hR - hL ) * dhr(ij)
733	fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
734	ELSE
735	farea = 0._wp
736	fvol(ij) = 0._wp
737	ENDIF
738
739	! Compute the fraction of rafted ice area and volume going to thickness category jl2
740	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
741	zswitch(ij) = 1._wp
742	ELSE
743	zswitch(ij) = 0._wp
744	ENDIF
745
746	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) )
747	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) )
748	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
749	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
750	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij) + &
751	& vsrft (ij) * rn_fsnowrft * zswitch(ij) )
752	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij) + &
753	& esrft (ij) * rn_fsnowrft * zswitch(ij) )
754
755	END DO
756
757	! Transfer ice energy to category jl2 by ridging
758	DO jk = 1, nlay_i
759	DO ij = 1, icells
760	ji = indxi(ij) ; jj = indxj(ij)
761	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)
762	END DO
763	END DO
764	!
765	END DO ! jl2
766
767	END DO ! jl1 (deforming categories)
768
769	!
770	CALL wrk_dealloc( jpij, indxi, indxj )
771	CALL wrk_dealloc( jpij, zswitch, fvol )
772	CALL wrk_dealloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
773	CALL wrk_dealloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
774	CALL wrk_dealloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
775	CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
776	!
777	END SUBROUTINE lim_itd_me_ridgeshift
778
779	SUBROUTINE lim_itd_me_icestrength( kstrngth )
780	!!----------------------------------------------------------------------
781	!! * ROUTINE lim_itd_me_icestrength *
782	!!
783	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
784	!!
785	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
786	!! dissipated per unit area removed from the ice pack under compression,
787	!! and assumed proportional to the change in potential energy caused
788	!! by ridging. Note that only Hibler's formulation is stable and that
789	!! ice strength has to be smoothed
790	!!
791	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
792	!!----------------------------------------------------------------------
793	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
794	INTEGER :: ji,jj, jl ! dummy loop indices
795	INTEGER :: ksmooth ! smoothing the resistance to deformation
796	INTEGER :: numts_rm ! number of time steps for the P smoothing
797	REAL(wp) :: zp, z1_3 ! local scalars
798	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
799	!!----------------------------------------------------------------------
800
801	CALL wrk_alloc( jpi, jpj, zworka )
802
803	!------------------------------------------------------------------------------!
804	! 1) Initialize
805	!------------------------------------------------------------------------------!
806	strength(:,:) = 0._wp
807
808	!------------------------------------------------------------------------------!
809	! 2) Compute thickness distribution of ridging and ridged ice
810	!------------------------------------------------------------------------------!
811	CALL lim_itd_me_ridgeprep
812
813	!------------------------------------------------------------------------------!
814	! 3) Rothrock(1975)'s method
815	!------------------------------------------------------------------------------!
816	IF( kstrngth == 1 ) THEN
817	z1_3 = 1._wp / 3._wp
818	DO jl = 1, jpl
819	DO jj= 1, jpj
820	DO ji = 1, jpi
821	!
822	IF( athorn(ji,jj,jl) > 0._wp ) THEN
823	!----------------------------
824	! PE loss from deforming ice
825	!----------------------------
826	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
827
828	!--------------------------
829	! PE gain from rafting ice
830	!--------------------------
831	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
832
833	!----------------------------
834	! PE gain from ridging ice
835	!----------------------------
836	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 * &
837	& ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) + &
838	& hrmin(ji,jj,jl) * hrmin(ji,jj,jl) + &
839	& hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
840	!!(a3-b3)/(a-b) = aa+ab+bb
841	ENDIF
842	!
843	END DO
844	END DO
845	END DO
846
847	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
848	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
849	ksmooth = 1
850
851	!------------------------------------------------------------------------------!
852	! 4) Hibler (1979)' method
853	!------------------------------------------------------------------------------!
854	ELSE ! kstrngth ne 1: Hibler (1979) form
855	!
856	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1)
857	!
858	ksmooth = 1
859	!
860	ENDIF ! kstrngth
861	!
862	!------------------------------------------------------------------------------!
863	! 5) Impact of brine volume
864	!------------------------------------------------------------------------------!
865	! CAN BE REMOVED
866	IF( ln_icestr_bvf ) THEN
867	DO jj = 1, jpj
868	DO ji = 1, jpi
869	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
870	END DO
871	END DO
872	ENDIF
873	!
874	!------------------------------------------------------------------------------!
875	! 6) Smoothing ice strength
876	!------------------------------------------------------------------------------!
877	!
878	!-------------------
879	! Spatial smoothing
880	!-------------------
881	IF ( ksmooth == 1 ) THEN
882
883	CALL lbc_lnk( strength, 'T', 1. )
884
885	DO jj = 2, jpjm1
886	DO ji = 2, jpim1
887	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
888	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
889	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
890	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
891	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
892	ELSE
893	zworka(ji,jj) = 0._wp
894	ENDIF
895	END DO
896	END DO
897
898	DO jj = 2, jpjm1
899	DO ji = 2, jpim1
900	strength(ji,jj) = zworka(ji,jj)
901	END DO
902	END DO
903	CALL lbc_lnk( strength, 'T', 1. )
904
905	ENDIF
906
907	!--------------------
908	! Temporal smoothing
909	!--------------------
910	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
911	strp1(:,:) = 0.0
912	strp2(:,:) = 0.0
913	ENDIF
914
915	IF ( ksmooth == 2 ) THEN
916
917	CALL lbc_lnk( strength, 'T', 1. )
918
919	DO jj = 1, jpj - 1
920	DO ji = 1, jpi - 1
921	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
922	numts_rm = 1 ! number of time steps for the running mean
923	IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
924	IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
925	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
926	strp2(ji,jj) = strp1(ji,jj)
927	strp1(ji,jj) = strength(ji,jj)
928	strength(ji,jj) = zp
929
930	ENDIF
931	END DO
932	END DO
933
934	ENDIF ! ksmooth
935
936	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
937
938	CALL wrk_dealloc( jpi, jpj, zworka )
939	!
940	END SUBROUTINE lim_itd_me_icestrength
941
942	SUBROUTINE lim_itd_me_init
943	!!-------------------------------------------------------------------
944	!! * ROUTINE lim_itd_me_init *
945	!!
946	!! ** Purpose : Physical constants and parameters linked
947	!! to the mechanical ice redistribution
948	!!
949	!! ** Method : Read the namiceitdme namelist
950	!! and check the parameters values
951	!! called at the first timestep (nit000)
952	!!
953	!! ** input : Namelist namiceitdme
954	!!-------------------------------------------------------------------
955	INTEGER :: ios ! Local integer output status for namelist read
956	NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, &
957	& rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
958	& nn_partfun
959	!!-------------------------------------------------------------------
960	!
961	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
962	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
963	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
964
965	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
966	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
967	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
968	IF(lwm) WRITE ( numoni, namiceitdme )
969	!
970	IF (lwp) THEN ! control print
971	WRITE(numout,*)
972	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
973	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
974	WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs
975	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg
976	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft
977	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar
978	WRITE(numout,)' Equivalent to G for an exponential part function rn_astar = ', rn_astar
979	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar
980	WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting
981	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
982	WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
983	WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg
984	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun
985	ENDIF
986	!
987	END SUBROUTINE lim_itd_me_init
988
989	#else
990	!!----------------------------------------------------------------------
991	!! Default option Empty module NO LIM-3 sea-ice model
992	!!----------------------------------------------------------------------
993	CONTAINS
994	SUBROUTINE lim_itd_me ! Empty routines
995	END SUBROUTINE lim_itd_me
996	SUBROUTINE lim_itd_me_icestrength
997	END SUBROUTINE lim_itd_me_icestrength
998	SUBROUTINE lim_itd_me_init
999	END SUBROUTINE lim_itd_me_init
1000	#endif
1001	!!======================================================================
1002	END MODULE limitd_me

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

Download in other formats: