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

Context Navigation

source: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 9199

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

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

Download in other formats: