New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limupdate.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 2293

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

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

Download in other formats: