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

source: NEMO/branches/UKMO/NEMO4_beta_mirror/src/OCE/FLO/floblk.F90 @ 9950

Last change on this file since 9950 was 9950, checked in by davestorkey, 6 years ago

UKMO/NEMO4_beta_mirror branch: remove SVN keywords

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