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

source: NEMO/branches/UKMO/NEMO_4.0_GO8_package_text_diagnostics/src/OCE/ICB/icbutl.F90 @ 10948

Last change on this file since 10948 was 10948, checked in by andmirek, 5 years ago

GMED 462 iceberg model

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