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/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/ICB – NEMO

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/ICB/icbutl.F90 @ 10126

Last change on this file since 10126 was 10126, checked in by jchanut, 6 years ago

Merge with trunk

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