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 branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

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

Download in other formats: