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

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8623

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

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

Download in other formats: