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/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB – NEMO

Context Navigation

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 3379

Last change on this file since 3379 was 3379, checked in by sga, 12 years ago
NEMO branch dev_r3337_NOCS10_ICB: remove deallocation of diagnostic arrays and add more comments
File size: 30.6 KB

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

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

Download in other formats: