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

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 4634

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

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

Download in other formats: