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/2011/dev_r2802_MERCATOR9_floats/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: branches/2011/dev_r2802_MERCATOR9_floats/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 @ 2839

Last change on this file since 2839 was 2839, checked in by cbricaud, 13 years ago

modified routine for netcdf output

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