source: NEMO/branches/2019/fix_ticket2238_solution1/src/OCE/ICB/icbutl.F90 @ 10680

Last change on this file since 10680 was 10680, checked in by mathiot, 20 months ago

cosmetic change (ticket #2238)

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