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 @ 247

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

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