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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4293

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

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

Download in other formats: