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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5183

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

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

Download in other formats: