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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/FLO – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/FLO/floblk.F90 @ 10970

Last change on this file since 10970 was 10970, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : CRS and FLO. Only tested compilation. Note that base code doesn't compile with key_floats (#2279), so changes to FLO not really tested at all.
  • Property svn:keywords set to Id
File size: 17.0 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 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) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * e3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl),Kmm)
116            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * e3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
117            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * e3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl),Kmm)
118            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * e3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl),Kmm)
119
120            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
121            zsurfz =          e1e2t(iiloc(jfl),ijloc(jfl))
122            zvol   = zsurfz * e3t(iiloc(jfl),ijloc(jfl),-ikl(jfl),Kmm)
123
124            !
125            zuinfl =( uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(1)
126            zuoutfl=( uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kbb) + uu(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl),Kmm) )/2.*zsurfx(2)
127            zvinfl =( vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl),Kmm) )/2.*zsurfy(1)
128            zvoutfl=( vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kbb) + vv(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl),Kmm) )/2.*zsurfy(2)
129            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
130               &   +  ww(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
131            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
132               &   +  ww(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
133           
134            ! interpolation of velocity field on the float initial position           
135            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
136            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
137            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
138           
139            ! faces of input and output
140            ! u-direction
141            IF( zufl(jfl) < 0. ) THEN
142               iioutfl(jfl) = iil(jfl) - 1.
143               iiinfl (jfl) = iil(jfl)
144               zind   = zuinfl
145               zuinfl = zuoutfl
146               zuoutfl= zind
147            ELSE
148               iioutfl(jfl) = iil(jfl)
149               iiinfl (jfl) = iil(jfl) - 1
150            ENDIF
151            ! v-direction       
152            IF( zvfl(jfl) < 0. ) THEN
153               ijoutfl(jfl) = ijl(jfl) - 1.
154               ijinfl (jfl) = ijl(jfl)
155               zind    = zvinfl
156               zvinfl  = zvoutfl
157               zvoutfl = zind
158            ELSE
159               ijoutfl(jfl) = ijl(jfl)
160               ijinfl (jfl) = ijl(jfl) - 1.
161            ENDIF
162            ! w-direction
163            IF( zwfl(jfl) < 0. ) THEN
164               ikoutfl(jfl) = ikl(jfl) - 1.
165               ikinfl (jfl) = ikl(jfl)
166               zind    = zwinfl
167               zwinfl  = zwoutfl
168               zwoutfl = zind
169            ELSE
170               ikoutfl(jfl) = ikl(jfl)
171               ikinfl (jfl) = ikl(jfl) - 1.
172            ENDIF
173           
174            ! compute the time to go out the mesh across a face
175            ! u-direction
176            zudfl (jfl) = zuoutfl - zuinfl
177            zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl))
178            IF( zufl(jfl)*zuoutfl <= 0. ) THEN
179               ztxfl(jfl) = 1.E99
180            ELSE
181               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
182                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
183               ELSE
184                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
185               ENDIF
186               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
187                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
188                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
189               ENDIF
190            ENDIF
191            ! v-direction
192            zvdfl (jfl) = zvoutfl - zvinfl
193            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
194            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
195               ztyfl(jfl) = 1.E99
196            ELSE
197               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
198                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
199               ELSE
200                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
201               ENDIF
202               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
203                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
204                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
205               ENDIF
206            ENDIF
207            ! w-direction       
208            IF( nisobfl(jfl) == 1. ) THEN
209               zwdfl (jfl) = zwoutfl - zwinfl
210               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
211               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
212                  ztzfl(jfl) = 1.E99
213               ELSE
214                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
215                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
216                  ELSE
217                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
218                  ENDIF
219                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
220                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
221                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
222                  ENDIF
223               ENDIF
224            ENDIF
225           
226            ! the time to go leave the mesh is the smallest time
227                   
228            IF( nisobfl(jfl) == 1. ) THEN
229               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
230            ELSE
231               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
232            ENDIF
233            ! new age of the FLOAT
234            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
235            ! test to know if the "age" of the float is not bigger than the
236            ! time step
237            IF( zagenewfl(jfl) > rdt ) THEN
238               zttfl(jfl) = (rdt-zagefl(jfl)) / zvol
239               zagenewfl(jfl) = rdt
240            ENDIF
241           
242            ! In the "minimal" direction we compute the index of new mesh
243            ! on i-direction
244            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
245               zgifl(jfl) = float(iioutfl(jfl))
246               ind = iioutfl(jfl)
247               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
248                  iioutfl(jfl) = iioutfl(jfl) + 1
249               ELSE
250                  iioutfl(jfl) = iioutfl(jfl) - 1
251               ENDIF
252               iiinfl(jfl) = ind
253            ELSE
254               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
255                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
256                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
257               ELSE
258                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
259               ENDIF
260            ENDIF
261            ! on j-direction
262            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
263               zgjfl(jfl) = float(ijoutfl(jfl))
264               ind = ijoutfl(jfl)
265               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
266                  ijoutfl(jfl) = ijoutfl(jfl) + 1
267               ELSE
268                  ijoutfl(jfl) = ijoutfl(jfl) - 1
269               ENDIF
270               ijinfl(jfl) = ind
271            ELSE
272               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
273                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
274                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
275               ELSE
276                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
277               ENDIF
278            ENDIF
279            ! on k-direction
280            IF( nisobfl(jfl) == 1. ) THEN
281               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
282                  zgkfl(jfl) = float(ikoutfl(jfl))
283                  ind = ikoutfl(jfl)
284                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
285                     ikoutfl(jfl) = ikoutfl(jfl)+1
286                  ELSE
287                     ikoutfl(jfl) = ikoutfl(jfl)-1
288                  ENDIF
289                  ikinfl(jfl) = ind
290               ELSE
291                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
292                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
293                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
294                  ELSE
295                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
296                  ENDIF
297               ENDIF
298            ENDIF
299           
300            ! coordinate of the new point on the temperature grid
301           
302            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
303            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
304            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
305!!Alexcadm   write(*,*)'PE ',narea,
306!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
307!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
308!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
309!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
310!!Alexcadm     .  zgjfl(jfl)
311!!Alexcadm  IF (jfl == 910) write(*,*)'Flotteur 910',
312!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
313!!Alexcadm     .     ,ijoutfl(jfl),ikinfl(jfl),
314!!Alexcadm     .    ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
315!!Alexcadm     .     ,ztzfl(jfl),zgifl(jfl),
316!!Alexcadm     .  zgjfl(jfl)
317            ! reinitialisation of the age of FLOAT
318            zagefl(jfl) = zagenewfl(jfl)
319# if   defined key_mpp_mpi
320         ELSE
321            ! we put zgifl, zgjfl, zgkfl, zagefl
322            zgifl (jfl) = 0.
323            zgjfl (jfl) = 0.
324            zgkfl (jfl) = 0.
325            zagefl(jfl) = 0.
326            iil(jfl) = 0
327            ijl(jfl) = 0
328         ENDIF
329# endif
330      END DO
331     
332      ! synchronisation
333      CALL mpp_sum( 'floblk', zgifl , jpnfl )   ! sums over the global domain
334      CALL mpp_sum( 'floblk', zgjfl , jpnfl )
335      CALL mpp_sum( 'floblk', zgkfl , jpnfl )
336      CALL mpp_sum( 'floblk', zagefl, jpnfl )
337      CALL mpp_sum( 'floblk', iil   , jpnfl )
338      CALL mpp_sum( 'floblk', ijl   , jpnfl )
339     
340      ! Test to know if a  float hasn't integrated enought time
341      IF( ln_argo ) THEN
342         ifin = 1
343         DO jfl = 1, jpnfl
344            IF( zagefl(jfl) < rdt )   ifin = 0
345            tpifl(jfl) = zgifl(jfl) + 0.5
346            tpjfl(jfl) = zgjfl(jfl) + 0.5
347         END DO
348      ELSE
349         ifin = 1
350         DO jfl = 1, jpnfl
351            IF( zagefl(jfl) < rdt )   ifin = 0
352            tpifl(jfl) = zgifl(jfl) + 0.5
353            tpjfl(jfl) = zgjfl(jfl) + 0.5
354            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
355         END DO
356      ENDIF
357!!Alexcadm  IF (lwp) write(numout,*) '---------'
358!!Alexcadm  IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
359!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
360!!Alexcadm  IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
361!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
362!!Alexcadm  IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
363!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
364      IF( ifin == 0 ) THEN
365         iloop = iloop + 1 
366         GO TO 222
367      ENDIF
368      !
369      !
370   END SUBROUTINE flo_blk
371
372#  else
373   !!----------------------------------------------------------------------
374   !!   Default option                                         Empty module
375   !!----------------------------------------------------------------------
376CONTAINS
377   SUBROUTINE flo_blk                  ! Empty routine
378   END SUBROUTINE flo_blk 
379#endif
380   
381   !!======================================================================
382END MODULE floblk 
Note: See TracBrowser for help on using the repository browser.