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/trunk/src/OCE/FLO – NEMO

source: NEMO/trunk/src/OCE/FLO/flo4rk.F90 @ 13891

Last change on this file since 13891 was 13237, checked in by smasson, 4 years ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 18.6 KB
RevLine 
[3]1MODULE flo4rk
2   !!======================================================================
3   !!                    ***  MODULE  flo4rk  ***
4   !! Ocean floats :   trajectory computation using a 4th order Runge-Kutta
5   !!======================================================================
[11536]6   !!
[3]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
[16]14   USE in_out_manager  ! I/O manager
[3]15
16   IMPLICIT NONE
17   PRIVATE
18
[2528]19   PUBLIC   flo_4rk    ! routine called by floats.F90
[3]20
[2528]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
[13237]28#  include "domzgr_substitute.h90"
[3]29   !!----------------------------------------------------------------------
[9598]30   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]31   !! $Id$
[10068]32   !! Software governed by the CeCILL license (see ./LICENSE)
[3]33   !!----------------------------------------------------------------------
34CONTAINS
35
[12377]36   SUBROUTINE flo_4rk( kt, Kbb, Kmm )
[3]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.
[2528]47      !!----------------------------------------------------------------------
[12377]48      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
49      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
[3]50      !!
51      INTEGER ::  jfl, jind           ! dummy loop indices
[3294]52      INTEGER ::  ierror              ! error value
53
[9125]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
[3]57      !!---------------------------------------------------------------------
[3294]58      !
59      IF( ierror /= 0 ) THEN
60         WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
61      ENDIF
62
[3]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       
[2528]126      DO  jind = 1, 4         
127     
[3]128         ! for each step we compute the compute the velocity with Lagrange interpolation
[12377]129         CALL flo_interp( Kbb, Kmm, zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
[3]130         
131         ! computation of Runge-Kutta factor
132         DO jfl = 1, jpnfl
[12489]133            zrkxfl(jfl,jind) = rn_Dt*zufl(jfl)
134            zrkyfl(jfl,jind) = rn_Dt*zvfl(jfl)
135            zrkzfl(jfl,jind) = rn_Dt*zwfl(jfl)
[3]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
[2528]152      !
[3294]153      !
[3]154   END SUBROUTINE flo_4rk
155
156
[12377]157   SUBROUTINE flo_interp( Kbb, Kmm,              &
158      &                   pxt , pyt , pzt ,      &
[2528]159      &                   pufl, pvfl, pwfl, ki )
[3]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.
[2528]169      !!----------------------------------------------------------------------
[12377]170      INTEGER                    , INTENT(in   ) ::   Kbb, Kmm           ! ocean time level indices
[2528]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                 !
[3]174      !!
[2528]175      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
176      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
[9125]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
[3]187      !!---------------------------------------------------------------------
[3294]188 
[3]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
[2528]202            IF( iilu(jfl) <= 2 ) THEN          ;   iidu(jfl,jind1) = jind1
[3]203            ELSE
[2528]204               IF( iilu(jfl) >= jpi-1 ) THEN   ;   iidu(jfl,jind1) = jpi       + jind1 - 4
205               ELSE                            ;   iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
[3]206               ENDIF
207            ENDIF
208            !  j-direction
[2528]209            IF( ijlu(jfl) <= 2 ) THEN          ;   ijdu(jfl,jind1) = jind1
[3]210            ELSE
[2528]211               IF( ijlu(jfl) >= jpj-1 ) THEN   ;   ijdu(jfl,jind1) = jpj       + jind1 - 4
212               ELSE                            ;   ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
[3]213               ENDIF
214            ENDIF
215            ! k-direction
[2528]216            IF( iklu(jfl) <= 2 ) THEN          ;   ikdu(jfl,jind1) = jind1
[3]217            ELSE
[2528]218               IF( iklu(jfl) >= jpk-1 ) THEN   ;   ikdu(jfl,jind1) = jpk + jind1 - 4
219               ELSE                            ;   ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
[3]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) =   &
[12377]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) )   &
[3]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
[2528]284            IF( iilv(jfl) <= 2 ) THEN          ;   iidv(jfl,jind1) = jind1
[3]285            ELSE
[2528]286               IF( iilv(jfl) >= jpi-1 ) THEN   ;   iidv(jfl,jind1) = jpi       + jind1 - 4
287               ELSE                            ;   iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
[3]288               ENDIF
289            ENDIF
290            ! j-direction
[2528]291            IF( ijlv(jfl) <= 2 ) THEN          ;   ijdv(jfl,jind1) = jind1
[3]292            ELSE
[2528]293               IF( ijlv(jfl) >= jpj-1 ) THEN   ;   ijdv(jfl,jind1) = jpj       + jind1 - 4
294               ELSE                            ;   ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
[3]295               ENDIF
296            ENDIF
297            ! k-direction
[2528]298            IF( iklv(jfl) <= 2 ) THEN          ;   ikdv(jfl,jind1) = jind1
[3]299            ELSE
[2528]300               IF( iklv(jfl) >= jpk-1 ) THEN   ;   ikdv(jfl,jind1) = jpk + jind1 - 4
301               ELSE                            ;   ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
[3]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
[2528]321                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))   ) )
[3]322                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
[2528]323                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))   ) )
[3]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)=   &
[12377]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) )    & 
[3]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
[2528]359         iilw(jfl) = INT( pxt(jfl)   )
360         ijlw(jfl) = INT( pyt(jfl)   )
361         iklw(jfl) = INT( pzt(jfl)+.5)
[3]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
[2528]368            IF( iilw(jfl) <= 2 ) THEN          ;   iidw(jfl,jind1) = jind1
[3]369            ELSE
[2528]370               IF( iilw(jfl) >= jpi-1 ) THEN   ;   iidw(jfl,jind1) = jpi       + jind1 - 4
371               ELSE                            ;   iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
[3]372               ENDIF
373            ENDIF
374            ! j-direction
[2528]375            IF( ijlw(jfl) <= 2 ) THEN          ;   ijdw(jfl,jind1) = jind1
[3]376            ELSE
[2528]377               IF( ijlw(jfl) >= jpj-1 ) THEN   ;   ijdw(jfl,jind1) = jpj       + jind1 - 4
378               ELSE                            ;   ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
[3]379               ENDIF
380            ENDIF
381            ! k-direction
[2528]382            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
[3]383            ELSE
[2528]384               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
385               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]386               ENDIF
387            ENDIF
388         END DO
389      END DO
390      DO jind1 = 1, 4
391         DO jfl = 1, jpnfl
[2528]392            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
[3]393            ELSE
[2528]394               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
395               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]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
[2528]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))   ) )
[3]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)=   &
[2528]427                     &   ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
[12377]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)
[3]430               END DO
431            END DO
432         END DO
433         
[2528]434         zsumw = 0.e0
[3]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
[2528]445      !   
[3294]446      !   
[3]447   END SUBROUTINE flo_interp
448
449   !!======================================================================
450END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.