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

Last change on this file since 1264 was 1156, checked in by rblod, 16 years ago
Update Id and licence information, see ticket #210
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
	26	USE ice_oce ! ice variables
[888]	27	USE sbc_oce ! Surface boundary condition: ocean fields
	28	USE sbc_ice ! Surface boundary condition: ice fields
[825]	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
[888]	38	USE limsbc
[825]	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
[863]	46	USE prtctl ! Print control
[825]	47
	48
	49	IMPLICIT NONE
	50	PRIVATE
	51
	52	!! * Accessibility
	53	PUBLIC lim_update ! routine called by ice_step
	54
[895]	55	!! * Substitutions
	56	# include "vectopt_loop_substitute.h90"
	57
[825]	58	!!----------------------------------------------------------------------
[1156]	59	!! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008)
	60	!! $Id$
	61	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]	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 :: &
[921]	84	ji, jj, & ! geographical indices
	85	jk, jl, jm ! layer, category and type indices
[825]	86	INTEGER :: &
[921]	87	jbnd1, jbnd2
[825]	88	INTEGER :: &
[921]	89	i_ice_switch
[825]	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
[921]	101	! in an ice-free cell
[825]	102	REAL(wp) :: & ! dummy switches and arguments
[834]	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, &
[825]	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
[834]	113
[825]	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
[921]	118	! internal melt
[825]	119	REAL(wp) :: &
	120	zweight, zesum
	121
[921]	122
[825]	123	!!-------------------------------------------------------------------
	124
[921]	125	IF( ln_nicep ) THEN
	126	WRITE(numout,*) ' lim_update '
	127	WRITE(numout,*) ' ~~~~~~~~~~ '
[825]	128
[921]	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
[825]	139
[921]	140	!------------------------------------------------------------------------------
	141	! 1. Update of Global variables \|
	142	!------------------------------------------------------------------------------
[825]	143
[921]	144	!---------------------
	145	! Ice dynamics
	146	!---------------------
[825]	147
[921]	148	u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
	149	v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
[825]	150
[921]	151	!-----------------------------
	152	! Update ice and snow volumes
	153	!-----------------------------
[825]	154
[921]	155	DO jl = 1, jpl
	156	DO jj = 1, jpj
	157	DO ji = 1, jpi
[825]	158
[921]	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
[825]	166
[921]	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)
[825]	176
[921]	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
[825]	201
[921]	202	ENDIF
	203	END DO
	204	END DO
	205	END DO
[825]	206
[921]	207	!--------------------
	208	! Excessive ablation
	209	!--------------------
[825]	210
[921]	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)
[825]	229
[921]	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))
[825]	242
[921]	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
[825]	247
[921]	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
[825]	258
[921]	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
[825]	268
[921]	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
[825]	280
[921]	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
[825]	291
[921]	292	!---------------------------------------------
	293	! Ice concentration and ice heat content
	294	!---------------------------------------------
[825]	295
[921]	296	a_i (:,:,:) = a_i (:,:,:) + d_a_i_trp(:,:,:) &
	297	+ d_a_i_thd(:,:,:)
	298	CALL lim_var_glo2eqv ! useless, just for debug
[1055]	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
[921]	304	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)
	305	CALL lim_var_glo2eqv ! useless, just for debug
[1055]	306	IF( ln_nicep) THEN
[921]	307	WRITE(numout,*) ' After transport update '
[1055]	308	DO jk = 1, nlay_i
	309	WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
	310	END DO
	311	ENDIF
[921]	312	e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)
	313	CALL lim_var_glo2eqv ! useless, just for debug
[1055]	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
[825]	320
[921]	321	at_i(:,:) = 0.0
	322	DO jl = 1, jpl
	323	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
	324	END DO
[825]	325
[921]	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
[825]	339
[921]	340	!------------------------------
	341	! Snow temperature and ice age
	342	!------------------------------
[825]	343
[921]	344	e_s(:,:,:,:) = e_s(:,:,:,:) + &
	345	d_e_s_trp(:,:,:,:) + &
	346	d_e_s_thd(:,:,:,:)
[825]	347
[921]	348	oa_i(:,:,:) = oa_i(:,:,:) + &
	349	d_oa_i_trp(:,:,:) + &
	350	d_oa_i_thd(:,:,:)
[825]	351
[921]	352	!--------------
	353	! Ice salinity
	354	!--------------
[825]	355
[921]	356	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
[825]	357
[921]	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
[825]	367
[921]	368	smv_i(:,:,:) = smv_i(:,:,:) + &
	369	d_smv_i_thd(:,:,:) + &
	370	d_smv_i_trp(:,:,:)
[825]	371
[921]	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
[825]	377
[921]	378	ENDIF ! num_sal .EQ. 2
[825]	379
[921]	380	CALL lim_var_glo2eqv
[825]	381
[921]	382	!--------------------------------------
	383	! 2. Review of all pathological cases
	384	!--------------------------------------
[825]	385
[921]	386	zrtt = 173.15 * rone
	387	zacrith = 1.0e-6
[825]	388
[921]	389	!-------------------------------------------
	390	! 2.1) Advection of ice in an ice-free cell
	391	!-------------------------------------------
	392	! should be removed since it is treated after dynamics now
[825]	393
[921]	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
[825]	410
[921]	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
[825]	422
[921]	423	!change this 14h44
	424	zhimax = 20.0 ! line added up
	425	! change this also 17 aug
	426	zhimax = 30.0 ! line added up
[825]	427
[921]	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,:))
[825]	449
[921]	450	END DO ! ji
	451	END DO !jj
	452	END DO !jl
[825]	453
[921]	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
[825]	465
[921]	466	at_i(:,:) = 0.0
	467	DO jl = 1, jpl
	468	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
	469	END DO
[825]	470
[921]	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
[825]	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
[921]	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
[825]	505
[921]	506	at_i(:,:) = 0.0
	507	DO jl = 1, jpl
	508	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
	509	END DO
[825]	510
[921]	511	!---------------------------------
	512	! 2.3) Melt of an internal layer
	513	!---------------------------------
	514	internal_melt(:,:,:) = .false.
[825]	515
[921]	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)
[825]	529	internal_melt(ji,jj,jl) = .true.
[921]	530	ENDIF
	531	END DO ! ji
	532	END DO ! jj
	533	END DO !jk
	534	END DO !jl
[825]	535
[921]	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)
[825]	547
[921]	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)
[825]	563
[921]	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
[825]	568
[921]	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
[825]	582
[921]	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)
[825]	587
[921]	588	! Redistributing energy on the new grid
	589	IF ( ind_im .GT. 0 ) THEN
[825]	590
[921]	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 )
[825]	598
[921]	599	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
	600	END DO !layer
	601	END DO ! jk
[825]	602
[921]	603	zesum = 0.0
	604	DO jk = 1, nlay_i
	605	zesum = zesum + e_i(ji,jj,jk,jl)
	606	END DO
[825]	607
[921]	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
[825]	611
[921]	612	ELSE ! ind_im .EQ. 0, total melt
	613	e_i(ji,jj,jk,jl) = 0.0
	614	ENDIF
[825]	615
[921]	616	ENDIF ! internal_melt
[825]	617
[921]	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
[825]	633
[921]	634	internal_melt(:,:,:) = .false.
[825]	635
[921]	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
[825]	644
[921]	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.
[825]	648
[921]	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
[825]	661
[921]	662	END DO
	663	END DO
	664	END DO
[825]	665
[921]	666	DO jl = 1, jpl
	667	DO jj = 1, jpj
	668	DO ji = 1, jpi
	669	IF ( internal_melt(ji,jj,jl) ) THEN
[825]	670	v_s(ji,jj,jl) = 0.0
	671	e_s(ji,jj,1,jl) = 0.0
[921]	672	! ! release heat
[825]	673	fheat_res(ji,jj) = fheat_res(ji,jj) &
[921]	674	+ ze_s * v_s(ji,jj,jl) / rdt_ice
[825]	675	! release mass
	676	rdmsnif(ji,jj) = rdmsnif(ji,jj) + rhosn * v_s(ji,jj,jl)
[921]	677	ENDIF
	678	END DO
	679	END DO
	680	END DO
[825]	681
[921]	682	zbigvalue = 1.0d+20
[825]	683
[921]	684	DO jl = 1, jpl
	685	DO jj = 1, jpj
	686	DO ji = 1, jpi
[825]	687
[921]	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
[825]	695
[921]	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)
[825]	705
[921]	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 ) )
[825]	713
[921]	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)
[825]	720
[921]	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
[825]	727
[921]	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
[825]	734
[921]	735	!-------------------------
	736	! 2.8) Snow heat content
	737	!-------------------------
[825]	738
[921]	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
[825]	742
[921]	743	END DO ! ji
	744	END DO ! jj
	745	END DO ! jl
[825]	746
[921]	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
[825]	758
[921]	759	!------------------------
	760	! 2.9) Ice heat content
	761	!------------------------
[825]	762
[921]	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
[825]	776
[921]	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)
[825]	783
[921]	784	WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
	785	ENDIF
[825]	786
[921]	787	!---------------------
	788	! 2.11) Ice salinity
	789	!---------------------
[825]	790
[921]	791	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
[825]	792
[921]	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
[825]	807
[921]	808	ENDIF
[825]	809
[921]	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
[825]	818
[921]	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
[825]	846
[921]	847	!--- 2.13 Total ice concentration should not exceed 1
	848	!-----------------------------------------------------
	849	zamax = amax
	850	! 2.13.1) individual concentrations cannot exceed zamax
	851	!------------------------------------------------------
[825]	852
[921]	853	at_i(:,:) = 0.0
	854	DO jl = 1, jpl
	855	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
	856	END DO
[825]	857
[921]	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
[825]	864
[921]	865	DO jj = 1, jpj
	866	DO ji = 1, jpi
[825]	867
[921]	868	! 0) Excessive area ?
	869	z_da_ex = MAX( at_i(ji,jj) - zamax , 0.0 )
[825]	870
[921]	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)
[825]	879
[921]	880	END DO
[825]	881
[921]	882	END DO !ji
	883	END DO !jj
[825]	884
[921]	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
[825]	893
[921]	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
[825]	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
[921]	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
[825]	948
[921]	949	at_i(:,:) = 0.0
	950	DO jl = 1, jpl
	951	at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
	952	END DO
[869]	953
[921]	954	!------------------------------------------------------------------------------
	955	! 2) Corrections to avoid wrong values \|
	956	!------------------------------------------------------------------------------
	957	! Ice drift
	958	!------------
	959
[894]	960	DO jj = 2, jpjm1
	961	DO ji = fs_2, fs_jpim1
[825]	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
[921]	977	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	978	! ALERTES
	979	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[825]	980
[921]	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
[825]	988
[921]	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
[825]	1005
[921]	1006	END DO !jj
	1007	END DO !ji
[825]	1008
[921]	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
[825]	1014
	1015
[863]	1016	IF(ln_ctl) THEN ! Control print
[867]	1017	CALL prt_ctl_info(' ')
[863]	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
[867]	1030	CALL prt_ctl_info(' ')
[863]	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)
[867]	1075	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update : t_i : ')
[863]	1076	END DO
	1077	END DO
	1078
[867]	1079	CALL prt_ctl_info(' ')
[863]	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 : ')
[888]	1083	CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ')
[863]	1084	CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ')
	1085
[867]	1086	CALL prt_ctl_info(' ')
[863]	1087	CALL prt_ctl_info(' - Stresses : ')
	1088	CALL prt_ctl_info(' ~~~~~~~~~~ ')
[888]	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 : ')
[863]	1092	ENDIF
	1093
[921]	1094	!---------------------
[863]	1095
[825]	1096	END SUBROUTINE lim_update
[834]	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
[825]	1107	END MODULE limupdate

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

Download in other formats: