Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

limitd_th.F90 @ 5418

Last change on this file since 5418 was 5418, checked in by deazer, 9 years ago
Removed SVN KEYWORDS ready for adding code changes before fcm merges
File size: 38.4 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
[4869]	8	!! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age
[2715]	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_rem :
	16	!! lim_itd_th_reb :
	17	!! lim_itd_fitline :
	18	!! lim_itd_shiftice :
[2528]	19	!!----------------------------------------------------------------------
[4161]	20	USE dom_ice ! LIM-3 domain
	21	USE par_oce ! ocean parameters
	22	USE dom_oce ! ocean domain
	23	USE phycst ! physical constants (ocean directory)
	24	USE thd_ice ! LIM-3 thermodynamic variables
	25	USE ice ! LIM-3 variables
	26	USE limvar ! LIM-3 variables
	27	USE prtctl ! Print control
	28	USE in_out_manager ! I/O manager
	29	USE lib_mpp ! MPP library
	30	USE wrk_nemo ! work arrays
	31	USE lib_fortran ! to use key_nosignedzero
[5123]	32	USE limcons ! conservation tests
[921]	33
[825]	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC lim_itd_th_rem
	38	PUBLIC lim_itd_th_reb
[825]	39
	40	!!----------------------------------------------------------------------
[4161]	41	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]	42	!! $Id$
[2715]	43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	44	!!----------------------------------------------------------------------
	45	CONTAINS
	46
[4869]	47	SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
[921]	48	!!------------------------------------------------------------------
	49	!! * ROUTINE lim_itd_th_rem *
	50	!!
[2715]	51	!! ** Purpose : computes the redistribution of ice thickness
	52	!! after thermodynamic growth of ice thickness
	53	!!
[921]	54	!! ** Method : Linear remapping
	55	!!
[2715]	56	!! References : W.H. Lipscomb, JGR 2001
[921]	57	!!------------------------------------------------------------------
[2715]	58	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	59	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	60	INTEGER , INTENT (in) :: kt ! Ocean time step
	61	!
	62	INTEGER :: ji, jj, jl ! dummy loop index
[4161]	63	INTEGER :: ii, ij ! 2D corresponding indices to ji
	64	INTEGER :: nd ! local integer
[2715]	65	REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars
[4161]	66	REAL(wp) :: zx2, zwk2, zda0, zetamax ! - -
[5123]	67	REAL(wp) :: zx3
[2715]	68	CHARACTER (len = 15) :: fieldid
[825]	69
[3294]	70	INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index
[825]	71
[3294]	72	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdhice ! ice thickness increment
	73	REAL(wp), POINTER, DIMENSION(:,:,:) :: g0 ! coefficients for fitting the line of the ITD
	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: g1 ! coefficients for fitting the line of the ITD
	75	REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness
	76	REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness
[4872]	77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness
[3294]	78	REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es
	79	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume
	80	REAL(wp), POINTER, DIMENSION(:) :: zvetamin, zvetamax ! maximum values for etas
	81	INTEGER , POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions
	82	INTEGER :: nbrem ! number of cells with ice to transfer
	83	REAL(wp) :: zslope ! used to compute local thermodynamic "speeds"
	84	REAL(wp), POINTER, DIMENSION(:,:) :: zhb0, zhb1 ! category boundaries for thinnes categories
	85	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
	86	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
	87	REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories
	88	REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init, et_s_final ! snow energy summed over categories
	89	INTEGER , POINTER, DIMENSION(:,:) :: zremap_flag ! compute remapping or not ????
	90	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhbnew ! new boundaries of ice categories
	91	!!------------------------------------------------------------------
[825]	92
[3294]	93	CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer
	94	CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer
[4872]	95	CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
[3294]	96	CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )
	97	CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )
	98	CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )
	99	CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	100	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]	101
[921]	102	!!----------------------------------------------------------------------------------------------
	103	!! 0) Conservation checkand changes in each ice category
	104	!!----------------------------------------------------------------------------------------------
[2715]	105	IF( con_i ) THEN
[825]	106	CALL lim_column_sum (jpl, v_i, vt_i_init)
	107	CALL lim_column_sum (jpl, v_s, vt_s_init)
	108	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
	109	dummy_es(:,:,:) = e_s(:,:,1,:)
	110	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
	111	ENDIF
	112
[921]	113	!!----------------------------------------------------------------------------------------------
	114	!! 1) Compute thickness and changes in each ice category
	115	!!----------------------------------------------------------------------------------------------
[2715]	116	IF( kt == nit000 .AND. lwp) THEN
[921]	117	WRITE(numout,*)
	118	WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'
	119	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	120	WRITE(numout,*) ' klbnd : ', klbnd
	121	WRITE(numout,*) ' kubnd : ', kubnd
	122	ENDIF
[825]	123
[2715]	124	zdhice(:,:,:) = 0._wp
[921]	125	DO jl = klbnd, kubnd
	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
[5134]	128	rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes
[4990]	129	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
[5134]	130	rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) !0 if no ice and 1 if yes
[4990]	131	zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
[5134]	132	zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)
[921]	133	END DO
	134	END DO
	135	END DO
	136
	137	!-----------------------------------------------------------------------------------------------
	138	! 2) Compute fractional ice area in each grid cell
	139	!-----------------------------------------------------------------------------------------------
[2715]	140	at_i(:,:) = 0._wp
[825]	141	DO jl = klbnd, kubnd
[2715]	142	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[825]	143	END DO
	144
[921]	145	!-----------------------------------------------------------------------------------------------
	146	! 3) Identify grid cells with ice
	147	!-----------------------------------------------------------------------------------------------
[825]	148	nbrem = 0
	149	DO jj = 1, jpj
	150	DO ji = 1, jpi
[5123]	151	IF ( at_i(ji,jj) > epsi10 ) THEN
[825]	152	nbrem = nbrem + 1
	153	nind_i(nbrem) = ji
	154	nind_j(nbrem) = jj
[3294]	155	zremap_flag(ji,jj) = 1
[825]	156	ELSE
[3294]	157	zremap_flag(ji,jj) = 0
[825]	158	ENDIF
[5123]	159	END DO
	160	END DO
[825]	161
[921]	162	!-----------------------------------------------------------------------------------------------
	163	! 4) Compute new category boundaries
	164	!-----------------------------------------------------------------------------------------------
[825]	165	!- 4.1 Compute category boundaries
[2715]	166	zhbnew(:,:,:) = 0._wp
[825]	167
	168	DO jl = klbnd, kubnd - 1
	169	DO ji = 1, nbrem
[4161]	170	ii = nind_i(ji)
	171	ij = nind_j(ji)
[825]	172	!
[4688]	173	zhbnew(ii,ij,jl) = hi_max(jl)
[4872]	174	IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
[825]	175	!interpolate between adjacent category growth rates
[4872]	176	zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
	177	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
	178	ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN
[4161]	179	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
[4872]	180	ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN
[4161]	181	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
[825]	182	ENDIF
[2715]	183	END DO
[825]	184
[921]	185	!- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
[825]	186	DO ji = 1, nbrem
[4161]	187	ii = nind_i(ji)
	188	ij = nind_j(ji)
[4688]	189	IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN
[4161]	190	zremap_flag(ii,ij) = 0
[4688]	191	ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN
[4161]	192	zremap_flag(ii,ij) = 0
[825]	193	ENDIF
	194
[921]	195	!- 4.3 Check that each zhbnew does not exceed maximal values hi_max
[4688]	196	IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
	197	IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
	198	END DO
	199
[5123]	200	END DO
[825]	201
[921]	202	!-----------------------------------------------------------------------------------------------
	203	! 5) Identify cells where ITD is to be remapped
	204	!-----------------------------------------------------------------------------------------------
	205	nbrem = 0
	206	DO jj = 1, jpj
	207	DO ji = 1, jpi
[4688]	208	IF( zremap_flag(ji,jj) == 1 ) THEN
[921]	209	nbrem = nbrem + 1
	210	nind_i(nbrem) = ji
	211	nind_j(nbrem) = jj
	212	ENDIF
[4688]	213	END DO
	214	END DO
[825]	215
[921]	216	!-----------------------------------------------------------------------------------------------
	217	! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
	218	!-----------------------------------------------------------------------------------------------
	219	DO jj = 1, jpj
	220	DO ji = 1, jpi
[4869]	221	zhb0(ji,jj) = hi_max(0) ! 0eme
	222	zhb1(ji,jj) = hi_max(1) ! 1er
[825]	223
[2715]	224	zhbnew(ji,jj,klbnd-1) = 0._wp
[825]	225
[2715]	226	IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
[5134]	227	zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
[921]	228	ELSE
[4293]	229	zhbnew(ji,jj,kubnd) = hi_max(kubnd)
	230	!!? clem bug: since hi_max(jpl)=99, this limit is very high
	231	!!? but I think it is erased in fitline subroutine
[921]	232	ENDIF
[825]	233
[5134]	234	END DO
	235	END DO
[825]	236
[921]	237	!-----------------------------------------------------------------------------------------------
	238	! 7) Compute g(h)
	239	!-----------------------------------------------------------------------------------------------
	240	!- 7.1 g(h) for category 1 at start of time step
[5123]	241	CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
[2715]	242	& hR(:,:,klbnd), zremap_flag )
[825]	243
[921]	244	!- 7.2 Area lost due to melting of thin ice (first category, klbnd)
	245	DO ji = 1, nbrem
[4161]	246	ii = nind_i(ji)
	247	ij = nind_j(ji)
[825]	248
[5123]	249	IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
[4161]	250	zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
[5123]	251	IF( zdh0 < 0.0 ) THEN !remove area from category 1
	252	zdh0 = MIN( -zdh0, hi_max(klbnd) )
[825]	253
[921]	254	!Integrate g(1) from 0 to dh0 to estimate area melted
[5123]	255	zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
	256	IF( zetamax > 0.0 ) THEN
[921]	257	zx1 = zetamax
[5123]	258	zx2 = 0.5 * zetamax * zetamax
[4161]	259	zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed
[921]	260	! Constrain new thickness <= ht_i
[5123]	261	zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0
[921]	262	!ice area lost due to melting of thin ice
[5123]	263	zda0 = MIN( zda0, zdamax )
[825]	264
[921]	265	! Remove area, conserving volume
[5123]	266	ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
[4161]	267	a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0
[5123]	268	v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
	269	ENDIF
[825]	270
[5134]	271	ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0
[5123]	272	zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )
[921]	273	! zhbnew was 0, and is shifted to the right to account for thin ice
	274	! growth in openwater (F0 = f1)
	275	ENDIF ! zdh0
[825]	276
[2715]	277	ENDIF ! a_i > epsi10
[825]	278
[5134]	279	END DO
[825]	280
[921]	281	!- 7.3 g(h) for each thickness category
	282	DO jl = klbnd, kubnd
[5123]	283	CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
	284	& g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
[921]	285	END DO
[825]	286
[921]	287	!-----------------------------------------------------------------------------------------------
	288	! 8) Compute area and volume to be shifted across each boundary
	289	!-----------------------------------------------------------------------------------------------
[825]	290
[921]	291	DO jl = klbnd, kubnd - 1
	292	DO jj = 1, jpj
	293	DO ji = 1, jpi
	294	zdonor(ji,jj,jl) = 0
	295	zdaice(ji,jj,jl) = 0.0
	296	zdvice(ji,jj,jl) = 0.0
	297	END DO
	298	END DO
[825]	299
[921]	300	DO ji = 1, nbrem
[4161]	301	ii = nind_i(ji)
	302	ij = nind_j(ji)
[825]	303
[5123]	304	IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
[825]	305
[921]	306	! left and right integration limits in eta space
[5123]	307	zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
	308	zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
[4161]	309	zdonor(ii,ij,jl) = jl
[825]	310
[921]	311	ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
[825]	312
[921]	313	! left and right integration limits in eta space
	314	zvetamin(ji) = 0.0
[5123]	315	zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
[4161]	316	zdonor(ii,ij,jl) = jl + 1
[825]	317
[921]	318	ENDIF ! zhbnew(jl) > hi_max(jl)
[825]	319
[5123]	320	zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
[921]	321	zetamin = zvetamin(ji)
[825]	322
[921]	323	zx1 = zetamax - zetamin
[5123]	324	zwk1 = zetamin * zetamin
	325	zwk2 = zetamax * zetamax
	326	zx2 = 0.5 * ( zwk2 - zwk1 )
[921]	327	zwk1 = zwk1 * zetamin
	328	zwk2 = zwk2 * zetamax
[5123]	329	zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 )
[4161]	330	nd = zdonor(ii,ij,jl)
	331	zdaice(ii,ij,jl) = g1(ii,ij,nd)zx2 + g0(ii,ij,nd)zx1
[4688]	332	zdvice(ii,ij,jl) = g1(ii,ij,nd)zx3 + g0(ii,ij,nd)zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
[921]	333
[5123]	334	END DO
[921]	335	END DO ! jl klbnd -> kubnd - 1
	336
	337	!!----------------------------------------------------------------------------------------------
	338	!! 9) Shift ice between categories
	339	!!----------------------------------------------------------------------------------------------
	340	CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
	341
	342	!!----------------------------------------------------------------------------------------------
	343	!! 10) Make sure ht_i >= minimum ice thickness hi_min
	344	!!----------------------------------------------------------------------------------------------
	345
	346	DO ji = 1, nbrem
[4161]	347	ii = nind_i(ji)
	348	ij = nind_j(ji)
[5123]	349	IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
[5134]	350	a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
[5123]	351	ht_i(ii,ij,1) = rn_himin
[921]	352	ENDIF
[5123]	353	END DO
[921]	354
	355	!!----------------------------------------------------------------------------------------------
	356	!! 11) Conservation check
	357	!!----------------------------------------------------------------------------------------------
[825]	358	IF ( con_i ) THEN
	359	CALL lim_column_sum (jpl, v_i, vt_i_final)
	360	fieldid = ' v_i : limitd_th '
	361	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
	362
	363	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_final)
	364	fieldid = ' e_i : limitd_th '
	365	CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)
	366
	367	CALL lim_column_sum (jpl, v_s, vt_s_final)
	368	fieldid = ' v_s : limitd_th '
	369	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	370
	371	dummy_es(:,:,:) = e_s(:,:,1,:)
	372	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
	373	fieldid = ' e_s : limitd_th '
	374	CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
	375	ENDIF
	376
[3294]	377	CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer
	378	CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer
[4872]	379	CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
[3294]	380	CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )
	381	CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )
	382	CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )
	383	CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	384	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 )
	385
[921]	386	END SUBROUTINE lim_itd_th_rem
[825]	387
	388
[5123]	389	SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
[921]	390	!!------------------------------------------------------------------
	391	!! * ROUTINE lim_itd_fitline *
	392	!!
[2715]	393	!! ** Purpose : fit g(h) with a line using area, volume constraints
[921]	394	!!
[2715]	395	!! ** Method : Fit g(h) with a line, satisfying area and volume constraints.
	396	!! To reduce roundoff errors caused by large values of g0 and g1,
	397	!! we actually compute g(eta), where eta = h - hL, and hL is the
	398	!! left boundary.
[921]	399	!!------------------------------------------------------------------
[2715]	400	INTEGER , INTENT(in ) :: num_cat ! category index
	401	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: HbL, HbR ! left and right category boundaries
	402	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hice ! ice thickness
	403	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: g0, g1 ! coefficients in linear equation for g(eta)
	404	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hL ! min value of range over which g(h) > 0
	405	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hR ! max value of range over which g(h) > 0
[3294]	406	INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag !
[2715]	407	!
	408	INTEGER :: ji,jj ! horizontal indices
	409	REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL)
	410	REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL)
	411	REAL(wp) :: zdhr ! 1 / (hR - hL)
	412	REAL(wp) :: zwk1, zwk2 ! temporary variables
	413	!!------------------------------------------------------------------
	414	!
	415	!
[825]	416	DO jj = 1, jpj
	417	DO ji = 1, jpi
[2715]	418	!
[4333]	419	IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 &
[3294]	420	& .AND. hice(ji,jj) > 0._wp ) THEN
[825]	421
[921]	422	! Initialize hL and hR
	423
[825]	424	hL(ji,jj) = HbL(ji,jj)
	425	hR(ji,jj) = HbR(ji,jj)
	426
[921]	427	! Change hL or hR if hice falls outside central third of range
[825]	428
[5123]	429	zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
	430	zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
[825]	431
[2715]	432	IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
	433	ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
[825]	434	ENDIF
	435
[921]	436	! Compute coefficients of g(eta) = g0 + g1*eta
	437
[2715]	438	zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
	439	zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
	440	zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
[5123]	441	g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
	442	g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
[2715]	443	!
	444	ELSE ! remap_flag = .false. or a_i < epsi10
	445	hL(ji,jj) = 0._wp
	446	hR(ji,jj) = 0._wp
	447	g0(ji,jj) = 0._wp
	448	g1(ji,jj) = 0._wp
	449	ENDIF ! a_i > epsi10
	450	!
	451	END DO
	452	END DO
	453	!
	454	END SUBROUTINE lim_itd_fitline
[825]	455
	456
[2715]	457	SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	458	!!------------------------------------------------------------------
	459	!! * ROUTINE lim_itd_shiftice *
[2715]	460	!!
	461	!! ** Purpose : shift ice across category boundaries, conserving everything
[921]	462	!! ( area, volume, energy, age*vol, and mass of salt )
	463	!!
	464	!! ** Method :
	465	!!------------------------------------------------------------------
[3294]	466	INTEGER , INTENT(in ) :: klbnd ! Start thickness category index point
	467	INTEGER , INTENT(in ) :: kubnd ! End point on which the the computation is applied
[2715]	468	INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in ) :: zdonor ! donor category index
	469	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdaice ! ice area transferred across boundary
	470	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdvice ! ice volume transferred across boundary
[825]	471
[2715]	472	INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices
[4161]	473	INTEGER :: ii, ij ! indices when changing from 2D-1D is done
[825]	474
[3294]	475	REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn
	476	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
[825]	477
[2715]	478	REAL(wp) :: zdvsnow, zdesnow ! snow volume and energy transferred
	479	REAL(wp) :: zdeice ! ice energy transferred
	480	REAL(wp) :: zdsm_vice ! ice salinity times volume transferred
	481	REAL(wp) :: zdo_aice ! ice age times volume transferred
	482	REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred
[825]	483
[3294]	484	INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions
[825]	485
[2715]	486	INTEGER :: nbrem ! number of cells with ice to transfer
[825]	487
[2715]	488	LOGICAL :: zdaice_negative ! true if daice < -puny
	489	LOGICAL :: zdvice_negative ! true if dvice < -puny
	490	LOGICAL :: zdaice_greater_aicen ! true if daice > aicen
	491	LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen
	492	!!------------------------------------------------------------------
[825]	493
[3294]	494	CALL wrk_alloc( jpi,jpj,jpl, zaTsfn )
	495	CALL wrk_alloc( jpi,jpj, zworka )
	496	CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	497
[921]	498	!----------------------------------------------------------------------------------------------
	499	! 1) Define a variable equal to a_i*T_su
	500	!----------------------------------------------------------------------------------------------
[825]	501
	502	DO jl = klbnd, kubnd
[5123]	503	zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
[2715]	504	END DO
[825]	505
[5134]	506	!clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10)
	507	! !----------------------------------------------------------------------------------------------
	508	! ! 2) Check for daice or dvice out of range, allowing for roundoff error
	509	! !----------------------------------------------------------------------------------------------
	510	! ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
	511	! ! has a small area, with h(n) very close to a boundary. Then
	512	! ! the coefficients of g(h) are large, and the computed daice and
	513	! ! dvice can be in error. If this happens, it is best to transfer
	514	! ! either the entire category or nothing at all, depending on which
	515	! ! side of the boundary hice(n) lies.
	516	! !-----------------------------------------------------------------
	517	! DO jl = klbnd, kubnd-1
	518	!
	519	! zdaice_negative = .false.
	520	! zdvice_negative = .false.
	521	! zdaice_greater_aicen = .false.
	522	! zdvice_greater_vicen = .false.
	523	!
	524	! DO jj = 1, jpj
	525	! DO ji = 1, jpi
	526	!
	527	! IF (zdonor(ji,jj,jl) > 0) THEN
	528	! jl1 = zdonor(ji,jj,jl)
	529	!
	530	! IF (zdaice(ji,jj,jl) < 0.0) THEN
	531	! IF (zdaice(ji,jj,jl) > -epsi10) THEN
	532	! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &
	533	! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
	534	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
	535	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	536	! ELSE
	537	! zdaice(ji,jj,jl) = 0.0 ! shift no ice
	538	! zdvice(ji,jj,jl) = 0.0
	539	! ENDIF
	540	! ELSE
	541	! zdaice_negative = .true.
	542	! ENDIF
	543	! ENDIF
	544	!
	545	! IF (zdvice(ji,jj,jl) < 0.0) THEN
	546	! IF (zdvice(ji,jj,jl) > -epsi10 ) THEN
	547	! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &
	548	! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
	549	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
	550	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	551	! ELSE
	552	! zdaice(ji,jj,jl) = 0.0 ! shift no ice
	553	! zdvice(ji,jj,jl) = 0.0
	554	! ENDIF
	555	! ELSE
	556	! zdvice_negative = .true.
	557	! ENDIF
	558	! ENDIF
	559	!
	560	! ! If daice is close to aicen, set daice = aicen.
	561	! IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN
	562	! IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN
	563	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
	564	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	565	! ELSE
	566	! zdaice_greater_aicen = .true.
	567	! ENDIF
	568	! ENDIF
	569	!
	570	! IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN
	571	! IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN
	572	! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
	573	! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
	574	! ELSE
	575	! zdvice_greater_vicen = .true.
	576	! ENDIF
	577	! ENDIF
	578	!
	579	! ENDIF ! donor > 0
	580	! END DO
	581	! END DO
	582	!
	583	! END DO
	584	!clem
[921]	585	!-------------------------------------------------------------------------------
	586	! 3) Transfer volume and energy between categories
	587	!-------------------------------------------------------------------------------
[825]	588
	589	DO jl = klbnd, kubnd - 1
	590	nbrem = 0
	591	DO jj = 1, jpj
	592	DO ji = 1, jpi
[5123]	593	IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
[825]	594	nbrem = nbrem + 1
	595	nind_i(nbrem) = ji
	596	nind_j(nbrem) = jj
[5123]	597	ENDIF
[825]	598	END DO
	599	END DO
	600
	601	DO ji = 1, nbrem
[4161]	602	ii = nind_i(ji)
	603	ij = nind_j(ji)
[825]	604
[4161]	605	jl1 = zdonor(ii,ij,jl)
[5134]	606	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) )
	607	zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch
[2715]	608	IF( jl1 == jl) THEN ; jl2 = jl1+1
[5123]	609	ELSE ; jl2 = jl
[825]	610	ENDIF
	611
	612	!--------------
	613	! Ice areas
	614	!--------------
	615
[4161]	616	a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
	617	a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
[825]	618
	619	!--------------
	620	! Ice volumes
	621	!--------------
	622
[4161]	623	v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)
	624	v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
[825]	625
	626	!--------------
	627	! Snow volumes
	628	!--------------
	629
[4688]	630	zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij)
[4161]	631	v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
	632	v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow
[825]	633
	634	!--------------------
	635	! Snow heat content
	636	!--------------------
	637
[4688]	638	zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij)
[4161]	639	e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow
	640	e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow
[825]	641
	642	!--------------
	643	! Ice age
	644	!--------------
	645
[4688]	646	zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	647	oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice
	648	oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice
[825]	649
	650	!--------------
	651	! Ice salinity
	652	!--------------
	653
[4688]	654	zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij)
[4161]	655	smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice
	656	smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice
[825]	657
	658	!---------------------
	659	! Surface temperature
	660	!---------------------
	661
[4688]	662	zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	663	zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf
	664	zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf
[825]	665
[5123]	666	END DO
[825]	667
	668	!------------------
	669	! Ice heat content
	670	!------------------
	671
	672	DO jk = 1, nlay_i
	673	DO ji = 1, nbrem
[4161]	674	ii = nind_i(ji)
	675	ij = nind_j(ji)
[825]	676
[4161]	677	jl1 = zdonor(ii,ij,jl)
[5123]	678	IF (jl1 == jl) THEN
[825]	679	jl2 = jl+1
	680	ELSE ! n1 = n+1
	681	jl2 = jl
	682	ENDIF
	683
[4161]	684	zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
	685	e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice
	686	e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice
[5123]	687	END DO
	688	END DO
[825]	689
	690	END DO ! boundaries, 1 to ncat-1
	691
	692	!-----------------------------------------------------------------
	693	! Update ice thickness and temperature
	694	!-----------------------------------------------------------------
	695
	696	DO jl = klbnd, kubnd
	697	DO jj = 1, jpj
[921]	698	DO ji = 1, jpi
[2715]	699	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
	700	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
[921]	701	t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)
	702	ELSE
[2715]	703	ht_i(ji,jj,jl) = 0._wp
[5123]	704	t_su(ji,jj,jl) = rt0
[921]	705	ENDIF
[5123]	706	END DO
	707	END DO
	708	END DO
[2715]	709	!
[3294]	710	CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn )
	711	CALL wrk_dealloc( jpi,jpj, zworka )
	712	CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer
	713	!
[921]	714	END SUBROUTINE lim_itd_shiftice
[2715]	715
[825]	716
[4869]	717	SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
[921]	718	!!------------------------------------------------------------------
	719	!! * ROUTINE lim_itd_th_reb *
[2715]	720	!!
[921]	721	!! ** Purpose : rebin - rebins thicknesses into defined categories
	722	!!
	723	!! ** Method :
	724	!!------------------------------------------------------------------
[2715]	725	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	726	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	727	!
	728	INTEGER :: ji,jj, jl ! dummy loop indices
	729	INTEGER :: zshiftflag ! = .true. if ice must be shifted
[921]	730	CHARACTER (len = 15) :: fieldid
[825]	731
[3294]	732	INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index
	733	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! ice area and volume transferred
[825]	734
[3294]	735	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
	736	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
[2715]	737	!!------------------------------------------------------------------
[3294]	738
	739	CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger
	740	CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )
	741	CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
[2715]	742	!
	743	IF( con_i ) THEN ! conservation check
[834]	744	CALL lim_column_sum (jpl, v_i, vt_i_init)
	745	CALL lim_column_sum (jpl, v_s, vt_s_init)
	746	ENDIF
[825]	747
[921]	748	!
	749	!------------------------------------------------------------------------------
	750	! 1) Compute ice thickness.
	751	!------------------------------------------------------------------------------
[825]	752	DO jl = klbnd, kubnd
	753	DO jj = 1, jpj
[921]	754	DO ji = 1, jpi
[5134]	755	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
	756	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
[2715]	757	END DO
	758	END DO
	759	END DO
[825]	760
[921]	761	!------------------------------------------------------------------------------
[5134]	762	! 2) If a category thickness is not in bounds, shift the
[921]	763	! entire area, volume, and energy to the neighboring category
	764	!------------------------------------------------------------------------------
[825]	765	!-------------------------
	766	! Initialize shift arrays
	767	!-------------------------
	768	DO jl = klbnd, kubnd
[2715]	769	zdonor(:,:,jl) = 0
	770	zdaice(:,:,jl) = 0._wp
	771	zdvice(:,:,jl) = 0._wp
[825]	772	END DO
	773
	774	!-------------------------
	775	! Move thin categories up
	776	!-------------------------
	777
	778	DO jl = klbnd, kubnd - 1 ! loop over category boundaries
	779
[921]	780	!---------------------------------------
	781	! identify thicknesses that are too big
	782	!---------------------------------------
[869]	783	zshiftflag = 0
[825]	784
	785	DO jj = 1, jpj
	786	DO ji = 1, jpi
[2715]	787	IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN
[869]	788	zshiftflag = 1
[825]	789	zdonor(ji,jj,jl) = jl
[4161]	790	! begin TECLIM change
[4293]	791	!zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp
	792	!zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)(hi_max(jl)+hi_max(jl-1)) 0.5_wp
[4161]	793	! end TECLIM change
[4293]	794	! clem: how much of a_i you send in cat sup is somewhat arbitrary
[5134]	795	zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)
	796	zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
[825]	797	ENDIF
[5123]	798	END DO
	799	END DO
[2715]	800	IF(lk_mpp) CALL mpp_max( zshiftflag )
[825]	801
[2715]	802	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	803	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	804	! Reset shift parameters
[2715]	805	zdonor(:,:,jl) = 0
	806	zdaice(:,:,jl) = 0._wp
	807	zdvice(:,:,jl) = 0._wp
	808	ENDIF
	809	!
[5134]	810	END DO
[825]	811
	812	!----------------------------
	813	! Move thick categories down
	814	!----------------------------
	815
	816	DO jl = kubnd - 1, 1, -1 ! loop over category boundaries
	817
[921]	818	!-----------------------------------------
	819	! Identify thicknesses that are too small
	820	!-----------------------------------------
[869]	821	zshiftflag = 0
[825]	822
[4688]	823	DO jj = 1, jpj
	824	DO ji = 1, jpi
	825	IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
	826	!
	827	zshiftflag = 1
	828	zdonor(ji,jj,jl) = jl + 1
	829	zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)
	830	zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
	831	ENDIF
[5123]	832	END DO
	833	END DO
[4688]	834
	835	IF(lk_mpp) CALL mpp_max( zshiftflag )
	836
	837	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	838	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
	839	! Reset shift parameters
	840	zdonor(:,:,jl) = 0
	841	zdaice(:,:,jl) = 0._wp
	842	zdvice(:,:,jl) = 0._wp
	843	ENDIF
	844
[5123]	845	END DO
[825]	846
[921]	847	!------------------------------------------------------------------------------
[5134]	848	! 3) Conservation check
[921]	849	!------------------------------------------------------------------------------
[825]	850
[2715]	851	IF( con_i ) THEN
[921]	852	CALL lim_column_sum (jpl, v_i, vt_i_final)
	853	fieldid = ' v_i : limitd_reb '
	854	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
[825]	855
[921]	856	CALL lim_column_sum (jpl, v_s, vt_s_final)
	857	fieldid = ' v_s : limitd_reb '
	858	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	859	ENDIF
[2715]	860	!
[3294]	861	CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger
	862	CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )
	863	CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
	864
[921]	865	END SUBROUTINE lim_itd_th_reb
[825]	866
	867	#else
[2715]	868	!!----------------------------------------------------------------------
	869	!! Default option Dummy module NO LIM sea-ice model
	870	!!----------------------------------------------------------------------
[825]	871	CONTAINS
	872	SUBROUTINE lim_itd_th_rem
	873	END SUBROUTINE lim_itd_th_rem
	874	SUBROUTINE lim_itd_fitline
	875	END SUBROUTINE lim_itd_fitline
	876	SUBROUTINE lim_itd_shiftice
	877	END SUBROUTINE lim_itd_shiftice
	878	SUBROUTINE lim_itd_th_reb
	879	END SUBROUTINE lim_itd_th_reb
	880	#endif
[2715]	881	!!======================================================================
[921]	882	END MODULE limitd_th

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

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

Context Navigation

source: branches/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5418

Download in other formats: