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.

limupdate.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limupdate.F90 @ 913

Last change on this file since 913 was 895, checked in by rblod, 16 years ago
Fix small bugs in sea-ice (previous commit) and old flyspray: 189, 185, 161
File size: 45.8 KB

Line
1	MODULE limupdate
2	!!======================================================================
3	!! * MODULE limupdate *
4	!! Update of sea-ice global variables
5	!! at the end of the time step
6	!!
7	!! Ice speed from ice dynamics
8	!! Ice thickness, Snow thickness, Temperatures, Lead fraction
9	!! from advection and ice thermodynamics
10	!!======================================================================
11	#if defined key_lim3
12	!!----------------------------------------------------------------------
13	!! 'key_lim3' LIM3 sea-ice model
14	!!----------------------------------------------------------------------
15	!! lim_update : computes update of sea-ice global variables
16	!! from trend terms
17	!!----------------------------------------------------------------------
18	!! * Modules used
19	USE limistate
20	USE limrhg ! ice rheology
21	USE lbclnk
22
23	USE dom_oce
24	USE oce ! dynamics and tracers variables
25	USE in_out_manager
26	USE ice_oce ! ice variables
27	USE sbc_oce ! Surface boundary condition: ocean fields
28	USE sbc_ice ! Surface boundary condition: ice fields
29	USE dom_ice
30	USE daymod
31	USE phycst ! Define parameters for the routines
32	USE ice
33	USE iceini
34	USE lbclnk
35	USE limdyn
36	USE limtrp
37	USE limthd
38	USE limsbc
39	USE limdia
40	USE limwri
41	USE limrst
42	USE thd_ice ! LIM thermodynamic sea-ice variables
43	USE par_ice
44	USE limitd_th
45	USE limvar
46	USE prtctl ! Print control
47
48
49	IMPLICIT NONE
50	PRIVATE
51
52	!! * Accessibility
53	PUBLIC lim_update ! routine called by ice_step
54
55	!! * Substitutions
56	# include "vectopt_loop_substitute.h90"
57
58	!! * Module variables
59
60	!!----------------------------------------------------------------------
61	!! LIM 2.1 , UCL-ASTR-LODYC-IPSL (2004)
62	!!----------------------------------------------------------------------
63
64	CONTAINS
65
66	SUBROUTINE lim_update
67	!!-------------------------------------------------------------------
68	!! * ROUTINE lim_update *
69	!!
70	!! ** Purpose : Computes update of sea-ice global variables at
71	!! the end of the time step.
72	!! Address pathological cases
73	!! This place is very important
74	!!
75	!! ** Method : Mathematical
76	!!
77	!! ** Action : -
78	!!
79	!! History : This routine was new for LIM 3.0
80	!! 3.0 ! 04-06 (M. Vancoppenolle) Tendencies
81	!!---------------------------------------------------------------------
82	!! * Local variables
83	INTEGER :: &
84	ji, jj, & ! geographical indices
85	jk, jl, jm ! layer, category and type indices
86	INTEGER :: &
87	jbnd1, jbnd2
88	INTEGER :: &
89	i_ice_switch
90
91	REAL(wp) :: & ! constant values
92	epsi06 = 1.e-06 , &
93	epsi03 = 1.e-03 , &
94	epsi16 = 1.e-16 , &
95	epsi20 = 1.e-20 , &
96	epsi04 = 1.e-04 , &
97	epsi10 = 1.e-10 , &
98	rzero = 0.e0 , &
99	rone = 1.e0 , &
100	zhimax ! maximum thickness tolerated for advection of
101	! in an ice-free cell
102	REAL(wp) :: & ! dummy switches and arguments
103	zindb, zindsn, zindic, zacrith, &
104	zrtt, zindg, zh, zdvres, zviold, &
105	zbigvalue, zvsold, z_da_ex, zamax, &
106	z_prescr_hi, zat_i_old, &
107	ztmelts, ze_s
108
109	REAL(wp), DIMENSION(jpl) :: z_da_i, z_dv_i
110
111	LOGICAL, DIMENSION(jpi,jpj,jpl) :: &
112	internal_melt
113
114	INTEGER :: &
115	ind_im, layer ! indices for internal melt
116	REAL(wp), DIMENSION(jkmax) :: &
117	zthick0, zqm0 ! thickness of the layers and heat contents for
118	! internal melt
119	REAL(wp) :: &
120	zweight, zesum
121
122
123	!!-------------------------------------------------------------------
124
125	WRITE(numout,*) ' lim_update '
126	WRITE(numout,*) ' ~~~~~~~~~~ '
127
128	!+++++ [
129	WRITE(numout,*) ' O) Initial values '
130	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
131	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
132	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
133	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
134	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
135	DO jk = 1, nlay_i
136	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
137	END DO
138	!+++++ ]
139
140	!------------------------------------------------------------------------------
141	! 1. Update of Global variables \|
142	!------------------------------------------------------------------------------
143
144	!---------------------
145	! Ice dynamics
146	!---------------------
147
148	u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
149	v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
150
151	!-----------------------------
152	! Update ice and snow volumes
153	!-----------------------------
154
155	DO jl = 1, jpl
156	DO jj = 1, jpj
157	DO ji = 1, jpi
158
159	v_i(ji,jj,jl) = v_i(ji,jj,jl) + d_v_i_trp(ji,jj,jl) &
160	+ d_v_i_thd(ji,jj,jl)
161	v_s(ji,jj,jl) = v_s(ji,jj,jl) + d_v_s_trp(ji,jj,jl) &
162	+ d_v_s_thd(ji,jj,jl)
163	END DO
164	END DO
165	END DO
166
167	!---------------------------------
168	! Classify the pathological cases
169	!---------------------------------
170	! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
171	! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
172	! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
173	! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation
174	! with negative advection, very pathological )
175	! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
176
177	DO jl = 1, jpl
178	DO jj = 1, jpj
179	DO ji = 1, jpi
180	patho_case(ji,jj,jl) = 1
181	IF ( v_i(ji,jj,jl) .GE. 0.0 ) THEN
182	IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
183	patho_case(ji,jj,jl) = 2
184	ENDIF
185	ELSE
186	patho_case(ji,jj,jl) = 3
187	IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN
188	patho_case(ji,jj,jl) = 4
189	ENDIF
190	ENDIF
191	IF ( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. &
192	( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
193	patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free
194	! cell
195	WRITE(numout,*) ' ALERTE patho_case still equal to 5 '
196	WRITE(numout,*) ' ji , jj : ', ji, jj
197	WRITE(numout,*) ' old_v_i : ', old_v_i(ji,jj,jl)
198	WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl)
199
200	ENDIF
201	END DO
202	END DO
203	END DO
204
205	!--------------------
206	! Excessive ablation
207	!--------------------
208
209	DO jl = 1, jpl
210	DO jj = 1, jpj
211	DO ji = 1, jpi
212	IF ( ( patho_case(ji,jj,jl) .EQ. 3 ) &
213	.OR. ( patho_case(ji,jj,jl) .EQ. 4 ) ) THEN
214	zviold = old_v_i(ji,jj,jl)
215	zvsold = old_v_s(ji,jj,jl)
216	! in cases 3 ( combined total ablation )
217	! and 4 ( total ablation with negative advection )
218	! there is excessive total ablation
219	! advection is chosen to be prioritary in order to conserve mass.
220	! dv_i_thd is computed as a residual
221	! negative energy has to be kept in memory and to be given to the ocean
222	! equivalent salt flux is given to the ocean
223	!
224	! This was the best solution found. Otherwise, mass conservation in advection
225	! scheme should have been revised, which could have been a big problem
226	! Martin Vancoppenolle (2006, updated 2007)
227
228	! is there any ice left ?
229	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
230	!=1 if hi > 1e-3 and 0 if not
231	zdvres = MAX(0.0,-v_i(ji,jj,jl)) !residual volume if too much ice was molten
232	!this quantity is positive
233	v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl) !ice volume cannot be negative
234	!correct thermodynamic ablation
235	d_v_i_thd(ji,jj,jl) = zindic * d_v_i_thd(ji,jj,jl) + &
236	(1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl))
237	! THIS IS NEW
238	d_a_i_thd(ji,jj,jl) = zindic * d_a_i_thd(ji,jj,jl) + &
239	(1.0-zindic) * (-old_a_i(ji,jj,jl))
240
241	!residual salt flux if ice is over-molten
242	fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &
243	( rhoic * zdvres / rdt_ice )
244	! fheat_res(ji,jj) = fheat_res(ji,jj) + rhoic * lfus * zdvres / rdt_ice
245
246	! is there any snow left ?
247	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
248	zvsold = v_s(ji,jj,jl)
249	zdvres = MAX(0.0,-v_s(ji,jj,jl)) !residual volume if too much ice was molten
250	!this quantity is positive
251	v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl) !snow volume cannot be negative
252	!correct thermodynamic ablation
253	d_v_s_thd(ji,jj,jl) = zindsn * d_v_s_thd(ji,jj,jl) + &
254	(1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl))
255	!unsure correction on salt flux.... maybe future will tell it was not that right
256
257	!residual salt flux if snow is over-molten
258	fsalt_res(ji,jj) = fsalt_res(ji,jj) + sss_m(ji,jj) * &
259	( rhosn * zdvres / rdt_ice )
260	!this flux will be positive if snow was over-molten
261	! fheat_res(ji,jj) = fheat_res(ji,jj) + rhosn * lfus * zdvres / rdt_ice
262	ENDIF
263	END DO !ji
264	END DO !jj
265	END DO !jl
266
267	!+++++ [
268	DO jj = 1, jpj
269	DO ji = 1, jpi
270	IF ( ABS(fsalt_res(ji,jj)) .GT. 1.0 ) THEN
271	WRITE(numout,*) ' ALERTE 1000 : residual salt flux of -> ', &
272	fsalt_res(ji,jj)
273	WRITE(numout,*) ' ji, jj : ', ji, jj, ' gphit, glamt : ', &
274	gphit(ji,jj), glamt(ji,jj)
275	ENDIF
276	END DO
277	END DO
278
279	WRITE(numout,*) ' 1. Before update of Global variables '
280	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
281	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
282	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
283	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
284	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
285	DO jk = 1, nlay_i
286	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
287	END DO
288	!+++++ ]
289
290	!---------------------------------------------
291	! Ice concentration and ice heat content
292	!---------------------------------------------
293
294	a_i (:,:,:) = a_i (:,:,:) + d_a_i_trp(:,:,:) &
295	+ d_a_i_thd(:,:,:)
296	CALL lim_var_glo2eqv ! useless, just for debug
297	DO jk = 1, nlay_i
298	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
299	END DO
300	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)
301	CALL lim_var_glo2eqv ! useless, just for debug
302	WRITE(numout,*) ' After transport update '
303	DO jk = 1, nlay_i
304	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
305	END DO
306	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)
307	CALL lim_var_glo2eqv ! useless, just for debug
308	WRITE(numout,*) ' After thermodyn update '
309	DO jk = 1, nlay_i
310	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
311	END DO
312
313	at_i(:,:) = 0.0
314	DO jl = 1, jpl
315	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
316	END DO
317
318	!+++++ [
319	WRITE(numout,*) ' 1. After update of Global variables (2) '
320	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
321	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
322	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
323	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
324	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
325	WRITE(numout,*) ' oa_i : ', oa_i(jiindx, jjindx, 1:jpl)
326	WRITE(numout,*) ' e_s : ', e_s(jiindx, jjindx, 1, 1:jpl)
327	DO jk = 1, nlay_i
328	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
329	END DO
330	!+++++ ]
331
332	!------------------------------
333	! Snow temperature and ice age
334	!------------------------------
335
336	e_s(:,:,:,:) = e_s(:,:,:,:) + &
337	d_e_s_trp(:,:,:,:) + &
338	d_e_s_thd(:,:,:,:)
339
340	oa_i(:,:,:) = oa_i(:,:,:) + &
341	d_oa_i_trp(:,:,:) + &
342	d_oa_i_thd(:,:,:)
343
344	!--------------
345	! Ice salinity
346	!--------------
347
348	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
349
350	!+++++
351	WRITE(numout,*) ' Before everything '
352	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
353	WRITE(numout,*) ' oa_i: ', oa_i(jiindx, jjindx, 1:jpl)
354	DO jk = 1, nlay_i
355	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
356	END DO
357	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
358	!+++++
359
360	smv_i(:,:,:) = smv_i(:,:,:) + &
361	d_smv_i_thd(:,:,:) + &
362	d_smv_i_trp(:,:,:)
363
364	!+++++
365	WRITE(numout,*) ' After advection '
366	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
367	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
368	!+++++
369
370	ENDIF ! num_sal .EQ. 2
371
372	CALL lim_var_glo2eqv
373
374	!--------------------------------------
375	! 2. Review of all pathological cases
376	!--------------------------------------
377
378	zrtt = 173.15 * rone
379	zacrith = 1.0e-6
380
381	!-------------------------------------------
382	! 2.1) Advection of ice in an ice-free cell
383	!-------------------------------------------
384	! should be removed since it is treated after dynamics now
385
386	zhimax = 5.0
387	! first category
388	DO jj = 1, jpj
389	DO ji = 1, jpi
390	!--- the thickness of such an ice is often out of bounds
391	!--- thus we recompute a new area while conserving ice volume
392	zat_i_old = SUM(old_a_i(ji,jj,:))
393	zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,1)) - epsi10 ) )
394	IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb.GT.zhimax) &
395	.AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb).GT.zhimax ) &
396	.AND.( zat_i_old.LT.zacrith ) ) THEN ! new line
397	z_prescr_hi = hi_max(1) / 2.0
398	a_i(ji,jj,1) = v_i(ji,jj,1) / z_prescr_hi
399	ENDIF
400	END DO
401	END DO
402
403	!+++++ [
404	WRITE(numout,*) ' 2.1 '
405	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
406	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
407	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
408	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
409	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
410	DO jk = 1, nlay_i
411	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
412	END DO
413	!+++++ ]
414
415	!change this 14h44
416	zhimax = 20.0 ! line added up
417	! change this also 17 aug
418	zhimax = 30.0 ! line added up
419
420	DO jl = 2, jpl
421	jm = ice_types(jl)
422	DO jj = 1, jpj
423	DO ji = 1, jpi
424	zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) )
425	! this correction is very tricky... sometimes, advection gets wrong i don't know why
426	! it makes problems when the advected volume and concentration do not seem to be
427	! related with each other
428	! the new thickness is sometimes very big!
429	! and sometimes d_a_i_trp and d_v_i_trp have different sign
430	! which of course is plausible
431	! but fuck! it fucks everything up :)
432	IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb.GT.zhimax) &
433	.AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb).GT.zhimax ) THEN
434	z_prescr_hi = ( hi_max_typ(jl-ice_cat_bounds(jm,1) ,jm) + &
435	hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / &
436	2.0
437	a_i(ji,jj,jl) = v_i(ji,jj,jl) / z_prescr_hi
438	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
439	ENDIF
440	zat_i_old = SUM(old_a_i(ji,jj,:))
441
442	END DO ! ji
443	END DO !jj
444	END DO !jl
445
446	!+++++ [
447	WRITE(numout,*) ' 2.1 initial '
448	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
449	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
450	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
451	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
452	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
453	DO jk = 1, nlay_i
454	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
455	END DO
456	!+++++ ]
457
458	at_i(:,:) = 0.0
459	DO jl = 1, jpl
460	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
461	END DO
462
463	!----------------------------------------------------
464	! 2.2) Rebin categories with thickness out of bounds
465	!----------------------------------------------------
466	!+++++ [
467	WRITE(numout,*) ' 2.1 before rebinning '
468	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
469	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
470	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
471	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
472	DO jk = 1, nlay_i
473	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
474	END DO
475	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
476	!+++++ ]
477
478	DO jm = 1, jpm
479	jbnd1 = ice_cat_bounds(jm,1)
480	jbnd2 = ice_cat_bounds(jm,2)
481	IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
482	END DO
483
484
485	!+++++ [
486	WRITE(numout,*) ' 2.1 after rebinning'
487	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
488	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
489	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
490	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
491	DO jk = 1, nlay_i
492	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
493	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
494	END DO
495	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
496	!+++++ ]
497
498	at_i(:,:) = 0.0
499	DO jl = 1, jpl
500	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
501	END DO
502
503	!---------------------------------
504	! 2.3) Melt of an internal layer
505	!---------------------------------
506	internal_melt(:,:,:) = .false.
507
508	DO jl = 1, jpl
509	DO jk = 1, nlay_i
510	DO jj = 1, jpj
511	DO ji = 1, jpi
512	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
513	IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
514	( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) .AND. &
515	( v_i(ji,jj,jl) .GT. 0.0 ) .AND. &
516	( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
517	! WRITE(numout,*) ' Internal layer melt : '
518	! WRITE(numout,*) ' ji, jj, jk, jl : ', ji,jj,jk,jl
519	! WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
520	! WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
521	internal_melt(ji,jj,jl) = .true.
522	ENDIF
523	END DO ! ji
524	END DO ! jj
525	END DO !jk
526	END DO !jl
527
528	DO jl = 1, jpl
529	DO jj = 1, jpj
530	DO ji = 1, jpi
531	IF ( internal_melt(ji,jj,jl) ) THEN
532	! initial ice thickness
533	!-----------------------
534	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
535	! WRITE(numout,*) ' ji,jj,jl : ', ji,jj,jl
536	! WRITE(numout,*) ' old ht_i: ', ht_i(ji,jj,jl)
537	! WRITE(numout,*) ' Enthalpy at the beg : ', e_i(ji,jj,1:nlay_i,jl)
538	! WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl)
539
540	! reduce ice thickness
541	!-----------------------
542	ind_im = 0
543	zesum = 0.0
544	DO jk = 1, nlay_i
545	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
546	IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
547	( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
548	ind_im = ind_im + 1
549	zesum = zesum + e_i(ji,jj,jk,jl)
550	END DO
551	IF (ind_im .LT.nlay_i ) smv_i(ji,jj,jl)= smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * &
552	( ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i )
553	ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i
554	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
555
556	! WRITE(numout,*) ' ind_im : ', ind_im
557	! WRITE(numout,*) ' new ht_i: ', ht_i(ji,jj,jl)
558	! WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl)
559	! WRITE(numout,*) ' zesum : ', zesum
560
561	! redistribute heat
562	!-----------------------
563	! old thicknesses and enthalpies
564	ind_im = 0
565	DO jk = 1, nlay_i
566	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
567	IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND. &
568	( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
569	ind_im = ind_im + 1
570	zthick0(ind_im) = ht_i(ji,jj,jl) * ind_im / nlay_i
571	zqm0 (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
572	ENDIF
573	END DO
574
575	! WRITE(numout,*) ' Old thickness, enthalpy '
576	! WRITE(numout,*) ' Number of layer : ind_im ', ind_im
577	! WRITE(numout,*) ' zthick0 : ', zthick0(1:ind_im)
578	! WRITE(numout,*) ' zqm0 : ', zqm0(1:ind_im)
579
580	! Redistributing energy on the new grid
581	IF ( ind_im .GT. 0 ) THEN
582
583	DO jk = 1, nlay_i
584	e_i(ji,jj,jk,jl) = 0.0
585	DO layer = 1, ind_im
586	zweight = MAX ( &
587	MIN( ht_i(ji,jj,jl) * layer / ind_im , ht_i(ji,jj,jl) * jk / nlay_i ) - &
588	MAX( ht_i(ji,jj,jl) * (layer-1) / ind_im , ht_i(ji,jj,jl) * (jk-1) / nlay_i ) , 0.0 ) &
589	/ ( ht_i(ji,jj,jl) / ind_im )
590
591	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
592	END DO !layer
593	END DO ! jk
594
595	zesum = 0.0
596	DO jk = 1, nlay_i
597	zesum = zesum + e_i(ji,jj,jk,jl)
598	END DO
599
600	! WRITE(numout,*) ' Enthalpy at the end : ', e_i(ji,jj,1:nlay_i,jl)
601	! WRITE(numout,*) ' Volume at the end : ', v_i(ji,jj,jl)
602	! WRITE(numout,*) ' zesum : ', zesum
603
604	ELSE ! ind_im .EQ. 0, total melt
605	e_i(ji,jj,jk,jl) = 0.0
606	ENDIF
607
608	ENDIF ! internal_melt
609
610	END DO ! ji
611	END DO !jj
612	END DO !jl
613	!+++++ [
614	WRITE(numout,*) ' 2.3 after melt of an internal ice layer '
615	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
616	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
617	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
618	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
619	DO jk = 1, nlay_i
620	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
621	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
622	END DO
623	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
624	!+++++ ]
625
626	internal_melt(:,:,:) = .false.
627
628	! Melt of snow
629	!--------------
630	DO jl = 1, jpl
631	DO jj = 1, jpj
632	DO ji = 1, jpi
633	! snow energy of melting
634	ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / &
635	MAX( v_s(ji,jj,jl), 1.0e-6 ) ! snow energy of melting
636
637	! If snow energy of melting smaller then Lf
638	! Then all snow melts and meltwater, heat go to the ocean
639	IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = .true.
640
641	!++++++
642	IF ( (ji.eq.jiindx) .AND. (jj.eq.jjindx) ) THEN
643	WRITE(numout,*) ' jl : ', jl
644	WRITE(numout,*) ' ze_s : ', ze_s
645	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl)
646	WRITE(numout,*) ' rhosn : ', rhosn
647	WRITE(numout,*) ' rhosn : ', lfus
648	WRITE(numout,*) ' area : ', area(ji,jj)
649	WRITE(numout,) ' rhosn lfus : ', rhosn * lfus
650	WRITE(numout,*) ' internal_melt : ', internal_melt(ji,jj,jl)
651	ENDIF
652	!++++++
653
654	END DO
655	END DO
656	END DO
657
658	DO jl = 1, jpl
659	DO jj = 1, jpj
660	DO ji = 1, jpi
661	IF ( internal_melt(ji,jj,jl) ) THEN
662	v_s(ji,jj,jl) = 0.0
663	e_s(ji,jj,1,jl) = 0.0
664	! ! release heat
665	fheat_res(ji,jj) = fheat_res(ji,jj) &
666	+ ze_s * v_s(ji,jj,jl) / rdt_ice
667	! release mass
668	rdmsnif(ji,jj) = rdmsnif(ji,jj) + rhosn * v_s(ji,jj,jl)
669	ENDIF
670	END DO
671	END DO
672	END DO
673
674	zbigvalue = 1.0d+20
675
676	DO jl = 1, jpl
677	DO jj = 1, jpj
678	DO ji = 1, jpi
679
680	!switches
681	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
682	!switch = 1 if a_i > 1e-06 and 0 if not
683	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
684	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
685	! bug fix 25 avril 2007
686	zindb = zindb*zindic
687
688	!--- 2.3 Correction to ice age
689	!------------------------------
690	! IF ((o_i(ji,jj,jl)-1.0)86400.0.gt.(rdt_icefloat(numit))) THEN
691	! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
692	! ENDIF
693	IF ((oa_i(ji,jj,jl)-1.0)86400.0.gt.(rdt_icenumit*a_i(ji,jj,jl))) THEN
694	oa_i(ji,jj,jl) = rdt_icenumit/86400.0a_i(ji,jj,jl)
695	ENDIF
696	oa_i(ji,jj,jl) = zindbzindicoa_i(ji,jj,jl)
697
698	!--- 2.4 Correction to snow thickness
699	!-------------------------------------
700	! ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
701	! v_s(ji,jj,jl) = MAX( zindb * v_s(ji,jj,jl), 0.0)
702	! snow thickness cannot be smaller than 1e-6
703	v_s(ji,jj,jl) = zindsnv_s(ji,jj,jl)zindb
704	v_s(ji,jj,jl) = v_s(ji,jj,jl) * MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) )
705
706	!--- 2.5 Correction to ice thickness
707	!-------------------------------------
708	! ice thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hi = 0
709	v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl), 0.0)
710	! ice thickness cannot be smaller than 1e-3
711	v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)
712
713	!--- 2.6 Snow is transformed into ice if the original ice cover disappears
714	!----------------------------------------------------------------------------
715	zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -v_i(ji,jj,jl) ) )
716	v_i(ji,jj,jl) = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0
717	v_s(ji,jj,jl) = ( rone - zindg ) * v_s(ji,jj,jl) + &
718	zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
719
720	!--- 2.7 Correction to ice concentrations
721	!--------------------------------------------
722	! if greater than 0, ice concentration cannot be smaller than 1e-10
723	a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
724	! then ice volume has to be corrected too...
725	! instead, zap small areas
726
727	!-------------------------
728	! 2.8) Snow heat content
729	!-------------------------
730
731	e_s(ji,jj,1,jl) = zindsn * &
732	( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) ) + &
733	( 1.0 - zindsn ) * 0.0
734
735	END DO ! ji
736	END DO ! jj
737	END DO ! jl
738
739	!+++++ [
740	WRITE(numout,*) ' 2.8 '
741	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
742	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
743	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
744	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
745	DO jk = 1, nlay_i
746	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
747	END DO
748	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
749	!+++++ ]
750
751	!------------------------
752	! 2.9) Ice heat content
753	!------------------------
754
755	DO jl = 1, jpl
756	DO jk = 1, nlay_i
757	DO jj = 1, jpj
758	DO ji = 1, jpi
759	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi06 ) )
760	! =1 if v_i > 1e-6 and 0 if not
761	e_i(ji,jj,jk,jl)= zindic * &
762	( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) + &
763	( 1.0 - zindic ) * 0.0
764	END DO ! ji
765	END DO ! jj
766	END DO !jk
767	END DO !jl
768
769	WRITE(numout,*) ' 2.9 '
770	DO jk = 1, nlay_i
771	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
772	END DO
773	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
774
775	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
776
777	!---------------------
778	! 2.11) Ice salinity
779	!---------------------
780
781	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
782
783	DO jl = 1, jpl
784	DO jk = 1, nlay_i
785	DO jj = 1, jpj
786	DO ji = 1, jpi
787	! salinity stays in bounds
788	smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), &
789	0.1 * v_i(ji,jj,jl) )
790	i_ice_switch = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
791	smv_i(ji,jj,jl) = i_ice_switch*smv_i(ji,jj,jl) + &
792	0.1(1.0-i_ice_switch)v_i(ji,jj,jl)
793	END DO ! ji
794	END DO ! jj
795	END DO !jk
796	END DO !jl
797
798	ENDIF
799
800	!+++++ [
801	WRITE(numout,*) ' 2.11 '
802	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
803	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
804	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
805	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
806	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
807	!+++++ ]
808
809	DO jm = 1, jpm
810	DO jj = 1, jpj
811	DO ji = 1, jpi
812	jl = ice_cat_bounds(jm,1)
813	!--- 2.12 Constrain the thickness of the smallest category above 5 cm
814	!----------------------------------------------------------------------
815	! the ice thickness of the smallest category should be higher than 5 cm
816	! we changed hiclim to 10
817	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
818	ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
819	zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi20 ) )
820	ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
821	! v_s(ji,jj,jl) = v_s(ji,jj,jl) * zh
822	ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
823	! v_i(ji,jj,jl) = v_i(ji,jj,jl) * zh
824	a_i (ji,jj,jl) = a_i(ji,jj,jl) /zh
825	END DO !ji
826	END DO !jj
827	END DO !jm
828	!+++++ [
829	WRITE(numout,*) ' 2.12 '
830	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
831	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
832	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
833	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
834	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
835	!+++++ ]
836
837	!--- 2.13 Total ice concentration should not exceed 1
838	!-----------------------------------------------------
839	zamax = amax
840	! 2.13.1) individual concentrations cannot exceed zamax
841	!------------------------------------------------------
842
843	at_i(:,:) = 0.0
844	DO jl = 1, jpl
845	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
846	END DO
847
848	! 2.13.2) Total ice concentration cannot exceed zamax
849	!----------------------------------------------------
850	at_i(:,:) = 0.0
851	DO jl = 1, jpl
852	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
853	END DO
854
855	DO jj = 1, jpj
856	DO ji = 1, jpi
857
858	! 0) Excessive area ?
859	z_da_ex = MAX( at_i(ji,jj) - zamax , 0.0 )
860
861	! 1) Count the number of existing categories
862	DO jl = 1, jpl
863	zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) )
864	zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) )
865	z_da_i(jl) = a_i(ji,jj,jl)zindbz_da_ex/MAX(at_i(ji,jj),epsi06)
866	z_dv_i(jl) = v_i(ji,jj,jl)*z_da_i(jl)/MAX(at_i(ji,jj),epsi06)
867	a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i(jl)
868	v_i(ji,jj,jl) = v_i(ji,jj,jl) + z_dv_i(jl)
869
870	END DO
871
872	END DO !ji
873	END DO !jj
874
875	!+++++ [
876	WRITE(numout,*) ' 2.13 '
877	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
878	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
879	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
880	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
881	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
882	!+++++ ]
883
884	at_i(:,:) = 0.0
885	DO jl = 1, jpl
886	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
887	END DO
888
889	DO jj = 1, jpj
890	DO ji = 1, jpi
891	IF (at_i(ji,jj).GT.1.0) THEN
892	WRITE(numout,*) ' lim_update ! : at_i > 1 -> PAS BIEN -> ALERTE '
893	WRITE(numout,*) ' ~~~~~~~~~~ at_i ', at_i(ji,jj)
894	WRITE(numout,*) ' Point ', ji, jj
895	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
896	DO jl = 1, jpl
897	WRITE(numout,) ' a_i ** ', a_i(ji,jj,jl), ' CAT no ', jl
898	WRITE(numout,) ' a_i_old ** ', old_a_i(ji,jj,jl), ' CAT no ', jl
899	WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
900	END DO
901	! WRITE(numout,*) ' CORRECTION BARBARE '
902	! z_da_ex = MAX( at_i(ji,jj) - zamax , 0.0 )
903	ENDIF
904	END DO
905	END DO
906
907	! Final thickness distribution rebinning
908	! --------------------------------------
909	!+++++ [
910	WRITE(numout,*) ' rebinning before'
911	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
912	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
913	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
914	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
915	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
916	!+++++ ]
917	!old version
918	! CALL lim_itd_th_reb(1,jpl)
919
920	DO jm = 1, jpm
921	jbnd1 = ice_cat_bounds(jm,1)
922	jbnd2 = ice_cat_bounds(jm,2)
923	IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
924	IF (ice_ncat_types(jm) .EQ. 1 ) THEN
925	ENDIF
926	END DO
927	!+++++ [
928	WRITE(numout,*) ' rebinning final'
929	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
930	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
931	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
932	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
933	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
934	!+++++ ]
935
936	at_i(:,:) = 0.0
937	DO jl = 1, jpl
938	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
939	END DO
940
941	!------------------------------------------------------------------------------
942	! 2) Corrections to avoid wrong values \|
943	!------------------------------------------------------------------------------
944	! Ice drift
945	!------------
946
947	DO jj = 2, jpjm1
948	DO ji = fs_2, fs_jpim1
949	IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
950	IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj) = 0.0 ! right side
951	IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
952	IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj) = 0.0 ! upper side
953	IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji-1,jj) = 0.0 ! bottom side
954	ENDIF
955	END DO
956	END DO
957	!mask velocities
958	u_ice(:,:) = u_ice(:,:) * tmu(:,:)
959	v_ice(:,:) = v_ice(:,:) * tmv(:,:)
960	!lateral boundary conditions
961	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
962	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
963
964	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
965	! ALERTES
966	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
967
968	DO jj = 1, jpj
969	DO ji = 1, jpi
970	DO jl = 1, jpl
971	! IF ((v_i(ji,jj,jl).NE.0.0).AND.(a_i(ji,jj,jl).EQ.0.0)) THEN
972	! WRITE(numout,*) ' lim_update : incompatible volume and concentration '
973	END DO ! jl
974
975	DO jl = 1, jpl
976	IF ( (a_i(ji,jj,jl).GT.1.0).OR.(at_i(ji,jj).GT.1.0) ) THEN
977	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
978	WRITE(numout,*) ' lim_update : a_i > 1 '
979	WRITE(numout,*) ' PAS BIEN ----> ALERTE !!! '
980	WRITE(numout,*) ' ~~~~~~~~~~ at_i ', at_i(ji,jj)
981	WRITE(numout,*) ' Point - category', ji, jj, jl
982	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
983	WRITE(numout,) ' a_i ** a_i_old ', a_i(ji,jj,jl), old_a_i(ji,jj,jl)
984	WRITE(numout,) ' v_i ** v_i_old ', v_i(ji,jj,jl), old_v_i(ji,jj,jl)
985	WRITE(numout,) ' ht_i ** ', v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi06)*zindb
986	WRITE(numout,*) ' hi_max(jl), hi_max(jl-1) ', hi_max(jl), hi_max(jl-1)
987	WRITE(numout,*) ' d_v_i_thd / trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
988	WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
989	ENDIF
990	END DO
991
992	END DO !jj
993	END DO !ji
994
995	WRITE(numout,*) ' TESTOSC1 ', tio_u(jiindx,jjindx), tio_v(jiindx,jjindx)
996	WRITE(numout,*) ' TESTOSC2 ', u_ice(jiindx,jjindx), v_ice(jiindx,jjindx)
997	WRITE(numout,*) ' TESTOSC3 ', u_oce(jiindx,jjindx), v_oce(jiindx,jjindx)
998	WRITE(numout,*) ' TESTOSC4 ', utau (jiindx,jjindx), vtau (jiindx,jjindx)
999
1000
1001	IF(ln_ctl) THEN ! Control print
1002	CALL prt_ctl_info(' ')
1003	CALL prt_ctl_info(' - Cell values : ')
1004	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
1005	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_update : cell area :')
1006	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update : at_i :')
1007	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update : vt_i :')
1008	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_update : vt_s :')
1009	CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update : strength :')
1010	CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
1011	CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
1012	CALL prt_ctl(tab2d_1=old_u_ice , clinfo1=' lim_update : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice :')
1013
1014	DO jl = 1, jpl
1015	CALL prt_ctl_info(' ')
1016	CALL prt_ctl_info(' - Category : ', ivar1=jl)
1017	CALL prt_ctl_info(' ~~~~~~~~~~')
1018	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_update : ht_i : ')
1019	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_update : ht_s : ')
1020	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_update : t_su : ')
1021	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_update : t_snow : ')
1022	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_update : sm_i : ')
1023	CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update : o_i : ')
1024	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update : a_i : ')
1025	CALL prt_ctl(tab2d_1=old_a_i (:,:,jl) , clinfo1= ' lim_update : old_a_i : ')
1026	CALL prt_ctl(tab2d_1=d_a_i_trp (:,:,jl) , clinfo1= ' lim_update : d_a_i_trp : ')
1027	CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update : d_a_i_thd : ')
1028	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update : v_i : ')
1029	CALL prt_ctl(tab2d_1=old_v_i (:,:,jl) , clinfo1= ' lim_update : old_v_i : ')
1030	CALL prt_ctl(tab2d_1=d_v_i_trp (:,:,jl) , clinfo1= ' lim_update : d_v_i_trp : ')
1031	CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update : d_v_i_thd : ')
1032	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update : v_s : ')
1033	CALL prt_ctl(tab2d_1=old_v_s (:,:,jl) , clinfo1= ' lim_update : old_v_s : ')
1034	CALL prt_ctl(tab2d_1=d_v_s_trp (:,:,jl) , clinfo1= ' lim_update : d_v_s_trp : ')
1035	CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update : d_v_s_thd : ')
1036	CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : e_i1 : ')
1037	CALL prt_ctl(tab2d_1=old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : old_e_i1 : ')
1038	CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : de_i1_trp : ')
1039	CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : de_i1_thd : ')
1040	CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : e_i2 : ')
1041	CALL prt_ctl(tab2d_1=old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : old_e_i2 : ')
1042	CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : de_i2_trp : ')
1043	CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : de_i2_thd : ')
1044	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update : e_snow : ')
1045	CALL prt_ctl(tab2d_1=old_e_s (:,:,1,jl) , clinfo1= ' lim_update : old_e_snow : ')
1046	CALL prt_ctl(tab2d_1=d_e_s_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : d_e_s_trp : ')
1047	CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : d_e_s_thd : ')
1048	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update : smv_i : ')
1049	CALL prt_ctl(tab2d_1=old_smv_i (:,:,jl) , clinfo1= ' lim_update : old_smv_i : ')
1050	CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl) , clinfo1= ' lim_update : d_smv_i_trp : ')
1051	CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update : d_smv_i_thd : ')
1052	CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update : oa_i : ')
1053	CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update : old_oa_i : ')
1054	CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update : d_oa_i_trp : ')
1055	CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update : d_oa_i_thd : ')
1056	CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update : Path. case : ')
1057
1058	DO jk = 1, nlay_i
1059	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
1060	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update : t_i : ')
1061	END DO
1062	END DO
1063
1064	CALL prt_ctl_info(' ')
1065	CALL prt_ctl_info(' - Heat / FW fluxes : ')
1066	CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ')
1067	CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')
1068	CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ')
1069	CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
1070
1071	CALL prt_ctl_info(' ')
1072	CALL prt_ctl_info(' - Stresses : ')
1073	CALL prt_ctl_info(' ~~~~~~~~~~ ')
1074	CALL prt_ctl(tab2d_1=utau , clinfo1= ' lim_update : utau : ', tab2d_2=vtau , clinfo2= ' vtau : ')
1075	CALL prt_ctl(tab2d_1=utaui_ice , clinfo1= ' lim_update : utaui_ice : ', tab2d_2=vtaui_ice , clinfo2= ' vtaui_ice : ')
1076	CALL prt_ctl(tab2d_1=u_oce , clinfo1= ' lim_update : u_oce : ', tab2d_2=v_oce , clinfo2= ' v_oce : ')
1077	ENDIF
1078
1079	!---------------------
1080
1081	END SUBROUTINE lim_update
1082	#else
1083	!!----------------------------------------------------------------------
1084	!! Default option Empty Module No sea-ice model
1085	!!----------------------------------------------------------------------
1086	CONTAINS
1087	SUBROUTINE lim_update ! Empty routine
1088	END SUBROUTINE lim_update
1089
1090	#endif
1091
1092	END MODULE limupdate

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

Download in other formats: