Context Navigation

icedyn_rdgrft.F90

Last change on this file was 12937, checked in by dancopsey, 4 years ago

Merge in Clem's branch. It was originally here:

svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch

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

Download in other formats: