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 NEMO/trunk/src/OCE/FLO – NEMO

source: NEMO/trunk/src/OCE/FLO/floblk.F90 @ 13286

Last change on this file since 13286 was 13286, checked in by smasson, 4 years ago

trunk: merge extra halos branch in trunk, see #2366
  • Property svn:keywords set to Id
File size: 16.8 KB
Line 
1MODULE floblk
2   !!======================================================================
3   !!                     ***  MODULE  floblk  ***
4   !! Ocean floats :   trajectory computation
5   !!======================================================================
6   !!
7   !!----------------------------------------------------------------------
8   !!    flotblk     : compute float trajectories with Blanke algorithme
9   !!----------------------------------------------------------------------
10   USE flo_oce         ! ocean drifting floats
11   USE oce             ! ocean dynamics and tracers
12   USE dom_oce         ! ocean space and time domain
13   USE phycst          ! physical constants
14   USE in_out_manager  ! I/O manager
15   USE lib_mpp         ! distribued memory computing library
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC   flo_blk    ! routine called by floats.F90
21
22#  include "domzgr_substitute.h90"
23
24   !!----------------------------------------------------------------------
25   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
26   !! $Id$
27   !! Software governed by the CeCILL license (see ./LICENSE)
28   !!----------------------------------------------------------------------
29CONTAINS
30
31   SUBROUTINE flo_blk( kt, Kbb, Kmm )
32      !!---------------------------------------------------------------------
33      !!                  ***  ROUTINE flo_blk  ***
34      !!           
35      !! ** Purpose :   Compute the geographical position,latitude, longitude
36      !!      and depth of each float at each time step.
37      !!
38      !! ** Method  :   The position of a float is computed with Bruno Blanke
39      !!      algorithm. We need to know the velocity field, the old positions
40      !!      of the floats and the grid defined on the domain.
41      !!----------------------------------------------------------------------
42      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
43      INTEGER, INTENT( in  ) ::   Kbb, Kmm ! ocean time level indices
44      !!
45#ifndef key_agrif
46
47!RB super quick fix to compile with agrif
48
49      INTEGER :: jfl              ! dummy loop arguments
50      INTEGER :: ind, ifin, iloop
51      REAL(wp)   ::       &
52         zuinfl,zvinfl,zwinfl,      &     ! transport across the input face
53         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
54         zvol,                      &     ! volume of the mesh
55         zsurfz,                    &     ! surface of the face of the mesh
56         zind
57
58      REAL(wp), DIMENSION ( 2 )  ::   zsurfx, zsurfy   ! surface of the face of the mesh
59
60      INTEGER  , DIMENSION ( jpnfl )  ::   iil, ijl, ikl                   ! index of nearest mesh
61      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc             
62      INTEGER  , DIMENSION ( jpnfl )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float.
63      INTEGER  , DIMENSION ( jpnfl )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float.
64      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on
65      !                                                                         ! velocity mesh.
66      REAL(wp) , DIMENSION ( jpnfl )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh
67      !                                                                         ! across one of the face x,y and z
68      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh
69      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of
70      !                                                                         ! the float has been computed
71      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation
72      !                                                                         ! of new position
73      REAL(wp) , DIMENSION ( jpnfl )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position
74      REAL(wp) , DIMENSION ( jpnfl )  ::   zudfl, zvdfl, zwdfl             ! velocity diff input/output of mesh
75      REAL(wp) , DIMENSION ( jpnfl )  ::   zgidfl, zgjdfl, zgkdfl          ! direction index of float
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_mpi
108         IF( iil(jfl) >= mig(Nis0) .AND. iil(jfl) <= mig(Nie0) .AND.   &
109             ijl(jfl) >= mjg(Njs0) .AND. ijl(jfl) <= mjg(Nje0)   ) THEN
110            iiloc(jfl) = iil(jfl) - mig(1) + 1
111            ijloc(jfl) = ijl(jfl) - mjg(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!!bug (gm) change e3t into e3. but never checked
119            zsurfx(1) =   &
120            &   e2u(iiloc(jfl)-1,ijloc(jfl)  )    &
121            & * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
122            zsurfx(2) =   &
123            &   e2u(iiloc(jfl)  ,ijloc(jfl)  )    &
124            & * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
125            zsurfy(1) =   &
126            &   e1v(iiloc(jfl)  ,ijloc(jfl)-1)    &
127            & * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
128            zsurfy(2) =   &
129            &   e1v(iiloc(jfl)  ,ijloc(jfl)  )    &
130            & * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
131
132            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
133            zsurfz =          e1e2t(iiloc(jfl),ijloc(jfl))
134            zvol   = zsurfz * e3t(iiloc(jfl),ijloc(jfl),-ikl(jfl),Kmm)
135
136            !
137            zuinfl =( uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(1)
138            zuoutfl=( uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(2)
139            zvinfl =( vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kmm) )/2.*zsurfy(1)
140            zvoutfl=( vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kmm) )/2.*zsurfy(2)
141            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
142               &   +  ww(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
143            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
144               &   +  ww(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
145           
146            ! interpolation of velocity field on the float initial position           
147            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
148            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
149            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
150           
151            ! faces of input and output
152            ! u-direction
153            IF( zufl(jfl) < 0. ) THEN
154               iioutfl(jfl) = iil(jfl) - 1.
155               iiinfl (jfl) = iil(jfl)
156               zind   = zuinfl
157               zuinfl = zuoutfl
158               zuoutfl= zind
159            ELSE
160               iioutfl(jfl) = iil(jfl)
161               iiinfl (jfl) = iil(jfl) - 1
162            ENDIF
163            ! v-direction       
164            IF( zvfl(jfl) < 0. ) THEN
165               ijoutfl(jfl) = ijl(jfl) - 1.
166               ijinfl (jfl) = ijl(jfl)
167               zind    = zvinfl
168               zvinfl  = zvoutfl
169               zvoutfl = zind
170            ELSE
171               ijoutfl(jfl) = ijl(jfl)
172               ijinfl (jfl) = ijl(jfl) - 1.
173            ENDIF
174            ! w-direction
175            IF( zwfl(jfl) < 0. ) THEN
176               ikoutfl(jfl) = ikl(jfl) - 1.
177               ikinfl (jfl) = ikl(jfl)
178               zind    = zwinfl
179               zwinfl  = zwoutfl
180               zwoutfl = zind
181            ELSE
182               ikoutfl(jfl) = ikl(jfl)
183               ikinfl (jfl) = ikl(jfl) - 1.
184            ENDIF
185           
186            ! compute the time to go out the mesh across a face
187            ! u-direction
188            zudfl (jfl) = zuoutfl - zuinfl
189            zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl))
190            IF( zufl(jfl)*zuoutfl <= 0. ) THEN
191               ztxfl(jfl) = HUGE(1._wp)
192            ELSE
193               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
194                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
195               ELSE
196                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
197               ENDIF
198               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
199                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
200                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
201               ENDIF
202            ENDIF
203            ! v-direction
204            zvdfl (jfl) = zvoutfl - zvinfl
205            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
206            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
207               ztyfl(jfl) = HUGE(1._wp)
208            ELSE
209               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
210                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
211               ELSE
212                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
213               ENDIF
214               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
215                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
216                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
217               ENDIF
218            ENDIF
219            ! w-direction       
220            IF( nisobfl(jfl) == 1. ) THEN
221               zwdfl (jfl) = zwoutfl - zwinfl
222               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
223               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
224                  ztzfl(jfl) = HUGE(1._wp)
225               ELSE
226                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
227                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
228                  ELSE
229                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
230                  ENDIF
231                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
232                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
233                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
234                  ENDIF
235               ENDIF
236            ENDIF
237           
238            ! the time to go leave the mesh is the smallest time
239                   
240            IF( nisobfl(jfl) == 1. ) THEN
241               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
242            ELSE
243               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
244            ENDIF
245            ! new age of the FLOAT
246            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
247            ! test to know if the "age" of the float is not bigger than the
248            ! time step
249            IF( zagenewfl(jfl) > rn_Dt ) THEN
250               zttfl(jfl) = (rn_Dt-zagefl(jfl)) / zvol
251               zagenewfl(jfl) = rn_Dt
252            ENDIF
253           
254            ! In the "minimal" direction we compute the index of new mesh
255            ! on i-direction
256            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
257               zgifl(jfl) = float(iioutfl(jfl))
258               ind = iioutfl(jfl)
259               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
260                  iioutfl(jfl) = iioutfl(jfl) + 1
261               ELSE
262                  iioutfl(jfl) = iioutfl(jfl) - 1
263               ENDIF
264               iiinfl(jfl) = ind
265            ELSE
266               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
267                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
268                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
269               ELSE
270                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
271               ENDIF
272            ENDIF
273            ! on j-direction
274            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
275               zgjfl(jfl) = float(ijoutfl(jfl))
276               ind = ijoutfl(jfl)
277               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
278                  ijoutfl(jfl) = ijoutfl(jfl) + 1
279               ELSE
280                  ijoutfl(jfl) = ijoutfl(jfl) - 1
281               ENDIF
282               ijinfl(jfl) = ind
283            ELSE
284               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
285                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
286                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
287               ELSE
288                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
289               ENDIF
290            ENDIF
291            ! on k-direction
292            IF( nisobfl(jfl) == 1. ) THEN
293               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
294                  zgkfl(jfl) = float(ikoutfl(jfl))
295                  ind = ikoutfl(jfl)
296                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
297                     ikoutfl(jfl) = ikoutfl(jfl)+1
298                  ELSE
299                     ikoutfl(jfl) = ikoutfl(jfl)-1
300                  ENDIF
301                  ikinfl(jfl) = ind
302               ELSE
303                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
304                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
305                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
306                  ELSE
307                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
308                  ENDIF
309               ENDIF
310            ENDIF
311           
312            ! coordinate of the new point on the temperature grid
313           
314            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
315            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
316            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
317!!Alexcadm   write(*,*)'PE ',narea,
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!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
324!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
325!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
326!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
327!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
328!!Alexcadm     .  zgjfl(jfl)
329            ! reinitialisation of the age of FLOAT
330            zagefl(jfl) = zagenewfl(jfl)
331# if   defined key_mpp_mpi
332         ELSE
333            ! we put zgifl, zgjfl, zgkfl, zagefl
334            zgifl (jfl) = 0.
335            zgjfl (jfl) = 0.
336            zgkfl (jfl) = 0.
337            zagefl(jfl) = 0.
338            iil(jfl) = 0
339            ijl(jfl) = 0
340         ENDIF
341# endif
342      END DO
343     
344      ! synchronisation
345      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
346      CALL mpp_sum( 'floblk', zgjfl , jpnfl )
347      CALL mpp_sum( 'floblk', zgkfl , jpnfl )
348      CALL mpp_sum( 'floblk', zagefl, jpnfl )
349      CALL mpp_sum( 'floblk', iil   , jpnfl )
350      CALL mpp_sum( 'floblk', ijl   , jpnfl )
351     
352      ! Test to know if a  float hasn't integrated enought time
353      IF( ln_argo ) THEN
354         ifin = 1
355         DO jfl = 1, jpnfl
356            IF( zagefl(jfl) < rn_Dt )   ifin = 0
357            tpifl(jfl) = zgifl(jfl) + 0.5
358            tpjfl(jfl) = zgjfl(jfl) + 0.5
359         END DO
360      ELSE
361         ifin = 1
362         DO jfl = 1, jpnfl
363            IF( zagefl(jfl) < rn_Dt )   ifin = 0
364            tpifl(jfl) = zgifl(jfl) + 0.5
365            tpjfl(jfl) = zgjfl(jfl) + 0.5
366            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
367         END DO
368      ENDIF
369!!Alexcadm  IF (lwp) write(numout,*) '---------'
370!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
371!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
372!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
373!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
374!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
375!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
376      IF( ifin == 0 ) THEN
377         iloop = iloop + 1 
378         GO TO 222
379      ENDIF
380#endif
381      !
382      !
383   END SUBROUTINE flo_blk
384
385   !!======================================================================
386END MODULE floblk 
Note: See TracBrowser for help on using the repository browser.