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.
flo4rk.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/FLO – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/FLO/flo4rk.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago

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