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

Last change on this file since 15372 was 15372, checked in by davestorkey, 12 months ago

trunk: fix for wind forcing of icebergs #2728

  • Property svn:keywords set to Id
File size: 43.4 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   !!           4.2   !  2020-07  (P. Mathiot)     simplification of interpolation routine
11   !!                 !                            and add Nacho Merino work
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   icb_utl_interp   :
16   !!   icb_utl_pos      : compute bottom left corner indice, weight and mask
17   !!   icb_utl_bilin_h  : interpolation field to icb position
18   !!   icb_utl_bilin_e  : interpolation of scale factor to icb position
19   !!----------------------------------------------------------------------
20   USE par_oce                             ! ocean parameters
21   USE oce,    ONLY: ts, uu, vv
22   USE dom_oce                             ! ocean domain
23   USE in_out_manager                      ! IO parameters
24   USE lbclnk                              ! lateral boundary condition
25   USE lib_mpp                             ! MPI code and lk_mpp in particular
26   USE icb_oce                             ! define iceberg arrays
27   USE sbc_oce                             ! ocean surface boundary conditions
28#if defined key_si3
29   USE ice,    ONLY: u_ice, v_ice, hm_i    ! SI3 variables
30   USE icevar                              ! ice_var_sshdyn
31   USE sbc_ice, ONLY: snwice_mass, snwice_mass_b
32#endif
33
34   IMPLICIT NONE
35   PRIVATE
36
37   INTERFACE icb_utl_bilin_h
38      MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h
39   END INTERFACE
40
41   PUBLIC   icb_utl_copy          ! routine called in icbstp module
42   PUBLIC   icb_utl_getkb         ! routine called in icbdyn and icbthm modules
43   PUBLIC   test_icb_utl_getkb    ! routine called in icbdyn and icbthm modules
44   PUBLIC   icb_utl_zavg          ! routine called in icbdyn and icbthm modules
45   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules
46   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module
47   PUBLIC   icb_utl_add           ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
48   PUBLIC   icb_utl_delete        ! routine called in icblbc, icbthm modules
49   PUBLIC   icb_utl_destroy       ! routine called in icbstp module
50   PUBLIC   icb_utl_track         ! routine not currently used, retain just in case
51   PUBLIC   icb_utl_print_berg    ! routine called in icbthm module
52   PUBLIC   icb_utl_print         ! routine called in icbini, icbstp module
53   PUBLIC   icb_utl_count         ! routine called in icbdia, icbini, icblbc, icbrst modules
54   PUBLIC   icb_utl_incr          ! routine called in icbini, icbclv modules
55   PUBLIC   icb_utl_yearday       ! routine called in icbclv, icbstp module
56   PUBLIC   icb_utl_mass          ! routine called in icbdia module
57   PUBLIC   icb_utl_heat          ! routine called in icbdia module
58
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
63   !! $Id$
64   !! Software governed by the CeCILL license (see ./LICENSE)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE icb_utl_copy( Kmm )
69      !!----------------------------------------------------------------------
70      !!                  ***  ROUTINE icb_utl_copy  ***
71      !!
72      !! ** Purpose :   iceberg initialization.
73      !!
74      !! ** Method  : - blah blah
75      !!----------------------------------------------------------------------
76      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp
77#if defined key_si3
78      REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m    !    ocean surface (ssh_m) if ice is not embedded
79      !                                              !    ocean surface in leads if ice is embedded   
80#endif
81      INTEGER :: jk   ! vertical loop index
82      INTEGER :: Kmm  ! ocean time levelindex
83      !
84      ! copy nemo forcing arrays into iceberg versions with extra halo
85      ! only necessary for variables not on T points
86      ! and ssh which is used to calculate gradients
87      !
88      ! surface forcing
89      !
90      ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
91      ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
92      sst_e(1:jpi,1:jpj) = sst_m(:,:)
93      sss_e(1:jpi,1:jpj) = sss_m(:,:)
94      fr_e (1:jpi,1:jpj) = fr_i (:,:)
95      ua_e (1:jpi,1:jpj) = utau_icb (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
96      va_e (1:jpi,1:jpj) = vtau_icb (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
97      ff_e(1:jpi,1:jpj) = ff_f (:,:) 
98      !
99      CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 )
100      CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -1._wp, 1, 1 )
101      CALL lbc_lnk_icb( 'icbutl', ua_e , 'U', -1._wp, 1, 1 )
102      CALL lbc_lnk_icb( 'icbutl', va_e , 'V', -1._wp, 1, 1 )
103#if defined key_si3
104      hi_e(1:jpi, 1:jpj) = hm_i (:,:) 
105      ui_e(1:jpi, 1:jpj) = u_ice(:,:)
106      vi_e(1:jpi, 1:jpj) = v_ice(:,:)
107      !     
108      ! compute ssh slope using ssh_lead if embedded
109      zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b)
110      ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)
111      !
112      CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 )
113      CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 )
114#else
115      ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)         
116#endif
117      !
118      ! (PM) could be improve with a 3d lbclnk gathering both variables
119      ! should be done once extra haloe generalised
120      IF ( ln_M2016 ) THEN
121         DO jk = 1,jpk
122            ! uoce
123            ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm)
124            CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 )
125            uoce_e(:,:,jk) = ztmp(:,:)
126            !
127            ! voce
128            ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm)
129            CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 )
130            voce_e(:,:,jk) = ztmp(:,:)
131            !
132            e3t_e(1:jpi,1:jpj,jk) = e3t(:,:,jk,Kmm)
133         END DO
134         toce_e(1:jpi,1:jpj,:) = ts(:,:,:,1,Kmm)
135      END IF
136      !
137   END SUBROUTINE icb_utl_copy
138
139
140   SUBROUTINE icb_utl_interp( pi, pj, pe1 , pssu, pui, pua, pssh_i,         &
141      &                               pe2 , pssv, pvi, pva, pssh_j,         &
142      &                               psst, psss, pcn, phi, pff   ,         &
143      &                               plon, plat, ptoce, puoce, pvoce, pe3t )
144      !!----------------------------------------------------------------------
145      !!                  ***  ROUTINE icb_utl_interp  ***
146      !!
147      !! ** Purpose :   interpolation
148      !!
149      !! ** Method  : - interpolate from various ocean arrays onto iceberg position
150      !!
151      !!       !!gm  CAUTION here I do not care of the slip/no-slip conditions
152      !!             this can be done later (not that easy to do...)
153      !!             right now, U is 0 in land so that the coastal value of velocity parallel to the coast
154      !!             is half the off shore value, wile the normal-to-the-coast value is zero.
155      !!             This is OK as a starting point.
156      !!       !!pm  HARD CODED: - rho_air now computed in sbcblk (what are the effect ?)
157      !!                         - drag coefficient (should it be namelist parameter ?)
158      !!
159      !!----------------------------------------------------------------------
160      REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential
161      REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors
162      REAL(wp), INTENT(  out), OPTIONAL ::   pssu, pssv, pui, pvi, pua, pva ! ocean, ice and wind speeds
163      REAL(wp), INTENT(  out), OPTIONAL ::   pssh_i, pssh_j                 ! ssh i- & j-gradients
164      REAL(wp), INTENT(  out), OPTIONAL ::   psst, psss, pcn, phi, pff      ! SST, SSS, ice concentration, ice thickness, Coriolis
165      REAL(wp), INTENT(  out), OPTIONAL ::   plat, plon                     ! position
166      REAL(wp), DIMENSION(jpk), INTENT(  out), OPTIONAL ::   ptoce, puoce, pvoce, pe3t   ! 3D variables
167      !
168      REAL(wp), DIMENSION(4) :: zwT  , zwU  , zwV  , zwF   ! interpolation weight
169      REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask
170      REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm
171      REAL(wp), DIMENSION(4,jpk) :: zw1d
172      INTEGER                :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner
173      INTEGER                :: iiTp, iiTm, ijTp, ijTm
174      REAL(wp) ::   zcd, zmod       ! local scalars
175      !!----------------------------------------------------------------------
176      !
177      ! get position, weight and mask
178      CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT )
179      CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU )
180      CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV )
181      CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF )
182      !
183      ! metrics and coordinates
184      IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors
185      IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
186      IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true.  )
187      IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. )
188      !
189      IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU        , .false. ) ! ocean velocities
190      IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV        , .false. ) !
191      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst
192      IF ( PRESENT(psss) ) psss = icb_utl_bilin_h( sss_e, iiT, ijT, zwT * zmskT, .false. ) ! sss
193      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration
194      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e , iiF, ijF, zwF        , .false. ) ! Coriolis parameter
195      !
196      IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN
197         pua  = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind
198         pva  = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! here (ua,va) are stress => rough conversion from stress to speed
199         zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient
200         zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  )
201         pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0
202         pva  = pva * zmod
203      END IF
204      !
205#if defined key_si3
206      IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU        , .false. ) ! sea-ice velocities
207      IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV        , .false. )
208      IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness
209#else
210      IF ( PRESENT(pui) ) pui = 0._wp
211      IF ( PRESENT(pvi) ) pvi = 0._wp
212      IF ( PRESENT(phi) ) phi = 0._wp
213#endif
214      !
215      ! Estimate SSH gradient in i- and j-direction (centred evaluation)
216      IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN
217         CALL icb_utl_pos( pi+0.1, pj    , 'T', iiTp, ijTp, zwTp, zmskTp )
218         CALL icb_utl_pos( pi-0.1, pj    , 'T', iiTm, ijTm, zwTm, zmskTm )
219         !
220         IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )
221         pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
222            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe1 )
223         !
224         CALL icb_utl_pos( pi    , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp )
225         CALL icb_utl_pos( pi    , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm )
226         !
227         IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
228         pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   &
229            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe2 )
230      END IF
231      !
232      ! 3d interpolation
233      IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN
234         ! no need to mask as 0 is a valid data for land
235         zw1d(1,:) = zwU(1) ; zw1d(2,:) = zwU(2) ; zw1d(3,:) = zwU(3) ; zw1d(4,:) = zwU(4) ;
236         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zw1d )
237
238         zw1d(1,:) = zwV(1) ; zw1d(2,:) = zwV(2) ; zw1d(3,:) = zwV(3) ; zw1d(4,:) = zwV(4) ;
239         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zw1d )
240      END IF
241
242      IF ( PRESENT(ptoce) ) THEN
243         ! for temperature we need to mask the weight properly
244         ! no need of extra halo as it is a T point variable
245         zw1d(1,:) = tmask(iiT  ,ijT  ,:) * zwT(1) * zmskT(1)
246         zw1d(2,:) = tmask(iiT+1,ijT  ,:) * zwT(2) * zmskT(2)
247         zw1d(3,:) = tmask(iiT  ,ijT+1,:) * zwT(3) * zmskT(3)
248         zw1d(4,:) = tmask(iiT+1,ijT+1,:) * zwT(4) * zmskT(4)
249         ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zw1d )
250      END IF
251      !
252      IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results
253      !
254   END SUBROUTINE icb_utl_interp
255
256   SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk )
257      !!----------------------------------------------------------------------
258      !!                  ***  FUNCTION icb_utl_bilin  ***
259      !!
260      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
261      !!                this version deals with extra halo points
262      !!
263      !!       !!gm  CAUTION an optional argument should be added to handle
264      !!             the slip/no-slip conditions  ==>>> to be done later
265      !!
266      !!----------------------------------------------------------------------
267      REAL(wp)              , INTENT(IN)  ::   pi, pj    ! targeted coordinates in (i,j) referential
268      CHARACTER(len=1)      , INTENT(IN)  ::   cd_type   ! point type
269      REAL(wp), DIMENSION(4), INTENT(OUT) ::   pw, pmsk  ! weight and mask
270      INTEGER ,               INTENT(OUT) ::   kii, kij  ! bottom left corner position in local domain
271      !
272      REAL(wp) :: zwi, zwj ! distance to bottom left corner
273      INTEGER  :: ierr 
274      !
275      !!----------------------------------------------------------------------
276      !
277      SELECT CASE ( cd_type )
278      CASE ( 'T' )
279         ! note that here there is no +0.5 added
280         ! since we're looking for four T points containing quadrant we're in of
281         ! current T cell
282         kii = MAX(0, INT( pi        ))
283         kij = MAX(0, INT( pj        ))    ! T-point
284         zwi = pi - REAL(kii,wp)
285         zwj = pj - REAL(kij,wp)
286      CASE ( 'U' )
287         kii = MAX(0, INT( pi-0.5_wp ))
288         kij = MAX(0, INT( pj        ))    ! U-point
289         zwi = pi - 0.5_wp - REAL(kii,wp)
290         zwj = pj - REAL(kij,wp)
291      CASE ( 'V' )
292         kii = MAX(0, INT( pi        ))
293         kij = MAX(0, INT( pj-0.5_wp ))    ! V-point
294         zwi = pi - REAL(kii,wp)
295         zwj = pj - 0.5_wp - REAL(kij,wp)
296      CASE ( 'F' )
297         kii = MAX(0, INT( pi-0.5_wp ))
298         kij = MAX(0, INT( pj-0.5_wp ))    ! F-point
299         zwi = pi - 0.5_wp - REAL(kii,wp)
300         zwj = pj - 0.5_wp - REAL(kij,wp)
301      END SELECT
302      kii = kii + (nn_hls-1)
303      kij = kij + (nn_hls-1)
304      !
305      ! compute weight
306      pw(1) = (1._wp-zwi) * (1._wp-zwj)
307      pw(2) =        zwi  * (1._wp-zwj)
308      pw(3) = (1._wp-zwi) *        zwj
309      pw(4) =        zwi  *        zwj
310      !
311      ! find position in this processor. Prevent near edge problems (see #1389)
312      !
313      IF (TRIM(cd_type) == 'T' ) THEN
314         ierr = 0
315         IF    ( kii <  mig( 1 ) ) THEN   ;  ierr = ierr + 1
316         ELSEIF( kii >= mig(jpi) ) THEN   ;  ierr = ierr + 1
317         ENDIF
318         !
319         IF    ( kij <  mjg( 1 ) ) THEN   ;   ierr = ierr + 1
320         ELSEIF( kij >= mjg(jpj) ) THEN   ;   ierr = ierr + 1
321         ENDIF
322         !
323         IF ( ierr > 0 ) THEN
324            WRITE(numicb,*) 'bottom left corner T point out of bound'
325            WRITE(numicb,*) pi, kii, mig( 1 ), mig(jpi)
326            WRITE(numicb,*) pj, kij, mjg( 1 ), mjg(jpj)
327            WRITE(numicb,*) pmsk
328            CALL FLUSH(numicb)
329            CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error).'       , &
330                 &                                'This can be fixed using rn_speed_limit=0.4 in &namberg.'                   , &
331                 &                                'More details in the corresponding iceberg.stat file (nn_verbose_level > 0).' )
332         END IF
333      END IF
334      !
335      ! find position in this processor. Prevent near edge problems (see #1389)
336      ! (PM) will be useless if extra halo is used in NEMO
337      !
338      IF    ( kii <= mig(1)-1 ) THEN   ;   kii = 0
339      ELSEIF( kii  > mig(jpi) ) THEN   ;   kii = jpi
340      ELSE                             ;   kii = mi1(kii)
341      ENDIF
342      IF    ( kij <= mjg(1)-1 ) THEN   ;   kij = 0
343      ELSEIF( kij  > mjg(jpj) ) THEN   ;   kij = jpj
344      ELSE                             ;   kij = mj1(kij)
345      ENDIF
346      !
347      ! define mask array
348      ! land value is not used in the interpolation
349      SELECT CASE ( cd_type )
350      CASE ( 'T' )
351         pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/)
352      CASE ( 'U' )
353         pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/)
354      CASE ( 'V' )
355         pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/)
356      CASE ( 'F' )
357         ! F case only used for coriolis, ff_f is not mask so zmask = 1
358         pmsk = 1.
359      END SELECT
360   END SUBROUTINE icb_utl_pos
361
362   REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon )
363      !!----------------------------------------------------------------------
364      !!                  ***  FUNCTION icb_utl_bilin  ***
365      !!
366      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
367      !!                this version deals with extra halo points
368      !!
369      !!       !!gm  CAUTION an optional argument should be added to handle
370      !!             the slip/no-slip conditions  ==>>> to be done later
371      !!
372      !!----------------------------------------------------------------------
373      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated
374      REAL(wp), DIMENSION(4)              , INTENT(in) ::   pw        ! weight
375      LOGICAL                             , INTENT(in) ::   pllon     ! input data is a longitude
376      INTEGER ,                             INTENT(in) ::   pii, pij  ! bottom left corner
377      !
378      REAL(wp), DIMENSION(4) :: zdat ! input data
379      !!----------------------------------------------------------------------
380      !
381      ! data
382      zdat(1) = pfld(pii  ,pij  )
383      zdat(2) = pfld(pii+1,pij  )
384      zdat(3) = pfld(pii  ,pij+1)
385      zdat(4) = pfld(pii+1,pij+1)
386      !
387      IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN
388         WHERE( zdat < 0._wp ) zdat = zdat + 360._wp
389      ENDIF
390      !
391      ! compute interpolated value
392      icb_utl_bilin_2d_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4)) 
393      !
394      IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp
395      !
396   END FUNCTION icb_utl_bilin_2d_h
397
398   FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw )
399      !!----------------------------------------------------------------------
400      !!                  ***  FUNCTION icb_utl_bilin  ***
401      !!
402      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
403      !!                this version deals with extra halo points
404      !!
405      !!       !!gm  CAUTION an optional argument should be added to handle
406      !!             the slip/no-slip conditions  ==>>> to be done later
407      !!
408      !!----------------------------------------------------------------------
409      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated
410      REAL(wp), DIMENSION(4,jpk)               , INTENT(in) ::   pw        ! weight
411      INTEGER ,                                  INTENT(in) ::   pii, pij  ! bottom left corner
412      REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h
413      !
414      REAL(wp), DIMENSION(4,jpk) :: zdat ! input data
415      INTEGER :: jk
416      !!----------------------------------------------------------------------
417      !
418      ! data
419      zdat(1,:) = pfld(pii  ,pij  ,:)
420      zdat(2,:) = pfld(pii+1,pij  ,:)
421      zdat(3,:) = pfld(pii  ,pij+1,:)
422      zdat(4,:) = pfld(pii+1,pij+1,:)
423      !
424      ! compute interpolated value
425      DO jk=1,jpk
426         icb_utl_bilin_3d_h(jk) =   ( zdat(1,jk)*pw(1,jk) + zdat(2,jk)*pw(2,jk) + zdat(3,jk)*pw(3,jk) + zdat(4,jk)*pw(4,jk) ) &
427            &                     /   MAX(1.e-20, pw(1,jk)+pw(2,jk)+pw(3,jk)+pw(4,jk)) 
428      END DO
429      !
430   END FUNCTION icb_utl_bilin_3d_h
431
432   REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
433      !!----------------------------------------------------------------------
434      !!                  ***  FUNCTION dom_init  ***
435      !!
436      !! ** Purpose :   bilinear interpolation at berg location of horizontal scale factor
437      !! ** Method  :   interpolation done using the 4 nearest grid points among
438      !!                t-, u-, v-, and f-points.
439      !!----------------------------------------------------------------------
440      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
441      REAL(wp)                , INTENT(IN) ::   pi , pj              ! iceberg position
442      !
443      ! weights corresponding to corner points of a T cell quadrant
444      REAL(wp) ::   zi, zj          ! local real
445      INTEGER  ::   ii, ij          ! bottom left corner coordinate in local domain
446      !
447      ! values at corner points of a T cell quadrant
448      ! 00 = bottom left, 10 = bottom right, 01 = top left, 11 = top right
449      REAL(wp) ::   ze00, ze10, ze01, ze11
450      !!----------------------------------------------------------------------
451      !
452      ! cannot used iiT because need ii/ij reltaive to global indices not local one
453      ii = MAX(1, INT( pi ))   ;   ij = MAX(1, INT( pj ))            ! left bottom T-point (i,j) indices
454      !
455      ! fractional box spacing
456      ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell
457      ! 0.5 <= zi < 1  , 0   <= zj < 0.5   -->  NE quadrant
458      ! 0   <= zi < 0.5, 0.5 <= zj < 1     -->  SE quadrant
459      ! 0.5 <= zi < 1  , 0.5 <= zj < 1     -->  SW quadrant
460
461      zi = pi - REAL(ii,wp)          !!gm use here mig, mjg arrays
462      zj = pj - REAL(ij,wp)
463
464      ! conversion to local domain (no need to do a sanity check already done in icbpos)
465      ii = mi1(ii) + (nn_hls-1)
466      ij = mj1(ij) + (nn_hls-1)
467      !
468      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
469         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NE quadrant
470            !                                                      !             i=I       i=I+1/2
471            ze01 = pev(ii  ,ij  )   ;   ze11 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
472            ze00 = pet(ii  ,ij  )   ;   ze10 = peu(ii  ,ij  )      !   j=J        T ------- U
473            zi = 2._wp * zi
474            zj = 2._wp * zj
475         ELSE                                                      !  SE quadrant
476            !                                                                    !             i=I       i=I+1/2
477            ze01 = pet(ii  ,ij+1)   ;   ze11 = peu(ii  ,ij+1)      !   j=J+1      T ------- U
478            ze00 = pev(ii  ,ij  )   ;   ze10 = pef(ii  ,ij  )      !   j=J+1/2    V ------- F
479            zi = 2._wp *  zi
480            zj = 2._wp * (zj-0.5_wp)
481         ENDIF
482      ELSE
483         IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NW quadrant
484            !                                                                    !             i=I       i=I+1/2
485            ze01 = pef(ii  ,ij  )   ;   ze11 = pev(ii+1,ij)        !   j=J+1/2    F ------- V
486            ze00 = peu(ii  ,ij  )   ;   ze10 = pet(ii+1,ij)        !   j=J        U ------- T
487            zi = 2._wp * (zi-0.5_wp)
488            zj = 2._wp *  zj
489         ELSE                                                      !  SW quadrant
490            !                                                                    !             i=I+1/2   i=I+1
491            ze01 = peu(ii  ,ij+1)   ;   ze11 = pet(ii+1,ij+1)      !   j=J+1      U ------- T
492            ze00 = pef(ii  ,ij  )   ;   ze10 = pev(ii+1,ij  )      !   j=J+1/2    F ------- V
493            zi = 2._wp * (zi-0.5_wp)
494            zj = 2._wp * (zj-0.5_wp)
495         ENDIF
496      ENDIF
497      !
498      icb_utl_bilin_e = ( ze01 * (1._wp-zi) + ze11 * zi ) *        zj    &
499         &            + ( ze00 * (1._wp-zi) + ze10 * zi ) * (1._wp-zj)
500      !
501   END FUNCTION icb_utl_bilin_e
502
503   SUBROUTINE icb_utl_getkb( kb, pe3, pD )
504      !!----------------------------------------------------------------------
505      !!                ***  ROUTINE icb_utl_getkb         ***
506      !!
507      !! ** Purpose :   compute the latest level affected by icb
508      !!
509      !!----------------------------------------------------------------------
510      INTEGER,                INTENT(out):: kb
511      REAL(wp), DIMENSION(:), INTENT(in) :: pe3
512      REAL(wp),               INTENT(in) :: pD
513      !!
514      INTEGER  :: jk
515      REAL(wp) :: zdepw
516      !!----------------------------------------------------------------------
517      !!
518      zdepw = pe3(1) ; kb = 2
519      DO WHILE ( zdepw <  pD)
520         zdepw = zdepw + pe3(kb)
521         kb = kb + 1
522      END DO
523      kb = MIN(kb - 1,jpk)
524   END SUBROUTINE
525
526   SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb )
527      !!----------------------------------------------------------------------
528      !!                ***  ROUTINE icb_utl_getkb         ***
529      !!
530      !! ** Purpose :   compute the vertical average of ocean properties affected by icb
531      !!
532      !!----------------------------------------------------------------------
533      INTEGER,                INTENT(in ) :: kb        ! deepest level affected by icb
534      REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile
535      REAL(wp),               INTENT(in ) :: pD        ! draft
536      REAL(wp),               INTENT(out) :: pzavg     ! z average
537      !!----------------------------------------------------------------------
538      INTEGER  :: jk
539      REAL(wp) :: zdep
540      !!----------------------------------------------------------------------
541      pzavg = 0.0 ; zdep = 0.0
542      DO jk = 1,kb-1
543         pzavg = pzavg + pe3(jk)*pdat(jk)
544         zdep  = zdep  + pe3(jk)
545      END DO
546      ! if kb is limited by mbkt  => bottom value is used between bathy and icb tail
547      ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v)
548      pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD
549   END SUBROUTINE
550
551   SUBROUTINE icb_utl_add( bergvals, ptvals )
552      !!----------------------------------------------------------------------
553      !!                ***  ROUTINE icb_utl_add           ***
554      !!
555      !! ** Purpose :   add a new berg to the iceberg list
556      !!
557      !!----------------------------------------------------------------------
558      TYPE(iceberg), INTENT(in)           ::   bergvals
559      TYPE(point)  , INTENT(in)           ::   ptvals
560      !
561      TYPE(iceberg), POINTER ::   new => NULL()
562      !!----------------------------------------------------------------------
563      !
564      new => NULL()
565      CALL icb_utl_create( new, bergvals, ptvals )
566      CALL icb_utl_insert( new )
567      new => NULL()     ! Clear new
568      !
569   END SUBROUTINE icb_utl_add         
570
571
572   SUBROUTINE icb_utl_create( berg, bergvals, ptvals )
573      !!----------------------------------------------------------------------
574      !!                ***  ROUTINE icb_utl_create  ***
575      !!
576      !! ** Purpose :   add a new berg to the iceberg list
577      !!
578      !!----------------------------------------------------------------------
579      TYPE(iceberg), INTENT(in) ::   bergvals
580      TYPE(point)  , INTENT(in) ::   ptvals
581      TYPE(iceberg), POINTER    ::   berg
582      !
583      TYPE(point)  , POINTER    ::   pt
584      INTEGER                   ::   istat
585      !!----------------------------------------------------------------------
586      !
587      IF( ASSOCIATED(berg) )   CALL ctl_stop( 'icebergs, icb_utl_create: berg already associated' )
588      ALLOCATE(berg, STAT=istat)
589      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
590      berg%number(:) = bergvals%number(:)
591      berg%mass_scaling = bergvals%mass_scaling
592      berg%prev => NULL()
593      berg%next => NULL()
594      !
595      ALLOCATE(pt, STAT=istat)
596      IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
597      pt = ptvals
598      berg%current_point => pt
599      !
600   END SUBROUTINE icb_utl_create
601
602
603   SUBROUTINE icb_utl_insert( newberg )
604      !!----------------------------------------------------------------------
605      !!                 ***  ROUTINE icb_utl_insert  ***
606      !!
607      !! ** Purpose :   add a new berg to the iceberg list
608      !!
609      !!----------------------------------------------------------------------
610      TYPE(iceberg), POINTER  ::   newberg
611      !
612      TYPE(iceberg), POINTER  ::   this, prev, last
613      !!----------------------------------------------------------------------
614      !
615      IF( ASSOCIATED( first_berg ) ) THEN
616         last => first_berg
617         DO WHILE (ASSOCIATED(last%next))
618            last => last%next
619         ENDDO
620         newberg%prev => last
621         last%next    => newberg
622         last         => newberg
623      ELSE                       ! list is empty so create it
624         first_berg => newberg
625      ENDIF
626      !
627   END SUBROUTINE icb_utl_insert
628
629
630   REAL(wp) FUNCTION icb_utl_yearday(kmon, kday, khr, kmin, ksec)
631      !!----------------------------------------------------------------------
632      !!                 ***  FUNCTION icb_utl_yearday  ***
633      !!
634      !! ** Purpose :   
635      !!
636      ! sga - improved but still only applies to 365 day year, need to do this properly
637      !
638      !!gm  all these info are already known in daymod, no???
639      !!
640      !!----------------------------------------------------------------------
641      INTEGER, INTENT(in)     :: kmon, kday, khr, kmin, ksec
642      !
643      INTEGER, DIMENSION(12)  :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
644      !!----------------------------------------------------------------------
645      !
646      icb_utl_yearday = REAL( SUM( imonths(1:kmon) ), wp )
647      icb_utl_yearday = icb_utl_yearday + REAL(kday-1,wp) + (REAL(khr,wp) + (REAL(kmin,wp) + REAL(ksec,wp)/60.)/60.)/24.
648      !
649   END FUNCTION icb_utl_yearday
650
651   !!-------------------------------------------------------------------------
652
653   SUBROUTINE icb_utl_delete( first, berg )
654      !!----------------------------------------------------------------------
655      !!                 ***  ROUTINE icb_utl_delete  ***
656      !!
657      !! ** Purpose :   
658      !!
659      !!----------------------------------------------------------------------
660      TYPE(iceberg), POINTER :: first, berg
661      !!----------------------------------------------------------------------
662      ! Connect neighbors to each other
663      IF ( ASSOCIATED(berg%prev) ) THEN
664        berg%prev%next => berg%next
665      ELSE
666        first => berg%next
667      ENDIF
668      IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
669      !
670      CALL icb_utl_destroy(berg)
671      !
672   END SUBROUTINE icb_utl_delete
673
674
675   SUBROUTINE icb_utl_destroy( berg )
676      !!----------------------------------------------------------------------
677      !!                 ***  ROUTINE icb_utl_destroy  ***
678      !!
679      !! ** Purpose :   remove a single iceberg instance
680      !!
681      !!----------------------------------------------------------------------
682      TYPE(iceberg), POINTER :: berg
683      !!----------------------------------------------------------------------
684      !
685      ! Remove any points
686      IF( ASSOCIATED( berg%current_point ) )   DEALLOCATE( berg%current_point )
687      !
688      DEALLOCATE(berg)
689      !
690   END SUBROUTINE icb_utl_destroy
691
692
693   SUBROUTINE icb_utl_track( knum, cd_label, kt )
694      !!----------------------------------------------------------------------
695      !!                 ***  ROUTINE icb_utl_track  ***
696      !!
697      !! ** Purpose :   
698      !!
699      !!----------------------------------------------------------------------
700      INTEGER, DIMENSION(nkounts)    :: knum       ! iceberg number
701      CHARACTER(len=*)               :: cd_label   !
702      INTEGER                        :: kt         ! timestep number
703      !
704      TYPE(iceberg), POINTER         :: this
705      LOGICAL                        :: match
706      INTEGER                        :: k
707      !!----------------------------------------------------------------------
708      !
709      this => first_berg
710      DO WHILE( ASSOCIATED(this) )
711         match = .TRUE.
712         DO k = 1, nkounts
713            IF( this%number(k) /= knum(k) ) match = .FALSE.
714         END DO
715         IF( match )   CALL icb_utl_print_berg(this, kt)
716         this => this%next
717      END DO
718      !
719   END SUBROUTINE icb_utl_track
720
721
722   SUBROUTINE icb_utl_print_berg( berg, kt )
723      !!----------------------------------------------------------------------
724      !!                 ***  ROUTINE icb_utl_print_berg  ***
725      !!
726      !! ** Purpose :   print one
727      !!
728      !!----------------------------------------------------------------------
729      TYPE(iceberg), POINTER :: berg
730      TYPE(point)  , POINTER :: pt
731      INTEGER                :: kt      ! timestep number
732      !!----------------------------------------------------------------------
733      !
734      IF (nn_verbose_level == 0) RETURN
735      pt => berg%current_point
736      WRITE(numicb, 9200) kt, berg%number(1), &
737                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
738                   pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi
739      CALL flush( numicb )
740 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
741      !
742   END SUBROUTINE icb_utl_print_berg
743
744
745   SUBROUTINE icb_utl_print( cd_label, kt )
746      !!----------------------------------------------------------------------
747      !!                 ***  ROUTINE icb_utl_print  ***
748      !!
749      !! ** Purpose :   print many
750      !!
751      !!----------------------------------------------------------------------
752      CHARACTER(len=*)       :: cd_label
753      INTEGER                :: kt             ! timestep number
754      !
755      INTEGER                :: ibergs, inbergs
756      TYPE(iceberg), POINTER :: this
757      !!----------------------------------------------------------------------
758      !
759      IF (nn_verbose_level == 0) RETURN
760      this => first_berg
761      IF( ASSOCIATED(this) ) THEN
762         WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
763         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
764            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi'
765      ENDIF
766      DO WHILE( ASSOCIATED(this) )
767        CALL icb_utl_print_berg(this, kt)
768        this => this%next
769      END DO
770      ibergs = icb_utl_count()
771      inbergs = ibergs
772      CALL mpp_sum('icbutl', inbergs)
773      IF( ibergs > 0 )   WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)')   &
774         &                                  cd_label, ibergs, inbergs, narea
775      !
776   END SUBROUTINE icb_utl_print
777
778
779   SUBROUTINE icb_utl_incr()
780      !!----------------------------------------------------------------------
781      !!                 ***  ROUTINE icb_utl_incr  ***
782      !!
783      !! ** Purpose :   
784      !!
785      ! Small routine for coping with very large integer values labelling icebergs
786      ! num_bergs is a array of integers
787      ! the first member is incremented in steps of jpnij starting from narea
788      ! this means each iceberg is labelled with a unique number
789      ! when this gets to the maximum allowed integer the second and subsequent members are
790      ! used to count how many times the member before cycles
791      !!----------------------------------------------------------------------
792      INTEGER ::   ii, ibig
793      !!----------------------------------------------------------------------
794
795      ibig = HUGE(num_bergs(1))
796      IF( ibig-jpnij < num_bergs(1) ) THEN
797         num_bergs(1) = narea
798         DO ii = 2,nkounts
799            IF( num_bergs(ii) == ibig ) THEN
800               num_bergs(ii) = 0
801               IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
802            ELSE
803               num_bergs(ii) = num_bergs(ii) + 1
804               EXIT
805            ENDIF
806         END DO
807      ELSE
808         num_bergs(1) = num_bergs(1) + jpnij
809      ENDIF
810      !
811   END SUBROUTINE icb_utl_incr
812
813
814   INTEGER FUNCTION icb_utl_count()
815      !!----------------------------------------------------------------------
816      !!                 ***  FUNCTION icb_utl_count  ***
817      !!
818      !! ** Purpose :   
819      !!----------------------------------------------------------------------
820      TYPE(iceberg), POINTER :: this
821      !!----------------------------------------------------------------------
822      !
823      icb_utl_count = 0
824      this => first_berg
825      DO WHILE( ASSOCIATED(this) )
826         icb_utl_count = icb_utl_count+1
827         this => this%next
828      END DO
829      !
830   END FUNCTION icb_utl_count
831
832
833   REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
834      !!----------------------------------------------------------------------
835      !!                 ***  FUNCTION icb_utl_mass  ***
836      !!
837      !! ** Purpose :   compute the mass all iceberg, all berg bits or all bergs.
838      !!----------------------------------------------------------------------
839      TYPE(iceberg)      , POINTER  ::   first
840      TYPE(point)        , POINTER  ::   pt
841      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
842      !
843      TYPE(iceberg), POINTER ::   this
844      !!----------------------------------------------------------------------
845      icb_utl_mass = 0._wp
846      this => first
847      !
848      IF( PRESENT( justbergs  ) ) THEN
849         DO WHILE( ASSOCIATED( this ) )
850            pt => this%current_point
851            icb_utl_mass = icb_utl_mass + pt%mass         * this%mass_scaling
852            this => this%next
853         END DO
854      ELSEIF( PRESENT(justbits) ) THEN
855         DO WHILE( ASSOCIATED( this ) )
856            pt => this%current_point
857            icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
858            this => this%next
859         END DO
860      ELSE
861         DO WHILE( ASSOCIATED( this ) )
862            pt => this%current_point
863            icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
864            this => this%next
865         END DO
866      ENDIF
867      !
868   END FUNCTION icb_utl_mass
869
870
871   REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
872      !!----------------------------------------------------------------------
873      !!                 ***  FUNCTION icb_utl_heat  ***
874      !!
875      !! ** Purpose :   compute the heat in all iceberg, all bergies or all bergs.
876      !!----------------------------------------------------------------------
877      TYPE(iceberg)      , POINTER  ::   first
878      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
879      !
880      TYPE(iceberg)      , POINTER  ::   this
881      TYPE(point)        , POINTER  ::   pt
882      !!----------------------------------------------------------------------
883      icb_utl_heat = 0._wp
884      this => first
885      !
886      IF( PRESENT( justbergs  ) ) THEN
887         DO WHILE( ASSOCIATED( this ) )
888            pt => this%current_point
889            icb_utl_heat = icb_utl_heat + pt%mass         * this%mass_scaling * pt%heat_density
890            this => this%next
891         END DO
892      ELSEIF( PRESENT(justbits) ) THEN
893         DO WHILE( ASSOCIATED( this ) )
894            pt => this%current_point
895            icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
896            this => this%next
897         END DO
898      ELSE
899         DO WHILE( ASSOCIATED( this ) )
900            pt => this%current_point
901            icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
902            this => this%next
903         END DO
904      ENDIF
905      !
906   END FUNCTION icb_utl_heat
907
908   SUBROUTINE test_icb_utl_getkb
909      !!----------------------------------------------------------------------
910      !!                 ***  FUNCTION test_icb_utl_getkb  ***
911      !!
912      !! ** Purpose : Test routine icb_utl_getkb, icb_utl_zavg
913      !! ** Methode : Call each subroutine with specific input data
914      !!              What should be the output is easy to determined and check
915      !!              if NEMO return the correct answer.
916      !! ** Comments : not called, if needed a CALL test_icb_utl_getkb need to be added in icb_step
917      !!----------------------------------------------------------------------
918      INTEGER :: ikb
919      REAL(wp) :: zD, zout
920      REAL(wp), DIMENSION(jpk) :: ze3, zin
921      WRITE(numout,*) 'Test icb_utl_getkb : '
922      zD = 0.0 ; ze3= 20.0
923      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
924      CALL icb_utl_getkb(ikb, ze3, zD)
925      WRITE(numout,*) 'OUTPUT : kb = ',ikb
926
927      zD = 8000000.0 ; ze3= 20.0
928      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
929      CALL icb_utl_getkb(ikb, ze3, zD)
930      WRITE(numout,*) 'OUTPUT : kb = ',ikb
931
932      zD = 80.0 ; ze3= 20.0
933      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
934      CALL icb_utl_getkb(ikb, ze3, zD)
935      WRITE(numout,*) 'OUTPUT : kb = ',ikb
936
937      zD = 85.0 ; ze3= 20.0
938      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
939      CALL icb_utl_getkb(ikb, ze3, zD)
940      WRITE(numout,*) 'OUTPUT : kb = ',ikb
941
942      zD = 75.0 ; ze3= 20.0
943      WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1)
944      CALL icb_utl_getkb(ikb, ze3, zD)
945      WRITE(numout,*) 'OUTPUT : kb = ',ikb
946
947      WRITE(numout,*) '=================================='
948      WRITE(numout,*) 'Test icb_utl_zavg'
949      zD = 0.0 ; ze3= 20.0 ; zin=1.0
950      CALL icb_utl_getkb(ikb, ze3, zD)
951      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
952      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
953      WRITE(numout,*) 'OUTPUT : zout = ',zout
954
955      zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
956      CALL icb_utl_getkb(ikb, ze3, zD)
957      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
958      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
959      WRITE(numout,*) 'OUTPUT : zout = ',zout
960      CALL FLUSH(numout)
961
962      zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0
963      CALL icb_utl_getkb(ikb, ze3, zD)
964      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
965      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
966      WRITE(numout,*) 'OUTPUT : zout = ',zout
967
968      zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0
969      CALL icb_utl_getkb(ikb, ze3, zD)
970      ikb = 2
971      CALL icb_utl_zavg(zout, zin, ze3, zD, ikb)
972      WRITE(numout,*) 'INPUT  : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb
973      WRITE(numout,*) 'OUTPUT : zout = ',zout
974
975      CALL FLUSH(numout)
976
977   END SUBROUTINE test_icb_utl_getkb
978
979   !!======================================================================
980END MODULE icbutl
Note: See TracBrowser for help on using the repository browser.