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

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

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

Download in other formats: