ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/FLO/flodom.F90

MODULE flodom
   !!======================================================================
   !!                       ***  MODULE  flodom  ***
   !! Ocean floats :   domain
   !!======================================================================
   !! History :  OPA          ! 1998-07 (Y.Drillet, CLIPPER)  Original code
   !!            NEMO_3.3.1   ! 2011-09 (C.Bricaud,S.Law-Chune Mercator-Ocean): 
                              ! add Ariane convention, Comsecitc changes
   !!----------------------------------------------------------------------
#if   defined key_floats   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_floats'                                     float trajectories
   !!----------------------------------------------------------------------
   !!   flo_dom               : initialization of floats
   !!   add_new_floats        : add new floats (long/lat/depth)
   !!   add_new_ariane_floats : add new floats with araine convention (i/j/k)
   !!   findmesh              : compute index of position 
   !!   dstnce                : compute distance between face mesh and floats 
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE flo_oce         ! ocean drifting floats
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! distribued memory computing library

   IMPLICIT NONE
   PRIVATE

   PUBLIC   flo_dom         ! routine called by floats.F90
   PUBLIC   flo_dom_alloc   ! Routine called in floats.F90

   CHARACTER (len=21) ::  clname1 = 'init_float'              ! floats initialisation filename
   CHARACTER (len=21) ::  clname2 = 'init_float_ariane'       ! ariane floats initialisation filename


   INTEGER , ALLOCATABLE, DIMENSION(:) ::   iimfl, ijmfl, ikmfl       ! index mesh of floats
   INTEGER , ALLOCATABLE, DIMENSION(:) ::   idomfl, ivtest, ihtest    !   -     
   REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zgifl, zgjfl,  zgkfl      ! distances in indexes

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: flodom.F90 3294 2012-01-28 16:44:18Z rblod $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE flo_dom
      !! ---------------------------------------------------------------------
      !!                  ***  ROUTINE flo_dom  ***
      !!                 
      !!  ** Purpose :   Initialisation of floats
      !!
      !!  ** Method  :   We put the floats  in the domain with the latitude,
      !!               the longitude (degree) and the depth (m).
      !!----------------------------------------------------------------------      
      INTEGER            ::   jfl    ! dummy loop  
      INTEGER            ::   inum   ! logical unit for file read
      !!---------------------------------------------------------------------
      
      ! Initialisation with the geographical position or restart
      
      IF(lwp) WRITE(numout,*) 'flo_dom : compute initial position of floats'
      IF(lwp) WRITE(numout,*) '~~~~~~~~'
      IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl
      
      !-------------------------!
      ! FLOAT RESTART FILE READ !
      !-------------------------!
      IF( ln_rstflo )THEN

         IF(lwp) WRITE(numout,*) '        float restart file read'
         
         ! open the restart file 
         !----------------------
         CALL ctl_opn( inum, 'restart_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )

         ! read of the restart file
         READ(inum,*)   ( tpifl  (jfl), jfl=1, jpnrstflo),   & 
                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   &
                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   &
                        ( nisobfl(jfl), jfl=1, jpnrstflo),   &
                        ( ngrpfl (jfl), jfl=1, jpnrstflo)    
         CLOSE(inum)

         ! if we want a  surface drift  ( like PROVOR floats )
         IF( ln_argo ) nisobfl(1:jpnrstflo) = 0
         
         ! It is possible to add new floats.          
         !---------------------------------
         IF( jpnfl > jpnrstflo )THEN

            IF(lwp) WRITE(numout,*) '        add new floats'

            IF( ln_ariane )THEN  !Add new floats with ariane convention
                CALL flo_add_new_ariane_floats(jpnrstflo+1,jpnfl) 
            ELSE                 !Add new floats with long/lat convention
                CALL flo_add_new_floats(jpnrstflo+1,jpnfl)
            ENDIF
         ENDIF

      !--------------------------------------!
      ! FLOAT INITILISATION: NO RESTART FILE !
      !--------------------------------------!
      ELSE    !ln_rstflo

         IF( ln_ariane )THEN       !Add new floats with ariane convention
            CALL flo_add_new_ariane_floats(1,jpnfl)
         ELSE                      !Add new floats with long/lat convention
            CALL flo_add_new_floats(1,jpnfl)
         ENDIF

      ENDIF
            
   END SUBROUTINE flo_dom

   SUBROUTINE flo_add_new_floats(kfl_start, kfl_end)
      !! -------------------------------------------------------------
      !!                 ***  SUBROUTINE add_new_arianefloats  ***
      !!          
      !! ** Purpose :   
      !!
      !!       First initialisation of floats
      !!       the initials positions of floats are written in a file
      !!       with a variable to know if it is a isobar float a number 
      !!       to identified who want the trajectories of this float and 
      !!       an index for the number of the float         
      !!       open the init file 
      !!               
      !! ** Method  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) :: kfl_start, kfl_end
      !!
      INTEGER           :: inum ! file unit
      INTEGER           :: jfl,ji, jj, jk ! dummy loop indices
      INTEGER           :: itrash         ! trash var for reading
      INTEGER           :: ifl            ! number of floats to read
      REAL(wp)          :: zdxab, zdyad
      LOGICAL           :: llinmesh
      CHARACTER(len=80) :: cltmp
      !!---------------------------------------------------------------------
      ifl = kfl_end-kfl_start+1

      ! we get the init values 
      !-----------------------
      CALL ctl_opn( inum , clname1, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
      DO jfl = kfl_start,kfl_end
         READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash
         if(lwp)write(numout,*)'read:',jfl,flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash ; call flush(numout)
      END DO
      CLOSE(inum)
            
      ! Test to find the grid point coordonate with the geographical position            
      !----------------------------------------------------------------------
      DO jfl = kfl_start,kfl_end
         ihtest(jfl) = 0
         ivtest(jfl) = 0
         ikmfl(jfl) = 0
# if   defined key_mpp_mpi
         DO ji = MAX(nldi,2), nlei
            DO jj = MAX(nldj,2), nlej   ! NO vector opt.
# else         
         DO ji = 2, jpi
            DO jj = 2, jpj   ! NO vector opt.
# endif                     
               ! For each float we find the indexes of the mesh                      
               CALL flo_findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1),   &
                                 glamf(ji-1,jj  ),gphif(ji-1,jj  ),   &
                                 glamf(ji  ,jj  ),gphif(ji  ,jj  ),   &
                                 glamf(ji  ,jj-1),gphif(ji  ,jj-1),   &
                                 flxx(jfl)       ,flyy(jfl)       ,   &
                                 glamt(ji  ,jj  ),gphit(ji  ,jj  ), llinmesh)
               IF( llinmesh )THEN
                  iimfl(jfl) = ji
                  ijmfl(jfl) = jj
                  ihtest(jfl) = ihtest(jfl)+1
                  DO jk = 1, jpk-1
                     IF( (fsdepw(ji,jj,jk) <= flzz(jfl)) .AND. (fsdepw(ji,jj,jk+1) > flzz(jfl)) ) THEN
                        ikmfl(jfl) = jk
                        ivtest(jfl) = ivtest(jfl) + 1
                     ENDIF
                  END DO
               ENDIF
            END DO
         END DO

         ! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1               
         IF( ihtest(jfl) ==  0 ) THEN
            iimfl(jfl) = -1
            ijmfl(jfl) = -1
         ENDIF
      END DO

      !Test if each float is in one and only one proc
      !----------------------------------------------
      IF( lk_mpp )   THEN 
         CALL mpp_sum(ihtest,jpnfl)
         CALL mpp_sum(ivtest,jpnfl)
      ENDIF
      DO jfl = kfl_start,kfl_end

         IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
             CALL ctl_stop('STOP',TRIM(cltmp) )
         ENDIF
         IF( (ihtest(jfl) == 0) ) THEN
             WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS IN NO MESH'
             CALL ctl_stop('STOP',TRIM(cltmp) )
         ENDIF
      END DO

      ! We compute the distance between the float and the face of the mesh            
      !-------------------------------------------------------------------
      DO jfl = kfl_start,kfl_end

         ! Made only if the float is in the domain of the processor              
         IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN

            ! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
            idomfl(jfl) = 0
            IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1

            ! Computation of the distance between the float and the faces of the mesh
            !            zdxab
            !             .
            !        B----.---------C
            !        |    .         |
            !        |<------>flo   |
            !        |        ^     |
            !        |        |.....|....zdyad
            !        |        |     |
            !        A--------|-----D
            !
            zdxab = flo_dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) )
            zdyad = flo_dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) )

            ! Translation of this distances (in meter) in indexes
            zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
            zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
            zgkfl(jfl) = (( fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))   &
               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) )                             &
               &                 + (( flzz(jfl)-fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1))   &
               &                 / (  fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1)                              &
               &                    - fsdepw(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) )
         ELSE
            zgifl(jfl) = 0.e0
            zgjfl(jfl) = 0.e0
            zgkfl(jfl) = 0.e0
         ENDIF

      END DO
                  
      ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats.
      IF( lk_mpp )   THEN 
         CALL mpp_sum( zgjfl, ifl )   ! sums over the global domain
         CALL mpp_sum( zgkfl, ifl )
      ENDIF
            
      DO jfl = kfl_start,kfl_end
         tpifl(jfl) = zgifl(jfl)
         tpjfl(jfl) = zgjfl(jfl)
         tpkfl(jfl) = zgkfl(jfl)
      END DO

      ! WARNING : initial position not in the sea         
      IF( .NOT. ln_rstflo ) THEN 
         DO jfl =  kfl_start,kfl_end
            IF( idomfl(jfl) == 1 ) THEN
               IF(lwp) WRITE(numout,*)'*****************************'
               IF(lwp) WRITE(numout,*)'!!!!!!!  WARNING   !!!!!!!!!!'
               IF(lwp) WRITE(numout,*)'*****************************'
               IF(lwp) WRITE(numout,*)'The float number',jfl,'is out of the sea.'
               IF(lwp) WRITE(numout,*)'geographical position',flxx(jfl),flyy(jfl),flzz(jfl)
               IF(lwp) WRITE(numout,*)'index position',tpifl(jfl),tpjfl(jfl),tpkfl(jfl)
            ENDIF
         END DO
      ENDIF

   END SUBROUTINE flo_add_new_floats

   SUBROUTINE flo_add_new_ariane_floats(kfl_start, kfl_end)
      !! -------------------------------------------------------------
      !!                 ***  SUBROUTINE add_new_arianefloats  ***
      !!          
      !! ** Purpose :   
      !!       First initialisation of floats with ariane convention
      !!       
      !!       The indexes are read directly from file (warning ariane
      !!       convention, are refered to 
      !!       U,V,W grids - and not T-) 
      !!       The isobar advection is managed with the sign of tpkfl ( >0 -> 3D
      !!       advection, <0 -> 2D) 
      !!       Some variables are not read, as - gl         : time index; 4th
      !!       column        
      !!                                       - transport  : transport ; 5th
      !!                                       column
      !!       and paste in the jtrash var
      !!       At the end, ones need to replace the indexes on T grid
      !!       RMQ : there is no float groups identification !
      !!
      !!               
      !! ** Method  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) :: kfl_start, kfl_end
      !!
      INTEGER  :: inum         ! file unit
      INTEGER  :: ierr, ifl
      INTEGER  :: jfl, jfl1    ! dummy loop indices
      INTEGER  :: itrash       ! trash var for reading  
      CHARACTER(len=80) :: cltmp

      !!----------------------------------------------------------------------
      nisobfl(kfl_start:kfl_end) = 1 ! we assume that by default we want 3D advection

      ifl = kfl_end - kfl_start + 1  ! number of floats to read  

      ! we check that the number of floats in the init_file are consistant with the namelist
      IF( lwp ) THEN

         jfl1=0
         ierr=0
         CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL',  1, numout, .TRUE., 1 )
         DO WHILE (ierr .EQ. 0)
            jfl1=jfl1+1
            READ(inum,*, iostat=ierr)
         END DO
         CLOSE(inum)
         IF( (jfl1-1) .NE. ifl )THEN 
            WRITE(cltmp,'(A25,A20,A3,i4.4,A10,i4.4)')"the number of floats in ",TRIM(clname2), &
                                                     " = ",jfl1," is not equal to jfl= ",ifl
            CALL ctl_stop('STOP',TRIM(cltmp) )
         ENDIF

      ENDIF
            
      ! we get the init values 
      CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .TRUE., 1 )
      DO jfl = kfl_start, kfl_end
          READ(inum,*) tpifl(jfl),tpjfl(jfl),tpkfl(jfl),itrash, itrash
              
          IF ( tpkfl(jfl) .LT. 0. ) nisobfl(jfl) = 0 !set the 2D advection according to init_float
          ngrpfl(jfl)=jfl
      END DO

      ! conversion from ariane index to T grid index
      tpkfl(kfl_start:kfl_end) = abs(tpkfl)-0.5 ! reversed vertical axis
      tpifl(kfl_start:kfl_end) = tpifl+0.5
      tpjfl(kfl_start:kfl_end) = tpjfl+0.5


   END SUBROUTINE flo_add_new_ariane_floats


   SUBROUTINE flo_findmesh( pax, pay, pbx, pby,   &
                            pcx, pcy, pdx, pdy,   &
                            px  ,py  ,ptx, pty, ldinmesh )
      !! -------------------------------------------------------------
      !!                ***  ROUTINE findmesh  ***
      !!      
      !! ** Purpose :   Find the index of mesh for the point spx spy
      !!
      !! ** Method  : 
      !!----------------------------------------------------------------------
      REAL(wp) ::   &
         pax, pay, pbx, pby,    &     ! ???
         pcx, pcy, pdx, pdy,    &     ! ???
         px, py,                &     ! longitude and latitude
         ptx, pty                     ! ???
      LOGICAL ::  ldinmesh            ! ???
      !!
      REAL(wp) ::   zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
      !!---------------------------------------------------------------------
      !! Statement function
      REAL(wp) ::   fsline
      REAL(wp) ::   psax, psay, psbx, psby, psx, psy
      fsline( psax, psay, psbx, psby, psx, psy ) = psy  * ( psbx - psax )   &
         &                                       - psx  * ( psby - psay )   &
         &                                       + psax *   psby - psay * psbx
      !!---------------------------------------------------------------------
      
      ! 4 semi plane defined by the 4 points and including the T point
      zabt = fsline(pax,pay,pbx,pby,ptx,pty)
      zbct = fsline(pbx,pby,pcx,pcy,ptx,pty)
      zcdt = fsline(pcx,pcy,pdx,pdy,ptx,pty)
      zdat = fsline(pdx,pdy,pax,pay,ptx,pty)
      
      ! 4 semi plane defined by the 4 points and including the extrememity
      zabpt = fsline(pax,pay,pbx,pby,px,py)
      zbcpt = fsline(pbx,pby,pcx,pcy,px,py)
      zcdpt = fsline(pcx,pcy,pdx,pdy,px,py)
      zdapt = fsline(pdx,pdy,pax,pay,px,py)
       
      ! We compare the semi plane T with the semi plane including the point
      ! to know if it is in this  mesh.
      ! For numerical reasons it is possible that for a point which is on
      ! the line we don't have exactly zero with fsline function. We want 
      ! that a point can't be in 2 mesh in the same time, so we put the 
      ! coefficient to zero if it is smaller than 1.E-12
      
      IF( ABS(zabpt) <= 1.E-12 ) zabpt = 0.
      IF( ABS(zbcpt) <= 1.E-12 ) zbcpt = 0.
      IF( ABS(zcdpt) <= 1.E-12 ) zcdpt = 0.
      IF( ABS(zdapt) <= 1.E-12 ) zdapt = 0.
      IF( (zabt*zabpt >  0.) .AND. (zbct*zbcpt >= 0. ) .AND. ( zcdt*zcdpt >= 0. ) .AND. ( zdat*zdapt > 0. )   &
         .AND. ( px <= MAX(pcx,pdx) ) .AND. ( px >= MIN(pax,pbx) )    &
         .AND. ( py <= MAX(pby,pcy) ) .AND. ( py >= MIN(pay,pdy) ) ) THEN
         ldinmesh=.TRUE.
      ELSE
         ldinmesh=.FALSE.
      ENDIF
      !
   END SUBROUTINE flo_findmesh


   FUNCTION flo_dstnce( pla1, phi1, pla2, phi2 )
      !! -------------------------------------------------------------
      !!                 ***  Function dstnce  ***
      !!          
      !! ** Purpose :   returns distance (in m) between two geographical
      !!                points
      !! ** Method  : 
      !!----------------------------------------------------------------------
      REAL(wp), INTENT(in) ::   pla1, phi1, pla2, phi2   ! ???
      !!
      REAL(wp) :: dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi
      REAL(wp) :: flo_dstnce
      !!---------------------------------------------------------------------
      !
      dpi  = 2._wp * ASIN(1._wp)
      dls  = dpi / 180._wp
      dly1 = phi1 * dls
      dly2 = phi2 * dls
      dlx1 = pla1 * dls
      dlx2 = pla2 * dls
      !
      dlx = SIN(dly1) * SIN(dly2) + COS(dly1) * COS(dly2) * COS(dlx2-dlx1)
      !
      IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp
      !
      dld = ATAN(DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls
      flo_dstnce = dld * 1000._wp
      !
   END FUNCTION flo_dstnce

   INTEGER FUNCTION flo_dom_alloc()
      !!----------------------------------------------------------------------
      !!                 ***  FUNCTION flo_dom_alloc  ***
      !!----------------------------------------------------------------------

      ALLOCATE( iimfl(jpnfl) , ijmfl(jpnfl) , ikmfl(jpnfl) ,                          &  
                idomfl(jpnfl), ivtest(jpnfl), ihtest(jpnfl),                 &
                zgifl(jpnfl) , zgjfl(jpnfl) , zgkfl(jpnfl)   , STAT=flo_dom_alloc )
      !
      IF( lk_mpp             )   CALL mpp_sum ( flo_dom_alloc )
      IF( flo_dom_alloc /= 0 )   CALL ctl_warn('flo_dom_alloc: failed to allocate arrays')
   END FUNCTION flo_dom_alloc


#else
   !!----------------------------------------------------------------------
   !!   Default option                                         Empty module
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE flo_dom                 ! Empty routine
         WRITE(*,*) 'flo_dom: : You should not have seen this print! error?'
   END SUBROUTINE flo_dom
#endif

   !!======================================================================
END MODULE flodom