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

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

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