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

Last change on this file since 6140 was 6140, checked in by timgraham, 5 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

  • Property svn:keywords set to Id
File size: 19.7 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
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   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE flo_4rk( kt )
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE flo_4rk  ***
41      !!
42      !!  ** Purpose :   Compute the geographical position (lat,lon,depth)
43      !!       of each float at each time step.
44      !!
45      !!  ** Method  :   The position of a float is computed with a 4th order
46      !!       Runge-Kutta scheme and and Lagrange interpolation.
47      !!         We need to know the velocity field, the old positions of the
48      !!       floats and the grid defined on the domain.
49      !!----------------------------------------------------------------------
50      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
51      !!
52      INTEGER ::  jfl, jind           ! dummy loop indices
53      INTEGER ::  ierror              ! error value
54
55      REAL(wp), POINTER, DIMENSION(:)   ::   zgifl , zgjfl , zgkfl    ! index RK  positions
56      REAL(wp), POINTER, DIMENSION(:)   ::   zufl  , zvfl  , zwfl     ! interpolated velocity at the float position
57      REAL(wp), POINTER, DIMENSION(:,:) ::   zrkxfl, zrkyfl, zrkzfl   ! RK coefficients
58      !!---------------------------------------------------------------------
59      CALL wrk_alloc( jpnfl,    zgifl , zgjfl , zgkfl  , zufl, zvfl, zwfl)
60      CALL wrk_alloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
61      !
62      IF( ierror /= 0 ) THEN
63         WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
64      ENDIF
65
66   
67      IF( kt == nit000 ) THEN
68         IF(lwp) WRITE(numout,*)
69         IF(lwp) WRITE(numout,*) 'flo_4rk : compute Runge Kutta trajectories for floats '
70         IF(lwp) WRITE(numout,*) '~~~~~~~'
71      ENDIF
72
73      ! Verification of the floats positions. If one of them leave the domain
74      ! domain we replace the float near the border.
75      DO jfl = 1, jpnfl
76         ! i-direction
77         IF( tpifl(jfl) <= 1.5 ) THEN
78            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
79            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the WEST border.'
80            tpifl(jfl) = tpifl(jfl) + 1.
81            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=',tpifl(jfl)
82         ENDIF
83         
84         IF( tpifl(jfl) >= jpi-.5 ) THEN
85            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
86            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the EAST border.'
87            tpifl(jfl) = tpifl(jfl) - 1.
88            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=', tpifl(jfl)
89         ENDIF
90         ! j-direction
91         IF( tpjfl(jfl) <= 1.5 ) THEN
92            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
93            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the SOUTH border.'
94            tpjfl(jfl) = tpjfl(jfl) + 1.
95            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
96         ENDIF
97           
98         IF( tpjfl(jfl) >= jpj-.5 ) THEN
99            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
100            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the NORTH border.'
101            tpjfl(jfl) = tpjfl(jfl) - 1.
102            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
103         ENDIF
104         ! k-direction
105         IF( tpkfl(jfl) <= .5 ) THEN
106            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
107            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the TOP border.'
108            tpkfl(jfl) = tpkfl(jfl) + 1.
109            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
110         ENDIF
111         
112         IF( tpkfl(jfl) >= jpk-.5 )  THEN
113            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
114            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the BOTTOM border.'
115            tpkfl(jfl) = tpkfl(jfl) - 1.
116            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
117         ENDIF
118      END DO
119     
120      ! 4 steps of Runge-Kutta algorithme
121      ! initialisation of the positions
122     
123      DO jfl = 1, jpnfl
124         zgifl(jfl) = tpifl(jfl)
125         zgjfl(jfl) = tpjfl(jfl)
126         zgkfl(jfl) = tpkfl(jfl)
127      END DO
128       
129      DO  jind = 1, 4         
130     
131         ! for each step we compute the compute the velocity with Lagrange interpolation
132         CALL flo_interp( zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
133         
134         ! computation of Runge-Kutta factor
135         DO jfl = 1, jpnfl
136            zrkxfl(jfl,jind) = rdt*zufl(jfl)
137            zrkyfl(jfl,jind) = rdt*zvfl(jfl)
138            zrkzfl(jfl,jind) = rdt*zwfl(jfl)
139         END DO
140         IF( jind /= 4 ) THEN
141            DO jfl = 1, jpnfl
142               zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
143               zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
144               zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
145            END DO
146         ENDIF
147      END DO
148      DO jind = 1, 4
149         DO jfl = 1, jpnfl
150            tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
151            tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
152            tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
153         END DO
154      END DO
155      !
156      CALL wrk_dealloc( jpnfl,    zgifl , zgjfl , zgkfl  , zufl, zvfl, zwfl)
157      CALL wrk_dealloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
158      !
159   END SUBROUTINE flo_4rk
160
161
162   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
163      &                   pufl, pvfl, pwfl, ki )
164      !!----------------------------------------------------------------------
165      !!                ***  ROUTINE flointerp  ***
166      !!
167      !! ** Purpose :   Interpolation of the velocity on the float position
168      !!
169      !! ** Method  :   Lagrange interpolation with the 64 neighboring
170      !!      points. This routine is call 4 time at each time step to
171      !!      compute velocity at the date and the position we need to
172      !!      integrated with RK method.
173      !!----------------------------------------------------------------------
174      REAL(wp) , DIMENSION(jpnfl), INTENT(in   ) ::   pxt , pyt , pzt    ! position of the float
175      REAL(wp) , DIMENSION(jpnfl), INTENT(  out) ::   pufl, pvfl, pwfl   ! velocity at this position
176      INTEGER                    , INTENT(in   ) ::   ki                 !
177      !!
178      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
179      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
180      INTEGER  , POINTER, DIMENSION(:)       ::   iilu, ijlu, iklu   ! nearest neighbour INDEX-u
181      INTEGER  , POINTER, DIMENSION(:)       ::   iilv, ijlv, iklv   ! nearest neighbour INDEX-v
182      INTEGER  , POINTER, DIMENSION(:)       ::   iilw, ijlw, iklw   ! nearest neighbour INDEX-w
183      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidu, ijdu, ikdu   ! 64 nearest neighbour INDEX-u
184      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidv, ijdv, ikdv   ! 64 nearest neighbour INDEX-v
185      INTEGER  , POINTER, DIMENSION(:,:)     ::   iidw, ijdw, ikdw   ! 64 nearest neighbour INDEX-w
186      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxu, zlagyu, zlagzu   ! Lagrange  coefficients
187      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxv, zlagyv, zlagzv   !    -           -
188      REAL(wp) , POINTER, DIMENSION(:,:)     ::   zlagxw, zlagyw, zlagzw   !    -           -
189      REAL(wp) , POINTER, DIMENSION(:,:,:,:) ::   ztufl , ztvfl , ztwfl   ! velocity at choosen time step
190      !!---------------------------------------------------------------------
191      CALL wrk_alloc( jpnfl,    iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
192      CALL wrk_alloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
193      CALL wrk_alloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
194      CALL wrk_alloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
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                     &   / e3w_n(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      CALL wrk_dealloc( jpnfl,    iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
454      CALL wrk_dealloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
455      CALL wrk_dealloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
456      CALL wrk_dealloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
457      !   
458   END SUBROUTINE flo_interp
459
460#  else
461   !!----------------------------------------------------------------------
462   !!   No floats                                              Dummy module
463   !!----------------------------------------------------------------------
464#endif
465   
466   !!======================================================================
467END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.