source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/FLO/floblk.F90 @ 12603

Last change on this file since 12603 was 12603, checked in by orioltp, 6 months ago

Adding several interfaces to work with both single and double precision

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