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 @ 3372

Last change on this file since 3372 was 3372, checked in by sga, 12 years ago
NEMO branch dev_r3337_NOCS10_ICB: change all routine names and create more Gurvanistic havoc
File size: 30.8 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 icbrun 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 icbrun 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, icbrun 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, icbrun 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 : iceberg initialization.
117	!!
118	!! ** Method : - blah blah
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 = ii - nimpp + 1
218	ij = ij - njmpp + 1
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 = ii - nimpp + 1
256	ij = ij - njmpp + 1
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 point 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 = ii - nimpp + 1
306	ij = ij - njmpp + 1
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	!! ** method : - ?
351	!!----------------------------------------------------------------------
352	TYPE(iceberg), INTENT(in) :: bergvals
353	TYPE(point) , INTENT(in) :: ptvals
354	!
355	TYPE(iceberg), POINTER :: new => NULL()
356	!!----------------------------------------------------------------------
357	!
358	new => NULL()
359	CALL icb_utl_create( new, bergvals, ptvals )
360	CALL icb_utl_insert( new )
361	new => NULL() ! Clear new
362	!
363	END SUBROUTINE icb_utl_add
364
365
366	SUBROUTINE icb_utl_create( berg, bergvals, ptvals )
367	!!----------------------------------------------------------------------
368	!! * ROUTINE icb_utl_create *
369	!!
370	!! ** Purpose : add a new berg to the iceberg list
371	!!
372	!! ** method : - ?
373	!!----------------------------------------------------------------------
374	TYPE(iceberg), INTENT(in) :: bergvals
375	TYPE(point) , INTENT(in) :: ptvals
376	TYPE(iceberg), POINTER :: berg
377	!
378	TYPE(point) , POINTER :: pt
379	INTEGER :: istat
380	!!----------------------------------------------------------------------
381	!
382	IF( ASSOCIATED(berg) ) CALL ctl_stop( 'icebergs, icb_utl_create: berg already associated' )
383	ALLOCATE(berg, STAT=istat)
384	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
385	berg%number(:) = bergvals%number(:)
386	berg%mass_scaling = bergvals%mass_scaling
387	berg%prev => NULL()
388	berg%next => NULL()
389	!
390	ALLOCATE(pt, STAT=istat)
391	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
392	pt = ptvals
393	berg%current_point => pt
394	!
395	END SUBROUTINE icb_utl_create
396
397
398	SUBROUTINE icb_utl_insert( newberg )
399	!!----------------------------------------------------------------------
400	!! * ROUTINE icb_utl_insert *
401	!!
402	!! ** Purpose : add a new berg to the iceberg list
403	!!
404	!! ** method : - ?
405	!!----------------------------------------------------------------------
406	TYPE(iceberg), POINTER :: newberg
407	!
408	TYPE(iceberg), POINTER :: this, prev, last
409	!!----------------------------------------------------------------------
410	!
411	IF( ASSOCIATED( first_berg ) ) THEN
412	! last = last_berg()
413	last=>first_berg
414	DO WHILE (ASSOCIATED(last%next))
415	last=>last%next
416	ENDDO
417	newberg%prev => last
418	last%next => newberg
419	last => newberg
420	ELSE ! list is empty so create it
421	first_berg => newberg
422	ENDIF
423	!
424	END SUBROUTINE icb_utl_insert
425
426
427	REAL(wp) FUNCTION icb_utl_yearday(kmon, kday, khr, kmin, ksec)
428	!!----------------------------------------------------------------------
429	!! * FUNCTION icb_utl_yearday *
430	!!
431	!! ** Purpose :
432	!!
433	! sga - improved but still only applies to 365 day year, need to do this properly
434	!
435	!!gm all these info are already known in daymod, no???
436	!!
437	!! ** method : - ?
438	!!----------------------------------------------------------------------
439	INTEGER, INTENT(in) :: kmon, kday, khr, kmin, ksec
440	!
441	INTEGER, DIMENSION(12) :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
442	!!----------------------------------------------------------------------
443	!
444	icb_utl_yearday = FLOAT( SUM( imonths(1:kmon) ) )
445	icb_utl_yearday = icb_utl_yearday + FLOAT(kday-1) + (FLOAT(khr) + (FLOAT(kmin) + FLOAT(ksec)/60.)/60.)/24.
446	!
447	END FUNCTION icb_utl_yearday
448
449	!!-------------------------------------------------------------------------
450
451	SUBROUTINE icb_utl_delete( first, berg )
452	!!----------------------------------------------------------------------
453	!! * ROUTINE icb_utl_delete *
454	!!
455	!! ** Purpose :
456	!!
457	!!----------------------------------------------------------------------
458	TYPE(iceberg), POINTER :: first, berg
459	!!----------------------------------------------------------------------
460	! Connect neighbors to each other
461	IF ( ASSOCIATED(berg%prev) ) THEN
462	berg%prev%next => berg%next
463	ELSE
464	first => berg%next
465	ENDIF
466	IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
467	!
468	! Bye-bye berg
469	CALL icb_utl_destroy(berg)
470	!
471	END SUBROUTINE icb_utl_delete
472
473
474	SUBROUTINE icb_utl_destroy( berg )
475	!!----------------------------------------------------------------------
476	!! * ROUTINE icb_utl_destroy *
477	!!
478	!! ** Purpose :
479	!!
480	!!----------------------------------------------------------------------
481	TYPE(iceberg), POINTER :: berg
482	!!----------------------------------------------------------------------
483	!
484	! Remove any points
485	IF( ASSOCIATED( berg%current_point ) ) DEALLOCATE( berg%current_point )
486	!
487	! Bye-bye berg
488	DEALLOCATE(berg)
489	!
490	END SUBROUTINE icb_utl_destroy
491
492
493	SUBROUTINE icb_utl_track( knum, cd_label, kt )
494	!!----------------------------------------------------------------------
495	!! * ROUTINE icb_utl_track *
496	!!
497	!! ** Purpose :
498	!!
499	!!----------------------------------------------------------------------
500	INTEGER, DIMENSION(nkounts) :: knum ! iceberg number
501	CHARACTER(len=*) :: cd_label !
502	INTEGER :: kt ! timestep number
503	!
504	TYPE(iceberg), POINTER :: this
505	LOGICAL :: match
506	INTEGER :: k
507	!!----------------------------------------------------------------------
508	!
509	this => first_berg
510	DO WHILE( ASSOCIATED(this) )
511	match = .TRUE.
512	DO k = 1, nkounts
513	IF( this%number(k) /= knum(k) ) match = .FALSE.
514	END DO
515	IF( match ) CALL icb_utl_print_berg(this, kt)
516	this => this%next
517	END DO
518	!
519	END SUBROUTINE icb_utl_track
520
521
522	SUBROUTINE icb_utl_print_berg( berg, kt )
523	!!----------------------------------------------------------------------
524	!! * ROUTINE icb_utl_print_berg *
525	!!
526	!! ** Purpose :
527	!!
528	!!----------------------------------------------------------------------
529	TYPE(iceberg), POINTER :: berg
530	TYPE(point) , POINTER :: pt
531	INTEGER :: kt ! timestep number
532	!!----------------------------------------------------------------------
533	!
534	pt => berg%current_point
535	WRITE(numicb, 9200) kt, berg%number(1), &
536	pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel, &
537	pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
538	CALL flush( numicb )
539	9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
540	!
541	END SUBROUTINE icb_utl_print_berg
542
543
544	SUBROUTINE icb_utl_print( cd_label, kt )
545	!!----------------------------------------------------------------------
546	!! * ROUTINE icb_utl_print *
547	!!
548	!! ** Purpose :
549	!!
550	!!----------------------------------------------------------------------
551	CHARACTER(len=*) :: cd_label
552	INTEGER :: kt ! timestep number
553	!
554	INTEGER :: ibergs, inbergs
555	TYPE(iceberg), POINTER :: this
556	!!----------------------------------------------------------------------
557	!
558	this => first_berg
559	IF( ASSOCIATED(this) ) THEN
560	WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
561	WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) &
562	& 'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
563	ENDIF
564	DO WHILE( ASSOCIATED(this) )
565	CALL icb_utl_print_berg(this, kt)
566	this => this%next
567	END DO
568	ibergs = icb_utl_count()
569	inbergs = ibergs
570	IF( lk_mpp ) CALL mpp_sum(inbergs)
571	IF( ibergs > 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &
572	& cd_label, ibergs, inbergs, narea
573	!
574	END SUBROUTINE icb_utl_print
575
576
577	SUBROUTINE icb_utl_incr()
578	!!----------------------------------------------------------------------
579	!! * ROUTINE icb_utl_incr *
580	!!
581	!! ** Purpose :
582	!!
583	! Small routine for coping with very large integer values labelling icebergs
584	! num_bergs is a array of integers
585	! the first member is incremented in steps of jpnij starting from narea
586	! this means each iceberg is labelled with a unique number
587	! when this gets to the maximum allowed integer the second and subsequent members are
588	! used to count how many times the member before cycles
589	!!----------------------------------------------------------------------
590	INTEGER :: ii, ibig
591	!!----------------------------------------------------------------------
592
593	ibig = HUGE(num_bergs(1))
594	IF( ibig-jpnij < num_bergs(1) ) THEN
595	num_bergs(1) = narea
596	DO ii = 2,nkounts
597	IF( num_bergs(ii) == ibig ) THEN
598	num_bergs(ii) = 0
599	IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
600	ELSE
601	num_bergs(ii) = num_bergs(ii) + 1
602	EXIT
603	ENDIF
604	END DO
605	ELSE
606	num_bergs(1) = num_bergs(1) + jpnij
607	ENDIF
608	!
609	END SUBROUTINE icb_utl_incr
610
611
612	INTEGER FUNCTION icb_utl_count()
613	!!----------------------------------------------------------------------
614	!! * FUNCTION icb_utl_count *
615	!!
616	!! ** Purpose :
617	!!----------------------------------------------------------------------
618	TYPE(iceberg), POINTER :: this
619	!!----------------------------------------------------------------------
620	!
621	icb_utl_count = 0
622	this => first_berg
623	DO WHILE( ASSOCIATED(this) )
624	icb_utl_count = icb_utl_count+1
625	this => this%next
626	END DO
627	!
628	END FUNCTION icb_utl_count
629
630
631	REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
632	!!----------------------------------------------------------------------
633	!! * FUNCTION icb_utl_mass *
634	!!
635	!! ** Purpose : compute the mass all iceberg, all bergies or all bergs.
636	!!----------------------------------------------------------------------
637	TYPE(iceberg) , POINTER :: first
638	TYPE(point) , POINTER :: pt
639	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
640	!
641	TYPE(iceberg), POINTER :: this
642	!!----------------------------------------------------------------------
643	icb_utl_mass = 0._wp
644	this => first
645	!
646	IF( PRESENT( justbergs ) ) THEN
647	DO WHILE( ASSOCIATED( this ) )
648	pt => this%current_point
649	icb_utl_mass = icb_utl_mass + pt%mass * this%mass_scaling
650	this => this%next
651	END DO
652	ELSEIF( PRESENT(justbits) ) THEN
653	DO WHILE( ASSOCIATED( this ) )
654	pt => this%current_point
655	icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
656	this => this%next
657	END DO
658	ELSE
659	DO WHILE( ASSOCIATED( this ) )
660	pt => this%current_point
661	icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
662	this => this%next
663	END DO
664	ENDIF
665	!
666	END FUNCTION icb_utl_mass
667
668
669	REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
670	!!----------------------------------------------------------------------
671	!! * FUNCTION icb_utl_heat *
672	!!
673	!! ** Purpose : compute the heat in all iceberg, all bergies or all bergs.
674	!!----------------------------------------------------------------------
675	TYPE(iceberg) , POINTER :: first
676	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
677	!
678	TYPE(iceberg) , POINTER :: this
679	TYPE(point) , POINTER :: pt
680	!!----------------------------------------------------------------------
681	icb_utl_heat = 0._wp
682	this => first
683	!
684	IF( PRESENT( justbergs ) ) THEN
685	DO WHILE( ASSOCIATED( this ) )
686	pt => this%current_point
687	icb_utl_heat = icb_utl_heat + pt%mass * this%mass_scaling * pt%heat_density
688	this => this%next
689	END DO
690	ELSEIF( PRESENT(justbits) ) THEN
691	DO WHILE( ASSOCIATED( this ) )
692	pt => this%current_point
693	icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
694	this => this%next
695	END DO
696	ELSE
697	DO WHILE( ASSOCIATED( this ) )
698	pt => this%current_point
699	icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
700	this => this%next
701	END DO
702	ENDIF
703	!
704	END FUNCTION icb_utl_heat
705
706	!!======================================================================
707	END MODULE icbutl

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

Download in other formats: