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

source: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90 @ 13359

Last change on this file since 13359 was 13359, checked in by mathiot, 4 years ago

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