source: NEMO/trunk/src/OCE/FLO/flo4rk.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 7 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 18.5 KB
Line 
1MODULE flo4rk
2   !!======================================================================
3   !!                    ***  MODULE  flo4rk  ***
4   !! Ocean floats :   trajectory computation using a 4th order Runge-Kutta
5   !!======================================================================
6   !!
7   !!----------------------------------------------------------------------
8   !!   flo_4rk        : Compute the geographical position of floats
9   !!   flo_interp     : interpolation
10   !!----------------------------------------------------------------------
11   USE flo_oce         ! ocean drifting floats
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE in_out_manager  ! I/O manager
15
16   IMPLICIT NONE
17   PRIVATE
18
19   PUBLIC   flo_4rk    ! routine called by floats.F90
20
21   !                                   ! RK4 and Lagrange interpolation coefficients
22   REAL(wp), DIMENSION (4) ::   tcoef1 = (/  1.0  ,  0.5  ,  0.5  ,  0.0  /)   !
23   REAL(wp), DIMENSION (4) ::   tcoef2 = (/  0.0  ,  0.5  ,  0.5  ,  1.0  /)   !
24   REAL(wp), DIMENSION (4) ::   scoef2 = (/  1.0  ,  2.0  ,  2.0  ,  1.0  /)   !
25   REAL(wp), DIMENSION (4) ::   rcoef  = (/-1./6. , 1./2. ,-1./2. , 1./6. /)   !
26   REAL(wp), DIMENSION (3) ::   scoef1 = (/  0.5  ,  0.5  ,  1.0  /)           !
27
28   !!----------------------------------------------------------------------
29   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
30   !! $Id$
31   !! Software governed by the CeCILL license (see ./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   SUBROUTINE flo_4rk( kt, Kbb, Kmm )
36      !!----------------------------------------------------------------------
37      !!                  ***  ROUTINE flo_4rk  ***
38      !!
39      !!  ** Purpose :   Compute the geographical position (lat,lon,depth)
40      !!       of each float at each time step.
41      !!
42      !!  ** Method  :   The position of a float is computed with a 4th order
43      !!       Runge-Kutta scheme and and Lagrange interpolation.
44      !!         We need to know the velocity field, the old positions of the
45      !!       floats and the grid defined on the domain.
46      !!----------------------------------------------------------------------
47      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
48      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
49      !!
50      INTEGER ::  jfl, jind           ! dummy loop indices
51      INTEGER ::  ierror              ! error value
52
53      REAL(wp), DIMENSION(jpnfl)   ::   zgifl , zgjfl , zgkfl    ! index RK  positions
54      REAL(wp), DIMENSION(jpnfl)   ::   zufl  , zvfl  , zwfl     ! interpolated velocity at the float position
55      REAL(wp), DIMENSION(jpnfl,4) ::   zrkxfl, zrkyfl, zrkzfl   ! RK coefficients
56      !!---------------------------------------------------------------------
57      !
58      IF( ierror /= 0 ) THEN
59         WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
60      ENDIF
61
62   
63      IF( kt == nit000 ) THEN
64         IF(lwp) WRITE(numout,*)
65         IF(lwp) WRITE(numout,*) 'flo_4rk : compute Runge Kutta trajectories for floats '
66         IF(lwp) WRITE(numout,*) '~~~~~~~'
67      ENDIF
68
69      ! Verification of the floats positions. If one of them leave the domain
70      ! domain we replace the float near the border.
71      DO jfl = 1, jpnfl
72         ! i-direction
73         IF( tpifl(jfl) <= 1.5 ) THEN
74            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
75            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the WEST border.'
76            tpifl(jfl) = tpifl(jfl) + 1.
77            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=',tpifl(jfl)
78         ENDIF
79         
80         IF( tpifl(jfl) >= jpi-.5 ) THEN
81            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
82            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the EAST border.'
83            tpifl(jfl) = tpifl(jfl) - 1.
84            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=', tpifl(jfl)
85         ENDIF
86         ! j-direction
87         IF( tpjfl(jfl) <= 1.5 ) THEN
88            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
89            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the SOUTH border.'
90            tpjfl(jfl) = tpjfl(jfl) + 1.
91            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
92         ENDIF
93           
94         IF( tpjfl(jfl) >= jpj-.5 ) THEN
95            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
96            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the NORTH border.'
97            tpjfl(jfl) = tpjfl(jfl) - 1.
98            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
99         ENDIF
100         ! k-direction
101         IF( tpkfl(jfl) <= .5 ) THEN
102            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
103            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the TOP border.'
104            tpkfl(jfl) = tpkfl(jfl) + 1.
105            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
106         ENDIF
107         
108         IF( tpkfl(jfl) >= jpk-.5 )  THEN
109            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
110            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the BOTTOM border.'
111            tpkfl(jfl) = tpkfl(jfl) - 1.
112            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
113         ENDIF
114      END DO
115     
116      ! 4 steps of Runge-Kutta algorithme
117      ! initialisation of the positions
118     
119      DO jfl = 1, jpnfl
120         zgifl(jfl) = tpifl(jfl)
121         zgjfl(jfl) = tpjfl(jfl)
122         zgkfl(jfl) = tpkfl(jfl)
123      END DO
124       
125      DO  jind = 1, 4         
126     
127         ! for each step we compute the compute the velocity with Lagrange interpolation
128         CALL flo_interp( Kbb, Kmm, zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
129         
130         ! computation of Runge-Kutta factor
131         DO jfl = 1, jpnfl
132            zrkxfl(jfl,jind) = rn_Dt*zufl(jfl)
133            zrkyfl(jfl,jind) = rn_Dt*zvfl(jfl)
134            zrkzfl(jfl,jind) = rn_Dt*zwfl(jfl)
135         END DO
136         IF( jind /= 4 ) THEN
137            DO jfl = 1, jpnfl
138               zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
139               zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
140               zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
141            END DO
142         ENDIF
143      END DO
144      DO jind = 1, 4
145         DO jfl = 1, jpnfl
146            tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
147            tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
148            tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
149         END DO
150      END DO
151      !
152      !
153   END SUBROUTINE flo_4rk
154
155
156   SUBROUTINE flo_interp( Kbb, Kmm,              &
157      &                   pxt , pyt , pzt ,      &
158      &                   pufl, pvfl, pwfl, ki )
159      !!----------------------------------------------------------------------
160      !!                ***  ROUTINE flointerp  ***
161      !!
162      !! ** Purpose :   Interpolation of the velocity on the float position
163      !!
164      !! ** Method  :   Lagrange interpolation with the 64 neighboring
165      !!      points. This routine is call 4 time at each time step to
166      !!      compute velocity at the date and the position we need to
167      !!      integrated with RK method.
168      !!----------------------------------------------------------------------
169      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm           ! ocean time level indices
170      REAL(wp) , DIMENSION(jpnfl), INTENT(in   ) ::   pxt , pyt , pzt    ! position of the float
171      REAL(wp) , DIMENSION(jpnfl), INTENT(  out) ::   pufl, pvfl, pwfl   ! velocity at this position
172      INTEGER                    , INTENT(in   ) ::   ki                 !
173      !!
174      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
175      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
176      INTEGER  , DIMENSION(jpnfl)       ::   iilu, ijlu, iklu   ! nearest neighbour INDEX-u
177      INTEGER  , DIMENSION(jpnfl)       ::   iilv, ijlv, iklv   ! nearest neighbour INDEX-v
178      INTEGER  , DIMENSION(jpnfl)       ::   iilw, ijlw, iklw   ! nearest neighbour INDEX-w
179      INTEGER  , DIMENSION(jpnfl,4)     ::   iidu, ijdu, ikdu   ! 64 nearest neighbour INDEX-u
180      INTEGER  , DIMENSION(jpnfl,4)     ::   iidv, ijdv, ikdv   ! 64 nearest neighbour INDEX-v
181      INTEGER  , DIMENSION(jpnfl,4)     ::   iidw, ijdw, ikdw   ! 64 nearest neighbour INDEX-w
182      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxu, zlagyu, zlagzu   ! Lagrange  coefficients
183      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxv, zlagyv, zlagzv   !    -           -
184      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxw, zlagyw, zlagzw   !    -           -
185      REAL(wp) , DIMENSION(jpnfl,4,4,4) ::   ztufl , ztvfl , ztwfl   ! velocity at choosen time step
186      !!---------------------------------------------------------------------
187 
188      ! Interpolation of U velocity
189
190      ! nearest neightboring point for computation of u       
191      DO jfl = 1, jpnfl
192         iilu(jfl) = INT(pxt(jfl)-.5)
193         ijlu(jfl) = INT(pyt(jfl)-.5)
194         iklu(jfl) = INT(pzt(jfl))
195      END DO
196     
197      !  64 neightboring points for computation of u
198      DO jind1 = 1, 4
199         DO jfl = 1, jpnfl
200            !  i-direction
201            IF( iilu(jfl) <= 2 ) THEN          ;   iidu(jfl,jind1) = jind1
202            ELSE
203               IF( iilu(jfl) >= jpi-1 ) THEN   ;   iidu(jfl,jind1) = jpi       + jind1 - 4
204               ELSE                            ;   iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
205               ENDIF
206            ENDIF
207            !  j-direction
208            IF( ijlu(jfl) <= 2 ) THEN          ;   ijdu(jfl,jind1) = jind1
209            ELSE
210               IF( ijlu(jfl) >= jpj-1 ) THEN   ;   ijdu(jfl,jind1) = jpj       + jind1 - 4
211               ELSE                            ;   ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
212               ENDIF
213            ENDIF
214            ! k-direction
215            IF( iklu(jfl) <= 2 ) THEN          ;   ikdu(jfl,jind1) = jind1
216            ELSE
217               IF( iklu(jfl) >= jpk-1 ) THEN   ;   ikdu(jfl,jind1) = jpk + jind1 - 4
218               ELSE                            ;   ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
219               ENDIF
220            ENDIF
221         END DO
222      END DO
223     
224      ! Lagrange coefficients
225      DO jfl = 1, jpnfl
226         DO jind1 = 1, 4
227            zlagxu(jfl,jind1) = 1.
228            zlagyu(jfl,jind1) = 1.
229            zlagzu(jfl,jind1) = 1.
230         END DO
231      END DO
232      DO jind1 = 1, 4
233         DO jind2 = 1, 4
234            DO jfl= 1, jpnfl
235               IF( jind1 /= jind2 ) THEN
236                  zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
237                  zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
238                  zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
239               ENDIF
240            END DO
241         END DO
242      END DO
243     
244      ! velocity when we compute at middle time step
245     
246      DO jfl = 1, jpnfl
247         DO jind1 = 1, 4
248            DO jind2 = 1, 4
249               DO jind3 = 1, 4
250                  ztufl(jfl,jind1,jind2,jind3) =   &
251                     &   (  tcoef1(ki) * uu(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3),Kbb) +   &
252                     &      tcoef2(ki) * uu(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3),Kmm) )   &
253                     &      / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 
254               END DO
255            END DO
256         END DO
257         
258         zsumu = 0.
259         DO jind1 = 1, 4
260            DO jind2 = 1, 4
261               DO jind3 = 1, 4
262                  zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2)   &
263                     &  * zlagzu(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3) 
264               END DO
265            END DO
266         END DO
267         pufl(jfl) = zsumu
268      END DO
269     
270      ! Interpolation of V velocity
271
272      ! nearest neightboring point for computation of v
273      DO jfl = 1, jpnfl
274         iilv(jfl) = INT(pxt(jfl)-.5)
275         ijlv(jfl) = INT(pyt(jfl)-.5)
276         iklv(jfl) = INT(pzt(jfl))
277      END DO
278     
279      ! 64 neightboring points for computation of v
280      DO jind1 = 1, 4
281         DO jfl = 1, jpnfl
282            ! i-direction
283            IF( iilv(jfl) <= 2 ) THEN          ;   iidv(jfl,jind1) = jind1
284            ELSE
285               IF( iilv(jfl) >= jpi-1 ) THEN   ;   iidv(jfl,jind1) = jpi       + jind1 - 4
286               ELSE                            ;   iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
287               ENDIF
288            ENDIF
289            ! j-direction
290            IF( ijlv(jfl) <= 2 ) THEN          ;   ijdv(jfl,jind1) = jind1
291            ELSE
292               IF( ijlv(jfl) >= jpj-1 ) THEN   ;   ijdv(jfl,jind1) = jpj       + jind1 - 4
293               ELSE                            ;   ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
294               ENDIF
295            ENDIF
296            ! k-direction
297            IF( iklv(jfl) <= 2 ) THEN          ;   ikdv(jfl,jind1) = jind1
298            ELSE
299               IF( iklv(jfl) >= jpk-1 ) THEN   ;   ikdv(jfl,jind1) = jpk + jind1 - 4
300               ELSE                            ;   ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
301               ENDIF
302            ENDIF
303         END DO
304      END DO
305     
306      ! Lagrange coefficients
307     
308      DO jfl = 1, jpnfl
309         DO jind1 = 1, 4
310            zlagxv(jfl,jind1) = 1.
311            zlagyv(jfl,jind1) = 1.
312            zlagzv(jfl,jind1) = 1.
313         END DO
314      END DO
315     
316      DO jind1 = 1, 4
317         DO jind2 = 1, 4
318            DO jfl = 1, jpnfl
319               IF( jind1 /= jind2 ) THEN
320                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))   ) )
321                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
322                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))   ) )
323               ENDIF
324            END DO
325         END DO
326      END DO
327     
328      ! velocity when we compute at middle time step
329     
330      DO jfl = 1, jpnfl
331         DO jind1 = 1, 4
332            DO jind2 = 1, 4
333               DO jind3 = 1 ,4
334                  ztvfl(jfl,jind1,jind2,jind3)=   &
335                     &   ( tcoef1(ki) * vv(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3),Kbb)  +   &
336                     &     tcoef2(ki) * vv(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3),Kmm) )    & 
337                     &     / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
338               END DO
339            END DO
340         END DO
341         
342         zsumv=0.
343         DO jind1 = 1, 4
344            DO jind2 = 1, 4
345               DO jind3 = 1, 4
346                  zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2)   &
347                     &  * zlagzv(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
348               END DO
349            END DO
350         END DO
351         pvfl(jfl) = zsumv
352      END DO
353     
354      ! Interpolation of W velocity
355
356      ! nearest neightboring point for computation of w
357      DO jfl = 1, jpnfl
358         iilw(jfl) = INT( pxt(jfl)   )
359         ijlw(jfl) = INT( pyt(jfl)   )
360         iklw(jfl) = INT( pzt(jfl)+.5)
361      END DO
362     
363      ! 64 neightboring points for computation of w
364      DO jind1 = 1, 4
365         DO jfl = 1, jpnfl
366            ! i-direction
367            IF( iilw(jfl) <= 2 ) THEN          ;   iidw(jfl,jind1) = jind1
368            ELSE
369               IF( iilw(jfl) >= jpi-1 ) THEN   ;   iidw(jfl,jind1) = jpi       + jind1 - 4
370               ELSE                            ;   iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
371               ENDIF
372            ENDIF
373            ! j-direction
374            IF( ijlw(jfl) <= 2 ) THEN          ;   ijdw(jfl,jind1) = jind1
375            ELSE
376               IF( ijlw(jfl) >= jpj-1 ) THEN   ;   ijdw(jfl,jind1) = jpj       + jind1 - 4
377               ELSE                            ;   ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
378               ENDIF
379            ENDIF
380            ! k-direction
381            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
382            ELSE
383               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
384               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
385               ENDIF
386            ENDIF
387         END DO
388      END DO
389      DO jind1 = 1, 4
390         DO jfl = 1, jpnfl
391            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
392            ELSE
393               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
394               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
395               ENDIF
396            ENDIF
397         END DO
398      END DO
399     
400      ! Lagrange coefficients  for w interpolation
401      DO jfl = 1, jpnfl
402         DO jind1 = 1, 4
403            zlagxw(jfl,jind1) = 1.
404            zlagyw(jfl,jind1) = 1.
405            zlagzw(jfl,jind1) = 1.
406         END DO
407      END DO
408      DO jind1 = 1, 4
409         DO jind2 = 1, 4
410            DO jfl = 1, jpnfl
411               IF( jind1 /= jind2 ) THEN
412                  zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2))   ) )
413                  zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2))   ) )
414                  zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
415               ENDIF
416            END DO
417         END DO
418      END DO
419     
420      ! velocity w  when we compute at middle time step
421      DO jfl = 1, jpnfl
422         DO jind1 = 1, 4
423            DO jind2 = 1, 4
424               DO jind3 = 1, 4
425                  ztwfl(jfl,jind1,jind2,jind3)=   &
426                     &   ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
427                     &     tcoef2(ki) * ww(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) )  &
428                     &   / e3w(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3),Kmm)
429               END DO
430            END DO
431         END DO
432         
433         zsumw = 0.e0
434         DO jind1 = 1, 4
435            DO jind2 = 1, 4
436               DO jind3 = 1, 4
437                  zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2)   &
438                     &  * zlagzw(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
439               END DO
440            END DO
441         END DO
442         pwfl(jfl) = zsumw
443      END DO
444      !   
445      !   
446   END SUBROUTINE flo_interp
447
448   !!======================================================================
449END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.