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

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

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

Download in other formats: