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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.6 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 flo_oce         ! ocean drifting floats
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE lib_mpp         ! distribued memory computing library
19
20   IMPLICIT NONE
21
22   !! * Accessibility
23   PRIVATE  dstnce
24   PUBLIC flo_dom     ! routine called by floats.F90
25
26   !! * Substitutions
27#  include "domzgr_substitute.h90"
28   !!----------------------------------------------------------------------
29   !!   OPA 9.0 , LODYC-IPSL  (2003)
30   !!----------------------------------------------------------------------
31
32CONTAINS
33
34   SUBROUTINE flo_dom
35      !! ---------------------------------------------------------------------
36      !!                  ***  ROUTINE flo_dom  ***
37      !!                 
38      !!  ** Purpose :   Initialisation of floats
39      !!
40      !!  ** Method  :   We put the floats  in the domain with the latitude,
41      !!       the longitude (degree) and the depth (m).
42      !!
43      !!----------------------------------------------------------------------     
44      !! * Local declarations
45      LOGICAL   :: llinmesh
46      CHARACTER (len=21) ::  clname
47      INTEGER   :: ji, jj, jk               ! DO loop index on 3 directions
48      INTEGER   :: jfl, jfl1                ! number of floats   
49      INTEGER   :: inum = 11                ! logical unit for file read
50      INTEGER , DIMENSION ( jpnfl    )  ::   &
51         iimfl, ijmfl, ikmfl,    &          ! index mesh of floats
52         idomfl,  ivtest, ihtest
53      REAL(wp)  :: zdxab,zdyad
54      REAL(wp) , DIMENSION ( jpnnewfl )  :: zgifl, zgjfl,  zgkfl
55      !!---------------------------------------------------------------------
56     
57      ! Initialisation with the geographical position or restart
58     
59      IF(lwp) WRITE(numout,*) 'flo_dom : compute initial position of floats'
60      IF(lwp) WRITE(numout,*) '~~~~~~~~'
61      IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl
62     
63      IF(ln_rstarfl) THEN
64         IF(lwp) WRITE(numout,*) '        float restart file read'
65         
66         ! open the restart file
67         clname='restart_float'
68         OPEN (inum,FILE=clname,FORM='UNFORMATTED')
69         REWIND inum
70
71         ! read of the restart file
72         READ(inum) ( tpifl  (jfl), jfl=1, jpnrstarfl),   & 
73                        ( tpjfl  (jfl), jfl=1, jpnrstarfl),   &
74                        ( tpkfl  (jfl), jfl=1, jpnrstarfl),   &
75                        ( nisobfl(jfl), jfl=1, jpnrstarfl),   &
76                        ( ngrpfl (jfl), jfl=1, jpnrstarfl)   
77         CLOSE(inum)
78
79         ! if we want a  surface drift  ( like PROVOR floats )
80         IF( ln_argo ) THEN
81            DO jfl = 1, jpnrstarfl
82               nisobfl(jfl) = 0
83            END DO
84         ENDIF
85
86         IF(lwp) WRITE(numout,*)' flo_dom: END of florstlec'
87         
88         ! It is possible to add new floats.         
89         IF(lwp) WRITE(numout,*)' flo_dom:jpnfl jpnrstarfl ',jpnfl,jpnrstarfl
90         IF( jpnfl > jpnrstarfl ) THEN
91            ! open the init file
92            clname='init_float'
93            OPEN(inum,FILE=clname,FORM='FORMATTED')
94            DO jfl = jpnrstarfl+1, jpnfl
95               READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),jfl1
96            END DO
97            CLOSE(inum)
98            IF(lwp) WRITE(numout,*)' flodom: END reading init_float file'
99           
100            ! Test to find the grid point coordonate with the geographical position           
101            DO jfl = jpnrstarfl+1, jpnfl
102               ihtest(jfl) = 0
103               ivtest(jfl) = 0
104               ikmfl(jfl) = 0
105# if defined key_mpp
106               DO ji = MAX(nldi,2), nlei
107                  DO jj = MAX(nldj,2), nlej
108# else
109               DO ji = 2, jpi
110                  DO jj = 2, jpj
111# endif                     
112                     ! For each float we find the indexes of the mesh                     
113                     CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
114                                   glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
115                                   glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
116                                   glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
117                                   flxx(jfl)       ,flyy(jfl)       ,   &
118                                   glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
119                     IF(llinmesh) THEN
120                        iimfl(jfl) = ji
121                        ijmfl(jfl) = jj
122                        ihtest(jfl) = ihtest(jfl)+1
123                        DO jk = 1, jpk-1
124                           IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN
125                              ikmfl(jfl) = jk
126                              ivtest(jfl) = ivtest(jfl) + 1
127                           ENDIF
128                        END DO
129                     ENDIF
130                  END DO
131               END DO
132               IF(lwp) WRITE(numout,*)'   flo_dom: END findmesh'
133               
134               ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1               
135               IF( ihtest(jfl) ==  0 ) THEN
136                  iimfl(jfl) = -1
137                  ijmfl(jfl) = -1
138               ENDIF
139            END DO
140           
141            ! A zero in the sum of the arrays "ihtest" and "ivtest"             
142# if defined key_mpp
143            CALL mpp_sum(ihtest,jpnfl,iwork)
144            CALL mpp_sum(ivtest,jpnfl,iwork)
145# endif
146            DO jfl = jpnrstarfl+1, jpnfl
147               IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
148                  IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
149                  STOP
150               ENDIF
151               IF( (ihtest(jfl) == 0) ) THEN
152                  IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH'
153                  STOP
154               ENDIF
155            END DO
156           
157            ! We compute the distance between the float and the face of the mesh           
158            DO jfl = jpnrstarfl+1, jpnfl               
159               ! Made only if the float is in the domain of the processor             
160               IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN
161                 
162                  ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
163                 
164                  idomfl(jfl) = 0
165                  IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1
166                                           
167                  ! Computation of the distance between the float and the faces of the mesh
168                  !            zdxab
169                  !             .
170                  !        B----.---------C
171                  !        |    .         |
172                  !        |<------>flo   |
173                  !        |        ^     |
174                  !        |        |.....|....zdyad
175                  !        |        |     |
176                  !        A--------|-----D
177                  !
178             
179                  zdxab = dstnce(flxx(jfl),flyy(jfl),glamf(iimfl(jfl)-1,ijmfl(jfl)-1),flyy(jfl))
180                  zdyad = dstnce(flxx(jfl),flyy(jfl),flxx(jfl),gphif(iimfl(jfl)-1,ijmfl(jfl)-1))
181                 
182                  ! Translation of this distances (in meter) in indexes
183                 
184                  zgifl(jfl-jpnrstarfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
185                  zgjfl(jfl-jpnrstarfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
186                  zgkfl(jfl-jpnrstarfl) = (( fsdepw(ji,jj,ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))  &
187                                        / (  fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl) ) )   &
188                                        + (( flzz(jfl)-fsdepw(ji,jj,ikmfl(jfl)) ) *(ikmfl(jfl)+1))   &
189                                        / (  fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl)) ) 
190               ELSE
191                  zgifl(jfl-jpnrstarfl) = 0.
192                  zgjfl(jfl-jpnrstarfl) = 0.
193                  zgkfl(jfl-jpnrstarfl) = 0.
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.
198# if defined key_mpp
199
200            CALL mpp_sum( zgjfl, jpnnewfl )
201            CALL mpp_sum( zgkfl, jpnnewfl )
202            IF(lwp) WRITE(numout,*) (zgifl(jfl),jfl=1,jpnnewfl)
203            IF(lwp) WRITE(numout,*) (zgjfl(jfl),jfl=1,jpnnewfl)
204            IF(lwp) WRITE(numout,*) (zgkfl(jfl),jfl=1,jpnnewfl) 
205# endif
206           
207            DO jfl = jpnrstarfl+1, jpnfl
208               tpifl(jfl) = zgifl(jfl-jpnrstarfl)
209               tpjfl(jfl) = zgjfl(jfl-jpnrstarfl)
210               tpkfl(jfl) = zgkfl(jfl-jpnrstarfl)
211            END DO
212         ENDIF
213      ELSE
214         IF(lwp) WRITE(numout,*) '                     init_float read '
215         
216         ! First initialisation of floats
217         ! the initials positions of floats are written in a file
218         ! with a variable to know if it is a isobar float a number
219         ! to identified who want the trajectories of this float and
220         ! an index for the number of the float         
221         ! open the init file
222         clname='init_float'
223         OPEN(inum,FILE=clname,FORM='FORMATTED')
224         READ(inum) (flxx(jfl)   , jfl=1, jpnfl),   &
225                    (flyy(jfl)   , jfl=1, jpnfl),   &
226                    (flzz(jfl)   , jfl=1, jpnfl),   &
227                    (nisobfl(jfl), jfl=1, jpnfl),   &
228                    (ngrpfl(jfl) , jfl=1, jpnfl)
229         CLOSE(inum)
230           
231         ! Test to find the grid point coordonate with the geographical position         
232         DO jfl = 1, jpnfl
233            ihtest(jfl) = 0
234            ivtest(jfl) = 0
235            ikmfl(jfl) = 0
236# if defined key_mpp
237            DO ji = MAX(nldi,2), nlei
238               DO jj = MAX(nldj,2), nlej
239# else
240            DO ji = 2, jpi
241               DO jj = 2, jpj
242# endif                 
243                  ! for each float we find the indexes of the mesh
244                 
245                  CALL findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
246                                glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
247                                glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
248                                glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
249                                flxx(jfl)       ,flyy(jfl)       ,   &
250                                glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
251                  IF(llinmesh) THEN
252                     iimfl(jfl)  = ji
253                     ijmfl(jfl)  = jj
254                     ihtest(jfl) = ihtest(jfl)+1
255                     DO jk = 1, jpk-1
256                        IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) >  flzz(jfl)) ) THEN
257                           ikmfl(jfl)  = jk
258                           ivtest(jfl) = ivtest(jfl) + 1
259                        ENDIF
260                     END DO
261                  ENDIF
262               END DO
263            END DO
264           
265            ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1           
266            IF( ihtest(jfl) == 0 ) THEN
267               iimfl(jfl) = -1
268               ijmfl(jfl) = -1
269            ENDIF
270         END DO
271         
272         ! A zero in the sum of the arrays "ihtest" and "ivtest"         
273# if defined key_mpp
274         CALL mpp_sum(ihtest,jpnfl,iwork)
275         CALL mpp_sum(ivtest,jpnfl,iwork)
276# endif
277
278         DO jfl = 1, jpnfl
279            IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1 )) THEN
280               IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
281            ENDIF
282            IF( ihtest(jfl) == 0 ) THEN
283               IF(lwp) WRITE(numout,*)'THE FLOAT',jfl,' IS IN NO MESH'
284            ENDIF
285         END DO
286       
287         ! We compute the distance between the float and the face of  the mesh         
288         DO jfl = 1, jpnfl
289            ! Made only if the float is in the domain of the processor
290            IF( (iimfl(jfl) >= 0 ) .AND. ( ijmfl(jfl) >= 0 ) ) THEN
291               
292               ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
293               
294               idomfl(jfl) = 0
295               IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl)=1
296               
297               ! Computation of the distance between the float
298               ! and the faces of the mesh
299               !            zdxab
300               !             .
301               !        B----.---------C
302               !        |    .         |
303               !        |<------>flo   |
304               !        |        ^     |
305               !        |        |.....|....zdyad
306               !        |        |     |
307               !        A--------|-----D
308               
309               zdxab = dstnce(flxx(jfl),flyy(jfl),glamf(iimfl(jfl)-1,ijmfl(jfl)-1),flyy(jfl))               
310               zdyad = dstnce(flxx(jfl),flyy(jfl),flxx(jfl),gphif(iimfl(jfl)-1,ijmfl(jfl)-1))
311               
312               ! Translation of this distances (in meter) in indexes
313               
314               tpifl(jfl) = (iimfl(jfl)-0.5)+zdxab/ e1u(iimfl(jfl)-1,ijmfl(jfl))+(mig(1)-jpizoom)
315               tpjfl(jfl) = (ijmfl(jfl)-0.5)+zdyad/ e2v(iimfl(jfl),ijmfl(jfl)-1)+(mjg(1)-jpjzoom)
316               tpkfl(jfl) = (fsdepw(ji,jj,ikmfl(jfl)+1) - flzz(jfl))*(ikmfl(jfl))   &
317                          / (fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl)))   &
318                          + (flzz(jfl) - fsdepw(ji,jj,ikmfl(jfl)))*(ikmfl(jfl)+1)   &
319                          / (fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl)))
320            ELSE
321               tpifl (jfl) = 0.e0
322               tpjfl (jfl) = 0.e0
323               tpkfl (jfl) = 0.e0
324               idomfl(jfl) = 0
325            ENDIF
326         END DO
327         
328         ! The sum of all the arrays tpifl, tpjfl, tpkfl give 3 arrays with the positions of all the floats.
329# if defined key_mpp
330         CALL mpp_sum( tpifl , jpnfl )
331         CALL mpp_sum( tpjfl , jpnfl )
332         CALL mpp_sum( tpkfl , jpnfl )
333         CALL mpp_sum( idomfl, jpnfl )
334# endif
335      ENDIF
336           
337      ! Print the initial positions of the floats
338      IF( .NOT. ln_rstarfl ) 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.