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

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

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

Download in other formats: