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.

limitd_th.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5146

Last change on this file since 5146 was 5134, checked in by clem, 9 years ago
LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)
Property svn:keywords set to `Id`
File size: 39.1 KB

Rev	Line
[825]	1	MODULE limitd_th
	2	!!======================================================================
	3	!! * MODULE limitd_th *
[3625]	4	!! LIM3 ice model : ice thickness distribution: Thermodynamics
[825]	5	!!======================================================================
[2715]	6	!! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
	7	!! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3
[4869]	8	!! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age
[2715]	9	!! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked
	10	!!----------------------------------------------------------------------
[2528]	11	#if defined key_lim3
[825]	12	!!----------------------------------------------------------------------
[2528]	13	!! 'key_lim3' : LIM3 sea-ice model
	14	!!----------------------------------------------------------------------
[2715]	15	!! lim_itd_th_rem :
	16	!! lim_itd_th_reb :
	17	!! lim_itd_fitline :
	18	!! lim_itd_shiftice :
[2528]	19	!!----------------------------------------------------------------------
[4161]	20	USE dom_ice ! LIM-3 domain
	21	USE par_oce ! ocean parameters
	22	USE dom_oce ! ocean domain
	23	USE phycst ! physical constants (ocean directory)
	24	USE thd_ice ! LIM-3 thermodynamic variables
	25	USE ice ! LIM-3 variables
	26	USE limvar ! LIM-3 variables
	27	USE prtctl ! Print control
	28	USE in_out_manager ! I/O manager
	29	USE lib_mpp ! MPP library
	30	USE wrk_nemo ! work arrays
	31	USE lib_fortran ! to use key_nosignedzero
[5123]	32	USE limcons ! conservation tests
[921]	33
[825]	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC lim_itd_th_rem
	38	PUBLIC lim_itd_th_reb
[825]	39
	40	!!----------------------------------------------------------------------
[4161]	41	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]	42	!! $Id$
[2715]	43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	44	!!----------------------------------------------------------------------
	45	CONTAINS
	46
[4869]	47	SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
[921]	48	!!------------------------------------------------------------------
	49	!! * ROUTINE lim_itd_th_rem *
	50	!!
[2715]	51	!! ** Purpose : computes the redistribution of ice thickness
	52	!! after thermodynamic growth of ice thickness
	53	!!
[921]	54	!! ** Method : Linear remapping
	55	!!
[2715]	56	!! References : W.H. Lipscomb, JGR 2001
[921]	57	!!------------------------------------------------------------------
[2715]	58	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	59	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	60	INTEGER , INTENT (in) :: kt ! Ocean time step
	61	!
	62	INTEGER :: ji, jj, jl ! dummy loop index
[4161]	63	INTEGER :: ii, ij ! 2D corresponding indices to ji
	64	INTEGER :: nd ! local integer
[2715]	65	REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars
[4161]	66	REAL(wp) :: zx2, zwk2, zda0, zetamax ! - -
[5123]	67	REAL(wp) :: zx3
[2715]	68	CHARACTER (len = 15) :: fieldid
[825]	69
[3294]	70	INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index
[825]	71
[3294]	72	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdhice ! ice thickness increment
	73	REAL(wp), POINTER, DIMENSION(:,:,:) :: g0 ! coefficients for fitting the line of the ITD
	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: g1 ! coefficients for fitting the line of the ITD
	75	REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness
	76	REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness
[4872]	77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness
[3294]	78	REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es
	79	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume
	80	REAL(wp), POINTER, DIMENSION(:) :: zvetamin, zvetamax ! maximum values for etas
	81	INTEGER , POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions
	82	INTEGER :: nbrem ! number of cells with ice to transfer
	83	REAL(wp) :: zslope ! used to compute local thermodynamic "speeds"
	84	REAL(wp), POINTER, DIMENSION(:,:) :: zhb0, zhb1 ! category boundaries for thinnes categories
	85	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
	86	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
	87	REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories
	88	REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init, et_s_final ! snow energy summed over categories
	89	INTEGER , POINTER, DIMENSION(:,:) :: zremap_flag ! compute remapping or not ????
	90	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhbnew ! new boundaries of ice categories
	91	!!------------------------------------------------------------------
[825]	92
[3294]	93	CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer
	94	CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer
[4872]	95	CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
[3294]	96	CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )
	97	CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )
	98	CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )
	99	CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	100	CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
[825]	101
[921]	102	!!----------------------------------------------------------------------------------------------
	103	!! 0) Conservation checkand changes in each ice category
	104	!!----------------------------------------------------------------------------------------------
[2715]	105	IF( con_i ) THEN
[825]	106	CALL lim_column_sum (jpl, v_i, vt_i_init)
	107	CALL lim_column_sum (jpl, v_s, vt_s_init)
	108	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
	109	dummy_es(:,:,:) = e_s(:,:,1,:)
	110	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
	111	ENDIF
	112
[921]	113	!!----------------------------------------------------------------------------------------------
	114	!! 1) Compute thickness and changes in each ice category
	115	!!----------------------------------------------------------------------------------------------
[2715]	116	IF( kt == nit000 .AND. lwp) THEN
[921]	117	WRITE(numout,*)
	118	WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'
	119	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	120	WRITE(numout,*) ' klbnd : ', klbnd
	121	WRITE(numout,*) ' kubnd : ', kubnd
	122	ENDIF
[825]	123
[2715]	124	zdhice(:,:,:) = 0._wp
[921]	125	DO jl = klbnd, kubnd
	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
[5134]	128	rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes
[4990]	129	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
[5134]	130	rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) !0 if no ice and 1 if yes
[4990]	131	zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
[5134]	132	!clem IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)
	133	zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)
[921]	134	END DO
	135	END DO
	136	END DO
	137
	138	!-----------------------------------------------------------------------------------------------
	139	! 2) Compute fractional ice area in each grid cell
	140	!-----------------------------------------------------------------------------------------------
[2715]	141	at_i(:,:) = 0._wp
[825]	142	DO jl = klbnd, kubnd
[2715]	143	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[825]	144	END DO
	145
[921]	146	!-----------------------------------------------------------------------------------------------
	147	! 3) Identify grid cells with ice
	148	!-----------------------------------------------------------------------------------------------
[825]	149	nbrem = 0
	150	DO jj = 1, jpj
	151	DO ji = 1, jpi
[5123]	152	IF ( at_i(ji,jj) > epsi10 ) THEN
[825]	153	nbrem = nbrem + 1
	154	nind_i(nbrem) = ji
	155	nind_j(nbrem) = jj
[3294]	156	zremap_flag(ji,jj) = 1
[825]	157	ELSE
[3294]	158	zremap_flag(ji,jj) = 0
[825]	159	ENDIF
[5123]	160	END DO
	161	END DO
[825]	162
[921]	163	!-----------------------------------------------------------------------------------------------
	164	! 4) Compute new category boundaries
	165	!-----------------------------------------------------------------------------------------------
[825]	166	!- 4.1 Compute category boundaries
[2715]	167	zhbnew(:,:,:) = 0._wp
[825]	168
	169	DO jl = klbnd, kubnd - 1
	170	DO ji = 1, nbrem
[4161]	171	ii = nind_i(ji)
	172	ij = nind_j(ji)
[825]	173	!
[4688]	174	zhbnew(ii,ij,jl) = hi_max(jl)
[4872]	175	IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
[825]	176	!interpolate between adjacent category growth rates
[4872]	177	zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
	178	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
	179	ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN
[4161]	180	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
[4872]	181	ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN
[4161]	182	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
[825]	183	ENDIF
[2715]	184	END DO
[825]	185
[921]	186	!- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
[825]	187	DO ji = 1, nbrem
[4161]	188	ii = nind_i(ji)
	189	ij = nind_j(ji)
[4688]	190	IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN
[4161]	191	zremap_flag(ii,ij) = 0
[4688]	192	ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN
[4161]	193	zremap_flag(ii,ij) = 0
[825]	194	ENDIF
	195
[921]	196	!- 4.3 Check that each zhbnew does not exceed maximal values hi_max
[4688]	197	IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
	198	IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
	199	END DO
	200
[5123]	201	END DO
[825]	202
[921]	203	!-----------------------------------------------------------------------------------------------
	204	! 5) Identify cells where ITD is to be remapped
	205	!-----------------------------------------------------------------------------------------------
	206	nbrem = 0
	207	DO jj = 1, jpj
	208	DO ji = 1, jpi
[4688]	209	IF( zremap_flag(ji,jj) == 1 ) THEN
[921]	210	nbrem = nbrem + 1
	211	nind_i(nbrem) = ji
	212	nind_j(nbrem) = jj
	213	ENDIF
[4688]	214	END DO
	215	END DO
[825]	216
[921]	217	!-----------------------------------------------------------------------------------------------
	218	! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
	219	!-----------------------------------------------------------------------------------------------
	220	DO jj = 1, jpj
	221	DO ji = 1, jpi
[4869]	222	zhb0(ji,jj) = hi_max(0) ! 0eme
	223	zhb1(ji,jj) = hi_max(1) ! 1er
[825]	224
[2715]	225	zhbnew(ji,jj,klbnd-1) = 0._wp
[825]	226
[2715]	227	IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
[5134]	228	zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
[921]	229	ELSE
[4293]	230	zhbnew(ji,jj,kubnd) = hi_max(kubnd)
	231	!!? clem bug: since hi_max(jpl)=99, this limit is very high
	232	!!? but I think it is erased in fitline subroutine
[921]	233	ENDIF
[825]	234
[5134]	235	END DO
	236	END DO
[825]	237
[921]	238	!-----------------------------------------------------------------------------------------------
	239	! 7) Compute g(h)
	240	!-----------------------------------------------------------------------------------------------
	241	!- 7.1 g(h) for category 1 at start of time step
[5123]	242	CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
[2715]	243	& hR(:,:,klbnd), zremap_flag )
[825]	244
[921]	245	!- 7.2 Area lost due to melting of thin ice (first category, klbnd)
	246	DO ji = 1, nbrem
[4161]	247	ii = nind_i(ji)
	248	ij = nind_j(ji)
[825]	249
[5123]	250	IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
[4161]	251	zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
[5123]	252	IF( zdh0 < 0.0 ) THEN !remove area from category 1
	253	zdh0 = MIN( -zdh0, hi_max(klbnd) )
[825]	254
[921]	255	!Integrate g(1) from 0 to dh0 to estimate area melted
[5123]	256	zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
	257	IF( zetamax > 0.0 ) THEN
[921]	258	zx1 = zetamax
[5123]	259	zx2 = 0.5 * zetamax * zetamax
[4161]	260	zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed
[921]	261	! Constrain new thickness <= ht_i
[5123]	262	zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0
[921]	263	!ice area lost due to melting of thin ice
[5123]	264	zda0 = MIN( zda0, zdamax )
[825]	265
[921]	266	! Remove area, conserving volume
[5123]	267	ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
[4161]	268	a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0
[5123]	269	v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
	270	ENDIF
[825]	271
[5134]	272	ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0
[5123]	273	zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )
[921]	274	! zhbnew was 0, and is shifted to the right to account for thin ice
	275	! growth in openwater (F0 = f1)
	276	ENDIF ! zdh0
[825]	277
[2715]	278	ENDIF ! a_i > epsi10
[825]	279
[5134]	280	END DO
[825]	281
[921]	282	!- 7.3 g(h) for each thickness category
	283	DO jl = klbnd, kubnd
[5123]	284	CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
	285	& g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
[921]	286	END DO
[825]	287
[921]	288	!-----------------------------------------------------------------------------------------------
	289	! 8) Compute area and volume to be shifted across each boundary
	290	!-----------------------------------------------------------------------------------------------
[825]	291
[921]	292	DO jl = klbnd, kubnd - 1
	293	DO jj = 1, jpj
	294	DO ji = 1, jpi
	295	zdonor(ji,jj,jl) = 0
	296	zdaice(ji,jj,jl) = 0.0
	297	zdvice(ji,jj,jl) = 0.0
	298	END DO
	299	END DO
[825]	300
[921]	301	DO ji = 1, nbrem
[4161]	302	ii = nind_i(ji)
	303	ij = nind_j(ji)
[825]	304
[5123]	305	IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
[825]	306
[921]	307	! left and right integration limits in eta space
[5123]	308	zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
	309	zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
[4161]	310	zdonor(ii,ij,jl) = jl
[825]	311
[921]	312	ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
[825]	313
[921]	314	! left and right integration limits in eta space
	315	zvetamin(ji) = 0.0
[5123]	316	zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
[4161]	317	zdonor(ii,ij,jl) = jl + 1
[825]	318
[921]	319	ENDIF ! zhbnew(jl) > hi_max(jl)
[825]	320
[5123]	321	zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
[921]	322	zetamin = zvetamin(ji)
[825]	323
[921]	324	zx1 = zetamax - zetamin
[5123]	325	zwk1 = zetamin * zetamin
	326	zwk2 = zetamax * zetamax
	327	zx2 = 0.5 * ( zwk2 - zwk1 )
[921]	328	zwk1 = zwk1 * zetamin
	329	zwk2 = zwk2 * zetamax
[5123]	330	zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 )
[4161]	331	nd = zdonor(ii,ij,jl)
	332	zdaice(ii,ij,jl) = g1(ii,ij,nd)zx2 + g0(ii,ij,nd)zx1
[4688]	333	zdvice(ii,ij,jl) = g1(ii,ij,nd)zx3 + g0(ii,ij,nd)zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
[921]	334
[5123]	335	END DO
[921]	336	END DO ! jl klbnd -> kubnd - 1
	337
	338	!!----------------------------------------------------------------------------------------------
	339	!! 9) Shift ice between categories
	340	!!----------------------------------------------------------------------------------------------
	341	CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
	342
	343	!!----------------------------------------------------------------------------------------------
	344	!! 10) Make sure ht_i >= minimum ice thickness hi_min
	345	!!----------------------------------------------------------------------------------------------
	346
	347	DO ji = 1, nbrem
[4161]	348	ii = nind_i(ji)
	349	ij = nind_j(ji)
[5123]	350	IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
[5134]	351	a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
[5123]	352	ht_i(ii,ij,1) = rn_himin
[921]	353	ENDIF
[5123]	354	END DO
[921]	355
	356	!!----------------------------------------------------------------------------------------------
	357	!! 11) Conservation check
	358	!!----------------------------------------------------------------------------------------------
[825]	359	IF ( con_i ) THEN
	360	CALL lim_column_sum (jpl, v_i, vt_i_final)
	361	fieldid = ' v_i : limitd_th '
	362	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
	363
	364	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_final)
	365	fieldid = ' e_i : limitd_th '
	366	CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)
	367
	368	CALL lim_column_sum (jpl, v_s, vt_s_final)
	369	fieldid = ' v_s : limitd_th '
	370	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	371
	372	dummy_es(:,:,:) = e_s(:,:,1,:)
	373	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
	374	fieldid = ' e_s : limitd_th '
	375	CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
	376	ENDIF
	377
[3294]	378	CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer
	379	CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer
[4872]	380	CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
[3294]	381	CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )
	382	CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )
	383	CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )
	384	CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	385	CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
	386
[921]	387	END SUBROUTINE lim_itd_th_rem
[825]	388
	389
[5123]	390	SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
[921]	391	!!------------------------------------------------------------------
	392	!! * ROUTINE lim_itd_fitline *
	393	!!
[2715]	394	!! ** Purpose : fit g(h) with a line using area, volume constraints
[921]	395	!!
[2715]	396	!! ** Method : Fit g(h) with a line, satisfying area and volume constraints.
	397	!! To reduce roundoff errors caused by large values of g0 and g1,
	398	!! we actually compute g(eta), where eta = h - hL, and hL is the
	399	!! left boundary.
[921]	400	!!------------------------------------------------------------------
[2715]	401	INTEGER , INTENT(in ) :: num_cat ! category index
	402	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: HbL, HbR ! left and right category boundaries
	403	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hice ! ice thickness
	404	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: g0, g1 ! coefficients in linear equation for g(eta)
	405	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hL ! min value of range over which g(h) > 0
	406	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hR ! max value of range over which g(h) > 0
[3294]	407	INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag !
[2715]	408	!
	409	INTEGER :: ji,jj ! horizontal indices
	410	REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL)
	411	REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL)
	412	REAL(wp) :: zdhr ! 1 / (hR - hL)
	413	REAL(wp) :: zwk1, zwk2 ! temporary variables
	414	!!------------------------------------------------------------------
	415	!
	416	!
[825]	417	DO jj = 1, jpj
	418	DO ji = 1, jpi
[2715]	419	!
[4333]	420	IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 &
[3294]	421	& .AND. hice(ji,jj) > 0._wp ) THEN
[825]	422
[921]	423	! Initialize hL and hR
	424
[825]	425	hL(ji,jj) = HbL(ji,jj)
	426	hR(ji,jj) = HbR(ji,jj)
	427
[921]	428	! Change hL or hR if hice falls outside central third of range
[825]	429
[5123]	430	zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
	431	zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
[825]	432
[2715]	433	IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
	434	ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
[825]	435	ENDIF
	436
[921]	437	! Compute coefficients of g(eta) = g0 + g1*eta
	438
[2715]	439	zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
	440	zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
	441	zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
[5123]	442	g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
	443	g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
[2715]	444	!
	445	ELSE ! remap_flag = .false. or a_i < epsi10
	446	hL(ji,jj) = 0._wp
	447	hR(ji,jj) = 0._wp
	448	g0(ji,jj) = 0._wp
	449	g1(ji,jj) = 0._wp
	450	ENDIF ! a_i > epsi10
	451	!
	452	END DO
	453	END DO
	454	!
	455	END SUBROUTINE lim_itd_fitline
[825]	456
	457
[2715]	458	SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	459	!!------------------------------------------------------------------
	460	!! * ROUTINE lim_itd_shiftice *
[2715]	461	!!
	462	!! ** Purpose : shift ice across category boundaries, conserving everything
[921]	463	!! ( area, volume, energy, age*vol, and mass of salt )
	464	!!
	465	!! ** Method :
	466	!!------------------------------------------------------------------
[3294]	467	INTEGER , INTENT(in ) :: klbnd ! Start thickness category index point
	468	INTEGER , INTENT(in ) :: kubnd ! End point on which the the computation is applied
[2715]	469	INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in ) :: zdonor ! donor category index
	470	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdaice ! ice area transferred across boundary
	471	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdvice ! ice volume transferred across boundary
[825]	472
[2715]	473	INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices
[4161]	474	INTEGER :: ii, ij ! indices when changing from 2D-1D is done
[825]	475
[3294]	476	REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn
	477	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
[825]	478
[2715]	479	REAL(wp) :: zdvsnow, zdesnow ! snow volume and energy transferred
	480	REAL(wp) :: zdeice ! ice energy transferred
	481	REAL(wp) :: zdsm_vice ! ice salinity times volume transferred
	482	REAL(wp) :: zdo_aice ! ice age times volume transferred
	483	REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred
[825]	484
[3294]	485	INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions
[825]	486
[2715]	487	INTEGER :: nbrem ! number of cells with ice to transfer
[825]	488
[2715]	489	LOGICAL :: zdaice_negative ! true if daice < -puny
	490	LOGICAL :: zdvice_negative ! true if dvice < -puny
	491	LOGICAL :: zdaice_greater_aicen ! true if daice > aicen
	492	LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen
	493	!!------------------------------------------------------------------
[825]	494
[3294]	495	CALL wrk_alloc( jpi,jpj,jpl, zaTsfn )
	496	CALL wrk_alloc( jpi,jpj, zworka )
	497	CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	498
[921]	499	!----------------------------------------------------------------------------------------------
	500	! 1) Define a variable equal to a_i*T_su
	501	!----------------------------------------------------------------------------------------------
[825]	502
	503	DO jl = klbnd, kubnd
[5123]	504	zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
[2715]	505	END DO
[825]	506
[5134]	507	!clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10)
	508	! !----------------------------------------------------------------------------------------------
	509	! ! 2) Check for daice or dvice out of range, allowing for roundoff error
	510	! !----------------------------------------------------------------------------------------------
	511	! ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
	512	! ! has a small area, with h(n) very close to a boundary. Then
	513	! ! the coefficients of g(h) are large, and the computed daice and
	514	! ! dvice can be in error. If this happens, it is best to transfer
	515	! ! either the entire category or nothing at all, depending on which
	516	! ! side of the boundary hice(n) lies.
	517	! !-----------------------------------------------------------------
	518	! DO jl = klbnd, kubnd-1
	519	!
	520	! zdaice_negative = .false.
	521	! zdvice_negative = .false.
	522	! zdaice_greater_aicen = .false.
	523	! zdvice_greater_vicen = .false.
	524	!
	525	! DO jj = 1, jpj
	526	! DO ji = 1, jpi
	527	!
	528	! IF (zdonor(ji,jj,jl) > 0) THEN
	529	! jl1 = zdonor(ji,jj,jl)
	530	!
	531	! IF (zdaice(ji,jj,jl) < 0.0) THEN
	532	! IF (zdaice(ji,jj,jl) > -epsi10) THEN
	533	! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &
	534	! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
	535	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
	536	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	537	! ELSE
	538	! zdaice(ji,jj,jl) = 0.0 ! shift no ice
	539	! zdvice(ji,jj,jl) = 0.0
	540	! ENDIF
	541	! ELSE
	542	! zdaice_negative = .true.
	543	! ENDIF
	544	! ENDIF
	545	!
	546	! IF (zdvice(ji,jj,jl) < 0.0) THEN
	547	! IF (zdvice(ji,jj,jl) > -epsi10 ) THEN
	548	! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &
	549	! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
	550	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
	551	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	552	! ELSE
	553	! zdaice(ji,jj,jl) = 0.0 ! shift no ice
	554	! zdvice(ji,jj,jl) = 0.0
	555	! ENDIF
	556	! ELSE
	557	! zdvice_negative = .true.
	558	! ENDIF
	559	! ENDIF
	560	!
	561	! ! If daice is close to aicen, set daice = aicen.
	562	! IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN
	563	! IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN
	564	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
	565	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	566	! ELSE
	567	! zdaice_greater_aicen = .true.
	568	! ENDIF
	569	! ENDIF
	570	!
	571	! IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN
	572	! IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN
	573	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
	574	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	575	! ELSE
	576	! zdvice_greater_vicen = .true.
	577	! ENDIF
	578	! ENDIF
	579	!
	580	! ENDIF ! donor > 0
	581	! END DO
	582	! END DO
	583	!
	584	! END DO
	585	!clem
[921]	586	!-------------------------------------------------------------------------------
	587	! 3) Transfer volume and energy between categories
	588	!-------------------------------------------------------------------------------
[825]	589
	590	DO jl = klbnd, kubnd - 1
	591	nbrem = 0
	592	DO jj = 1, jpj
	593	DO ji = 1, jpi
[5123]	594	IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
[825]	595	nbrem = nbrem + 1
	596	nind_i(nbrem) = ji
	597	nind_j(nbrem) = jj
[5123]	598	ENDIF
[825]	599	END DO
	600	END DO
	601
	602	DO ji = 1, nbrem
[4161]	603	ii = nind_i(ji)
	604	ij = nind_j(ji)
[825]	605
[4161]	606	jl1 = zdonor(ii,ij,jl)
[5134]	607	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) )
	608	zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch
[2715]	609	IF( jl1 == jl) THEN ; jl2 = jl1+1
[5123]	610	ELSE ; jl2 = jl
[825]	611	ENDIF
	612
	613	!--------------
	614	! Ice areas
	615	!--------------
	616
[4161]	617	a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
	618	a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
[825]	619
	620	!--------------
	621	! Ice volumes
	622	!--------------
	623
[4161]	624	v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)
	625	v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
[825]	626
	627	!--------------
	628	! Snow volumes
	629	!--------------
	630
[4688]	631	zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij)
[4161]	632	v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
	633	v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow
[825]	634
	635	!--------------------
	636	! Snow heat content
	637	!--------------------
	638
[4688]	639	zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij)
[4161]	640	e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow
	641	e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow
[825]	642
	643	!--------------
	644	! Ice age
	645	!--------------
	646
[4688]	647	zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	648	oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice
	649	oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice
[825]	650
	651	!--------------
	652	! Ice salinity
	653	!--------------
	654
[4688]	655	zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij)
[4161]	656	smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice
	657	smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice
[825]	658
	659	!---------------------
	660	! Surface temperature
	661	!---------------------
	662
[4688]	663	zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	664	zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf
	665	zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf
[825]	666
[5123]	667	END DO
[825]	668
	669	!------------------
	670	! Ice heat content
	671	!------------------
	672
	673	DO jk = 1, nlay_i
	674	DO ji = 1, nbrem
[4161]	675	ii = nind_i(ji)
	676	ij = nind_j(ji)
[825]	677
[4161]	678	jl1 = zdonor(ii,ij,jl)
[5123]	679	IF (jl1 == jl) THEN
[825]	680	jl2 = jl+1
	681	ELSE ! n1 = n+1
	682	jl2 = jl
	683	ENDIF
	684
[4161]	685	zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
	686	e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice
	687	e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice
[5123]	688	END DO
	689	END DO
[825]	690
	691	END DO ! boundaries, 1 to ncat-1
	692
	693	!-----------------------------------------------------------------
	694	! Update ice thickness and temperature
	695	!-----------------------------------------------------------------
	696
	697	DO jl = klbnd, kubnd
	698	DO jj = 1, jpj
[921]	699	DO ji = 1, jpi
[2715]	700	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
	701	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
[921]	702	t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)
	703	ELSE
[2715]	704	ht_i(ji,jj,jl) = 0._wp
[5123]	705	t_su(ji,jj,jl) = rt0
[921]	706	ENDIF
[5123]	707	END DO
	708	END DO
	709	END DO
[2715]	710	!
[3294]	711	CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn )
	712	CALL wrk_dealloc( jpi,jpj, zworka )
	713	CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	714	!
[921]	715	END SUBROUTINE lim_itd_shiftice
[2715]	716
[825]	717
[4869]	718	SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
[921]	719	!!------------------------------------------------------------------
	720	!! * ROUTINE lim_itd_th_reb *
[2715]	721	!!
[921]	722	!! ** Purpose : rebin - rebins thicknesses into defined categories
	723	!!
	724	!! ** Method :
	725	!!------------------------------------------------------------------
[2715]	726	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	727	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	728	!
	729	INTEGER :: ji,jj, jl ! dummy loop indices
	730	INTEGER :: zshiftflag ! = .true. if ice must be shifted
[921]	731	CHARACTER (len = 15) :: fieldid
[825]	732
[3294]	733	INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index
	734	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! ice area and volume transferred
[825]	735
[3294]	736	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
	737	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
[2715]	738	!!------------------------------------------------------------------
[4688]	739	!! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate
[3294]	740
	741	CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger
	742	CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )
	743	CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
[2715]	744	!
	745	IF( con_i ) THEN ! conservation check
[834]	746	CALL lim_column_sum (jpl, v_i, vt_i_init)
	747	CALL lim_column_sum (jpl, v_s, vt_s_init)
	748	ENDIF
[825]	749
[921]	750	!
	751	!------------------------------------------------------------------------------
	752	! 1) Compute ice thickness.
	753	!------------------------------------------------------------------------------
[825]	754	DO jl = klbnd, kubnd
	755	DO jj = 1, jpj
[921]	756	DO ji = 1, jpi
[5134]	757	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
	758	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
[2715]	759	END DO
	760	END DO
	761	END DO
[825]	762
[921]	763	!------------------------------------------------------------------------------
[5134]	764	! 2) If a category thickness is not in bounds, shift the
[921]	765	! entire area, volume, and energy to the neighboring category
	766	!------------------------------------------------------------------------------
[825]	767	!-------------------------
	768	! Initialize shift arrays
	769	!-------------------------
	770	DO jl = klbnd, kubnd
[2715]	771	zdonor(:,:,jl) = 0
	772	zdaice(:,:,jl) = 0._wp
	773	zdvice(:,:,jl) = 0._wp
[825]	774	END DO
	775
	776	!-------------------------
	777	! Move thin categories up
	778	!-------------------------
	779
	780	DO jl = klbnd, kubnd - 1 ! loop over category boundaries
	781
[921]	782	!---------------------------------------
	783	! identify thicknesses that are too big
	784	!---------------------------------------
[869]	785	zshiftflag = 0
[825]	786
	787	DO jj = 1, jpj
	788	DO ji = 1, jpi
[2715]	789	IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN
[869]	790	zshiftflag = 1
[825]	791	zdonor(ji,jj,jl) = jl
[4161]	792	! begin TECLIM change
[4293]	793	!zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp
	794	!zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)(hi_max(jl)+hi_max(jl-1)) 0.5_wp
[4161]	795	! end TECLIM change
[4293]	796	! clem: how much of a_i you send in cat sup is somewhat arbitrary
[5134]	797	zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)
	798	zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
[825]	799	ENDIF
[5123]	800	END DO
	801	END DO
[2715]	802	IF(lk_mpp) CALL mpp_max( zshiftflag )
[825]	803
[2715]	804	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	805	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	806	! Reset shift parameters
[2715]	807	zdonor(:,:,jl) = 0
	808	zdaice(:,:,jl) = 0._wp
	809	zdvice(:,:,jl) = 0._wp
	810	ENDIF
	811	!
[5134]	812	END DO
[825]	813
	814	!----------------------------
	815	! Move thick categories down
	816	!----------------------------
	817
	818	DO jl = kubnd - 1, 1, -1 ! loop over category boundaries
	819
[921]	820	!-----------------------------------------
	821	! Identify thicknesses that are too small
	822	!-----------------------------------------
[869]	823	zshiftflag = 0
[825]	824
[4688]	825	DO jj = 1, jpj
	826	DO ji = 1, jpi
	827	IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
	828	!
	829	zshiftflag = 1
	830	zdonor(ji,jj,jl) = jl + 1
	831	zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)
	832	zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
	833	ENDIF
[5123]	834	END DO
	835	END DO
[4688]	836
	837	IF(lk_mpp) CALL mpp_max( zshiftflag )
	838
	839	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	840	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
	841	! Reset shift parameters
	842	zdonor(:,:,jl) = 0
	843	zdaice(:,:,jl) = 0._wp
	844	zdvice(:,:,jl) = 0._wp
	845	ENDIF
	846
	847	! ! clem-change begin: why not doing that?
[4333]	848	! DO jj = 1, jpj
	849	! DO ji = 1, jpi
[4688]	850	! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
	851	! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10
	852	! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)
[4333]	853	! ENDIF
[5123]	854	! END DO
	855	! END DO
[4333]	856	! clem-change end
[825]	857
[5123]	858	END DO
[825]	859
[921]	860	!------------------------------------------------------------------------------
[5134]	861	! 3) Conservation check
[921]	862	!------------------------------------------------------------------------------
[825]	863
[2715]	864	IF( con_i ) THEN
[921]	865	CALL lim_column_sum (jpl, v_i, vt_i_final)
	866	fieldid = ' v_i : limitd_reb '
	867	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
[825]	868
[921]	869	CALL lim_column_sum (jpl, v_s, vt_s_final)
	870	fieldid = ' v_s : limitd_reb '
	871	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	872	ENDIF
[2715]	873	!
[3294]	874	CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger
	875	CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )
	876	CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
	877
[921]	878	END SUBROUTINE lim_itd_th_reb
[825]	879
	880	#else
[2715]	881	!!----------------------------------------------------------------------
	882	!! Default option Dummy module NO LIM sea-ice model
	883	!!----------------------------------------------------------------------
[825]	884	CONTAINS
	885	SUBROUTINE lim_itd_th_rem
	886	END SUBROUTINE lim_itd_th_rem
	887	SUBROUTINE lim_itd_fitline
	888	END SUBROUTINE lim_itd_fitline
	889	SUBROUTINE lim_itd_shiftice
	890	END SUBROUTINE lim_itd_shiftice
	891	SUBROUTINE lim_itd_th_reb
	892	END SUBROUTINE lim_itd_th_reb
	893	#endif
[2715]	894	!!======================================================================
[921]	895	END MODULE limitd_th

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

Download in other formats: