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

Context Navigation

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8738

Last change on this file since 8738 was 8738, checked in by dancopsey, 6 years ago
Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) up to revision 8588
File size: 32.0 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	! MV MP 2016
100	IF ( ln_pnd ) THEN ! Melt pond
101	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
102	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
103	ENDIF
104	! END MP 2016
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	!
127	tm_i(:,:) = 0._wp
128	sm_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	sm_i(:,:) = sm_i(:,:) + r1_nlay_i * sz_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
133	END DO
134	END DO
135	!
136	! ! put rt0 where there is no ice
137	WHERE( at_i(:,:)<=epsi20 )
138	tm_su(:,:) = rt0
139	tm_si(:,:) = rt0
140	tm_i (:,:) = rt0
141	END WHERE
142
143	DEALLOCATE( z1_at_i , z1_vt_i )
144	ENDIF
145	!
146	END SUBROUTINE ice_var_agg
147
148
149	SUBROUTINE ice_var_glo2eqv
150	!!-------------------------------------------------------------------
151	!! * ROUTINE ice_var_glo2eqv *
152	!!
153	!! ** Purpose : computes equivalent variables as function of
154	!! global variables, i.e. it turns VGLO into VEQV
155	!!-------------------------------------------------------------------
156	INTEGER :: ji, jj, jk, jl ! dummy loop indices
157	REAL(wp) :: ze_i ! local scalars
158	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
159	REAL(wp) :: zhmax, z1_zhmax ! - -
160	REAL(wp) :: zlay_i, zlay_s ! - -
161	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
162	!!-------------------------------------------------------------------
163
164	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
165	!! ice area and ice volume ?
166	!! the idea is to be able to define one for all at the begining of this routine
167	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
168
169	!-------------------------------------------------------
170	! Ice thickness, snow thickness, ice salinity, ice age
171	!-------------------------------------------------------
172	! !--- inverse of the ice area
173	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
174	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
175	END WHERE
176	!
177	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
178	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
179	END WHERE
180	!
181	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness
182
183	zhmax = hi_max(jpl)
184	z1_zhmax = 1._wp / hi_max(jpl)
185	WHERE( h_i(:,:,jpl) > zhmax ) !--- bound h_i by hi_max (i.e. 99 m) with associated update of ice area
186	h_i (:,:,jpl) = zhmax
187	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
188	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) ! NB: v_i always /=0 as h_i > hi_max
189	END WHERE
190
191	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness
192
193	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age
194
195	IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere)
196	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
197	ELSEWHERE ; s_i(:,:,:) = rn_simin
198	END WHERE
199	ENDIF
200
201	CALL ice_var_salprof ! salinity profile
202
203	!-------------------
204	! Ice temperature [K] (with a minimum value (rt0 - 100.))
205	!-------------------
206	zlay_i = REAL( nlay_i , wp ) ! number of layers
207	DO jl = 1, jpl
208	DO jk = 1, nlay_i
209	DO jj = 1, jpj
210	DO ji = 1, jpi
211	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
212	!
213	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]
214	ztmelts = - sz_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C]
215	! Conversion q(S,T) -> T (second order equation)
216	zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
217	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
218	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
219	!
220	ELSE !--- no ice
221	t_i(ji,jj,jk,jl) = rt0
222	ENDIF
223	END DO
224	END DO
225	END DO
226	END DO
227
228	!--------------------
229	! Snow temperature [K] (with a minimum value (rt0 - 100.))
230	!--------------------
231	zlay_s = REAL( nlay_s , wp )
232	DO jk = 1, nlay_s
233	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
234	t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
235	ELSEWHERE !--- no ice
236	t_s(:,:,jk,:) = rt0
237	END WHERE
238	END DO
239
240	! integrated values
241	vt_i (:,:) = SUM( v_i, dim=3 )
242	vt_s (:,:) = SUM( v_s, dim=3 )
243	at_i (:,:) = SUM( a_i, dim=3 )
244
245	! MV MP 2016
246	! probably should resum for melt ponds ???
247
248	!
249	END SUBROUTINE ice_var_glo2eqv
250
251
252	SUBROUTINE ice_var_eqv2glo
253	!!-------------------------------------------------------------------
254	!! * ROUTINE ice_var_eqv2glo *
255	!!
256	!! ** Purpose : computes global variables as function of
257	!! equivalent variables, i.e. it turns VEQV into VGLO
258	!!-------------------------------------------------------------------
259	!
260	v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:)
261	v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:)
262	sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:)
263	!
264	END SUBROUTINE ice_var_eqv2glo
265
266
267	SUBROUTINE ice_var_salprof
268	!!-------------------------------------------------------------------
269	!! * ROUTINE ice_var_salprof *
270	!!
271	!! ** Purpose : computes salinity profile in function of bulk salinity
272	!!
273	!! ** Method : If bulk salinity greater than zsi1,
274	!! the profile is assumed to be constant (S_inf)
275	!! If bulk salinity lower than zsi0,
276	!! the profile is linear with 0 at the surface (S_zero)
277	!! If it is between zsi0 and zsi1, it is a
278	!! alpha-weighted linear combination of s_inf and s_zero
279	!!
280	!! ** References : Vancoppenolle et al., 2007
281	!!-------------------------------------------------------------------
282	INTEGER :: ji, jj, jk, jl ! dummy loop index
283	REAL(wp) :: zsal, z1_dS
284	REAL(wp) :: zargtemp , zs0, zs
285	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
286	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
287	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
288	!!-------------------------------------------------------------------
289
290	!!gm Question: Remove the option 3 ? How many years since it last use ?
291
292	SELECT CASE ( nn_icesal )
293	!
294	! !---------------------------------------!
295	CASE( 1 ) ! constant salinity in time and space !
296	! !---------------------------------------!
297	sz_i(:,:,:,:) = rn_icesal
298	s_i(:,:,:) = rn_icesal
299	!
300	! !---------------------------------------------!
301	CASE( 2 ) ! time varying salinity with linear profile !
302	! !---------------------------------------------!
303	!
304	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
305	!
306	DO jl = 1, jpl
307	DO jk = 1, nlay_i
308	sz_i(:,:,jk,jl) = s_i(:,:,jl)
309	END DO
310	END DO
311	! ! Slope of the linear profile
312	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
313	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
314	END WHERE
315	!
316	z1_dS = 1._wp / ( zsi1 - zsi0 )
317	DO jl = 1, jpl
318	DO jj = 1, jpj
319	DO ji = 1, jpi
320	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
321	! ! force a constant profile when SSS too low (Baltic Sea)
322	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
323	END DO
324	END DO
325	END DO
326
327	! Computation of the profile
328	DO jl = 1, jpl
329	DO jk = 1, nlay_i
330	DO jj = 1, jpj
331	DO ji = 1, jpi
332	! ! linear profile with 0 surface value
333	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
334	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
335	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
336	END DO
337	END DO
338	END DO
339	END DO
340	!
341	DEALLOCATE( z_slope_s , zalpha )
342	!
343	! !-------------------------------------------!
344	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
345	! !-------------------------------------------! (mean = 2.30)
346	!
347	s_i(:,:,:) = 2.30_wp
348	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
349	!! a array S_prof(1:nlay_i) containing the calculation and just do:
350	! DO jk = 1, nlay_i
351	! sz_i(:,:,jk,:) = S_prof(jk)
352	! END DO
353	!!gm end
354	!
355	DO jl = 1, jpl
356	DO jk = 1, nlay_i
357	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
358	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
359	END DO
360	END DO
361	!
362	END SELECT
363	!
364	END SUBROUTINE ice_var_salprof
365
366	SUBROUTINE ice_var_salprof1d
367	!!-------------------------------------------------------------------
368	!! * ROUTINE ice_var_salprof1d *
369	!!
370	!! ** Purpose : 1d computation of the sea ice salinity profile
371	!! Works with 1d vectors and is used by thermodynamic modules
372	!!-------------------------------------------------------------------
373	INTEGER :: ji, jk ! dummy loop indices
374	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
375	REAL(wp) :: zs, zs0 ! - -
376	!
377	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
378	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
379	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
380	!!-------------------------------------------------------------------
381	!
382	SELECT CASE ( nn_icesal )
383	!
384	! !---------------------------------------!
385	CASE( 1 ) ! constant salinity in time and space !
386	! !---------------------------------------!
387	sz_i_1d(1:npti,:) = rn_icesal
388	!
389	! !---------------------------------------------!
390	CASE( 2 ) ! time varying salinity with linear profile !
391	! !---------------------------------------------!
392	!
393	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
394	!
395	! ! Slope of the linear profile
396	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
397	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
398	END WHERE
399
400	z1_dS = 1._wp / ( zsi1 - zsi0 )
401	DO ji = 1, npti
402	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
403	! ! force a constant profile when SSS too low (Baltic Sea)
404	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
405	END DO
406	!
407	! Computation of the profile
408	DO jk = 1, nlay_i
409	DO ji = 1, npti
410	! ! linear profile with 0 surface value
411	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
412	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
413	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
414	END DO
415	END DO
416	!
417	DEALLOCATE( z_slope_s, zalpha )
418
419	! !-------------------------------------------!
420	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
421	! !-------------------------------------------! (mean = 2.30)
422	!
423	s_i_1d(1:npti) = 2.30_wp
424	!
425	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
426	DO jk = 1, nlay_i
427	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
428	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
429	DO ji = 1, npti
430	sz_i_1d(ji,jk) = zsal
431	END DO
432	END DO
433	!
434	END SELECT
435	!
436	END SUBROUTINE ice_var_salprof1d
437
438
439	SUBROUTINE ice_var_zapsmall
440	!!-------------------------------------------------------------------
441	!! * ROUTINE ice_var_zapsmall *
442	!!
443	!! ** Purpose : Remove too small sea ice areas and correct fluxes
444	!!-------------------------------------------------------------------
445	INTEGER :: ji, jj, jl, jk ! dummy loop indices
446	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
447	!!-------------------------------------------------------------------
448	!
449	DO jl = 1, jpl !== loop over the categories ==!
450	!
451	!-----------------------------------------------------------------
452	! Zap ice energy and use ocean heat to melt ice
453	!-----------------------------------------------------------------
454	WHERE( a_i(:,:,jl) > epsi20 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
455	ELSEWHERE ; h_i(:,:,jl) = 0._wp
456	END WHERE
457	!
458	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
459	ELSEWHERE ; zswitch(:,:) = 1._wp
460	END WHERE
461
462	DO jk = 1, nlay_i
463	DO jj = 1 , jpj
464	DO ji = 1 , jpi
465	! update exchanges with ocean
466	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
467	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
468	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
469	END DO
470	END DO
471	END DO
472
473	DO jj = 1 , jpj
474	DO ji = 1 , jpi
475	! update exchanges with ocean
476	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_rdtice
477	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice
478	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice
479	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
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	END DO
501	END DO
502
503	IF( ln_pnd ) THEN
504	DO jj = 1 , jpj
505	DO ji = 1 , jpi
506	IF( ln_pnd_fw ) &
507	& wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
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	END DO
511	END DO
512	ENDIF
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 :: ijpij, 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	ijpij = SIZE( zhti , 1 )
571	zh_i(1:ijpij,1:jpl) = 0._wp
572	zh_s(1:ijpij,1:jpl) = 0._wp
573	za_i (1:ijpij,1:jpl) = 0._wp
574
575	DO ji = 1, ijpij
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	! clem: 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	ENDDO
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, ijpij
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_bv
674	!!-------------------------------------------------------------------
675	!! * ROUTINE ice_var_bv *
676	!!
677	!! ** Purpose : computes mean brine volume (%) in sea ice
678	!!
679	!! ** Method : e = - 0.054 * S (ppt) / T (C)
680	!!
681	!! References : Vancoppenolle et al., JGR, 2007
682	!!-------------------------------------------------------------------
683	INTEGER :: ji, jj, jk, jl ! dummy loop indices
684	!!-------------------------------------------------------------------
685	!
686	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
687	!! instead of setting everything to zero as just below
688	bv_i (:,:,:) = 0._wp
689	DO jl = 1, jpl
690	DO jk = 1, nlay_i
691	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
692	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
693	END WHERE
694	END DO
695	END DO
696	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
697	ELSEWHERE ; bvm_i(:,:) = 0._wp
698	END WHERE
699	!
700	END SUBROUTINE ice_var_bv
701
702
703	#else
704	!!----------------------------------------------------------------------
705	!! Default option Dummy module NO ESIM sea-ice model
706	!!----------------------------------------------------------------------
707	#endif
708
709	!!======================================================================
710	END MODULE icevar

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

Download in other formats: