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

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 3634

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

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

Download in other formats: