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

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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