Context Navigation

icevar.F90 @ 8984

Last change on this file since 8984 was 8984, checked in by clem, 6 years ago
clarify vertical heat diffusion. There are still issues with heat conservation in case nice_jules=2 but it will be solved in due time. Sette tests are all OK except ISOMIP
File size: 38.9 KB

Line
1	MODULE icevar
2	!!======================================================================
3	!! * MODULE icevar *
4	!! sea-ice: Different sets of ice model variables
5	!! how to switch from one to another
6	!!
7	!! There are three sets of variables
8	!! VGLO : global variables of the model
9	!! - v_i (jpi,jpj,jpl)
10	!! - v_s (jpi,jpj,jpl)
11	!! - a_i (jpi,jpj,jpl)
12	!! - t_s (jpi,jpj,jpl)
13	!! - e_i (jpi,jpj,nlay_i,jpl)
14	!! - sv_i(jpi,jpj,jpl)
15	!! - oa_i(jpi,jpj,jpl)
16	!! VEQV : equivalent variables sometimes used in the model
17	!! - h_i(jpi,jpj,jpl)
18	!! - h_s(jpi,jpj,jpl)
19	!! - t_i(jpi,jpj,nlay_i,jpl)
20	!! ...
21	!! VAGG : aggregate variables, averaged/summed over all
22	!! thickness categories
23	!! - vt_i(jpi,jpj)
24	!! - vt_s(jpi,jpj)
25	!! - at_i(jpi,jpj)
26	!! - et_s(jpi,jpj) !total snow heat content
27	!! - et_i(jpi,jpj) !total ice thermal content
28	!! - sm_i(jpi,jpj) !mean ice salinity
29	!! - tm_i (jpi,jpj) !mean ice temperature
30	!!======================================================================
31	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
32	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
33	!! 3.5 ! 2012 (M. Vancoppenolle) add ice_var_itd
34	!! 3.6 ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
35	!!----------------------------------------------------------------------
36	#if defined key_lim3
37	!!----------------------------------------------------------------------
38	!! 'key_lim3' ESIM sea-ice model
39	!!----------------------------------------------------------------------
40	!! ice_var_agg : integrate variables over layers and categories
41	!! ice_var_glo2eqv : transform from VGLO to VEQV
42	!! ice_var_eqv2glo : transform from VEQV to VGLO
43	!! ice_var_salprof : salinity profile in the ice
44	!! ice_var_salprof1d : salinity profile in the ice 1D
45	!! ice_var_zapsmall : remove very small area and volume
46	!! ice_var_itd : convert 1-cat to jpl-cat
47	!! ice_var_itd2 : convert N-cat to jpl-cat
48	!! ice_var_bv : brine volume
49	!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
50	!!----------------------------------------------------------------------
51	USE dom_oce ! ocean space and time domain
52	USE phycst ! physical constants (ocean directory)
53	USE sbc_oce , ONLY : sss_m
54	USE ice ! sea-ice: variables
55	USE ice1D ! sea-ice: thermodynamics variables
56	!
57	USE in_out_manager ! I/O manager
58	USE lib_mpp ! MPP library
59	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
60
61	IMPLICIT NONE
62	PRIVATE
63
64	PUBLIC ice_var_agg
65	PUBLIC ice_var_glo2eqv
66	PUBLIC ice_var_eqv2glo
67	PUBLIC ice_var_salprof
68	PUBLIC ice_var_salprof1d
69	PUBLIC ice_var_zapsmall
70	PUBLIC ice_var_itd
71	PUBLIC ice_var_itd2
72	PUBLIC ice_var_bv
73	PUBLIC ice_var_enthalpy
74
75	!!----------------------------------------------------------------------
76	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
77	!! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
78	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
79	!!----------------------------------------------------------------------
80	CONTAINS
81
82	SUBROUTINE ice_var_agg( kn )
83	!!-------------------------------------------------------------------
84	!! * ROUTINE ice_var_agg *
85	!!
86	!! ** Purpose : aggregates ice-thickness-category variables to
87	!! all-ice variables, i.e. it turns VGLO into VAGG
88	!!-------------------------------------------------------------------
89	INTEGER, INTENT( in ) :: kn ! =1 state variables only
90	! ! >1 state variables + others
91	!
92	INTEGER :: ji, jj, jk, jl ! dummy loop indices
93	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i
94	!!-------------------------------------------------------------------
95	!
96	! ! integrated values
97	vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 )
98	vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 )
99	at_i(:,:) = SUM( a_i(:,:,:) , dim=3 )
100	et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
101	et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
102	!
103	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
104	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
105	!
106	ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction
107
108	IF( kn > 1 ) THEN
109	!
110	ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
111	WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:)
112	ELSEWHERE ; z1_at_i(:,:) = 0._wp
113	END WHERE
114	WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:)
115	ELSEWHERE ; z1_vt_i(:,:) = 0._wp
116	END WHERE
117	!
118	! ! mean ice/snow thickness
119	hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
120	hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
121	!
122	! ! mean temperature (K), salinity and age
123	tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
124	tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
125	om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
126	!
127	tm_i(:,:) = 0._wp
128	sm_i(:,:) = 0._wp
129	DO jl = 1, jpl
130	DO jk = 1, nlay_i
131	tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
132	sm_i(:,:) = sm_i(:,:) + r1_nlay_i * sz_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
133	END DO
134	END DO
135	!
136	! ! put rt0 where there is no ice
137	WHERE( at_i(:,:)<=epsi20 )
138	tm_su(:,:) = rt0
139	tm_si(:,:) = rt0
140	tm_i (:,:) = rt0
141	END WHERE
142
143	DEALLOCATE( z1_at_i , z1_vt_i )
144	ENDIF
145	!
146	END SUBROUTINE ice_var_agg
147
148
149	SUBROUTINE ice_var_glo2eqv
150	!!-------------------------------------------------------------------
151	!! * ROUTINE ice_var_glo2eqv *
152	!!
153	!! ** Purpose : computes equivalent variables as function of
154	!! global variables, i.e. it turns VGLO into VEQV
155	!!-------------------------------------------------------------------
156	INTEGER :: ji, jj, jk, jl ! dummy loop indices
157	REAL(wp) :: ze_i ! local scalars
158	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
159	REAL(wp) :: zhmax, z1_zhmax ! - -
160	REAL(wp) :: zlay_i, zlay_s ! - -
161	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
162	!!-------------------------------------------------------------------
163
164	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
165	!! ice area and ice volume ?
166	!! the idea is to be able to define one for all at the begining of this routine
167	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
168
169	!---------------------------------------------------------------
170	! Ice thickness, snow thickness, ice salinity, ice age and ponds
171	!---------------------------------------------------------------
172	! !--- inverse of the ice area
173	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
174	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
175	END WHERE
176	!
177	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
178	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
179	END WHERE
180	! !--- ice thickness
181	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
182
183	zhmax = hi_max(jpl)
184	z1_zhmax = 1._wp / hi_max(jpl)
185	WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
186	h_i (:,:,jpl) = zhmax
187	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
188	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
189	END WHERE
190	! !--- snow thickness
191	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
192	! !--- ice age
193	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
194	! !--- pond fraction and thickness
195	a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
196	WHERE( a_ip_frac(:,:,:) > epsi20 ) ; h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
197	ELSEWHERE ; h_ip(:,:,:) = 0._wp
198	END WHERE
199	!
200	! !--- salinity (with a minimum value imposed everywhere)
201	IF( nn_icesal == 2 ) THEN
202	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
203	ELSEWHERE ; s_i(:,:,:) = rn_simin
204	END WHERE
205	ENDIF
206	CALL ice_var_salprof ! salinity profile
207
208	!-------------------
209	! Ice temperature [K] (with a minimum value (rt0 - 100.))
210	!-------------------
211	zlay_i = REAL( nlay_i , wp ) ! number of layers
212	DO jl = 1, jpl
213	DO jk = 1, nlay_i
214	DO jj = 1, jpj
215	DO ji = 1, jpi
216	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
217	!
218	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]
219	ztmelts = - sz_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C]
220	! Conversion q(S,T) -> T (second order equation)
221	zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
222	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
223	t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
224	!
225	ELSE !--- no ice
226	t_i(ji,jj,jk,jl) = rt0
227	ENDIF
228	END DO
229	END DO
230	END DO
231	END DO
232
233	!--------------------
234	! Snow temperature [K] (with a minimum value (rt0 - 100.))
235	!--------------------
236	zlay_s = REAL( nlay_s , wp )
237	DO jk = 1, nlay_s
238	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
239	t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
240	ELSEWHERE !--- no ice
241	t_s(:,:,jk,:) = rt0
242	END WHERE
243	END DO
244	!
245	! integrated values
246	vt_i (:,:) = SUM( v_i, dim=3 )
247	vt_s (:,:) = SUM( v_s, dim=3 )
248	at_i (:,:) = SUM( a_i, dim=3 )
249	!
250	END SUBROUTINE ice_var_glo2eqv
251
252
253	SUBROUTINE ice_var_eqv2glo
254	!!-------------------------------------------------------------------
255	!! * ROUTINE ice_var_eqv2glo *
256	!!
257	!! ** Purpose : computes global variables as function of
258	!! equivalent variables, i.e. it turns VEQV into VGLO
259	!!-------------------------------------------------------------------
260	!
261	v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
262	v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
263	sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
264	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
265	!
266	END SUBROUTINE ice_var_eqv2glo
267
268
269	SUBROUTINE ice_var_salprof
270	!!-------------------------------------------------------------------
271	!! * ROUTINE ice_var_salprof *
272	!!
273	!! ** Purpose : computes salinity profile in function of bulk salinity
274	!!
275	!! ** Method : If bulk salinity greater than zsi1,
276	!! the profile is assumed to be constant (S_inf)
277	!! If bulk salinity lower than zsi0,
278	!! the profile is linear with 0 at the surface (S_zero)
279	!! If it is between zsi0 and zsi1, it is a
280	!! alpha-weighted linear combination of s_inf and s_zero
281	!!
282	!! ** References : Vancoppenolle et al., 2007
283	!!-------------------------------------------------------------------
284	INTEGER :: ji, jj, jk, jl ! dummy loop index
285	REAL(wp) :: zsal, z1_dS
286	REAL(wp) :: zargtemp , zs0, zs
287	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
288	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
289	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
290	!!-------------------------------------------------------------------
291
292	!!gm Question: Remove the option 3 ? How many years since it last use ?
293
294	SELECT CASE ( nn_icesal )
295	!
296	! !---------------------------------------!
297	CASE( 1 ) ! constant salinity in time and space !
298	! !---------------------------------------!
299	sz_i(:,:,:,:) = rn_icesal
300	s_i(:,:,:) = rn_icesal
301	!
302	! !---------------------------------------------!
303	CASE( 2 ) ! time varying salinity with linear profile !
304	! !---------------------------------------------!
305	!
306	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
307	!
308	DO jl = 1, jpl
309	DO jk = 1, nlay_i
310	sz_i(:,:,jk,jl) = s_i(:,:,jl)
311	END DO
312	END DO
313	! ! Slope of the linear profile
314	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
315	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
316	END WHERE
317	!
318	z1_dS = 1._wp / ( zsi1 - zsi0 )
319	DO jl = 1, jpl
320	DO jj = 1, jpj
321	DO ji = 1, jpi
322	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
323	! ! force a constant profile when SSS too low (Baltic Sea)
324	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
325	END DO
326	END DO
327	END DO
328	!
329	! Computation of the profile
330	DO jl = 1, jpl
331	DO jk = 1, nlay_i
332	DO jj = 1, jpj
333	DO ji = 1, jpi
334	! ! linear profile with 0 surface value
335	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
336	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
337	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
338	END DO
339	END DO
340	END DO
341	END DO
342	!
343	DEALLOCATE( z_slope_s , zalpha )
344	!
345	! !-------------------------------------------!
346	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
347	! !-------------------------------------------! (mean = 2.30)
348	!
349	s_i(:,:,:) = 2.30_wp
350	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
351	!! a array S_prof(1:nlay_i) containing the calculation and just do:
352	! DO jk = 1, nlay_i
353	! sz_i(:,:,jk,:) = S_prof(jk)
354	! END DO
355	!!gm end
356	!
357	DO jl = 1, jpl
358	DO jk = 1, nlay_i
359	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
360	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
361	END DO
362	END DO
363	!
364	END SELECT
365	!
366	END SUBROUTINE ice_var_salprof
367
368
369	SUBROUTINE ice_var_salprof1d
370	!!-------------------------------------------------------------------
371	!! * ROUTINE ice_var_salprof1d *
372	!!
373	!! ** Purpose : 1d computation of the sea ice salinity profile
374	!! Works with 1d vectors and is used by thermodynamic modules
375	!!-------------------------------------------------------------------
376	INTEGER :: ji, jk ! dummy loop indices
377	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
378	REAL(wp) :: zs, zs0 ! - -
379	!
380	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
381	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
382	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
383	!!-------------------------------------------------------------------
384	!
385	SELECT CASE ( nn_icesal )
386	!
387	! !---------------------------------------!
388	CASE( 1 ) ! constant salinity in time and space !
389	! !---------------------------------------!
390	sz_i_1d(1:npti,:) = rn_icesal
391	!
392	! !---------------------------------------------!
393	CASE( 2 ) ! time varying salinity with linear profile !
394	! !---------------------------------------------!
395	!
396	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
397	!
398	! ! Slope of the linear profile
399	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
400	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
401	END WHERE
402
403	z1_dS = 1._wp / ( zsi1 - zsi0 )
404	DO ji = 1, npti
405	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
406	! ! force a constant profile when SSS too low (Baltic Sea)
407	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
408	END DO
409	!
410	! Computation of the profile
411	DO jk = 1, nlay_i
412	DO ji = 1, npti
413	! ! linear profile with 0 surface value
414	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
415	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
416	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
417	END DO
418	END DO
419	!
420	DEALLOCATE( z_slope_s, zalpha )
421
422	! !-------------------------------------------!
423	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
424	! !-------------------------------------------! (mean = 2.30)
425	!
426	s_i_1d(1:npti) = 2.30_wp
427	!
428	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
429	DO jk = 1, nlay_i
430	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
431	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
432	DO ji = 1, npti
433	sz_i_1d(ji,jk) = zsal
434	END DO
435	END DO
436	!
437	END SELECT
438	!
439	END SUBROUTINE ice_var_salprof1d
440
441
442	SUBROUTINE ice_var_zapsmall
443	!!-------------------------------------------------------------------
444	!! * ROUTINE ice_var_zapsmall *
445	!!
446	!! ** Purpose : Remove too small sea ice areas and correct fluxes
447	!!-------------------------------------------------------------------
448	INTEGER :: ji, jj, jl, jk ! dummy loop indices
449	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
450	!!-------------------------------------------------------------------
451	!
452	DO jl = 1, jpl !== loop over the categories ==!
453	!
454	!-----------------------------------------------------------------
455	! Zap ice energy and use ocean heat to melt ice
456	!-----------------------------------------------------------------
457	WHERE( a_i(:,:,jl) > epsi20 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
458	ELSEWHERE ; h_i(:,:,jl) = 0._wp
459	END WHERE
460	!
461	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
462	ELSEWHERE ; zswitch(:,:) = 1._wp
463	END WHERE
464	!
465	DO jk = 1, nlay_i
466	DO jj = 1 , jpj
467	DO ji = 1 , jpi
468	! update exchanges with ocean
469	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
470	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
471	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
472	END DO
473	END DO
474	END DO
475	!
476	DO jj = 1 , jpj
477	DO ji = 1 , jpi
478	! update exchanges with ocean
479	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_rdtice
480	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
481	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
482	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
483	!-----------------------------------------------------------------
484	! Zap snow energy
485	!-----------------------------------------------------------------
486	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
487	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
488	!
489	!-----------------------------------------------------------------
490	! zap ice and snow volume, add water and salt to ocean
491	!-----------------------------------------------------------------
492	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
493	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
494	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
495	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
496	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
497	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
498	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
499	!
500	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
501	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
502	!
503	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
504	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
505	!
506	END DO
507	END DO
508	!
509	END DO
510
511	! to be sure that at_i is the sum of a_i(jl)
512	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
513	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
514
515	! open water = 1 if at_i=0
516	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
517	!
518	END SUBROUTINE ice_var_zapsmall
519
520
521	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
522	!!-------------------------------------------------------------------
523	!! * ROUTINE ice_var_itd *
524	!!
525	!! ** Purpose : converting 1-cat ice to multiple ice categories
526	!!
527	!! ice thickness distribution follows a gaussian law
528	!! around the concentration of the most likely ice thickness
529	!! (similar as iceistate.F90)
530	!!
531	!! ** Method: Iterative procedure
532	!!
533	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
534	!!
535	!! 2) Check whether the distribution conserves area and volume, positivity and
536	!! category boundaries
537	!!
538	!! 3) If not (input ice is too thin), the last category is empty and
539	!! the number of categories is reduced (jpl-1)
540	!!
541	!! 4) Iterate until ok (SUM(itest(:) = 4)
542	!!
543	!! ** Arguments : zhti: 1-cat ice thickness
544	!! zhts: 1-cat snow depth
545	!! zati: 1-cat ice concentration
546	!!
547	!! ** Output : jpl-cat
548	!!
549	!! (Example of application: BDY forcings when input are cell averaged)
550	!!-------------------------------------------------------------------
551	INTEGER :: ji, jk, jl ! dummy loop indices
552	INTEGER :: idim, i_fill, jl0
553	REAL(wp) :: zarg, zV, zconv, zdh, zdv
554	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
555	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
556	INTEGER , DIMENSION(4) :: itest
557	!!-------------------------------------------------------------------
558	!
559	! ----------------------------------------
560	! distribution over the jpl ice categories
561	! ----------------------------------------
562	! a gaussian distribution for ice concentration is used
563	! then we check whether the distribution fullfills
564	! volume and area conservation, positivity and ice categories bounds
565	idim = SIZE( zhti , 1 )
566	zh_i(1:idim,1:jpl) = 0._wp
567	zh_s(1:idim,1:jpl) = 0._wp
568	za_i(1:idim,1:jpl) = 0._wp
569	!
570	DO ji = 1, idim
571	!
572	IF( zhti(ji) > 0._wp ) THEN
573	!
574	! find which category (jl0) the input ice thickness falls into
575	jl0 = jpl
576	DO jl = 1, jpl
577	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
578	jl0 = jl
579	CYCLE
580	ENDIF
581	END DO
582	!
583	itest(:) = 0
584	i_fill = jpl + 1 !------------------------------------
585	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
586	! !------------------------------------
587	i_fill = i_fill - 1
588	!
589	zh_i(ji,1:jpl) = 0._wp
590	za_i(ji,1:jpl) = 0._wp
591	itest(:) = 0
592	!
593	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
594	zh_i(ji,1) = zhti(ji)
595	za_i (ji,1) = zati (ji)
596	ELSE !-- case ice is thicker: fill categories >1
597	! thickness
598	DO jl = 1, i_fill - 1
599	zh_i(ji,jl) = hi_mean(jl)
600	END DO
601	!
602	! concentration
603	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
604	DO jl = 1, i_fill - 1
605	IF ( jl /= jl0 ) THEN
606	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
607	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
608	ENDIF
609	END DO
610	!
611	! last category
612	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
613	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
614	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
615	!
616	! correction if concentration of upper cat is greater than lower cat
617	! (it should be a gaussian around jl0 but sometimes it is not)
618	IF ( jl0 /= jpl ) THEN
619	DO jl = jpl, jl0+1, -1
620	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
621	zdv = zh_i(ji,jl) * za_i(ji,jl)
622	zh_i(ji,jl ) = 0._wp
623	za_i (ji,jl ) = 0._wp
624	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
625	END IF
626	END DO
627	ENDIF
628	!
629	ENDIF
630	!
631	! Compatibility tests
632	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
633	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
634	!
635	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
636	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
637	!
638	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
639	!
640	itest(4) = 1
641	DO jl = 1, i_fill
642	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
643	END DO
644	! !----------------------------
645	END DO ! end iteration on categories
646	! !----------------------------
647	ENDIF
648	END DO
649
650	! Add Snow in each category where za_i is not 0
651	DO jl = 1, jpl
652	DO ji = 1, idim
653	IF( za_i(ji,jl) > 0._wp ) THEN
654	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
655	! In case snow load is in excess that would lead to transformation from snow to ice
656	! Then, transfer the snow excess into the ice (different from icethd_dh)
657	zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
658	! recompute h_i, h_s avoiding out of bounds values
659	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
660	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
661	ENDIF
662	END DO
663	END DO
664	!
665	END SUBROUTINE ice_var_itd
666
667
668	SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
669	!!-------------------------------------------------------------------
670	!! * ROUTINE ice_var_itd2 *
671	!!
672	!! ** Purpose : converting N-cat ice to jpl ice categories
673	!!
674	!! ice thickness distribution follows a gaussian law
675	!! around the concentration of the most likely ice thickness
676	!! (similar as iceistate.F90)
677	!!
678	!! ** Method: Iterative procedure
679	!!
680	!! 1) Fill ice cat that correspond to input thicknesses
681	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
682	!!
683	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
684	!! by removing 25% ice area from jlmin and jlmax (resp.)
685	!!
686	!! 3) Expand the filling to the empty cat between jlmin and jlmax
687	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
688	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
689	!!
690	!! ** Arguments : zhti: N-cat ice thickness
691	!! zhts: N-cat snow depth
692	!! zati: N-cat ice concentration
693	!!
694	!! ** Output : jpl-cat
695	!!
696	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
697	!!-------------------------------------------------------------------
698	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
699	INTEGER :: idim, icat
700	INTEGER, PARAMETER :: ztrans = 0.25_wp
701	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
702	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
703	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
704	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
705	!!-------------------------------------------------------------------
706	!
707	idim = SIZE( zhti, 1 )
708	icat = SIZE( zhti, 2 )
709	!
710	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
711	ALLOCATE( jlmin(idim), jlmax(idim) )
712
713	! --- initialize output fields to 0 --- !
714	zh_i(1:idim,1:jpl) = 0._wp
715	zh_s(1:idim,1:jpl) = 0._wp
716	za_i(1:idim,1:jpl) = 0._wp
717	!
718	! --- fill the categories --- !
719	! find where cat-input = cat-output and fill cat-output fields
720	jlmax(:) = 0
721	jlmin(:) = 999
722	jlfil(:,:) = 0
723	DO jl1 = 1, jpl
724	DO jl2 = 1, icat
725	DO ji = 1, idim
726	IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
727	! fill the right category
728	zh_i(ji,jl1) = zhti(ji,jl2)
729	zh_s(ji,jl1) = zhts(ji,jl2)
730	za_i(ji,jl1) = zati(ji,jl2)
731	! record categories that are filled
732	jlmax(ji) = MAX( jlmax(ji), jl1 )
733	jlmin(ji) = MIN( jlmin(ji), jl1 )
734	jlfil(ji,jl1) = jl1
735	ENDIF
736	END DO
737	END DO
738	END DO
739	!
740	! --- fill the gaps between categories --- !
741	! transfer from categories filled at the previous step to the empty ones in between
742	DO ji = 1, idim
743	jl1 = jlmin(ji)
744	jl2 = jlmax(ji)
745	IF( jl1 > 1 ) THEN
746	! fill the lower cat (jl1-1)
747	za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
748	zh_i(ji,jl1-1) = hi_mean(jl1-1)
749	! remove from cat jl1
750	za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
751	ENDIF
752	IF( jl2 < jpl ) THEN
753	! fill the upper cat (jl2+1)
754	za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
755	zh_i(ji,jl2+1) = hi_mean(jl2+1)
756	! remove from cat jl2
757	za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
758	ENDIF
759	END DO
760	!
761	jlfil2(:,:) = jlfil(:,:)
762	! fill categories from low to high
763	DO jl = 2, jpl-1
764	DO ji = 1, idim
765	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
766	! fill high
767	za_i(ji,jl) = ztrans * za_i(ji,jl-1)
768	zh_i(ji,jl) = hi_mean(jl)
769	jlfil(ji,jl) = jl
770	! remove low
771	za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
772	ENDIF
773	END DO
774	END DO
775	!
776	! fill categories from high to low
777	DO jl = jpl-1, 2, -1
778	DO ji = 1, idim
779	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
780	! fill low
781	za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
782	zh_i(ji,jl) = hi_mean(jl)
783	jlfil2(ji,jl) = jl
784	! remove high
785	za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
786	ENDIF
787	END DO
788	END DO
789	!
790	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
791	DEALLOCATE( jlmin, jlmax )
792	!
793	END SUBROUTINE ice_var_itd2
794
795
796	SUBROUTINE ice_var_bv
797	!!-------------------------------------------------------------------
798	!! * ROUTINE ice_var_bv *
799	!!
800	!! ** Purpose : computes mean brine volume (%) in sea ice
801	!!
802	!! ** Method : e = - 0.054 * S (ppt) / T (C)
803	!!
804	!! References : Vancoppenolle et al., JGR, 2007
805	!!-------------------------------------------------------------------
806	INTEGER :: ji, jj, jk, jl ! dummy loop indices
807	!!-------------------------------------------------------------------
808	!
809	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
810	!! instead of setting everything to zero as just below
811	bv_i (:,:,:) = 0._wp
812	DO jl = 1, jpl
813	DO jk = 1, nlay_i
814	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
815	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
816	END WHERE
817	END DO
818	END DO
819	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
820	ELSEWHERE ; bvm_i(:,:) = 0._wp
821	END WHERE
822	!
823	END SUBROUTINE ice_var_bv
824
825
826	SUBROUTINE ice_var_enthalpy
827	!!-------------------------------------------------------------------
828	!! * ROUTINE ice_var_enthalpy *
829	!!
830	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
831	!!
832	!! ** Method : Formula (Bitz and Lipscomb, 1999)
833	!!-------------------------------------------------------------------
834	INTEGER :: ji, jk ! dummy loop indices
835	REAL(wp) :: ztmelts ! local scalar
836	!!-------------------------------------------------------------------
837	!
838	DO jk = 1, nlay_i ! Sea ice energy of melting
839	DO ji = 1, npti
840	ztmelts = - tmut * sz_i_1d(ji,jk)
841	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
842	! (sometimes zdf scheme produces abnormally high temperatures)
843	e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
844	& + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
845	& - rcp * ztmelts )
846	END DO
847	END DO
848	DO jk = 1, nlay_s ! Snow energy of melting
849	DO ji = 1, npti
850	e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
851	END DO
852	END DO
853	!
854	END SUBROUTINE ice_var_enthalpy
855
856	#else
857	!!----------------------------------------------------------------------
858	!! Default option Dummy module NO ESIM sea-ice model
859	!!----------------------------------------------------------------------
860	#endif
861
862	!!======================================================================
863	END MODULE icevar

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

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

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8984

Download in other formats: