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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 3821

Last change on this file since 3821 was 3821, checked in by acc, 11 years ago

Branch dev_MERGE_2012. #1060. Minor alterations to icbutl.F90 and icbclv.F90 (ICeBerg trajectory component)

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