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 branches/2011/dev_r2802_MERCATOR9_floats/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: branches/2011/dev_r2802_MERCATOR9_floats/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90 @ 2907

Last change on this file since 2907 was 2907, checked in by cbricaud, 13 years ago

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