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

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

NEMO branch dev_r3337_NOCS10_ICB: remove unnecessary deallocation of arrays at end of iceberg run

File size: 30.6 KB
Line 
1MODULE 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
58CONTAINS
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( 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 = 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 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 = 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      !! ** 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      CALL icb_utl_destroy(berg)
469      !
470   END SUBROUTINE icb_utl_delete
471
472
473   SUBROUTINE icb_utl_destroy( berg )
474      !!----------------------------------------------------------------------
475      !!                 ***  ROUTINE icb_utl_destroy  ***
476      !!
477      !! ** Purpose :   
478      !!
479      !!----------------------------------------------------------------------
480      TYPE(iceberg), POINTER :: berg
481      !!----------------------------------------------------------------------
482      !
483      ! Remove any points
484      IF( ASSOCIATED( berg%current_point ) )   DEALLOCATE( berg%current_point )
485      !
486      DEALLOCATE(berg)
487      !
488   END SUBROUTINE icb_utl_destroy
489
490
491   SUBROUTINE icb_utl_track( knum, cd_label, kt )
492      !!----------------------------------------------------------------------
493      !!                 ***  ROUTINE icb_utl_track  ***
494      !!
495      !! ** Purpose :   
496      !!
497      !!----------------------------------------------------------------------
498      INTEGER, DIMENSION(nkounts)    :: knum       ! iceberg number
499      CHARACTER(len=*)               :: cd_label   !
500      INTEGER                        :: kt         ! timestep number
501      !
502      TYPE(iceberg), POINTER         :: this
503      LOGICAL                        :: match
504      INTEGER                        :: k
505      !!----------------------------------------------------------------------
506      !
507      this => first_berg
508      DO WHILE( ASSOCIATED(this) )
509         match = .TRUE.
510         DO k = 1, nkounts
511            IF( this%number(k) /= knum(k) ) match = .FALSE.
512         END DO
513         IF( match )   CALL icb_utl_print_berg(this, kt)
514         this => this%next
515      END DO
516      !
517   END SUBROUTINE icb_utl_track
518
519
520   SUBROUTINE icb_utl_print_berg( berg, kt )
521      !!----------------------------------------------------------------------
522      !!                 ***  ROUTINE icb_utl_print_berg  ***
523      !!
524      !! ** Purpose :   
525      !!
526      !!----------------------------------------------------------------------
527      TYPE(iceberg), POINTER :: berg
528      TYPE(point)  , POINTER :: pt
529      INTEGER                :: kt      ! timestep number
530      !!----------------------------------------------------------------------
531      !
532      pt => berg%current_point
533      WRITE(numicb, 9200) kt, berg%number(1), &
534                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
535                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
536      CALL flush( numicb )
537 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
538      !
539   END SUBROUTINE icb_utl_print_berg
540
541
542   SUBROUTINE icb_utl_print( cd_label, kt )
543      !!----------------------------------------------------------------------
544      !!                 ***  ROUTINE icb_utl_print  ***
545      !!
546      !! ** Purpose :   
547      !!
548      !!----------------------------------------------------------------------
549      CHARACTER(len=*)       :: cd_label
550      INTEGER                :: kt             ! timestep number
551      !
552      INTEGER                :: ibergs, inbergs
553      TYPE(iceberg), POINTER :: this
554      !!----------------------------------------------------------------------
555      !
556      this => first_berg
557      IF( ASSOCIATED(this) ) THEN
558         WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
559         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
560            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
561      ENDIF
562      DO WHILE( ASSOCIATED(this) )
563        CALL icb_utl_print_berg(this, kt)
564        this => this%next
565      END DO
566      ibergs = icb_utl_count()
567      inbergs = ibergs
568      IF( lk_mpp )   CALL mpp_sum(inbergs)
569      IF( ibergs > 0 )   WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)')   &
570         &                                  cd_label, ibergs, inbergs, narea
571      !
572   END SUBROUTINE icb_utl_print
573
574
575   SUBROUTINE icb_utl_incr()
576      !!----------------------------------------------------------------------
577      !!                 ***  ROUTINE icb_utl_incr  ***
578      !!
579      !! ** Purpose :   
580      !!
581      ! Small routine for coping with very large integer values labelling icebergs
582      ! num_bergs is a array of integers
583      ! the first member is incremented in steps of jpnij starting from narea
584      ! this means each iceberg is labelled with a unique number
585      ! when this gets to the maximum allowed integer the second and subsequent members are
586      ! used to count how many times the member before cycles
587      !!----------------------------------------------------------------------
588      INTEGER ::   ii, ibig
589      !!----------------------------------------------------------------------
590
591      ibig = HUGE(num_bergs(1))
592      IF( ibig-jpnij < num_bergs(1) ) THEN
593         num_bergs(1) = narea
594         DO ii = 2,nkounts
595            IF( num_bergs(ii) == ibig ) THEN
596               num_bergs(ii) = 0
597               IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
598            ELSE
599               num_bergs(ii) = num_bergs(ii) + 1
600               EXIT
601            ENDIF
602         END DO
603      ELSE
604         num_bergs(1) = num_bergs(1) + jpnij
605      ENDIF
606      !
607   END SUBROUTINE icb_utl_incr
608
609
610   INTEGER FUNCTION icb_utl_count()
611      !!----------------------------------------------------------------------
612      !!                 ***  FUNCTION icb_utl_count  ***
613      !!
614      !! ** Purpose :   
615      !!----------------------------------------------------------------------
616      TYPE(iceberg), POINTER :: this
617      !!----------------------------------------------------------------------
618      !
619      icb_utl_count = 0
620      this => first_berg
621      DO WHILE( ASSOCIATED(this) )
622         icb_utl_count = icb_utl_count+1
623         this => this%next
624      END DO
625      !
626   END FUNCTION icb_utl_count
627
628
629   REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
630      !!----------------------------------------------------------------------
631      !!                 ***  FUNCTION icb_utl_mass  ***
632      !!
633      !! ** Purpose :   compute the mass all iceberg, all bergies or all bergs.
634      !!----------------------------------------------------------------------
635      TYPE(iceberg)      , POINTER  ::   first
636      TYPE(point)        , POINTER  ::   pt
637      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
638      !
639      TYPE(iceberg), POINTER ::   this
640      !!----------------------------------------------------------------------
641      icb_utl_mass = 0._wp
642      this => first
643      !
644      IF( PRESENT( justbergs  ) ) THEN
645         DO WHILE( ASSOCIATED( this ) )
646            pt => this%current_point
647            icb_utl_mass = icb_utl_mass + pt%mass         * this%mass_scaling
648            this => this%next
649         END DO
650      ELSEIF( PRESENT(justbits) ) THEN
651         DO WHILE( ASSOCIATED( this ) )
652            pt => this%current_point
653            icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
654            this => this%next
655         END DO
656      ELSE
657         DO WHILE( ASSOCIATED( this ) )
658            pt => this%current_point
659            icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
660            this => this%next
661         END DO
662      ENDIF
663      !
664   END FUNCTION icb_utl_mass
665
666
667   REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
668      !!----------------------------------------------------------------------
669      !!                 ***  FUNCTION icb_utl_heat  ***
670      !!
671      !! ** Purpose :   compute the heat in all iceberg, all bergies or all bergs.
672      !!----------------------------------------------------------------------
673      TYPE(iceberg)      , POINTER  ::   first
674      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
675      !
676      TYPE(iceberg)      , POINTER  ::   this
677      TYPE(point)        , POINTER  ::   pt
678      !!----------------------------------------------------------------------
679      icb_utl_heat = 0._wp
680      this => first
681      !
682      IF( PRESENT( justbergs  ) ) THEN
683         DO WHILE( ASSOCIATED( this ) )
684            pt => this%current_point
685            icb_utl_heat = icb_utl_heat + pt%mass         * this%mass_scaling * pt%heat_density
686            this => this%next
687         END DO
688      ELSEIF( PRESENT(justbits) ) THEN
689         DO WHILE( ASSOCIATED( this ) )
690            pt => this%current_point
691            icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
692            this => this%next
693         END DO
694      ELSE
695         DO WHILE( ASSOCIATED( this ) )
696            pt => this%current_point
697            icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
698            this => this%next
699         END DO
700      ENDIF
701      !
702   END FUNCTION icb_utl_heat
703
704   !!======================================================================
705END MODULE icbutl
Note: See TracBrowser for help on using the repository browser.