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/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 3817

Last change on this file since 3817 was 3817, checked in by clevy, 11 years ago
Bugfix on articifial ice source, see ticket:#1054
Property svn:keywords set to `Id`
File size: 47.3 KB

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

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

Download in other formats: