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

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

Last change on this file since 3339 was 3339, checked in by sga, 12 years ago

NEMO branch dev_r3337_NOCS10_ICB: add new iceberg sub-directory ICB
File size: 27.6 KB
Line 
1MODULE icbutl
2
3   !!======================================================================
4   !!                       ***  MODULE  icbutl  ***
5   !! Ocean physics:  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   !!   interp_flds   :
14   !!   bilin         :
15   !!   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   copy_flds                ! routine called in xxx.F90 module
37   PUBLIC   interp_flds              ! routine called in xxx.F90 module
38   PUBLIC   bilin
39   PUBLIC   bilin_x
40   PUBLIC   bilin_e
41   PUBLIC   add_new_berg_to_list
42   PUBLIC   insert_berg_into_list
43   PUBLIC   delete_iceberg_from_list
44   PUBLIC   destroy_iceberg
45   PUBLIC   track_berg
46   PUBLIC   print_berg
47   PUBLIC   print_bergs
48   PUBLIC   count_bergs
49   PUBLIC   increment_kounter
50   PUBLIC   yearday
51   PUBLIC   sum_mass
52   PUBLIC   sum_heat
53
54   PRIVATE  create_iceberg
55
56   !!-------------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE copy_flds()
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE copy_flds  ***
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 copy_flds
108
109   !!-------------------------------------------------------------------------
110
111   SUBROUTINE interp_flds( pi, pe1, puo, pui, pua, pssh_i,                     &
112      &                    pj, pe2, pvo, pvi, pva, pssh_j, psst, pcn, phi, pff )
113      !!----------------------------------------------------------------------
114      !!                  ***  ROUTINE interp_flds  ***
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 = bilin_e( e1t, e1u, e1v, e1f, pi, pj )                 ! scale factors
137      pe2 = bilin_e( e2t, e2u, e2v, e2f, pi, pj )
138      !
139      puo = bilin( uo_e, pi, pj, 'U', 1, 1 )                      ! ocean velocities
140      pvo = bilin( vo_e, pi, pj, 'V', 1, 1 )
141      psst = bilin( sst_m, pi, pj, 'T', 0, 0 )                    ! SST
142      pcn  = bilin( fr_i , pi, pj, 'T', 0, 0 )                    ! ice concentration
143      pff  = bilin( ff_e , pi, pj, 'F', 1, 1 )                    ! Coriolis parameter
144      !
145      pua = bilin( ua_e , pi, pj, 'U', 1, 1 )                     ! 10m wind
146      pva = 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( pua*pua + pva*pva)  )  )
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 = bilin( ui_e, pi, pj, 'U', 1, 1 )                      ! sea-ice velocities
154      pvi = bilin( vi_e, pi, pj, 'V', 1, 1 )
155      phi = 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 = ( bilin( ssh_e, pi+0.1_wp, pj, 'T', 1, 1 ) - bilin( ssh_e, pi-0.1_wp, pj, 'T', 1, 1 )  ) / ( 0.2_wp * pe1 )
164      pssh_j = ( bilin( ssh_e, pi, pj+0.1_wp, 'T', 1, 1 ) - bilin( ssh_e, pi, pj-0.1_wp, 'T', 1, 1 )  ) / ( 0.2_wp * pe2 )
165      !
166   END SUBROUTINE interp_flds
167
168   !!-------------------------------------------------------------------------
169
170   REAL(wp) FUNCTION bilin( pfld, pi, pj, cd_type, jdi, jdj )
171      !!----------------------------------------------------------------------
172      !!                  ***  FUNCTION bilin  ***
173      !!
174      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
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      INTEGER                                         , INTENT(in) ::   jdi, jdj  ! extra halo on grid
181      REAL(wp), DIMENSION(1-jdi:jpi+jdi,1-jdj:jpj+jdj), INTENT(in) ::   pfld      ! field to be interpolated
182      REAL(wp)                                        , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
183      CHARACTER(len=1)                                , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
184      !
185      INTEGER  ::   ii, ij   ! local integer
186      REAL(wp) ::   zi, zj   ! local real
187      !!----------------------------------------------------------------------
188      !
189      SELECT CASE ( cd_type )
190         CASE ( 'T' )
191            ! note that here there is no +0.5 added
192            ! since we're looking for four T points containing quadrant we're in of
193            ! current T cell
194            ii = INT( pi     )
195            ij = INT( pj      )    ! T-point
196            zi = pi - REAL(ii,wp)
197            zj = pj - REAL(ij,wp)
198         CASE ( 'U' )
199            ii = INT( pi-0.5 )
200            ij = INT( pj      )    ! U-point
201            zi = pi - 0.5 - REAL(ii,wp)
202            zj = pj - REAL(ij,wp)
203         CASE ( 'V' )
204            ii = INT( pi     )
205            ij = INT( pj -0.5 )    ! V-point
206            zi = pi - REAL(ii,wp)
207            zj = pj - 0.5 - REAL(ij,wp)
208         CASE ( 'F' )
209            ii = INT( pi-0.5 )
210            ij = INT( pj -0.5 )    ! F-point
211            zi = pi - 0.5 - REAL(ii,wp)
212            zj = pj - 0.5 - REAL(ij,wp)
213      END SELECT
214      !
215      ! find position in this processor
216      ii = ii - nimpp + 1
217      ij = ij - njmpp + 1
218      !
219      bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
220         &  + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
221      !
222   END FUNCTION bilin
223
224   !!-------------------------------------------------------------------------
225
226   REAL(wp) FUNCTION bilin_x( pfld, pi, pj )
227      !!----------------------------------------------------------------------
228      !!                  ***  FUNCTION 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), DIMENSION(4)                   ::   z4
243      !!----------------------------------------------------------------------
244      !
245      ! note that here there is no +0.5 added
246      ! since we're looking for four T points containing quadrant we're in of
247      ! current T cell
248      ii = INT( pi     )
249      ij = INT( pj      )    ! T-point
250      zi = pi - REAL(ii,wp)
251      zj = pj - REAL(ij,wp)
252      !
253      ! find position in this processor
254      ii = ii - nimpp + 1
255      ij = ij - njmpp + 1
256      z4(1) = pfld(ii  ,ij  )
257      z4(2) = pfld(ii+1,ij  )
258      z4(3) = pfld(ii  ,ij+1)
259      z4(4) = pfld(ii+1,ij+1)
260      IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
261         WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
262      ENDIF
263      !
264      bilin_x = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
265      IF( bilin_x > 180._wp ) bilin_x = bilin_x - 360._wp
266      !
267   END FUNCTION bilin_x
268
269   !!-------------------------------------------------------------------------
270
271   REAL(wp) FUNCTION bilin_e( pet, peu, pev, pef, pi, pj )
272      !!----------------------------------------------------------------------
273      !!                  ***  FUNCTION dom_init  ***
274      !!
275      !! ** Purpose :   bilinear interpolation at berg location of horizontal scale factor
276      !! ** Method  :   interpolation done using the 4 nearest grid point among
277      !!                t-, u-, v-, and f-points.
278      !!----------------------------------------------------------------------
279      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
280      REAL(wp)                , INTENT(in) ::   pi, pj               ! targeted coordinates in (i,j) referential
281      !
282      INTEGER  ::   ii, ij, icase   ! local integer
283      !
284      ! weights corresponding to corner points of a T cell quadrant
285      REAL(wp) ::   zi, zj          ! local real
286      !
287      ! values at corner points of a T cell quadrant
288      ! 00 = bottom left, 10 = bottom right, 01 = top left, 11 = top right
289      REAL(wp) ::   ze00, ze10, ze01, ze11
290      !!----------------------------------------------------------------------
291      !
292      ii = INT( pi )   ;   ij = INT( pj )    ! left bottom T-point (i,j) indices
293
294      ! fractional box spacing
295      ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell
296      ! 0.5 <= zi < 1  , 0   <= zj < 0.5   -->  NE quadrant
297      ! 0   <= zi < 0.5, 0.5 <= zj < 1     -->  SE quadrant
298      ! 0.5 <= zi < 1  , 0.5 <= zj < 1     -->  SW quadrant
299
300      zi = pi - REAL(ii,wp)
301      zj = pj - REAL(ij,wp)
302
303      ! find position in this processor
304      ii = ii - nimpp + 1
305      ij = ij - njmpp + 1
306
307      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
308         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NE quadrant
309            !                                                      !             i=I       i=I+1/2
310            ze01 = pev(ii  ,ij  )   ;   ze11 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
311            ze00 = pet(ii  ,ij  )   ;   ze10 = peu(ii  ,ij  )      !   j=J        T ------- U
312            zi = 2._wp * zi
313            zj = 2._wp * zj
314         ELSE                                                      !  SE quadrant
315            !                                                                    !             i=I       i=I+1/2
316            ze01 = pet(ii  ,ij+1)   ;   ze11 = peu(ii  ,ij+1)      !   j=J+1      T ------- U
317            ze00 = pev(ii  ,ij  )   ;   ze10 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
318            zi = 2._wp *  zi
319            zj = 2._wp * (zj-0.5_wp)
320         ENDIF
321      ELSE
322         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NW quadrant
323            !                                                                    !             i=I       i=I+1/2
324            ze01 = pef(ii  ,ij  )   ;   ze11 = pev(ii+1,ij)        !   j=J+1/2    F ------- V
325            ze00 = peu(ii  ,ij  )   ;   ze10 = pet(ii+1,ij)        !   j=J        U ------- T
326            zi = 2._wp * (zi-0.5_wp)
327            zj = 2._wp *  zj
328         ELSE                                                      !  SW quadrant
329            !                                                                    !             i=I+1/2   i=I+1
330            ze01 = peu(ii  ,ij+1)   ;   ze11 = pet(ii+1,ij+1)      !   j=J+1      U ------- T
331            ze00 = pef(ii  ,ij  )   ;   ze10 = pev(ii+1,ij  )      !   j=J+1/2    F ------- V
332            zi = 2._wp * (zi-0.5_wp)
333            zj = 2._wp * (zj-0.5_wp)
334         ENDIF
335      ENDIF
336      !
337      bilin_e = ( ze01 * (1.-zi) + ze11 * zi ) *     zj    &
338         &    + ( ze00 * (1.-zi) + ze10 * zi ) * (1.-zj)
339      !
340   END FUNCTION bilin_e
341
342   !!-------------------------------------------------------------------------
343
344
345   SUBROUTINE add_new_berg_to_list( bergvals, ptvals )
346      !!----------------------------------------------------------------------
347      !!                ***  ROUTINE add_new_berg_to_list  ***
348      !!
349      !! ** Purpose :   add a new berg to the iceberg list
350      !!
351      !! ** method  : - ?
352      !!----------------------------------------------------------------------
353      TYPE(iceberg), INTENT(in)           ::   bergvals
354      TYPE(point)  , INTENT(in)           ::   ptvals
355      !
356      TYPE(iceberg), POINTER ::   new => NULL()
357      !!----------------------------------------------------------------------
358      !
359      new => NULL()
360      CALL create_iceberg( new, bergvals, ptvals )
361      CALL insert_berg_into_list( new )
362      new => NULL()     ! Clear new
363      !
364   END SUBROUTINE add_new_berg_to_list
365
366   !!-------------------------------------------------------------------------
367
368   SUBROUTINE create_iceberg(berg, bergvals, ptvals )
369      ! Arguments
370      TYPE(iceberg), INTENT(in) :: bergvals
371      TYPE(point)  , INTENT(in) :: ptvals
372      TYPE(iceberg), POINTER    :: berg
373      ! local variables
374      TYPE(point)  , POINTER    :: pt
375      INTEGER                   :: istat
376
377      IF ( ASSOCIATED(berg) ) then
378         CALL ctl_stop( 'icebergs, create_iceberg: berg already associated' )
379      ENDIF
380      ALLOCATE(berg, STAT=istat)
381      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
382      berg%number(:) = bergvals%number(:)
383      berg%mass_scaling = bergvals%mass_scaling
384      berg%prev => NULL()
385      berg%next => NULL()
386
387      ALLOCATE(pt, STAT=istat)
388      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
389      pt = ptvals
390      berg%current_point => pt
391
392   END SUBROUTINE create_iceberg
393
394   !!-------------------------------------------------------------------------
395
396   SUBROUTINE insert_berg_into_list( newberg )
397      !!----------------------------------------------------------------------
398      !!                 ***  ROUTINE insert_berg_into_list  ***
399      !!
400      !! ** Purpose :   add a new berg to the iceberg list
401      !!
402      !! ** method  : - ?
403      !!----------------------------------------------------------------------
404      TYPE(iceberg), POINTER  ::   newberg
405      !
406      TYPE(iceberg), POINTER  ::   this, prev, last
407      !!----------------------------------------------------------------------
408
409      IF( ASSOCIATED( first_berg ) ) THEN
410!        last = last_berg()
411         last=>first_berg
412         DO WHILE (ASSOCIATED(last%next))
413            last=>last%next
414         ENDDO
415         newberg%prev => last
416         last%next    => newberg
417         last         => newberg
418      ELSE                       ! list is empty so create it
419         first_berg => newberg
420      ENDIF
421      !
422   END SUBROUTINE insert_berg_into_list
423
424   !!-------------------------------------------------------------------------
425
426   REAL(wp) FUNCTION yearday(imon, iday, ihr, imin, isec)
427      ! sga - improved but still only applies to 365 day year, need to do this properly
428      ! Arguments
429      INTEGER, intent(in)     :: imon, iday, ihr, imin, isec
430      ! Local variables
431      INTEGER                 :: i
432      INTEGER, DIMENSION(12)  :: months = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
433
434      yearday = FLOAT( SUM( months(1:imon) ) )
435      yearday = yearday + FLOAT(iday-1) + (FLOAT(ihr) + (FLOAT(imin) + FLOAT(isec)/60.)/60.)/24.
436
437   END FUNCTION yearday
438
439   !!-------------------------------------------------------------------------
440
441   SUBROUTINE delete_iceberg_from_list( first, berg )
442      ! Arguments
443      TYPE(iceberg), POINTER :: first, berg
444
445      ! Connect neighbors to each other
446      IF ( ASSOCIATED(berg%prev) ) THEN
447        berg%prev%next => berg%next
448      ELSE
449        first => berg%next
450      ENDIF
451      IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
452
453      ! Bye-bye berg
454      CALL destroy_iceberg(berg)
455
456   END SUBROUTINE delete_iceberg_from_list
457
458   !!-------------------------------------------------------------------------
459
460   SUBROUTINE destroy_iceberg(berg)
461      ! Arguments
462      TYPE(iceberg), POINTER :: berg
463
464      ! Remove any points
465      IF( ASSOCIATED( berg%current_point ) ) DEALLOCATE( berg%current_point )
466
467      ! Bye-bye berg
468      DEALLOCATE(berg)
469
470   END SUBROUTINE destroy_iceberg
471
472   !!-------------------------------------------------------------------------
473
474   SUBROUTINE track_berg( num, label, kt )
475      ! Arguments
476      INTEGER, DIMENSION(nkounts)    :: num            ! iceberg number
477      CHARACTER(len=*)               :: label
478      INTEGER                        :: kt             ! timestep number
479      ! Local variables
480      TYPE(iceberg), POINTER         :: this
481      LOGICAL                        :: match
482      INTEGER                        :: k
483
484      this => first_berg
485      DO WHILE( ASSOCIATED(this) )
486         match = .TRUE.
487         DO k=1,nkounts
488            IF( this%number(k) /= num(k) ) match = .FALSE.
489         END DO
490         IF( match ) CALL print_berg(this, kt)
491         this => this%next
492      ENDDO
493
494   END SUBROUTINE track_berg
495
496   !!-------------------------------------------------------------------------
497
498   SUBROUTINE print_berg( berg, kt )
499      ! Arguments
500      TYPE(iceberg), POINTER :: berg
501      TYPE(point)  , POINTER :: pt
502      INTEGER                :: kt             ! timestep number
503
504      pt => berg%current_point
505      WRITE(numicb, 9200) kt, berg%number(1), &
506                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
507                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
508      CALL flush( numicb )
509 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
510
511   END SUBROUTINE print_berg
512
513   !!-------------------------------------------------------------------------
514
515   SUBROUTINE print_bergs( label, kt )
516      ! Arguments
517      CHARACTER(len=*)       :: label
518      INTEGER                :: kt             ! timestep number
519      ! Local variables
520      INTEGER                :: nbergs, nnbergs
521      TYPE(iceberg), POINTER :: this
522
523      this => first_berg
524      IF( ASSOCIATED(this) ) THEN
525         WRITE(numicb,'(a," pe=(",i3,")")' ) label, narea
526         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
527                      'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
528      ENDIF
529      DO WHILE( ASSOCIATED(this) )
530        CALL print_berg(this, kt)
531        this => this%next
532      ENDDO
533      nbergs = count_bergs()
534      nnbergs = nbergs
535      IF( lk_mpp ) CALL mpp_sum(nnbergs)
536      IF ( nbergs .GT. 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &
537                                           label, nbergs, nnbergs, narea
538
539   END SUBROUTINE print_bergs
540
541   !!-------------------------------------------------------------------------
542
543   SUBROUTINE increment_kounter()
544      ! Small routine for coping with very large integer values labelling icebergs
545      ! kount_bergs is a array of integers
546      ! the first member is incremented in steps of jpnij starting from narea
547      ! this means each iceberg is labelled with a unique number
548      ! when this gets to the maximum allowed integer the second and subsequent members are
549      ! used to count how many times the member before cycles
550
551      ! Local variables
552      INTEGER                :: i, ibig
553
554      ibig = HUGE(kount_bergs(1))
555      IF( ibig-jpnij < kount_bergs(1) ) THEN
556         kount_bergs(1) = narea
557         DO i = 2,nkounts
558            IF( kount_bergs(i) == ibig ) THEN
559               kount_bergs(i) = 0
560               IF( i == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
561            ELSE
562               kount_bergs(i) = kount_bergs(i) + 1
563               EXIT
564            ENDIF
565         END DO
566      ELSE
567         kount_bergs(1) = kount_bergs(1) + jpnij
568      ENDIF
569
570   END SUBROUTINE increment_kounter
571
572   !!-------------------------------------------------------------------------
573
574   INTEGER function count_bergs()
575      ! Local variables
576      TYPE(iceberg), POINTER :: this
577
578      count_bergs = 0
579      this => first_berg
580      DO WHILE( ASSOCIATED(this) )
581         count_bergs = count_bergs+1
582         this => this%next
583      ENDDO
584
585   END FUNCTION count_bergs
586
587   !!-------------------------------------------------------------------------
588
589   REAL(wp) FUNCTION sum_mass( first, justbits, justbergs )
590      !!----------------------------------------------------------------------
591      !!                 ***  FUNCTION sum_mass  ***
592      !!
593      !! ** Purpose :   compute the mass all iceberg, all bergies or all bergs.
594      !!----------------------------------------------------------------------
595      TYPE(iceberg)      , POINTER  ::   first
596      TYPE(point)        , POINTER  ::   pt
597      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
598      !
599      TYPE(iceberg), POINTER ::   this
600      !!----------------------------------------------------------------------
601      sum_mass = 0._wp
602      this => first
603      !
604      IF( PRESENT( justbergs  ) ) THEN
605         DO WHILE( ASSOCIATED( this ) )
606            pt => this%current_point
607            sum_mass = sum_mass + pt%mass         * this%mass_scaling
608            this => this%next
609         END DO
610      ELSEIF( PRESENT(justbits) ) THEN
611         DO WHILE( ASSOCIATED( this ) )
612            pt => this%current_point
613            sum_mass = sum_mass + pt%mass_of_bits * this%mass_scaling
614            this => this%next
615         END DO
616      ELSE
617         DO WHILE( ASSOCIATED( this ) )
618            pt => this%current_point
619            sum_mass = sum_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
620            this => this%next
621         END DO
622      ENDIF
623      !
624   END FUNCTION sum_mass
625
626   !!-------------------------------------------------------------------------
627
628   REAL(wp) FUNCTION sum_heat( first, justbits, justbergs )
629      !!----------------------------------------------------------------------
630      !!                 ***  FUNCTION sum_heat  ***
631      !!
632      !! ** Purpose :   compute the heat in all iceberg, all bergies or all bergs.
633      !!----------------------------------------------------------------------
634      TYPE(iceberg)      , POINTER  ::   first
635      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
636      !
637      TYPE(iceberg)      , POINTER  ::   this
638      TYPE(point)        , POINTER  ::   pt
639      !!----------------------------------------------------------------------
640      sum_heat = 0._wp
641      this => first
642      !
643      IF( PRESENT( justbergs  ) ) THEN
644         DO WHILE( ASSOCIATED( this ) )
645            pt => this%current_point
646            sum_heat = sum_heat + pt%mass         * this%mass_scaling * pt%heat_density
647            this => this%next
648         END DO
649      ELSEIF( PRESENT(justbits) ) THEN
650         DO WHILE( ASSOCIATED( this ) )
651            pt => this%current_point
652            sum_heat = sum_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
653            this => this%next
654         END DO
655      ELSE
656         DO WHILE( ASSOCIATED( this ) )
657            pt => this%current_point
658            sum_heat = sum_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
659            this => this%next
660         END DO
661      ENDIF
662      !
663   END FUNCTION sum_heat
664
665   !!-------------------------------------------------------------------------
666
667END MODULE icbutl
Note: See TracBrowser for help on using the repository browser.