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 @ 6140

Last change on this file since 6140 was 6140, checked in by timgraham, 8 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
RevLine 
[3]1MODULE flo4rk
2   !!======================================================================
3   !!                    ***  MODULE  flo4rk  ***
4   !! Ocean floats :   trajectory computation using a 4th order Runge-Kutta
5   !!======================================================================
[5836]6#if   defined key_floats
[3]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]16   USE in_out_manager  ! I/O manager
[3294]17   USE wrk_nemo        ! working array
[3]18
19   IMPLICIT NONE
20   PRIVATE
21
[2528]22   PUBLIC   flo_4rk    ! routine called by floats.F90
[3]23
[2528]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
[3]31   !!----------------------------------------------------------------------
[2528]32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]33   !! $Id$
[2528]34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]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.
[2528]49      !!----------------------------------------------------------------------
50      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[3]51      !!
52      INTEGER ::  jfl, jind           ! dummy loop indices
[3294]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
[3]58      !!---------------------------------------------------------------------
[3294]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
[3]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       
[2528]129      DO  jind = 1, 4         
130     
[3]131         ! for each step we compute the compute the velocity with Lagrange interpolation
[2528]132         CALL flo_interp( zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
[3]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
[2528]155      !
[3294]156      CALL wrk_dealloc( jpnfl,    zgifl , zgjfl , zgkfl  , zufl, zvfl, zwfl)
157      CALL wrk_dealloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
158      !
[3]159   END SUBROUTINE flo_4rk
160
161
162   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
[2528]163      &                   pufl, pvfl, pwfl, ki )
[3]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.
[2528]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                 !
[3]177      !!
[2528]178      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
179      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
[3294]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
[3]190      !!---------------------------------------------------------------------
[3294]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 
[3]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
[2528]209            IF( iilu(jfl) <= 2 ) THEN          ;   iidu(jfl,jind1) = jind1
[3]210            ELSE
[2528]211               IF( iilu(jfl) >= jpi-1 ) THEN   ;   iidu(jfl,jind1) = jpi       + jind1 - 4
212               ELSE                            ;   iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
[3]213               ENDIF
214            ENDIF
215            !  j-direction
[2528]216            IF( ijlu(jfl) <= 2 ) THEN          ;   ijdu(jfl,jind1) = jind1
[3]217            ELSE
[2528]218               IF( ijlu(jfl) >= jpj-1 ) THEN   ;   ijdu(jfl,jind1) = jpj       + jind1 - 4
219               ELSE                            ;   ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
[3]220               ENDIF
221            ENDIF
222            ! k-direction
[2528]223            IF( iklu(jfl) <= 2 ) THEN          ;   ikdu(jfl,jind1) = jind1
[3]224            ELSE
[2528]225               IF( iklu(jfl) >= jpk-1 ) THEN   ;   ikdu(jfl,jind1) = jpk + jind1 - 4
226               ELSE                            ;   ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
[3]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) =   &
[2528]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)) )   &
[3]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
[2528]291            IF( iilv(jfl) <= 2 ) THEN          ;   iidv(jfl,jind1) = jind1
[3]292            ELSE
[2528]293               IF( iilv(jfl) >= jpi-1 ) THEN   ;   iidv(jfl,jind1) = jpi       + jind1 - 4
294               ELSE                            ;   iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
[3]295               ENDIF
296            ENDIF
297            ! j-direction
[2528]298            IF( ijlv(jfl) <= 2 ) THEN          ;   ijdv(jfl,jind1) = jind1
[3]299            ELSE
[2528]300               IF( ijlv(jfl) >= jpj-1 ) THEN   ;   ijdv(jfl,jind1) = jpj       + jind1 - 4
301               ELSE                            ;   ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
[3]302               ENDIF
303            ENDIF
304            ! k-direction
[2528]305            IF( iklv(jfl) <= 2 ) THEN          ;   ikdv(jfl,jind1) = jind1
[3]306            ELSE
[2528]307               IF( iklv(jfl) >= jpk-1 ) THEN   ;   ikdv(jfl,jind1) = jpk + jind1 - 4
308               ELSE                            ;   ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
[3]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
[2528]328                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))   ) )
[3]329                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
[2528]330                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))   ) )
[3]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)=   &
[2528]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)) )    & 
[3]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
[2528]366         iilw(jfl) = INT( pxt(jfl)   )
367         ijlw(jfl) = INT( pyt(jfl)   )
368         iklw(jfl) = INT( pzt(jfl)+.5)
[3]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
[2528]375            IF( iilw(jfl) <= 2 ) THEN          ;   iidw(jfl,jind1) = jind1
[3]376            ELSE
[2528]377               IF( iilw(jfl) >= jpi-1 ) THEN   ;   iidw(jfl,jind1) = jpi       + jind1 - 4
378               ELSE                            ;   iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
[3]379               ENDIF
380            ENDIF
381            ! j-direction
[2528]382            IF( ijlw(jfl) <= 2 ) THEN          ;   ijdw(jfl,jind1) = jind1
[3]383            ELSE
[2528]384               IF( ijlw(jfl) >= jpj-1 ) THEN   ;   ijdw(jfl,jind1) = jpj       + jind1 - 4
385               ELSE                            ;   ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
[3]386               ENDIF
387            ENDIF
388            ! k-direction
[2528]389            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
[3]390            ELSE
[2528]391               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
392               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]393               ENDIF
394            ENDIF
395         END DO
396      END DO
397      DO jind1 = 1, 4
398         DO jfl = 1, jpnfl
[2528]399            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
[3]400            ELSE
[2528]401               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
402               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]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
[2528]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))   ) )
[3]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)=   &
[2528]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)) )  &
[6140]436                     &   / e3w_n(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))
[3]437               END DO
438            END DO
439         END DO
440         
[2528]441         zsumw = 0.e0
[3]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
[2528]452      !   
[3294]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      !   
[3]458   END SUBROUTINE flo_interp
459
460#  else
461   !!----------------------------------------------------------------------
[2528]462   !!   No floats                                              Dummy module
[3]463   !!----------------------------------------------------------------------
464#endif
465   
466   !!======================================================================
467END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.