source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/FLO/flodom.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 19 months ago

GMED 450 add flush after prints

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