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

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

limitd_th.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4317

Last change on this file since 4317 was 4293, checked in by clem, 10 years ago
small corrections for ice categories hi_max and for shifting ice between categories
Property svn:keywords set to `Id`
File size: 49.1 KB

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

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

Download in other formats: