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 @ 4332

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

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

Download in other formats: