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

Last change on this file since 11536 was 11536, checked in by smasson, 20 months ago

trunk: merge dev_r10984_HPC-13 into the trunk

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