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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/FLO/flodom.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

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