ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/FLO/flo4rk.F90

MODULE flo4rk
   !!======================================================================
   !!                    ***  MODULE  flo4rk  ***
   !! Ocean floats :   trajectory computation using a 4th order Runge-Kutta
   !!======================================================================
#if   defined key_floats   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_floats'                                     float trajectories
   !!----------------------------------------------------------------------
   !!   flo_4rk        : Compute the geographical position of floats
   !!   flo_interp     : interpolation
   !!----------------------------------------------------------------------
   USE flo_oce         ! ocean drifting floats
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE in_out_manager  ! I/O manager
   USE wrk_nemo        ! working array

   IMPLICIT NONE
   PRIVATE

   PUBLIC   flo_4rk    ! routine called by floats.F90

   !                                   ! RK4 and Lagrange interpolation coefficients
   REAL(wp), DIMENSION (4) ::   tcoef1 = (/  1.0  ,  0.5  ,  0.5  ,  0.0  /)   ! 
   REAL(wp), DIMENSION (4) ::   tcoef2 = (/  0.0  ,  0.5  ,  0.5  ,  1.0  /)   !
   REAL(wp), DIMENSION (4) ::   scoef2 = (/  1.0  ,  2.0  ,  2.0  ,  1.0  /)   !
   REAL(wp), DIMENSION (4) ::   rcoef  = (/-1./6. , 1./2. ,-1./2. , 1./6. /)   !
   REAL(wp), DIMENSION (3) ::   scoef1 = (/  0.5  ,  0.5  ,  1.0  /)           !

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: flo4rk.F90 3294 2012-01-28 16:44:18Z rblod $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE flo_4rk( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE flo_4rk  ***
      !!
      !!  ** Purpose :   Compute the geographical position (lat,lon,depth)
      !!       of each float at each time step.
      !! 
      !!  ** Method  :   The position of a float is computed with a 4th order
      !!       Runge-Kutta scheme and and Lagrange interpolation.
      !!         We need to know the velocity field, the old positions of the
      !!       floats and the grid defined on the domain.
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      !!
      INTEGER ::  jfl, jind           ! dummy loop indices
      INTEGER ::  ierror              ! error value

      REAL(wp), POINTER, DIMENSION(:)   ::   zgifl , zgjfl , zgkfl    ! index RK  positions
      REAL(wp), POINTER, DIMENSION(:)   ::   zufl  , zvfl  , zwfl     ! interpolated velocity at the float position 
      REAL(wp), POINTER, DIMENSION(:,:) ::   zrkxfl, zrkyfl, zrkzfl   ! RK coefficients
      !!---------------------------------------------------------------------
      CALL wrk_alloc( jpnfl,    zgifl , zgjfl , zgkfl  , zufl, zvfl, zwfl)
      CALL wrk_alloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
      !
      IF( ierror /= 0 ) THEN
         WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
      ENDIF

    
      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'flo_4rk : compute Runge Kutta trajectories for floats '
         IF(lwp) WRITE(numout,*) '~~~~~~~'
      ENDIF

      ! Verification of the floats positions. If one of them leave the domain
      ! domain we replace the float near the border.
      DO jfl = 1, jpnfl
         ! i-direction
         IF( tpifl(jfl) <= 1.5 ) THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the WEST border.'
            tpifl(jfl) = tpifl(jfl) + 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=',tpifl(jfl)
         ENDIF
          
         IF( tpifl(jfl) >= jpi-.5 ) THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the EAST border.'
            tpifl(jfl) = tpifl(jfl) - 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=', tpifl(jfl)
         ENDIF
         ! j-direction
         IF( tpjfl(jfl) <= 1.5 ) THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the SOUTH border.'
            tpjfl(jfl) = tpjfl(jfl) + 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
         ENDIF
           
         IF( tpjfl(jfl) >= jpj-.5 ) THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the NORTH border.'
            tpjfl(jfl) = tpjfl(jfl) - 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
         ENDIF
         ! k-direction
         IF( tpkfl(jfl) <= .5 ) THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the TOP border.'
            tpkfl(jfl) = tpkfl(jfl) + 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
         ENDIF
         
         IF( tpkfl(jfl) >= jpk-.5 )  THEN
            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the BOTTOM border.'
            tpkfl(jfl) = tpkfl(jfl) - 1.
            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
         ENDIF
      END DO
      
      ! 4 steps of Runge-Kutta algorithme
      ! initialisation of the positions 
      
      DO jfl = 1, jpnfl
         zgifl(jfl) = tpifl(jfl)
         zgjfl(jfl) = tpjfl(jfl)
         zgkfl(jfl) = tpkfl(jfl)
      END DO
       
      DO  jind = 1, 4         
      
         ! for each step we compute the compute the velocity with Lagrange interpolation
         CALL flo_interp( zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
         
         ! computation of Runge-Kutta factor
         DO jfl = 1, jpnfl
            zrkxfl(jfl,jind) = rdt*zufl(jfl)
            zrkyfl(jfl,jind) = rdt*zvfl(jfl)
            zrkzfl(jfl,jind) = rdt*zwfl(jfl)
         END DO
         IF( jind /= 4 ) THEN
            DO jfl = 1, jpnfl
               zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
               zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
               zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
            END DO
         ENDIF
      END DO
      DO jind = 1, 4
         DO jfl = 1, jpnfl
            tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
            tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
            tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
         END DO
      END DO
      !
      CALL wrk_dealloc( jpnfl,    zgifl , zgjfl , zgkfl  , zufl, zvfl, zwfl)
      CALL wrk_dealloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
      !
   END SUBROUTINE flo_4rk


   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
      &                   pufl, pvfl, pwfl, ki )
      !!----------------------------------------------------------------------
      !!                ***  ROUTINE flointerp  ***
      !!
      !! ** Purpose :   Interpolation of the velocity on the float position
      !! 
      !! ** Method  :   Lagrange interpolation with the 64 neighboring
      !!      points. This routine is call 4 time at each time step to
      !!      compute velocity at the date and the position we need to
      !!      integrated with RK method.
      !!----------------------------------------------------------------------
      REAL(wp) , DIMENSION(jpnfl), INTENT(in   ) ::   pxt , pyt , pzt    ! position of the float
      REAL(wp) , DIMENSION(jpnfl), INTENT(  out) ::   pufl, pvfl, pwfl   ! velocity at this position
      INTEGER                    , INTENT(in   ) ::   ki                 !
      !!
      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
      INTEGER  , POINTER, DIMENSION(:)       ::   iilu, ijlu, iklu   ! nearest neighbour INDEX-u
      INTEGER  , POINTER, DIMENSION(:)       ::   iilv, ijlv, iklv   ! nearest neighbour INDEX-v
      INTEGER  , POINTER, DIMENSION(:)       ::   iilw, ijlw, iklw   ! nearest neighbour INDEX-w
      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidu, ijdu, ikdu   ! 64 nearest neighbour INDEX-u
      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidv, ijdv, ikdv   ! 64 nearest neighbour INDEX-v
      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidw, ijdw, ikdw   ! 64 nearest neighbour INDEX-w
      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxu, zlagyu, zlagzu   ! Lagrange  coefficients
      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxv, zlagyv, zlagzv   !    -           -
      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxw, zlagyw, zlagzw   !    -           -
      REAL(wp) , POINTER, DIMENSION(:,:,:,:) ::   ztufl , ztvfl , ztwfl   ! velocity at choosen time step
      !!---------------------------------------------------------------------
      CALL wrk_alloc( jpnfl,    iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
      CALL wrk_alloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
      CALL wrk_alloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
      CALL wrk_alloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
 
      ! Interpolation of U velocity

      ! nearest neightboring point for computation of u       
      DO jfl = 1, jpnfl
         iilu(jfl) = INT(pxt(jfl)-.5)
         ijlu(jfl) = INT(pyt(jfl)-.5)
         iklu(jfl) = INT(pzt(jfl))
      END DO
      
      !  64 neightboring points for computation of u 
      DO jind1 = 1, 4
         DO jfl = 1, jpnfl
            !  i-direction
            IF( iilu(jfl) <= 2 ) THEN          ;   iidu(jfl,jind1) = jind1
            ELSE
               IF( iilu(jfl) >= jpi-1 ) THEN   ;   iidu(jfl,jind1) = jpi       + jind1 - 4
               ELSE                            ;   iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
               ENDIF
            ENDIF
            !  j-direction
            IF( ijlu(jfl) <= 2 ) THEN          ;   ijdu(jfl,jind1) = jind1
            ELSE
               IF( ijlu(jfl) >= jpj-1 ) THEN   ;   ijdu(jfl,jind1) = jpj       + jind1 - 4
               ELSE                            ;   ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
               ENDIF
            ENDIF
            ! k-direction
            IF( iklu(jfl) <= 2 ) THEN          ;   ikdu(jfl,jind1) = jind1
            ELSE
               IF( iklu(jfl) >= jpk-1 ) THEN   ;   ikdu(jfl,jind1) = jpk + jind1 - 4
               ELSE                            ;   ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
               ENDIF
            ENDIF
         END DO
      END DO
      
      ! Lagrange coefficients
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            zlagxu(jfl,jind1) = 1.
            zlagyu(jfl,jind1) = 1.
            zlagzu(jfl,jind1) = 1.
         END DO
      END DO
      DO jind1 = 1, 4
         DO jind2 = 1, 4
            DO jfl= 1, jpnfl
               IF( jind1 /= jind2 ) THEN
                  zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
                  zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
                  zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
               ENDIF
            END DO
         END DO
      END DO
      
      ! velocity when we compute at middle time step
      
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1, 4
                  ztufl(jfl,jind1,jind2,jind3) =   &
                     &   (  tcoef1(ki) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) +   &
                     &      tcoef2(ki) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) )   &
                     &      / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 
               END DO
            END DO
         END DO
         
         zsumu = 0.
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1, 4
                  zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2)   &
                     &  * zlagzu(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3) 
               END DO
            END DO
         END DO
         pufl(jfl) = zsumu
      END DO
      
      ! Interpolation of V velocity 

      ! nearest neightboring point for computation of v 
      DO jfl = 1, jpnfl
         iilv(jfl) = INT(pxt(jfl)-.5)
         ijlv(jfl) = INT(pyt(jfl)-.5)
         iklv(jfl) = INT(pzt(jfl))
      END DO
      
      ! 64 neightboring points for computation of v 
      DO jind1 = 1, 4
         DO jfl = 1, jpnfl
            ! i-direction
            IF( iilv(jfl) <= 2 ) THEN          ;   iidv(jfl,jind1) = jind1
            ELSE
               IF( iilv(jfl) >= jpi-1 ) THEN   ;   iidv(jfl,jind1) = jpi       + jind1 - 4
               ELSE                            ;   iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
               ENDIF
            ENDIF
            ! j-direction
            IF( ijlv(jfl) <= 2 ) THEN          ;   ijdv(jfl,jind1) = jind1
            ELSE
               IF( ijlv(jfl) >= jpj-1 ) THEN   ;   ijdv(jfl,jind1) = jpj       + jind1 - 4
               ELSE                            ;   ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
               ENDIF
            ENDIF
            ! k-direction
            IF( iklv(jfl) <= 2 ) THEN          ;   ikdv(jfl,jind1) = jind1
            ELSE
               IF( iklv(jfl) >= jpk-1 ) THEN   ;   ikdv(jfl,jind1) = jpk + jind1 - 4
               ELSE                            ;   ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
               ENDIF
            ENDIF
         END DO
      END DO
      
      ! Lagrange coefficients
      
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            zlagxv(jfl,jind1) = 1.
            zlagyv(jfl,jind1) = 1.
            zlagzv(jfl,jind1) = 1.
         END DO
      END DO
      
      DO jind1 = 1, 4
         DO jind2 = 1, 4
            DO jfl = 1, jpnfl
               IF( jind1 /= jind2 ) THEN
                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))   ) )
                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))   ) )
               ENDIF
            END DO
         END DO
      END DO
      
      ! velocity when we compute at middle time step
      
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1 ,4
                  ztvfl(jfl,jind1,jind2,jind3)=   &
                     &   ( tcoef1(ki) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3))  +   &
                     &     tcoef2(ki) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) )    & 
                     &     / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
               END DO
            END DO
         END DO
         
         zsumv=0.
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1, 4
                  zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2)   &
                     &  * zlagzv(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
               END DO
            END DO
         END DO
         pvfl(jfl) = zsumv
      END DO
      
      ! Interpolation of W velocity

      ! nearest neightboring point for computation of w 
      DO jfl = 1, jpnfl
         iilw(jfl) = INT( pxt(jfl)   )
         ijlw(jfl) = INT( pyt(jfl)   )
         iklw(jfl) = INT( pzt(jfl)+.5)
      END DO
      
      ! 64 neightboring points for computation of w 
      DO jind1 = 1, 4
         DO jfl = 1, jpnfl
            ! i-direction
            IF( iilw(jfl) <= 2 ) THEN          ;   iidw(jfl,jind1) = jind1
            ELSE
               IF( iilw(jfl) >= jpi-1 ) THEN   ;   iidw(jfl,jind1) = jpi       + jind1 - 4
               ELSE                            ;   iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
               ENDIF
            ENDIF
            ! j-direction
            IF( ijlw(jfl) <= 2 ) THEN          ;   ijdw(jfl,jind1) = jind1
            ELSE
               IF( ijlw(jfl) >= jpj-1 ) THEN   ;   ijdw(jfl,jind1) = jpj       + jind1 - 4
               ELSE                            ;   ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
               ENDIF
            ENDIF
            ! k-direction
            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
            ELSE
               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
               ENDIF
            ENDIF
         END DO
      END DO
      DO jind1 = 1, 4
         DO jfl = 1, jpnfl
            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
            ELSE
               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
               ENDIF
            ENDIF
         END DO
      END DO
      
      ! Lagrange coefficients  for w interpolation
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            zlagxw(jfl,jind1) = 1.
            zlagyw(jfl,jind1) = 1.
            zlagzw(jfl,jind1) = 1.
         END DO
      END DO
      DO jind1 = 1, 4
         DO jind2 = 1, 4
            DO jfl = 1, jpnfl
               IF( jind1 /= jind2 ) THEN
                  zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2))   ) )
                  zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2))   ) )
                  zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
               ENDIF
            END DO
         END DO
      END DO
      
      ! velocity w  when we compute at middle time step
      DO jfl = 1, jpnfl
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1, 4
                  ztwfl(jfl,jind1,jind2,jind3)=   &
                     &   ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
                     &     tcoef2(ki) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) )  &
                     &   / fse3w(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))
               END DO
            END DO
         END DO
         
         zsumw = 0.e0
         DO jind1 = 1, 4
            DO jind2 = 1, 4
               DO jind3 = 1, 4
                  zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2)   &
                     &  * zlagzw(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
               END DO
            END DO
         END DO
         pwfl(jfl) = zsumw
      END DO
      !   
      CALL wrk_dealloc( jpnfl,    iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
      CALL wrk_dealloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
      CALL wrk_dealloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
      CALL wrk_dealloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
      !   
   END SUBROUTINE flo_interp

#  else
   !!----------------------------------------------------------------------
   !!   No floats                                              Dummy module
   !!----------------------------------------------------------------------
#endif
   
   !!======================================================================
END MODULE flo4rk