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

Last change on this file since 8531 was 8522, checked in by clem, 7 years ago
changes in style - part6 - ice diffusion (hope the scheme still converges)
File size: 32.4 KB

Line
1	MODULE icevar
2	!!======================================================================
3	!! * MODULE icevar *
4	!! 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	!! - smv_i(jpi,jpj,jpl)
15	!! - oa_i (jpi,jpj,jpl)
16	!! VEQV : equivalent variables sometimes used in the model
17	!! - ht_i(jpi,jpj,jpl)
18	!! - ht_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	!! - smt_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' LIM3 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_bv : brine volume
45	!! ice_var_salprof1d : salinity profile in the ice 1D
46	!! ice_var_zapsmall : remove very small area and volume
47	!! ice_var_itd : convert 1-cat to multiple cat
48	!!----------------------------------------------------------------------
49	USE par_oce ! ocean parameters
50	USE phycst ! physical constants (ocean directory)
51	USE sbc_oce , ONLY : sss_m
52	USE ice ! ice variables
53	USE ice1D ! ice variables (thermodynamics)
54	!
55	USE in_out_manager ! I/O manager
56	USE lib_mpp ! MPP library
57	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
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_bv
67	PUBLIC ice_var_salprof1d
68	PUBLIC ice_var_zapsmall
69	PUBLIC ice_var_itd
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	htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
120	htm_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	smt_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	smt_i(:,:) = smt_i(:,:) + r1_nlay_i * s_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
133	END DO
134	END DO
135	!
136	!!gm QUESTION 1 : why salinity is named smt_i and not just sm_i ? since the 4D field is named s_i. (NB for temp: tm_i, t_i)
137	!
138	DEALLOCATE( z1_at_i , z1_vt_i )
139	ENDIF
140	!
141	END SUBROUTINE ice_var_agg
142
143
144	SUBROUTINE ice_var_glo2eqv
145	!!------------------------------------------------------------------
146	!! * ROUTINE ice_var_glo2eqv *
147	!!
148	!! ** Purpose : computes equivalent variables as function of
149	!! global variables, i.e. it turns VGLO into VEQV
150	!!------------------------------------------------------------------
151	INTEGER :: ji, jj, jk, jl ! dummy loop indices
152	REAL(wp) :: ze_i ! local scalars
153	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
154	REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - -
155	REAL(wp) :: zlay_i, zlay_s ! - -
156	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
157	!!------------------------------------------------------------------
158
159	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
160	!! ice area and ice volume ?
161	!! the idea is to be able to define one for all at the begining of this routine
162	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
163
164	!-------------------------------------------------------
165	! Ice thickness, snow thickness, ice salinity, ice age
166	!-------------------------------------------------------
167	! !--- inverse of the ice area
168	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
169	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
170	END WHERE
171	!
172	ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness
173
174	zhmax = hi_max(jpl)
175	z1_zhmax = 1._wp / hi_max(jpl)
176	WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
177	ht_i (:,:,jpl) = zhmax
178	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
179	z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max
180	END WHERE
181
182	ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness
183
184	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age
185
186	IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere)
187	WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )
188	ELSEWHERE ; sm_i(:,:,:) = rn_simin
189	END WHERE
190	ENDIF
191
192	CALL ice_var_salprof ! salinity profile
193
194	!-------------------
195	! Ice temperature [K] (with a minimum value (rt0 - 100.))
196	!-------------------
197	zlay_i = REAL( nlay_i , wp ) ! number of layers
198	DO jl = 1, jpl
199	DO jk = 1, nlay_i
200	DO jj = 1, jpj
201	DO ji = 1, jpi
202	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
203	!
204	ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
205	ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C]
206	! Conversion q(S,T) -> T (second order equation)
207	zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
208	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
209	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
210	!
211	ELSE !--- no ice
212	t_i(ji,jj,jk,jl) = rt0
213	ENDIF
214	END DO
215	END DO
216	END DO
217	END DO
218
219	!--------------------
220	! Snow temperature [K] (with a minimum value (rt0 - 100.))
221	!--------------------
222	zlay_s = REAL( nlay_s , wp )
223	DO jk = 1, nlay_s
224	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
225	t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
226	ELSEWHERE !--- no ice
227	t_s(:,:,jk,:) = rt0
228	END WHERE
229	END DO
230
231	! integrated values
232	vt_i (:,:) = SUM( v_i, dim=3 )
233	vt_s (:,:) = SUM( v_s, dim=3 )
234	at_i (:,:) = SUM( a_i, dim=3 )
235
236	! MV MP 2016
237	! probably should resum for melt ponds ???
238
239	!
240	END SUBROUTINE ice_var_glo2eqv
241
242
243	SUBROUTINE ice_var_eqv2glo
244	!!------------------------------------------------------------------
245	!! * ROUTINE ice_var_eqv2glo *
246	!!
247	!! ** Purpose : computes global variables as function of
248	!! equivalent variables, i.e. it turns VEQV into VGLO
249	!!------------------------------------------------------------------
250	!
251	v_i (:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
252	v_s (:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
253	smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
254	!
255	END SUBROUTINE ice_var_eqv2glo
256
257
258	SUBROUTINE ice_var_salprof
259	!!------------------------------------------------------------------
260	!! * ROUTINE ice_var_salprof *
261	!!
262	!! ** Purpose : computes salinity profile in function of bulk salinity
263	!!
264	!! ** Method : If bulk salinity greater than zsi1,
265	!! the profile is assumed to be constant (S_inf)
266	!! If bulk salinity lower than zsi0,
267	!! the profile is linear with 0 at the surface (S_zero)
268	!! If it is between zsi0 and zsi1, it is a
269	!! alpha-weighted linear combination of s_inf and s_zero
270	!!
271	!! ** References : Vancoppenolle et al., 2007
272	!!------------------------------------------------------------------
273	INTEGER :: ji, jj, jk, jl ! dummy loop index
274	REAL(wp) :: zsal, z1_dS
275	REAL(wp) :: zargtemp , zs0, zs
276	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
277	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
278	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
279	!!------------------------------------------------------------------
280
281	!!gm Question: Remove the option 3 ? How many years since it last use ?
282
283	SELECT CASE ( nn_icesal )
284	!
285	! !---------------------------------------!
286	CASE( 1 ) ! constant salinity in time and space !
287	! !---------------------------------------!
288	s_i (:,:,:,:) = rn_icesal
289	sm_i(:,:,:) = rn_icesal
290	!
291	! !---------------------------------------------!
292	CASE( 2 ) ! time varying salinity with linear profile !
293	! !---------------------------------------------!
294	!
295	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
296	!
297	DO jk = 1, nlay_i
298	s_i(:,:,jk,:) = sm_i(:,:,:)
299	END DO
300	! ! Slope of the linear profile
301	WHERE( ht_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:)
302	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
303	END WHERE
304	!
305	z1_dS = 1._wp / ( zsi1 - zsi0 )
306	DO jl = 1, jpl
307	DO jj = 1, jpj
308	DO ji = 1, jpi
309	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
310	! ! force a constant profile when SSS too low (Baltic Sea)
311	IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
312	END DO
313	END DO
314	END DO
315
316	! Computation of the profile
317	DO jl = 1, jpl
318	DO jk = 1, nlay_i
319	DO jj = 1, jpj
320	DO ji = 1, jpi
321	! ! linear profile with 0 surface value
322	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
323	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) ! weighting the profile
324	s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
325	END DO
326	END DO
327	END DO
328	END DO
329	!
330	DEALLOCATE( z_slope_s , zalpha )
331	!
332	! !-------------------------------------------!
333	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
334	! !-------------------------------------------! (mean = 2.30)
335	!
336	sm_i(:,:,:) = 2.30_wp
337	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
338	!! a array S_prof(1:nlay_i) containing the calculation and just do:
339	! DO jk = 1, nlay_i
340	! s_i(:,:,jk,:) = S_prof(jk)
341	! END DO
342	!!gm end
343	!
344	DO jl = 1, jpl
345	DO jk = 1, nlay_i
346	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
347	s_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
348	END DO
349	END DO
350	!
351	END SELECT
352	!
353	END SUBROUTINE ice_var_salprof
354
355
356	SUBROUTINE ice_var_bv
357	!!------------------------------------------------------------------
358	!! * ROUTINE ice_var_bv *
359	!!
360	!! ** Purpose : computes mean brine volume (%) in sea ice
361	!!
362	!! ** Method : e = - 0.054 * S (ppt) / T (C)
363	!!
364	!! References : Vancoppenolle et al., JGR, 2007
365	!!------------------------------------------------------------------
366	INTEGER :: ji, jj, jk, jl ! dummy loop indices
367	!!------------------------------------------------------------------
368	!
369	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
370	!! instead of setting everything to zero as just below
371	bv_i (:,:,:) = 0._wp
372	DO jl = 1, jpl
373	DO jk = 1, nlay_i
374	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
375	bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
376	END WHERE
377	END DO
378	END DO
379	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
380	ELSEWHERE ; bvm_i(:,:) = 0._wp
381	END WHERE
382	!
383	END SUBROUTINE ice_var_bv
384
385
386	SUBROUTINE ice_var_salprof1d
387	!!-------------------------------------------------------------------
388	!! * ROUTINE ice_var_salprof1d *
389	!!
390	!! ** Purpose : 1d computation of the sea ice salinity profile
391	!! Works with 1d vectors and is used by thermodynamic modules
392	!!-------------------------------------------------------------------
393	INTEGER :: ji, jk ! dummy loop indices
394	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
395	REAL(wp) :: zalpha, zs, zs0 ! - -
396	!
397	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s !
398	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
399	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
400	!!---------------------------------------------------------------------
401	!
402	SELECT CASE ( nn_icesal )
403	!
404	! !---------------------------------------!
405	CASE( 1 ) ! constant salinity in time and space !
406	! !---------------------------------------!
407	s_i_1d(:,:) = rn_icesal
408	!
409	! !---------------------------------------------!
410	CASE( 2 ) ! time varying salinity with linear profile !
411	! !---------------------------------------------!
412	!
413	ALLOCATE( z_slope_s(jpij) )
414	!
415	DO ji = 1, nidx ! Slope of the linear profile
416	rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
417	z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
418	END DO
419
420	z1_dS = 1._wp / ( zsi1 - zsi0 )
421	DO jk = 1, nlay_i
422	DO ji = 1, nidx
423	zalpha = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) )
424	IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp ! cst profile when SSS too low (Baltic Sea)
425	!
426	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i ! linear profile with 0 surface value
427	zs = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji) ! weighting the profile
428	s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) ! bounding salinity
429	END DO
430	END DO
431	!
432	DEALLOCATE( z_slope_s )
433
434	! !-------------------------------------------!
435	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
436	! !-------------------------------------------! (mean = 2.30)
437	!
438	sm_i_1d(:) = 2.30_wp
439	!
440	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
441	DO jk = 1, nlay_i
442	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
443	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
444	DO ji = 1, nidx
445	s_i_1d(ji,jk) = zsal
446	END DO
447	END DO
448	!
449	END SELECT
450	!
451	END SUBROUTINE ice_var_salprof1d
452
453
454	SUBROUTINE ice_var_zapsmall
455	!!-------------------------------------------------------------------
456	!! * ROUTINE ice_var_zapsmall *
457	!!
458	!! ** Purpose : Remove too small sea ice areas and correct fluxes
459	!!-------------------------------------------------------------------
460	INTEGER :: ji, jj, jl, jk ! dummy loop indices
461	REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp
462	!!-------------------------------------------------------------------
463	!
464	DO jl = 1, jpl !== loop over the categories ==!
465	!
466	!-----------------------------------------------------------------
467	! Zap ice energy and use ocean heat to melt ice
468	!-----------------------------------------------------------------
469	DO jk = 1, nlay_i
470	DO jj = 1 , jpj
471	DO ji = 1 , jpi
472	!!gm I think we can do better/faster : to be discussed
473	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
474	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
475	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
476	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
477	zei = e_i(ji,jj,jk,jl)
478	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
479	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
480	! update exchanges with ocean
481	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
482	END DO
483	END DO
484	END DO
485
486	DO jj = 1 , jpj
487	DO ji = 1 , jpi
488	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
489	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
490	rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &
491	& / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
492	zsal = smv_i(ji,jj, jl)
493	zvi = v_i (ji,jj, jl)
494	zvs = v_s (ji,jj, jl)
495	zes = e_s (ji,jj,1,jl)
496	IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl)
497	!-----------------------------------------------------------------
498	! Zap snow energy
499	!-----------------------------------------------------------------
500	t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
501	e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
502
503	!-----------------------------------------------------------------
504	! zap ice and snow volume, add water and salt to ocean
505	!-----------------------------------------------------------------
506	ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
507	a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch
508	v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch
509	v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch
510	t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
511	oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
512	smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
513
514	ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
515	ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
516
517	! MV MP 2016
518	IF ( ln_pnd ) THEN
519	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
520	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
521	IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice
522	ENDIF
523	! END MV MP 2016
524
525	! update exchanges with ocean
526	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
527	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
528	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
529	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
530	END DO
531	END DO
532	END DO
533
534	! to be sure that at_i is the sum of a_i(jl)
535	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
536	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
537
538	! open water = 1 if at_i=0
539	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
540	!
541	END SUBROUTINE ice_var_zapsmall
542
543
544	SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
545	!!------------------------------------------------------------------
546	!! * ROUTINE ice_var_itd *
547	!!
548	!! ** Purpose : converting 1-cat ice to multiple ice categories
549	!!
550	!! ice thickness distribution follows a gaussian law
551	!! around the concentration of the most likely ice thickness
552	!! (similar as iceistate.F90)
553	!!
554	!! ** Method: Iterative procedure
555	!!
556	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
557	!!
558	!! 2) Check whether the distribution conserves area and volume, positivity and
559	!! category boundaries
560	!!
561	!! 3) If not (input ice is too thin), the last category is empty and
562	!! the number of categories is reduced (jpl-1)
563	!!
564	!! 4) Iterate until ok (SUM(itest(:) = 4)
565	!!
566	!! ** Arguments : zhti: 1-cat ice thickness
567	!! zhts: 1-cat snow depth
568	!! zai : 1-cat ice concentration
569	!!
570	!! ** Output : jpl-cat
571	!!
572	!! (Example of application: BDY forcings when input are cell averaged)
573	!!-------------------------------------------------------------------
574	INTEGER :: ji, jk, jl ! dummy loop indices
575	INTEGER :: ijpij, i_fill, jl0
576	REAL(wp) :: zarg, zV, zconv, zdh, zdv
577	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables
578	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables
579	INTEGER , DIMENSION(4) :: itest
580	!!-------------------------------------------------------------------
581
582	!--------------------------------------------------------------------
583	! initialisation of variables
584	!--------------------------------------------------------------------
585	ijpij = SIZE( zhti , 1 )
586	zht_i(1:ijpij,1:jpl) = 0._wp
587	zht_s(1:ijpij,1:jpl) = 0._wp
588	za_i (1:ijpij,1:jpl) = 0._wp
589
590	! ----------------------------------------
591	! distribution over the jpl ice categories
592	! ----------------------------------------
593	DO ji = 1, ijpij
594
595	IF( zhti(ji) > 0._wp ) THEN
596
597	! find which category (jl0) the input ice thickness falls into
598	jl0 = jpl
599	DO jl = 1, jpl
600	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
601	jl0 = jl
602	CYCLE
603	ENDIF
604	END DO
605
606	! initialisation of tests
607	itest(:) = 0
608
609	i_fill = jpl + 1 !====================================
610	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
611	! iteration !====================================
612	i_fill = i_fill - 1
613
614	! initialisation of ice variables for each try
615	zht_i(ji,1:jpl) = 0._wp
616	za_i (ji,1:jpl) = 0._wp
617	itest(:) = 0
618
619	! *** case very thin ice: fill only category 1
620	IF ( i_fill == 1 ) THEN
621	zht_i(ji,1) = zhti(ji)
622	za_i (ji,1) = zai (ji)
623
624	! *** case ice is thicker: fill categories >1
625	ELSE
626
627	! Fill ice thicknesses in the (i_fill-1) cat by hmean
628	DO jl = 1, i_fill - 1
629	zht_i(ji,jl) = hi_mean(jl)
630	END DO
631
632	! Concentrations in the (i_fill-1) categories
633	za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
634	DO jl = 1, i_fill - 1
635	IF ( jl /= jl0 ) THEN
636	zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
637	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
638	ENDIF
639	END DO
640
641	! Concentration in the last (i_fill) category
642	za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
643
644	! Ice thickness in the last (i_fill) category
645	zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
646	zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
647
648	! clem: correction if concentration of upper cat is greater than lower cat
649	! (it should be a gaussian around jl0 but sometimes it is not)
650	IF ( jl0 /= jpl ) THEN
651	DO jl = jpl, jl0+1, -1
652	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
653	zdv = zht_i(ji,jl) * za_i(ji,jl)
654	zht_i(ji,jl ) = 0._wp
655	za_i (ji,jl ) = 0._wp
656	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
657	END IF
658	ENDDO
659	ENDIF
660
661	ENDIF ! case ice is thick or thin
662
663	!---------------------
664	! Compatibility tests
665	!---------------------
666	! Test 1: area conservation
667	zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
668	IF ( zconv < epsi06 ) itest(1) = 1
669
670	! Test 2: volume conservation
671	zconv = ABS( zhti(ji)zai(ji) - SUM( za_i(ji,1:jpl)zht_i(ji,1:jpl) ) )
672	IF ( zconv < epsi06 ) itest(2) = 1
673
674	! Test 3: thickness of the last category is in-bounds ?
675	IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
676
677	! Test 4: positivity of ice concentrations
678	itest(4) = 1
679	DO jl = 1, i_fill
680	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
681	END DO
682	! !============================
683	END DO ! end iteration on categories
684	! !============================
685	ENDIF ! if zhti > 0
686	END DO ! i loop
687
688	! ------------------------------------------------
689	! Adding Snow in each category where za_i is not 0
690	! ------------------------------------------------
691	DO jl = 1, jpl
692	DO ji = 1, ijpij
693	IF( za_i(ji,jl) > 0._wp ) THEN
694	zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
695	! In case snow load is in excess that would lead to transformation from snow to ice
696	! Then, transfer the snow excess into the ice (different from icethd_dh)
697	zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )
698	! recompute ht_i, ht_s avoiding out of bounds values
699	zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
700	zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
701	ENDIF
702	END DO
703	END DO
704	!
705	END SUBROUTINE ice_var_itd
706
707	#else
708	!!----------------------------------------------------------------------
709	!! Default option Dummy module NO LIM3 sea-ice model
710	!!----------------------------------------------------------------------
711	#endif
712
713	!!======================================================================
714	END MODULE icevar

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

Download in other formats: