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 trunk/NEMO/OPA_SRC/FLO – NEMO

source: trunk/NEMO/OPA_SRC/FLO/flodom.F90 @ 171

Last change on this file since 171 was 171, checked in by opalod, 19 years ago

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