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.

icbutl.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 8306

Last change on this file since 8306 was 8306, checked in by clem, 7 years ago
step1: remove LIM2 from the code
Property svn:keywords set to `Id`
File size: 33.5 KB

Line
1	MODULE icbutl
2	!!======================================================================
3	!! * MODULE icbutl *
4	!! Icebergs: various iceberg utility routines
5	!!======================================================================
6	!! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code
7	!! - ! 2011-03 (Madec) Part conversion to NEMO form
8	!! - ! Removal of mapping from another grid
9	!! - ! 2011-04 (Alderson) Split into separate modules
10	!!----------------------------------------------------------------------
11	!!----------------------------------------------------------------------
12	!! icb_utl_interp :
13	!! icb_utl_bilin :
14	!! icb_utl_bilin_e :
15	!!----------------------------------------------------------------------
16	USE par_oce ! ocean parameters
17	USE dom_oce ! ocean domain
18	USE in_out_manager ! IO parameters
19	USE lbclnk ! lateral boundary condition
20	USE lib_mpp ! MPI code and lk_mpp in particular
21	USE icb_oce ! define iceberg arrays
22	USE sbc_oce ! ocean surface boundary conditions
23	#if defined key_lim3
24	USE ice, ONLY: u_ice, v_ice, htm_i ! LIM-3 variables
25	#endif
26
27	IMPLICIT NONE
28	PRIVATE
29
30	PUBLIC icb_utl_copy ! routine called in icbstp module
31	PUBLIC icb_utl_interp ! routine called in icbdyn, icbthm modules
32	PUBLIC icb_utl_bilin ! routine called in icbini, icbdyn modules
33	PUBLIC icb_utl_bilin_x ! routine called in icbdyn module
34	PUBLIC icb_utl_add ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
35	PUBLIC icb_utl_delete ! routine called in icblbc, icbthm modules
36	PUBLIC icb_utl_destroy ! routine called in icbstp module
37	PUBLIC icb_utl_track ! routine not currently used, retain just in case
38	PUBLIC icb_utl_print_berg ! routine called in icbthm module
39	PUBLIC icb_utl_print ! routine called in icbini, icbstp module
40	PUBLIC icb_utl_count ! routine called in icbdia, icbini, icblbc, icbrst modules
41	PUBLIC icb_utl_incr ! routine called in icbini, icbclv modules
42	PUBLIC icb_utl_yearday ! routine called in icbclv, icbstp module
43	PUBLIC icb_utl_mass ! routine called in icbdia module
44	PUBLIC icb_utl_heat ! routine called in icbdia module
45
46	!!----------------------------------------------------------------------
47	!! NEMO/OPA 3.3 , NEMO Consortium (2011)
48	!! $Id$
49	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
50	!!-------------------------------------------------------------------------
51
52	CONTAINS
53
54	SUBROUTINE icb_utl_copy()
55	!!----------------------------------------------------------------------
56	!! * ROUTINE icb_utl_copy *
57	!!
58	!! ** Purpose : iceberg initialization.
59	!!
60	!! ** Method : - blah blah
61	!!----------------------------------------------------------------------
62
63	! copy nemo forcing arrays into iceberg versions with extra halo
64	! only necessary for variables not on T points
65	! and ssh which is used to calculate gradients
66
67	uo_e(:,:) = 0._wp ; uo_e(1:jpi, 1:jpj) = ssu_m(:,:) * umask(:,:,1)
68	vo_e(:,:) = 0._wp ; vo_e(1:jpi, 1:jpj) = ssv_m(:,:) * vmask(:,:,1)
69	ff_e(:,:) = 0._wp ; ff_e(1:jpi, 1:jpj) = ff_f (:,:)
70	tt_e(:,:) = 0._wp ; tt_e(1:jpi, 1:jpj) = sst_m(:,:)
71	fr_e(:,:) = 0._wp ; fr_e(1:jpi, 1:jpj) = fr_i (:,:)
72	ua_e(:,:) = 0._wp ; ua_e(1:jpi, 1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
73	va_e(:,:) = 0._wp ; va_e(1:jpi, 1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
74
75	CALL lbc_lnk_icb( uo_e, 'U', -1._wp, 1, 1 )
76	CALL lbc_lnk_icb( vo_e, 'V', -1._wp, 1, 1 )
77	CALL lbc_lnk_icb( ff_e, 'F', +1._wp, 1, 1 )
78	CALL lbc_lnk_icb( ua_e, 'U', -1._wp, 1, 1 )
79	CALL lbc_lnk_icb( va_e, 'V', -1._wp, 1, 1 )
80	CALL lbc_lnk_icb( fr_e, 'T', +1._wp, 1, 1 )
81	CALL lbc_lnk_icb( tt_e, 'T', +1._wp, 1, 1 )
82	#if defined key_lim3
83	hicth(:,:) = 0._wp ; hicth(1:jpi,1:jpj) = htm_i(:,:)
84	ui_e(:,:) = 0._wp ; ui_e(1:jpi, 1:jpj) = u_ice(:,:)
85	vi_e(:,:) = 0._wp ; vi_e(1:jpi, 1:jpj) = v_ice(:,:)
86	CALL lbc_lnk_icb(hicth, 'T', +1._wp, 1, 1 )
87	CALL lbc_lnk_icb( ui_e, 'U', -1._wp, 1, 1 )
88	CALL lbc_lnk_icb( vi_e, 'V', -1._wp, 1, 1 )
89	#endif
90
91	!! special for ssh which is used to calculate slope
92	!! so fudge some numbers all the way around the boundary
93
94	ssh_e(:,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
95	ssh_e(0 , :) = ssh_e(1 , :)
96	ssh_e(jpi+1, :) = ssh_e(jpi, :)
97	ssh_e(: , 0) = ssh_e(: , 1)
98	ssh_e(: ,jpj+1) = ssh_e(: ,jpj)
99	ssh_e(0,0) = ssh_e(1,1)
100	ssh_e(jpi+1,0) = ssh_e(jpi,1)
101	ssh_e(0,jpj+1) = ssh_e(1,jpj)
102	ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj)
103	CALL lbc_lnk_icb( ssh_e, 'T', +1._wp, 1, 1 )
104	!
105	END SUBROUTINE icb_utl_copy
106
107
108	SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i, &
109	& pj, pe2, pvo, pvi, pva, pssh_j, &
110	& psst, pcn, phi, pff )
111	!!----------------------------------------------------------------------
112	!! * ROUTINE icb_utl_interp *
113	!!
114	!! ** Purpose : interpolation
115	!!
116	!! ** Method : - interpolate from various ocean arrays onto iceberg position
117	!!
118	!! !!gm CAUTION here I do not care of the slip/no-slip conditions
119	!! this can be done later (not that easy to do...)
120	!! right now, U is 0 in land so that the coastal value of velocity parallel to the coast
121	!! is half the off shore value, wile the normal-to-the-coast value is zero.
122	!! This is OK as a starting point.
123	!!
124	!!----------------------------------------------------------------------
125	REAL(wp), INTENT(in ) :: pi , pj ! position in (i,j) referential
126	REAL(wp), INTENT( out) :: pe1, pe2 ! i- and j scale factors
127	REAL(wp), INTENT( out) :: puo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds
128	REAL(wp), INTENT( out) :: pssh_i, pssh_j ! ssh i- & j-gradients
129	REAL(wp), INTENT( out) :: psst, pcn, phi, pff ! SST, ice concentration, ice thickness, Coriolis
130	!
131	REAL(wp) :: zcd, zmod ! local scalars
132	!!----------------------------------------------------------------------
133
134	pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors
135	pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
136	!
137	puo = icb_utl_bilin_h( uo_e, pi, pj, 'U' ) ! ocean velocities
138	pvo = icb_utl_bilin_h( vo_e, pi, pj, 'V' )
139	psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' ) ! SST
140	pcn = icb_utl_bilin_h( fr_e , pi, pj, 'T' ) ! ice concentration
141	pff = icb_utl_bilin_h( ff_e , pi, pj, 'F' ) ! Coriolis parameter
142	!
143	pua = icb_utl_bilin_h( ua_e , pi, pj, 'U' ) ! 10m wind
144	pva = icb_utl_bilin_h( va_e , pi, pj, 'V' ) ! here (ua,va) are stress => rough conversion from stress to speed
145	zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient
146	zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( puapua + pvapva) ) )
147	pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0
148	pva = pva * zmod
149
150	#if defined key_lim3
151	pui = icb_utl_bilin_h( ui_e, pi, pj, 'U' ) ! sea-ice velocities
152	pvi = icb_utl_bilin_h( vi_e, pi, pj, 'V' )
153	phi = icb_utl_bilin_h(hicth, pi, pj, 'T' ) ! ice thickness
154	#else
155	pui = 0._wp
156	pvi = 0._wp
157	phi = 0._wp
158	#endif
159
160	! Estimate SSH gradient in i- and j-direction (centred evaluation)
161	pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) - &
162	& icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 )
163	pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) - &
164	& icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 )
165	!
166	END SUBROUTINE icb_utl_interp
167
168
169	REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )
170	!!----------------------------------------------------------------------
171	!! * FUNCTION icb_utl_bilin *
172	!!
173	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
174	!! this version deals with extra halo points
175	!!
176	!! !!gm CAUTION an optional argument should be added to handle
177	!! the slip/no-slip conditions ==>>> to be done later
178	!!
179	!!----------------------------------------------------------------------
180	REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated
181	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
182	CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points
183	!
184	INTEGER :: ii, ij ! local integer
185	REAL(wp) :: zi, zj ! local real
186	!!----------------------------------------------------------------------
187	!
188	SELECT CASE ( cd_type )
189	CASE ( 'T' )
190	! note that here there is no +0.5 added
191	! since we're looking for four T points containing quadrant we're in of
192	! current T cell
193	ii = MAX(1, INT( pi ))
194	ij = MAX(1, INT( pj )) ! T-point
195	zi = pi - REAL(ii,wp)
196	zj = pj - REAL(ij,wp)
197	CASE ( 'U' )
198	ii = MAX(1, INT( pi-0.5 ))
199	ij = MAX(1, INT( pj )) ! U-point
200	zi = pi - 0.5 - REAL(ii,wp)
201	zj = pj - REAL(ij,wp)
202	CASE ( 'V' )
203	ii = MAX(1, INT( pi ))
204	ij = MAX(1, INT( pj-0.5 )) ! V-point
205	zi = pi - REAL(ii,wp)
206	zj = pj - 0.5 - REAL(ij,wp)
207	CASE ( 'F' )
208	ii = MAX(1, INT( pi-0.5 ))
209	ij = MAX(1, INT( pj-0.5 )) ! F-point
210	zi = pi - 0.5 - REAL(ii,wp)
211	zj = pj - 0.5 - REAL(ij,wp)
212	END SELECT
213	!
214	! find position in this processor. Prevent near edge problems (see #1389)
215
216	if (ii.lt.mig(1)) then
217	ii = 1
218	else if (ii.gt.mig(jpi)) then
219	ii = jpi
220	else
221	ii = mi1( ii )
222	end if
223
224	if (ij.lt.mjg(1)) then
225	ij = 1
226	else if (ij.gt.mjg(jpj)) then
227	ij = jpj
228	else
229	ij = mj1( ij )
230	end if
231
232	if (ij.eq.jpj) ij=ij-1
233	if (ii.eq.jpi) ii=ii-1
234
235	!
236	icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) &
237	& + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj
238	!
239	END FUNCTION icb_utl_bilin_h
240
241
242	REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
243	!!----------------------------------------------------------------------
244	!! * FUNCTION icb_utl_bilin *
245	!!
246	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
247	!!
248	!! !!gm CAUTION an optional argument should be added to handle
249	!! the slip/no-slip conditions ==>>> to be done later
250	!!
251	!!----------------------------------------------------------------------
252	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated
253	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
254	CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points
255	!
256	INTEGER :: ii, ij ! local integer
257	REAL(wp) :: zi, zj ! local real
258	!!----------------------------------------------------------------------
259	!
260	SELECT CASE ( cd_type )
261	CASE ( 'T' )
262	! note that here there is no +0.5 added
263	! since we're looking for four T points containing quadrant we're in of
264	! current T cell
265	ii = MAX(1, INT( pi ))
266	ij = MAX(1, INT( pj )) ! T-point
267	zi = pi - REAL(ii,wp)
268	zj = pj - REAL(ij,wp)
269	CASE ( 'U' )
270	ii = MAX(1, INT( pi-0.5 ))
271	ij = MAX(1, INT( pj )) ! U-point
272	zi = pi - 0.5 - REAL(ii,wp)
273	zj = pj - REAL(ij,wp)
274	CASE ( 'V' )
275	ii = MAX(1, INT( pi ))
276	ij = MAX(1, INT( pj-0.5 )) ! V-point
277	zi = pi - REAL(ii,wp)
278	zj = pj - 0.5 - REAL(ij,wp)
279	CASE ( 'F' )
280	ii = MAX(1, INT( pi-0.5 ))
281	ij = MAX(1, INT( pj-0.5 )) ! F-point
282	zi = pi - 0.5 - REAL(ii,wp)
283	zj = pj - 0.5 - REAL(ij,wp)
284	END SELECT
285	!
286	! find position in this processor. Prevent near edge problems (see #1389)
287
288	if (ii.lt.mig(1)) then
289	ii = 1
290	else if (ii.gt.mig(jpi)) then
291	ii = jpi
292	else
293	ii = mi1( ii )
294	end if
295
296	if (ij.lt.mjg(1)) then
297	ij = 1
298	else if (ij.gt.mjg(jpj)) then
299	ij = jpj
300	else
301	ij = mj1( ij )
302	end if
303
304	if (ij.eq.jpj) ij=ij-1
305	if (ii.eq.jpi) ii=ii-1
306
307	icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) &
308	& + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj
309	!
310	END FUNCTION icb_utl_bilin
311
312
313	REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj )
314	!!----------------------------------------------------------------------
315	!! * FUNCTION icb_utl_bilin_x *
316	!!
317	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
318	!! Special case for interpolating longitude
319	!!
320	!! !!gm CAUTION an optional argument should be added to handle
321	!! the slip/no-slip conditions ==>>> to be done later
322	!!
323	!!----------------------------------------------------------------------
324	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated
325	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
326	!
327	INTEGER :: ii, ij ! local integer
328	REAL(wp) :: zi, zj ! local real
329	REAL(wp) :: zret ! local real
330	REAL(wp), DIMENSION(4) :: z4
331	!!----------------------------------------------------------------------
332	!
333	! note that here there is no +0.5 added
334	! since we're looking for four T points containing quadrant we're in of
335	! current T cell
336	ii = MAX(1, INT( pi ))
337	ij = MAX(1, INT( pj )) ! T-point
338	zi = pi - REAL(ii,wp)
339	zj = pj - REAL(ij,wp)
340	!
341	! find position in this processor. Prevent near edge problems (see #1389)
342
343	if (ii.lt.mig(1)) then
344	ii = 1
345	else if (ii.gt.mig(jpi)) then
346	ii = jpi
347	else
348	ii = mi1( ii )
349	end if
350
351	if (ij.lt.mjg(1)) then
352	ij = 1
353	else if (ij.gt.mjg(jpj)) then
354	ij = jpj
355	else
356	ij = mj1( ij )
357	end if
358
359	if (ij.eq.jpj) ij=ij-1
360	if (ii.eq.jpi) ii=ii-1
361
362	z4(1) = pfld(ii ,ij )
363	z4(2) = pfld(ii+1,ij )
364	z4(3) = pfld(ii ,ij+1)
365	z4(4) = pfld(ii+1,ij+1)
366	IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
367	WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
368	ENDIF
369	!
370	zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
371	IF( zret > 180._wp ) zret = zret - 360._wp
372	icb_utl_bilin_x = zret
373	!
374	END FUNCTION icb_utl_bilin_x
375
376
377	REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
378	!!----------------------------------------------------------------------
379	!! * FUNCTION dom_init *
380	!!
381	!! ** Purpose : bilinear interpolation at berg location of horizontal scale factor
382	!! ** Method : interpolation done using the 4 nearest grid points among
383	!! t-, u-, v-, and f-points.
384	!!----------------------------------------------------------------------
385	REAL(wp), DIMENSION(:,:), INTENT(in) :: pet, peu, pev, pef ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
386	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
387	!
388	INTEGER :: ii, ij, icase ! local integer
389	!
390	! weights corresponding to corner points of a T cell quadrant
391	REAL(wp) :: zi, zj ! local real
392	!
393	! values at corner points of a T cell quadrant
394	! 00 = bottom left, 10 = bottom right, 01 = top left, 11 = top right
395	REAL(wp) :: ze00, ze10, ze01, ze11
396	!!----------------------------------------------------------------------
397	!
398	ii = MAX(1, INT( pi )) ; ij = MAX(1, INT( pj )) ! left bottom T-point (i,j) indices
399
400	! fractional box spacing
401	! 0 <= zi < 0.5, 0 <= zj < 0.5 --> NW quadrant of current T cell
402	! 0.5 <= zi < 1 , 0 <= zj < 0.5 --> NE quadrant
403	! 0 <= zi < 0.5, 0.5 <= zj < 1 --> SE quadrant
404	! 0.5 <= zi < 1 , 0.5 <= zj < 1 --> SW quadrant
405
406	zi = pi - REAL(ii,wp) !!gm use here mig, mjg arrays
407	zj = pj - REAL(ij,wp)
408
409	! find position in this processor. Prevent near edge problems (see #1389)
410
411	if (ii.lt.mig(1)) then
412	ii = 1
413	else if (ii.gt.mig(jpi)) then
414	ii = jpi
415	else
416	ii = mi1( ii )
417	end if
418
419	if (ij.lt.mjg(1)) then
420	ij = 1
421	else if (ij.gt.mjg(jpj)) then
422	ij = jpj
423	else
424	ij = mj1( ij )
425	end if
426
427	if (ij.eq.jpj) ij=ij-1
428	if (ii.eq.jpi) ii=ii-1
429
430	IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN
431	IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NE quadrant
432	! ! i=I i=I+1/2
433	ze01 = pev(ii ,ij ) ; ze11 = pef(ii ,ij ) ! j=J+1/2 V ------- F
434	ze00 = pet(ii ,ij ) ; ze10 = peu(ii ,ij ) ! j=J T ------- U
435	zi = 2._wp * zi
436	zj = 2._wp * zj
437	ELSE ! SE quadrant
438	! ! i=I i=I+1/2
439	ze01 = pet(ii ,ij+1) ; ze11 = peu(ii ,ij+1) ! j=J+1 T ------- U
440	ze00 = pev(ii ,ij ) ; ze10 = pef(ii ,ij ) ! j=J+1/2 V ------- F
441	zi = 2._wp * zi
442	zj = 2._wp * (zj-0.5_wp)
443	ENDIF
444	ELSE
445	IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NW quadrant
446	! ! i=I i=I+1/2
447	ze01 = pef(ii ,ij ) ; ze11 = pev(ii+1,ij) ! j=J+1/2 F ------- V
448	ze00 = peu(ii ,ij ) ; ze10 = pet(ii+1,ij) ! j=J U ------- T
449	zi = 2._wp * (zi-0.5_wp)
450	zj = 2._wp * zj
451	ELSE ! SW quadrant
452	! ! i=I+1/2 i=I+1
453	ze01 = peu(ii ,ij+1) ; ze11 = pet(ii+1,ij+1) ! j=J+1 U ------- T
454	ze00 = pef(ii ,ij ) ; ze10 = pev(ii+1,ij ) ! j=J+1/2 F ------- V
455	zi = 2._wp * (zi-0.5_wp)
456	zj = 2._wp * (zj-0.5_wp)
457	ENDIF
458	ENDIF
459	!
460	icb_utl_bilin_e = ( ze01 * (1.-zi) + ze11 * zi ) * zj &
461	& + ( ze00 * (1.-zi) + ze10 * zi ) * (1.-zj)
462	!
463	END FUNCTION icb_utl_bilin_e
464
465
466	SUBROUTINE icb_utl_add( bergvals, ptvals )
467	!!----------------------------------------------------------------------
468	!! * ROUTINE icb_utl_add *
469	!!
470	!! ** Purpose : add a new berg to the iceberg list
471	!!
472	!!----------------------------------------------------------------------
473	TYPE(iceberg), INTENT(in) :: bergvals
474	TYPE(point) , INTENT(in) :: ptvals
475	!
476	TYPE(iceberg), POINTER :: new => NULL()
477	!!----------------------------------------------------------------------
478	!
479	new => NULL()
480	CALL icb_utl_create( new, bergvals, ptvals )
481	CALL icb_utl_insert( new )
482	new => NULL() ! Clear new
483	!
484	END SUBROUTINE icb_utl_add
485
486
487	SUBROUTINE icb_utl_create( berg, bergvals, ptvals )
488	!!----------------------------------------------------------------------
489	!! * ROUTINE icb_utl_create *
490	!!
491	!! ** Purpose : add a new berg to the iceberg list
492	!!
493	!!----------------------------------------------------------------------
494	TYPE(iceberg), INTENT(in) :: bergvals
495	TYPE(point) , INTENT(in) :: ptvals
496	TYPE(iceberg), POINTER :: berg
497	!
498	TYPE(point) , POINTER :: pt
499	INTEGER :: istat
500	!!----------------------------------------------------------------------
501	!
502	IF( ASSOCIATED(berg) ) CALL ctl_stop( 'icebergs, icb_utl_create: berg already associated' )
503	ALLOCATE(berg, STAT=istat)
504	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
505	berg%number(:) = bergvals%number(:)
506	berg%mass_scaling = bergvals%mass_scaling
507	berg%prev => NULL()
508	berg%next => NULL()
509	!
510	ALLOCATE(pt, STAT=istat)
511	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
512	pt = ptvals
513	berg%current_point => pt
514	!
515	END SUBROUTINE icb_utl_create
516
517
518	SUBROUTINE icb_utl_insert( newberg )
519	!!----------------------------------------------------------------------
520	!! * ROUTINE icb_utl_insert *
521	!!
522	!! ** Purpose : add a new berg to the iceberg list
523	!!
524	!!----------------------------------------------------------------------
525	TYPE(iceberg), POINTER :: newberg
526	!
527	TYPE(iceberg), POINTER :: this, prev, last
528	!!----------------------------------------------------------------------
529	!
530	IF( ASSOCIATED( first_berg ) ) THEN
531	last => first_berg
532	DO WHILE (ASSOCIATED(last%next))
533	last => last%next
534	ENDDO
535	newberg%prev => last
536	last%next => newberg
537	last => newberg
538	ELSE ! list is empty so create it
539	first_berg => newberg
540	ENDIF
541	!
542	END SUBROUTINE icb_utl_insert
543
544
545	REAL(wp) FUNCTION icb_utl_yearday(kmon, kday, khr, kmin, ksec)
546	!!----------------------------------------------------------------------
547	!! * FUNCTION icb_utl_yearday *
548	!!
549	!! ** Purpose :
550	!!
551	! sga - improved but still only applies to 365 day year, need to do this properly
552	!
553	!!gm all these info are already known in daymod, no???
554	!!
555	!!----------------------------------------------------------------------
556	INTEGER, INTENT(in) :: kmon, kday, khr, kmin, ksec
557	!
558	INTEGER, DIMENSION(12) :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
559	!!----------------------------------------------------------------------
560	!
561	icb_utl_yearday = REAL( SUM( imonths(1:kmon) ), wp )
562	icb_utl_yearday = icb_utl_yearday + REAL(kday-1,wp) + (REAL(khr,wp) + (REAL(kmin,wp) + REAL(ksec,wp)/60.)/60.)/24.
563	!
564	END FUNCTION icb_utl_yearday
565
566	!!-------------------------------------------------------------------------
567
568	SUBROUTINE icb_utl_delete( first, berg )
569	!!----------------------------------------------------------------------
570	!! * ROUTINE icb_utl_delete *
571	!!
572	!! ** Purpose :
573	!!
574	!!----------------------------------------------------------------------
575	TYPE(iceberg), POINTER :: first, berg
576	!!----------------------------------------------------------------------
577	! Connect neighbors to each other
578	IF ( ASSOCIATED(berg%prev) ) THEN
579	berg%prev%next => berg%next
580	ELSE
581	first => berg%next
582	ENDIF
583	IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
584	!
585	CALL icb_utl_destroy(berg)
586	!
587	END SUBROUTINE icb_utl_delete
588
589
590	SUBROUTINE icb_utl_destroy( berg )
591	!!----------------------------------------------------------------------
592	!! * ROUTINE icb_utl_destroy *
593	!!
594	!! ** Purpose : remove a single iceberg instance
595	!!
596	!!----------------------------------------------------------------------
597	TYPE(iceberg), POINTER :: berg
598	!!----------------------------------------------------------------------
599	!
600	! Remove any points
601	IF( ASSOCIATED( berg%current_point ) ) DEALLOCATE( berg%current_point )
602	!
603	DEALLOCATE(berg)
604	!
605	END SUBROUTINE icb_utl_destroy
606
607
608	SUBROUTINE icb_utl_track( knum, cd_label, kt )
609	!!----------------------------------------------------------------------
610	!! * ROUTINE icb_utl_track *
611	!!
612	!! ** Purpose :
613	!!
614	!!----------------------------------------------------------------------
615	INTEGER, DIMENSION(nkounts) :: knum ! iceberg number
616	CHARACTER(len=*) :: cd_label !
617	INTEGER :: kt ! timestep number
618	!
619	TYPE(iceberg), POINTER :: this
620	LOGICAL :: match
621	INTEGER :: k
622	!!----------------------------------------------------------------------
623	!
624	this => first_berg
625	DO WHILE( ASSOCIATED(this) )
626	match = .TRUE.
627	DO k = 1, nkounts
628	IF( this%number(k) /= knum(k) ) match = .FALSE.
629	END DO
630	IF( match ) CALL icb_utl_print_berg(this, kt)
631	this => this%next
632	END DO
633	!
634	END SUBROUTINE icb_utl_track
635
636
637	SUBROUTINE icb_utl_print_berg( berg, kt )
638	!!----------------------------------------------------------------------
639	!! * ROUTINE icb_utl_print_berg *
640	!!
641	!! ** Purpose : print one
642	!!
643	!!----------------------------------------------------------------------
644	TYPE(iceberg), POINTER :: berg
645	TYPE(point) , POINTER :: pt
646	INTEGER :: kt ! timestep number
647	!!----------------------------------------------------------------------
648	!
649	pt => berg%current_point
650	WRITE(numicb, 9200) kt, berg%number(1), &
651	pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel, &
652	pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
653	CALL flush( numicb )
654	9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
655	!
656	END SUBROUTINE icb_utl_print_berg
657
658
659	SUBROUTINE icb_utl_print( cd_label, kt )
660	!!----------------------------------------------------------------------
661	!! * ROUTINE icb_utl_print *
662	!!
663	!! ** Purpose : print many
664	!!
665	!!----------------------------------------------------------------------
666	CHARACTER(len=*) :: cd_label
667	INTEGER :: kt ! timestep number
668	!
669	INTEGER :: ibergs, inbergs
670	TYPE(iceberg), POINTER :: this
671	!!----------------------------------------------------------------------
672	!
673	this => first_berg
674	IF( ASSOCIATED(this) ) THEN
675	WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
676	WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) &
677	& 'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
678	ENDIF
679	DO WHILE( ASSOCIATED(this) )
680	CALL icb_utl_print_berg(this, kt)
681	this => this%next
682	END DO
683	ibergs = icb_utl_count()
684	inbergs = ibergs
685	IF( lk_mpp ) CALL mpp_sum(inbergs)
686	IF( ibergs > 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &
687	& cd_label, ibergs, inbergs, narea
688	!
689	END SUBROUTINE icb_utl_print
690
691
692	SUBROUTINE icb_utl_incr()
693	!!----------------------------------------------------------------------
694	!! * ROUTINE icb_utl_incr *
695	!!
696	!! ** Purpose :
697	!!
698	! Small routine for coping with very large integer values labelling icebergs
699	! num_bergs is a array of integers
700	! the first member is incremented in steps of jpnij starting from narea
701	! this means each iceberg is labelled with a unique number
702	! when this gets to the maximum allowed integer the second and subsequent members are
703	! used to count how many times the member before cycles
704	!!----------------------------------------------------------------------
705	INTEGER :: ii, ibig
706	!!----------------------------------------------------------------------
707
708	ibig = HUGE(num_bergs(1))
709	IF( ibig-jpnij < num_bergs(1) ) THEN
710	num_bergs(1) = narea
711	DO ii = 2,nkounts
712	IF( num_bergs(ii) == ibig ) THEN
713	num_bergs(ii) = 0
714	IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
715	ELSE
716	num_bergs(ii) = num_bergs(ii) + 1
717	EXIT
718	ENDIF
719	END DO
720	ELSE
721	num_bergs(1) = num_bergs(1) + jpnij
722	ENDIF
723	!
724	END SUBROUTINE icb_utl_incr
725
726
727	INTEGER FUNCTION icb_utl_count()
728	!!----------------------------------------------------------------------
729	!! * FUNCTION icb_utl_count *
730	!!
731	!! ** Purpose :
732	!!----------------------------------------------------------------------
733	TYPE(iceberg), POINTER :: this
734	!!----------------------------------------------------------------------
735	!
736	icb_utl_count = 0
737	this => first_berg
738	DO WHILE( ASSOCIATED(this) )
739	icb_utl_count = icb_utl_count+1
740	this => this%next
741	END DO
742	!
743	END FUNCTION icb_utl_count
744
745
746	REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
747	!!----------------------------------------------------------------------
748	!! * FUNCTION icb_utl_mass *
749	!!
750	!! ** Purpose : compute the mass all iceberg, all berg bits or all bergs.
751	!!----------------------------------------------------------------------
752	TYPE(iceberg) , POINTER :: first
753	TYPE(point) , POINTER :: pt
754	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
755	!
756	TYPE(iceberg), POINTER :: this
757	!!----------------------------------------------------------------------
758	icb_utl_mass = 0._wp
759	this => first
760	!
761	IF( PRESENT( justbergs ) ) THEN
762	DO WHILE( ASSOCIATED( this ) )
763	pt => this%current_point
764	icb_utl_mass = icb_utl_mass + pt%mass * this%mass_scaling
765	this => this%next
766	END DO
767	ELSEIF( PRESENT(justbits) ) THEN
768	DO WHILE( ASSOCIATED( this ) )
769	pt => this%current_point
770	icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
771	this => this%next
772	END DO
773	ELSE
774	DO WHILE( ASSOCIATED( this ) )
775	pt => this%current_point
776	icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
777	this => this%next
778	END DO
779	ENDIF
780	!
781	END FUNCTION icb_utl_mass
782
783
784	REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
785	!!----------------------------------------------------------------------
786	!! * FUNCTION icb_utl_heat *
787	!!
788	!! ** Purpose : compute the heat in all iceberg, all bergies or all bergs.
789	!!----------------------------------------------------------------------
790	TYPE(iceberg) , POINTER :: first
791	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
792	!
793	TYPE(iceberg) , POINTER :: this
794	TYPE(point) , POINTER :: pt
795	!!----------------------------------------------------------------------
796	icb_utl_heat = 0._wp
797	this => first
798	!
799	IF( PRESENT( justbergs ) ) THEN
800	DO WHILE( ASSOCIATED( this ) )
801	pt => this%current_point
802	icb_utl_heat = icb_utl_heat + pt%mass * this%mass_scaling * pt%heat_density
803	this => this%next
804	END DO
805	ELSEIF( PRESENT(justbits) ) THEN
806	DO WHILE( ASSOCIATED( this ) )
807	pt => this%current_point
808	icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
809	this => this%next
810	END DO
811	ELSE
812	DO WHILE( ASSOCIATED( this ) )
813	pt => this%current_point
814	icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
815	this => this%next
816	END DO
817	ENDIF
818	!
819	END FUNCTION icb_utl_heat
820
821	!!======================================================================
822	END MODULE icbutl

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

Download in other formats: