Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icedyn_rdgrft.F90

Last change on this file was 13442, checked in by dancopsey, 4 years ago
Put majority of code in to get melt pond freezing to feed back into ice thickness.
File size: 53.4 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, lhprdg, svprdg ! area etc of new ridges
506	REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft, svprft ! 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	lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg
581	svprdg(ji) = sv_ip_2d(ji,jl1) * afrdg
582	aprft1 = a_ip_2d(ji,jl1) * afrft
583	aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
584	vprft (ji) = v_ip_2d(ji,jl1) * afrft
585	lhprft(ji) = lh_ip_2d(ji,jl1) * afrft
586	svprft(ji) = sv_ip_2d(ji,jl1) * afrft
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	lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji)
617	sv_ip_2d(ji,jl1) = sv_ip_2d(ji,jl1) - svprdg(ji) - svprft(ji)
618	ENDIF
619	ENDIF
620
621	END DO ! ji
622
623	! special loop for e_s because of layers jk
624	DO jk = 1, nlay_s
625	DO ji = 1, npti
626	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
627	! Compute ridging /rafting fractions
628	afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
629	afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
630	! Compute ridging /rafting ice and new ridges for es
631	esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
632	esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
633	! Put the snow lost by ridging into the ocean
634	hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2)
635	& - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
636	!
637	! Remove energy of new ridge to each category jl1
638	!-------------------------------------------------
639	ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
640	ENDIF
641	END DO
642	END DO
643
644	! special loop for e_i because of layers jk
645	DO jk = 1, nlay_i
646	DO ji = 1, npti
647	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
648	! Compute ridging /rafting fractions
649	afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
650	afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
651	! Compute ridging ice and new ridges for ei
652	eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
653	eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
654	!
655	! Remove energy of new ridge to each category jl1
656	!-------------------------------------------------
657	ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
658	ENDIF
659	END DO
660	END DO
661
662	! 3) compute categories in which ice is added (jl2)
663	!--------------------------------------------------
664	itest_rdg(1:npti) = 0
665	itest_rft(1:npti) = 0
666	DO jl2 = 1, jpl
667	!
668	DO ji = 1, npti
669
670	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
671
672	! Compute the fraction of ridged ice area and volume going to thickness category jl2
673	IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
674	hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
675	hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
676	farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
677	fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
678	!
679	itest_rdg(ji) = 1 ! test for conservation
680	ELSE
681	farea = 0._wp
682	fvol(ji) = 0._wp
683	ENDIF
684
685	! Compute the fraction of rafted ice area and volume going to thickness category jl2
686	IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN
687	zswitch(ji) = 1._wp
688	!
689	itest_rft(ji) = 1 ! test for conservation
690	ELSE
691	zswitch(ji) = 0._wp
692	ENDIF
693	!
694	! Patch to ensure perfect conservation if ice thickness goes mad
695	! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
696	! Then ice volume is removed from one category but the ridging/rafting scheme
697	! does not know where to move it, leading to a conservation issue.
698	IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF
699	IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp
700	!
701	! Add area, volume of new ridge to category jl2
702	!----------------------------------------------
703	a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) )
704	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) )
705	v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
706	sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
707	v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + &
708	& vsrft (ji) * rn_fsnwrft * zswitch(ji) )
709	IF ( ln_pnd_H12 ) THEN
710	v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) &
711	& + vprft (ji) * rn_fpndrft * zswitch(ji) )
712	a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea &
713	& + aprft2(ji) * rn_fpndrft * zswitch(ji) )
714	lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + ( lhprdg (ji) * rn_fpndrdg * fvol (ji) &
715	& + lhprft (ji) * rn_fpndrft * zswitch(ji) )
716	sv_ip_2d (ji,jl2) = sv_ip_2d(ji,jl2) + ( svprdg (ji) * rn_fpndrdg * fvol (ji) &
717	& + svprft (ji) * rn_fpndrft * zswitch(ji) )
718	ENDIF
719
720	ENDIF
721
722	END DO
723	! Add snow energy of new ridge to category jl2
724	!---------------------------------------------
725	DO jk = 1, nlay_s
726	DO ji = 1, npti
727	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) &
728	& ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji) + &
729	& esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
730	END DO
731	END DO
732	! Add ice energy of new ridge to category jl2
733	!--------------------------------------------
734	DO jk = 1, nlay_i
735	DO ji = 1, npti
736	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) &
737	& ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)
738	END DO
739	END DO
740	!
741	END DO ! jl2
742	!
743	END DO ! jl1
744	!
745	! roundoff errors
746	!----------------
747	! In case ridging/rafting lead to very small negative values (sometimes it happens)
748	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, lh_ip_2d, ,sv_ip_2d, ze_s_2d, ze_i_2d )
749	!
750	END SUBROUTINE rdgrft_shift
751
752
753	SUBROUTINE ice_strength
754	!!----------------------------------------------------------------------
755	!! * ROUTINE ice_strength *
756	!!
757	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
758	!!
759	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
760	!! dissipated per unit area removed from the ice pack under compression,
761	!! and assumed proportional to the change in potential energy caused
762	!! by ridging. Note that only Hibler's formulation is stable and that
763	!! ice strength has to be smoothed
764	!!----------------------------------------------------------------------
765	INTEGER :: ji, jj, jl ! dummy loop indices
766	INTEGER :: ismooth ! smoothing the resistance to deformation
767	INTEGER :: itframe ! number of time steps for the P smoothing
768	REAL(wp) :: zp, z1_3 ! local scalars
769	REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
770	REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps
771	!!----------------------------------------------------------------------
772	! !--------------------------------------------------!
773	IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method !
774	! !--------------------------------------------------!
775	strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
776	ismooth = 1
777	! !--------------------------------------------------!
778	ELSE ! Zero strength !
779	! !--------------------------------------------------!
780	strength(:,:) = 0._wp
781	ismooth = 0
782	ENDIF
783	! !--------------------------------------------------!
784	SELECT CASE( ismooth ) ! Smoothing ice strength !
785	! !--------------------------------------------------!
786	CASE( 1 ) !--- Spatial smoothing
787	DO jj = 2, jpjm1
788	DO ji = 2, jpim1
789	IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
790	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
791	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
792	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
793	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
794	ELSE
795	zworka(ji,jj) = 0._wp
796	ENDIF
797	END DO
798	END DO
799
800	DO jj = 2, jpjm1
801	DO ji = 2, jpim1
802	strength(ji,jj) = zworka(ji,jj)
803	END DO
804	END DO
805	CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
806	!
807	CASE( 2 ) !--- Temporal smoothing
808	IF ( kt_ice == nit000 ) THEN
809	zstrp1(:,:) = 0._wp
810	zstrp2(:,:) = 0._wp
811	ENDIF
812	!
813	DO jj = 2, jpjm1
814	DO ji = 2, jpim1
815	IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
816	itframe = 1 ! number of time steps for the running mean
817	IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
818	IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
819	zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
820	zstrp2 (ji,jj) = zstrp1 (ji,jj)
821	zstrp1 (ji,jj) = strength(ji,jj)
822	strength(ji,jj) = zp
823	ENDIF
824	END DO
825	END DO
826	CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
827	!
828	END SELECT
829	!
830	END SUBROUTINE ice_strength
831
832
833	SUBROUTINE ice_dyn_1d2d( kn )
834	!!-----------------------------------------------------------------------
835	!! * ROUTINE ice_dyn_1d2d *
836	!!
837	!! ** Purpose : move arrays from 1d to 2d and the reverse
838	!!-----------------------------------------------------------------------
839	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
840	!
841	INTEGER :: jl, jk ! dummy loop indices
842	!!-----------------------------------------------------------------------
843	!
844	SELECT CASE( kn )
845	! !---------------------!
846	CASE( 1 ) !== from 2D to 1D ==!
847	! !---------------------!
848	! fields used but not modified
849	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
850	CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
851	! the following fields are modified in this routine
852	!!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
853	!!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
854	!!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
855	CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
856	CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
857	CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
858	CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
859	CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
860	CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) )
861	CALL tab_3d_2d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) )
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	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd (:,:) )
877	!
878	! !---------------------!
879	CASE( 2 ) !== from 1D to 2D ==!
880	! !---------------------!
881	CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
882	CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
883	CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
884	CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
885	CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
886	CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
887	CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
888	CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
889	CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) )
890	CALL tab_2d_3d( npti, nptidx(1:npti), sv_ip_2d(1:npti,1:jpl), sv_ip(:,:,:) )
891	DO jl = 1, jpl
892	DO jk = 1, nlay_s
893	CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
894	END DO
895	DO jk = 1, nlay_i
896	CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
897	END DO
898	END DO
899	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
900	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
901	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
902	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) )
903	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
904	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
905	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd (:,:) )
906	!
907	END SELECT
908	!
909	END SUBROUTINE ice_dyn_1d2d
910
911
912	SUBROUTINE ice_dyn_rdgrft_init
913	!!-------------------------------------------------------------------
914	!! * ROUTINE ice_dyn_rdgrft_init *
915	!!
916	!! ** Purpose : Physical constants and parameters linked
917	!! to the mechanical ice redistribution
918	!!
919	!! ** Method : Read the namdyn_rdgrft namelist
920	!! and check the parameters values
921	!! called at the first timestep (nit000)
922	!!
923	!! ** input : Namelist namdyn_rdgrft
924	!!-------------------------------------------------------------------
925	INTEGER :: ios ! Local integer output status for namelist read
926	!!
927	NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
928	& rn_csrdg , &
929	& ln_partf_lin, rn_gstar, &
930	& ln_partf_exp, rn_astar, &
931	& ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, &
932	& ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft
933	!!-------------------------------------------------------------------
934	!
935	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
936	READ ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
937	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
938	REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
939	READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
940	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' )
941	IF(lwm) WRITE ( numoni, namdyn_rdgrft )
942	!
943	IF (lwp) THEN ! control print
944	WRITE(numout,*)
945	WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
946	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
947	WRITE(numout,*) ' Namelist namdyn_rdgrft:'
948	WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79
949	WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar
950	WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg
951	WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg
952	WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
953	WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar
954	WRITE(numout,*) ' Exponential ridging participation function ln_partf_exp = ', ln_partf_exp
955	WRITE(numout,) ' Equivalent to G for an exponential function rn_astar = ', rn_astar
956	WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging
957	WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar
958	WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg
959	WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg
960	WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg
961	WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting
962	WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
963	WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
964	WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft
965	WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft
966	ENDIF
967	!
968	IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
969	CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
970	ENDIF
971	!
972	IF( .NOT. ln_icethd ) THEN
973	rn_porordg = 0._wp
974	rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
975	rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
976	IF( lwp ) THEN
977	WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed'
978	WRITE(numout,*) ' rn_porordg = ', rn_porordg
979	WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg
980	WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg
981	WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft
982	WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft
983	ENDIF
984	ENDIF
985	! ! allocate arrays
986	IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
987	!
988	END SUBROUTINE ice_dyn_rdgrft_init
989
990	#else
991	!!----------------------------------------------------------------------
992	!! Default option Empty module NO SI3 sea-ice model
993	!!----------------------------------------------------------------------
994	#endif
995
996	!!======================================================================
997	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_retain_melt/src/ICE/icedyn_rdgrft.F90

Download in other formats: