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_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/FLO – NEMO

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

Last change on this file since 13135 was 13135, checked in by orioltp, 4 years ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: merge with trunk@13134, see #2364
  • Property svn:keywords set to Id
File size: 16.7 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<<<<<<< .working
178               ztxfl(jfl) = HUGE(0.0_wp)
179||||||| .merge-left.r12512
180               ztxfl(jfl) = 1.E99
181=======
182               ztxfl(jfl) = HUGE(1._wp)
183>>>>>>> .merge-right.r13134
184            ELSE
185               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
186                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
187               ELSE
188                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
189               ENDIF
190               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
191                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
192                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
193               ENDIF
194            ENDIF
195            ! v-direction
196            zvdfl (jfl) = zvoutfl - zvinfl
197            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
198            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
199               ztyfl(jfl) = HUGE(1._wp)
200            ELSE
201               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
202                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
203               ELSE
204                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
205               ENDIF
206               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
207                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
208                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
209               ENDIF
210            ENDIF
211            ! w-direction       
212            IF( nisobfl(jfl) == 1. ) THEN
213               zwdfl (jfl) = zwoutfl - zwinfl
214               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
215               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
216                  ztzfl(jfl) = HUGE(1._wp)
217               ELSE
218                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
219                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
220                  ELSE
221                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
222                  ENDIF
223                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
224                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
225                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
226                  ENDIF
227               ENDIF
228            ENDIF
229           
230            ! the time to go leave the mesh is the smallest time
231                   
232            IF( nisobfl(jfl) == 1. ) THEN
233               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
234            ELSE
235               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
236            ENDIF
237            ! new age of the FLOAT
238            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
239            ! test to know if the "age" of the float is not bigger than the
240            ! time step
241            IF( zagenewfl(jfl) > rn_Dt ) THEN
242               zttfl(jfl) = (rn_Dt-zagefl(jfl)) / zvol
243               zagenewfl(jfl) = rn_Dt
244            ENDIF
245           
246            ! In the "minimal" direction we compute the index of new mesh
247            ! on i-direction
248            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
249               zgifl(jfl) = float(iioutfl(jfl))
250               ind = iioutfl(jfl)
251               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
252                  iioutfl(jfl) = iioutfl(jfl) + 1
253               ELSE
254                  iioutfl(jfl) = iioutfl(jfl) - 1
255               ENDIF
256               iiinfl(jfl) = ind
257            ELSE
258               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
259                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
260                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
261               ELSE
262                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
263               ENDIF
264            ENDIF
265            ! on j-direction
266            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
267               zgjfl(jfl) = float(ijoutfl(jfl))
268               ind = ijoutfl(jfl)
269               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
270                  ijoutfl(jfl) = ijoutfl(jfl) + 1
271               ELSE
272                  ijoutfl(jfl) = ijoutfl(jfl) - 1
273               ENDIF
274               ijinfl(jfl) = ind
275            ELSE
276               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
277                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
278                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
279               ELSE
280                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
281               ENDIF
282            ENDIF
283            ! on k-direction
284            IF( nisobfl(jfl) == 1. ) THEN
285               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
286                  zgkfl(jfl) = float(ikoutfl(jfl))
287                  ind = ikoutfl(jfl)
288                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
289                     ikoutfl(jfl) = ikoutfl(jfl)+1
290                  ELSE
291                     ikoutfl(jfl) = ikoutfl(jfl)-1
292                  ENDIF
293                  ikinfl(jfl) = ind
294               ELSE
295                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
296                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
297                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
298                  ELSE
299                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
300                  ENDIF
301               ENDIF
302            ENDIF
303           
304            ! coordinate of the new point on the temperature grid
305           
306            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
307            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
308            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
309!!Alexcadm   write(*,*)'PE ',narea,
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!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
316!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
317!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
318!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
319!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
320!!Alexcadm     .  zgjfl(jfl)
321            ! reinitialisation of the age of FLOAT
322            zagefl(jfl) = zagenewfl(jfl)
323# if   defined key_mpp_mpi
324         ELSE
325            ! we put zgifl, zgjfl, zgkfl, zagefl
326            zgifl (jfl) = 0.
327            zgjfl (jfl) = 0.
328            zgkfl (jfl) = 0.
329            zagefl(jfl) = 0.
330            iil(jfl) = 0
331            ijl(jfl) = 0
332         ENDIF
333# endif
334      END DO
335     
336      ! synchronisation
337      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
338      CALL mpp_sum( 'floblk', zgjfl , jpnfl )
339      CALL mpp_sum( 'floblk', zgkfl , jpnfl )
340      CALL mpp_sum( 'floblk', zagefl, jpnfl )
341      CALL mpp_sum( 'floblk', iil   , jpnfl )
342      CALL mpp_sum( 'floblk', ijl   , jpnfl )
343     
344      ! Test to know if a  float hasn't integrated enought time
345      IF( ln_argo ) THEN
346         ifin = 1
347         DO jfl = 1, jpnfl
348            IF( zagefl(jfl) < rn_Dt )   ifin = 0
349            tpifl(jfl) = zgifl(jfl) + 0.5
350            tpjfl(jfl) = zgjfl(jfl) + 0.5
351         END DO
352      ELSE
353         ifin = 1
354         DO jfl = 1, jpnfl
355            IF( zagefl(jfl) < rn_Dt )   ifin = 0
356            tpifl(jfl) = zgifl(jfl) + 0.5
357            tpjfl(jfl) = zgjfl(jfl) + 0.5
358            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
359         END DO
360      ENDIF
361!!Alexcadm  IF (lwp) write(numout,*) '---------'
362!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
363!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
364!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
365!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
366!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
367!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
368      IF( ifin == 0 ) THEN
369         iloop = iloop + 1 
370         GO TO 222
371      ENDIF
372      !
373      !
374   END SUBROUTINE flo_blk
375
376   !!======================================================================
377END MODULE floblk 
Note: See TracBrowser for help on using the repository browser.