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/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90 @ 2528

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

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

Download in other formats: