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.
floblk.F90 in trunk/NEMO/OPA_SRC/FLO – NEMO

source: trunk/NEMO/OPA_SRC/FLO/floblk.F90 @ 84

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

CT : UPDATE057 : # General syntax, alignement, comments corrections

# l_ctl alone replace the set (l_ctl .AND. lwp)
# Add of diagnostics which are activated when using l_ctl logical
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.1 KB
Line 
1MODULE floblk
2   !!======================================================================
3   !!                     ***  MODULE  floblk  ***
4   !! Ocean floats :   trajectory computation
5   !!======================================================================
6#if   defined key_floats   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_floats'                                     float trajectories
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!    flotblk     : compute float trajectories with Blanke algorithme
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE flo_oce         ! ocean drifting floats
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distribued memory computing library
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Accessibility
26   PUBLIC flo_blk      ! routine called by floats.F90
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LODYC-IPSL  (2003)
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE flo_blk( kt )
34      !!---------------------------------------------------------------------
35      !!                  ***  ROUTINE flo_blk  ***
36      !!           
37      !! ** Purpose :   Compute the geographical position,latitude, longitude
38      !!      and depth of each float at each time step.
39      !!
40      !! ** Method  :   The position of a float is computed with Bruno Blanke
41      !!      algorithm. We need to know the velocity field, the old positions
42      !!      of the floats and the grid defined on the domain.
43      !!
44      !!----------------------------------------------------------------------
45      !! * arguments
46      INTEGER, INTENT( in  ) ::   kt ! ocean time step
47
48      !! * Local declarations
49      INTEGER :: jfl              ! dummy loop arguments
50      INTEGER :: ind, ifin, iloop
51      INTEGER , DIMENSION ( jpnfl )  ::   &
52         iil, ijl, ikl,             &     ! index of nearest mesh
53         iiloc , ijloc,             &
54         iiinfl, ijinfl, ikinfl,    &     ! index of input mesh of the float.
55         iioutfl, ijoutfl, ikoutfl        ! index of output mesh of the float.
56      REAL(wp) , DIMENSION ( jpnfl )  ::    &
57         zgifl, zgjfl, zgkfl,       &     ! position of floats, index on
58                                          ! velocity mesh.
59         ztxfl, ztyfl, ztzfl,       &     ! time for a float to quit the mesh
60                                          ! across one of the face x,y and z
61         zttfl,                     &     ! time for a float to quit the mesh
62         zagefl,                    &     ! time during which, trajectorie of
63                                          ! the float has been computed
64         zagenewfl,                 &     ! new age of float after calculation
65                                          ! of new position
66         zufl, zvfl, zwfl,          &     ! interpolated vel. at float position
67         zudfl, zvdfl, zwdfl,       &     ! velocity diff input/output of mesh
68         zgidfl, zgjdfl, zgkdfl           ! direction index of float
69      REAL(wp)   ::       &
70         zuinfl,zvinfl,zwinfl,      &     ! transport across the input face
71         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
72         zvol,                      &     ! volume of the mesh
73         zsurfz,                    &     ! surface of the face of the mesh
74         zind
75      REAL(wp), DIMENSION ( 2 )  ::   &
76         zsurfx, zsurfy                   ! surface of the face of the mesh
77      !!---------------------------------------------------------------------
78     
79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'flo_blk : compute Blanke trajectories for floats '
82         IF(lwp) WRITE(numout,*) '~~~~~~~ '
83      ENDIF
84
85      ! Initialisation of parameters
86     
87      DO jfl = 1, jpnfl
88         ! ages of floats are put at zero
89         zagefl(jfl) = 0.
90         ! index on the velocity grid
91         ! We considere k coordinate negative, with this transformation
92         ! the computation in the 3 direction is the same.
93         zgifl(jfl) = tpifl(jfl) - 0.5
94         zgjfl(jfl) = tpjfl(jfl) - 0.5
95         zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl)))
96         ! surface drift every 10 days
97         IF( ln_argo ) THEN
98            IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 )  zgkfl(jfl) = -1.
99         ENDIF
100         ! index of T mesh
101         iil(jfl) = 1 + INT(zgifl(jfl))
102         ijl(jfl) = 1 + INT(zgjfl(jfl))
103         ikl(jfl) =     INT(zgkfl(jfl))
104      END DO
105       
106      iloop = 0
107222   DO jfl = 1, jpnfl
108# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
109         IF( (iil(jfl) >= (mig(nldi)-jpizoom+1)) .AND. (iil(jfl) <= (mig(nlei)-jpizoom+1)) .AND.   &
110             (ijl(jfl) >= (mjg(nldj)-jpjzoom+1)) .AND. (ijl(jfl) <= (mjg(nlej)-jpjzoom+1)) ) THEN
111            iiloc(jfl) = iil(jfl) - (mig(1)-jpizoom+1) + 1
112            ijloc(jfl) = ijl(jfl) - (mjg(1)-jpjzoom+1) + 1
113# else
114            iiloc(jfl) = iil(jfl)
115            ijloc(jfl) = ijl(jfl)
116# endif
117           
118            ! compute the transport across the mesh where the float is.           
119            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * e3t(-ikl(jfl))
120            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * e3t(-ikl(jfl))
121            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * e3t(-ikl(jfl))
122            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * e3t(-ikl(jfl))
123
124            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
125            zsurfz=  e1t(iiloc(jfl),ijloc(jfl)) * e2t(iiloc(jfl),ijloc(jfl))
126            zvol  =( e1t(iiloc(jfl),ijloc(jfl)) * e2t(iiloc(jfl),ijloc(jfl)) * e3t(-ikl(jfl)) )
127
128            !
129            zuinfl =( ub(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(1)
130            zuoutfl=( ub(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(2)
131            zvinfl =( vb(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) )/2.*zsurfy(1)
132            zvoutfl=( vb(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl)) )/2.*zsurfy(2)
133            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
134               &   +  wn(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
135            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
136               &   +  wn(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
137           
138            ! interpolation of velocity field on the float initial position           
139            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
140            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
141            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
142           
143            ! faces of input and output
144            ! u-direction
145            IF( zufl(jfl) < 0. ) THEN
146               iioutfl(jfl) = iil(jfl) - 1.
147               iiinfl (jfl) = iil(jfl)
148               zind   = zuinfl
149               zuinfl = zuoutfl
150               zuoutfl= zind
151            ELSE
152               iioutfl(jfl) = iil(jfl)
153               iiinfl (jfl) = iil(jfl) - 1
154            ENDIF
155            ! v-direction       
156            IF( zvfl(jfl) < 0. ) THEN
157               ijoutfl(jfl) = ijl(jfl) - 1.
158               ijinfl (jfl) = ijl(jfl)
159               zind    = zvinfl
160               zvinfl  = zvoutfl
161               zvoutfl = zind
162            ELSE
163               ijoutfl(jfl) = ijl(jfl)
164               ijinfl (jfl) = ijl(jfl) - 1.
165            ENDIF
166            ! w-direction
167            IF( zwfl(jfl) < 0. ) THEN
168               ikoutfl(jfl) = ikl(jfl) - 1.
169               ikinfl (jfl) = ikl(jfl)
170               zind    = zwinfl
171               zwinfl  = zwoutfl
172               zwoutfl = zind
173            ELSE
174               ikoutfl(jfl) = ikl(jfl)
175               ikinfl (jfl) = ikl(jfl) - 1.
176            ENDIF
177           
178            ! compute the time to go out the mesh across a face
179            ! u-direction
180            zudfl (jfl) = zuoutfl - zuinfl
181            zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl))
182            IF( zufl(jfl)*zuoutfl <= 0. ) THEN
183               ztxfl(jfl) = 1.E99
184            ELSE
185               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
186                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
187               ELSE
188                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
189               ENDIF
190               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
191                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
192                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
193               ENDIF
194            ENDIF
195            ! v-direction
196            zvdfl (jfl) = zvoutfl - zvinfl
197            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
198            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
199               ztyfl(jfl) = 1.E99
200            ELSE
201               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
202                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
203               ELSE
204                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
205               ENDIF
206               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
207                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
208                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
209               ENDIF
210            ENDIF
211            ! w-direction       
212            IF( nisobfl(jfl) == 1. ) THEN
213               zwdfl (jfl) = zwoutfl - zwinfl
214               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
215               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
216                  ztzfl(jfl) = 1.E99
217               ELSE
218                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
219                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
220                  ELSE
221                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
222                  ENDIF
223                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
224                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
225                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
226                  ENDIF
227               ENDIF
228            ENDIF
229           
230            ! the time to go leave the mesh is the smallest time
231                   
232            IF( nisobfl(jfl) == 1. ) THEN
233               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
234            ELSE
235               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
236            ENDIF
237            ! new age of the FLOAT
238            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
239            ! test to know if the "age" of the float is not bigger than the
240            ! time step
241            IF( zagenewfl(jfl) > rdt ) THEN
242               zttfl(jfl) = (rdt-zagefl(jfl)) / zvol
243               zagenewfl(jfl) = rdt
244            ENDIF
245           
246            ! In the "minimal" direction we compute the index of new mesh
247            ! on i-direction
248            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
249               zgifl(jfl) = float(iioutfl(jfl))
250               ind = iioutfl(jfl)
251               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
252                  iioutfl(jfl) = iioutfl(jfl) + 1
253               ELSE
254                  iioutfl(jfl) = iioutfl(jfl) - 1
255               ENDIF
256               iiinfl(jfl) = ind
257            ELSE
258               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
259                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
260                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
261               ELSE
262                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
263               ENDIF
264            ENDIF
265            ! on j-direction
266            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
267               zgjfl(jfl) = float(ijoutfl(jfl))
268               ind = ijoutfl(jfl)
269               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
270                  ijoutfl(jfl) = ijoutfl(jfl) + 1
271               ELSE
272                  ijoutfl(jfl) = ijoutfl(jfl) - 1
273               ENDIF
274               ijinfl(jfl) = ind
275            ELSE
276               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
277                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
278                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
279               ELSE
280                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
281               ENDIF
282            ENDIF
283            ! on k-direction
284            IF( nisobfl(jfl) == 1. ) THEN
285               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
286                  zgkfl(jfl) = float(ikoutfl(jfl))
287                  ind = ikoutfl(jfl)
288                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
289                     ikoutfl(jfl) = ikoutfl(jfl)+1
290                  ELSE
291                     ikoutfl(jfl) = ikoutfl(jfl)-1
292                  ENDIF
293                  ikinfl(jfl) = ind
294               ELSE
295                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
296                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
297                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
298                  ELSE
299                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
300                  ENDIF
301               ENDIF
302            ENDIF
303           
304            ! coordinate of the new point on the temperature grid
305           
306            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
307            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
308            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
309!!Alexcadm   write(*,*)'PE ',narea,
310!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
311!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
312!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
313!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
314!!Alexcadm     .  zgjfl(jfl)
315!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
316!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
317!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
318!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
319!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
320!!Alexcadm     .  zgjfl(jfl)
321            ! reinitialisation of the age of FLOAT
322            zagefl(jfl) = zagenewfl(jfl)
323# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
324         ELSE
325            ! we put zgifl, zgjfl, zgkfl, zagefl
326            zgifl (jfl) = 0.
327            zgjfl (jfl) = 0.
328            zgkfl (jfl) = 0.
329            zagefl(jfl) = 0.
330            iil(jfl) = 0
331            ijl(jfl) = 0
332         ENDIF
333# endif
334      END DO
335     
336      ! synchronisation
337      IF( lk_mpp )   CALL mpp_sum( zgifl , jpnfl )   ! sums over the global domain
338      IF( lk_mpp )   CALL mpp_sum( zgjfl , jpnfl )
339      IF( lk_mpp )   CALL mpp_sum( zgkfl , jpnfl )
340      IF( lk_mpp )   CALL mpp_sum( zagefl, jpnfl )
341      IF( lk_mpp )   CALL mpp_sum( iil   , jpnfl )
342      IF( lk_mpp )   CALL mpp_sum( ijl   , jpnfl )
343     
344      ! in the case of open boundaries we need to test if the floats don't
345      ! go out of the domain. If it goes out, the float is put at the
346      ! middle of the mesh in the domain but the trajectory isn't compute
347      ! more time.     
348# if defined key_obc
349      DO jfl = 1, jpnfl
350         IF( lp_obc_east ) THEN
351            IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN
352               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
353               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
354               zagefl(jfl) = rdt
355            END IF
356         END IF
357         IF( lp_obc_west ) THEN
358            IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN
359               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
360               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
361               zagefl(jfl) = rdt
362            END IF
363         END IF
364         IF( lp_obc_north ) THEN
365            IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN
366               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
367               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
368               zagefl(jfl) = rdt
369            END IF
370         END IF
371         IF( lp_obc_south ) THEN
372            IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN
373               zgifl (jfl) = INT(zgifl(jfl)) + 0.5
374               zgjfl (jfl) = INT(zgjfl(jfl)) + 0.5
375               zagefl(jfl) = rdt
376            END IF
377         END IF
378      END DO
379#endif
380
381      ! Test to know if a  float hasn't integrated enought time
382      IF( ln_argo ) THEN
383         ifin = 1
384         DO jfl = 1, jpnfl
385            IF( zagefl(jfl) < rdt )   ifin = 0
386            tpifl(jfl) = zgifl(jfl) + 0.5
387            tpjfl(jfl) = zgjfl(jfl) + 0.5
388         END DO
389      ELSE
390         ifin = 1
391         DO jfl = 1, jpnfl
392            IF( zagefl(jfl) < rdt )   ifin = 0
393            tpifl(jfl) = zgifl(jfl) + 0.5
394            tpjfl(jfl) = zgjfl(jfl) + 0.5
395            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
396         END DO
397      ENDIF
398!!Alexcadm  IF (lwp) write(numout,*) '---------'
399!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
400!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
401!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
402!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
403!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
404!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
405      IF( ifin == 0 ) THEN
406         iloop = iloop + 1 
407         GO TO 222
408      ENDIF
409
410   END SUBROUTINE flo_blk
411
412#  else
413   !!----------------------------------------------------------------------
414   !!   Default option                                         Empty module
415   !!----------------------------------------------------------------------
416CONTAINS
417   SUBROUTINE flo_blk                  ! Empty routine
418   END SUBROUTINE flo_blk 
419#endif
420   
421   !!======================================================================
422END MODULE floblk 
Note: See TracBrowser for help on using the repository browser.