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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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