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

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

icevar.F90 in branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

Last change on this file since 8906 was 8906, checked in by clem, 6 years ago
make melt ponds from Holland2012 and associated freshwater flux conservative
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
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	!-----------------------------------------------------------------
482	! Zap snow energy
483	!-----------------------------------------------------------------
484	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
485	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
486	!
487	!-----------------------------------------------------------------
488	! zap ice and snow volume, add water and salt to ocean
489	!-----------------------------------------------------------------
490	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
491	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
492	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
493	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
494	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
495	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
496	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
497	!
498	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
499	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
500	!
501	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
502	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
503	!
504	END DO
505	END DO
506	!
507	END DO
508
509	! to be sure that at_i is the sum of a_i(jl)
510	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
511	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
512
513	! open water = 1 if at_i=0
514	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
515	!
516	END SUBROUTINE ice_var_zapsmall
517
518
519	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
520	!!-------------------------------------------------------------------
521	!! * ROUTINE ice_var_itd *
522	!!
523	!! ** Purpose : converting 1-cat ice to multiple ice categories
524	!!
525	!! ice thickness distribution follows a gaussian law
526	!! around the concentration of the most likely ice thickness
527	!! (similar as iceistate.F90)
528	!!
529	!! ** Method: Iterative procedure
530	!!
531	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
532	!!
533	!! 2) Check whether the distribution conserves area and volume, positivity and
534	!! category boundaries
535	!!
536	!! 3) If not (input ice is too thin), the last category is empty and
537	!! the number of categories is reduced (jpl-1)
538	!!
539	!! 4) Iterate until ok (SUM(itest(:) = 4)
540	!!
541	!! ** Arguments : zhti: 1-cat ice thickness
542	!! zhts: 1-cat snow depth
543	!! zati: 1-cat ice concentration
544	!!
545	!! ** Output : jpl-cat
546	!!
547	!! (Example of application: BDY forcings when input are cell averaged)
548	!!-------------------------------------------------------------------
549	INTEGER :: ji, jk, jl ! dummy loop indices
550	INTEGER :: idim, i_fill, jl0
551	REAL(wp) :: zarg, zV, zconv, zdh, zdv
552	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
553	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
554	INTEGER , DIMENSION(4) :: itest
555	!!-------------------------------------------------------------------
556	!
557	! ----------------------------------------
558	! distribution over the jpl ice categories
559	! ----------------------------------------
560	! a gaussian distribution for ice concentration is used
561	! then we check whether the distribution fullfills
562	! volume and area conservation, positivity and ice categories bounds
563	idim = SIZE( zhti , 1 )
564	zh_i(1:idim,1:jpl) = 0._wp
565	zh_s(1:idim,1:jpl) = 0._wp
566	za_i(1:idim,1:jpl) = 0._wp
567	!
568	DO ji = 1, idim
569	!
570	IF( zhti(ji) > 0._wp ) THEN
571	!
572	! find which category (jl0) the input ice thickness falls into
573	jl0 = jpl
574	DO jl = 1, jpl
575	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
576	jl0 = jl
577	CYCLE
578	ENDIF
579	END DO
580	!
581	itest(:) = 0
582	i_fill = jpl + 1 !------------------------------------
583	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
584	! !------------------------------------
585	i_fill = i_fill - 1
586	!
587	zh_i(ji,1:jpl) = 0._wp
588	za_i(ji,1:jpl) = 0._wp
589	itest(:) = 0
590	!
591	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
592	zh_i(ji,1) = zhti(ji)
593	za_i (ji,1) = zati (ji)
594	ELSE !-- case ice is thicker: fill categories >1
595	! thickness
596	DO jl = 1, i_fill - 1
597	zh_i(ji,jl) = hi_mean(jl)
598	END DO
599	!
600	! concentration
601	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
602	DO jl = 1, i_fill - 1
603	IF ( jl /= jl0 ) THEN
604	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
605	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
606	ENDIF
607	END DO
608	!
609	! last category
610	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
611	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
612	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
613	!
614	! correction if concentration of upper cat is greater than lower cat
615	! (it should be a gaussian around jl0 but sometimes it is not)
616	IF ( jl0 /= jpl ) THEN
617	DO jl = jpl, jl0+1, -1
618	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
619	zdv = zh_i(ji,jl) * za_i(ji,jl)
620	zh_i(ji,jl ) = 0._wp
621	za_i (ji,jl ) = 0._wp
622	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
623	END IF
624	END DO
625	ENDIF
626	!
627	ENDIF
628	!
629	! Compatibility tests
630	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
631	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
632	!
633	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
634	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
635	!
636	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
637	!
638	itest(4) = 1
639	DO jl = 1, i_fill
640	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
641	END DO
642	! !----------------------------
643	END DO ! end iteration on categories
644	! !----------------------------
645	ENDIF
646	END DO
647
648	! Add Snow in each category where za_i is not 0
649	DO jl = 1, jpl
650	DO ji = 1, idim
651	IF( za_i(ji,jl) > 0._wp ) THEN
652	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
653	! In case snow load is in excess that would lead to transformation from snow to ice
654	! Then, transfer the snow excess into the ice (different from icethd_dh)
655	zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
656	! recompute h_i, h_s avoiding out of bounds values
657	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
658	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
659	ENDIF
660	END DO
661	END DO
662	!
663	END SUBROUTINE ice_var_itd
664
665
666	SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
667	!!-------------------------------------------------------------------
668	!! * ROUTINE ice_var_itd2 *
669	!!
670	!! ** Purpose : converting N-cat ice to jpl ice categories
671	!!
672	!! ice thickness distribution follows a gaussian law
673	!! around the concentration of the most likely ice thickness
674	!! (similar as iceistate.F90)
675	!!
676	!! ** Method: Iterative procedure
677	!!
678	!! 1) Fill ice cat that correspond to input thicknesses
679	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
680	!!
681	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
682	!! by removing 25% ice area from jlmin and jlmax (resp.)
683	!!
684	!! 3) Expand the filling to the empty cat between jlmin and jlmax
685	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
686	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
687	!!
688	!! ** Arguments : zhti: N-cat ice thickness
689	!! zhts: N-cat snow depth
690	!! zati: N-cat ice concentration
691	!!
692	!! ** Output : jpl-cat
693	!!
694	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
695	!!-------------------------------------------------------------------
696	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
697	INTEGER :: idim, icat
698	INTEGER, PARAMETER :: ztrans = 0.25_wp
699	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
700	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
701	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
702	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
703	!!-------------------------------------------------------------------
704	!
705	idim = SIZE( zhti, 1 )
706	icat = SIZE( zhti, 2 )
707	!
708	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
709	ALLOCATE( jlmin(idim), jlmax(idim) )
710
711	! --- initialize output fields to 0 --- !
712	zh_i(1:idim,1:jpl) = 0._wp
713	zh_s(1:idim,1:jpl) = 0._wp
714	za_i(1:idim,1:jpl) = 0._wp
715	!
716	! --- fill the categories --- !
717	! find where cat-input = cat-output and fill cat-output fields
718	jlmax(:) = 0
719	jlmin(:) = 999
720	jlfil(:,:) = 0
721	DO jl1 = 1, jpl
722	DO jl2 = 1, icat
723	DO ji = 1, idim
724	IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
725	! fill the right category
726	zh_i(ji,jl1) = zhti(ji,jl2)
727	zh_s(ji,jl1) = zhts(ji,jl2)
728	za_i(ji,jl1) = zati(ji,jl2)
729	! record categories that are filled
730	jlmax(ji) = MAX( jlmax(ji), jl1 )
731	jlmin(ji) = MIN( jlmin(ji), jl1 )
732	jlfil(ji,jl1) = jl1
733	ENDIF
734	END DO
735	END DO
736	END DO
737	!
738	! --- fill the gaps between categories --- !
739	! transfer from categories filled at the previous step to the empty ones in between
740	DO ji = 1, idim
741	jl1 = jlmin(ji)
742	jl2 = jlmax(ji)
743	IF( jl1 > 1 ) THEN
744	! fill the lower cat (jl1-1)
745	za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
746	zh_i(ji,jl1-1) = hi_mean(jl1-1)
747	! remove from cat jl1
748	za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
749	ENDIF
750	IF( jl2 < jpl ) THEN
751	! fill the upper cat (jl2+1)
752	za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
753	zh_i(ji,jl2+1) = hi_mean(jl2+1)
754	! remove from cat jl2
755	za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
756	ENDIF
757	END DO
758	!
759	jlfil2(:,:) = jlfil(:,:)
760	! fill categories from low to high
761	DO jl = 2, jpl-1
762	DO ji = 1, idim
763	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
764	! fill high
765	za_i(ji,jl) = ztrans * za_i(ji,jl-1)
766	zh_i(ji,jl) = hi_mean(jl)
767	jlfil(ji,jl) = jl
768	! remove low
769	za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
770	ENDIF
771	END DO
772	END DO
773	!
774	! fill categories from high to low
775	DO jl = jpl-1, 2, -1
776	DO ji = 1, idim
777	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
778	! fill low
779	za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
780	zh_i(ji,jl) = hi_mean(jl)
781	jlfil2(ji,jl) = jl
782	! remove high
783	za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
784	ENDIF
785	END DO
786	END DO
787	!
788	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
789	DEALLOCATE( jlmin, jlmax )
790	!
791	END SUBROUTINE ice_var_itd2
792
793
794	SUBROUTINE ice_var_bv
795	!!-------------------------------------------------------------------
796	!! * ROUTINE ice_var_bv *
797	!!
798	!! ** Purpose : computes mean brine volume (%) in sea ice
799	!!
800	!! ** Method : e = - 0.054 * S (ppt) / T (C)
801	!!
802	!! References : Vancoppenolle et al., JGR, 2007
803	!!-------------------------------------------------------------------
804	INTEGER :: ji, jj, jk, jl ! dummy loop indices
805	!!-------------------------------------------------------------------
806	!
807	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
808	!! instead of setting everything to zero as just below
809	bv_i (:,:,:) = 0._wp
810	DO jl = 1, jpl
811	DO jk = 1, nlay_i
812	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
813	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
814	END WHERE
815	END DO
816	END DO
817	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
818	ELSEWHERE ; bvm_i(:,:) = 0._wp
819	END WHERE
820	!
821	END SUBROUTINE ice_var_bv
822
823
824	#else
825	!!----------------------------------------------------------------------
826	!! Default option Dummy module NO ESIM sea-ice model
827	!!----------------------------------------------------------------------
828	#endif
829
830	!!======================================================================
831	END MODULE icevar

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

Download in other formats: