source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 6486

Last change on this file since 6486 was 6486, checked in by davestorkey, 5 years ago

Remove SVN keywords from UKMO/dev_r5518_GO6_package branch.

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