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

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago
merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113
File size: 46.4 KB

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

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

Download in other formats: