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

Context Navigation

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8637

Last change on this file since 8637 was 8637, checked in by gm, 7 years ago
#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626
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	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 :: ijpij, 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	ijpij = SIZE( zhti , 1 )
564	zh_i(1:ijpij,1:jpl) = 0._wp
565	zh_s(1:ijpij,1:jpl) = 0._wp
566	za_i (1:ijpij,1:jpl) = 0._wp
567
568	DO ji = 1, ijpij
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	! clem: 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	ENDDO
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, ijpij
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_bv
667	!!-------------------------------------------------------------------
668	!! * ROUTINE ice_var_bv *
669	!!
670	!! ** Purpose : computes mean brine volume (%) in sea ice
671	!!
672	!! ** Method : e = - 0.054 * S (ppt) / T (C)
673	!!
674	!! References : Vancoppenolle et al., JGR, 2007
675	!!-------------------------------------------------------------------
676	INTEGER :: ji, jj, jk, jl ! dummy loop indices
677	!!-------------------------------------------------------------------
678	!
679	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
680	!! instead of setting everything to zero as just below
681	bv_i (:,:,:) = 0._wp
682	DO jl = 1, jpl
683	DO jk = 1, nlay_i
684	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
685	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
686	END WHERE
687	END DO
688	END DO
689	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
690	ELSEWHERE ; bvm_i(:,:) = 0._wp
691	END WHERE
692	!
693	END SUBROUTINE ice_var_bv
694
695
696	#else
697	!!----------------------------------------------------------------------
698	!! Default option Dummy module NO ESIM sea-ice model
699	!!----------------------------------------------------------------------
700	#endif
701
702	!!======================================================================
703	END MODULE icevar

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

Download in other formats: