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 NEMO/trunk/src/OCE/FLO – NEMO

source: NEMO/trunk/src/OCE/FLO/flodom.F90

Last change on this file was 15235, checked in by clem, 3 years ago

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