source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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