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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90 @ 2335

Last change on this file since 2335 was 2335, checked in by gm, 13 years ago

v3.3beta: Suppress obsolete key_mpp_shmem

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