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.
flodom.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 @ 2335

Last change on this file since 2335 was 2335, checked in by gm, 14 years ago

v3.3beta: Suppress obsolete key_mpp_shmem

  • Property svn:keywords set to Id
File size: 20.3 KB
Line 
1MODULE flodom
2   !!======================================================================
3   !!                       ***  MODULE  flodom  ***
4   !! Ocean floats :   domain
5   !!======================================================================
6   !! History :  OPA  ! 1998-07 (Y.Drillet, CLIPPER)  Original code
7   !!----------------------------------------------------------------------
8#if   defined key_floats   ||   defined key_esopa
9   !!----------------------------------------------------------------------
10   !!   'key_floats'                                     float trajectories
11   !!----------------------------------------------------------------------
12   !!   flo_dom        : initialization of floats
13   !!   findmesh       : compute index of position
14   !!   dstnce         : compute distance between face mesh and floats
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE flo_oce         ! ocean drifting floats
19   USE in_out_manager  ! I/O manager
20   USE lib_mpp         ! distribued memory computing library
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   flo_dom    ! routine called by floats.F90
26
27   !! * Substitutions
28#  include "domzgr_substitute.h90"
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE flo_dom
37      !! ---------------------------------------------------------------------
38      !!                  ***  ROUTINE flo_dom  ***
39      !!                 
40      !!  ** Purpose :   Initialisation of floats
41      !!
42      !!  ** Method  :   We put the floats  in the domain with the latitude,
43      !!               the longitude (degree) and the depth (m).
44      !!----------------------------------------------------------------------     
45      LOGICAL  ::   llinmesh
46      INTEGER  ::   ji, jj, jk   ! DO loop index on 3 directions
47      INTEGER  ::   jfl, jfl1    ! number of floats   
48      INTEGER  ::   inum         ! logical unit for file read
49      INTEGER, DIMENSION(jpnfl) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats
50      INTEGER, DIMENSION(jpnfl) ::   idomfl,  ivtest, ihtest   !   -             -
51      REAL(wp) ::   zdxab, zdyad
52      REAL(wp), DIMENSION(jpnnewflo+1)  :: zgifl, zgjfl,  zgkfl
53      !!---------------------------------------------------------------------
54     
55      ! Initialisation with the geographical position or restart
56     
57      IF(lwp) WRITE(numout,*) 'flo_dom : compute initial position of floats'
58      IF(lwp) WRITE(numout,*) '~~~~~~~~'
59      IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl
60     
61      IF(ln_rstflo) THEN
62         IF(lwp) WRITE(numout,*) '        float restart file read'
63         
64         ! open the restart file
65         CALL ctl_opn( inum, 'restart_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
66
67         ! read of the restart file
68         READ(inum) ( tpifl  (jfl), jfl=1, jpnrstflo),   & 
69                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   &
70                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   &
71                        ( nisobfl(jfl), jfl=1, jpnrstflo),   &
72                        ( ngrpfl (jfl), jfl=1, jpnrstflo)   
73         CLOSE(inum)
74
75         ! if we want a  surface drift  ( like PROVOR floats )
76         IF( ln_argo ) THEN
77            DO jfl = 1, jpnrstflo
78               nisobfl(jfl) = 0
79            END DO
80         ENDIF
81
82         IF(lwp) WRITE(numout,*)' flo_dom: END of florstlec'
83         
84         ! It is possible to add new floats.         
85         IF(lwp) WRITE(numout,*)' flo_dom:jpnfl jpnrstflo ',jpnfl,jpnrstflo
86         IF( jpnfl > jpnrstflo ) THEN
87            ! open the init file
88            CALL ctl_opn( inum, 'init_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
89            DO jfl = jpnrstflo+1, jpnfl
90               READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),jfl1
91            END DO
92            CLOSE(inum)
93            IF(lwp) WRITE(numout,*)' flodom: END reading init_float file'
94           
95            ! Test to find the grid point coordonate with the geographical position           
96            DO jfl = jpnrstflo+1, jpnfl
97               ihtest(jfl) = 0
98               ivtest(jfl) = 0
99               ikmfl(jfl) = 0
100# if   defined key_mpp_mpi
101               DO ji = MAX(nldi,2), nlei
102                  DO jj = MAX(nldj,2), nlej   ! NO vector opt.
103# else
104               DO ji = 2, jpi
105                  DO jj = 2, jpj   ! NO vector opt.
106# endif                     
107                     ! For each float we find the indexes of the mesh                     
108                     CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
109                                   glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
110                                   glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
111                                   glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
112                                   flxx(jfl)       ,flyy(jfl)       ,   &
113                                   glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
114                     IF(llinmesh) THEN
115                        iimfl(jfl) = ji
116                        ijmfl(jfl) = jj
117                        ihtest(jfl) = ihtest(jfl)+1
118                        DO jk = 1, jpk-1
119                           IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN
120                              ikmfl(jfl) = jk
121                              ivtest(jfl) = ivtest(jfl) + 1
122                           ENDIF
123                        END DO
124                     ENDIF
125                  END DO
126               END DO
127               IF(lwp) WRITE(numout,*)'   flo_dom: END findmesh'
128               
129               ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1               
130               IF( ihtest(jfl) ==  0 ) THEN
131                  iimfl(jfl) = -1
132                  ijmfl(jfl) = -1
133               ENDIF
134            END DO
135           
136            ! A zero in the sum of the arrays "ihtest" and "ivtest"             
137# if   defined key_mpp_mpi
138            CALL mpp_sum(ihtest,jpnfl)
139            CALL mpp_sum(ivtest,jpnfl)
140# endif
141            DO jfl = jpnrstflo+1, jpnfl
142               IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
143                  IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
144                  STOP
145               ENDIF
146               IF( (ihtest(jfl) == 0) ) THEN
147                  IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH'
148                  STOP
149               ENDIF
150            END DO
151           
152            ! We compute the distance between the float and the face of the mesh           
153            DO jfl = jpnrstflo+1, jpnfl               
154               ! Made only if the float is in the domain of the processor             
155               IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN
156                 
157                  ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
158                 
159                  idomfl(jfl) = 0
160                  IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1
161                                           
162                  ! Computation of the distance between the float and the faces of the mesh
163                  !            zdxab
164                  !             .
165                  !        B----.---------C
166                  !        |    .         |
167                  !        |<------>flo   |
168                  !        |        ^     |
169                  !        |        |.....|....zdyad
170                  !        |        |     |
171                  !        A--------|-----D
172                  !
173             
174                  zdxab = dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) )
175                  zdyad = dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) )
176                 
177                  ! Translation of this distances (in meter) in indexes
178                 
179                  zgifl(jfl-jpnrstflo)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
180                  zgjfl(jfl-jpnrstflo)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
181                  zgkfl(jfl-jpnrstflo) = (( fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))   &
182                     &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
183                     &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) )                             &
184                     &                 + (( flzz(jfl)-fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1))   &
185                     &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
186                     &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) )
187               ELSE
188                  zgifl(jfl-jpnrstflo) = 0.e0
189                  zgjfl(jfl-jpnrstflo) = 0.e0
190                  zgkfl(jfl-jpnrstflo) = 0.e0
191               ENDIF
192            END DO
193           
194            ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats.
195            IF( lk_mpp )   THEN
196               CALL mpp_sum( zgjfl, jpnnewflo )   ! sums over the global domain
197               CALL mpp_sum( zgkfl, jpnnewflo )
198               IF(lwp) WRITE(numout,*) (zgifl(jfl),jfl=1,jpnnewflo)
199               IF(lwp) WRITE(numout,*) (zgjfl(jfl),jfl=1,jpnnewflo)
200               IF(lwp) WRITE(numout,*) (zgkfl(jfl),jfl=1,jpnnewflo) 
201            ENDIF
202           
203            DO jfl = jpnrstflo+1, jpnfl
204               tpifl(jfl) = zgifl(jfl-jpnrstflo)
205               tpjfl(jfl) = zgjfl(jfl-jpnrstflo)
206               tpkfl(jfl) = zgkfl(jfl-jpnrstflo)
207            END DO
208         ENDIF
209      ELSE
210         IF(lwp) WRITE(numout,*) '                     init_float read '
211         
212         ! First initialisation of floats
213         ! the initials positions of floats are written in a file
214         ! with a variable to know if it is a isobar float a number
215         ! to identified who want the trajectories of this float and
216         ! an index for the number of the float         
217         ! open the init file
218         CALL ctl_opn( inum, 'init_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
219         READ(inum) (flxx(jfl)   , jfl=1, jpnfl),   &
220                    (flyy(jfl)   , jfl=1, jpnfl),   &
221                    (flzz(jfl)   , jfl=1, jpnfl),   &
222                    (nisobfl(jfl), jfl=1, jpnfl),   &
223                    (ngrpfl(jfl) , jfl=1, jpnfl)
224         CLOSE(inum)
225           
226         ! Test to find the grid point coordonate with the geographical position         
227         DO jfl = 1, jpnfl
228            ihtest(jfl) = 0
229            ivtest(jfl) = 0
230            ikmfl(jfl) = 0
231# if   defined key_mpp_mpi
232            DO ji = MAX(nldi,2), nlei
233               DO jj = MAX(nldj,2), nlej   ! NO vector opt.
234# else
235            DO ji = 2, jpi
236               DO jj = 2, jpj   ! NO vector opt.
237# endif                 
238                  ! for each float we find the indexes of the mesh
239                 
240                  CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
241                                glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
242                                glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
243                                glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
244                                flxx(jfl)       ,flyy(jfl)       ,   &
245                                glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
246                  IF(llinmesh) THEN
247                     iimfl(jfl)  = ji
248                     ijmfl(jfl)  = jj
249                     ihtest(jfl) = ihtest(jfl)+1
250                     DO jk = 1, jpk-1
251                        IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) >  flzz(jfl)) ) THEN
252                           ikmfl(jfl)  = jk
253                           ivtest(jfl) = ivtest(jfl) + 1
254                        ENDIF
255                     END DO
256                  ENDIF
257               END DO
258            END DO
259           
260            ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1           
261            IF( ihtest(jfl) == 0 ) THEN
262               iimfl(jfl) = -1
263               ijmfl(jfl) = -1
264            ENDIF
265         END DO
266         
267         ! A zero in the sum of the arrays "ihtest" and "ivtest"         
268         IF( lk_mpp )   CALL mpp_sum(ihtest,jpnfl)   ! sums over the global domain
269         IF( lk_mpp )   CALL mpp_sum(ivtest,jpnfl)
270
271         DO jfl = 1, jpnfl
272            IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1 )) THEN
273               IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
274            ENDIF
275            IF( ihtest(jfl) == 0 ) THEN
276               IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH'
277            ENDIF
278         END DO
279       
280         ! We compute the distance between the float and the face of  the mesh         
281         DO jfl = 1, jpnfl
282            ! Made only if the float is in the domain of the processor
283            IF( (iimfl(jfl) >= 0 ) .AND. ( ijmfl(jfl) >= 0 ) ) THEN
284               
285               ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
286               
287               idomfl(jfl) = 0
288               IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl)=1
289               
290               ! Computation of the distance between the float
291               ! and the faces of the mesh
292               !            zdxab
293               !             .
294               !        B----.---------C
295               !        |    .         |
296               !        |<------>flo   |
297               !        |        ^     |
298               !        |        |.....|....zdyad
299               !        |        |     |
300               !        A--------|-----D
301               
302               zdxab = dstnce(flxx(jfl),flyy(jfl),glamf(iimfl(jfl)-1,ijmfl(jfl)-1),flyy(jfl))               
303               zdyad = dstnce(flxx(jfl),flyy(jfl),flxx(jfl),gphif(iimfl(jfl)-1,ijmfl(jfl)-1))
304               
305               ! Translation of this distances (in meter) in indexes
306               
307               tpifl(jfl) = (iimfl(jfl)-0.5)+zdxab/ e1u(iimfl(jfl)-1,ijmfl(jfl))+(mig(1)-jpizoom)
308               tpjfl(jfl) = (ijmfl(jfl)-0.5)+zdyad/ e2v(iimfl(jfl),ijmfl(jfl)-1)+(mjg(1)-jpjzoom)
309               tpkfl(jfl) = (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl))*(ikmfl(jfl))                     &
310                          / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))   &
311                          + (flzz(jfl) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))*(ikmfl(jfl)+1)                     &
312                          / (fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)))
313            ELSE
314               tpifl (jfl) = 0.e0
315               tpjfl (jfl) = 0.e0
316               tpkfl (jfl) = 0.e0
317               idomfl(jfl) = 0
318            ENDIF
319         END DO
320         
321         ! The sum of all the arrays tpifl, tpjfl, tpkfl give 3 arrays with the positions of all the floats.
322         IF( lk_mpp )   CALL mpp_sum( tpifl , jpnfl )   ! sums over the global domain
323         IF( lk_mpp )   CALL mpp_sum( tpjfl , jpnfl )
324         IF( lk_mpp )   CALL mpp_sum( tpkfl , jpnfl )
325         IF( lk_mpp )   CALL mpp_sum( idomfl, jpnfl )
326      ENDIF
327           
328      ! Print the initial positions of the floats
329      IF( .NOT. ln_rstflo ) THEN 
330         ! WARNING : initial position not in the sea         
331         DO jfl = 1, jpnfl
332            IF( idomfl(jfl) == 1 ) THEN
333               IF(lwp) WRITE(numout,*)'*****************************'
334               IF(lwp) WRITE(numout,*)'!!!!!!!  WARNING   !!!!!!!!!!'
335               IF(lwp) WRITE(numout,*)'*****************************'
336               IF(lwp) WRITE(numout,*)'The float number',jfl,'is out of the sea.'
337               IF(lwp) WRITE(numout,*)'geographical position',flxx(jfl),flyy(jfl),flzz(jfl)
338               IF(lwp) WRITE(numout,*)'index position',tpifl(jfl),tpjfl(jfl),tpkfl(jfl)
339            ENDIF
340         END DO
341      ENDIF
342
343   END SUBROUTINE flo_dom
344
345
346   SUBROUTINE findmesh( pax, pay, pbx, pby,   &
347                        pcx, pcy, pdx, pdy,   &
348                        px  ,py  ,ptx, pty, ldinmesh )
349      !! -------------------------------------------------------------
350      !!                ***  ROUTINE findmesh  ***
351      !!     
352      !! ** Purpose :   Find the index of mesh for the point spx spy
353      !!
354      !! ** Method  :
355      !!----------------------------------------------------------------------
356      REAL(wp) ::   &
357         pax, pay, pbx, pby,    &     ! ???
358         pcx, pcy, pdx, pdy,    &     ! ???
359         px, py,                &     ! longitude and latitude
360         ptx, pty                     ! ???
361      LOGICAL ::  ldinmesh            ! ???
362      !!
363      REAL(wp) ::   zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
364      !!---------------------------------------------------------------------
365      !! Statement function
366      REAL(wp) ::   fsline
367      REAL(wp) ::   psax, psay, psbx, psby, psx, psy
368      fsline( psax, psay, psbx, psby, psx, psy ) = psy  * ( psbx - psax )   &
369         &                                       - psx  * ( psby - psay )   &
370         &                                       + psax *   psby - psay * psbx
371      !!---------------------------------------------------------------------
372     
373      ! 4 semi plane defined by the 4 points and including the T point
374      zabt = fsline(pax,pay,pbx,pby,ptx,pty)
375      zbct = fsline(pbx,pby,pcx,pcy,ptx,pty)
376      zcdt = fsline(pcx,pcy,pdx,pdy,ptx,pty)
377      zdat = fsline(pdx,pdy,pax,pay,ptx,pty)
378     
379      ! 4 semi plane defined by the 4 points and including the extrememity
380      zabpt = fsline(pax,pay,pbx,pby,px,py)
381      zbcpt = fsline(pbx,pby,pcx,pcy,px,py)
382      zcdpt = fsline(pcx,pcy,pdx,pdy,px,py)
383      zdapt = fsline(pdx,pdy,pax,pay,px,py)
384       
385      ! We compare the semi plane T with the semi plane including the point
386      ! to know if it is in this  mesh.
387      ! For numerical reasons it is possible that for a point which is on
388      ! the line we don't have exactly zero with fsline function. We want
389      ! that a point can't be in 2 mesh in the same time, so we put the
390      ! coefficient to zero if it is smaller than 1.E-12
391     
392      IF( ABS(zabpt) <= 1.E-12 ) zabpt = 0.
393      IF( ABS(zbcpt) <= 1.E-12 ) zbcpt = 0.
394      IF( ABS(zcdpt) <= 1.E-12 ) zcdpt = 0.
395      IF( ABS(zdapt) <= 1.E-12 ) zdapt = 0.
396      IF( (zabt*zabpt >  0.) .AND. (zbct*zbcpt >= 0. ) .AND. ( zcdt*zcdpt >= 0. ) .AND. ( zdat*zdapt > 0. )   &
397         .AND. ( px <= MAX(pcx,pdx) ) .AND. ( px >= MIN(pax,pbx) )    &
398         .AND. ( py <= MAX(pby,pcy) ) .AND. ( py >= MIN(pay,pdy) ) ) THEN
399         ldinmesh=.TRUE.
400      ELSE
401         ldinmesh=.FALSE.
402      ENDIF
403      !
404   END SUBROUTINE findmesh
405
406
407   FUNCTION dstnce( pla1, phi1, pla2, phi2 )
408      !! -------------------------------------------------------------
409      !!                 ***  Function dstnce  ***
410      !!         
411      !! ** Purpose :   returns distance (in m) between two geographical
412      !!                points
413      !! ** Method  :
414      !!----------------------------------------------------------------------
415      REAL(wp), INTENT(in) ::   pla1, phi1, pla2, phi2   ! ???
416      !!
417      REAL(wp) ::   dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi
418      REAL(wp) ::   dstnce
419      !!---------------------------------------------------------------------
420      !
421      dpi  = 2.* ASIN(1.)
422      dls  = dpi / 180.
423      dly1 = phi1 * dls
424      dly2 = phi2 * dls
425      dlx1 = pla1 * dls
426      dlx2 = pla2 * dls
427      !
428      dlx = SIN(dly1) * SIN(dly2) + COS(dly1) * COS(dly2) * COS(dlx2-dlx1)
429      !
430      IF( ABS(dlx) > 1.0 ) dlx = 1.0
431      !
432      dld = ATAN(DSQRT( ( 1-dlx )/( 1+dlx ) )) * 222.24 / dls
433      dstnce = dld * 1000.
434      !
435   END FUNCTION dstnce
436
437#  else
438   !!----------------------------------------------------------------------
439   !!   Default option                                         Empty module
440   !!----------------------------------------------------------------------
441CONTAINS
442   SUBROUTINE flo_dom                 ! Empty routine
443   END SUBROUTINE flo_dom
444#endif
445
446   !!======================================================================
447END MODULE flodom
Note: See TracBrowser for help on using the repository browser.