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

source: trunk/NEMO/OPA_SRC/FLO/flo4rk.F90 @ 16

Last change on this file since 16 was 16, checked in by opalod, 20 years ago

CT : UPDATE001 : First major NEMO update
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.4 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   !! * Modules used
14   USE flo_oce         ! ocean drifting floats
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE in_out_manager  ! I/O manager
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Accessibility
23   PUBLIC flo_4rk      ! routine called by floats.F90
24
25   !! * Module variables
26   REAL(wp), DIMENSION (4) ::   &   ! RK4 and Lagrange interpolation
27      tcoef1 = (/  1.0  ,  0.5  ,  0.5  ,  0.0  /) ,  &  ! coeffients  for
28      tcoef2 = (/  0.0  ,  0.5  ,  0.5  ,  1.0  /) ,  &  ! lagrangian interp.
29      scoef2 = (/  1.0  ,  2.0  ,  2.0  ,  1.0  /) ,  &  ! RK4 coefficients
30      rcoef  = (/-1./6. , 1./2. ,-1./2. , 1./6. /)       ! ???
31   REAL(wp), DIMENSION (3) ::   &
32      scoef1 = (/ .5, .5, 1. /)       ! compute position with interpolated
33   !!----------------------------------------------------------------------
34   !!   OPA 9.0 , LODYC-IPSL  (2003)
35   !!----------------------------------------------------------------------
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      !!----------------------------------------------------------------------
52      !! * Arguments
53      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
54
55      !! * Local declarations
56      INTEGER ::  jfl, jind           ! dummy loop indices
57      REAL(wp), DIMENSION ( jpnfl)  ::   &
58         zgifl, zgjfl, zgkfl,   &     ! index RK  positions
59         zufl, zvfl, zwfl             ! interpolated velocity at the
60      ! float position
61      REAL(wp), DIMENSION ( jpnfl, 4 )  ::    &
62         zrkxfl, zrkyfl, zrkzfl       ! RK coefficients
63      !!---------------------------------------------------------------------
64   
65      IF( kt == nit000 ) THEN
66         IF(lwp) WRITE(numout,*)
67         IF(lwp) WRITE(numout,*) 'flo_4rk : compute Runge Kutta trajectories for floats '
68         IF(lwp) WRITE(numout,*) '~~~~~~~'
69      ENDIF
70
71      ! Verification of the floats positions. If one of them leave the domain
72      ! domain we replace the float near the border.
73      DO jfl = 1, jpnfl
74         ! i-direction
75         IF( tpifl(jfl) <= 1.5 ) THEN
76            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
77            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the WEST border.'
78            tpifl(jfl) = tpifl(jfl) + 1.
79            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=',tpifl(jfl)
80         ENDIF
81         
82         IF( tpifl(jfl) >= jpi-.5 ) THEN
83            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
84            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the EAST border.'
85            tpifl(jfl) = tpifl(jfl) - 1.
86            IF(lwp)WRITE(numout,*)'New initialisation for this float at i=', tpifl(jfl)
87         ENDIF
88         ! j-direction
89         IF( tpjfl(jfl) <= 1.5 ) THEN
90            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
91            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the SOUTH border.'
92            tpjfl(jfl) = tpjfl(jfl) + 1.
93            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
94         ENDIF
95           
96         IF( tpjfl(jfl) >= jpj-.5 ) THEN
97            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
98            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the NORTH border.'
99            tpjfl(jfl) = tpjfl(jfl) - 1.
100            IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
101         ENDIF
102         ! k-direction
103         IF( tpkfl(jfl) <= .5 ) THEN
104            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
105            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the TOP border.'
106            tpkfl(jfl) = tpkfl(jfl) + 1.
107            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
108         ENDIF
109         
110         IF( tpkfl(jfl) >= jpk-.5 )  THEN
111            IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
112            IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the BOTTOM border.'
113            tpkfl(jfl) = tpkfl(jfl) - 1.
114            IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
115         ENDIF
116      END DO
117     
118      ! 4 steps of Runge-Kutta algorithme
119      ! initialisation of the positions
120     
121      DO jfl = 1, jpnfl
122         zgifl(jfl) = tpifl(jfl)
123         zgjfl(jfl) = tpjfl(jfl)
124         zgkfl(jfl) = tpkfl(jfl)
125      END DO
126       
127      DO  jind = 1, 4
128         
129         ! for each step we compute the compute the velocity with Lagrange interpolation
130         CALL flo_interp(zgifl,zgjfl,zgkfl,zufl,zvfl,zwfl,jind)
131         
132         ! computation of Runge-Kutta factor
133         
134         DO jfl = 1, jpnfl
135            zrkxfl(jfl,jind) = rdt*zufl(jfl)
136            zrkyfl(jfl,jind) = rdt*zvfl(jfl)
137            zrkzfl(jfl,jind) = rdt*zwfl(jfl)
138         END DO
139         IF( jind /= 4 ) THEN
140            DO jfl = 1, jpnfl
141               zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
142               zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
143               zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
144            END DO
145         ENDIF
146      END DO
147      DO jind = 1, 4
148         DO jfl = 1, jpnfl
149            tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
150            tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
151            tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
152         END DO
153      END DO
154
155   END SUBROUTINE flo_4rk
156
157
158   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
159      &                   pufl, pvfl, pwfl, kind )
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.
169      !!
170      !!----------------------------------------------------------------------
171      !! * Local declarations
172      INTEGER ::   &
173         kind,   &
174         jfl, jind1, jind2, jind3,   &
175         zsumu, zsumv, zsumw
176      INTEGER , DIMENSION  ( jpnfl  ) ::   &
177         iilu, ijlu, iklu,     &           ! nearest neighbour INDEX-u
178         iilv, ijlv, iklv,     &           ! nearest neighbour INDEX-v
179         iilw, ijlw, iklw                  ! nearest neighbour INDEX-w
180      INTEGER , DIMENSION  ( jpnfl, 4  ) ::   &
181         iidu, ijdu, ikdu,    &            ! 64 nearest neighbour INDEX-u
182         iidv, ijdv, ikdv,    &            ! 64 nearest neighbour INDEX-v
183         iidw, ijdw, ikdw                  ! 64 nearest neighbour INDEX-w
184      REAL(wp) , DIMENSION  ( jpnfl  ) ::   &
185         pxt , pyt , pzt,     &            ! position of the float
186         pufl, pvfl, pwfl                  ! velocity at this position
187      REAL(wp) , DIMENSION  ( jpnfl, 4, 4, 4 ) ::   &
188         ztufl, ztvfl, ztwfl               ! velocity at choosen time step
189      REAL(wp) , DIMENSION  ( jpnfl, 4 ) ::   &
190         zlagxu, zlagyu, zlagzu,   &       ! Lagrange  coefficients
191         zlagxv, zlagyv, zlagzv,   &
192         zlagxw, zlagyw, zlagzw
193      !!---------------------------------------------------------------------
194     
195      ! Interpolation of U velocity
196
197      ! nearest neightboring point for computation of u       
198      DO jfl = 1, jpnfl
199         iilu(jfl) = INT(pxt(jfl)-.5)
200         ijlu(jfl) = INT(pyt(jfl)-.5)
201         iklu(jfl) = INT(pzt(jfl))
202      END DO
203     
204      !  64 neightboring points for computation of u
205      DO jind1 = 1, 4
206         DO jfl = 1, jpnfl
207            !  i-direction
208            IF( iilu(jfl) <= 2 ) THEN
209               iidu(jfl,jind1) = jind1
210            ELSE
211               IF( iilu(jfl) >= jpi-1 ) THEN
212                  iidu(jfl,jind1) = jpi + jind1 - 4
213               ELSE
214                  iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
215               ENDIF
216            ENDIF
217            !  j-direction
218            IF( ijlu(jfl) <= 2 ) THEN
219               ijdu(jfl,jind1) = jind1
220            ELSE
221               IF( ijlu(jfl) >= jpj-1 ) THEN
222                  ijdu(jfl,jind1) = jpj + jind1 - 4
223               ELSE
224                  ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
225               ENDIF
226            ENDIF
227            ! k-direction
228            IF( iklu(jfl) <= 2 ) THEN
229               ikdu(jfl,jind1) = jind1
230            ELSE
231               IF( iklu(jfl) >= jpk-1 ) THEN
232                  ikdu(jfl,jind1) = jpk + jind1 - 4
233               ELSE
234                  ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
235               ENDIF
236            ENDIF
237         END DO
238      END DO
239     
240      ! Lagrange coefficients
241     
242      DO jfl = 1, jpnfl
243         DO jind1 = 1, 4
244            zlagxu(jfl,jind1) = 1.
245            zlagyu(jfl,jind1) = 1.
246            zlagzu(jfl,jind1) = 1.
247         END DO
248      END DO
249     
250      DO jind1 = 1, 4
251         DO jind2 = 1, 4
252            DO jfl= 1, jpnfl
253               IF( jind1 /= jind2 ) THEN
254                  zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
255                  zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
256                  zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
257               ENDIF
258            END DO
259         END DO
260      END DO
261     
262      ! velocity when we compute at middle time step
263     
264      DO jfl = 1, jpnfl
265         DO jind1 = 1, 4
266            DO jind2 = 1, 4
267               DO jind3 = 1, 4
268                  ztufl(jfl,jind1,jind2,jind3) =   &
269                     &   (  tcoef1(kind) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) +   &
270                     &      tcoef2(kind) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) )   &
271                     &      / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 
272               END DO
273            END DO
274         END DO
275         
276         zsumu = 0.
277         DO jind1 = 1, 4
278            DO jind2 = 1, 4
279               DO jind3 = 1, 4
280                  zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2)   &
281                     &  * zlagzu(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3) 
282               END DO
283            END DO
284         END DO
285         pufl(jfl) = zsumu
286      END DO
287     
288      ! Interpolation of V velocity
289
290      ! nearest neightboring point for computation of v
291      DO jfl = 1, jpnfl
292         iilv(jfl) = INT(pxt(jfl)-.5)
293         ijlv(jfl) = INT(pyt(jfl)-.5)
294         iklv(jfl) = INT(pzt(jfl))
295      END DO
296     
297      ! 64 neightboring points for computation of v
298      DO jind1 = 1, 4
299         DO jfl = 1, jpnfl
300            ! i-direction
301            IF( iilv(jfl) <= 2 ) THEN
302               iidv(jfl,jind1) = jind1
303            ELSE
304               IF( iilv(jfl) >= jpi-1 ) THEN
305                  iidv(jfl,jind1) = jpi + jind1 - 4
306               ELSE
307                  iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
308               ENDIF
309            ENDIF
310            ! j-direction
311            IF( ijlv(jfl) <= 2 ) THEN
312               ijdv(jfl,jind1) = jind1
313            ELSE
314               IF( ijlv(jfl) >= jpj-1 ) THEN
315                  ijdv(jfl,jind1) = jpj + jind1 - 4
316               ELSE
317                  ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
318               ENDIF
319            ENDIF
320            ! k-direction
321            IF( iklv(jfl) <= 2 ) THEN
322               ikdv(jfl,jind1) = jind1
323            ELSE
324               IF( iklv(jfl) >= jpk-1 ) THEN
325                  ikdv(jfl,jind1) = jpk + jind1 - 4
326               ELSE
327                  ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
328               ENDIF
329            ENDIF
330         END DO
331      END DO
332     
333      ! Lagrange coefficients
334     
335      DO jfl = 1, jpnfl
336         DO jind1 = 1, 4
337            zlagxv(jfl,jind1) = 1.
338            zlagyv(jfl,jind1) = 1.
339            zlagzv(jfl,jind1) = 1.
340         END DO
341      END DO
342     
343      DO jind1 = 1, 4
344         DO jind2 = 1, 4
345            DO jfl = 1, jpnfl
346               IF( jind1 /= jind2 ) THEN
347                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))) )
348                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
349                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))) )
350               ENDIF
351            END DO
352         END DO
353      END DO
354     
355      ! velocity when we compute at middle time step
356     
357      DO jfl = 1, jpnfl
358         DO jind1 = 1, 4
359            DO jind2 = 1, 4
360               DO jind3 = 1 ,4
361                  ztvfl(jfl,jind1,jind2,jind3)=   &
362                     &   ( tcoef1(kind) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3))  +   &
363                     &     tcoef2(kind) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) )    & 
364                     &     / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
365               END DO
366            END DO
367         END DO
368         
369         zsumv=0.
370         DO jind1 = 1, 4
371            DO jind2 = 1, 4
372               DO jind3 = 1, 4
373                  zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2)   &
374                     &  * zlagzv(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
375               END DO
376            END DO
377         END DO
378         pvfl(jfl) = zsumv
379      END DO
380     
381      ! Interpolation of W velocity
382
383      ! nearest neightboring point for computation of w
384      DO jfl = 1, jpnfl
385         iilw(jfl) = INT(pxt(jfl))
386         ijlw(jfl) = INT(pyt(jfl))
387         iklw(jfl) = INT(pzt(jfl)+.5)
388      END DO
389     
390      ! 64 neightboring points for computation of w
391      DO jind1 = 1, 4
392         DO jfl = 1, jpnfl
393            ! i-direction
394            IF( iilw(jfl) <= 2 ) THEN
395               iidw(jfl,jind1) = jind1
396            ELSE
397               IF( iilw(jfl) >= jpi-1 ) THEN
398                  iidw(jfl,jind1) = jpi + jind1 - 4
399               ELSE
400                  iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
401               ENDIF
402            ENDIF
403            ! j-direction
404            IF( ijlw(jfl) <= 2 ) THEN
405               ijdw(jfl,jind1) = jind1
406            ELSE
407               IF( ijlw(jfl) >= jpj-1 ) THEN
408                  ijdw(jfl,jind1) = jpj + jind1 - 4
409               ELSE
410                  ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
411               ENDIF
412            ENDIF
413            ! k-direction
414            IF( iklw(jfl) <= 2 ) THEN
415               ikdw(jfl,jind1) = jind1
416            ELSE
417               IF( iklw(jfl) >= jpk-1 ) THEN
418                  ikdw(jfl,jind1) = jpk + jind1 - 4
419               ELSE
420                  ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
421               ENDIF
422            ENDIF
423         END DO
424      END DO
425      DO jind1 = 1, 4
426         DO jfl = 1, jpnfl
427            IF( iklw(jfl) <= 2 ) THEN
428               ikdw(jfl,jind1) = jind1
429            ELSE
430               IF( iklw(jfl) >= jpk-1 ) THEN
431                  ikdw(jfl,jind1) = jpk + jind1 - 4
432               ELSE
433                  ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
434               ENDIF
435            ENDIF
436         END DO
437      END DO
438     
439      ! Lagrange coefficients  for w interpolation
440     
441      DO jfl = 1, jpnfl
442         DO jind1 = 1, 4
443            zlagxw(jfl,jind1) = 1.
444            zlagyw(jfl,jind1) = 1.
445            zlagzw(jfl,jind1) = 1.
446         END DO
447      END DO
448       
449      DO jind1 = 1, 4
450         DO jind2 = 1, 4
451            DO jfl = 1, jpnfl
452               IF( jind1 /= jind2 ) THEN
453                  zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2))) )
454                  zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2))) )
455                  zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
456               ENDIF
457            END DO
458         END DO
459      END DO
460     
461      ! velocity w  when we compute at middle time step
462     
463      DO jfl = 1, jpnfl
464         DO jind1 = 1, 4
465            DO jind2 = 1, 4
466               DO jind3 = 1, 4
467                  ztwfl(jfl,jind1,jind2,jind3)=   &
468                     &   ( tcoef1(kind) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
469                     &     tcoef2(kind) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) )  &
470!!bug e3w instead of fse3
471                     &   / e3w(ikdw(jfl,jind3))
472               END DO
473            END DO
474         END DO
475         
476         zsumw=0.
477         DO jind1 = 1, 4
478            DO jind2 = 1, 4
479               DO jind3 = 1, 4
480                  zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2)   &
481                     &  * zlagzw(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
482               END DO
483            END DO
484         END DO
485         pwfl(jfl) = zsumw
486      END DO
487     
488   END SUBROUTINE flo_interp
489
490#  else
491   !!----------------------------------------------------------------------
492   !!   Default option                                         Empty module
493   !!----------------------------------------------------------------------
494CONTAINS
495   SUBROUTINE flo_4rk                 ! Empty routine
496   END SUBROUTINE flo_4rk
497#endif
498   
499   !!======================================================================
500END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.