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

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

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