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/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/FLO – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/FLO/flodom.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

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