Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

limitd_th.F90 @ 7910

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

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 7910

Download in other formats: