Context Navigation

limitd_th.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

Property svn:keywords set to Id

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

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

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

Context Navigation

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

Download in other formats: