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

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 8355

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

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

Download in other formats: