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

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

icevar.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icevar.F90 @ 12182

Last change on this file since 12182 was 12182, checked in by davestorkey, 4 years ago
2019/dev_r11943_MERGE_2019: Merge in dev_ASINTER-01-05_merge.
Property svn:keywords set to `Id`
File size: 57.6 KB

Line
1	MODULE icevar
2	!!======================================================================
3	!! * MODULE icevar *
4	!! sea-ice: series of functions to transform or compute ice variables
5	!!======================================================================
6	!! History : - ! 2006-01 (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	!!
14	!! There are three sets of variables
15	!! VGLO : global variables of the model
16	!! - v_i (jpi,jpj,jpl)
17	!! - v_s (jpi,jpj,jpl)
18	!! - a_i (jpi,jpj,jpl)
19	!! - t_s (jpi,jpj,jpl)
20	!! - e_i (jpi,jpj,nlay_i,jpl)
21	!! - e_s (jpi,jpj,nlay_s,jpl)
22	!! - sv_i(jpi,jpj,jpl)
23	!! - oa_i(jpi,jpj,jpl)
24	!! VEQV : equivalent variables sometimes used in the model
25	!! - h_i(jpi,jpj,jpl)
26	!! - h_s(jpi,jpj,jpl)
27	!! - t_i(jpi,jpj,nlay_i,jpl)
28	!! ...
29	!! VAGG : aggregate variables, averaged/summed over all
30	!! thickness categories
31	!! - vt_i(jpi,jpj)
32	!! - vt_s(jpi,jpj)
33	!! - at_i(jpi,jpj)
34	!! - st_i(jpi,jpj)
35	!! - et_s(jpi,jpj) total snow heat content
36	!! - et_i(jpi,jpj) total ice thermal content
37	!! - sm_i(jpi,jpj) mean ice salinity
38	!! - tm_i(jpi,jpj) mean ice temperature
39	!! - tm_s(jpi,jpj) mean snw temperature
40	!!----------------------------------------------------------------------
41	!! ice_var_agg : integrate variables over layers and categories
42	!! ice_var_glo2eqv : transform from VGLO to VEQV
43	!! ice_var_eqv2glo : transform from VEQV to VGLO
44	!! ice_var_salprof : salinity profile in the ice
45	!! ice_var_salprof1d : salinity profile in the ice 1D
46	!! ice_var_zapsmall : remove very small area and volume
47	!! ice_var_zapneg : remove negative ice fields
48	!! ice_var_roundoff : remove negative values arising from roundoff erros
49	!! ice_var_bv : brine volume
50	!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
51	!! ice_var_sshdyn : compute equivalent ssh in lead
52	!! ice_var_itd : convert N-cat to M-cat
53	!!----------------------------------------------------------------------
54	USE dom_oce ! ocean space and time domain
55	USE phycst ! physical constants (ocean directory)
56	USE sbc_oce , ONLY : sss_m, ln_ice_embd, nn_fsbc
57	USE ice ! sea-ice: variables
58	USE ice1D ! sea-ice: thermodynamics variables
59	!
60	USE in_out_manager ! I/O manager
61	USE lib_mpp ! MPP library
62	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
63
64	IMPLICIT NONE
65	PRIVATE
66
67	PUBLIC ice_var_agg
68	PUBLIC ice_var_glo2eqv
69	PUBLIC ice_var_eqv2glo
70	PUBLIC ice_var_salprof
71	PUBLIC ice_var_salprof1d
72	PUBLIC ice_var_zapsmall
73	PUBLIC ice_var_zapneg
74	PUBLIC ice_var_roundoff
75	PUBLIC ice_var_bv
76	PUBLIC ice_var_enthalpy
77	PUBLIC ice_var_sshdyn
78	PUBLIC ice_var_itd
79
80	INTERFACE ice_var_itd
81	MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
82	END INTERFACE
83
84	!!----------------------------------------------------------------------
85	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
86	!! $Id$
87	!! Software governed by the CeCILL license (see ./LICENSE)
88	!!----------------------------------------------------------------------
89	CONTAINS
90
91	SUBROUTINE ice_var_agg( kn )
92	!!-------------------------------------------------------------------
93	!! * ROUTINE ice_var_agg *
94	!!
95	!! ** Purpose : aggregates ice-thickness-category variables to
96	!! all-ice variables, i.e. it turns VGLO into VAGG
97	!!-------------------------------------------------------------------
98	INTEGER, INTENT( in ) :: kn ! =1 state variables only
99	! ! >1 state variables + others
100	!
101	INTEGER :: ji, jj, jk, jl ! dummy loop indices
102	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i, z1_vt_s
103	!!-------------------------------------------------------------------
104	!
105	! ! integrated values
106	vt_i(:,:) = SUM( v_i (:,:,:) , dim=3 )
107	vt_s(:,:) = SUM( v_s (:,:,:) , dim=3 )
108	st_i(:,:) = SUM( sv_i(:,:,:) , dim=3 )
109	at_i(:,:) = SUM( a_i (:,:,:) , dim=3 )
110	et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
111	et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
112	!
113	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
114	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
115	!
116	ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction
117	!
118	!!GS: tm_su always needed by ABL over sea-ice
119	ALLOCATE( z1_at_i(jpi,jpj) )
120	WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:)
121	ELSEWHERE ; z1_at_i(:,:) = 0._wp
122	END WHERE
123	tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
124	WHERE( at_i(:,:)<=epsi20 ) tm_su(:,:) = rt0
125	!
126	! The following fields are calculated for diagnostics and outputs only
127	! ==> Do not use them for other purposes
128	IF( kn > 1 ) THEN
129	!
130	ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
131	WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:)
132	ELSEWHERE ; z1_vt_i(:,:) = 0._wp
133	END WHERE
134	WHERE( vt_s(:,:) > epsi20 ) ; z1_vt_s(:,:) = 1._wp / vt_s(:,:)
135	ELSEWHERE ; z1_vt_s(:,:) = 0._wp
136	END WHERE
137	!
138	! ! mean ice/snow thickness
139	hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
140	hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
141	!
142	! ! mean temperature (K), salinity and age
143	tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
144	om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
145	sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:)
146	!
147	tm_i(:,:) = 0._wp
148	tm_s(:,:) = 0._wp
149	DO jl = 1, jpl
150	DO jk = 1, nlay_i
151	tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
152	END DO
153	DO jk = 1, nlay_s
154	tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
155	END DO
156	END DO
157	!
158	! ! put rt0 where there is no ice
159	WHERE( at_i(:,:)<=epsi20 )
160	tm_si(:,:) = rt0
161	tm_i (:,:) = rt0
162	tm_s (:,:) = rt0
163	END WHERE
164	!
165	! ! mean melt pond depth
166	WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
167	ELSEWHERE ; hm_ip(:,:) = 0._wp
168	END WHERE
169	!
170	DEALLOCATE( z1_vt_i , z1_vt_s )
171	!
172	ENDIF
173	!
174	DEALLOCATE( z1_at_i )
175	!
176	END SUBROUTINE ice_var_agg
177
178
179	SUBROUTINE ice_var_glo2eqv
180	!!-------------------------------------------------------------------
181	!! * ROUTINE ice_var_glo2eqv *
182	!!
183	!! ** Purpose : computes equivalent variables as function of
184	!! global variables, i.e. it turns VGLO into VEQV
185	!!-------------------------------------------------------------------
186	INTEGER :: ji, jj, jk, jl ! dummy loop indices
187	REAL(wp) :: ze_i ! local scalars
188	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
189	REAL(wp) :: zhmax, z1_zhmax ! - -
190	REAL(wp) :: zlay_i, zlay_s ! - -
191	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
192	!!-------------------------------------------------------------------
193
194	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
195	!! ice area and ice volume ?
196	!! the idea is to be able to define one for all at the begining of this routine
197	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
198
199	!---------------------------------------------------------------
200	! Ice thickness, snow thickness, ice salinity, ice age and ponds
201	!---------------------------------------------------------------
202	! !--- inverse of the ice area
203	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
204	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
205	END WHERE
206	!
207	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
208	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
209	END WHERE
210	! !--- ice thickness
211	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
212
213	zhmax = hi_max(jpl)
214	z1_zhmax = 1._wp / hi_max(jpl)
215	WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
216	h_i (:,:,jpl) = zhmax
217	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
218	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
219	END WHERE
220	! !--- snow thickness
221	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
222	! !--- ice age
223	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
224	! !--- pond fraction and thickness
225	a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
226	WHERE( a_ip_frac(:,:,:) > epsi20 ) ; h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
227	ELSEWHERE ; h_ip(:,:,:) = 0._wp
228	END WHERE
229	!
230	! !--- salinity (with a minimum value imposed everywhere)
231	IF( nn_icesal == 2 ) THEN
232	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
233	ELSEWHERE ; s_i(:,:,:) = rn_simin
234	END WHERE
235	ENDIF
236	CALL ice_var_salprof ! salinity profile
237
238	!-------------------
239	! Ice temperature [K] (with a minimum value (rt0 - 100.))
240	!-------------------
241	zlay_i = REAL( nlay_i , wp ) ! number of layers
242	DO jl = 1, jpl
243	DO jk = 1, nlay_i
244	DO jj = 1, jpj
245	DO ji = 1, jpi
246	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
247	!
248	ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
249	ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
250	! Conversion q(S,T) -> T (second order equation)
251	zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
252	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
253	t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
254	!
255	ELSE !--- no ice
256	t_i(ji,jj,jk,jl) = rt0
257	ENDIF
258	END DO
259	END DO
260	END DO
261	END DO
262
263	!--------------------
264	! Snow temperature [K] (with a minimum value (rt0 - 100.))
265	!--------------------
266	zlay_s = REAL( nlay_s , wp )
267	DO jk = 1, nlay_s
268	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
269	t_s(:,:,jk,:) = rt0 + MAX( -100._wp , &
270	& MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
271	ELSEWHERE !--- no ice
272	t_s(:,:,jk,:) = rt0
273	END WHERE
274	END DO
275	!
276	! integrated values
277	vt_i (:,:) = SUM( v_i , dim=3 )
278	vt_s (:,:) = SUM( v_s , dim=3 )
279	at_i (:,:) = SUM( a_i , dim=3 )
280	!
281	END SUBROUTINE ice_var_glo2eqv
282
283
284	SUBROUTINE ice_var_eqv2glo
285	!!-------------------------------------------------------------------
286	!! * ROUTINE ice_var_eqv2glo *
287	!!
288	!! ** Purpose : computes global variables as function of
289	!! equivalent variables, i.e. it turns VEQV into VGLO
290	!!-------------------------------------------------------------------
291	!
292	v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
293	v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
294	sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
295	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
296	!
297	END SUBROUTINE ice_var_eqv2glo
298
299
300	SUBROUTINE ice_var_salprof
301	!!-------------------------------------------------------------------
302	!! * ROUTINE ice_var_salprof *
303	!!
304	!! ** Purpose : computes salinity profile in function of bulk salinity
305	!!
306	!! ** Method : If bulk salinity greater than zsi1,
307	!! the profile is assumed to be constant (S_inf)
308	!! If bulk salinity lower than zsi0,
309	!! the profile is linear with 0 at the surface (S_zero)
310	!! If it is between zsi0 and zsi1, it is a
311	!! alpha-weighted linear combination of s_inf and s_zero
312	!!
313	!! ** References : Vancoppenolle et al., 2007
314	!!-------------------------------------------------------------------
315	INTEGER :: ji, jj, jk, jl ! dummy loop index
316	REAL(wp) :: zsal, z1_dS
317	REAL(wp) :: zargtemp , zs0, zs
318	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
319	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
320	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
321	!!-------------------------------------------------------------------
322
323	!!gm Question: Remove the option 3 ? How many years since it last use ?
324
325	SELECT CASE ( nn_icesal )
326	!
327	! !---------------------------------------!
328	CASE( 1 ) ! constant salinity in time and space !
329	! !---------------------------------------!
330	sz_i(:,:,:,:) = rn_icesal
331	s_i (:,:,:) = rn_icesal
332	!
333	! !---------------------------------------------!
334	CASE( 2 ) ! time varying salinity with linear profile !
335	! !---------------------------------------------!
336	!
337	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
338	!
339	DO jl = 1, jpl
340	DO jk = 1, nlay_i
341	sz_i(:,:,jk,jl) = s_i(:,:,jl)
342	END DO
343	END DO
344	! ! Slope of the linear profile
345	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
346	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
347	END WHERE
348	!
349	z1_dS = 1._wp / ( zsi1 - zsi0 )
350	DO jl = 1, jpl
351	DO jj = 1, jpj
352	DO ji = 1, jpi
353	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
354	! ! force a constant profile when SSS too low (Baltic Sea)
355	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
356	END DO
357	END DO
358	END DO
359	!
360	! Computation of the profile
361	DO jl = 1, jpl
362	DO jk = 1, nlay_i
363	DO jj = 1, jpj
364	DO ji = 1, jpi
365	! ! linear profile with 0 surface value
366	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
367	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
368	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
369	END DO
370	END DO
371	END DO
372	END DO
373	!
374	DEALLOCATE( z_slope_s , zalpha )
375	!
376	! !-------------------------------------------!
377	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
378	! !-------------------------------------------! (mean = 2.30)
379	!
380	s_i(:,:,:) = 2.30_wp
381	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
382	!! a array S_prof(1:nlay_i) containing the calculation and just do:
383	! DO jk = 1, nlay_i
384	! sz_i(:,:,jk,:) = S_prof(jk)
385	! END DO
386	!!gm end
387	!
388	DO jl = 1, jpl
389	DO jk = 1, nlay_i
390	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
391	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
392	END DO
393	END DO
394	!
395	END SELECT
396	!
397	END SUBROUTINE ice_var_salprof
398
399
400	SUBROUTINE ice_var_salprof1d
401	!!-------------------------------------------------------------------
402	!! * ROUTINE ice_var_salprof1d *
403	!!
404	!! ** Purpose : 1d computation of the sea ice salinity profile
405	!! Works with 1d vectors and is used by thermodynamic modules
406	!!-------------------------------------------------------------------
407	INTEGER :: ji, jk ! dummy loop indices
408	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
409	REAL(wp) :: zs, zs0 ! - -
410	!
411	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
412	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
413	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
414	!!-------------------------------------------------------------------
415	!
416	SELECT CASE ( nn_icesal )
417	!
418	! !---------------------------------------!
419	CASE( 1 ) ! constant salinity in time and space !
420	! !---------------------------------------!
421	sz_i_1d(1:npti,:) = rn_icesal
422	!
423	! !---------------------------------------------!
424	CASE( 2 ) ! time varying salinity with linear profile !
425	! !---------------------------------------------!
426	!
427	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
428	!
429	! ! Slope of the linear profile
430	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
431	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
432	END WHERE
433
434	z1_dS = 1._wp / ( zsi1 - zsi0 )
435	DO ji = 1, npti
436	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
437	! ! force a constant profile when SSS too low (Baltic Sea)
438	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
439	END DO
440	!
441	! Computation of the profile
442	DO jk = 1, nlay_i
443	DO ji = 1, npti
444	! ! linear profile with 0 surface value
445	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
446	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
447	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
448	END DO
449	END DO
450	!
451	DEALLOCATE( z_slope_s, zalpha )
452
453	! !-------------------------------------------!
454	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
455	! !-------------------------------------------! (mean = 2.30)
456	!
457	s_i_1d(1:npti) = 2.30_wp
458	!
459	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
460	DO jk = 1, nlay_i
461	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
462	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
463	DO ji = 1, npti
464	sz_i_1d(ji,jk) = zsal
465	END DO
466	END DO
467	!
468	END SELECT
469	!
470	END SUBROUTINE ice_var_salprof1d
471
472
473	SUBROUTINE ice_var_zapsmall
474	!!-------------------------------------------------------------------
475	!! * ROUTINE ice_var_zapsmall *
476	!!
477	!! ** Purpose : Remove too small sea ice areas and correct fluxes
478	!!-------------------------------------------------------------------
479	INTEGER :: ji, jj, jl, jk ! dummy loop indices
480	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
481	!!-------------------------------------------------------------------
482	!
483	DO jl = 1, jpl !== loop over the categories ==!
484	!
485	WHERE( a_i(:,:,jl) > epsi10 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
486	ELSEWHERE ; h_i(:,:,jl) = 0._wp
487	END WHERE
488	!
489	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
490	ELSEWHERE ; zswitch(:,:) = 1._wp
491	END WHERE
492	!
493	!-----------------------------------------------------------------
494	! Zap ice energy and use ocean heat to melt ice
495	!-----------------------------------------------------------------
496	DO jk = 1, nlay_i
497	DO jj = 1 , jpj
498	DO ji = 1 , jpi
499	! update exchanges with ocean
500	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
501	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
502	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
503	END DO
504	END DO
505	END DO
506	!
507	DO jk = 1, nlay_s
508	DO jj = 1 , jpj
509	DO ji = 1 , jpi
510	! update exchanges with ocean
511	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
512	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
513	t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
514	END DO
515	END DO
516	END DO
517	!
518	!-----------------------------------------------------------------
519	! zap ice and snow volume, add water and salt to ocean
520	!-----------------------------------------------------------------
521	DO jj = 1 , jpj
522	DO ji = 1 , jpi
523	! update exchanges with ocean
524	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice
525	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice
526	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice
527	!
528	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
529	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
530	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
531	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
532	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
533	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
534	!
535	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
536	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
537	!
538	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
539	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
540	!
541	END DO
542	END DO
543	!
544	END DO
545
546	! to be sure that at_i is the sum of a_i(jl)
547	at_i (:,:) = SUM( a_i (:,:,:), dim=3 )
548	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 )
549	!!clem add?
550	! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
551	! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
552	! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
553	! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
554	!!clem
555
556	! open water = 1 if at_i=0
557	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
558	!
559	END SUBROUTINE ice_var_zapsmall
560
561
562	SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
563	!!-------------------------------------------------------------------
564	!! * ROUTINE ice_var_zapneg *
565	!!
566	!! ** Purpose : Remove negative sea ice fields and correct fluxes
567	!!-------------------------------------------------------------------
568	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
569	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
570	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
571	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
572	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
573	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
574	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
575	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
576	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
577	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
578	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
579	!
580	INTEGER :: ji, jj, jl, jk ! dummy loop indices
581	REAL(wp) :: z1_dt
582	!!-------------------------------------------------------------------
583	!
584	z1_dt = 1._wp / pdt
585	!
586	DO jl = 1, jpl !== loop over the categories ==!
587	!
588	! make sure a_i=0 where v_i<=0
589	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
590
591	!----------------------------------------
592	! zap ice energy and send it to the ocean
593	!----------------------------------------
594	DO jk = 1, nlay_i
595	DO jj = 1 , jpj
596	DO ji = 1 , jpi
597	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
598	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
599	pe_i(ji,jj,jk,jl) = 0._wp
600	ENDIF
601	END DO
602	END DO
603	END DO
604	!
605	DO jk = 1, nlay_s
606	DO jj = 1 , jpj
607	DO ji = 1 , jpi
608	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
609	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
610	pe_s(ji,jj,jk,jl) = 0._wp
611	ENDIF
612	END DO
613	END DO
614	END DO
615	!
616	!-----------------------------------------------------
617	! zap ice and snow volume, add water and salt to ocean
618	!-----------------------------------------------------
619	DO jj = 1 , jpj
620	DO ji = 1 , jpi
621	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
622	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
623	pv_i (ji,jj,jl) = 0._wp
624	ENDIF
625	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
626	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
627	pv_s (ji,jj,jl) = 0._wp
628	ENDIF
629	IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
630	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
631	psv_i (ji,jj,jl) = 0._wp
632	ENDIF
633	END DO
634	END DO
635	!
636	END DO
637	!
638	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
639	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
640	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
641	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
642	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
643	! but it does not change conservation, so keep it this way is ok
644	!
645	END SUBROUTINE ice_var_zapneg
646
647
648	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
649	!!-------------------------------------------------------------------
650	!! * ROUTINE ice_var_roundoff *
651	!!
652	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
653	!!-------------------------------------------------------------------
654	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
655	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
656	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
657	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
658	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
659	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
660	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
661	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
662	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
663	!!-------------------------------------------------------------------
664	!
665	WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
666	WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
667	WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
668	WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
669	WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
670	WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
671	WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
672	IF( ln_pnd_H12 ) THEN
673	WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
674	WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
675	ENDIF
676	!
677	END SUBROUTINE ice_var_roundoff
678
679
680	SUBROUTINE ice_var_bv
681	!!-------------------------------------------------------------------
682	!! * ROUTINE ice_var_bv *
683	!!
684	!! ** Purpose : computes mean brine volume (%) in sea ice
685	!!
686	!! ** Method : e = - 0.054 * S (ppt) / T (C)
687	!!
688	!! References : Vancoppenolle et al., JGR, 2007
689	!!-------------------------------------------------------------------
690	INTEGER :: ji, jj, jk, jl ! dummy loop indices
691	!!-------------------------------------------------------------------
692	!
693	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
694	!! instead of setting everything to zero as just below
695	bv_i (:,:,:) = 0._wp
696	DO jl = 1, jpl
697	DO jk = 1, nlay_i
698	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
699	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
700	END WHERE
701	END DO
702	END DO
703	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
704	ELSEWHERE ; bvm_i(:,:) = 0._wp
705	END WHERE
706	!
707	END SUBROUTINE ice_var_bv
708
709
710	SUBROUTINE ice_var_enthalpy
711	!!-------------------------------------------------------------------
712	!! * ROUTINE ice_var_enthalpy *
713	!!
714	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
715	!!
716	!! ** Method : Formula (Bitz and Lipscomb, 1999)
717	!!-------------------------------------------------------------------
718	INTEGER :: ji, jk ! dummy loop indices
719	REAL(wp) :: ztmelts ! local scalar
720	!!-------------------------------------------------------------------
721	!
722	DO jk = 1, nlay_i ! Sea ice energy of melting
723	DO ji = 1, npti
724	ztmelts = - rTmlt * sz_i_1d(ji,jk)
725	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
726	! (sometimes zdf scheme produces abnormally high temperatures)
727	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
728	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
729	& - rcp * ztmelts )
730	END DO
731	END DO
732	DO jk = 1, nlay_s ! Snow energy of melting
733	DO ji = 1, npti
734	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
735	END DO
736	END DO
737	!
738	END SUBROUTINE ice_var_enthalpy
739
740
741	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
742	!!---------------------------------------------------------------------
743	!! * ROUTINE ice_var_sshdyn *
744	!!
745	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
746	!!
747	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
748	!!
749	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
750	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
751	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
752	!!----------------------------------------------------------------------
753	!
754	! input
755	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
756	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
757	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
758	!
759	! output
760	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
761	!
762	! temporary
763	REAL(wp) :: zintn, zintb ! time interpolation weights []
764	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
765	!
766	! compute ice load used to define the equivalent ssh in lead
767	IF( ln_ice_embd ) THEN
768	!
769	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
770	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
771	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
772	!
773	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
774	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
775	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
776	!
777	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
778	!
779	ELSE
780	zsnwiceload(:,:) = 0.0_wp
781	ENDIF
782	! compute equivalent ssh in lead
783	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
784	!
785	END FUNCTION ice_var_sshdyn
786
787
788	!!-------------------------------------------------------------------
789	!! * INTERFACE ice_var_itd *
790	!!
791	!! ** Purpose : converting N-cat ice to jpl ice categories
792	!!-------------------------------------------------------------------
793	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
794	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
795	!!-------------------------------------------------------------------
796	!! ** Purpose : converting 1-cat ice to 1 ice category
797	!!-------------------------------------------------------------------
798	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
799	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
800	REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
801	REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
802	!!-------------------------------------------------------------------
803	! == thickness and concentration == !
804	ph_i(:) = phti(:)
805	ph_s(:) = phts(:)
806	pa_i(:) = pati(:)
807	!
808	! == temperature and salinity and ponds == !
809	pt_i (:) = ptmi (:)
810	pt_s (:) = ptms (:)
811	pt_su(:) = ptmsu(:)
812	ps_i (:) = psmi (:)
813	pa_ip(:) = patip(:)
814	ph_ip(:) = phtip(:)
815
816	END SUBROUTINE ice_var_itd_1c1c
817
818	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
819	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
820	!!-------------------------------------------------------------------
821	!! ** Purpose : converting N-cat ice to 1 ice category
822	!!-------------------------------------------------------------------
823	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
824	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
825	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
826	REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
827	!
828	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
829	!
830	INTEGER :: idim
831	!!-------------------------------------------------------------------
832	!
833	idim = SIZE( phti, 1 )
834	!
835	! == thickness and concentration == !
836	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) )
837	!
838	pa_i(:) = SUM( pati(:,:), dim=2 )
839
840	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
841	ELSEWHERE ; z1_ai(:) = 0._wp
842	END WHERE
843
844	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
845	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
846	!
847	! == temperature and salinity == !
848	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
849	ELSEWHERE ; z1_vi(:) = 0._wp
850	END WHERE
851	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
852	ELSEWHERE ; z1_vs(:) = 0._wp
853	END WHERE
854	pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
855	pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
856	pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
857	ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
858
859	! == ponds == !
860	pa_ip(:) = SUM( patip(:,:), dim=2 )
861	WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
862	ELSEWHERE ; ph_ip(:) = 0._wp
863	END WHERE
864	!
865	DEALLOCATE( z1_ai, z1_vi, z1_vs )
866	!
867	END SUBROUTINE ice_var_itd_Nc1c
868
869	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
870	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
871	!!-------------------------------------------------------------------
872	!!
873	!! ** Purpose : converting 1-cat ice to jpl ice categories
874	!!
875	!!
876	!! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015)
877	!! it has the property of conserving total concentration and volume
878	!!
879	!!
880	!! ** Arguments : phti: 1-cat ice thickness
881	!! phts: 1-cat snow depth
882	!! pati: 1-cat ice concentration
883	!!
884	!! ** Output : jpl-cat
885	!!
886	!! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
887	!! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
888	!! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614
889	!!-------------------------------------------------------------------
890	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
891	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
892	REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
893	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
894	!
895	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti
896	INTEGER :: ji, jk, jl
897	INTEGER :: idim
898	REAL(wp) :: zv, zdh
899	!!-------------------------------------------------------------------
900	!
901	idim = SIZE( phti , 1 )
902	!
903	ph_i(1:idim,1:jpl) = 0._wp
904	ph_s(1:idim,1:jpl) = 0._wp
905	pa_i(1:idim,1:jpl) = 0._wp
906	!
907	ALLOCATE( z1_hti(idim) )
908	WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:)
909	ELSEWHERE ; z1_hti(:) = 0._wp
910	END WHERE
911	!
912	! == thickness and concentration == !
913	! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
914	DO jl = 1, jpl-1
915	DO ji = 1, idim
916	!
917	IF( phti(ji) > 0._wp ) THEN
918	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
919	pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.hi_max(jl-1) ) EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
920	& - ( phti(ji) + 2.hi_max(jl ) ) EXP( -2.hi_max(jl )z1_hti(ji) ) )
921	!
922	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
923	zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl-1) + 2.hi_max(jl-1)*hi_max(jl-1) ) &
924	& * EXP( -2.hi_max(jl-1)z1_hti(ji) ) &
925	& - ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jl) + 2.hi_max(jl)*hi_max(jl) ) &
926	& * EXP(-2.hi_max(jl)z1_hti(ji)) )
927	! thickness
928	IF( pa_i(ji,jl) > epsi06 ) THEN
929	ph_i(ji,jl) = zv / pa_i(ji,jl)
930	ELSE
931	ph_i(ji,jl) = 0.
932	pa_i(ji,jl) = 0.
933	ENDIF
934	ENDIF
935	!
936	ENDDO
937	ENDDO
938	!
939	! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
940	DO ji = 1, idim
941	!
942	IF( phti(ji) > 0._wp ) THEN
943	! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
944	pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.hi_max(jpl-1) ) EXP( -2.hi_max(jpl-1)z1_hti(ji) )
945
946	! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
947	zv = pati(ji) * z1_hti(ji) * ( phti(ji)phti(ji) + 2.phti(ji)hi_max(jpl-1) + 2.hi_max(jpl-1)*hi_max(jpl-1) ) &
948	& * EXP( -2.hi_max(jpl-1)z1_hti(ji) )
949	! thickness
950	IF( pa_i(ji,jpl) > epsi06 ) THEN
951	ph_i(ji,jpl) = zv / pa_i(ji,jpl)
952	else
953	ph_i(ji,jpl) = 0.
954	pa_i(ji,jpl) = 0.
955	ENDIF
956	ENDIF
957	!
958	ENDDO
959	!
960	! Add Snow in each category where pa_i is not 0
961	DO jl = 1, jpl
962	DO ji = 1, idim
963	IF( pa_i(ji,jl) > 0._wp ) THEN
964	ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
965	! In case snow load is in excess that would lead to transformation from snow to ice
966	! Then, transfer the snow excess into the ice (different from icethd_dh)
967	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )
968	! recompute h_i, h_s avoiding out of bounds values
969	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
970	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
971	ENDIF
972	END DO
973	END DO
974	!
975	DEALLOCATE( z1_hti )
976	!
977	! == temperature and salinity == !
978	DO jl = 1, jpl
979	pt_i (:,jl) = ptmi (:)
980	pt_s (:,jl) = ptms (:)
981	pt_su(:,jl) = ptmsu(:)
982	ps_i (:,jl) = psmi (:)
983	ps_i (:,jl) = psmi (:)
984	END DO
985	!
986	! == ponds == !
987	ALLOCATE( zfra(idim) )
988	! keep the same pond fraction atip/ati for each category
989	WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:)
990	ELSEWHERE ; zfra(:) = 0._wp
991	END WHERE
992	DO jl = 1, jpl
993	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
994	END DO
995	! keep the same v_ip/v_i ratio for each category
996	WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
997	ELSEWHERE ; zfra(:) = 0._wp
998	END WHERE
999	DO jl = 1, jpl
1000	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1001	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1002	END WHERE
1003	END DO
1004	DEALLOCATE( zfra )
1005	!
1006	END SUBROUTINE ice_var_itd_1cMc
1007
1008	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
1009	& ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
1010	!!-------------------------------------------------------------------
1011	!!
1012	!! ** Purpose : converting N-cat ice to jpl ice categories
1013	!!
1014	!! ice thickness distribution follows a gaussian law
1015	!! around the concentration of the most likely ice thickness
1016	!! (similar as iceistate.F90)
1017	!!
1018	!! ** Method: Iterative procedure
1019	!!
1020	!! 1) Fill ice cat that correspond to input thicknesses
1021	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
1022	!!
1023	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
1024	!! by removing 25% ice area from jlmin and jlmax (resp.)
1025	!!
1026	!! 3) Expand the filling to the empty cat between jlmin and jlmax
1027	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
1028	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
1029	!!
1030	!! ** Arguments : phti: N-cat ice thickness
1031	!! phts: N-cat snow depth
1032	!! pati: N-cat ice concentration
1033	!!
1034	!! ** Output : jpl-cat
1035	!!
1036	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
1037	!!-------------------------------------------------------------------
1038	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
1039	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
1040	REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds
1041	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds
1042	!
1043	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
1044	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
1045	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra
1046	!
1047	REAL(wp), PARAMETER :: ztrans = 0.25_wp
1048	INTEGER :: ji, jl, jl1, jl2
1049	INTEGER :: idim, icat
1050	!!-------------------------------------------------------------------
1051	!
1052	idim = SIZE( phti, 1 )
1053	icat = SIZE( phti, 2 )
1054	!
1055	! == thickness and concentration == !
1056	! ! ---------------------- !
1057	IF( icat == jpl ) THEN ! input cat = output cat !
1058	! ! ---------------------- !
1059	ph_i(:,:) = phti(:,:)
1060	ph_s(:,:) = phts(:,:)
1061	pa_i(:,:) = pati(:,:)
1062	!
1063	! == temperature and salinity and ponds == !
1064	pt_i (:,:) = ptmi (:,:)
1065	pt_s (:,:) = ptms (:,:)
1066	pt_su(:,:) = ptmsu(:,:)
1067	ps_i (:,:) = psmi (:,:)
1068	pa_ip(:,:) = patip(:,:)
1069	ph_ip(:,:) = phtip(:,:)
1070	! ! ---------------------- !
1071	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
1072	! ! ---------------------- !
1073	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
1074	& ph_i(:,:), ph_s(:,:), pa_i (:,:), &
1075	& ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
1076	& pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) )
1077	! ! ---------------------- !
1078	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
1079	! ! ---------------------- !
1080	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
1081	& ph_i(:,1), ph_s(:,1), pa_i (:,1), &
1082	& ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
1083	& pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) )
1084	! ! ----------------------- !
1085	ELSE ! input cat /= output cat !
1086	! ! ----------------------- !
1087
1088	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
1089	ALLOCATE( jlmin(idim), jlmax(idim) )
1090
1091	! --- initialize output fields to 0 --- !
1092	ph_i(1:idim,1:jpl) = 0._wp
1093	ph_s(1:idim,1:jpl) = 0._wp
1094	pa_i(1:idim,1:jpl) = 0._wp
1095	!
1096	! --- fill the categories --- !
1097	! find where cat-input = cat-output and fill cat-output fields
1098	jlmax(:) = 0
1099	jlmin(:) = 999
1100	jlfil(:,:) = 0
1101	DO jl1 = 1, jpl
1102	DO jl2 = 1, icat
1103	DO ji = 1, idim
1104	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
1105	! fill the right category
1106	ph_i(ji,jl1) = phti(ji,jl2)
1107	ph_s(ji,jl1) = phts(ji,jl2)
1108	pa_i(ji,jl1) = pati(ji,jl2)
1109	! record categories that are filled
1110	jlmax(ji) = MAX( jlmax(ji), jl1 )
1111	jlmin(ji) = MIN( jlmin(ji), jl1 )
1112	jlfil(ji,jl1) = jl1
1113	ENDIF
1114	END DO
1115	END DO
1116	END DO
1117	!
1118	! --- fill the gaps between categories --- !
1119	! transfer from categories filled at the previous step to the empty ones in between
1120	DO ji = 1, idim
1121	jl1 = jlmin(ji)
1122	jl2 = jlmax(ji)
1123	IF( jl1 > 1 ) THEN
1124	! fill the lower cat (jl1-1)
1125	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
1126	ph_i(ji,jl1-1) = hi_mean(jl1-1)
1127	! remove from cat jl1
1128	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
1129	ENDIF
1130	IF( jl2 < jpl ) THEN
1131	! fill the upper cat (jl2+1)
1132	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
1133	ph_i(ji,jl2+1) = hi_mean(jl2+1)
1134	! remove from cat jl2
1135	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
1136	ENDIF
1137	END DO
1138	!
1139	jlfil2(:,:) = jlfil(:,:)
1140	! fill categories from low to high
1141	DO jl = 2, jpl-1
1142	DO ji = 1, idim
1143	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
1144	! fill high
1145	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
1146	ph_i(ji,jl) = hi_mean(jl)
1147	jlfil(ji,jl) = jl
1148	! remove low
1149	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
1150	ENDIF
1151	END DO
1152	END DO
1153	!
1154	! fill categories from high to low
1155	DO jl = jpl-1, 2, -1
1156	DO ji = 1, idim
1157	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
1158	! fill low
1159	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
1160	ph_i(ji,jl) = hi_mean(jl)
1161	jlfil2(ji,jl) = jl
1162	! remove high
1163	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
1164	ENDIF
1165	END DO
1166	END DO
1167	!
1168	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
1169	DEALLOCATE( jlmin, jlmax )
1170	!
1171	! == temperature and salinity == !
1172	!
1173	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
1174	!
1175	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
1176	ELSEWHERE ; z1_ai(:) = 0._wp
1177	END WHERE
1178	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
1179	ELSEWHERE ; z1_vi(:) = 0._wp
1180	END WHERE
1181	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
1182	ELSEWHERE ; z1_vs(:) = 0._wp
1183	END WHERE
1184	!
1185	! fill all the categories with the same value
1186	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1187	DO jl = 1, jpl
1188	pt_i (:,jl) = ztmp(:)
1189	END DO
1190	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
1191	DO jl = 1, jpl
1192	pt_s (:,jl) = ztmp(:)
1193	END DO
1194	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
1195	DO jl = 1, jpl
1196	pt_su(:,jl) = ztmp(:)
1197	END DO
1198	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
1199	DO jl = 1, jpl
1200	ps_i (:,jl) = ztmp(:)
1201	END DO
1202	!
1203	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
1204	!
1205	! == ponds == !
1206	ALLOCATE( zfra(idim) )
1207	! keep the same pond fraction atip/ati for each category
1208	WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
1209	ELSEWHERE ; zfra(:) = 0._wp
1210	END WHERE
1211	DO jl = 1, jpl
1212	pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
1213	END DO
1214	! keep the same v_ip/v_i ratio for each category
1215	WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
1216	zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
1217	ELSEWHERE
1218	zfra(:) = 0._wp
1219	END WHERE
1220	DO jl = 1, jpl
1221	WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
1222	ELSEWHERE ; ph_ip(:,jl) = 0._wp
1223	END WHERE
1224	END DO
1225	DEALLOCATE( zfra )
1226	!
1227	ENDIF
1228	!
1229	END SUBROUTINE ice_var_itd_NcMc
1230
1231	#else
1232	!!----------------------------------------------------------------------
1233	!! Default option Dummy module NO SI3 sea-ice model
1234	!!----------------------------------------------------------------------
1235	#endif
1236
1237	!!======================================================================
1238	END MODULE icevar

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

Download in other formats: