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/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/FLO – NEMO

source: NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/FLO/flodom.F90 @ 11692

Last change on this file since 11692 was 11692, checked in by francesca, 4 years ago

Update branch to integrate the development starting from the current v4.01 ready trunk

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