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

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

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

Download in other formats: