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

Context Navigation

source: branches/CMIP5_IPSL/NEMO/LIM_SRC_3/limupdate.F90 @ 8506

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

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

Download in other formats: