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/UKMO/NEMO_4.0.4_bouncing_icebergs/src/OCE/ICB – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_bouncing_icebergs/src/OCE/ICB/icbutl.F90 @ 15421

Last change on this file since 15421 was 15421, checked in by dancopsey, 11 months ago

Reverted all my efforts to generate top melt out of this branch. This branch is going back to simply do the bouncing icebergs thing which is what it was designed for.

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