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

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

Initial revision
  • 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
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Accessibility
22   PUBLIC flo_4rk      ! routine called by floats.F90
23
24   !! * Module variables
25   REAL(wp), DIMENSION (4) ::   &   ! RK4 and Lagrange interpolation
26      tcoef1 = /  1.0  ,  0.5  ,  0.5  ,  0.0  / ,  &  ! coeffients  for
27      tcoef2 = /  0.0  ,  0.5  ,  0.5  ,  1.0  / ,  &  ! lagrangian interp.
28      scoef2 = /  1.0  ,  2.0  ,  2.0  ,  1.0  / ,  &  ! RK4 coefficients
29      rcoef  = /-1./6. , 1./2. ,-1./2. , 1./6. /       ! ???
30   REAL(wp), DIMENSION (3) ::   &
31      scoef1 = / .5, .5, 1. /       ! compute position with interpolated
32   !!----------------------------------------------------------------------
33   !!   OPA 9.0 , LODYC-IPSL  (2003)
34   !!----------------------------------------------------------------------
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      !!----------------------------------------------------------------------
51      !! * Arguments
52      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
53
54      !! * Local declarations
55      INTEGER ::  jfl, jind           ! dummy loop indices
56      REAL(wp), DIMENSION ( jpnfl)  ::   &
57         zgifl, zgjfl, zgkfl,   &     ! index RK  positions
58         zufl, zvfl, zwfl             ! interpolated velocity at the
59      ! float position
60      REAL(wp), DIMENSION ( jpnfl, 4 )  ::    &
61         zrkxfl, zrkyfl, zrkzfl       ! RK coefficients
62      !!---------------------------------------------------------------------
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       
126      DO  jind = 1, 4
127         
128         ! for each step we compute the compute the velocity with Lagrange interpolation
129         CALL flo_interp(zgifl,zgjfl,zgkfl,zufl,zvfl,zwfl,jind)
130         
131         ! computation of Runge-Kutta factor
132         
133         DO jfl = 1, jpnfl
134            zrkxfl(jfl,jind) = rdt*zufl(jfl)
135            zrkyfl(jfl,jind) = rdt*zvfl(jfl)
136            zrkzfl(jfl,jind) = rdt*zwfl(jfl)
137         END DO
138         IF( jind /= 4 ) THEN
139            DO jfl = 1, jpnfl
140               zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
141               zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
142               zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
143            END DO
144         ENDIF
145      END DO
146      DO jind = 1, 4
147         DO jfl = 1, jpnfl
148            tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
149            tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
150            tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
151         END DO
152      END DO
153
154   END SUBROUTINE flo_4rk
155
156
157   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
158      &                   pufl, pvfl, pwfl, kind )
159      !!----------------------------------------------------------------------
160      !!                ***  ROUTINE flointerp  ***
161      !!
162      !! ** Purpose :   Interpolation of the velocity on the float position
163      !!
164      !! ** Method  :   Lagrange interpolation with the 64 neighboring
165      !!      points. This routine is call 4 time at each time step to
166      !!      compute velocity at the date and the position we need to
167      !!      integrated with RK method.
168      !!
169      !!----------------------------------------------------------------------
170      !! * Local declarations
171      INTEGER ::   &
172         kind,   &
173         jfl, jind1, jind2, jind3,   &
174         zsumu, zsumv, zsumw
175      INTEGER , DIMENSION  ( jpnfl  ) ::   &
176         iilu, ijlu, iklu,     &           ! nearest neighbour INDEX-u
177         iilv, ijlv, iklv,     &           ! nearest neighbour INDEX-v
178         iilw, ijlw, iklw                  ! nearest neighbour INDEX-w
179      INTEGER , DIMENSION  ( jpnfl, 4  ) ::   &
180         iidu, ijdu, ikdu,    &            ! 64 nearest neighbour INDEX-u
181         iidv, ijdv, ikdv,    &            ! 64 nearest neighbour INDEX-v
182         iidw, ijdw, ikdw                  ! 64 nearest neighbour INDEX-w
183      REAL(wp) , DIMENSION  ( jpnfl  ) ::   &
184         pxt , pyt , pzt,     &            ! position of the float
185         pufl, pvfl, pwfl                  ! velocity at this position
186      REAL(wp) , DIMENSION  ( jpnfl, 4, 4, 4 ) ::   &
187         ztufl, ztvfl, ztwfl               ! velocity at choosen time step
188      REAL(wp) , DIMENSION  ( jpnfl, 4 ) ::   &
189         zlagxu, zlagyu, zlagzu,   &       ! Lagrange  coefficients
190         zlagxv, zlagyv, zlagzv,   &
191         zlagxw, zlagyw, zlagzw
192      !!---------------------------------------------------------------------
193     
194      ! Interpolation of U velocity
195
196      ! nearest neightboring point for computation of u       
197      DO jfl = 1, jpnfl
198         iilu(jfl) = INT(pxt(jfl)-.5)
199         ijlu(jfl) = INT(pyt(jfl)-.5)
200         iklu(jfl) = INT(pzt(jfl))
201      END DO
202     
203      !  64 neightboring points for computation of u
204      DO jind1 = 1, 4
205         DO jfl = 1, jpnfl
206            !  i-direction
207            IF( iilu(jfl) <= 2 ) THEN
208               iidu(jfl,jind1) = jind1
209            ELSE
210               IF( iilu(jfl) >= jpi-1 ) THEN
211                  iidu(jfl,jind1) = jpi + jind1 - 4
212               ELSE
213                  iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
214               ENDIF
215            ENDIF
216            !  j-direction
217            IF( ijlu(jfl) <= 2 ) THEN
218               ijdu(jfl,jind1) = jind1
219            ELSE
220               IF( ijlu(jfl) >= jpj-1 ) THEN
221                  ijdu(jfl,jind1) = jpj + jind1 - 4
222               ELSE
223                  ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
224               ENDIF
225            ENDIF
226            ! k-direction
227            IF( iklu(jfl) <= 2 ) THEN
228               ikdu(jfl,jind1) = jind1
229            ELSE
230               IF( iklu(jfl) >= jpk-1 ) THEN
231                  ikdu(jfl,jind1) = jpk + jind1 - 4
232               ELSE
233                  ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
234               ENDIF
235            ENDIF
236         END DO
237      END DO
238     
239      ! Lagrange coefficients
240     
241      DO jfl = 1, jpnfl
242         DO jind1 = 1, 4
243            zlagxu(jfl,jind1) = 1.
244            zlagyu(jfl,jind1) = 1.
245            zlagzu(jfl,jind1) = 1.
246         END DO
247      END DO
248     
249      DO jind1 = 1, 4
250         DO jind2 = 1, 4
251            DO jfl= 1, jpnfl
252               IF( jind1 /= jind2 ) THEN
253                  zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
254                  zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
255                  zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
256               ENDIF
257            END DO
258         END DO
259      END DO
260     
261      ! velocity when we compute at middle time step
262     
263      DO jfl = 1, jpnfl
264         DO jind1 = 1, 4
265            DO jind2 = 1, 4
266               DO jind3 = 1, 4
267                  ztufl(jfl,jind1,jind2,jind3) =   &
268                     &   (  tcoef1(kind) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) +   &
269                     &      tcoef2(kind) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) )   &
270                     &      / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 
271               END DO
272            END DO
273         END DO
274         
275         zsumu = 0.
276         DO jind1 = 1, 4
277            DO jind2 = 1, 4
278               DO jind3 = 1, 4
279                  zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2)   &
280                     &  * zlagzu(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3) 
281               END DO
282            END DO
283         END DO
284         pufl(jfl) = zsumu
285      END DO
286     
287      ! Interpolation of V velocity
288
289      ! nearest neightboring point for computation of v
290      DO jfl = 1, jpnfl
291         iilv(jfl) = INT(pxt(jfl)-.5)
292         ijlv(jfl) = INT(pyt(jfl)-.5)
293         iklv(jfl) = INT(pzt(jfl))
294      END DO
295     
296      ! 64 neightboring points for computation of v
297      DO jind1 = 1, 4
298         DO jfl = 1, jpnfl
299            ! i-direction
300            IF( iilv(jfl) <= 2 ) THEN
301               iidv(jfl,jind1) = jind1
302            ELSE
303               IF( iilv(jfl) >= jpi-1 ) THEN
304                  iidv(jfl,jind1) = jpi + jind1 - 4
305               ELSE
306                  iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
307               ENDIF
308            ENDIF
309            ! j-direction
310            IF( ijlv(jfl) <= 2 ) THEN
311               ijdv(jfl,jind1) = jind1
312            ELSE
313               IF( ijlv(jfl) >= jpj-1 ) THEN
314                  ijdv(jfl,jind1) = jpj + jind1 - 4
315               ELSE
316                  ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
317               ENDIF
318            ENDIF
319            ! k-direction
320            IF( iklv(jfl) <= 2 ) THEN
321               ikdv(jfl,jind1) = jind1
322            ELSE
323               IF( iklv(jfl) >= jpk-1 ) THEN
324                  ikdv(jfl,jind1) = jpk + jind1 - 4
325               ELSE
326                  ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
327               ENDIF
328            ENDIF
329         END DO
330      END DO
331     
332      ! Lagrange coefficients
333     
334      DO jfl = 1, jpnfl
335         DO jind1 = 1, 4
336            zlagxv(jfl,jind1) = 1.
337            zlagyv(jfl,jind1) = 1.
338            zlagzv(jfl,jind1) = 1.
339         END DO
340      END DO
341     
342      DO jind1 = 1, 4
343         DO jind2 = 1, 4
344            DO jfl = 1, jpnfl
345               IF( jind1 /= jind2 ) THEN
346                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))) )
347                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
348                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))) )
349               ENDIF
350            END DO
351         END DO
352      END DO
353     
354      ! velocity when we compute at middle time step
355     
356      DO jfl = 1, jpnfl
357         DO jind1 = 1, 4
358            DO jind2 = 1, 4
359               DO jind3 = 1 ,4
360                  ztvfl(jfl,jind1,jind2,jind3)=   &
361                     &   ( tcoef1(kind) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3))  +   &
362                     &     tcoef2(kind) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) )    & 
363                     &     / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
364               END DO
365            END DO
366         END DO
367         
368         zsumv=0.
369         DO jind1 = 1, 4
370            DO jind2 = 1, 4
371               DO jind3 = 1, 4
372                  zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2)   &
373                     &  * zlagzv(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
374               END DO
375            END DO
376         END DO
377         pvfl(jfl) = zsumv
378      END DO
379     
380      ! Interpolation of W velocity
381
382      ! nearest neightboring point for computation of w
383      DO jfl = 1, jpnfl
384         iilw(jfl) = INT(pxt(jfl))
385         ijlw(jfl) = INT(pyt(jfl))
386         iklw(jfl) = INT(pzt(jfl)+.5)
387      END DO
388     
389      ! 64 neightboring points for computation of w
390      DO jind1 = 1, 4
391         DO jfl = 1, jpnfl
392            ! i-direction
393            IF( iilw(jfl) <= 2 ) THEN
394               iidw(jfl,jind1) = jind1
395            ELSE
396               IF( iilw(jfl) >= jpi-1 ) THEN
397                  iidw(jfl,jind1) = jpi + jind1 - 4
398               ELSE
399                  iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
400               ENDIF
401            ENDIF
402            ! j-direction
403            IF( ijlw(jfl) <= 2 ) THEN
404               ijdw(jfl,jind1) = jind1
405            ELSE
406               IF( ijlw(jfl) >= jpj-1 ) THEN
407                  ijdw(jfl,jind1) = jpj + jind1 - 4
408               ELSE
409                  ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
410               ENDIF
411            ENDIF
412            ! k-direction
413            IF( iklw(jfl) <= 2 ) THEN
414               ikdw(jfl,jind1) = jind1
415            ELSE
416               IF( iklw(jfl) >= jpk-1 ) THEN
417                  ikdw(jfl,jind1) = jpk + jind1 - 4
418               ELSE
419                  ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
420               ENDIF
421            ENDIF
422         END DO
423      END DO
424      DO jind1 = 1, 4
425         DO jfl = 1, jpnfl
426            IF( iklw(jfl) <= 2 ) THEN
427               ikdw(jfl,jind1) = jind1
428            ELSE
429               IF( iklw(jfl) >= jpk-1 ) THEN
430                  ikdw(jfl,jind1) = jpk + jind1 - 4
431               ELSE
432                  ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
433               ENDIF
434            ENDIF
435         END DO
436      END DO
437     
438      ! Lagrange coefficients  for w interpolation
439     
440      DO jfl = 1, jpnfl
441         DO jind1 = 1, 4
442            zlagxw(jfl,jind1) = 1.
443            zlagyw(jfl,jind1) = 1.
444            zlagzw(jfl,jind1) = 1.
445         END DO
446      END DO
447       
448      DO jind1 = 1, 4
449         DO jind2 = 1, 4
450            DO jfl = 1, jpnfl
451               IF( jind1 /= jind2 ) THEN
452                  zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2))) )
453                  zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2))) )
454                  zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
455               ENDIF
456            END DO
457         END DO
458      END DO
459     
460      ! velocity w  when we compute at middle time step
461     
462      DO jfl = 1, jpnfl
463         DO jind1 = 1, 4
464            DO jind2 = 1, 4
465               DO jind3 = 1, 4
466                  ztwfl(jfl,jind1,jind2,jind3)=   &
467                     &   ( tcoef1(kind) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
468                     &     tcoef2(kind) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) )  &
469!!bug e3w instead of fse3
470                     &   / e3w(ikdw(jfl,jind3))
471               END DO
472            END DO
473         END DO
474         
475         zsumw=0.
476         DO jind1 = 1, 4
477            DO jind2 = 1, 4
478               DO jind3 = 1, 4
479                  zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2)   &
480                     &  * zlagzw(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
481               END DO
482            END DO
483         END DO
484         pwfl(jfl) = zsumw
485      END DO
486     
487   END SUBROUTINE flo_interp
488
489#  else
490   !!----------------------------------------------------------------------
491   !!   Default option                                         Empty module
492   !!----------------------------------------------------------------------
493CONTAINS
494   SUBROUTINE flo_4rk                 ! Empty routine
495   END SUBROUTINE flo_4rk
496#endif
497   
498   !!======================================================================
499END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.