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 trunk/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90 @ 4528

Last change on this file since 4528 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

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