Context Navigation

icevar.F90 @ 8813

Last change on this file since 8813 was 8813, checked in by gm, 6 years ago
#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787
File size: 37.5 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
367	SUBROUTINE ice_var_salprof1d
368	!!-------------------------------------------------------------------
369	!! * ROUTINE ice_var_salprof1d *
370	!!
371	!! ** Purpose : 1d computation of the sea ice salinity profile
372	!! Works with 1d vectors and is used by thermodynamic modules
373	!!-------------------------------------------------------------------
374	INTEGER :: ji, jk ! dummy loop indices
375	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
376	REAL(wp) :: zs, zs0 ! - -
377	!
378	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
379	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
380	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
381	!!-------------------------------------------------------------------
382	!
383	SELECT CASE ( nn_icesal )
384	!
385	! !---------------------------------------!
386	CASE( 1 ) ! constant salinity in time and space !
387	! !---------------------------------------!
388	sz_i_1d(1:npti,:) = rn_icesal
389	!
390	! !---------------------------------------------!
391	CASE( 2 ) ! time varying salinity with linear profile !
392	! !---------------------------------------------!
393	!
394	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
395	!
396	! ! Slope of the linear profile
397	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
398	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
399	END WHERE
400
401	z1_dS = 1._wp / ( zsi1 - zsi0 )
402	DO ji = 1, npti
403	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
404	! ! force a constant profile when SSS too low (Baltic Sea)
405	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
406	END DO
407	!
408	! Computation of the profile
409	DO jk = 1, nlay_i
410	DO ji = 1, npti
411	! ! linear profile with 0 surface value
412	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
413	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
414	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
415	END DO
416	END DO
417	!
418	DEALLOCATE( z_slope_s, zalpha )
419
420	! !-------------------------------------------!
421	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
422	! !-------------------------------------------! (mean = 2.30)
423	!
424	s_i_1d(1:npti) = 2.30_wp
425	!
426	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
427	DO jk = 1, nlay_i
428	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
429	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
430	DO ji = 1, npti
431	sz_i_1d(ji,jk) = zsal
432	END DO
433	END DO
434	!
435	END SELECT
436	!
437	END SUBROUTINE ice_var_salprof1d
438
439
440	SUBROUTINE ice_var_zapsmall
441	!!-------------------------------------------------------------------
442	!! * ROUTINE ice_var_zapsmall *
443	!!
444	!! ** Purpose : Remove too small sea ice areas and correct fluxes
445	!!-------------------------------------------------------------------
446	INTEGER :: ji, jj, jl, jk ! dummy loop indices
447	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
448	!!-------------------------------------------------------------------
449	!
450	DO jl = 1, jpl !== loop over the categories ==!
451	!
452	!-----------------------------------------------------------------
453	! Zap ice energy and use ocean heat to melt ice
454	!-----------------------------------------------------------------
455	WHERE( a_i(:,:,jl) > epsi20 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
456	ELSEWHERE ; h_i(:,:,jl) = 0._wp
457	END WHERE
458	!
459	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
460	ELSEWHERE ; zswitch(:,:) = 1._wp
461	END WHERE
462	!
463	DO jk = 1, nlay_i
464	DO jj = 1 , jpj
465	DO ji = 1 , jpi
466	! update exchanges with ocean
467	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
468	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
469	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
470	END DO
471	END DO
472	END DO
473	!
474	DO jj = 1 , jpj
475	DO ji = 1 , jpi
476	! update exchanges with ocean
477	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_rdtice
478	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
479	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
480	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
481	IF( ln_pnd_fwb ) THEN
482	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
483	ENDIF
484	!-----------------------------------------------------------------
485	! Zap snow energy
486	!-----------------------------------------------------------------
487	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
488	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
489	!
490	!-----------------------------------------------------------------
491	! zap ice and snow volume, add water and salt to ocean
492	!-----------------------------------------------------------------
493	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
494	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
495	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
496	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
497	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
498	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
499	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
500	!
501	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
502	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
503	!
504	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
505	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
506	!
507	END DO
508	END DO
509	!
510	END DO
511
512	! to be sure that at_i is the sum of a_i(jl)
513	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
514	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
515
516	! open water = 1 if at_i=0
517	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
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 :: idim, 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	idim = SIZE( zhti , 1 )
567	zh_i(1:idim,1:jpl) = 0._wp
568	zh_s(1:idim,1:jpl) = 0._wp
569	za_i(1:idim,1:jpl) = 0._wp
570	!
571	DO ji = 1, idim
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	END DO
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, idim
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
669	SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
670	!!-------------------------------------------------------------------
671	!! * ROUTINE ice_var_itd2 *
672	!!
673	!! ** Purpose : converting N-cat ice to jpl ice categories
674	!!
675	!! ice thickness distribution follows a gaussian law
676	!! around the concentration of the most likely ice thickness
677	!! (similar as iceistate.F90)
678	!!
679	!! ** Method: Iterative procedure
680	!!
681	!! 1) Fill ice cat that correspond to input thicknesses
682	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
683	!!
684	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
685	!! by removing 25% ice area from jlmin and jlmax (resp.)
686	!!
687	!! 3) Expand the filling to the empty cat between jlmin and jlmax
688	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
689	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
690	!!
691	!! ** Arguments : zhti: N-cat ice thickness
692	!! zhts: N-cat snow depth
693	!! zati: N-cat ice concentration
694	!!
695	!! ** Output : jpl-cat
696	!!
697	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
698	!!-------------------------------------------------------------------
699	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
700	INTEGER :: idim, icat
701	INTEGER, PARAMETER :: ztrans = 0.25_wp
702	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
703	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
704	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
705	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
706	!!-------------------------------------------------------------------
707	!
708	idim = SIZE( zhti, 1 )
709	icat = SIZE( zhti, 2 )
710	!
711	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
712	ALLOCATE( jlmin(idim), jlmax(idim) )
713
714	! --- initialize output fields to 0 --- !
715	zh_i(1:idim,1:jpl) = 0._wp
716	zh_s(1:idim,1:jpl) = 0._wp
717	za_i(1:idim,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, icat
726	DO ji = 1, idim
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, idim
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, idim
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, idim
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( jlfil, jlfil2 ) ! deallocate arrays
792	DEALLOCATE( jlmin, jlmax )
793	!
794	END SUBROUTINE ice_var_itd2
795
796
797	SUBROUTINE ice_var_bv
798	!!-------------------------------------------------------------------
799	!! * ROUTINE ice_var_bv *
800	!!
801	!! ** Purpose : computes mean brine volume (%) in sea ice
802	!!
803	!! ** Method : e = - 0.054 * S (ppt) / T (C)
804	!!
805	!! References : Vancoppenolle et al., JGR, 2007
806	!!-------------------------------------------------------------------
807	INTEGER :: ji, jj, jk, jl ! dummy loop indices
808	!!-------------------------------------------------------------------
809	!
810	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
811	!! instead of setting everything to zero as just below
812	bv_i (:,:,:) = 0._wp
813	DO jl = 1, jpl
814	DO jk = 1, nlay_i
815	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
816	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
817	END WHERE
818	END DO
819	END DO
820	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
821	ELSEWHERE ; bvm_i(:,:) = 0._wp
822	END WHERE
823	!
824	END SUBROUTINE ice_var_bv
825
826
827	#else
828	!!----------------------------------------------------------------------
829	!! Default option Dummy module NO ESIM sea-ice model
830	!!----------------------------------------------------------------------
831	#endif
832
833	!!======================================================================
834	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_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8813

Download in other formats: