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

MODULE floblk
   !!======================================================================
   !!                     ***  MODULE  floblk  ***
   !! Ocean floats :   trajectory computation
   !!======================================================================
#if   defined key_floats   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_floats'                                     float trajectories
   !!----------------------------------------------------------------------
   !!    flotblk     : compute float trajectories with Blanke algorithme
   !!----------------------------------------------------------------------
   USE flo_oce         ! ocean drifting floats
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! distribued memory computing library
   USE wrk_nemo        ! working array

   IMPLICIT NONE
   PRIVATE

   PUBLIC   flo_blk    ! routine called by floats.F90

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: floblk.F90 4328 2013-12-06 10:25:13Z davestorkey $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE flo_blk( kt )
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE flo_blk  ***
      !!           
      !! ** Purpose :   Compute the geographical position,latitude, longitude
      !!      and depth of each float at each time step.
      !! 
      !! ** Method  :   The position of a float is computed with Bruno Blanke
      !!      algorithm. We need to know the velocity field, the old positions
      !!      of the floats and the grid defined on the domain.
      !!----------------------------------------------------------------------
      INTEGER, INTENT( in  ) ::   kt ! ocean time step
      !!
      INTEGER :: jfl              ! dummy loop arguments
      INTEGER :: ind, ifin, iloop
      REAL(wp)   ::       &
         zuinfl,zvinfl,zwinfl,      &     ! transport across the input face
         zuoutfl,zvoutfl,zwoutfl,   &     ! transport across the ouput face
         zvol,                      &     ! volume of the mesh
         zsurfz,                    &     ! surface of the face of the mesh 
         zind

      REAL(wp), DIMENSION ( 2 )  ::   zsurfx, zsurfy   ! surface of the face of the mesh

      INTEGER  , POINTER, DIMENSION ( : )  ::   iil, ijl, ikl                   ! index of nearest mesh
      INTEGER  , POINTER, DIMENSION ( : )  ::   iiloc , ijloc              
      INTEGER  , POINTER, DIMENSION ( : )  ::   iiinfl, ijinfl, ikinfl          ! index of input mesh of the float.
      INTEGER  , POINTER, DIMENSION ( : )  ::   iioutfl, ijoutfl, ikoutfl       ! index of output mesh of the float.
      REAL(wp) , POINTER, DIMENSION ( : )  ::   zgifl, zgjfl, zgkfl             ! position of floats, index on 
      !                                                                         ! velocity mesh.
      REAL(wp) , POINTER, DIMENSION ( : )  ::    ztxfl, ztyfl, ztzfl            ! time for a float to quit the mesh
      !                                                                         ! across one of the face x,y and z 
      REAL(wp) , POINTER, DIMENSION ( : )  ::    zttfl                          ! time for a float to quit the mesh 
      REAL(wp) , POINTER, DIMENSION ( : )  ::    zagefl                         ! time during which, trajectorie of 
      !                                                                         ! the float has been computed
      REAL(wp) , POINTER, DIMENSION ( : )  ::   zagenewfl                       ! new age of float after calculation 
      !                                                                         ! of new position
      REAL(wp) , POINTER, DIMENSION ( : )  ::   zufl, zvfl, zwfl                ! interpolated vel. at float position
      REAL(wp) , POINTER, DIMENSION ( : )  ::   zudfl, zvdfl, zwdfl             ! velocity diff input/output of mesh
      REAL(wp) , POINTER, DIMENSION ( : )  ::   zgidfl, zgjdfl, zgkdfl          ! direction index of float
      !!---------------------------------------------------------------------
      CALL wrk_alloc( jpnfl , iil   , ijl   , ikl   , iiloc  ,  ijloc           )
      CALL wrk_alloc( jpnfl , iiinfl, ijinfl, ikinfl, iioutfl, ijoutfl, ikoutfl )
      CALL wrk_alloc( jpnfl , zgifl , zgjfl , zgkfl , ztxfl  , ztyfl  , ztzfl   , zttfl , zagefl, zagenewfl) 
      CALL wrk_alloc( jpnfl , zufl  , zvfl  , zwfl  , zudfl  , zvdfl  , zwdfl   , zgidfl, zgjdfl, zgkdfl   )

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'flo_blk : compute Blanke trajectories for floats '
         IF(lwp) WRITE(numout,*) '~~~~~~~ '
      ENDIF

      ! Initialisation of parameters
      
      DO jfl = 1, jpnfl
         ! ages of floats are put at zero
         zagefl(jfl) = 0.
         ! index on the velocity grid 
         ! We considere k coordinate negative, with this transformation 
         ! the computation in the 3 direction is the same. 
         zgifl(jfl) = tpifl(jfl) - 0.5
         zgjfl(jfl) = tpjfl(jfl) - 0.5
         zgkfl(jfl) = MIN(-1.,-(tpkfl(jfl)))
         ! surface drift every 10 days 
         IF( ln_argo ) THEN
            IF( MOD(kt,150) >= 146 .OR. MOD(kt,150) == 0 )  zgkfl(jfl) = -1.
         ENDIF
         ! index of T mesh
         iil(jfl) = 1 + INT(zgifl(jfl))
         ijl(jfl) = 1 + INT(zgjfl(jfl))
         ikl(jfl) =     INT(zgkfl(jfl))
      END DO
       
      iloop = 0
222   DO jfl = 1, jpnfl
# if   defined key_mpp_mpi
         IF( (iil(jfl) >= (mig(nldi)-jpizoom+1)) .AND. (iil(jfl) <= (mig(nlei)-jpizoom+1)) .AND.   &
             (ijl(jfl) >= (mjg(nldj)-jpjzoom+1)) .AND. (ijl(jfl) <= (mjg(nlej)-jpjzoom+1)) ) THEN
            iiloc(jfl) = iil(jfl) - (mig(1)-jpizoom+1) + 1
            ijloc(jfl) = ijl(jfl) - (mjg(1)-jpjzoom+1) + 1
# else 
            iiloc(jfl) = iil(jfl)
            ijloc(jfl) = ijl(jfl)
# endif
            
            ! compute the transport across the mesh where the float is.            
!!bug (gm) change e3t into fse3. but never checked 
            zsurfx(1) = e2u(iiloc(jfl)-1,ijloc(jfl)  ) * fse3u(iiloc(jfl)-1,ijloc(jfl)  ,-ikl(jfl))
            zsurfx(2) = e2u(iiloc(jfl)  ,ijloc(jfl)  ) * fse3u(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl))
            zsurfy(1) = e1v(iiloc(jfl)  ,ijloc(jfl)-1) * fse3v(iiloc(jfl)  ,ijloc(jfl)-1,-ikl(jfl))
            zsurfy(2) = e1v(iiloc(jfl)  ,ijloc(jfl)  ) * fse3v(iiloc(jfl)  ,ijloc(jfl)  ,-ikl(jfl))

            ! for a isobar float zsurfz is put to zero. The vertical velocity will be zero too.
            zsurfz = e1t(iiloc(jfl),ijloc(jfl)) * e2t(iiloc(jfl),ijloc(jfl))
            zvol   = zsurfz * fse3t(iiloc(jfl),ijloc(jfl),-ikl(jfl))

            !
            zuinfl =( ub(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl)-1,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(1)
            zuoutfl=( ub(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl)) + un(iiloc(jfl)  ,ijloc(jfl),-ikl(jfl)) )/2.*zsurfx(2)
            zvinfl =( vb(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl)-1,-ikl(jfl)) )/2.*zsurfy(1)
            zvoutfl=( vb(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl)) + vn(iiloc(jfl),ijloc(jfl)  ,-ikl(jfl)) )/2.*zsurfy(2)
            zwinfl =-(wb(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1))    &
               &   +  wn(iiloc(jfl),ijloc(jfl),-(ikl(jfl)-1)) )/2. *  zsurfz*nisobfl(jfl)
            zwoutfl=-(wb(iiloc(jfl),ijloc(jfl),- ikl(jfl)   )   &
               &   +  wn(iiloc(jfl),ijloc(jfl),- ikl(jfl)   ) )/2. *  zsurfz*nisobfl(jfl)
            
            ! interpolation of velocity field on the float initial position            
            zufl(jfl)=  zuinfl  + ( zgifl(jfl) - float(iil(jfl)-1) ) * ( zuoutfl - zuinfl)
            zvfl(jfl)=  zvinfl  + ( zgjfl(jfl) - float(ijl(jfl)-1) ) * ( zvoutfl - zvinfl)
            zwfl(jfl)=  zwinfl  + ( zgkfl(jfl) - float(ikl(jfl)-1) ) * ( zwoutfl - zwinfl)
            
            ! faces of input and output
            ! u-direction
            IF( zufl(jfl) < 0. ) THEN
               iioutfl(jfl) = iil(jfl) - 1.
               iiinfl (jfl) = iil(jfl)
               zind   = zuinfl
               zuinfl = zuoutfl
               zuoutfl= zind
            ELSE
               iioutfl(jfl) = iil(jfl)
               iiinfl (jfl) = iil(jfl) - 1
            ENDIF
            ! v-direction       
            IF( zvfl(jfl) < 0. ) THEN
               ijoutfl(jfl) = ijl(jfl) - 1.
               ijinfl (jfl) = ijl(jfl)
               zind    = zvinfl
               zvinfl  = zvoutfl
               zvoutfl = zind
            ELSE
               ijoutfl(jfl) = ijl(jfl)
               ijinfl (jfl) = ijl(jfl) - 1.
            ENDIF
            ! w-direction
            IF( zwfl(jfl) < 0. ) THEN
               ikoutfl(jfl) = ikl(jfl) - 1.
               ikinfl (jfl) = ikl(jfl)
               zind    = zwinfl
               zwinfl  = zwoutfl
               zwoutfl = zind
            ELSE
               ikoutfl(jfl) = ikl(jfl)
               ikinfl (jfl) = ikl(jfl) - 1.
            ENDIF
            
            ! compute the time to go out the mesh across a face
            ! u-direction
            zudfl (jfl) = zuoutfl - zuinfl
            zgidfl(jfl) = float(iioutfl(jfl) - iiinfl(jfl))
            IF( zufl(jfl)*zuoutfl <= 0. ) THEN
               ztxfl(jfl) = 1.E99
            ELSE
               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN
                  ztxfl(jfl)= zgidfl(jfl)/zudfl(jfl) * LOG(zuoutfl/zufl (jfl))
               ELSE
                  ztxfl(jfl)=(float(iioutfl(jfl))-zgifl(jfl))/zufl(jfl)
               ENDIF
               IF( (ABS(zgifl(jfl)-float(iiinfl (jfl))) <=  1.E-7) .OR.   &
                   (ABS(zgifl(jfl)-float(iioutfl(jfl))) <=  1.E-7) ) THEN
                  ztxfl(jfl)=(zgidfl(jfl))/zufl(jfl)
               ENDIF
            ENDIF
            ! v-direction
            zvdfl (jfl) = zvoutfl - zvinfl
            zgjdfl(jfl) = float(ijoutfl(jfl)-ijinfl(jfl))
            IF( zvfl(jfl)*zvoutfl <= 0. ) THEN
               ztyfl(jfl) = 1.E99
            ELSE
               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN
                  ztyfl(jfl) = zgjdfl(jfl)/zvdfl(jfl) * LOG(zvoutfl/zvfl (jfl))
               ELSE
                  ztyfl(jfl) = (float(ijoutfl(jfl)) - zgjfl(jfl))/zvfl(jfl)
               ENDIF
               IF( (ABS(zgjfl(jfl)-float(ijinfl (jfl))) <= 1.E-7) .OR.   &
                   (ABS(zgjfl(jfl)-float(ijoutfl(jfl))) <=  1.E-7) ) THEN
                  ztyfl(jfl) = (zgjdfl(jfl)) / zvfl(jfl)
               ENDIF
            ENDIF
            ! w-direction        
            IF( nisobfl(jfl) == 1. ) THEN 
               zwdfl (jfl) = zwoutfl - zwinfl
               zgkdfl(jfl) = float(ikoutfl(jfl) - ikinfl(jfl))
               IF( zwfl(jfl)*zwoutfl <= 0. ) THEN
                  ztzfl(jfl) = 1.E99
               ELSE
                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN
                     ztzfl(jfl) = zgkdfl(jfl)/zwdfl(jfl) * LOG(zwoutfl/zwfl (jfl))
                  ELSE
                     ztzfl(jfl) = (float(ikoutfl(jfl)) - zgkfl(jfl))/zwfl(jfl)
                  ENDIF
                  IF( (ABS(zgkfl(jfl)-float(ikinfl (jfl))) <=  1.E-7) .OR.   &
                      (ABS(zgkfl(jfl)-float(ikoutfl(jfl))) <= 1.E-7) ) THEN
                     ztzfl(jfl) = (zgkdfl(jfl)) / zwfl(jfl)
                  ENDIF
               ENDIF
            ENDIF
            
            ! the time to go leave the mesh is the smallest time
                   
            IF( nisobfl(jfl) == 1. ) THEN 
               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl),ztzfl(jfl))
            ELSE
               zttfl(jfl) = MIN(ztxfl(jfl),ztyfl(jfl))
            ENDIF
            ! new age of the FLOAT
            zagenewfl(jfl) = zagefl(jfl) + zttfl(jfl)*zvol
            ! test to know if the "age" of the float is not bigger than the 
            ! time step
            IF( zagenewfl(jfl) > rdt ) THEN
               zttfl(jfl) = (rdt-zagefl(jfl)) / zvol
               zagenewfl(jfl) = rdt
            ENDIF
            
            ! In the "minimal" direction we compute the index of new mesh
            ! on i-direction
            IF( ztxfl(jfl) <=  zttfl(jfl) ) THEN
               zgifl(jfl) = float(iioutfl(jfl))
               ind = iioutfl(jfl)
               IF( iioutfl(jfl) >= iiinfl(jfl) ) THEN
                  iioutfl(jfl) = iioutfl(jfl) + 1
               ELSE
                  iioutfl(jfl) = iioutfl(jfl) - 1
               ENDIF
               iiinfl(jfl) = ind
            ELSE
               IF( ABS(zudfl(jfl)) >= 1.E-5 ) THEN 
                  zgifl(jfl) = zgifl(jfl) + zgidfl(jfl)*zufl(jfl)    &
                     &       * ( EXP( zudfl(jfl)/zgidfl(jfl)*zttfl(jfl) ) - 1. ) /  zudfl(jfl)
               ELSE
                  zgifl(jfl) = zgifl(jfl) + zufl(jfl) * zttfl(jfl)
               ENDIF
            ENDIF
            ! on j-direction
            IF( ztyfl(jfl) <= zttfl(jfl) ) THEN
               zgjfl(jfl) = float(ijoutfl(jfl))
               ind = ijoutfl(jfl)
               IF( ijoutfl(jfl) >= ijinfl(jfl) ) THEN
                  ijoutfl(jfl) = ijoutfl(jfl) + 1
               ELSE
                  ijoutfl(jfl) = ijoutfl(jfl) - 1
               ENDIF
               ijinfl(jfl) = ind
            ELSE
               IF( ABS(zvdfl(jfl)) >= 1.E-5 ) THEN 
                  zgjfl(jfl) = zgjfl(jfl)+zgjdfl(jfl)*zvfl(jfl)   &
                     &       * ( EXP(zvdfl(jfl)/zgjdfl(jfl)*zttfl(jfl)) - 1. ) /  zvdfl(jfl)
               ELSE
                  zgjfl(jfl) = zgjfl(jfl)+zvfl(jfl)*zttfl(jfl)
               ENDIF
            ENDIF
            ! on k-direction
            IF( nisobfl(jfl) == 1. ) THEN
               IF( ztzfl(jfl) <= zttfl(jfl) ) THEN
                  zgkfl(jfl) = float(ikoutfl(jfl))
                  ind = ikoutfl(jfl)
                  IF( ikoutfl(jfl) >= ikinfl(jfl) ) THEN
                     ikoutfl(jfl) = ikoutfl(jfl)+1
                  ELSE
                     ikoutfl(jfl) = ikoutfl(jfl)-1
                  ENDIF
                  ikinfl(jfl) = ind
               ELSE
                  IF( ABS(zwdfl(jfl)) >= 1.E-5 ) THEN 
                     zgkfl(jfl) = zgkfl(jfl)+zgkdfl(jfl)*zwfl(jfl)    &
                        &       * ( EXP(zwdfl(jfl)/zgkdfl(jfl)*zttfl(jfl)) - 1. ) /  zwdfl(jfl)
                  ELSE
                     zgkfl(jfl) = zgkfl(jfl)+zwfl(jfl)*zttfl(jfl)
                  ENDIF
               ENDIF
            ENDIF
            
            ! coordinate of the new point on the temperature grid
            
            iil(jfl) = MAX(iiinfl(jfl),iioutfl(jfl))
            ijl(jfl) = MAX(ijinfl(jfl),ijoutfl(jfl))
            IF( nisobfl(jfl) ==  1 ) ikl(jfl) = MAX(ikinfl(jfl),ikoutfl(jfl))
!!Alexcadm	 write(*,*)'PE ',narea,
!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
!!Alexcadm     .		,ijoutfl(jfl),ikinfl(jfl),
!!Alexcadm     .	  ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
!!Alexcadm     .		,ztzfl(jfl),zgifl(jfl),
!!Alexcadm     .  zgjfl(jfl) 
!!Alexcadm	IF (jfl == 910) write(*,*)'Flotteur 910',
!!Alexcadm     .    iiinfl(jfl),iioutfl(jfl),ijinfl(jfl)
!!Alexcadm     .		,ijoutfl(jfl),ikinfl(jfl),
!!Alexcadm     .	  ikoutfl(jfl),ztxfl(jfl),ztyfl(jfl)
!!Alexcadm     .		,ztzfl(jfl),zgifl(jfl),
!!Alexcadm     .  zgjfl(jfl) 
            ! reinitialisation of the age of FLOAT
            zagefl(jfl) = zagenewfl(jfl)
# if   defined key_mpp_mpi
         ELSE
            ! we put zgifl, zgjfl, zgkfl, zagefl
            zgifl (jfl) = 0.
            zgjfl (jfl) = 0.
            zgkfl (jfl) = 0.
            zagefl(jfl) = 0.
            iil(jfl) = 0
            ijl(jfl) = 0
         ENDIF
# endif
      END DO
      
      ! synchronisation
      IF( lk_mpp )   CALL mpp_sum( zgifl , jpnfl )   ! sums over the global domain
      IF( lk_mpp )   CALL mpp_sum( zgjfl , jpnfl )
      IF( lk_mpp )   CALL mpp_sum( zgkfl , jpnfl )
      IF( lk_mpp )   CALL mpp_sum( zagefl, jpnfl )
      IF( lk_mpp )   CALL mpp_sum( iil   , jpnfl )
      IF( lk_mpp )   CALL mpp_sum( ijl   , jpnfl )
      
      ! Test to know if a  float hasn't integrated enought time
      IF( ln_argo ) THEN
         ifin = 1
         DO jfl = 1, jpnfl
            IF( zagefl(jfl) < rdt )   ifin = 0
            tpifl(jfl) = zgifl(jfl) + 0.5
            tpjfl(jfl) = zgjfl(jfl) + 0.5
         END DO
      ELSE
         ifin = 1
         DO jfl = 1, jpnfl
            IF( zagefl(jfl) < rdt )   ifin = 0
            tpifl(jfl) = zgifl(jfl) + 0.5
            tpjfl(jfl) = zgjfl(jfl) + 0.5
            IF( nisobfl(jfl) == 1 ) tpkfl(jfl) = -(zgkfl(jfl))
         END DO
      ENDIF
!!Alexcadm	IF (lwp) write(numout,*) '---------'
!!Alexcadm	IF (lwp) write(numout,*) 'before Erika:',tpifl(880),tpjfl(880),
!!Alexcadm     .       tpkfl(880),zufl(880),zvfl(880),zwfl(880)
!!Alexcadm	IF (lwp) write(numout,*) 'first Erika:',tpifl(900),tpjfl(900),
!!Alexcadm     .       tpkfl(900),zufl(900),zvfl(900),zwfl(900)
!!Alexcadm	IF (lwp) write(numout,*) 'last Erika:',tpifl(jpnfl),tpjfl(jpnfl),
!!Alexcadm     .       tpkfl(jpnfl),zufl(jpnfl),zvfl(jpnfl),zwfl(jpnfl)
      IF( ifin == 0 ) THEN
         iloop = iloop + 1 
         GO TO 222
      ENDIF
      !
      CALL wrk_dealloc( jpnfl , iil   , ijl   , ikl   , iiloc  ,  ijloc           )
      CALL wrk_dealloc( jpnfl , iiinfl, ijinfl, ikinfl, iioutfl, ijoutfl, ikoutfl )
      CALL wrk_dealloc( jpnfl , zgifl , zgjfl , zgkfl , ztxfl  , ztyfl  , ztzfl   , zttfl , zagefl, zagenewfl) 
      CALL wrk_dealloc( jpnfl , zufl  , zvfl  , zwfl  , zudfl  , zvdfl  , zwdfl   , zgidfl, zgjdfl, zgkdfl   )
      !
   END SUBROUTINE flo_blk

#  else
   !!----------------------------------------------------------------------
   !!   Default option                                         Empty module
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE flo_blk                  ! Empty routine
   END SUBROUTINE flo_blk 
#endif
   
   !!======================================================================
END MODULE floblk