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

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

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

Download in other formats: