source: trunk/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 9 years ago

Merge of 3.4beta into the trunk

  • Property svn:keywords set to Id
File size: 20.3 KB
Line 
1MODULE flodom
2   !!======================================================================
3   !!                       ***  MODULE  flodom  ***
4   !! Ocean floats :   domain
5   !!======================================================================
6   !! History :  OPA          ! 1998-07 (Y.Drillet, CLIPPER)  Original code
7   !!            NEMO_3.3.1   ! 2011-09 (C.Bricaud,S.Law-Chune Mercator-Ocean):
8                              ! add Ariane convention, Comsecitc changes
9   !!----------------------------------------------------------------------
10#if   defined key_floats   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_floats'                                     float trajectories
13   !!----------------------------------------------------------------------
14   !!   flo_dom               : initialization of floats
15   !!   add_new_floats        : add new floats (long/lat/depth)
16   !!   add_new_ariane_floats : add new floats with araine convention (i/j/k)
17   !!   findmesh              : compute index of position
18   !!   dstnce                : compute distance between face mesh and floats
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE flo_oce         ! ocean drifting floats
23   USE in_out_manager  ! I/O manager
24   USE lib_mpp         ! distribued memory computing library
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   flo_dom         ! routine called by floats.F90
30   PUBLIC   flo_dom_alloc   ! Routine called in floats.F90
31
32   CHARACTER (len=21) ::  clname1 = 'init_float'              ! floats initialisation filename
33   CHARACTER (len=21) ::  clname2 = 'init_float_ariane'       ! ariane floats initialisation filename
34
35
36   INTEGER , ALLOCATABLE, DIMENSION(:) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats
37   INTEGER , ALLOCATABLE, DIMENSION(:) ::   idomfl, ivtest, ihtest    !   -     
38   REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zgifl, zgjfl,  zgkfl      ! distances in indexes
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE flo_dom
50      !! ---------------------------------------------------------------------
51      !!                  ***  ROUTINE flo_dom  ***
52      !!                 
53      !!  ** Purpose :   Initialisation of floats
54      !!
55      !!  ** Method  :   We put the floats  in the domain with the latitude,
56      !!               the longitude (degree) and the depth (m).
57      !!----------------------------------------------------------------------     
58      INTEGER            ::   jfl    ! dummy loop 
59      INTEGER            ::   inum   ! logical unit for file read
60      !!---------------------------------------------------------------------
61     
62      ! Initialisation with the geographical position or restart
63     
64      IF(lwp) WRITE(numout,*) 'flo_dom : compute initial position of floats'
65      IF(lwp) WRITE(numout,*) '~~~~~~~~'
66      IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl
67     
68      !-------------------------!
69      ! FLOAT RESTART FILE READ !
70      !-------------------------!
71      IF( ln_rstflo )THEN
72
73         IF(lwp) WRITE(numout,*) '        float restart file read'
74         
75         ! open the restart file
76         !----------------------
77         CALL ctl_opn( inum, 'restart_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
78
79         ! read of the restart file
80         READ(inum,*)   ( tpifl  (jfl), jfl=1, jpnrstflo),   & 
81                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   &
82                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   &
83                        ( nisobfl(jfl), jfl=1, jpnrstflo),   &
84                        ( ngrpfl (jfl), jfl=1, jpnrstflo)   
85         CLOSE(inum)
86
87         ! if we want a  surface drift  ( like PROVOR floats )
88         IF( ln_argo ) nisobfl(1:jpnrstflo) = 0
89         
90         ! It is possible to add new floats.         
91         !---------------------------------
92         IF( jpnfl > jpnrstflo )THEN
93
94            IF(lwp) WRITE(numout,*) '        add new floats'
95
96            IF( ln_ariane )THEN  !Add new floats with ariane convention
97                CALL flo_add_new_ariane_floats(jpnrstflo+1,jpnfl) 
98            ELSE                 !Add new floats with long/lat convention
99                CALL flo_add_new_floats(jpnrstflo+1,jpnfl)
100            ENDIF
101         ENDIF
102
103      !--------------------------------------!
104      ! FLOAT INITILISATION: NO RESTART FILE !
105      !--------------------------------------!
106      ELSE    !ln_rstflo
107
108         IF( ln_ariane )THEN       !Add new floats with ariane convention
109            CALL flo_add_new_ariane_floats(1,jpnfl)
110         ELSE                      !Add new floats with long/lat convention
111            CALL flo_add_new_floats(1,jpnfl)
112         ENDIF
113
114      ENDIF
115           
116   END SUBROUTINE flo_dom
117
118   SUBROUTINE flo_add_new_floats(kfl_start, kfl_end)
119      !! -------------------------------------------------------------
120      !!                 ***  SUBROUTINE add_new_arianefloats  ***
121      !!         
122      !! ** Purpose :   
123      !!
124      !!       First initialisation of floats
125      !!       the initials positions of floats are written in a file
126      !!       with a variable to know if it is a isobar float a number
127      !!       to identified who want the trajectories of this float and
128      !!       an index for the number of the float         
129      !!       open the init file
130      !!               
131      !! ** Method  :
132      !!----------------------------------------------------------------------
133      INTEGER, INTENT(in) :: kfl_start, kfl_end
134      !!
135      INTEGER           :: inum ! file unit
136      INTEGER           :: jfl,ji, jj, jk ! dummy loop indices
137      INTEGER           :: itrash         ! trash var for reading
138      INTEGER           :: ifl            ! number of floats to read
139      REAL(wp)          :: zdxab, zdyad
140      LOGICAL           :: llinmesh
141      CHARACTER(len=80) :: cltmp
142      !!---------------------------------------------------------------------
143      ifl = kfl_end-kfl_start+1
144
145      ! we get the init values
146      !-----------------------
147      CALL ctl_opn( inum , clname1, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
148      DO jfl = kfl_start,kfl_end
149         READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash
150         if(lwp)write(numout,*)'read:',jfl,flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash ; call flush(numout)
151      END DO
152      CLOSE(inum)
153           
154      ! Test to find the grid point coordonate with the geographical position           
155      !----------------------------------------------------------------------
156      DO jfl = kfl_start,kfl_end
157         ihtest(jfl) = 0
158         ivtest(jfl) = 0
159         ikmfl(jfl) = 0
160# if   defined key_mpp_mpi
161         DO ji = MAX(nldi,2), nlei
162            DO jj = MAX(nldj,2), nlej   ! NO vector opt.
163# else         
164         DO ji = 2, jpi
165            DO jj = 2, jpj   ! NO vector opt.
166# endif                     
167               ! For each float we find the indexes of the mesh                     
168               CALL flo_findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
169                                 glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
170                                 glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
171                                 glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
172                                 flxx(jfl)       ,flyy(jfl)       ,   &
173                                 glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
174               IF( llinmesh )THEN
175                  iimfl(jfl) = ji
176                  ijmfl(jfl) = jj
177                  ihtest(jfl) = ihtest(jfl)+1
178                  DO jk = 1, jpk-1
179                     IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN
180                        ikmfl(jfl) = jk
181                        ivtest(jfl) = ivtest(jfl) + 1
182                     ENDIF
183                  END DO
184               ENDIF
185            END DO
186         END DO
187
188         ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1               
189         IF( ihtest(jfl) ==  0 ) THEN
190            iimfl(jfl) = -1
191            ijmfl(jfl) = -1
192         ENDIF
193      END DO
194
195      !Test if each float is in one and only one proc
196      !----------------------------------------------
197      IF( lk_mpp )   THEN
198         CALL mpp_sum(ihtest,jpnfl)
199         CALL mpp_sum(ivtest,jpnfl)
200      ENDIF
201      DO jfl = kfl_start,kfl_end
202
203         IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
204             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
205             CALL ctl_stop('STOP',TRIM(cltmp) )
206         ENDIF
207         IF( (ihtest(jfl) == 0) ) THEN
208             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS IN NO MESH'
209             CALL ctl_stop('STOP',TRIM(cltmp) )
210         ENDIF
211      END DO
212
213      ! We compute the distance between the float and the face of the mesh           
214      !-------------------------------------------------------------------
215      DO jfl = kfl_start,kfl_end
216
217         ! Made only if the float is in the domain of the processor             
218         IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN
219
220            ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
221            idomfl(jfl) = 0
222            IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1
223
224            ! Computation of the distance between the float and the faces of the mesh
225            !            zdxab
226            !             .
227            !        B----.---------C
228            !        |    .         |
229            !        |<------>flo   |
230            !        |        ^     |
231            !        |        |.....|....zdyad
232            !        |        |     |
233            !        A--------|-----D
234            !
235            zdxab = flo_dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) )
236            zdyad = flo_dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) )
237
238            ! Translation of this distances (in meter) in indexes
239            zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
240            zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
241            zgkfl(jfl) = (( fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))   &
242               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
243               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) )                             &
244               &                 + (( flzz(jfl)-fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1))   &
245               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
246               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) )
247         ELSE
248            zgifl(jfl) = 0.e0
249            zgjfl(jfl) = 0.e0
250            zgkfl(jfl) = 0.e0
251         ENDIF
252
253      END DO
254                 
255      ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats.
256      IF( lk_mpp )   THEN
257         CALL mpp_sum( zgjfl, ifl )   ! sums over the global domain
258         CALL mpp_sum( zgkfl, ifl )
259      ENDIF
260           
261      DO jfl = kfl_start,kfl_end
262         tpifl(jfl) = zgifl(jfl)
263         tpjfl(jfl) = zgjfl(jfl)
264         tpkfl(jfl) = zgkfl(jfl)
265      END DO
266
267      ! WARNING : initial position not in the sea         
268      IF( .NOT. ln_rstflo ) THEN
269         DO jfl =  kfl_start,kfl_end
270            IF( idomfl(jfl) == 1 ) THEN
271               IF(lwp) WRITE(numout,*)'*****************************'
272               IF(lwp) WRITE(numout,*)'!!!!!!!  WARNING   !!!!!!!!!!'
273               IF(lwp) WRITE(numout,*)'*****************************'
274               IF(lwp) WRITE(numout,*)'The float number',jfl,'is out of the sea.'
275               IF(lwp) WRITE(numout,*)'geographical position',flxx(jfl),flyy(jfl),flzz(jfl)
276               IF(lwp) WRITE(numout,*)'index position',tpifl(jfl),tpjfl(jfl),tpkfl(jfl)
277            ENDIF
278         END DO
279      ENDIF
280
281   END SUBROUTINE flo_add_new_floats
282
283   SUBROUTINE flo_add_new_ariane_floats(kfl_start, kfl_end)
284      !! -------------------------------------------------------------
285      !!                 ***  SUBROUTINE add_new_arianefloats  ***
286      !!         
287      !! ** Purpose :   
288      !!       First initialisation of floats with ariane convention
289      !!       
290      !!       The indexes are read directly from file (warning ariane
291      !!       convention, are refered to
292      !!       U,V,W grids - and not T-)
293      !!       The isobar advection is managed with the sign of tpkfl ( >0 -> 3D
294      !!       advection, <0 -> 2D)
295      !!       Some variables are not read, as - gl         : time index; 4th
296      !!       column       
297      !!                                       - transport  : transport ; 5th
298      !!                                       column
299      !!       and paste in the jtrash var
300      !!       At the end, ones need to replace the indexes on T grid
301      !!       RMQ : there is no float groups identification !
302      !!
303      !!               
304      !! ** Method  :
305      !!----------------------------------------------------------------------
306      INTEGER, INTENT(in) :: kfl_start, kfl_end
307      !!
308      INTEGER  :: inum         ! file unit
309      INTEGER  :: ierr, ifl
310      INTEGER  :: jfl, jfl1    ! dummy loop indices
311      INTEGER  :: itrash       ! trash var for reading 
312      CHARACTER(len=80) :: cltmp
313
314      !!----------------------------------------------------------------------
315      nisobfl(kfl_start:kfl_end) = 1 ! we assume that by default we want 3D advection
316
317      ifl = kfl_end - kfl_start + 1  ! number of floats to read 
318
319      ! we check that the number of floats in the init_file are consistant with the namelist
320      IF( lwp ) THEN
321
322         jfl1=0
323         ierr=0
324         CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL',  1, numout, .TRUE., 1 )
325         DO WHILE (ierr .EQ. 0)
326            jfl1=jfl1+1
327            READ(inum,*, iostat=ierr)
328         END DO
329         CLOSE(inum)
330         IF( (jfl1-1) .NE. ifl )THEN
331            WRITE(cltmp,'(A25,A20,A3,i4.4,A10,i4.4)')"the number of floats in ",TRIM(clname2), &
332                                                     " = ",jfl1," is not equal to jfl= ",ifl
333            CALL ctl_stop('STOP',TRIM(cltmp) )
334         ENDIF
335
336      ENDIF
337           
338      ! we get the init values
339      CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .TRUE., 1 )
340      DO jfl = kfl_start, kfl_end
341          READ(inum,*) tpifl(jfl),tpjfl(jfl),tpkfl(jfl),itrash, itrash
342             
343          IF ( tpkfl(jfl) .LT. 0. ) nisobfl(jfl) = 0 !set the 2D advection according to init_float
344          ngrpfl(jfl)=jfl
345      END DO
346
347      ! conversion from ariane index to T grid index
348      tpkfl(kfl_start:kfl_end) = abs(tpkfl)-0.5 ! reversed vertical axis
349      tpifl(kfl_start:kfl_end) = tpifl+0.5
350      tpjfl(kfl_start:kfl_end) = tpjfl+0.5
351
352
353   END SUBROUTINE flo_add_new_ariane_floats
354
355
356   SUBROUTINE flo_findmesh( pax, pay, pbx, pby,   &
357                            pcx, pcy, pdx, pdy,   &
358                            px  ,py  ,ptx, pty, ldinmesh )
359      !! -------------------------------------------------------------
360      !!                ***  ROUTINE findmesh  ***
361      !!     
362      !! ** Purpose :   Find the index of mesh for the point spx spy
363      !!
364      !! ** Method  :
365      !!----------------------------------------------------------------------
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      REAL(wp) ::   zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
374      !!---------------------------------------------------------------------
375      !! Statement function
376      REAL(wp) ::   fsline
377      REAL(wp) ::   psax, psay, psbx, psby, psx, psy
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 flo_findmesh
415
416
417   FUNCTION flo_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      REAL(wp), INTENT(in) ::   pla1, phi1, pla2, phi2   ! ???
426      !!
427      REAL(wp) :: dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi
428      REAL(wp) :: flo_dstnce
429      !!---------------------------------------------------------------------
430      !
431      dpi  = 2._wp * ASIN(1._wp)
432      dls  = dpi / 180._wp
433      dly1 = phi1 * dls
434      dly2 = phi2 * dls
435      dlx1 = pla1 * dls
436      dlx2 = pla2 * dls
437      !
438      dlx = SIN(dly1) * SIN(dly2) + COS(dly1) * COS(dly2) * COS(dlx2-dlx1)
439      !
440      IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp
441      !
442      dld = ATAN(DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls
443      flo_dstnce = dld * 1000._wp
444      !
445   END FUNCTION flo_dstnce
446
447   INTEGER FUNCTION flo_dom_alloc()
448      !!----------------------------------------------------------------------
449      !!                 ***  FUNCTION flo_dom_alloc  ***
450      !!----------------------------------------------------------------------
451
452      ALLOCATE( iimfl(jpnfl) , ijmfl(jpnfl) , ikmfl(jpnfl) ,                          & 
453                idomfl(jpnfl), ivtest(jpnfl), ihtest(jpnfl),                 &
454                zgifl(jpnfl) , zgjfl(jpnfl) , zgkfl(jpnfl)   , STAT=flo_dom_alloc )
455      !
456      IF( lk_mpp             )   CALL mpp_sum ( flo_dom_alloc )
457      IF( flo_dom_alloc /= 0 )   CALL ctl_warn('flo_dom_alloc: failed to allocate arrays')
458   END FUNCTION flo_dom_alloc
459
460
461#else
462   !!----------------------------------------------------------------------
463   !!   Default option                                         Empty module
464   !!----------------------------------------------------------------------
465CONTAINS
466   SUBROUTINE flo_dom                 ! Empty routine
467         WRITE(*,*) 'flo_dom: : You should not have seen this print! error?'
468   END SUBROUTINE flo_dom
469#endif
470
471   !!======================================================================
472END MODULE flodom
Note: See TracBrowser for help on using the repository browser.