Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icedyn_rdgrft.F90

Last change on this file was 15549, checked in by clem, 2 years ago
commit ice namelist changes to be added to nemo4.2
Property svn:keywords set to `Id`
File size: 56.7 KB

Line
1	MODULE icedyn_rdgrft
2	!!======================================================================
3	!! * MODULE icedyn_rdgrft *
4	!! sea-ice : Mechanical impact on ice thickness distribution
5	!!======================================================================
6	!! History : ! 2006-02 (M. Vancoppenolle) Original code
7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
8	!!----------------------------------------------------------------------
9	#if defined key_si3
10	!!----------------------------------------------------------------------
11	!! 'key_si3' SI3 sea-ice model
12	!!----------------------------------------------------------------------
13	!! ice_dyn_rdgrft : ridging/rafting of sea ice
14	!! ice_dyn_rdgrft_init : initialization of ridging/rafting of sea ice
15	!! ice_strength : ice strength calculation
16	!!----------------------------------------------------------------------
17	USE dom_oce ! ocean domain
18	USE phycst ! physical constants (ocean directory)
19	USE sbc_oce , ONLY : sss_m, sst_m ! surface boundary condition: ocean fields
20	USE ice1D ! sea-ice: thermodynamics
21	USE ice ! sea-ice: variables
22	USE icetab ! sea-ice: 1D <==> 2D transformation
23	USE icevar ! sea-ice: operations
24	USE icectl ! sea-ice: control prints
25	!
26	USE in_out_manager ! I/O manager
27	USE iom ! I/O manager library
28	USE lib_mpp ! MPP library
29	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
30	USE lbclnk ! lateral boundary conditions (or mpp links)
31	USE timing ! Timing
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC ice_dyn_rdgrft ! called by icestp
37	PUBLIC ice_dyn_rdgrft_init ! called by icedyn
38	PUBLIC ice_strength ! called by icedyn_rhg_evp
39
40	INTEGER :: nice_str ! choice of the type of strength
41	! ! associated indices:
42	INTEGER, PARAMETER :: np_strh79 = 1 ! Hibler 79
43	INTEGER, PARAMETER :: np_strr75 = 2 ! Rothrock 75
44	INTEGER, PARAMETER :: np_strcst = 3 ! Constant value
45
46	! Variables shared among ridging subroutines
47	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_net ! net rate at which area is removed (1/s)
48	! ! (ridging ice area - area of new ridges) / dt
49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear
50	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges
51	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n
52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness
53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice
55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness
56	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging
57	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: araft ! participating ice rafting
58	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_i_2d
59	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d
60	!
61	REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi)
62	REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipliyer: (hi/hraft)
63	!
64	! namelist (namdyn_rdgrft)
65	LOGICAL :: ln_str_R75 ! ice strength parameterization: Rothrock 75
66	REAL(wp) :: rn_pe_rdg ! coef accounting for frictional dissipation
67	LOGICAL :: ln_str_CST ! ice strength parameterization: Constant
68	REAL(wp) :: rn_str ! constant value of ice strength
69	LOGICAL :: ln_str_smooth ! ice strength spatial smoothing
70	LOGICAL :: ln_distf_lin ! redistribution of ridged ice: linear (Hibler 1980)
71	LOGICAL :: ln_distf_exp ! redistribution of ridged ice: exponential
72	REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5)
73	REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging
74	LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975))
75	REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging
76	LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007))
77	REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function
78	LOGICAL :: ln_ridging ! ridging of ice or not
79	REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice
80	REAL(wp) :: rn_porordg ! initial porosity of ridges (0.3 regular value)
81	REAL(wp) :: rn_fsnwrdg ! fractional snow loss to the ocean during ridging
82	REAL(wp) :: rn_fpndrdg ! fractional pond loss to the ocean during ridging
83	LOGICAL :: ln_rafting ! rafting of ice or not
84	REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging
85	REAL(wp) :: rn_craft ! coefficient for smoothness of the hyperbolic tangent in rafting
86	REAL(wp) :: rn_fsnwrft ! fractional snow loss to the ocean during rafting
87	REAL(wp) :: rn_fpndrft ! fractional pond loss to the ocean during rafting
88	!
89	!! * Substitutions
90	# include "do_loop_substitute.h90"
91	!!----------------------------------------------------------------------
92	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
93	!! $Id$
94	!! Software governed by the CeCILL licence (./LICENSE)
95	!!----------------------------------------------------------------------
96	CONTAINS
97
98	INTEGER FUNCTION ice_dyn_rdgrft_alloc()
99	!!-------------------------------------------------------------------
100	!! * ROUTINE ice_dyn_rdgrft_alloc *
101	!!-------------------------------------------------------------------
102	ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , &
103	& apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
104	& hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , &
105	& ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
106
107	CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
108	IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_alloc: failed to allocate arrays' )
109	!
110	END FUNCTION ice_dyn_rdgrft_alloc
111
112
113	SUBROUTINE ice_dyn_rdgrft( kt )
114	!!-------------------------------------------------------------------
115	!! * ROUTINE ice_dyn_rdgrft *
116	!!
117	!! ** Purpose : computes the mechanical redistribution of ice thickness
118	!!
119	!! ** Method : Steps :
120	!! 0) Identify grid cells with ice
121	!! 1) Calculate closing rate, divergence and opening
122	!! 2) Identify grid cells with ridging
123	!! 3) Start ridging iterations
124	!! - prep = ridged and rafted ice + closing_gross
125	!! - shift = move ice from one category to another
126	!!
127	!! ** Details
128	!! step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
129	!! The energy dissipation rate is equal to the net closing rate times the ice strength.
130	!!
131	!! step3: The gross closing rate is equal to the first two terms (open
132	!! water closing and thin ice ridging) without the third term
133	!! (thick, newly ridged ice).
134	!!
135	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
136	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
137	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
138	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
139	!! Bitz et al., JGR, 2001
140	!! Amundrud and Melling, JGR 2005
141	!! Babko et al., JGR 2002
142	!!
143	!! This routine is based on CICE code and authors William H. Lipscomb,
144	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
145	!!-------------------------------------------------------------------
146	INTEGER, INTENT(in) :: kt ! number of iteration
147	!!
148	INTEGER :: ji, jj, jk, jl ! dummy loop index
149	INTEGER :: iter, iterate_ridging ! local integer
150	INTEGER :: ipti ! local integer
151	REAL(wp) :: zfac ! local scalar
152	INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not
153	REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i
154	REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology)
155	!
156	INTEGER, PARAMETER :: jp_itermax = 20
157	!!-------------------------------------------------------------------
158	! controls
159	IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing
160	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
161	IF( ln_icediachk ) CALL ice_cons2D (0, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
162
163	IF( kt == nit000 ) THEN
164	IF(lwp) WRITE(numout,*)
165	IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
166	IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
167	ENDIF
168
169	!--------------------------------
170	! 0) Identify grid cells with ice
171	!--------------------------------
172	at_i(:,:) = SUM( a_i, dim=3 )
173	!
174	npti = 0 ; nptidx(:) = 0
175	ipti = 0 ; iptidx(:) = 0
176	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
177	IF ( at_i(ji,jj) > epsi10 ) THEN
178	npti = npti + 1
179	nptidx( npti ) = (jj - 1) * jpi + ji
180	ENDIF
181	END_2D
182
183	!--------------------------------------------------------
184	! 1) Dynamical inputs (closing rate, divergence, opening)
185	!--------------------------------------------------------
186	IF( npti > 0 ) THEN
187
188	! just needed here
189	CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i )
190	CALL tab_2d_1d( npti, nptidx(1:npti), zconv (1:npti) , rdg_conv )
191	! needed here and in the iteration loop
192	CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i)
193	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
194	CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i )
195	CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
196
197	DO ji = 1, npti
198	! closing_net = rate at which open water area is removed + ice area removed by ridging
199	! - ice area added in new ridges
200	IF( ln_rhg_EVP .OR. ln_rhg_VP ) &
201	& closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
202	IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji)
203	!
204	IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough
205	! ! to give asum = 1.0 after ridging
206	! Opening rate (non-negative) that will give asum = 1.0 after ridging.
207	opning(ji) = closing_net(ji) + zdivu(ji)
208	END DO
209	!
210	!------------------------------------
211	! 2) Identify grid cells with ridging
212	!------------------------------------
213	CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
214
215	DO ji = 1, npti
216	IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
217	ipti = ipti + 1
218	iptidx (ipti) = nptidx (ji)
219	! adjust to new indices
220	a_i_2d (ipti,:) = a_i_2d (ji,:)
221	v_i_2d (ipti,:) = v_i_2d (ji,:)
222	ato_i_1d (ipti) = ato_i_1d (ji)
223	closing_net(ipti) = closing_net(ji)
224	zdivu (ipti) = zdivu (ji)
225	opning (ipti) = opning (ji)
226	ENDIF
227	END DO
228
229	ENDIF
230
231	! grid cells with ridging
232	nptidx(:) = iptidx(:)
233	npti = ipti
234
235	!-----------------
236	! 3) Start ridging
237	!-----------------
238	IF( npti > 0 ) THEN
239
240	CALL ice_dyn_1d2d( 1 ) ! --- Move to 1D arrays --- !
241
242	iter = 1
243	iterate_ridging = 1
244	! !----------------------!
245	DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) ! ridging iterations !
246	! !----------------------!
247	! Calculate participation function (apartf)
248	! and transfer function
249	! and closing_gross (+correction on opening)
250	CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
251
252	! Redistribute area, volume, and energy between categories
253	CALL rdgrft_shift
254
255	! Do we keep on iterating?
256	!-------------------------
257	! Check whether a_i + ato_i = 0
258	! If not, because the closing and opening rates were reduced above, ridge again with new rates
259	iterate_ridging = 0
260	DO ji = 1, npti
261	zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
262	IF( ABS( zfac ) < epsi10 ) THEN
263	closing_net(ji) = 0._wp
264	opning (ji) = 0._wp
265	ato_i_1d (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
266	ELSE
267	iterate_ridging = 1
268	zdivu (ji) = zfac * r1_Dt_ice
269	closing_net(ji) = MAX( 0._wp, -zdivu(ji) )
270	opning (ji) = MAX( 0._wp, zdivu(ji) )
271	ENDIF
272	END DO
273	!
274	iter = iter + 1
275	IF( iter > jp_itermax ) CALL ctl_stop( 'STOP', 'icedyn_rdgrft: non-converging ridging scheme' )
276	!
277	END DO
278
279	CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- !
280
281	ENDIF
282
283	CALL ice_var_agg( 1 )
284
285	! controls
286	IF( sn_cfctl%l_prtctl ) CALL ice_prt3D('icedyn_rdgrft') ! prints
287	IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints
288	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
289	IF( ln_icediachk ) CALL ice_cons2D (1, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
290	IF( ln_timing ) CALL timing_stop ('icedyn_rdgrft') ! timing
291	!
292	END SUBROUTINE ice_dyn_rdgrft
293
294
295	SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
296	!!-------------------------------------------------------------------
297	!! * ROUTINE rdgrft_prep *
298	!!
299	!! ** Purpose : preparation for ridging calculations
300	!!
301	!! ** Method : Compute the thickness distribution of the ice and open water
302	!! participating in ridging and of the resulting ridges.
303	!!-------------------------------------------------------------------
304	REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net
305	REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i
306	!!
307	INTEGER :: ji, jl ! dummy loop indices
308	REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar
309	REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse
310	REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness
311	REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n
312	!--------------------------------------------------------------------
313
314	z1_gstar = 1._wp / rn_gstar
315	z1_astar = 1._wp / rn_astar
316
317	! ! Ice thickness needed for rafting
318	! In single precision there were floating point invalids due a sqrt of zhi which happens to have negative values
319	! To solve that an extra check about the value of pv_i was added.
320	! Although adding this condition is safe, the double definition (one for single other for double) has been kept to preserve the results of the sette test.
321	#if defined key_single
322
323	WHERE( pa_i(1:npti,:) > epsi10 .and. pv_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
324	#else
325	WHERE( pa_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
326	#endif
327	ELSEWHERE ; zhi(1:npti,:) = 0._wp
328	END WHERE
329
330	! 1) Participation function (apartf): a(h) = b(h).g(h)
331	!-----------------------------------------------------------------
332	! Compute the participation function = total area lost due to ridging/closing
333	! This is analogous to
334	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
335	! assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
336	!
337	! apartf = integrating b(h)g(h) between the category boundaries
338	! apartf is always >= 0 and SUM(apartf(0:jpl))=1
339	!-----------------------------------------------------------------
340	!
341	! Compute total area of ice plus open water.
342	! This is in general not equal to one because of divergence during transport
343	zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
344	!
345	WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)
346	ELSEWHERE ; z1_asum(1:npti) = 0._wp
347	END WHERE
348	!
349	! Compute cumulative thickness distribution function
350	! Compute the cumulative thickness distribution function zGsum,
351	! where zGsum(n) is the fractional area in categories 0 to n.
352	! initial value (in h = 0) = open water area
353	zGsum(1:npti,-1) = 0._wp
354	zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
355	DO jl = 1, jpl
356	zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok (and not jpl)
357	END DO
358	!
359	IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975)
360	DO jl = 0, jpl
361	DO ji = 1, npti
362	IF ( zGsum(ji,jl) < rn_gstar ) THEN
363	apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
364	& ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
365	ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
366	apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * &
367	& ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar )
368	ELSE
369	apartf(ji,jl) = 0._wp
370	ENDIF
371	END DO
372	END DO
373	!
374	ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007)
375	!
376	zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
377	DO jl = -1, jpl
378	DO ji = 1, npti
379	zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
380	END DO
381	END DO
382	DO jl = 0, jpl
383	DO ji = 1, npti
384	apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
385	END DO
386	END DO
387	!
388	ENDIF
389
390	! !--- Ridging and rafting participation concentrations
391	IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting
392	DO jl = 1, jpl
393	DO ji = 1, npti
394	aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
395	araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
396	END DO
397	END DO
398	ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone
399	DO jl = 1, jpl
400	DO ji = 1, npti
401	aridge(ji,jl) = apartf(ji,jl)
402	araft (ji,jl) = 0._wp
403	END DO
404	END DO
405	ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone
406	DO jl = 1, jpl
407	DO ji = 1, npti
408	aridge(ji,jl) = 0._wp
409	araft (ji,jl) = apartf(ji,jl)
410	END DO
411	END DO
412	ELSE !- no ridging & no rafting
413	DO jl = 1, jpl
414	DO ji = 1, npti
415	aridge(ji,jl) = 0._wp
416	araft (ji,jl) = 0._wp
417	END DO
418	END DO
419	ENDIF
420
421	! 2) Transfer function
422	!-----------------------------------------------------------------
423	! Compute max and min ridged ice thickness for each ridging category.
424	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
425	!
426	! This parameterization is a modified version of Hibler (1980).
427	! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
428	! and for very thick ridging ice must be >= hrdg_hi_min*hi
429	!
430	! The minimum ridging thickness, hrmin, is equal to 2*hi
431	! (i.e., rafting) and for very thick ridging ice is
432	! constrained by hrmin <= (zhmean + hi)/2.
433	!
434	! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
435	!
436	! These modifications have the effect of reducing the ice strength
437	! (relative to the Hibler formulation) when very thick ice is ridging.
438	!
439	! zaksum = net area removed/ total area removed
440	! where total area removed = area of ice that ridges
441	! net area removed = total area removed - area of new ridges
442	!-----------------------------------------------------------------
443	zfac = 1._wp / hi_hrft
444	zaksum(1:npti) = apartf(1:npti,0)
445	!
446	DO jl = 1, jpl
447	DO ji = 1, npti
448	IF ( apartf(ji,jl) > 0._wp ) THEN
449	zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
450	hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
451	hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
452	hraft (ji,jl) = zhi(ji,jl) * zfac
453	hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
454	!
455	! Normalization factor : zaksum, ensures mass conservation
456	zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) &
457	& + araft (ji,jl) * ( 1._wp - hi_hrft )
458	ELSE
459	hrmin (ji,jl) = 0._wp
460	hrmax (ji,jl) = 0._wp
461	hraft (ji,jl) = 0._wp
462	hi_hrdg(ji,jl) = 1._wp
463	ENDIF
464	END DO
465	END DO
466	!
467	! 3) closing_gross
468	!-----------------
469	! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
470	! NOTE: 0 < aksum <= 1
471	WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
472	ELSEWHERE ; closing_gross(1:npti) = 0._wp
473	END WHERE
474
475	! correction to closing rate if excessive ice removal
476	!----------------------------------------------------
477	! Reduce the closing rate if more than 100% of any ice category would be removed
478	! Reduce the opening rate in proportion
479	DO jl = 1, jpl
480	DO ji = 1, npti
481	zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
482	IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
483	closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
484	ENDIF
485	END DO
486	END DO
487
488	! 4) correction to opening if excessive open water removal
489	!---------------------------------------------------------
490	! Reduce the closing rate if more than 100% of the open water would be removed
491	! Reduce the opening rate in proportion
492	DO ji = 1, npti
493	zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
494	IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
495	opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
496	ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
497	opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
498	ENDIF
499	END DO
500	!
501	END SUBROUTINE rdgrft_prep
502
503
504	SUBROUTINE rdgrft_shift
505	!!-------------------------------------------------------------------
506	!! * ROUTINE rdgrft_shift *
507	!!
508	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
509	!!
510	!! ** Method : Remove area, volume, and energy from each ridging category
511	!! and add to thicker ice categories.
512	!!-------------------------------------------------------------------
513	!
514	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
515	REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2
516	REAL(wp) :: vsw ! vol of water trapped into ridges
517	REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted
518	REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1
519	REAL(wp) :: airft1, oirft1, aprft1
520	REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg ! area etc of new ridges
521	REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft ! area etc of rafted ice
522	!
523	REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges
524	REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2
525	REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a
526	REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i)
527	!
528	REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice
529	REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice
530	REAL(wp), DIMENSION(jpij,nlay_s) :: esrdg ! enth*volume of new ridges
531	REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges
532	!
533	INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation
534	LOGICAL , DIMENSION(jpij) :: ll_shift ! logical for doing calculation or not
535	!!-------------------------------------------------------------------
536	!
537	zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume
538	!
539	! 1) Change in open water area due to closing and opening
540	!--------------------------------------------------------
541	DO ji = 1, npti
542	ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice )
543	END DO
544
545	! 2) compute categories in which ice is removed (jl1)
546	!----------------------------------------------------
547	DO jl1 = 1, jpl
548
549	IF( nn_icesal /= 2 ) THEN
550	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
551	ENDIF
552
553	DO ji = 1, npti
554
555	! set logical to true when ridging
556	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ; ll_shift(ji) = .TRUE.
557	ELSE ; ll_shift(ji) = .FALSE.
558	ENDIF
559
560	IF( ll_shift(ji) ) THEN ! only if ice is ridging
561
562	IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
563	ELSE ; z1_ai(ji) = 0._wp
564	ENDIF
565
566	! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
567	airdg1 = aridge(ji,jl1) * closing_gross(ji) * rDt_ice
568	airft1 = araft (ji,jl1) * closing_gross(ji) * rDt_ice
569
570	airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
571	airft2(ji) = airft1 * hi_hrft
572
573	! ridging /rafting fractions
574	afrdg = airdg1 * z1_ai(ji)
575	afrft = airft1 * z1_ai(ji)
576
577	! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
578	IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg
579	ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0
580	ELSE ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
581	ENDIF
582	ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?)
583
584	! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
585	virdg1 = v_i_2d (ji,jl1) * afrdg
586	virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw
587	vsrdg(ji) = v_s_2d (ji,jl1) * afrdg
588	sirdg1 = sv_i_2d(ji,jl1) * afrdg
589	sirdg2(ji) = sv_i_2d(ji,jl1) * afrdg + vsw * sss_1d(ji)
590	oirdg1 = oa_i_2d(ji,jl1) * afrdg
591	oirdg2(ji) = oa_i_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
592
593	virft(ji) = v_i_2d (ji,jl1) * afrft
594	vsrft(ji) = v_s_2d (ji,jl1) * afrft
595	sirft(ji) = sv_i_2d(ji,jl1) * afrft
596	oirft1 = oa_i_2d(ji,jl1) * afrft
597	oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft
598
599	IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
600	aprdg1 = a_ip_2d(ji,jl1) * afrdg
601	aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
602	vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
603	aprft1 = a_ip_2d(ji,jl1) * afrft
604	aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
605	vprft (ji) = v_ip_2d(ji,jl1) * afrft
606	IF ( ln_pnd_lids ) THEN
607	vlrdg (ji) = v_il_2d(ji,jl1) * afrdg
608	vlrft (ji) = v_il_2d(ji,jl1) * afrft
609	ENDIF
610	ENDIF
611
612	! Ice-ocean exchanges associated with ice porosity
613	wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_Dt_ice ! increase in ice volume due to seawater frozen in voids
614	sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice
615	hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2]
616
617	! Put the snow lost by ridging into the ocean
618	! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
619	wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean
620	& + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
621
622	! virtual salt flux to keep salinity constant
623	IF( nn_icesal /= 2 ) THEN
624	sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i
625	sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice & ! put back sss_m into the ocean
626	& - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice ! and get s_i from the ocean
627	ENDIF
628
629	! Remove area, volume of new ridge to each category jl1
630	!------------------------------------------------------
631	a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1
632	v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji)
633	v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
634	sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1 - sirft(ji)
635	oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1 - oirft1
636	IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
637	a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1
638	v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
639	IF ( ln_pnd_lids ) THEN
640	v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji)
641	ENDIF
642	ENDIF
643	ENDIF
644
645	END DO ! ji
646
647	! special loop for e_s because of layers jk
648	DO jk = 1, nlay_s
649	DO ji = 1, npti
650	IF( ll_shift(ji) ) THEN
651	! Compute ridging /rafting fractions
652	afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
653	afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
654	! Compute ridging /rafting ice and new ridges for es
655	esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
656	esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
657	! Put the snow lost by ridging into the ocean
658	hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2)
659	& - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
660	!
661	! Remove energy of new ridge to each category jl1
662	!-------------------------------------------------
663	ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
664	ENDIF
665	END DO
666	END DO
667
668	! special loop for e_i because of layers jk
669	DO jk = 1, nlay_i
670	DO ji = 1, npti
671	IF( ll_shift(ji) ) THEN
672	! Compute ridging /rafting fractions
673	afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
674	afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
675	! Compute ridging ice and new ridges for ei
676	eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
677	eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
678	!
679	! Remove energy of new ridge to each category jl1
680	!-------------------------------------------------
681	ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
682	ENDIF
683	END DO
684	END DO
685
686	! 3) compute categories in which ice is added (jl2)
687	!--------------------------------------------------
688	itest_rdg(1:npti) = 0
689	itest_rft(1:npti) = 0
690	DO jl2 = 1, jpl
691	!
692	DO ji = 1, npti
693
694	IF( ll_shift(ji) ) THEN
695
696	! Compute the fraction of ridged ice area and volume going to thickness category jl2
697	IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
698	hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
699	hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
700	farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
701	fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
702	!
703	itest_rdg(ji) = 1 ! test for conservation
704	ELSE
705	farea = 0._wp
706	fvol(ji) = 0._wp
707	ENDIF
708
709	! Compute the fraction of rafted ice area and volume going to thickness category jl2
710	IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN
711	zswitch(ji) = 1._wp
712	!
713	itest_rft(ji) = 1 ! test for conservation
714	ELSE
715	zswitch(ji) = 0._wp
716	ENDIF
717	!
718	! Patch to ensure perfect conservation if ice thickness goes mad
719	! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
720	! Then ice volume is removed from one category but the ridging/rafting scheme
721	! does not know where to move it, leading to a conservation issue.
722	IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF
723	IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp
724	!
725	! Add area, volume of new ridge to category jl2
726	!----------------------------------------------
727	a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) )
728	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) )
729	v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
730	sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
731	v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + &
732	& vsrft (ji) * rn_fsnwrft * zswitch(ji) )
733	IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
734	v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) &
735	& + vprft (ji) * rn_fpndrft * zswitch(ji) )
736	a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea &
737	& + aprft2(ji) * rn_fpndrft * zswitch(ji) )
738	IF ( ln_pnd_lids ) THEN
739	v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + ( vlrdg(ji) * rn_fpndrdg * fvol (ji) &
740	& + vlrft(ji) * rn_fpndrft * zswitch(ji) )
741	ENDIF
742	ENDIF
743
744	ENDIF
745
746	END DO
747	! Add snow energy of new ridge to category jl2
748	!---------------------------------------------
749	DO jk = 1, nlay_s
750	DO ji = 1, npti
751	IF( ll_shift(ji) ) &
752	& ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji) + &
753	& esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
754	END DO
755	END DO
756	! Add ice energy of new ridge to category jl2
757	!--------------------------------------------
758	DO jk = 1, nlay_i
759	DO ji = 1, npti
760	IF( ll_shift(ji) ) &
761	& ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)
762	END DO
763	END DO
764	!
765	END DO ! jl2
766	!
767	END DO ! jl1
768	!
769	! roundoff errors
770	!----------------
771	! In case ridging/rafting lead to very small negative values (sometimes it happens)
772	CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
773	!
774	END SUBROUTINE rdgrft_shift
775
776
777	SUBROUTINE ice_strength
778	!!----------------------------------------------------------------------
779	!! * ROUTINE ice_strength *
780	!!
781	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
782	!!
783	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
784	!! dissipated per unit area removed from the ice pack under compression,
785	!! and assumed proportional to the change in potential energy caused
786	!! by ridging. Note that only Hibler's formulation is stable and that
787	!! ice strength has to be smoothed
788	!!----------------------------------------------------------------------
789	INTEGER :: ji, jj, jl ! dummy loop indices
790	REAL(wp) :: z1_3 ! local scalars
791	REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
792	!!clem
793	LOGICAL :: ln_str_R75
794	REAL(wp) :: zhi, zcp, rn_pe_rdg
795	REAL(wp), DIMENSION(jpij) :: zstrength, zaksum ! strength in 1D
796	!!----------------------------------------------------------------------
797	!
798	SELECT CASE( nice_str ) !--- Set which ice strength is chosen
799
800	CASE ( np_strr75 ) !== Rothrock(1975)'s method ==!
801
802	! these 2 param should be defined once for all at the 1st time step
803	zcp = 0.5_wp * grav * (rho0-rhoi) * rhoi * r1_rho0 ! proport const for PE
804	!
805	strength(:,:) = 0._wp
806	!
807	! Identify grid cells with ice
808	at_i(:,:) = SUM( a_i, dim=3 )
809	npti = 0 ; nptidx(:) = 0
810	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
811	IF ( at_i(ji,jj) > epsi10 ) THEN
812	npti = npti + 1
813	nptidx( npti ) = (jj - 1) * jpi + ji
814	ENDIF
815	END_2D
816
817	IF( npti > 0 ) THEN
818	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
819	CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i )
820	CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
821	CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti) , strength )
822
823	CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
824	!
825	zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module
826	DO jl = 1, jpl
827	DO ji = 1, npti
828	IF ( apartf(ji,jl) > 0._wp ) THEN
829	zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) &
830	& + araft (ji,jl) * ( 1._wp - hi_hrft )
831	ENDIF
832	END DO
833	END DO
834	!
835	z1_3 = 1._wp / 3._wp
836	DO jl = 1, jpl
837	DO ji = 1, npti
838	!
839	IF( apartf(ji,jl) > 0._wp ) THEN
840	!
841	IF( a_i_2d(ji,jl) > epsi10 ) THEN ; zhi = v_i_2d(ji,jl) / a_i_2d(ji,jl)
842	ELSE ; zhi = 0._wp
843	ENDIF
844	zstrength(ji) = zstrength(ji) - apartf(ji,jl) * zhi * zhi ! PE loss from deforming ice
845	zstrength(ji) = zstrength(ji) + 2._wp * araft (ji,jl) * zhi * zhi ! PE gain from rafting ice
846	zstrength(ji) = zstrength(ji) + aridge(ji,jl) * hi_hrdg(ji,jl) * z1_3 * & ! PE gain from ridging ice
847	& ( hrmax(ji,jl) * hrmax(ji,jl) + & ! (a3-b3)/(a-b) = aa+ab+bb
848	& hrmin(ji,jl) * hrmin(ji,jl) + &
849	& hrmax(ji,jl) * hrmin(ji,jl) )
850	ENDIF
851	!
852	END DO
853	END DO
854	!
855	zstrength(1:npti) = rn_pe_rdg * zcp * zstrength(1:npti) / zaksum(1:npti)
856	!
857	CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength )
858	!
859	ENDIF
860	!
861	CASE ( np_strh79 ) !== Hibler(1979)'s method ==!
862	strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
863	!
864	CASE ( np_strcst ) !== Constant strength ==!
865	strength(:,:) = rn_str
866	!
867	END SELECT
868	!
869	IF( ln_str_smooth ) THEN !--- Spatial smoothing
870	DO_2D( 0, 0, 0, 0 )
871	IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
872	zworka(ji,jj) = ( 4._wp * strength(ji,jj) &
873	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
874	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
875	& ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
876	ELSE
877	zworka(ji,jj) = 0._wp
878	ENDIF
879	END_2D
880
881	DO_2D( 0, 0, 0, 0 )
882	strength(ji,jj) = zworka(ji,jj)
883	END_2D
884	CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
885	!
886	ENDIF
887	!
888	END SUBROUTINE ice_strength
889
890
891	SUBROUTINE ice_dyn_1d2d( kn )
892	!!-----------------------------------------------------------------------
893	!! * ROUTINE ice_dyn_1d2d *
894	!!
895	!! ** Purpose : move arrays from 1d to 2d and the reverse
896	!!-----------------------------------------------------------------------
897	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
898	!
899	INTEGER :: jl, jk ! dummy loop indices
900	!!-----------------------------------------------------------------------
901	!
902	SELECT CASE( kn )
903	! !---------------------!
904	CASE( 1 ) !== from 2D to 1D ==!
905	! !---------------------!
906	! fields used but not modified
907	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
908	CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
909	! the following fields are modified in this routine
910	!!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
911	!!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
912	!!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
913	CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
914	CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
915	CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
916	CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
917	CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
918	CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
919	DO jl = 1, jpl
920	DO jk = 1, nlay_s
921	CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
922	END DO
923	DO jk = 1, nlay_i
924	CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
925	END DO
926	END DO
927	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
928	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
929	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
930	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) )
931	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
932	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
933	!
934	! !---------------------!
935	CASE( 2 ) !== from 1D to 2D ==!
936	! !---------------------!
937	CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
938	CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
939	CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
940	CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
941	CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
942	CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
943	CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
944	CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
945	CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) )
946	DO jl = 1, jpl
947	DO jk = 1, nlay_s
948	CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
949	END DO
950	DO jk = 1, nlay_i
951	CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
952	END DO
953	END DO
954	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
955	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
956	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
957	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) )
958	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
959	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
960	!
961	END SELECT
962	!
963	END SUBROUTINE ice_dyn_1d2d
964
965
966	SUBROUTINE ice_dyn_rdgrft_init
967	!!-------------------------------------------------------------------
968	!! * ROUTINE ice_dyn_rdgrft_init *
969	!!
970	!! ** Purpose : Physical constants and parameters linked
971	!! to the mechanical ice redistribution
972	!!
973	!! ** Method : Read the namdyn_rdgrft namelist
974	!! and check the parameters values
975	!! called at the first timestep (nit000)
976	!!
977	!! ** input : Namelist namdyn_rdgrft
978	!!-------------------------------------------------------------------
979	INTEGER :: ios, ioptio ! Local integer output status for namelist read
980	!!
981	NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, ln_str_R75, rn_pe_rdg, ln_str_CST, rn_str, ln_str_smooth, &
982	& ln_distf_lin, ln_distf_exp, rn_murdg, rn_csrdg, &
983	& ln_partf_lin, rn_gstar, ln_partf_exp, rn_astar, &
984	& ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, &
985	& ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft
986	!!-------------------------------------------------------------------
987	!
988	READ ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
989	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
990	READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
991	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' )
992	IF(lwm) WRITE ( numoni, namdyn_rdgrft )
993	!
994	IF (lwp) THEN ! control print
995	WRITE(numout,*)
996	WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
997	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
998	WRITE(numout,*) ' Namelist namdyn_rdgrft:'
999	WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79
1000	WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar
1001	WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg
1002	WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75
1003	WRITE(numout,*) ' coef accounting for frictional dissipation rn_pe_rdg = ', rn_pe_rdg
1004	WRITE(numout,*) ' ice strength parameterization Constant ln_str_CST = ', ln_str_CST
1005	WRITE(numout,*) ' ice strength value rn_str = ', rn_str
1006	WRITE(numout,*) ' spatial smoothing of the strength ln_str_smooth= ', ln_str_smooth
1007	WRITE(numout,*) ' redistribution of ridged ice: linear (Hibler 1980) ln_distf_lin = ', ln_distf_lin
1008	WRITE(numout,*) ' redistribution of ridged ice: exponential ln_distf_exp = ', ln_distf_exp
1009	WRITE(numout,*) ' e-folding scale of ridged ice rn_murdg = ', rn_murdg
1010	WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg
1011	WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
1012	WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar
1013	WRITE(numout,*) ' Exponential ridging participation function ln_partf_exp = ', ln_partf_exp
1014	WRITE(numout,) ' Equivalent to G for an exponential function rn_astar = ', rn_astar
1015	WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging
1016	WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar
1017	WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg
1018	WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg
1019	WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg
1020	WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting
1021	WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
1022	WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
1023	WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft
1024	WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft
1025	ENDIF
1026	!
1027	ioptio = 0
1028	IF( ln_str_H79 ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strh79 ; ENDIF
1029	IF( ln_str_R75 ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strr75 ; ENDIF
1030	IF( ln_str_CST ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strcst ; ENDIF
1031	IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_rdgrft_init: one and only one ice strength option has to be defined ' )
1032	!
1033	IF ( ( ln_distf_lin .AND. ln_distf_exp ) .OR. ( .NOT.ln_distf_lin .AND. .NOT.ln_distf_exp ) ) THEN
1034	CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one redistribution function (ln_distf_lin or ln_distf_exp)' )
1035	ENDIF
1036	!!clem
1037	IF( ln_distf_exp ) CALL ctl_stop( 'ice_dyn_rdgrft_init: exponential redistribution function not coded yet (ln_distf_exp)' )
1038	!
1039	IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
1040	CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
1041	ENDIF
1042	!
1043	IF( .NOT. ln_icethd ) THEN
1044	rn_porordg = 0._wp
1045	rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
1046	rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
1047	IF( lwp ) THEN
1048	WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed'
1049	WRITE(numout,*) ' rn_porordg = ', rn_porordg
1050	WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg
1051	WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg
1052	WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft
1053	WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft
1054	ENDIF
1055	ENDIF
1056	! ! allocate arrays
1057	IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
1058	!
1059	END SUBROUTINE ice_dyn_rdgrft_init
1060
1061	#else
1062	!!----------------------------------------------------------------------
1063	!! Default option Empty module NO SI3 sea-ice model
1064	!!----------------------------------------------------------------------
1065	#endif
1066
1067	!!======================================================================
1068	END MODULE icedyn_rdgrft

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

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

Context Navigation

source: NEMO/trunk/src/ICE/icedyn_rdgrft.F90

Download in other formats: