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 @ 8597

Last change on this file since 8597 was 8597, checked in by clem, 7 years ago
first step to make melt ponds compliant with the new code
File size: 32.1 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	!
261	END SUBROUTINE ice_var_eqv2glo
262
263
264	SUBROUTINE ice_var_salprof
265	!!-------------------------------------------------------------------
266	!! * ROUTINE ice_var_salprof *
267	!!
268	!! ** Purpose : computes salinity profile in function of bulk salinity
269	!!
270	!! ** Method : If bulk salinity greater than zsi1,
271	!! the profile is assumed to be constant (S_inf)
272	!! If bulk salinity lower than zsi0,
273	!! the profile is linear with 0 at the surface (S_zero)
274	!! If it is between zsi0 and zsi1, it is a
275	!! alpha-weighted linear combination of s_inf and s_zero
276	!!
277	!! ** References : Vancoppenolle et al., 2007
278	!!-------------------------------------------------------------------
279	INTEGER :: ji, jj, jk, jl ! dummy loop index
280	REAL(wp) :: zsal, z1_dS
281	REAL(wp) :: zargtemp , zs0, zs
282	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
283	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
284	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
285	!!-------------------------------------------------------------------
286
287	!!gm Question: Remove the option 3 ? How many years since it last use ?
288
289	SELECT CASE ( nn_icesal )
290	!
291	! !---------------------------------------!
292	CASE( 1 ) ! constant salinity in time and space !
293	! !---------------------------------------!
294	sz_i(:,:,:,:) = rn_icesal
295	s_i(:,:,:) = rn_icesal
296	!
297	! !---------------------------------------------!
298	CASE( 2 ) ! time varying salinity with linear profile !
299	! !---------------------------------------------!
300	!
301	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
302	!
303	DO jl = 1, jpl
304	DO jk = 1, nlay_i
305	sz_i(:,:,jk,jl) = s_i(:,:,jl)
306	END DO
307	END DO
308	! ! Slope of the linear profile
309	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
310	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
311	END WHERE
312	!
313	z1_dS = 1._wp / ( zsi1 - zsi0 )
314	DO jl = 1, jpl
315	DO jj = 1, jpj
316	DO ji = 1, jpi
317	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
318	! ! force a constant profile when SSS too low (Baltic Sea)
319	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
320	END DO
321	END DO
322	END DO
323
324	! Computation of the profile
325	DO jl = 1, jpl
326	DO jk = 1, nlay_i
327	DO jj = 1, jpj
328	DO ji = 1, jpi
329	! ! linear profile with 0 surface value
330	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
331	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
332	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
333	END DO
334	END DO
335	END DO
336	END DO
337	!
338	DEALLOCATE( z_slope_s , zalpha )
339	!
340	! !-------------------------------------------!
341	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
342	! !-------------------------------------------! (mean = 2.30)
343	!
344	s_i(:,:,:) = 2.30_wp
345	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
346	!! a array S_prof(1:nlay_i) containing the calculation and just do:
347	! DO jk = 1, nlay_i
348	! sz_i(:,:,jk,:) = S_prof(jk)
349	! END DO
350	!!gm end
351	!
352	DO jl = 1, jpl
353	DO jk = 1, nlay_i
354	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
355	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
356	END DO
357	END DO
358	!
359	END SELECT
360	!
361	END SUBROUTINE ice_var_salprof
362
363	SUBROUTINE ice_var_salprof1d
364	!!-------------------------------------------------------------------
365	!! * ROUTINE ice_var_salprof1d *
366	!!
367	!! ** Purpose : 1d computation of the sea ice salinity profile
368	!! Works with 1d vectors and is used by thermodynamic modules
369	!!-------------------------------------------------------------------
370	INTEGER :: ji, jk ! dummy loop indices
371	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
372	REAL(wp) :: zs, zs0 ! - -
373	!
374	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
375	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
376	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
377	!!-------------------------------------------------------------------
378	!
379	SELECT CASE ( nn_icesal )
380	!
381	! !---------------------------------------!
382	CASE( 1 ) ! constant salinity in time and space !
383	! !---------------------------------------!
384	sz_i_1d(1:npti,:) = rn_icesal
385	!
386	! !---------------------------------------------!
387	CASE( 2 ) ! time varying salinity with linear profile !
388	! !---------------------------------------------!
389	!
390	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
391	!
392	! ! Slope of the linear profile
393	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
394	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
395	END WHERE
396
397	z1_dS = 1._wp / ( zsi1 - zsi0 )
398	DO ji = 1, npti
399	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
400	! ! force a constant profile when SSS too low (Baltic Sea)
401	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
402	END DO
403	!
404	! Computation of the profile
405	DO jk = 1, nlay_i
406	DO ji = 1, npti
407	! ! linear profile with 0 surface value
408	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
409	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
410	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
411	END DO
412	END DO
413	!
414	DEALLOCATE( z_slope_s, zalpha )
415
416	! !-------------------------------------------!
417	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
418	! !-------------------------------------------! (mean = 2.30)
419	!
420	s_i_1d(1:npti) = 2.30_wp
421	!
422	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
423	DO jk = 1, nlay_i
424	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
425	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
426	DO ji = 1, npti
427	sz_i_1d(ji,jk) = zsal
428	END DO
429	END DO
430	!
431	END SELECT
432	!
433	END SUBROUTINE ice_var_salprof1d
434
435
436	SUBROUTINE ice_var_zapsmall
437	!!-------------------------------------------------------------------
438	!! * ROUTINE ice_var_zapsmall *
439	!!
440	!! ** Purpose : Remove too small sea ice areas and correct fluxes
441	!!-------------------------------------------------------------------
442	INTEGER :: ji, jj, jl, jk ! dummy loop indices
443	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
444	!!-------------------------------------------------------------------
445	!
446	DO jl = 1, jpl !== loop over the categories ==!
447	!
448	!-----------------------------------------------------------------
449	! Zap ice energy and use ocean heat to melt ice
450	!-----------------------------------------------------------------
451	WHERE( a_i(:,:,jl) > epsi20 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
452	ELSEWHERE ; h_i(:,:,jl) = 0._wp
453	END WHERE
454	!
455	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
456	ELSEWHERE ; zswitch(:,:) = 1._wp
457	END WHERE
458
459	DO jk = 1, nlay_i
460	DO jj = 1 , jpj
461	DO ji = 1 , jpi
462	! update exchanges with ocean
463	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
464	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
465	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
466	END DO
467	END DO
468	END DO
469
470	DO jj = 1 , jpj
471	DO ji = 1 , jpi
472	! update exchanges with ocean
473	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_rdtice
474	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
475	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
476	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
477	IF( ln_pnd_fwb ) THEN
478	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
479	ENDIF
480	!-----------------------------------------------------------------
481	! Zap snow energy
482	!-----------------------------------------------------------------
483	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
484	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
485
486	!-----------------------------------------------------------------
487	! zap ice and snow volume, add water and salt to ocean
488	!-----------------------------------------------------------------
489	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
490	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
491	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
492	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
493	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
494	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
495	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
496
497	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
498	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
499
500	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
501	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
502
503	END DO
504	END DO
505
506	END DO
507
508	! to be sure that at_i is the sum of a_i(jl)
509	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
510	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
511
512	! open water = 1 if at_i=0
513	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
514	!
515	END SUBROUTINE ice_var_zapsmall
516
517
518	SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
519	!!-------------------------------------------------------------------
520	!! * ROUTINE ice_var_itd *
521	!!
522	!! ** Purpose : converting 1-cat ice to multiple ice categories
523	!!
524	!! ice thickness distribution follows a gaussian law
525	!! around the concentration of the most likely ice thickness
526	!! (similar as iceistate.F90)
527	!!
528	!! ** Method: Iterative procedure
529	!!
530	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
531	!!
532	!! 2) Check whether the distribution conserves area and volume, positivity and
533	!! category boundaries
534	!!
535	!! 3) If not (input ice is too thin), the last category is empty and
536	!! the number of categories is reduced (jpl-1)
537	!!
538	!! 4) Iterate until ok (SUM(itest(:) = 4)
539	!!
540	!! ** Arguments : zhti: 1-cat ice thickness
541	!! zhts: 1-cat snow depth
542	!! zati: 1-cat ice concentration
543	!!
544	!! ** Output : jpl-cat
545	!!
546	!! (Example of application: BDY forcings when input are cell averaged)
547	!!-------------------------------------------------------------------
548	INTEGER :: ji, jk, jl ! dummy loop indices
549	INTEGER :: ijpij, i_fill, jl0
550	REAL(wp) :: zarg, zV, zconv, zdh, zdv
551	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
552	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
553	INTEGER , DIMENSION(4) :: itest
554	!!-------------------------------------------------------------------
555	!
556	! ----------------------------------------
557	! distribution over the jpl ice categories
558	! ----------------------------------------
559	! a gaussian distribution for ice concentration is used
560	! then we check whether the distribution fullfills
561	! volume and area conservation, positivity and ice categories bounds
562	ijpij = SIZE( zhti , 1 )
563	zh_i(1:ijpij,1:jpl) = 0._wp
564	zh_s(1:ijpij,1:jpl) = 0._wp
565	za_i (1:ijpij,1:jpl) = 0._wp
566
567	DO ji = 1, ijpij
568
569	IF( zhti(ji) > 0._wp ) THEN
570
571	! find which category (jl0) the input ice thickness falls into
572	jl0 = jpl
573	DO jl = 1, jpl
574	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
575	jl0 = jl
576	CYCLE
577	ENDIF
578	END DO
579
580	itest(:) = 0
581	i_fill = jpl + 1 !------------------------------------
582	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
583	! !------------------------------------
584	i_fill = i_fill - 1
585	!
586	zh_i(ji,1:jpl) = 0._wp
587	za_i (ji,1:jpl) = 0._wp
588	itest(:) = 0
589
590	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
591	zh_i(ji,1) = zhti(ji)
592	za_i (ji,1) = zati (ji)
593	ELSE !-- case ice is thicker: fill categories >1
594	! thickness
595	DO jl = 1, i_fill - 1
596	zh_i(ji,jl) = hi_mean(jl)
597	END DO
598
599	! concentration
600	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
601	DO jl = 1, i_fill - 1
602	IF ( jl /= jl0 ) THEN
603	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
604	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
605	ENDIF
606	END DO
607
608	! last category
609	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
610	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
611	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
612
613	! clem: correction if concentration of upper cat is greater than lower cat
614	! (it should be a gaussian around jl0 but sometimes it is not)
615	IF ( jl0 /= jpl ) THEN
616	DO jl = jpl, jl0+1, -1
617	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
618	zdv = zh_i(ji,jl) * za_i(ji,jl)
619	zh_i(ji,jl ) = 0._wp
620	za_i (ji,jl ) = 0._wp
621	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
622	END IF
623	ENDDO
624	ENDIF
625
626	ENDIF
627
628	! Compatibility tests
629	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
630	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
631
632	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
633	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
634
635	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
636
637	itest(4) = 1
638	DO jl = 1, i_fill
639	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
640	END DO
641	! !----------------------------
642	END DO ! end iteration on categories
643	! !----------------------------
644	ENDIF
645	END DO
646
647	! Add Snow in each category where za_i is not 0
648	DO jl = 1, jpl
649	DO ji = 1, ijpij
650	IF( za_i(ji,jl) > 0._wp ) THEN
651	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
652	! In case snow load is in excess that would lead to transformation from snow to ice
653	! Then, transfer the snow excess into the ice (different from icethd_dh)
654	zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
655	! recompute h_i, h_s avoiding out of bounds values
656	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
657	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
658	ENDIF
659	END DO
660	END DO
661	!
662	END SUBROUTINE ice_var_itd
663
664
665	SUBROUTINE ice_var_bv
666	!!-------------------------------------------------------------------
667	!! * ROUTINE ice_var_bv *
668	!!
669	!! ** Purpose : computes mean brine volume (%) in sea ice
670	!!
671	!! ** Method : e = - 0.054 * S (ppt) / T (C)
672	!!
673	!! References : Vancoppenolle et al., JGR, 2007
674	!!-------------------------------------------------------------------
675	INTEGER :: ji, jj, jk, jl ! dummy loop indices
676	!!-------------------------------------------------------------------
677	!
678	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
679	!! instead of setting everything to zero as just below
680	bv_i (:,:,:) = 0._wp
681	DO jl = 1, jpl
682	DO jk = 1, nlay_i
683	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
684	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
685	END WHERE
686	END DO
687	END DO
688	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
689	ELSEWHERE ; bvm_i(:,:) = 0._wp
690	END WHERE
691	!
692	END SUBROUTINE ice_var_bv
693
694
695	#else
696	!!----------------------------------------------------------------------
697	!! Default option Dummy module NO ESIM sea-ice model
698	!!----------------------------------------------------------------------
699	#endif
700
701	!!======================================================================
702	END MODULE icevar

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

Download in other formats: