Context Navigation

icedyn_rdgrft.F90 @ 15653

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

Download in other formats: