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

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

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

Download in other formats: