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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limitd_th.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limitd_th.F90 @ 2011

Last change on this file since 2011 was 1465, checked in by smasson, 15 years ago
supress ice_oce module, see ticket:448
Property svn:keywords set to `Id`
File size: 47.8 KB

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

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

Download in other formats: