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

Last change on this file since 13216 was 13216, checked in by rblod, 4 years ago

Merge dev_r12973_AGRIF_CMEMS
  • Property svn:keywords set to Id
File size: 16.6 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   !!----------------------------------------------------------------------
23   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
24   !! $Id$
25   !! Software governed by the CeCILL license (see ./LICENSE)
26   !!----------------------------------------------------------------------
27CONTAINS
28
29   SUBROUTINE flo_blk( kt, Kbb, Kmm )
30      !!---------------------------------------------------------------------
31      !!                  ***  ROUTINE flo_blk  ***
32      !!           
33      !! ** Purpose :   Compute the geographical position,latitude, longitude
34      !!      and depth of each float at each time step.
35      !!
36      !! ** Method  :   The position of a float is computed with Bruno Blanke
37      !!      algorithm. We need to know the velocity field, the old positions
38      !!      of the floats and the grid defined on the domain.
39      !!----------------------------------------------------------------------
40      INTEGER, INTENT( in  ) ::   kt       ! ocean time step
41      INTEGER, INTENT( in  ) ::   Kbb, Kmm ! ocean time level indices
42      !!
43#ifndef key_agrif
44
45!RB super quick fix to compile with agrif
46
47      INTEGER :: jfl              ! dummy loop arguments
48      INTEGER :: ind, ifin, iloop
49      REAL(wp)   ::       &
50         zuinfl,zvinfl,zwinfl,      &     ! transport across the input face
51         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
52         zvol,                      &     ! volume of the mesh
53         zsurfz,                    &     ! surface of the face of the mesh
54         zind
55
56      REAL(wp), DIMENSION ( 2 )  ::   zsurfx, zsurfy   ! surface of the face of the mesh
57
58      INTEGER  , DIMENSION ( jpnfl )  ::   iil, ijl, ikl                   ! index of nearest mesh
59      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc             
60      INTEGER  , DIMENSION ( jpnfl )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float.
61      INTEGER  , DIMENSION ( jpnfl )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float.
62      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on
63      !                                                                         ! velocity mesh.
64      REAL(wp) , DIMENSION ( jpnfl )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh
65      !                                                                         ! across one of the face x,y and z
66      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh
67      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of
68      !                                                                         ! the float has been computed
69      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation
70      !                                                                         ! of new position
71      REAL(wp) , DIMENSION ( jpnfl )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position
72      REAL(wp) , DIMENSION ( jpnfl )  ::   zudfl, zvdfl, zwdfl             ! velocity diff input/output of mesh
73      REAL(wp) , DIMENSION ( jpnfl )  ::   zgidfl, zgjdfl, zgkdfl          ! direction index of float
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) .AND. iil(jfl) <= mig(nlei) .AND.   &
107             ijl(jfl) >= mjg(nldj) .AND. ijl(jfl) <= mjg(nlej)   ) THEN
108            iiloc(jfl) = iil(jfl) - mig(1) + 1
109            ijloc(jfl) = ijl(jfl) - mjg(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 e3. but never checked
117            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
118            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
119            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
120            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
121
122            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
123            zsurfz =          e1e2t(iiloc(jfl),ijloc(jfl))
124            zvol   = zsurfz * e3t(iiloc(jfl),ijloc(jfl),-ikl(jfl),Kmm)
125
126            !
127            zuinfl =( uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(1)
128            zuoutfl=( uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(2)
129            zvinfl =( vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kmm) )/2.*zsurfy(1)
130            zvoutfl=( vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kmm) )/2.*zsurfy(2)
131            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
132               &   +  ww(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
133            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
134               &   +  ww(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) = HUGE(1._wp)
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) = HUGE(1._wp)
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) = HUGE(1._wp)
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) > rn_Dt ) THEN
240               zttfl(jfl) = (rn_Dt-zagefl(jfl)) / zvol
241               zagenewfl(jfl) = rn_Dt
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      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
336      CALL mpp_sum( 'floblk', zgjfl , jpnfl )
337      CALL mpp_sum( 'floblk', zgkfl , jpnfl )
338      CALL mpp_sum( 'floblk', zagefl, jpnfl )
339      CALL mpp_sum( 'floblk', iil   , jpnfl )
340      CALL mpp_sum( 'floblk', ijl   , jpnfl )
341     
342      ! Test to know if a  float hasn't integrated enought time
343      IF( ln_argo ) THEN
344         ifin = 1
345         DO jfl = 1, jpnfl
346            IF( zagefl(jfl) < rn_Dt )   ifin = 0
347            tpifl(jfl) = zgifl(jfl) + 0.5
348            tpjfl(jfl) = zgjfl(jfl) + 0.5
349         END DO
350      ELSE
351         ifin = 1
352         DO jfl = 1, jpnfl
353            IF( zagefl(jfl) < rn_Dt )   ifin = 0
354            tpifl(jfl) = zgifl(jfl) + 0.5
355            tpjfl(jfl) = zgjfl(jfl) + 0.5
356            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
357         END DO
358      ENDIF
359!!Alexcadm  IF (lwp) write(numout,*) '---------'
360!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
361!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
362!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
363!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
364!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
365!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
366      IF( ifin == 0 ) THEN
367         iloop = iloop + 1 
368         GO TO 222
369      ENDIF
370#endif
371      !
372      !
373   END SUBROUTINE flo_blk
374
375   !!======================================================================
376END MODULE floblk 
Note: See TracBrowser for help on using the repository browser.