ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/BDY/bdydyn2d.F90

MODULE bdydyn2d
   !!======================================================================
   !!                       ***  MODULE  bdydyn  ***
   !! Unstructured Open Boundary Cond. :   Apply boundary conditions to barotropic solution
   !!======================================================================
   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
   !!----------------------------------------------------------------------
#if defined key_bdy 
   !!----------------------------------------------------------------------
   !!   'key_bdy' :                    Unstructured Open Boundary Condition
   !!----------------------------------------------------------------------
   !!   bdy_dyn2d          : Apply open boundary conditions to barotropic variables.
   !!   bdy_dyn2d_frs      : Apply Flow Relaxation Scheme 
   !!   bdy_dyn2d_fla      : Apply Flather condition
   !!   bdy_dyn2d_orlanski : Orlanski Radiation
   !!   bdy_ssh            : Duplicate sea level across open boundaries
   !!----------------------------------------------------------------------
   USE timing          ! Timing
   USE oce             ! ocean dynamics and tracers 
   USE dom_oce         ! ocean space and time domain
   USE bdy_oce         ! ocean open boundary conditions
   USE bdylib          ! BDY library routines
   USE dynspg_oce      ! for barotropic variables
   USE phycst          ! physical constants
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE in_out_manager  !

   IMPLICIT NONE
   PRIVATE

   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn
   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: bdydyn2d.F90 2678 2015-11-26 09:59:07Z ufla $ 
   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE bdy_dyn2d( kt, pua2d, pva2d, pub2d, pvb2d, phur, phvr, pssh  )
      !!----------------------------------------------------------------------
      !!                  ***  SUBROUTINE bdy_dyn2d  ***
      !!
      !! ** Purpose : - Apply open boundary conditions for barotropic variables
      !!
      !!----------------------------------------------------------------------
      INTEGER,                      INTENT(in) ::   kt   ! Main time step counter
      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pua2d, pva2d 
      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pub2d, pvb2d
      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: phur, phvr
      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pssh
      !!
      INTEGER                                  ::   ib_bdy ! Loop counter

      DO ib_bdy=1, nb_bdy

         SELECT CASE( cn_dyn2d(ib_bdy) )
         CASE('none')
            CYCLE
         CASE('frs')
            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d )
         CASE('flather')
            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr )
         CASE('orlanski')
            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.false.)
         CASE('orlanski_npo')
            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.true. )
         CASE DEFAULT
            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
         END SELECT
      ENDDO

   END SUBROUTINE bdy_dyn2d

   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
      !!----------------------------------------------------------------------
      !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
      !!
      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
      !!                at open boundaries.
      !!
      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in 
      !!               a three-dimensional baroclinic ocean model with realistic
      !!               topography. Tellus, 365-382.
      !!----------------------------------------------------------------------
      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pua2d, pva2d 
      !!
      INTEGER  ::   jb, jk         ! dummy loop indices
      INTEGER  ::   ii, ij, igrd   ! local integers
      REAL(wp) ::   zwgt           ! boundary weight
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_frs')
      !
      igrd = 2                      ! Relaxation of zonal velocity
      DO jb = 1, idx%nblen(igrd)
         ii   = idx%nbi(jb,igrd)
         ij   = idx%nbj(jb,igrd)
         zwgt = idx%nbw(jb,igrd)
         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
      END DO
      !
      igrd = 3                      ! Relaxation of meridional velocity
      DO jb = 1, idx%nblen(igrd)
         ii   = idx%nbi(jb,igrd)
         ij   = idx%nbj(jb,igrd)
         zwgt = idx%nbw(jb,igrd)
         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
      END DO 
      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) 
      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
      !
      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
      !

   END SUBROUTINE bdy_dyn2d_frs


   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )
      !!----------------------------------------------------------------------
      !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
      !!             
      !!              - Apply Flather boundary conditions on normal barotropic velocities 
      !!
      !! ** WARNINGS about FLATHER implementation:
      !!1. According to Palma and Matano, 1998 "after ssh" is used. 
      !!   In ROMS and POM implementations, it is "now ssh". In the current 
      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable. 
      !!   So I use "before ssh" in the following.
      !!
      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of 
      !!   fact, the model ssh just inside the dynamical boundary is used (the outside  
      !!   ssh in the code is not updated).
      !!
      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
      !!----------------------------------------------------------------------
      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pua2d, pva2d
      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pssh, phur, phvr 

      INTEGER  ::   jb, igrd                         ! dummy loop indices
      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
      REAL(wp), POINTER :: flagu, flagv              ! short cuts
      REAL(wp) ::   zcorr                            ! Flather correction
      REAL(wp) ::   zforc                            ! temporary scalar
      REAL(wp) ::   zflag, z1_2                      !    "        "
      !!----------------------------------------------------------------------

      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')

      z1_2 = 0.5_wp

      ! ---------------------------------!
      ! Flather boundary conditions     :!
      ! ---------------------------------! 
     
!!! REPLACE spgu with nemo_wrk work space

      ! Fill temporary array with ssh data (here spgu):
      igrd = 1
      spgu(:,:) = 0.0
      DO jb = 1, idx%nblenrim(igrd)
         ii = idx%nbi(jb,igrd)
         ij = idx%nbj(jb,igrd)
         spgu(ii, ij) = dta%ssh(jb)
      END DO

      CALL lbc_bdy_lnk( spgu(:,:), 'T', 1., ib_bdy )
      !
      igrd = 2      ! Flather bc on u-velocity; 
      !             ! remember that flagu=-1 if normal velocity direction is outward
      !             ! I think we should rather use after ssh ?
      DO jb = 1, idx%nblenrim(igrd)
         ii  = idx%nbi(jb,igrd)
         ij  = idx%nbj(jb,igrd) 
         flagu => idx%flagu(jb,igrd)
         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary 
         !
         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )

         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
         ! Use characteristics method instead
         zflag = ABS(flagu)
         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
      END DO
      !
      igrd = 3      ! Flather bc on v-velocity
      !             ! remember that flagv=-1 if normal velocity direction is outward
      DO jb = 1, idx%nblenrim(igrd)
         ii  = idx%nbi(jb,igrd)
         ij  = idx%nbj(jb,igrd) 
         flagv => idx%flagv(jb,igrd)
         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary 
         !
         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )

         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
         ! Use characteristics method instead
         zflag = ABS(flagv)
         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
      END DO
      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
      !
      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
      !
   END SUBROUTINE bdy_dyn2d_fla


   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
      !!----------------------------------------------------------------------
      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
      !!             
      !!              - Apply Orlanski radiation condition adaptively:
      !!                  - radiation plus weak nudging at outflow points
      !!                  - no radiation and strong nudging at inflow points
      !! 
      !!
      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
      !!----------------------------------------------------------------------
      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
      REAL(wp), DIMENSION(:,:),     INTENT(inout) :: pua2d, pva2d
      REAL(wp), DIMENSION(:,:),     INTENT(in) :: pub2d, pvb2d 
      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version

      INTEGER  ::   ib, igrd                               ! dummy loop indices
      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
      !!----------------------------------------------------------------------

      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
      !
      igrd = 2      ! Orlanski bc on u-velocity; 
      !            
      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )

      igrd = 3      ! Orlanski bc on v-velocity
      !  
      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
      !
      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
      !
      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
      !
      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
      !
   END SUBROUTINE bdy_dyn2d_orlanski

   SUBROUTINE bdy_ssh( zssh )
      !!----------------------------------------------------------------------
      !!                  ***  SUBROUTINE bdy_ssh  ***
      !!
      !! ** Purpose : Duplicate sea level across open boundaries
      !!
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zssh ! Sea level
      !!
      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
      INTEGER  ::   ii, ij, zcoef, ip, jp   !   "       "

      igrd = 1                       ! Everything is at T-points here

      DO ib_bdy = 1, nb_bdy
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
            ! Set gradient direction:
            zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
            IF ( zcoef == 0 ) THEN
               zssh(ii,ij) = 0._wp
            ELSE
               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
            ENDIF
         END DO

         ! Boundary points should be updated
         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
      END DO

   END SUBROUTINE bdy_ssh

#else
   !!----------------------------------------------------------------------
   !!   Dummy module                   NO Unstruct Open Boundary Conditions
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
      INTEGER, intent(in) :: kt
      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt
   END SUBROUTINE bdy_dyn2d

#endif

   !!======================================================================
END MODULE bdydyn2d