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 branches/UKMO/r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

source: branches/UKMO/r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 @ 7144

Last change on this file since 7144 was 7144, checked in by jcastill, 7 years ago

Remove svn keywords

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