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

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

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

Download in other formats: