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.

limthd_dh.F90 in branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago
dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 48.4 KB

Line
1	MODULE limthd_dh
2	#if defined key_lim3
3	!!======================================================================
4	!! * MODULE limthd_dh *
5	!! thermodynamic growth and decay of the ice
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice
10	!! * Modules used
11
12	USE par_oce ! ocean parameters
13	USE phycst ! physical constants (OCE directory)
14	USE ice_oce ! ice variables
15	USE thd_ice
16	USE iceini
17	USE limistate
18	USE in_out_manager
19	USE ice
20	USE par_ice
21	USE limicepoints
22
23	IMPLICIT NONE
24	PRIVATE
25
26	!! * Routine accessibility
27	PUBLIC lim_thd_dh ! called by lim_thd
28
29	!! * Module variables
30	REAL(wp) :: & ! constant values
31	epsi20 = 1e-20 , &
32	epsi13 = 1e-13 , &
33	epsi16 = 1e-16 , &
34	zzero = 0.e0 , &
35	zone = 1.e0
36
37	!!----------------------------------------------------------------------
38	!! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005)
39	!!----------------------------------------------------------------------
40
41	CONTAINS
42
43	SUBROUTINE lim_thd_dh(kideb,kiut,jl)
44	!!------------------------------------------------------------------
45	!! * ROUTINE lim_thd_dh *
46	!!------------------------------------------------------------------
47	!! ** Purpose :
48	!! This routine determines variations of ice and snow thicknesses.
49	!! ** Method :
50	!! Ice/Snow surface melting arises from imbalance in surface fluxes
51	!! Bottom accretion/ablation arises from flux budget
52	!! Snow thickness can increase by precipitation and decrease by
53	!! sublimation
54	!! If snow load excesses Archmiede limit, snow-ice is formed by
55	!! the flooding of sea-water in the snow
56	!! ** Steps
57	!! 1) Compute available flux of heat for surface ablation
58	!! 2) Compute snow and sea ice enthalpies
59	!! 3) Surface ablation and sublimation
60	!! 4) Bottom accretion/ablation
61	!! 5) Case of Total ablation
62	!! 6) Snow ice formation
63	!!
64	!! ** Arguments
65	!!
66	!! ** Inputs / Outputs
67	!!
68	!! ** External
69	!!
70	!! ** References : Bitz and Lipscomb, JGR 99
71	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
72	!! Vancoppenolle Fichefet and Bitz, GRL 2005
73	!!
74	!! ** History :
75	!! original code 01-04 (LIM)
76	!! original routine
77	!! (05-2003) Martin Vancoppenolle, Louvain-La-Neuve, Belgium
78	!! (06/07-2005) Martin Vancoppenolle for the 3D version
79	!!
80	!!------------------------------------------------------------------
81	!! * Arguments
82	INTEGER , INTENT (IN) :: &
83	kideb , & !: Start point on which the the computation is applied
84	kiut , & !: End point on which the the computation is applied
85	jl !: Thickness cateogry number
86
87	!! * Local variables
88	INTEGER :: &
89	ji , & !: space index
90	jk , & !: ice layer index
91	isnow , & !: switch for presence (1) or absence (0) of snow
92	index , & !: ice point index (limicepoints)
93	zji , & !: 2D corresponding indices to ji
94	zjj , &
95	isnowic , & !: snow ice formation not
96	i_ice_switch , & !: ice thickness above a certain treshold or not
97	iter
98
99	REAL(wp) :: &
100	zhsnew , & !: new snow thickness
101	zihgnew , & !: switch for total ablation
102	ztmelts , & !: melting point
103	zhn , &
104	zdhcf , &
105	zdhbf , &
106	zhni , &
107	zhnfi , &
108	zihg , &
109	zdhgm , &
110	zdhnm , &
111	zhnnew , &
112	zdfseqv , &
113	zeps = 1.0e-13, &
114	zhisn , &
115	zfracs , & !: fractionation coefficient for bottom salt
116	!: entrapment
117	zds , & !: increment of bottom ice salinity
118	zoldsinew , & !: dummy argument for error diagnosis
119	zcoeff , & !: dummy argument for snowfall partitioning
120	!: over ice and leads
121	zjiindex , & !: zdhs_temp, &
122	zsm_snowice, & !: snow-ice salinity
123	zswi1 , & !: switch for computation of bottom salinity
124	zswi12 , & !: switch for computation of bottom salinity
125	zswi2 , & !: switch for computation of bottom salinity
126	zgrr , & !: bottom growth rate
127	zihic , & !:
128	ztform !: bottom formation temperature
129
130	REAL(wp) , DIMENSION(jpij) :: &
131	zh_i , & ! ice layer thickness
132	zh_s , & ! snow layer thickness
133	ztfs , & ! melting point
134	zhsold , & ! old snow thickness
135	zqprec , & !: energy of fallen snow
136	zqfont_su , & ! incoming, remaining surface energy
137	zqfont_bo ! incoming, bottom energy
138
139	REAL(wp) , DIMENSION(jpij) :: &
140	z_f_surf, & ! surface heat for ablation
141	zhgnew ! new ice thickness
142
143	REAL(wp), DIMENSION(jpij) :: &
144	zdh_s_mel , & ! snow melt
145	zdh_s_pre , & ! snow precipitation
146	zdh_s_sub , & ! snow sublimation
147	zfsalt_melt ! salt flux due to ice melt
148
149	REAL(wp) , DIMENSION(jpij,jkmax) :: &
150	zdeltah
151
152	! Pathological cases
153	REAL(wp), DIMENSION(jpij) :: &
154	zfdt_init , & !: total incoming heat for ice melt
155	zfdt_final , & !: total remaing heat for ice melt
156	zqt_i , & !: total ice heat content
157	zqt_s , & !: total snow heat content
158	zqt_dummy !: dummy heat content
159
160	REAL(wp), DIMENSION(jpij,jkmax) :: &
161	zqt_i_lay , & !: total ice heat content
162	zqt_s_lay !: total snow heat content
163
164	! Heat conservation
165	REAL(wp), DIMENSION(jpij) :: &
166	zfbase, &
167	zdq_i
168
169	INTEGER, DIMENSION(jpij) :: &
170	innermelt
171
172	REAL(wp) :: &
173	meance_dh
174
175	INTEGER :: &
176	num_iter_max, &
177	numce_dh
178
179	!!-----------------------------------------------------------------------------
180
181	WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation'
182	WRITE(numout,*) '~~~~~~~~~'
183	zfsalt_melt(:) = 0.0
184	ftotal_fin(:) = 0.0
185	zfdt_init(:) = 0.0
186	zfdt_final(:) = 0.0
187
188	DO ji = kideb, kiut
189	old_ht_i_b(ji) = ht_i_b(ji)
190	old_ht_s_b(ji) = ht_s_b(ji)
191	END DO
192	!
193	!------------------------------------------------------------------------------!
194	! 1) Calculate available heat for surface ablation !
195	!------------------------------------------------------------------------------!
196	!
197	DO ji = kideb,kiut
198	isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
199	ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt
200	z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * &
201	qsr_ice_1d(ji) - fc_su(ji)
202	z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * &
203	MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
204	zfdt_init(ji) = ( z_f_surf(ji) + &
205	MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) &
206	* rdt_ice
207	END DO ! ji
208
209	zqfont_su(:) = 0.0
210	zqfont_bo(:) = 0.0
211	dsm_i_se_1d(:) = 0.0
212	dsm_i_si_1d(:) = 0.0
213	! ! heat conservation
214	! IF ( con_i ) THEN
215	! IF (jiindex_1d .GT.0) WRITE(numout,*) ' zf_surf : ', z_f_surf(jiindex_1d)
216	! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qnsr_ice : ', qnsr_ice_1d(jiindex_1d)
217	! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qsr_ice : ', qsr_ice_1d(jiindex_1d)
218	! IF (jiindex_1d .GT.0) WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)
219	! IF (jiindex_1d .GT.0) WRITE(numout,*) ' i0 : ', i0(jiindex_1d)
220	! ENDIF
221
222	!------------------------------------------------------------------------------!
223	! 2) Computing layer thicknesses and snow and sea-ice enthalpies. !
224	!------------------------------------------------------------------------------!
225
226	!-----------------
227	! Layer thickness
228	!-----------------
229	DO ji = kideb,kiut
230	zh_i(ji) = ht_i_b(ji) / nlay_i
231	zh_s(ji) = ht_s_b(ji) / nlay_s
232	END DO
233
234	!--------------------------------------
235	! Snow energy of melting, heat content
236	!--------------------------------------
237	zqt_s(:) = 0.0
238	DO jk = 1, nlay_s
239	DO ji = kideb,kiut
240	! this line could be removed
241	! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
242	! heat conservation
243	zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
244	END DO
245	END DO
246
247	! ! heat conservation
248	! qt_s_in(:,jl) = zqt_s(:)
249	! IF (jiindex_1d.GT.0) &
250	! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / &
251	! rdt_ice
252
253	! IF (jiindex_1d.GT.0) THEN
254	! WRITE(numout,*) ' Vertical profile : '
255	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
256	! WRITE(numout,) ' qm0 : ', q_s_b(jiindex_1d,1:nlay_s) ht_s_b(jiindex_1d) / &
257	! rdt_ice
258	! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
259	! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
260	! WRITE(numout,*) ' zthick0: ', ht_s_b(jiindex_1d)
261	! ENDIF
262
263	!-------------------------------------
264	! Ice energy of melting, heat content
265	!-------------------------------------
266	zqt_i(:) = 0.0
267	DO jk = 1, nlay_i
268	DO ji = kideb,kiut
269	! Melting point in K
270	! ztmelts = - tmut * s_i_b(ji,jk) + rtt
271	! this line could be removed
272	! q_i_b(ji,jk) = rhoic * &
273	! ( cpic * ( ztmelts - t_i_b(ji,jk) ) &
274	! + lfus * ( 1.0 - (ztmelts-rtt) / &
275	! MIN( ( t_i_b(ji,jk) - rtt ) , - zeps ) ) &
276	! - rcp * ( ztmelts-rtt ) )
277	zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
278	zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
279	END DO
280	END DO
281	! IF (jiindex_1d.GT.0) &
282	! WRITE(numout,*) 'jl : ', jl, ' zqti : ', zqt_i(jiindex_1d) / &
283	! rdt_ice
284
285	! IF (jiindex_1d.GT.0) THEN
286	! WRITE(numout,*) ' --- Vertical profile : '
287	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
288	! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
289	! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)
290	! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)
291	! WRITE(numout,*) ' qm0 : ', zqt_i_lay(jiindex_1d,1:nlay_i) / rdt_ice
292	! WRITE(numout,*) ' zthick0: ', ht_i_b(jiindex_1d) / nlay_i
293	! ENDIF
294
295	!
296	!------------------------------------------------------------------------------\|
297	! 3) Surface ablation and sublimation \|
298	!------------------------------------------------------------------------------\|
299	!
300	!---------------------
301	! Snow precips / melt
302	!---------------------
303	! Snow accumulation in one thermodynamic time step
304	! snowfall is partitionned between leads and ice
305	! if snow fall was uniform, a fraction (1-at_i) would fall into leads
306	! but because of the winds, more snow falls on leads than on sea ice
307	! and a greater fraction (1-at_i)^beta of the total mass of snow
308	! (beta < 1) falls in leads
309	! In reality, beta depends on wind speed,
310	! and should decrease with increasing wind speed but here, it is
311	! considered as a constant
312	! an average value is 0.66
313	! Martin Vancoppenolle, December 2006
314
315	! Snow fall
316	DO ji = kideb, kiut
317	zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji)
318	zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
319	END DO
320	zdh_s_mel(:) = 0.0
321
322	! Melt of fallen snow
323	DO ji = kideb, kiut
324
325	! tatm_ice is now in K
326	zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )
327	zqfont_su(ji) = z_f_surf(ji) * rdt_ice
328	zdeltah(ji,1) = MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
329	zqfont_su(ji) = MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * &
330	zqprec(ji)
331	zdeltah(ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
332	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)
333	! heat conservation
334	qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
335	zqt_s(ji) = zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji)
336
337	END DO
338
339	! IF (jiindex_1d .GT. 0) THEN
340	! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)
341	! WRITE(numout,) ' zqprec : ', zqprec(jiindex_1d)( zdh_s_pre(jiindex_1d) + zdh_s_mel(jiindex_1d) ) / rdt_ice
342	! WRITE(numout,*) ' tatm : ', tatm_ice_1d(jiindex_1d)
343	! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / &
344	! rdt_ice
345	! ENDIF
346
347	! updated total snow heat content
348	zqt_s(ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 )
349
350	! IF (jiindex_1D .GT. 0) &
351	! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / &
352	! rdt_ice
353
354	! Snow melt due to surface heat imbalance
355	DO jk = 1, nlay_s
356	DO ji = kideb, kiut
357	zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
358	zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * &
359	q_s_b(ji,jk)
360	! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &
361	! zqfont_su(jiindex_1d)
362	zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) )
363	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow
364	END DO
365	END DO
366
367	! !+++++
368	! IF (jiindex_1d .GT. 0 ) THEN
369	! WRITE(numout,*) ' dhs_pre :', zdh_s_pre(jiindex_1d)
370	! WRITE(numout,*) ' dhs_mel :', zdh_s_mel(jiindex_1d)
371	! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)
372	! ENDIF
373	! !+++++
374
375	! Apply snow melt to snow depth
376	DO ji = kideb, kiut
377	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji)
378	! Old and new snow depths
379	zhsold(ji) = ht_s_b(ji)
380	zhsnew = ht_s_b(ji) + dh_s_tot(ji)
381	! If snow is still present zhn = 1, else zhn = 0
382	zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
383	ht_s_b(ji) = MAX( zzero , zhsnew )
384	! Volume and mass variations of snow
385	dvsbq_1d(ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) &
386	- zdh_s_mel(ji) )
387	dvsbq_1d(ji) = MIN( zzero, dvsbq_1d(ji) )
388	rdmsnif_1d(ji) = rhosn*dvsbq_1d(ji)
389	END DO ! ji
390
391	! !+++++
392	! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' dhs_tot :', dh_s_tot(jiindex_1d)
393	! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' zqfont_su :', zqfont_su(jiindex_1d)
394	! !+++++
395
396	!----------------------
397	! Surface ice ablation
398	!----------------------
399	DO ji = kideb, kiut
400	! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &
401	! zqfont_su(jiindex_1d)
402	dh_i_surf(ji) = 0.0
403	! Heat conservation test
404	z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test
405	zdq_i(ji) = 0.0
406	END DO ! ji
407
408	DO jk = 1, nlay_i
409	DO ji = kideb, kiut
410	zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)
411	zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
412	q_i_b(ji,jk)
413	zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )
414	dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk)
415	zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * &
416	q_i_b(ji,jk) / rdt_ice
417	! Contribution to salt flux -> to change
418	zji = MOD( npb(ji) - 1, jpi ) + 1
419	zjj = ( npb(ji) - 1 ) / jpi + 1
420	zfsalt_melt(ji) = zfsalt_melt(ji) + &
421	! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &
422	( sss_io(zji,zjj) - sm_i_b(ji) ) * &
423	a_i_b(ji) * &
424	MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice
425	!+++++++++++
426	! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN
427	! WRITE(numout,*) ' surf abl '
428	! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) / ( sss_io(zji,zjj)+epsi16)
429	! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)
430	! ENDIF
431	!+++++++++++
432	END DO ! ji
433	END DO ! jk
434
435	!-------------------
436	! Conservation test
437	!-------------------
438	IF ( con_i ) THEN
439	numce_dh = 0
440	meance_dh = 0.0
441	DO ji = kideb, kiut
442
443	IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
444	numce_dh = numce_dh + 1
445	meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
446	ENDIF
447
448	IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. &
449	! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN
450	WRITE(numout,*) ' ALERTE heat loss for surface melt '
451	WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
452	WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
453	WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji)
454	WRITE(numout,*) ' zdq_i : ', zdq_i(ji)
455	WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
456	WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
457	WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
458	WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
459	WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
460	WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)
461	ENDIF
462
463	END DO ! ji
464
465	IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
466	WRITE(numout,*) ' Error report - Category : ', jl
467	WRITE(numout,*) ' ~~~~~~~~~~~~ '
468	WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
469	WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
470
471	ENDIF ! con_i
472
473	!------------------
474	! Snow sublimation
475	!------------------
476	! !++++++
477	! zqt_dummy(:) = 0.0
478	! DO jk = 1, nlay_s
479	! DO ji = kideb,kiut
480	! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
481	! ! heat conservation
482	! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * old_ht_s_b(ji) / nlay_s
483	! END DO
484	! END DO
485
486	! IF (jiindex_1d .GT. 0) &
487	! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &
488	! rdt_ice
489	!++++++
490	! !++++++
491	! zqt_dummy(:) = 0.0
492	! DO jk = 1, nlay_s
493	! DO ji = kideb,kiut
494	! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
495	! ! heat conservation
496	! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
497	! END DO
498	! END DO
499	! IF (jiindex_1d .GT. 0) &
500	! WRITE(numout,*) 'jl : ', jl, ' New zqts : ', zqt_dummy(jiindex_1d) / &
501	! rdt_ice
502	! !++++++
503
504	DO ji = kideb, kiut
505	! if qla is positive (upwards), heat goes to the atmosphere, therefore
506	! snow sublimates, if qla is negative (downwards), snow condensates
507	zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
508	dh_s_tot(ji) = dh_s_tot(ji) + zdh_s_sub(ji)
509	zdhcf = ht_s_b(ji) + zdh_s_sub(ji)
510	ht_s_b(ji) = MAX( zzero , zdhcf )
511	! we recompute dh_s_tot
512	! which may be different in case of total snow melt
513	dh_s_tot(ji) = ht_s_b(ji) - zhsold(ji)
514	! IF ( ht_s_b(ji) .LE. 0.0 ) THEN
515	! dh_s_tot(ji) = MAX( 0.0 , dh_s_tot(ji) )
516	! ENDIF
517
518	qt_s_in(ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
519	END DO !ji
520
521	! IF ( jiindex_1d .GT. 0 ) THEN
522	! WRITE(numout,*) ' dh_s_sub : ', zdh_s_sub(jiindex_1d)
523	! WRITE(numout,*) ' dhs_tot : ', dh_s_tot(jiindex_1d)
524	! ENDIF
525
526	!++++++
527	zqt_dummy(:) = 0.0
528	DO jk = 1, nlay_s
529	DO ji = kideb,kiut
530	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
531	! heat conservation
532	zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
533	END DO
534	END DO
535	! IF (jiindex_1d .GT. 0) &
536	! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &
537	! rdt_ice
538	! !++++++
539
540	! IF (jiindex_1d .GT. 0) THEN
541	! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
542	! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
543	! ENDIF
544
545	DO jk = 1, nlay_s !n
546	DO ji = kideb, kiut !n
547	! In case of disparition of the snow, we have to update the snow
548	! temperatures
549	zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
550	t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
551	q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
552	END DO
553	END DO
554
555	! IF (jiindex_1d .GT. 0) THEN
556	! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)
557	! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)
558	! ENDIF
559
560	!
561	!------------------------------------------------------------------------------!
562	! 4) Basal growth / melt !
563	!------------------------------------------------------------------------------!
564	!
565	! Ice basal growth / melt is given by the ratio of heat budget over basal
566	! ice heat content. Basal heat budget is given by the difference between
567	! the inner conductive flux (fc_bo_i), from the open water heat flux
568	! (qlbbqb) and the turbulent ocean flux (fbif).
569	! fc_bo_i is positive downwards
570	! fbif and qlbbq are positive to the ice
571
572	!---------------------------------------------
573	! Basal growth - salinity not varying in time
574	!---------------------------------------------
575	IF ( num_sal .NE. 2 ) THEN
576	DO ji = kideb, kiut
577	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN
578	s_i_new(ji) = sm_i_b(ji)
579	! Melting point in K
580	ztmelts = - tmut * s_i_new(ji) + rtt
581	! New ice heat content (Bitz and Lipscomb, 1999)
582	ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the
583	! Baltic
584	q_i_b(ji,nlay_i+1) = rhoic * &
585	( cpic * ( ztmelts - ztform ) &
586	+ lfus * ( 1.0 - ( ztmelts - rtt ) / &
587	( ztform - rtt ) ) &
588	- rcp * ( ztmelts-rtt ) )
589	! Basal growth rate = - F*dt / q
590	dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
591	qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
592	! IF ( ji .EQ. jiindex_1d ) THEN
593	! WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
594	! WRITE(numout,*) ' ztmelts : ', ztmelts
595	! WRITE(numout,*) ' t_bo : ', t_bo_b(ji)
596	! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,nlay_i+1)
597	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
598	! ENDIF
599	ENDIF ! heat budget
600	END DO ! ji
601	ENDIF ! num_sal
602
603	! IF ( jiindex_1d .GT. 0 ) THEN
604	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
605	! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,nlay_i+1)
606	! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d)
607	! WRITE(numout,*) ' fbif_1d : ', fbif_1d(jiindex_1d)
608	! WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(jiindex_1d)
609	! ENDIF
610
611	!-----------------------------------------
612	! Basal growth - salinity varying in time
613	!-----------------------------------------
614	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
615	! the growth rate (dh_i_bott) is function of the new ice
616	! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice
617	! salinity (snewice). snewice depends on dh_i_bott
618	! it converges quickly, so, no problem
619
620	! Initial value (tested 1D, can be anything between 1 and 20)
621	num_iter_max = 4
622	s_i_new(:) = 4.0
623
624	! Iterative procedure
625	DO iter = 1, num_iter_max
626	DO ji = kideb, kiut
627	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN
628	zji = MOD( npb(ji) - 1, jpi ) + 1
629	zjj = ( npb(ji) - 1 ) / jpi + 1
630	! Melting point in K
631	ztmelts = - tmut * s_i_new(ji) + rtt
632	! New ice heat content (Bitz and Lipscomb, 1999)
633	q_i_b(ji,nlay_i+1) = rhoic * &
634	( cpic * ( ztmelts - t_bo_b(ji) ) &
635	+ lfus * ( 1.0 - ( ztmelts - rtt ) / &
636	( t_bo_b(ji) - rtt ) ) &
637	- rcp * ( ztmelts-rtt ) )
638	! !+++++++
639	! IF (ji.EQ.jiindex_1d) THEN
640	! WRITE(numout,*) ' ztmelts : ', ztmelts
641	! WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji)
642	! WRITE(numout,*) ' q_i_b(ji,nlay_i+1) : ', q_i_b(ji,nlay_i+1)
643	! ENDIF
644	! !+++++++
645	! Bottom growth rate = - F*dt / q
646	dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
647	+ qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
648	! New ice salinity ( Cox and Weeks, JGR, 1988 )
649	! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7
650	! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
651	! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8
652	! zgrr = MAX( dh_i_bott(ji) / rdt_ice, 0.0 )
653	zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) )
654	zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )
655	zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
656	zswi1 = 1. - zswi2 * zswi12
657	! IF (ji.EQ.jiindex_1d) THEN
658	! WRITE(numout,*) ' zgrr : ', zgrr
659	! WRITE(numout,*) ' zswi2 : ', zswi2, ' zswi12 :', zswi12, ' zswi1 :', zswi1
660	! ENDIF
661	zfracs = zswi1 * 0.12 + &
662	zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
663	zswi2 * 0.26 / &
664	( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )
665	zds = zfracs*sss_io(zji,zjj) - s_i_new(ji)
666	s_i_new(ji) = zfracs * sss_io(zji,zjj)
667	ENDIF ! fc_bo_i
668	END DO ! ji
669	END DO ! iter
670
671	! Final values
672	DO ji = kideb, kiut
673	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN
674	! New ice salinity must not exceed 15 psu
675	zoldsinew = s_i_new(ji)
676	s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
677	! Metling point in K
678	ztmelts = - tmut * s_i_new(ji) + rtt
679	! New ice heat content (Bitz and Lipscomb, 1999)
680	q_i_b(ji,nlay_i+1) = rhoic * &
681	( cpic * ( ztmelts - t_bo_b(ji) ) &
682	+ lfus * ( 1.0 - ( ztmelts - rtt ) / &
683	( t_bo_b(ji) - rtt ) ) &
684	- rcp * ( ztmelts-rtt ) )
685	! Basal growth rate = - F*dt / q
686	! sometimes, growth is very high because of very high bottom cond
687	! flux... this is dangerous
688	dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
689	qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
690	! new lines
691	!-----------------
692	! Salinity update
693	!-----------------
694	! entrapment during bottom growth
695	dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + &
696	sm_i_b(ji)*ht_i_b(ji) ) / &
697	MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps ) &
698	- sm_i_b(ji)
699	ENDIF ! heat budget
700	END DO ! ji
701	ENDIF ! num_sal
702
703	! IF ( jiindex_1d .GT. 0 ) THEN
704	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
705	! WRITE(numout,*) ' s_i_new : ', s_i_new(jiindex_1d)
706	! ENDIF
707
708	!------------
709	! Basal melt
710	!------------
711	!++++++
712	meance_dh = 0.0
713	numce_dh = 0
714	innermelt(:) = 0
715	!++++++
716
717	DO ji = kideb, kiut
718	! heat convergence at the surface > 0
719	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN
720
721	s_i_new(ji) = s_i_b(ji,nlay_i)
722	zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
723
724	zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test
725	zdq_i(ji) = 0.0
726
727	dh_i_bott(ji) = 0.0
728	ENDIF
729	END DO
730
731	! IF ( jiindex_1d .GT. 0 ) THEN
732	! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
733	! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d)
734	! WRITE(numout,*) ' fbif : ', fbif_1d(jiindex_1d)
735	! WRITE(numout,*) ' qlbbq : ', qlbbq_1d(jiindex_1d)
736	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
737	! ENDIF
738
739	DO jk = nlay_i, 1, -1
740	DO ji = kideb, kiut
741	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN
742	ztmelts = - tmut * s_i_b(ji,jk) + rtt
743	! sometimes internal temperature gt melting point
744	! the whole layer is removed
745	! IF (ji.EQ.jiindex_1d) THEN
746	! WRITE(numout,*) ' jk : ', jk
747	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
748	! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)
749	! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)
750	! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)
751	! ENDIF
752	IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN
753	zdeltah(ji,jk) = - zh_i(ji)
754	! zqfont_bo(ji) = zqfont_bo(ji)
755	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk)
756	!++++++
757	innermelt(ji) = 1
758	!++++++
759	! IF (ji.EQ.jiindex_1d) THEN
760	! WRITE(numout,*) ' Inner melt: ', innermelt(ji)
761	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
762	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
763	! ENDIF
764	ELSE ! normal ablation
765	zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk)
766	zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
767	q_i_b(ji,jk)
768	zdeltah(ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) )
769	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk)
770	zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * &
771	q_i_b(ji,jk) / rdt_ice
772	! contribution to salt flux
773	zji = MOD( npb(ji) - 1, jpi ) + 1
774	zjj = ( npb(ji) - 1 ) / jpi + 1
775	zfsalt_melt(ji) = zfsalt_melt(ji) + &
776	! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &
777	( sss_io(zji,zjj) - sm_i_b(ji) ) * &
778	a_i_b(ji) * &
779	MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice
780	! IF (ji.EQ.jiindex_1d) THEN
781	! WRITE(numout,*) ' No inner melt '
782	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
783	! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)
784	! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)
785	! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)
786	! WRITE(numout,*)
787	! ENDIF
788	ENDIF
789	! IF ( ji .EQ. jiindex_1d ) THEN
790	! WRITE(numout,*) ' basal melt '
791	! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)
792	! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji)
793	! WRITE(numout,*) ' zfsalt_melt good units : ', zfsalt_melt(ji) / ( sss_io(zji,zjj) + epsi16 )
794	! ENDIF
795	!+++++
796	ENDIF
797	END DO ! ji
798	END DO ! jk
799
800	IF ( con_i ) THEN
801	DO ji = kideb, kiut
802	IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN
803	!-------------------
804	! Conservation test
805	!-------------------
806	IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
807	numce_dh = numce_dh + 1
808	meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
809	ENDIF
810	IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. &
811	! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN
812	WRITE(numout,*) ' ALERTE heat loss for basal melt '
813	WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
814	WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
815	WRITE(numout,*) ' zfbase : ', zfbase(ji)
816	WRITE(numout,*) ' zdq_i : ', zdq_i(ji)
817	WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)
818	WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
819	WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)
820	WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)
821	WRITE(numout,*) ' s_i_new : ', s_i_new(ji)
822	WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)
823	WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
824	WRITE(numout,*) ' innermelt : ', innermelt(ji)
825	ENDIF
826	ENDIF ! heat convergence at the surface
827	END DO ! ji
828
829	IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh
830	WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
831	WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
832	WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
833
834	ENDIF ! con_i
835
836	! IF ( jiindex_1d .GT. 0 ) THEN
837	! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
838	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)
839	! ENDIF
840	!
841	!------------------------------------------------------------------------------!
842	! 5) Pathological cases !
843	!------------------------------------------------------------------------------!
844	!
845	!------------------------------------------
846	! Excessive ablation in a 1-category model
847	!------------------------------------------
848
849	DO ji = kideb, kiut
850	! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
851	zdhbf = dh_i_bott(ji)
852	IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )
853	! excessive energy is sent to lateral ablation
854	fsup(ji) = rhoiclfus at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &
855	* ( zdhbf - dh_i_bott(ji) ) / rdt_ice
856
857	dh_i_bott(ji) = zdhbf
858	!since ice volume is only used for outputs, we keep it global for all categories
859	dvbbq_1d(ji) = a_i_b(ji)*dh_i_bott(ji)
860	!new ice thickness
861	zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
862
863	! diagnostic ( bottom ice growth )
864	zji = MOD( npb(ji) - 1, jpi ) + 1
865	zjj = ( npb(ji) - 1 ) / jpi + 1
866	diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) &
867	/ rdt_ice
868	diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) &
869	/ rdt_ice
870	diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) &
871	/ rdt_ice
872	END DO
873
874	! IF (jiindex_1d .GT. 0 ) THEN
875	! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
876	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
877	! ENDIF
878
879	!--------------------------
880	! If too much ice melts
881	!--------------------------
882	! then heat applied minus heat content at previous time step
883	! should equal heat remaining
884	!
885	DO ji = kideb, kiut
886	! Adapt the remaining energy if too much ice melts
887	!--------------------------------------------------
888	zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
889	! 0 if no more ice
890	zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0
891	! remaining heat
892	zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) ! &
893	! + zihgnew * MAX ( zfdt_init(ji) - zqt_i(ji) - zqt_s(ji), 0.0 )
894
895	! !++++++++++
896	! IF (ji.EQ.jiindex_1d) THEN
897	! WRITE(numout,*) ' zfdt_init : ', zfdt_init(ji) / rdt_ice
898	! WRITE(numout,*) ' zqt_i : ', zqt_i(ji) / rdt_ice
899	! WRITE(numout,*) ' 1. zqt_s : ', zqt_s(ji) / rdt_ice
900	! WRITE(numout,*) ' zqt : ', ( zqt_i(ji) + zqt_s(ji) ) / rdt_ice
901	! WRITE(numout,*) ' zhgnew : ', zhgnew(ji)
902	! WRITE(numout,*) ' zihgnew : ', zihgnew
903	! WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
904	! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
905	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)
906	! WRITE(numout,*) ' 1. zfdt_final: ', zfdt_final(ji) / rdt_ice
907	! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) / rdt_ice
908	! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice
909	! ENDIF
910	! !++++++++++
911
912	! If snow remains, energy is used to melt snow
913	zhni = ht_s_b(ji) ! snow depth at previous time step
914	zihg = MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow
915
916	! energy of melting of remaining snow
917	zqt_s(ji) = ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps )
918	zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) / &
919	MAX( zqt_s(ji) , zeps ) )
920	zhnfi = zhni + zdhnm
921	zfdt_final(ji) = MAX ( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
922	ht_s_b(ji) = MAX( zzero , zhnfi )
923	zqt_s(ji) = zqt_s(ji) * ht_s_b(ji)
924
925	! IF (ji.EQ.jiindex_1d) THEN
926	! WRITE(numout,*) ' 2. zfdt_final : ', zfdt_final(ji) / rdt_ice
927	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)
928	! WRITE(numout,*) ' zhgnew : ', zhgnew(ji)
929	! WRITE(numout,*) ' 2. zqt_s : ', zqt_s(ji) / rdt_ice
930	! WRITE(numout,*) ' zdhnm : ', zdhnm
931	! WRITE(numout,*)
932	! ENDIF
933
934	! Mass variations of ice and snow
935	!---------------------------------
936	rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * &
937	(zhgnew(ji)-ht_i_b(ji))*rhoic ! good
938
939	rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * &
940	(ht_s_b(ji)-zhni)*rhosn ! good too
941
942	! Remaining heat to the ocean
943	!---------------------------------
944	! check the switches here
945	! focea is in W.m-2 * dt
946
947	focea(ji) = - zfdt_final(ji) / rdt_ice
948
949	! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN
950	! WRITE(numout,*) ' focea : ', focea(ji)
951	! ENDIF
952
953	END DO
954
955	! IF (jiindex_1d .GT. 0 ) THEN
956	! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
957	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
958	! ENDIF
959
960	ftotal_fin (:) = zfdt_final(:) / rdt_ice
961	! IF (jiindex_1d .GT. 1 ) WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)
962
963	!---------------------------
964	! Salt flux and heat fluxes
965	!---------------------------
966	DO ji = kideb, kiut
967	zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
968
969	! Salt flux
970	zji = MOD( npb(ji) - 1, jpi ) + 1
971	zjj = ( npb(ji) - 1 ) / jpi + 1
972	IF ( num_sal .NE. 4 ) &
973	fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + &
974	(1.0 - zihgnew) * rdmicif_1d(ji) * &
975	( sss_io(zji,zjj) - sm_i_b(ji) ) / rdt_ice
976	! new lines
977	IF ( num_sal .EQ. 4 ) &
978	fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + &
979	(1.0 - zihgnew) * rdmicif_1d(ji) * &
980	( sss_io(zji,zjj) - bulk_sal ) / rdt_ice
981	! Heat flux
982	! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
983	! excessive total ablation energy (focea) sent to the ocean
984	qfvbq_1d(ji) = qfvbq_1d(ji) + &
985	fsup(ji) + ( 1.0 - zihgnew ) * &
986	focea(ji) * a_i_b(ji) * rdt_ice
987
988	! ! astarojna
989	zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
990	! equals 0 if ht_i = 0, 1 if ht_i gt 0
991	! fstbif_1d est cumulatif merde
992	fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji)
993	! here there is a bug
994	qldif_1d(ji) = qldif_1d(ji) &
995	+ fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) &
996	* rdt_ice &
997	+ ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
998	END DO ! ji
999
1000	!-------------------------------------------
1001	! Correct temperature, energy and thickness
1002	!-------------------------------------------
1003	DO ji = kideb, kiut
1004	zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )
1005	t_su_b(ji) = zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
1006	END DO ! ji
1007
1008	DO jk = 1, nlay_i
1009	DO ji = kideb, kiut
1010	zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )
1011	t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
1012	q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk)
1013	! IF (ji.eq.jiindex_1d) THEN
1014	! WRITE(numout,*) ' jk : ', jk
1015	! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)
1016	! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)
1017	! WRITE(numout,*) ' zihgnew : ', zihgnew, ' zhgnew : ', &
1018	! zhgnew(ji)
1019	! ENDIF
1020	END DO
1021	END DO ! ji
1022
1023	! IF (jiindex_1d.GT.0) THEN
1024	! WRITE(numout,*) ' --- Vertical profile : '
1025	! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)
1026	! ENDIF
1027
1028	DO ji = kideb, kiut
1029	ht_i_b(ji) = zhgnew(ji)
1030	END DO ! ji
1031
1032	! IF (jiindex_1d .GT. 0 ) THEN
1033	! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
1034	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
1035	! ENDIF
1036	!
1037	!------------------------------------------------------------------------------\|
1038	! 6) Snow-Ice formation \|
1039	!------------------------------------------------------------------------------\|
1040	!
1041	! When snow load excesses Archimede's limit, snow-ice interface goes down
1042	! under sea-level, flooding of seawater transforms snow into ice
1043	! dh_snowice is positive for the ice
1044	DO ji = kideb, kiut
1045
1046	dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) &
1047	* ht_i_b(ji))/(rhosn+rau0-rhoic))
1048	zhgnew(ji) = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji))
1049	zhnnew = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji))
1050
1051	! Changes in ice volume and ice mass.
1052	dvnbq_1d(ji) = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji))
1053	dmgwi_1d(ji) = dmgwi_1d(ji) + a_i_b(ji) &
1054	(ht_s_b(ji)-zhnnew)rhosn
1055
1056	#if defined key_lim_fdd
1057	!(presently Activated)
1058	rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) &
1059	* ( zhgnew(ji) - ht_i_b(ji) )*rhoic
1060	rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) &
1061	* ( zhnnew - ht_s_b(ji) )*rhosn
1062
1063	! Equivalent salt flux (1) Snow-ice formation component
1064	! -----------------------------------------------------
1065	zji = MOD( npb(ji) - 1, jpi ) + 1
1066	zjj = ( npb(ji) - 1 ) / jpi + 1
1067	! zsm_snowice = MIN ( s_i_max, ( rhoic - rhosn ) / rhoic * &
1068	! sss_io(zji,zjj) )
1069
1070	zsm_snowice = ( rhoic - rhosn ) / rhoic * &
1071	sss_io(zji,zjj)
1072
1073	IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)
1074
1075	IF ( num_sal .NE. 4 ) &
1076	fseqv_1d(ji) = fseqv_1d(ji) + &
1077	( sss_io(zji,zjj) - zsm_snowice ) * &
1078	a_i_b(ji) * &
1079	( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
1080	! new lines
1081	IF ( num_sal .EQ. 4 ) &
1082	fseqv_1d(ji) = fseqv_1d(ji) + &
1083	( sss_io(zji,zjj) - bulk_sal ) * &
1084	a_i_b(ji) * &
1085	( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
1086
1087	! entrapment during snow ice formation
1088	i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0d-6 ) )
1089	isnowic = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
1090	i_ice_switch
1091	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
1092	dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) &
1093	+ sm_i_b(ji) * ht_i_b(ji) &
1094	/ MAX( ht_i_b(ji) + dh_snowice(ji), zeps) &
1095	- sm_i_b(ji) ) * isnowic
1096
1097	#else
1098	rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) &
1099	* ( zhgnew(ji) - ht_i_b(ji) ) * rhoic &
1100	+ ( zhnnew - ht_s_b(ji) ) * rhosn )
1101	#endif
1102	! Actualize new snow and ice thickness.
1103	ht_s_b(ji) = zhnnew
1104	ht_i_b(ji) = zhgnew(ji)
1105
1106	! Total ablation ! new lines added to debug
1107	IF( ht_i_b(ji).LE.0.0 ) a_i_b(ji) = 0.0
1108
1109	! diagnostic ( snow ice growth )
1110	zji = MOD( npb(ji) - 1, jpi ) + 1
1111	zjj = ( npb(ji) - 1 ) / jpi + 1
1112	diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
1113	rdt_ice
1114
1115	END DO !ji
1116
1117	! IF (jiindex_1d .GT. 1 ) THEN
1118	! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)
1119	! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)
1120	! ENDIF
1121
1122	END SUBROUTINE lim_thd_dh
1123	#else
1124	!!======================================================================
1125	!! * MODULE limthd_dh *
1126	!! no sea ice model
1127	!!======================================================================
1128	CONTAINS
1129	SUBROUTINE lim_thd_dh ! Empty routine
1130	END SUBROUTINE lim_thd_dh
1131	#endif
1132	END MODULE limthd_dh

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

Download in other formats: