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

source: trunk/NEMOGCM/NEMO/OPA_SRC/FLO/floblk.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

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