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 @ 14898

Last change on this file since 14898 was 14898, checked in by dancopsey, 16 months ago

Make iceberg top melt a combination of sea ice top melt and ocean solar.

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