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

Last change on this file since 941 was 921, checked in by rblod, 16 years ago
Correct indentation and print for debug in LIM3, see ticket #134, step I
File size: 47.7 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	DO jk = 1, nlay_i
300	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
301	END DO
302	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)
303	CALL lim_var_glo2eqv ! useless, just for debug
304	WRITE(numout,*) ' After transport update '
305	DO jk = 1, nlay_i
306	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
307	END DO
308	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)
309	CALL lim_var_glo2eqv ! useless, just for debug
310	WRITE(numout,*) ' After thermodyn update '
311	DO jk = 1, nlay_i
312	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
313	END DO
314
315	at_i(:,:) = 0.0
316	DO jl = 1, jpl
317	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
318	END DO
319
320	IF( ln_nicep ) THEN
321	WRITE(numout,*) ' 1. After update of Global variables (2) '
322	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
323	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
324	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
325	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
326	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
327	WRITE(numout,*) ' oa_i : ', oa_i(jiindx, jjindx, 1:jpl)
328	WRITE(numout,*) ' e_s : ', e_s(jiindx, jjindx, 1, 1:jpl)
329	DO jk = 1, nlay_i
330	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
331	END DO
332	ENDIF
333
334	!------------------------------
335	! Snow temperature and ice age
336	!------------------------------
337
338	e_s(:,:,:,:) = e_s(:,:,:,:) + &
339	d_e_s_trp(:,:,:,:) + &
340	d_e_s_thd(:,:,:,:)
341
342	oa_i(:,:,:) = oa_i(:,:,:) + &
343	d_oa_i_trp(:,:,:) + &
344	d_oa_i_thd(:,:,:)
345
346	!--------------
347	! Ice salinity
348	!--------------
349
350	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
351
352	IF( ln_nicep ) THEN
353	WRITE(numout,*) ' Before everything '
354	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
355	WRITE(numout,*) ' oa_i: ', oa_i(jiindx, jjindx, 1:jpl)
356	DO jk = 1, nlay_i
357	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
358	END DO
359	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
360	ENDIF
361
362	smv_i(:,:,:) = smv_i(:,:,:) + &
363	d_smv_i_thd(:,:,:) + &
364	d_smv_i_trp(:,:,:)
365
366	IF( ln_nicep ) THEN
367	WRITE(numout,*) ' After advection '
368	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
369	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
370	ENDIF
371
372	ENDIF ! num_sal .EQ. 2
373
374	CALL lim_var_glo2eqv
375
376	!--------------------------------------
377	! 2. Review of all pathological cases
378	!--------------------------------------
379
380	zrtt = 173.15 * rone
381	zacrith = 1.0e-6
382
383	!-------------------------------------------
384	! 2.1) Advection of ice in an ice-free cell
385	!-------------------------------------------
386	! should be removed since it is treated after dynamics now
387
388	zhimax = 5.0
389	! first category
390	DO jj = 1, jpj
391	DO ji = 1, jpi
392	!--- the thickness of such an ice is often out of bounds
393	!--- thus we recompute a new area while conserving ice volume
394	zat_i_old = SUM(old_a_i(ji,jj,:))
395	zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,1)) - epsi10 ) )
396	IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb.GT.zhimax) &
397	.AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb).GT.zhimax ) &
398	.AND.( zat_i_old.LT.zacrith ) ) THEN ! new line
399	z_prescr_hi = hi_max(1) / 2.0
400	a_i(ji,jj,1) = v_i(ji,jj,1) / z_prescr_hi
401	ENDIF
402	END DO
403	END DO
404
405	IF( ln_nicep ) THEN
406	WRITE(numout,*) ' 2.1 '
407	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
408	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
409	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
410	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
411	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
412	DO jk = 1, nlay_i
413	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
414	END DO
415	ENDIF
416
417	!change this 14h44
418	zhimax = 20.0 ! line added up
419	! change this also 17 aug
420	zhimax = 30.0 ! line added up
421
422	DO jl = 2, jpl
423	jm = ice_types(jl)
424	DO jj = 1, jpj
425	DO ji = 1, jpi
426	zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) )
427	! this correction is very tricky... sometimes, advection gets wrong i don't know why
428	! it makes problems when the advected volume and concentration do not seem to be
429	! related with each other
430	! the new thickness is sometimes very big!
431	! and sometimes d_a_i_trp and d_v_i_trp have different sign
432	! which of course is plausible
433	! but fuck! it fucks everything up :)
434	IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb.GT.zhimax) &
435	.AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb).GT.zhimax ) THEN
436	z_prescr_hi = ( hi_max_typ(jl-ice_cat_bounds(jm,1) ,jm) + &
437	hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / &
438	2.0
439	a_i(ji,jj,jl) = v_i(ji,jj,jl) / z_prescr_hi
440	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
441	ENDIF
442	zat_i_old = SUM(old_a_i(ji,jj,:))
443
444	END DO ! ji
445	END DO !jj
446	END DO !jl
447
448	IF( ln_nicep ) THEN
449	WRITE(numout,*) ' 2.1 initial '
450	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
451	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
452	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
453	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
454	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
455	DO jk = 1, nlay_i
456	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
457	END DO
458	ENDIF
459
460	at_i(:,:) = 0.0
461	DO jl = 1, jpl
462	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
463	END DO
464
465	!----------------------------------------------------
466	! 2.2) Rebin categories with thickness out of bounds
467	!----------------------------------------------------
468	IF( ln_nicep ) THEN
469	WRITE(numout,*) ' 2.1 before rebinning '
470	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
471	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
472	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
473	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
474	DO jk = 1, nlay_i
475	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
476	END DO
477	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
478	ENDIF
479
480	DO jm = 1, jpm
481	jbnd1 = ice_cat_bounds(jm,1)
482	jbnd2 = ice_cat_bounds(jm,2)
483	IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
484	END DO
485
486
487	IF( ln_nicep ) THEN
488	WRITE(numout,*) ' 2.1 after rebinning'
489	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
490	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
491	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
492	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
493	DO jk = 1, nlay_i
494	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
495	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
496	END DO
497	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
498	ENDIF
499
500	at_i(:,:) = 0.0
501	DO jl = 1, jpl
502	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
503	END DO
504
505	!---------------------------------
506	! 2.3) Melt of an internal layer
507	!---------------------------------
508	internal_melt(:,:,:) = .false.
509
510	DO jl = 1, jpl
511	DO jk = 1, nlay_i
512	DO jj = 1, jpj
513	DO ji = 1, jpi
514	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
515	IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
516	( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) .AND. &
517	( v_i(ji,jj,jl) .GT. 0.0 ) .AND. &
518	( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
519	! WRITE(numout,*) ' Internal layer melt : '
520	! WRITE(numout,*) ' ji, jj, jk, jl : ', ji,jj,jk,jl
521	! WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
522	! WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
523	internal_melt(ji,jj,jl) = .true.
524	ENDIF
525	END DO ! ji
526	END DO ! jj
527	END DO !jk
528	END DO !jl
529
530	DO jl = 1, jpl
531	DO jj = 1, jpj
532	DO ji = 1, jpi
533	IF ( internal_melt(ji,jj,jl) ) THEN
534	! initial ice thickness
535	!-----------------------
536	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
537	! WRITE(numout,*) ' ji,jj,jl : ', ji,jj,jl
538	! WRITE(numout,*) ' old ht_i: ', ht_i(ji,jj,jl)
539	! WRITE(numout,*) ' Enthalpy at the beg : ', e_i(ji,jj,1:nlay_i,jl)
540	! WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl)
541
542	! reduce ice thickness
543	!-----------------------
544	ind_im = 0
545	zesum = 0.0
546	DO jk = 1, nlay_i
547	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
548	IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
549	( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
550	ind_im = ind_im + 1
551	zesum = zesum + e_i(ji,jj,jk,jl)
552	END DO
553	IF (ind_im .LT.nlay_i ) smv_i(ji,jj,jl)= smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * &
554	( ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i )
555	ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i
556	v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
557
558	! WRITE(numout,*) ' ind_im : ', ind_im
559	! WRITE(numout,*) ' new ht_i: ', ht_i(ji,jj,jl)
560	! WRITE(numout,*) ' smv_i : ', smv_i(ji,jj,jl)
561	! WRITE(numout,*) ' zesum : ', zesum
562
563	! redistribute heat
564	!-----------------------
565	! old thicknesses and enthalpies
566	ind_im = 0
567	DO jk = 1, nlay_i
568	ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
569	IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND. &
570	( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
571	ind_im = ind_im + 1
572	zthick0(ind_im) = ht_i(ji,jj,jl) * ind_im / nlay_i
573	zqm0 (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
574	ENDIF
575	END DO
576
577	! WRITE(numout,*) ' Old thickness, enthalpy '
578	! WRITE(numout,*) ' Number of layer : ind_im ', ind_im
579	! WRITE(numout,*) ' zthick0 : ', zthick0(1:ind_im)
580	! WRITE(numout,*) ' zqm0 : ', zqm0(1:ind_im)
581
582	! Redistributing energy on the new grid
583	IF ( ind_im .GT. 0 ) THEN
584
585	DO jk = 1, nlay_i
586	e_i(ji,jj,jk,jl) = 0.0
587	DO layer = 1, ind_im
588	zweight = MAX ( &
589	MIN( ht_i(ji,jj,jl) * layer / ind_im , ht_i(ji,jj,jl) * jk / nlay_i ) - &
590	MAX( ht_i(ji,jj,jl) * (layer-1) / ind_im , ht_i(ji,jj,jl) * (jk-1) / nlay_i ) , 0.0 ) &
591	/ ( ht_i(ji,jj,jl) / ind_im )
592
593	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
594	END DO !layer
595	END DO ! jk
596
597	zesum = 0.0
598	DO jk = 1, nlay_i
599	zesum = zesum + e_i(ji,jj,jk,jl)
600	END DO
601
602	! WRITE(numout,*) ' Enthalpy at the end : ', e_i(ji,jj,1:nlay_i,jl)
603	! WRITE(numout,*) ' Volume at the end : ', v_i(ji,jj,jl)
604	! WRITE(numout,*) ' zesum : ', zesum
605
606	ELSE ! ind_im .EQ. 0, total melt
607	e_i(ji,jj,jk,jl) = 0.0
608	ENDIF
609
610	ENDIF ! internal_melt
611
612	END DO ! ji
613	END DO !jj
614	END DO !jl
615	IF( ln_nicep ) THEN
616	WRITE(numout,*) ' 2.3 after melt of an internal ice layer '
617	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
618	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
619	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
620	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
621	DO jk = 1, nlay_i
622	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
623	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
624	END DO
625	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
626	ENDIF
627
628	internal_melt(:,:,:) = .false.
629
630	! Melt of snow
631	!--------------
632	DO jl = 1, jpl
633	DO jj = 1, jpj
634	DO ji = 1, jpi
635	! snow energy of melting
636	ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / &
637	MAX( v_s(ji,jj,jl), 1.0e-6 ) ! snow energy of melting
638
639	! If snow energy of melting smaller then Lf
640	! Then all snow melts and meltwater, heat go to the ocean
641	IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = .true.
642
643	IF( ln_nicep ) THEN
644	IF ( (ji.eq.jiindx) .AND. (jj.eq.jjindx) ) THEN
645	WRITE(numout,*) ' jl : ', jl
646	WRITE(numout,*) ' ze_s : ', ze_s
647	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl)
648	WRITE(numout,*) ' rhosn : ', rhosn
649	WRITE(numout,*) ' rhosn : ', lfus
650	WRITE(numout,*) ' area : ', area(ji,jj)
651	WRITE(numout,) ' rhosn lfus : ', rhosn * lfus
652	WRITE(numout,*) ' internal_melt : ', internal_melt(ji,jj,jl)
653	ENDIF
654	ENDIF
655
656	END DO
657	END DO
658	END DO
659
660	DO jl = 1, jpl
661	DO jj = 1, jpj
662	DO ji = 1, jpi
663	IF ( internal_melt(ji,jj,jl) ) THEN
664	v_s(ji,jj,jl) = 0.0
665	e_s(ji,jj,1,jl) = 0.0
666	! ! release heat
667	fheat_res(ji,jj) = fheat_res(ji,jj) &
668	+ ze_s * v_s(ji,jj,jl) / rdt_ice
669	! release mass
670	rdmsnif(ji,jj) = rdmsnif(ji,jj) + rhosn * v_s(ji,jj,jl)
671	ENDIF
672	END DO
673	END DO
674	END DO
675
676	zbigvalue = 1.0d+20
677
678	DO jl = 1, jpl
679	DO jj = 1, jpj
680	DO ji = 1, jpi
681
682	!switches
683	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
684	!switch = 1 if a_i > 1e-06 and 0 if not
685	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
686	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
687	! bug fix 25 avril 2007
688	zindb = zindb*zindic
689
690	!--- 2.3 Correction to ice age
691	!------------------------------
692	! IF ((o_i(ji,jj,jl)-1.0)86400.0.gt.(rdt_icefloat(numit))) THEN
693	! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
694	! ENDIF
695	IF ((oa_i(ji,jj,jl)-1.0)86400.0.gt.(rdt_icenumit*a_i(ji,jj,jl))) THEN
696	oa_i(ji,jj,jl) = rdt_icenumit/86400.0a_i(ji,jj,jl)
697	ENDIF
698	oa_i(ji,jj,jl) = zindbzindicoa_i(ji,jj,jl)
699
700	!--- 2.4 Correction to snow thickness
701	!-------------------------------------
702	! ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
703	! v_s(ji,jj,jl) = MAX( zindb * v_s(ji,jj,jl), 0.0)
704	! snow thickness cannot be smaller than 1e-6
705	v_s(ji,jj,jl) = zindsnv_s(ji,jj,jl)zindb
706	v_s(ji,jj,jl) = v_s(ji,jj,jl) * MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) )
707
708	!--- 2.5 Correction to ice thickness
709	!-------------------------------------
710	! ice thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hi = 0
711	v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl), 0.0)
712	! ice thickness cannot be smaller than 1e-3
713	v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)
714
715	!--- 2.6 Snow is transformed into ice if the original ice cover disappears
716	!----------------------------------------------------------------------------
717	zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -v_i(ji,jj,jl) ) )
718	v_i(ji,jj,jl) = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0
719	v_s(ji,jj,jl) = ( rone - zindg ) * v_s(ji,jj,jl) + &
720	zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
721
722	!--- 2.7 Correction to ice concentrations
723	!--------------------------------------------
724	! if greater than 0, ice concentration cannot be smaller than 1e-10
725	a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
726	! then ice volume has to be corrected too...
727	! instead, zap small areas
728
729	!-------------------------
730	! 2.8) Snow heat content
731	!-------------------------
732
733	e_s(ji,jj,1,jl) = zindsn * &
734	( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) ) + &
735	( 1.0 - zindsn ) * 0.0
736
737	END DO ! ji
738	END DO ! jj
739	END DO ! jl
740
741	IF( ln_nicep ) THEN
742	WRITE(numout,*) ' 2.8 '
743	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
744	WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
745	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
746	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
747	DO jk = 1, nlay_i
748	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
749	END DO
750	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
751	ENDIF
752
753	!------------------------
754	! 2.9) Ice heat content
755	!------------------------
756
757	DO jl = 1, jpl
758	DO jk = 1, nlay_i
759	DO jj = 1, jpj
760	DO ji = 1, jpi
761	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi06 ) )
762	! =1 if v_i > 1e-6 and 0 if not
763	e_i(ji,jj,jk,jl)= zindic * &
764	( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) + &
765	( 1.0 - zindic ) * 0.0
766	END DO ! ji
767	END DO ! jj
768	END DO !jk
769	END DO !jl
770
771	IF( ln_nicep ) THEN
772	WRITE(numout,*) ' 2.9 '
773	DO jk = 1, nlay_i
774	WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
775	END DO
776	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
777
778	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
779	ENDIF
780
781	!---------------------
782	! 2.11) Ice salinity
783	!---------------------
784
785	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
786
787	DO jl = 1, jpl
788	DO jk = 1, nlay_i
789	DO jj = 1, jpj
790	DO ji = 1, jpi
791	! salinity stays in bounds
792	smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), &
793	0.1 * v_i(ji,jj,jl) )
794	i_ice_switch = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
795	smv_i(ji,jj,jl) = i_ice_switch*smv_i(ji,jj,jl) + &
796	0.1(1.0-i_ice_switch)v_i(ji,jj,jl)
797	END DO ! ji
798	END DO ! jj
799	END DO !jk
800	END DO !jl
801
802	ENDIF
803
804	IF( ln_nicep ) THEN
805	WRITE(numout,*) ' 2.11 '
806	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
807	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
808	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
809	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
810	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
811	ENDIF
812
813	DO jm = 1, jpm
814	DO jj = 1, jpj
815	DO ji = 1, jpi
816	jl = ice_cat_bounds(jm,1)
817	!--- 2.12 Constrain the thickness of the smallest category above 5 cm
818	!----------------------------------------------------------------------
819	! the ice thickness of the smallest category should be higher than 5 cm
820	! we changed hiclim to 10
821	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
822	ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
823	zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi20 ) )
824	ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
825	! v_s(ji,jj,jl) = v_s(ji,jj,jl) * zh
826	ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
827	! v_i(ji,jj,jl) = v_i(ji,jj,jl) * zh
828	a_i (ji,jj,jl) = a_i(ji,jj,jl) /zh
829	END DO !ji
830	END DO !jj
831	END DO !jm
832	IF( ln_nicep ) THEN
833	WRITE(numout,*) ' 2.12 '
834	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
835	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
836	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
837	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
838	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
839	ENDIF
840
841	!--- 2.13 Total ice concentration should not exceed 1
842	!-----------------------------------------------------
843	zamax = amax
844	! 2.13.1) individual concentrations cannot exceed zamax
845	!------------------------------------------------------
846
847	at_i(:,:) = 0.0
848	DO jl = 1, jpl
849	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
850	END DO
851
852	! 2.13.2) Total ice concentration cannot exceed zamax
853	!----------------------------------------------------
854	at_i(:,:) = 0.0
855	DO jl = 1, jpl
856	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
857	END DO
858
859	DO jj = 1, jpj
860	DO ji = 1, jpi
861
862	! 0) Excessive area ?
863	z_da_ex = MAX( at_i(ji,jj) - zamax , 0.0 )
864
865	! 1) Count the number of existing categories
866	DO jl = 1, jpl
867	zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) )
868	zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) )
869	z_da_i(jl) = a_i(ji,jj,jl)zindbz_da_ex/MAX(at_i(ji,jj),epsi06)
870	z_dv_i(jl) = v_i(ji,jj,jl)*z_da_i(jl)/MAX(at_i(ji,jj),epsi06)
871	a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i(jl)
872	v_i(ji,jj,jl) = v_i(ji,jj,jl) + z_dv_i(jl)
873
874	END DO
875
876	END DO !ji
877	END DO !jj
878
879	IF( ln_nicep ) THEN
880	WRITE(numout,*) ' 2.13 '
881	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
882	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
883	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
884	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
885	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
886	ENDIF
887
888	at_i(:,:) = 0.0
889	DO jl = 1, jpl
890	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
891	END DO
892
893	IF( ln_nicep ) THEN
894	DO jj = 1, jpj
895	DO ji = 1, jpi
896	IF (at_i(ji,jj).GT.1.0) THEN
897	WRITE(numout,*) ' lim_update ! : at_i > 1 -> PAS BIEN -> ALERTE '
898	WRITE(numout,*) ' ~~~~~~~~~~ at_i ', at_i(ji,jj)
899	WRITE(numout,*) ' Point ', ji, jj
900	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
901	DO jl = 1, jpl
902	WRITE(numout,) ' a_i ** ', a_i(ji,jj,jl), ' CAT no ', jl
903	WRITE(numout,) ' a_i_old ** ', old_a_i(ji,jj,jl), ' CAT no ', jl
904	WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
905	END DO
906	! WRITE(numout,*) ' CORRECTION BARBARE '
907	! z_da_ex = MAX( at_i(ji,jj) - zamax , 0.0 )
908	ENDIF
909	END DO
910	END DO
911	ENDIF
912
913	! Final thickness distribution rebinning
914	! --------------------------------------
915	IF( ln_nicep ) THEN
916	WRITE(numout,*) ' rebinning before'
917	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
918	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
919	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
920	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
921	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
922	ENDIF
923	!old version
924	! CALL lim_itd_th_reb(1,jpl)
925
926	DO jm = 1, jpm
927	jbnd1 = ice_cat_bounds(jm,1)
928	jbnd2 = ice_cat_bounds(jm,2)
929	IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
930	IF (ice_ncat_types(jm) .EQ. 1 ) THEN
931	ENDIF
932	END DO
933
934	IF( ln_nicep ) THEN
935	WRITE(numout,*) ' rebinning final'
936	WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
937	WRITE(numout,*) ' at_i ', at_i(jiindx,jjindx)
938	WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
939	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
940	WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
941	ENDIF
942
943	at_i(:,:) = 0.0
944	DO jl = 1, jpl
945	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
946	END DO
947
948	!------------------------------------------------------------------------------
949	! 2) Corrections to avoid wrong values \|
950	!------------------------------------------------------------------------------
951	! Ice drift
952	!------------
953
954	DO jj = 2, jpjm1
955	DO ji = fs_2, fs_jpim1
956	IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
957	IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj) = 0.0 ! right side
958	IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
959	IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj) = 0.0 ! upper side
960	IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji-1,jj) = 0.0 ! bottom side
961	ENDIF
962	END DO
963	END DO
964	!mask velocities
965	u_ice(:,:) = u_ice(:,:) * tmu(:,:)
966	v_ice(:,:) = v_ice(:,:) * tmv(:,:)
967	!lateral boundary conditions
968	CALL lbc_lnk( u_ice(:,:), 'U', -1. )
969	CALL lbc_lnk( v_ice(:,:), 'V', -1. )
970
971	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
972	! ALERTES
973	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
974
975	IF( ln_nicep ) THEN
976	DO jj = 1, jpj
977	DO ji = 1, jpi
978	DO jl = 1, jpl
979	! IF ((v_i(ji,jj,jl).NE.0.0).AND.(a_i(ji,jj,jl).EQ.0.0)) THEN
980	! WRITE(numout,*) ' lim_update : incompatible volume and concentration '
981	END DO ! jl
982
983	DO jl = 1, jpl
984	IF ( (a_i(ji,jj,jl).GT.1.0).OR.(at_i(ji,jj).GT.1.0) ) THEN
985	zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )
986	WRITE(numout,*) ' lim_update : a_i > 1 '
987	WRITE(numout,*) ' PAS BIEN ----> ALERTE !!! '
988	WRITE(numout,*) ' ~~~~~~~~~~ at_i ', at_i(ji,jj)
989	WRITE(numout,*) ' Point - category', ji, jj, jl
990	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
991	WRITE(numout,) ' a_i ** a_i_old ', a_i(ji,jj,jl), old_a_i(ji,jj,jl)
992	WRITE(numout,) ' v_i ** v_i_old ', v_i(ji,jj,jl), old_v_i(ji,jj,jl)
993	WRITE(numout,) ' ht_i ** ', v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi06)*zindb
994	WRITE(numout,*) ' hi_max(jl), hi_max(jl-1) ', hi_max(jl), hi_max(jl-1)
995	WRITE(numout,*) ' d_v_i_thd / trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
996	WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
997	ENDIF
998	END DO
999
1000	END DO !jj
1001	END DO !ji
1002
1003	WRITE(numout,*) ' TESTOSC1 ', tio_u(jiindx,jjindx), tio_v(jiindx,jjindx)
1004	WRITE(numout,*) ' TESTOSC2 ', u_ice(jiindx,jjindx), v_ice(jiindx,jjindx)
1005	WRITE(numout,*) ' TESTOSC3 ', u_oce(jiindx,jjindx), v_oce(jiindx,jjindx)
1006	WRITE(numout,*) ' TESTOSC4 ', utau (jiindx,jjindx), vtau (jiindx,jjindx)
1007	ENDIF
1008
1009
1010	IF(ln_ctl) THEN ! Control print
1011	CALL prt_ctl_info(' ')
1012	CALL prt_ctl_info(' - Cell values : ')
1013	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
1014	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_update : cell area :')
1015	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update : at_i :')
1016	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update : vt_i :')
1017	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_update : vt_s :')
1018	CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update : strength :')
1019	CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
1020	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 :')
1021	CALL prt_ctl(tab2d_1=old_u_ice , clinfo1=' lim_update : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice :')
1022
1023	DO jl = 1, jpl
1024	CALL prt_ctl_info(' ')
1025	CALL prt_ctl_info(' - Category : ', ivar1=jl)
1026	CALL prt_ctl_info(' ~~~~~~~~~~')
1027	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_update : ht_i : ')
1028	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_update : ht_s : ')
1029	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_update : t_su : ')
1030	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_update : t_snow : ')
1031	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_update : sm_i : ')
1032	CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update : o_i : ')
1033	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update : a_i : ')
1034	CALL prt_ctl(tab2d_1=old_a_i (:,:,jl) , clinfo1= ' lim_update : old_a_i : ')
1035	CALL prt_ctl(tab2d_1=d_a_i_trp (:,:,jl) , clinfo1= ' lim_update : d_a_i_trp : ')
1036	CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update : d_a_i_thd : ')
1037	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update : v_i : ')
1038	CALL prt_ctl(tab2d_1=old_v_i (:,:,jl) , clinfo1= ' lim_update : old_v_i : ')
1039	CALL prt_ctl(tab2d_1=d_v_i_trp (:,:,jl) , clinfo1= ' lim_update : d_v_i_trp : ')
1040	CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update : d_v_i_thd : ')
1041	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update : v_s : ')
1042	CALL prt_ctl(tab2d_1=old_v_s (:,:,jl) , clinfo1= ' lim_update : old_v_s : ')
1043	CALL prt_ctl(tab2d_1=d_v_s_trp (:,:,jl) , clinfo1= ' lim_update : d_v_s_trp : ')
1044	CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update : d_v_s_thd : ')
1045	CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : e_i1 : ')
1046	CALL prt_ctl(tab2d_1=old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : old_e_i1 : ')
1047	CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : de_i1_trp : ')
1048	CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : de_i1_thd : ')
1049	CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : e_i2 : ')
1050	CALL prt_ctl(tab2d_1=old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : old_e_i2 : ')
1051	CALL prt_ctl(tab2d_1=d_e_i_trp (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : de_i2_trp : ')
1052	CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update : de_i2_thd : ')
1053	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update : e_snow : ')
1054	CALL prt_ctl(tab2d_1=old_e_s (:,:,1,jl) , clinfo1= ' lim_update : old_e_snow : ')
1055	CALL prt_ctl(tab2d_1=d_e_s_trp (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : d_e_s_trp : ')
1056	CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update : d_e_s_thd : ')
1057	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update : smv_i : ')
1058	CALL prt_ctl(tab2d_1=old_smv_i (:,:,jl) , clinfo1= ' lim_update : old_smv_i : ')
1059	CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl) , clinfo1= ' lim_update : d_smv_i_trp : ')
1060	CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update : d_smv_i_thd : ')
1061	CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update : oa_i : ')
1062	CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update : old_oa_i : ')
1063	CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update : d_oa_i_trp : ')
1064	CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update : d_oa_i_thd : ')
1065	CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update : Path. case : ')
1066
1067	DO jk = 1, nlay_i
1068	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
1069	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update : t_i : ')
1070	END DO
1071	END DO
1072
1073	CALL prt_ctl_info(' ')
1074	CALL prt_ctl_info(' - Heat / FW fluxes : ')
1075	CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ')
1076	CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')
1077	CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ')
1078	CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
1079
1080	CALL prt_ctl_info(' ')
1081	CALL prt_ctl_info(' - Stresses : ')
1082	CALL prt_ctl_info(' ~~~~~~~~~~ ')
1083	CALL prt_ctl(tab2d_1=utau , clinfo1= ' lim_update : utau : ', tab2d_2=vtau , clinfo2= ' vtau : ')
1084	CALL prt_ctl(tab2d_1=utaui_ice , clinfo1= ' lim_update : utaui_ice : ', tab2d_2=vtaui_ice , clinfo2= ' vtaui_ice : ')
1085	CALL prt_ctl(tab2d_1=u_oce , clinfo1= ' lim_update : u_oce : ', tab2d_2=v_oce , clinfo2= ' v_oce : ')
1086	ENDIF
1087
1088	!---------------------
1089
1090	END SUBROUTINE lim_update
1091	#else
1092	!!----------------------------------------------------------------------
1093	!! Default option Empty Module No sea-ice model
1094	!!----------------------------------------------------------------------
1095	CONTAINS
1096	SUBROUTINE lim_update ! Empty routine
1097	END SUBROUTINE lim_update
1098
1099	#endif
1100
1101	END MODULE limupdate

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

Download in other formats: