rp_omega.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE rp_omega_mod
#if defined TL_IOMS && defined SOLVE3D
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine computes S-coordinate vertical velocity (m^3/s),       !
!                                                                      !
!                  W=[Hz/(m*n)]*omega,                                 !
!                                                                      !
!  diagnostically at horizontal RHO-points and vertical W-points.      !
!                                                                      !
!  BASIC STATE variables needed: W, z_w.                               !
!                                                                      !
!  NOTE: We need to recompute  basic state W in this routine since     !
!  ----  intermediate values of W are needed by the tangent linear     !
!        and adjoint routines.                                         !
!                                                                      !
!=======================================================================
!
      implicit none
!
      PRIVATE
      PUBLIC  :: rp_omega
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE rp_omega (ng, tile, model)
!***********************************************************************
!
      USE mod_param
      USE mod_grid
      USE mod_ocean
# if defined SEDIMENT && defined SED_MORPH
      USE mod_sedbed
      USE mod_stepping
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, model, 13, __line__, myfile)
# endif
      CALL rp_omega_tile (ng, tile, model,                              &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    imins, imaxs, jmins, jmaxs,                   &
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
     &                    nstp(ng), nnew(ng),                           &
     &                    grid(ng) % omn,                               &
     &                    sedbed(ng) % bed_thick,                       &
     &                    sedbed(ng) % tl_bed_thick,                    &
# endif
     &                    grid(ng) % Huon,                              &
     &                    grid(ng) % Hvom,                              &
     &                    grid(ng) % z_w,                               &
     &                    grid(ng) % tl_Huon,                           &
     &                    grid(ng) % tl_Hvom,                           &
     &                    grid(ng) % tl_z_w,                            &
     &                    ocean(ng) % W,                                &
     &                    ocean(ng) % tl_W)
# ifdef PROFILE
      CALL wclock_off (ng, model, 13, __line__, myfile)
# endif
!
      RETURN
      END SUBROUTINE rp_omega
!
!***********************************************************************
      SUBROUTINE rp_omega_tile (ng, tile, model,                        &
     &                          LBi, UBi, LBj, UBj,                     &
     &                          IminS, ImaxS, JminS, JmaxS,             &
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
     &                          nstp, nnew,                             &
     &                          omn,
     &                          bed_thick, tl_bed_thick,                &
# endif
     &                          Huon, Hvom, z_w,                        &
     &                          tl_Huon, tl_Hvom, tl_z_w,               &
     &                          W, tl_W)
!***********************************************************************
!
      USE mod_param
      USE mod_scalars
      USE mod_sources
!
      USE bc_3d_mod, ONLY : bc_w3d_tile
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange3d
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
      integer, intent(in) :: nstp, nnew
# endif
!
# ifdef ASSUMED_SHAPE
#  if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
      real(r8), intent(in) :: omn(LBi:,LBj:)
      real(r8), intent(in):: bed_thick(LBi:,LBj:,:)
      real(r8), intent(in):: tl_bed_thick(LBi:,LBj:,:)
#  endif
      real(r8), intent(in) :: Huon(LBi:,LBj:,:)
      real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
      real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
      real(r8), intent(in) :: tl_Huon(LBi:,LBj:,:)
      real(r8), intent(in) :: tl_Hvom(LBi:,LBj:,:)
      real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
 
      real(r8), intent(out) :: W(LBi:,LBj:,0:)
      real(r8), intent(out) :: tl_W(LBi:,LBj:,0:)
 
# else
 
#  if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
      real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
      real(r8), intent(in):: bed_thick(LBi:UBi,LBj:UBj,3)
      real(r8), intent(in):: tl_bed_thick(LBi:UBi,LBj:UBj,3)
#  endif
      real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(in) :: tl_Huon(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: tl_Hvom(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
 
      real(r8), intent(out) :: W(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(out) :: tl_W(LBi:UBi,LBj:UBj,0:N(ng))
# endif
!
!  Local variable declarations.
!
      integer :: i, ii, is, j, jj, k
      real(r8) :: cff, tl_cff
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
      real(r8) :: cff1
# endif
      real(r8), dimension(IminS:ImaxS) :: wrk
      real(r8), dimension(IminS:ImaxS) :: tl_wrk
 
# include "set_bounds.h"
!
!------------------------------------------------------------------------
!  Vertically integrate horizontal mass flux divergence.
!------------------------------------------------------------------------
!
!  Starting with zero vertical velocity at the bottom, integrate
!  from the bottom (k=0) to the free-surface (k=N).  The w(:,:,N(ng))
!  contains the vertical velocity at the free-surface, d(zeta)/d(t).
!  Notice that barotropic mass flux divergence is not used directly.
!
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
      cff1=1.0_r8/dt(ng)
# endif
      DO j=jstr,jend
        DO i=istr,iend
# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
          w(i,j,0)=-cff1*(bed_thick(i,j,nstp)-                          &
     &                    bed_thick(i,j,nnew))*omn(i,j)
          tl_w(i,j,0)=-cff1*(tl_bed_thick(i,j,nstp)-                    &
     &                       tl_bed_thick(i,j,nnew))*omn(i,j)
# else
          w(i,j,0)=0.0_r8
          tl_w(i,j,0)=0.0_r8
# endif
        END DO
!
!  Code added to clear tl_W to be consistent with adjoint.
!
        DO k=1,n(ng)
          DO i=istr,iend
           tl_w(i,j,k)=0.0_r8
          END DO
        END DO
        DO k=1,n(ng)
          DO i=istr,iend
            w(i,j,k)=w(i,j,k-1)-                                        &
     &               (huon(i+1,j,k)-huon(i,j,k)+                        &
     &                hvom(i,j+1,k)-hvom(i,j,k))
            tl_w(i,j,k)=tl_w(i,j,k-1)-                                  &
     &                  (tl_huon(i+1,j,k)-tl_huon(i,j,k)+               &
     &                   tl_hvom(i,j+1,k)-tl_hvom(i,j,k))
          END DO
        END DO
!
!  Apply mass point sources (volume vertical influx), if any.
!
!  Overwrite W(Isrc,Jsrc,k) with the same divergence of Huon,Hvom as
!  above but add in point source Qsrc(k) and reaccumulate the vertical
!  sum to obtain the correct net Qbar given in user input - J. Levin
!  (Jupiter Intelligence Inc.) and J. Wilkin
!
!    Dsrc(is) = 2,  flow across grid cell w-face (positive or negative)
!
        IF (lwsrc(ng)) THEN
          DO is=1,nsrc(ng)
            IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
              ii=sources(ng)%Isrc(is)
              jj=sources(ng)%Jsrc(is)
              IF (((istrr.le.ii).and.(ii.le.iendr)).and.                &
     &            ((jstrr.le.jj).and.(jj.le.jendr)).and.                &
     &            (j.eq.jj)) THEN
                DO k=1,n(ng)
                  w(ii,jj,k)=w(ii,jj,k-1)-                              &
     &                       (huon(ii+1,jj,k)-huon(ii,jj,k)+            &
     &                        hvom(ii,jj+1,k)-hvom(ii,jj,k))+           &
     &                       sources(ng)%Qsrc(is,k)
                  tl_w(ii,jj,k)=tl_w(ii,jj,k-1)-                        &
     &                          (tl_huon(ii+1,jj,k)-tl_huon(ii,jj,k)+   &
     &                           tl_hvom(ii,jj+1,k)-tl_hvom(ii,jj,k))+  &
     &                          sources(ng)%tl_Qsrc(is,k)
                END DO
              END IF
            END IF
          END DO
        END IF
!
        DO i=istr,iend
          cff=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
          tl_cff=-cff*cff*(tl_z_w(i,j,n(ng))-tl_z_w(i,j,0))+            &
# ifdef TL_IOMS
     &           2.0_r8*cff
# endif
          wrk(i)=cff*w(i,j,n(ng))
          tl_wrk(i)=tl_cff*w(i,j,n(ng))+cff*tl_w(i,j,n(ng))-            &
# ifdef TL_IOMS
     &              wrk(i)
# endif
        END DO
!
!  In order to insure zero vertical velocity at the free-surface,
!  subtract the vertical velocities of the moving S-coordinates
!  isosurfaces. These isosurfaces are proportional to d(zeta)/d(t).
!  The proportionally coefficients are a linear function of the
!  S-coordinate with zero value at the bottom (k=0) and unity at
!  the free-surface (k=N).
!
        DO k=n(ng)-1,1,-1
          DO i=istr,iend
            w(i,j,k)=w(i,j,k)-                                          &
     &               wrk(i)*(z_w(i,j,k)-z_w(i,j,0))
            tl_w(i,j,k)=tl_w(i,j,k)-                                    &
     &                  tl_wrk(i)*(z_w(i,j,k)-z_w(i,j,0))-              &
     &                  wrk(i)*(tl_z_w(i,j,k)-tl_z_w(i,j,0))+           &
# ifdef TL_IOMS
     &                  wrk(i)*(z_w(i,j,k)-z_w(i,j,0))
# endif
          END DO
        END DO
        DO i=istr,iend
          w(i,j,n(ng))=0.0_r8
          tl_w(i,j,n(ng))=0.0_r8
        END DO
      END DO
!
!  Set lateral boundary conditions.
!
      CALL bc_w3d_tile (ng, tile,                                       &
     &                  lbi, ubi, lbj, ubj, 0, n(ng),                   &
     &                  w)
      CALL bc_w3d_tile (ng, tile,                                       &
     &                  lbi, ubi, lbj, ubj, 0, n(ng),                   &
     &                  tl_w)
 
# ifdef DISTRIBUTE
      CALL mp_exchange3d (ng, tile, model, 1,                           &
     &                    lbi, ubi, lbj, ubj, 0, n(ng),                 &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    w)
      CALL mp_exchange3d (ng, tile, model, 1,                           &
     &                    lbi, ubi, lbj, ubj, 0, n(ng),                 &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    tl_w)
# endif
!
      RETURN
      END SUBROUTINE rp_omega_tile
#endif
      END MODULE rp_omega_mod