source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 18.7 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
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   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
32   !! $Id$
33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE flo_4rk( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE flo_4rk  ***
40      !!
41      !!  ** Purpose :   Compute the geographical position (lat,lon,depth)
42      !!       of each float at each time step.
43      !!
44      !!  ** Method  :   The position of a float is computed with a 4th order
45      !!       Runge-Kutta scheme and and Lagrange interpolation.
46      !!         We need to know the velocity field, the old positions of the
47      !!       floats and the grid defined on the domain.
48      !!----------------------------------------------------------------------
49      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
50      !!
51      INTEGER ::  jfl, jind           ! dummy loop indices
52      INTEGER ::  ierror              ! error value
53
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( ierror /= 0 ) THEN
60         WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
61      ENDIF
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         DO jfl = 1, jpnfl
133            zrkxfl(jfl,jind) = rdt*zufl(jfl)
134            zrkyfl(jfl,jind) = rdt*zvfl(jfl)
135            zrkzfl(jfl,jind) = rdt*zwfl(jfl)
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
152      !
153      !
154   END SUBROUTINE flo_4rk
155
156
157   SUBROUTINE flo_interp( pxt , pyt , pzt ,      &
158      &                   pufl, pvfl, pwfl, ki )
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      REAL(wp) , DIMENSION(jpnfl), INTENT(in   ) ::   pxt , pyt , pzt    ! position of the float
170      REAL(wp) , DIMENSION(jpnfl), INTENT(  out) ::   pufl, pvfl, pwfl   ! velocity at this position
171      INTEGER                    , INTENT(in   ) ::   ki                 !
172      !!
173      INTEGER  ::   jfl, jind1, jind2, jind3   ! dummy loop indices
174      REAL(wp) ::   zsumu, zsumv, zsumw        ! local scalar
175      INTEGER  , DIMENSION(jpnfl)       ::   iilu, ijlu, iklu   ! nearest neighbour INDEX-u
176      INTEGER  , DIMENSION(jpnfl)       ::   iilv, ijlv, iklv   ! nearest neighbour INDEX-v
177      INTEGER  , DIMENSION(jpnfl)       ::   iilw, ijlw, iklw   ! nearest neighbour INDEX-w
178      INTEGER  , DIMENSION(jpnfl,4)     ::   iidu, ijdu, ikdu   ! 64 nearest neighbour INDEX-u
179      INTEGER  , DIMENSION(jpnfl,4)     ::   iidv, ijdv, ikdv   ! 64 nearest neighbour INDEX-v
180      INTEGER  , DIMENSION(jpnfl,4)     ::   iidw, ijdw, ikdw   ! 64 nearest neighbour INDEX-w
181      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxu, zlagyu, zlagzu   ! Lagrange  coefficients
182      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxv, zlagyv, zlagzv   !    -           -
183      REAL(wp) , DIMENSION(jpnfl,4)     ::   zlagxw, zlagyw, zlagzw   !    -           -
184      REAL(wp) , DIMENSION(jpnfl,4,4,4) ::   ztufl , ztvfl , ztwfl   ! velocity at choosen time step
185      !!---------------------------------------------------------------------
186 
187      ! Interpolation of U velocity
188
189      ! nearest neightboring point for computation of u       
190      DO jfl = 1, jpnfl
191         iilu(jfl) = INT(pxt(jfl)-.5)
192         ijlu(jfl) = INT(pyt(jfl)-.5)
193         iklu(jfl) = INT(pzt(jfl))
194      END DO
195     
196      !  64 neightboring points for computation of u
197      DO jind1 = 1, 4
198         DO jfl = 1, jpnfl
199            !  i-direction
200            IF( iilu(jfl) <= 2 ) THEN          ;   iidu(jfl,jind1) = jind1
201            ELSE
202               IF( iilu(jfl) >= jpi-1 ) THEN   ;   iidu(jfl,jind1) = jpi       + jind1 - 4
203               ELSE                            ;   iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
204               ENDIF
205            ENDIF
206            !  j-direction
207            IF( ijlu(jfl) <= 2 ) THEN          ;   ijdu(jfl,jind1) = jind1
208            ELSE
209               IF( ijlu(jfl) >= jpj-1 ) THEN   ;   ijdu(jfl,jind1) = jpj       + jind1 - 4
210               ELSE                            ;   ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
211               ENDIF
212            ENDIF
213            ! k-direction
214            IF( iklu(jfl) <= 2 ) THEN          ;   ikdu(jfl,jind1) = jind1
215            ELSE
216               IF( iklu(jfl) >= jpk-1 ) THEN   ;   ikdu(jfl,jind1) = jpk + jind1 - 4
217               ELSE                            ;   ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
218               ENDIF
219            ENDIF
220         END DO
221      END DO
222     
223      ! Lagrange coefficients
224      DO jfl = 1, jpnfl
225         DO jind1 = 1, 4
226            zlagxu(jfl,jind1) = 1.
227            zlagyu(jfl,jind1) = 1.
228            zlagzu(jfl,jind1) = 1.
229         END DO
230      END DO
231      DO jind1 = 1, 4
232         DO jind2 = 1, 4
233            DO jfl= 1, jpnfl
234               IF( jind1 /= jind2 ) THEN
235                  zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
236                  zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
237                  zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
238               ENDIF
239            END DO
240         END DO
241      END DO
242     
243      ! velocity when we compute at middle time step
244     
245      DO jfl = 1, jpnfl
246         DO jind1 = 1, 4
247            DO jind2 = 1, 4
248               DO jind3 = 1, 4
249                  ztufl(jfl,jind1,jind2,jind3) =   &
250                     &   (  tcoef1(ki) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) +   &
251                     &      tcoef2(ki) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) )   &
252                     &      / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 
253               END DO
254            END DO
255         END DO
256         
257         zsumu = 0.
258         DO jind1 = 1, 4
259            DO jind2 = 1, 4
260               DO jind3 = 1, 4
261                  zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2)   &
262                     &  * zlagzu(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3) 
263               END DO
264            END DO
265         END DO
266         pufl(jfl) = zsumu
267      END DO
268     
269      ! Interpolation of V velocity
270
271      ! nearest neightboring point for computation of v
272      DO jfl = 1, jpnfl
273         iilv(jfl) = INT(pxt(jfl)-.5)
274         ijlv(jfl) = INT(pyt(jfl)-.5)
275         iklv(jfl) = INT(pzt(jfl))
276      END DO
277     
278      ! 64 neightboring points for computation of v
279      DO jind1 = 1, 4
280         DO jfl = 1, jpnfl
281            ! i-direction
282            IF( iilv(jfl) <= 2 ) THEN          ;   iidv(jfl,jind1) = jind1
283            ELSE
284               IF( iilv(jfl) >= jpi-1 ) THEN   ;   iidv(jfl,jind1) = jpi       + jind1 - 4
285               ELSE                            ;   iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
286               ENDIF
287            ENDIF
288            ! j-direction
289            IF( ijlv(jfl) <= 2 ) THEN          ;   ijdv(jfl,jind1) = jind1
290            ELSE
291               IF( ijlv(jfl) >= jpj-1 ) THEN   ;   ijdv(jfl,jind1) = jpj       + jind1 - 4
292               ELSE                            ;   ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
293               ENDIF
294            ENDIF
295            ! k-direction
296            IF( iklv(jfl) <= 2 ) THEN          ;   ikdv(jfl,jind1) = jind1
297            ELSE
298               IF( iklv(jfl) >= jpk-1 ) THEN   ;   ikdv(jfl,jind1) = jpk + jind1 - 4
299               ELSE                            ;   ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
300               ENDIF
301            ENDIF
302         END DO
303      END DO
304     
305      ! Lagrange coefficients
306     
307      DO jfl = 1, jpnfl
308         DO jind1 = 1, 4
309            zlagxv(jfl,jind1) = 1.
310            zlagyv(jfl,jind1) = 1.
311            zlagzv(jfl,jind1) = 1.
312         END DO
313      END DO
314     
315      DO jind1 = 1, 4
316         DO jind2 = 1, 4
317            DO jfl = 1, jpnfl
318               IF( jind1 /= jind2 ) THEN
319                  zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2))   ) )
320                  zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
321                  zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2))   ) )
322               ENDIF
323            END DO
324         END DO
325      END DO
326     
327      ! velocity when we compute at middle time step
328     
329      DO jfl = 1, jpnfl
330         DO jind1 = 1, 4
331            DO jind2 = 1, 4
332               DO jind3 = 1 ,4
333                  ztvfl(jfl,jind1,jind2,jind3)=   &
334                     &   ( tcoef1(ki) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3))  +   &
335                     &     tcoef2(ki) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) )    & 
336                     &     / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
337               END DO
338            END DO
339         END DO
340         
341         zsumv=0.
342         DO jind1 = 1, 4
343            DO jind2 = 1, 4
344               DO jind3 = 1, 4
345                  zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2)   &
346                     &  * zlagzv(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
347               END DO
348            END DO
349         END DO
350         pvfl(jfl) = zsumv
351      END DO
352     
353      ! Interpolation of W velocity
354
355      ! nearest neightboring point for computation of w
356      DO jfl = 1, jpnfl
357         iilw(jfl) = INT( pxt(jfl)   )
358         ijlw(jfl) = INT( pyt(jfl)   )
359         iklw(jfl) = INT( pzt(jfl)+.5)
360      END DO
361     
362      ! 64 neightboring points for computation of w
363      DO jind1 = 1, 4
364         DO jfl = 1, jpnfl
365            ! i-direction
366            IF( iilw(jfl) <= 2 ) THEN          ;   iidw(jfl,jind1) = jind1
367            ELSE
368               IF( iilw(jfl) >= jpi-1 ) THEN   ;   iidw(jfl,jind1) = jpi       + jind1 - 4
369               ELSE                            ;   iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
370               ENDIF
371            ENDIF
372            ! j-direction
373            IF( ijlw(jfl) <= 2 ) THEN          ;   ijdw(jfl,jind1) = jind1
374            ELSE
375               IF( ijlw(jfl) >= jpj-1 ) THEN   ;   ijdw(jfl,jind1) = jpj       + jind1 - 4
376               ELSE                            ;   ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
377               ENDIF
378            ENDIF
379            ! k-direction
380            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
381            ELSE
382               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
383               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
384               ENDIF
385            ENDIF
386         END DO
387      END DO
388      DO jind1 = 1, 4
389         DO jfl = 1, jpnfl
390            IF( iklw(jfl) <= 2 ) THEN          ;   ikdw(jfl,jind1) = jind1
391            ELSE
392               IF( iklw(jfl) >= jpk-1 ) THEN   ;   ikdw(jfl,jind1) = jpk       + jind1 - 4
393               ELSE                            ;   ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
394               ENDIF
395            ENDIF
396         END DO
397      END DO
398     
399      ! Lagrange coefficients  for w interpolation
400      DO jfl = 1, jpnfl
401         DO jind1 = 1, 4
402            zlagxw(jfl,jind1) = 1.
403            zlagyw(jfl,jind1) = 1.
404            zlagzw(jfl,jind1) = 1.
405         END DO
406      END DO
407      DO jind1 = 1, 4
408         DO jind2 = 1, 4
409            DO jfl = 1, jpnfl
410               IF( jind1 /= jind2 ) THEN
411                  zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2))   ) )
412                  zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2))   ) )
413                  zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
414               ENDIF
415            END DO
416         END DO
417      END DO
418     
419      ! velocity w  when we compute at middle time step
420      DO jfl = 1, jpnfl
421         DO jind1 = 1, 4
422            DO jind2 = 1, 4
423               DO jind3 = 1, 4
424                  ztwfl(jfl,jind1,jind2,jind3)=   &
425                     &   ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+   &
426                     &     tcoef2(ki) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) )  &
427                     &   / e3w_n(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))
428               END DO
429            END DO
430         END DO
431         
432         zsumw = 0.e0
433         DO jind1 = 1, 4
434            DO jind2 = 1, 4
435               DO jind3 = 1, 4
436                  zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2)   &
437                     &  * zlagzw(jfl,jind3) * rcoef(jind1)*rcoef(jind2)*rcoef(jind3)
438               END DO
439            END DO
440         END DO
441         pwfl(jfl) = zsumw
442      END DO
443      !   
444      !   
445   END SUBROUTINE flo_interp
446
447#  else
448   !!----------------------------------------------------------------------
449   !!   No floats                                              Dummy module
450   !!----------------------------------------------------------------------
451#endif
452   
453   !!======================================================================
454END MODULE flo4rk
Note: See TracBrowser for help on using the repository browser.