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_merge_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 9271

Last change on this file since 9271 was 9271, checked in by clem, 6 years ago
first steps for having more than 1 snow layers in the ice (in theory). There is still icethd_dh.F90 routine to change
File size: 38.9 KB

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

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

Download in other formats: