source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/FLO/floblk.F90 @ 11574

Last change on this file since 11574 was 11573, checked in by jchanut, 20 months ago

#2222, merged with trunk

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