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

Last change on this file since 13286 was 13286, checked in by smasson, 4 years ago

trunk: merge extra halos branch in trunk, see #2366

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