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

Context Navigation

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

Last change on this file since 3577 was 3558, checked in by rblod, 12 years ago
Fix issues when using key_nosignedzeo, see ticket #996
Property svn:keywords set to `Id`
File size: 46.3 KB

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

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

Download in other formats: