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

Last change on this file since 1656 was 1469, checked in by smasson, 15 years ago
[uv]taui_ice renamed [uv]tau_ice, see ticket:452
Property svn:keywords set to `Id`
File size: 47.9 KB

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

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

Download in other formats: