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

source: NEMO/trunk/src/OCE/ICB/icbutl.F90 @ 13286

Last change on this file since 13286 was 13281, checked in by ayoung, 4 years ago

Merging changes for ticket #2407. Patch for iceberg melting stability. Fully SETTE tested: this merge introduces a change to the results in ORCA2_ICE_PISCES.

  • Property svn:keywords set to Id
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.