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/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

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

Last change on this file since 10892 was 10892, checked in by andmirek, 3 years ago

GMED 450 reviewer comments

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