source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/FLO/flodom.F90 @ 10420

Last change on this file since 10420 was 10420, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: force STOP when fail to allocate array, see #2133

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