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 NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB – NEMO

source: NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90 @ 10679

Last change on this file since 10679 was 10679, checked in by mathiot, 5 years ago

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