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/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/FLO – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/FLO/floblk.F90 @ 12724

Last change on this file since 12724 was 12724, checked in by techene, 4 years ago

branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
  • Property svn:keywords set to Id
File size: 16.5 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      INTEGER :: jfl              ! dummy loop arguments
46      INTEGER :: ind, ifin, iloop
47      REAL(wp)   ::       &
48         zuinfl,zvinfl,zwinfl,      &     ! transport across the input face
49         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
50         zvol,                      &     ! volume of the mesh
51         zsurfz,                    &     ! surface of the face of the mesh
52         zind
53
54      REAL(wp), DIMENSION ( 2 )  ::   zsurfx, zsurfy   ! surface of the face of the mesh
55
56      INTEGER  , DIMENSION ( jpnfl )  ::   iil, ijl, ikl                   ! index of nearest mesh
57      INTEGER  , DIMENSION ( jpnfl )  ::   iiloc , ijloc
58      INTEGER  , DIMENSION ( jpnfl )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float.
59      INTEGER  , DIMENSION ( jpnfl )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float.
60      REAL(wp) , DIMENSION ( jpnfl )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on
61      !                                                                         ! velocity mesh.
62      REAL(wp) , DIMENSION ( jpnfl )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh
63      !                                                                         ! across one of the face x,y and z
64      REAL(wp) , DIMENSION ( jpnfl )  ::    zttfl                          ! time for a float to quit the mesh
65      REAL(wp) , DIMENSION ( jpnfl )  ::    zagefl                         ! time during which, trajectorie of
66      !                                                                         ! the float has been computed
67      REAL(wp) , DIMENSION ( jpnfl )  ::   zagenewfl                       ! new age of float after calculation
68      !                                                                         ! of new position
69      REAL(wp) , DIMENSION ( jpnfl )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position
70      REAL(wp) , DIMENSION ( jpnfl )  ::   zudfl, zvdfl, zwdfl             ! velocity diff input/output of mesh
71      REAL(wp) , DIMENSION ( jpnfl )  ::   zgidfl, zgjdfl, zgkdfl          ! direction index of float
72      !!---------------------------------------------------------------------
73
74      IF( kt == nit000 ) THEN
75         IF(lwp) WRITE(numout,*)
76         IF(lwp) WRITE(numout,*) 'flo_blk : compute Blanke trajectories for floats '
77         IF(lwp) WRITE(numout,*) '~~~~~~~ '
78      ENDIF
79
80      ! Initialisation of parameters
81
82      DO jfl = 1, jpnfl
83         ! ages of floats are put at zero
84         zagefl(jfl) = 0.
85         ! index on the velocity grid
86         ! We considere k coordinate negative, with this transformation
87         ! the computation in the 3 direction is the same.
88         zgifl(jfl) = tpifl(jfl) - 0.5
89         zgjfl(jfl) = tpjfl(jfl) - 0.5
90         zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl)))
91         ! surface drift every 10 days
92         IF( ln_argo ) THEN
93            IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 )  zgkfl(jfl) = -1.
94         ENDIF
95         ! index of T mesh
96         iil(jfl) = 1 + INT(zgifl(jfl))
97         ijl(jfl) = 1 + INT(zgjfl(jfl))
98         ikl(jfl) =     INT(zgkfl(jfl))
99      END DO
100
101      iloop = 0
102222   DO jfl = 1, jpnfl
103# if   defined key_mpp_mpi
104         IF( iil(jfl) >= mig(nldi) .AND. iil(jfl) <= mig(nlei) .AND.   &
105             ijl(jfl) >= mjg(nldj) .AND. ijl(jfl) <= mjg(nlej)   ) THEN
106            iiloc(jfl) = iil(jfl) - mig(1) + 1
107            ijloc(jfl) = ijl(jfl) - mjg(1) + 1
108# else
109            iiloc(jfl) = iil(jfl)
110            ijloc(jfl) = ijl(jfl)
111# endif
112
113            ! compute the transport across the mesh where the float is.
114!!bug (gm) change e3t into e3. but never checked
115            zsurfx(1) =   &
116            &   e2u(iiloc(jfl)-1,ijloc(jfl)  )    &
117            & * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
118            zsurfx(2) =   &
119            &   e2u(iiloc(jfl)  ,ijloc(jfl)  )    &
120            & * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
121            zsurfy(1) =   &
122            &   e1v(iiloc(jfl)  ,ijloc(jfl)-1)    &
123            & * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
124            zsurfy(2) =   &
125            &   e1v(iiloc(jfl)  ,ijloc(jfl)  )    &
126            & * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
127
128            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
129            zsurfz =          e1e2t(iiloc(jfl),ijloc(jfl))
130            zvol   = zsurfz * e3t(iiloc(jfl),ijloc(jfl),-ikl(jfl),Kmm)
131
132            !
133            zuinfl =( uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(1)
134            zuoutfl=( uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(2)
135            zvinfl =( vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kmm) )/2.*zsurfy(1)
136            zvoutfl=( vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kmm) )/2.*zsurfy(2)
137            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
138               &   +  ww(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
139            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
140               &   +  ww(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
141
142            ! interpolation of velocity field on the float initial position
143            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
144            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
145            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
146
147            ! faces of input and output
148            ! u-direction
149            IF( zufl(jfl) < 0. ) THEN
150               iioutfl(jfl) = iil(jfl) - 1.
151               iiinfl (jfl) = iil(jfl)
152               zind   = zuinfl
153               zuinfl = zuoutfl
154               zuoutfl= zind
155            ELSE
156               iioutfl(jfl) = iil(jfl)
157               iiinfl (jfl) = iil(jfl) - 1
158            ENDIF
159            ! v-direction
160            IF( zvfl(jfl) < 0. ) THEN
161               ijoutfl(jfl) = ijl(jfl) - 1.
162               ijinfl (jfl) = ijl(jfl)
163               zind    = zvinfl
164               zvinfl  = zvoutfl
165               zvoutfl = zind
166            ELSE
167               ijoutfl(jfl) = ijl(jfl)
168               ijinfl (jfl) = ijl(jfl) - 1.
169            ENDIF
170            ! w-direction
171            IF( zwfl(jfl) < 0. ) THEN
172               ikoutfl(jfl) = ikl(jfl) - 1.
173               ikinfl (jfl) = ikl(jfl)
174               zind    = zwinfl
175               zwinfl  = zwoutfl
176               zwoutfl = zind
177            ELSE
178               ikoutfl(jfl) = ikl(jfl)
179               ikinfl (jfl) = ikl(jfl) - 1.
180            ENDIF
181
182            ! compute the time to go out the mesh across a face
183            ! u-direction
184            zudfl (jfl) = zuoutfl - zuinfl
185            zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl))
186            IF( zufl(jfl)*zuoutfl <= 0. ) THEN
187               ztxfl(jfl) = HUGE(1._wp)
188            ELSE
189               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
190                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
191               ELSE
192                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
193               ENDIF
194               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
195                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
196                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
197               ENDIF
198            ENDIF
199            ! v-direction
200            zvdfl (jfl) = zvoutfl - zvinfl
201            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
202            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
203               ztyfl(jfl) = HUGE(1._wp)
204            ELSE
205               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
206                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
207               ELSE
208                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
209               ENDIF
210               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
211                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
212                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
213               ENDIF
214            ENDIF
215            ! w-direction
216            IF( nisobfl(jfl) == 1. ) THEN
217               zwdfl (jfl) = zwoutfl - zwinfl
218               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
219               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
220                  ztzfl(jfl) = HUGE(1._wp)
221               ELSE
222                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
223                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
224                  ELSE
225                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
226                  ENDIF
227                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
228                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
229                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
230                  ENDIF
231               ENDIF
232            ENDIF
233
234            ! the time to go leave the mesh is the smallest time
235
236            IF( nisobfl(jfl) == 1. ) THEN
237               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
238            ELSE
239               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
240            ENDIF
241            ! new age of the FLOAT
242            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
243            ! test to know if the "age" of the float is not bigger than the
244            ! time step
245            IF( zagenewfl(jfl) > rn_Dt ) THEN
246               zttfl(jfl) = (rn_Dt-zagefl(jfl)) / zvol
247               zagenewfl(jfl) = rn_Dt
248            ENDIF
249
250            ! In the "minimal" direction we compute the index of new mesh
251            ! on i-direction
252            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
253               zgifl(jfl) = float(iioutfl(jfl))
254               ind = iioutfl(jfl)
255               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
256                  iioutfl(jfl) = iioutfl(jfl) + 1
257               ELSE
258                  iioutfl(jfl) = iioutfl(jfl) - 1
259               ENDIF
260               iiinfl(jfl) = ind
261            ELSE
262               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
263                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
264                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
265               ELSE
266                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
267               ENDIF
268            ENDIF
269            ! on j-direction
270            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
271               zgjfl(jfl) = float(ijoutfl(jfl))
272               ind = ijoutfl(jfl)
273               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
274                  ijoutfl(jfl) = ijoutfl(jfl) + 1
275               ELSE
276                  ijoutfl(jfl) = ijoutfl(jfl) - 1
277               ENDIF
278               ijinfl(jfl) = ind
279            ELSE
280               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
281                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
282                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
283               ELSE
284                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
285               ENDIF
286            ENDIF
287            ! on k-direction
288            IF( nisobfl(jfl) == 1. ) THEN
289               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
290                  zgkfl(jfl) = float(ikoutfl(jfl))
291                  ind = ikoutfl(jfl)
292                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
293                     ikoutfl(jfl) = ikoutfl(jfl)+1
294                  ELSE
295                     ikoutfl(jfl) = ikoutfl(jfl)-1
296                  ENDIF
297                  ikinfl(jfl) = ind
298               ELSE
299                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
300                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
301                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
302                  ELSE
303                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
304                  ENDIF
305               ENDIF
306            ENDIF
307
308            ! coordinate of the new point on the temperature grid
309
310            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
311            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
312            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
313!!Alexcadm   write(*,*)'PE ',narea,
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!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
320!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
321!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
322!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
323!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
324!!Alexcadm     .  zgjfl(jfl)
325            ! reinitialisation of the age of FLOAT
326            zagefl(jfl) = zagenewfl(jfl)
327# if   defined key_mpp_mpi
328         ELSE
329            ! we put zgifl, zgjfl, zgkfl, zagefl
330            zgifl (jfl) = 0.
331            zgjfl (jfl) = 0.
332            zgkfl (jfl) = 0.
333            zagefl(jfl) = 0.
334            iil(jfl) = 0
335            ijl(jfl) = 0
336         ENDIF
337# endif
338      END DO
339
340      ! synchronisation
341      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
342      CALL mpp_sum( 'floblk', zgjfl , jpnfl )
343      CALL mpp_sum( 'floblk', zgkfl , jpnfl )
344      CALL mpp_sum( 'floblk', zagefl, jpnfl )
345      CALL mpp_sum( 'floblk', iil   , jpnfl )
346      CALL mpp_sum( 'floblk', ijl   , jpnfl )
347
348      ! Test to know if a  float hasn't integrated enought time
349      IF( ln_argo ) THEN
350         ifin = 1
351         DO jfl = 1, jpnfl
352            IF( zagefl(jfl) < rn_Dt )   ifin = 0
353            tpifl(jfl) = zgifl(jfl) + 0.5
354            tpjfl(jfl) = zgjfl(jfl) + 0.5
355         END DO
356      ELSE
357         ifin = 1
358         DO jfl = 1, jpnfl
359            IF( zagefl(jfl) < rn_Dt )   ifin = 0
360            tpifl(jfl) = zgifl(jfl) + 0.5
361            tpjfl(jfl) = zgjfl(jfl) + 0.5
362            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
363         END DO
364      ENDIF
365!!Alexcadm  IF (lwp) write(numout,*) '---------'
366!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
367!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
368!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
369!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
370!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
371!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
372      IF( ifin == 0 ) THEN
373         iloop = iloop + 1
374         GO TO 222
375      ENDIF
376      !
377      !
378   END SUBROUTINE flo_blk
379
380   !!======================================================================
381END MODULE floblk
Note: See TracBrowser for help on using the repository browser.