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

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

limitd_me.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

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

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

Download in other formats: