ana_sponge.h

Go to the documentation of this file.

!!
      SUBROUTINE ana_sponge (ng, tile, model)
!
!! git $Id$
!!================================================= Hernan G. Arango ===
!! Copyright (c) 2002-2025 The ROMS Group                              !
!!   Licensed under a MIT/X style license                              !
!!   See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This subroutine computes spatially varying horizontal mixing        !
!  coefficients in sponge areas. The horizontal viscosity and/or       !
!  diffusion are increased in sponge areas to supress noise due        !
!  open boundary conditions or nesting.                                !
!                                                                      !
!  There are two different ways to code the increased viscosity        !
!  and/or diffusion in the sponge areas:                               !
!                                                                      !
!  (1) Rescale the horizontal mixing coefficients computed earlier     !
!      in routine "ini_hmixcoef.F" with a nondimentional factor:       !
!                                                                      !
!      visc2_r(i,j) = ABS(factor(i,j)) * visc2_r(i,j)                  !
!                                                                      !
!      visc2_p(i,j) = 0.25_r8 * ABS(factor(i-1,j-1)+                   !
!                                   factor(i  ,j-1)+                   !
!                                   factor(i-1,j  )+                   !
!                                   factor(i  ,j  )) * visc2_p(i,j)    !
!                                                                      !
!      where factor(i,j) is defined at RHO-points and its values       !
!      can be ZERO (no mixing), ONE (same values), or linearly         !
!      greater than ONE (sponge are with larger mixing).               !
!                                                                      !
!      (See Adriatic Sea application below)                            !
!                                                                      !
!  (2) Overwrite the horizontal mixing coefficients computed earlier   !
!      in routine "ini_hmixcoef.F" with a new distribution:            !
!                                                                      !
!      visc2_r(i,j) = my_values(i,j)                                   !
!                                                                      !
!      visc2_p(i,j) = my_values(i,j)                                   !
!                                                                      !
!      (See US West Coast application below)                           !
!                                                                      !
!  However,                                                            !
!                                                                      !
!  It is HIGHLY recommended to write the nondimentional spatial        !
!  distribution arrays "visc_factor(i,j)" and "diff_factor(i,j)"       !
!  into the grid NetCDF file instead of using the analytical code      !
!  below. It is very easy to introduce parallel bugs.  Also, Users     !
!  can plot their spatial distribution and fine tune their values      !
!  during at the pre-proccessing stage for a particular application.   !
!                                                                      !
!  The Metadata for these scale factors in the Grid NetCDF file is     !
!  as follows (spherical grid case):                                   !
!                                                                      !
!    double visc_factor(eta_rho, xi_rho)                               !
!        visc_factor:long_name = "horizontal viscosity factor"         !
!        visc_factor:valid_min = 0.                                    !
!        visc_factor:coordinates = "lon_rho lat_rho"                   !
!                                                                      !
!    double diff_factor(eta_rho, xi_rho)                               !
!        diff_factor:long_name = "horizontal diffusivity factor"       !
!        diff_factor:valid_min = 0.                                    !
!        diff_factor:coordinates = "lon_rho lat_rho"                   !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_ncparam
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
#include "tile.h"
 
      CALL ana_sponge_tile (ng, tile, model,                            &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      imins, imaxs, jmins, jmaxs)
!
! Set analytical header file name used.
!
#ifdef DISTRIBUTE
      IF (lanafile) THEN
#else
      IF (lanafile.and.(tile.eq.0)) THEN
#endif
        ananame( 8)=myfile
      END IF
!
      RETURN
      END SUBROUTINE ana_sponge
!
!***********************************************************************
      SUBROUTINE ana_sponge_tile (ng, tile, model,                      &
     &                            LBi, UBi, LBj, UBj,                   &
     &                            IminS, ImaxS, JminS, JmaxS)
!***********************************************************************
!
      USE mod_param
      USE mod_grid
      USE mod_mixing
      USE mod_scalars
!
      USE exchange_2d_mod
#ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d
# ifdef SOLVE3D
      USE mp_exchange_mod, ONLY : mp_exchange3d
# endif
#endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
!
!  Local variable declarations.
!
      integer :: i, itrc, j
!
      real(r8) :: cff, innerF, outerF, val, width
!
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: factor
 
#include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Increase horizontal mixing in the sponge areas.
!-----------------------------------------------------------------------
 
#if defined ADRIA02
!
!  Adriatic Sea southern sponge areas.
!
      cff=4.0_r8                          ! quadruple horizontal mixing
      width=6.0_r8
!
!  Set horizontal mixing factor over 6 grid points in the southern
!  boundary and increase linearly its values by four. Otherwise, the
!  horizontal mixing is set to zero.
!
      factor(imins:imaxs,jmins:jmaxs)=0.0_r8               ! initialize
 
      DO i=istrt,iendt
        DO j=jstrt,min(int(width),jendt)
          factor(i,j)=1.0_r8+(cff-1.0_r8)*(width-real(j,r8))/width
        END DO
      END DO
 
# if defined UV_VIS2
      IF (luvsponge(ng)) THEN
        DO i=istrt,iendt
          DO j=jstrt,jendt
            mixing(ng) % visc2_r(i,j)=abs(factor(i,j))*                 &
     &                                mixing(ng) % visc2_r(i,j)
          END DO
        END DO
        DO j=jstrp,jendt
          DO i=istrp,iendt
            mixing(ng) % visc2_p(i,j)=0.25_r8*abs(factor(i-1,j-1)+      &
     &                                            factor(i  ,j-1)+      &
     &                                            factor(i-1,j  )+      &
     &                                            factor(i  ,j  ))*     &
     &                                mixing(ng) % visc2_p(i,j)
          END DO
        END DO
      END IF
# endif
 
# if defined TS_DIF2
      DO itrc=1,nt(ng)
        IF (ltracersponge(itrc,ng)) THEN
          DO i=istrt,iendt
            DO j=jstrt,jendt
              mixing(ng) % diff2(i,j,itrc)=abs(factor(i,j)*             &
     &                                     mixing(ng) % diff2(i,j,itrc)
            END DO
          END DO
        END IF
      END DO
# endif
 
#elif defined WC13
!
!  US West Coast sponge areas.
!
      width=user(1)                        ! sponge width in grid points
 
# if defined UV_VIS2
!
!  Momentum sponge regions:  sponge viscosities as in Marchesiello
!  et al 2003.
!
      IF (luvsponge(ng)) THEN
        innerf=visc2(ng)                   ! inner limit match value
        outerf=100.0_r8                    ! outer limit maximum value
!
!  Southern edge.
!
        DO j=jstrt,min(int(width),jendt)
          val=innerf+(outerf-innerf)*(width-real(j,r8))/width
          DO i=istrt,iendt
            mixing(ng)%visc2_r(i,j)=max(min(val,outerf),innerf)
            mixing(ng)%visc2_p(i,j)=max(min(val,outerf),innerf)
          END DO
        END DO
!
!  Northern edge.
!
        DO j=max(jstrt,mm(ng)+1-int(width)),jendt
          val=outerf+(innerf-outerf)*real(mm(ng)+1-j,r8)/width
          DO i=istrt,iendt
            mixing(ng) % visc2_r(i,j)=max(min(val,outerf),innerf)
            mixing(ng) % visc2_p(i,j)=max(min(val,outerf),innerf)
          END DO
        END DO
!
!  Western edge.
!
        DO i=istrt,min(int(width),iendt)
          DO j=max(jstrt,i),min(mm(ng)+1-i,jendt)
            val=innerf+(outerf-innerf)*(width-real(i,r8))/width
            mixing(ng) % visc2_r(i,j)=max(min(val,outerf),innerf)
            mixing(ng) % visc2_p(i,j)=max(min(val,outerf),innerf)
          END DO
        END DO
      END IF
# endif
 
# if defined TS_DIF2
!
!  Tracer sponge regions: sponge diffusivities as in Marchesiello
!  et al 2003.
!
      DO itrc=1,nt(ng)
        IF (ltracersponge(itrc,ng)) THEN
          innerf=tnu2(itrc,ng)             ! inner limit match value
          outerf=50.0_r8                   ! outer limit maximum value
!
!  Southern edge.
!
          DO j=jstrt,min(int(width),jendt)
            val=innerf+(outerf-innerf)*(width-real(j,r8))/width
            DO i=istrt,iendt
              mixing(ng) % diff2(i,j,itrc)=max(min(val,outerf),innerf)
            END DO
          END DO
!
!  Northern edge.
!
          DO j=max(jstrt,mm(ng)+1-int(width)),jendt
            val=outerf+(innerf-outerf)*real(mm(ng)+1-j,r8)/width
            DO i=istrt,iendt
              mixing(ng) % diff2(i,j,itrc)=max(min(val,outerf),innerf)
            END DO
          END DO
!
!  Western edge.
!
          DO i=istrt,min(int(width),iendt)
            DO j=max(jstrt,i),min(mm(ng)+1-i,jendt)
              val=innerf+(outerf-innerf)*(width-real(i,r8))/width
              mixing(ng) % diff2(i,j,itrc)=max(min(val,outerf),innerf)
            END DO
          END DO
        END IF
      END DO
# endif
#endif
!
!-----------------------------------------------------------------------
!  Exchange boundary data.
!-----------------------------------------------------------------------
!
!! WARNING:  This section is generic for all applications. Please do not
!!           change the code below.
!!
#ifdef UV_VIS2
      IF (luvsponge(ng)) THEN
        IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
          CALL exchange_r2d_tile (ng, tile,                             &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            mixing(ng) % visc2_r)
          CALL exchange_p2d_tile (ng, tile,                             &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            mixing(ng) % visc2_p)
        END IF
      END IF
#endif
 
#ifdef UV_VIS4
      IF (luvsponge(ng)) THEN
        IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
          CALL exchange_r2d_tile (ng, tile,                             &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            mixing(ng) % visc4_r)
          CALL exchange_p2d_tile (ng, tile,                             &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            mixing(ng) % visc4_p)
        END IF
      END IF
#endif
 
#ifdef SOLVE3D
# ifdef TS_DIF2
      IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
        DO itrc=1,nt(ng)
          IF (ltracersponge(itrc,ng)) THEN
            CALL exchange_r2d_tile (ng, tile,                           &
     &                              lbi, ubi, lbj, ubj,                 &
     &                              mixing(ng) % diff2(:,:,itrc))
          END IF
        END DO
      END IF
# endif
 
# ifdef TS_DIF4
      IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
        DO itrc=1,nt(ng)
          IF (ltracersponge(itrc,ng)) THEN
            CALL exchange_r2d_tile (ng, tile,                           &
     &                              lbi, ubi, lbj, ubj,                 &
     &                              mixing(ng) % diff4(:,:,itrc))
          END IF
        END DO
      END IF
# endif
#endif
 
#ifdef DISTRIBUTE
!
# ifdef UV_VIS2
      IF (luvsponge(ng)) THEN
        CALL mp_exchange2d (ng, tile, model, 2,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      mixing(ng) % visc2_r,                       &
     &                      mixing(ng) % visc2_p)
      END IF
# endif
 
# ifdef UV_VIS4
      IF (luvsponge(ng)) THEN
        CALL mp_exchange2d (ng, tile, model, 2,                         &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      mixing(ng) % visc4_r,                       &
     &                      mixing(ng) % visc4_p)
      END IF
# endif
 
# ifdef SOLVE3D
#  ifdef TS_DIF2
      IF (any(ltracersponge(:,ng))) THEN
        CALL mp_exchange3d (ng, tile, model, 1,                         &
     &                      lbi, ubi, lbj, ubj, 1, nt(ng),              &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      mixing(ng) % diff2)
      END IF
#  endif
 
#  ifdef TS_DIF4
      IF (any(ltracersponge(:,ng))) THEN
        CALL mp_exchange3d (ng, tile, model, 1,                         &
     &                      lbi, ubi, lbj, ubj, 1, nt(ng),              &
     &                      nghostpoints,                               &
     &                      ewperiodic(ng), nsperiodic(ng),             &
     &                      mixing(ng) % diff4)
      END IF
#  endif
# endif
#endif
!
      RETURN
      END SUBROUTINE ana_sponge_tile