Context Navigation

icevar.F90 @ 8787

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

Download in other formats: