MODULE prsgrd_mod
!
!git $Id$
!svn $Id: prsgrd31.h 1178 2023-07-11 17:50:57Z arango $
!=======================================================================
! Copyright (c) 2002-2023 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.md Hernan G. Arango !
!========================================== Alexander F. Shchepetkin ===
! !
! This subroutine evaluates the baroclinic hydrostatic pressure !
! gradient term using the STANDARD density Jacobian or WEIGHTED !
! density Jacobian scheme of Song (1998). Both of these approaches !
! compute horizontal differences of density before of the vertical !
! integration. !
! !
! The pressure gradient terms (m4/s2) are loaded into right-hand- !
! side arrays "ru" and "rv". !
! !
! Reference: !
! !
! Song, Y.T., 1998: A general pressure gradient formulation for !
! numerical ocean models. Part I: Scheme design and diagnostic !
! analysis, Monthly Weather Rev., 126, 3213-3230. !
! !
!=======================================================================
!
implicit none
!
PRIVATE
PUBLIC :: prsgrd
!
CONTAINS
!
!***********************************************************************
SUBROUTINE prsgrd (ng, tile)
!***********************************************************************
!
USE mod_param
#ifdef DIAGNOSTICS
USE mod_diags
#endif
#ifdef ATM_PRESS
USE mod_forces
#endif
USE mod_grid
USE mod_ocean
USE mod_stepping
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
!
! Local variable declarations.
!
character (len=*), parameter :: MyFile = &
& __FILE__
!
#include "tile.h"
!
#ifdef PROFILE
CALL wclock_on (ng, iNLM, 23, __LINE__, MyFile)
#endif
CALL prsgrd31_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs(ng), &
#ifdef WET_DRY
& GRID(ng)%umask_wet, &
& GRID(ng)%vmask_wet, &
#endif
& GRID(ng) % Hz, &
& GRID(ng) % om_v, &
& GRID(ng) % on_u, &
& GRID(ng) % z_r, &
& GRID(ng) % z_w, &
& OCEAN(ng) % rho, &
#ifdef TIDE_GENERATING_FORCES
& OCEAN(ng) % eq_tide, &
#endif
#ifdef WEC_VF
& OCEAN(ng) % zetat, &
#endif
#ifdef ATM_PRESS
& FORCES(ng) % Pair, &
#endif
#ifdef DIAGNOSTICS_UV
& DIAGS(ng) % DiaRU, &
& DIAGS(ng) % DiaRV, &
#endif
& OCEAN(ng) % ru, &
& OCEAN(ng) % rv)
#ifdef PROFILE
CALL wclock_off (ng, iNLM, 23, __LINE__, MyFile)
#endif
!
RETURN
END SUBROUTINE prsgrd
!
!***********************************************************************
SUBROUTINE prsgrd31_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs, &
#ifdef WET_DRY
& umask_wet, vmask_wet, &
#endif
& Hz, om_v, on_u, z_r, z_w, &
& rho, &
#ifdef TIDE_GENERATING_FORCES
& eq_tide, &
#endif
#ifdef WEC_VF
& zetat, &
#endif
#ifdef ATM_PRESS
& Pair, &
#endif
#ifdef DIAGNOSTICS_UV
& DiaRU, DiaRV, &
#endif
& ru, rv)
!***********************************************************************
!
USE mod_param
USE mod_scalars
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: nrhs
#ifdef ASSUMED_SHAPE
# ifdef WET_DRY
real(r8), intent(in) :: umask_wet(LBi:,LBj:)
real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
# endif
real(r8), intent(in) :: Hz(LBi:,LBj:,:)
real(r8), intent(in) :: om_v(LBi:,LBj:)
real(r8), intent(in) :: on_u(LBi:,LBj:)
real(r8), intent(in) :: z_r(LBi:,LBj:,:)
real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
real(r8), intent(in) :: rho(LBi:,LBj:,:)
# ifdef TIDE_GENERATING_FORCES
real(r8), intent(in) :: eq_tide(LBi:,LBj:)
# endif
# ifdef WEC_VF
real(r8), intent(in) :: zetat(LBi:,LBj:)
# endif
# ifdef ATM_PRESS
real(r8), intent(in) :: Pair(LBi:,LBj:)
# endif
# ifdef DIAGNOSTICS_UV
real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
# endif
real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
#else
# ifdef WET_DRY
real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
# ifdef TIDE_GENERATING_FORCES
real(r8), intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
# endif
# ifdef WEC_VF
real(r8), intent(in) :: zetat(LBi:UBi,LBj:UBj)
# endif
# ifdef ATM_PRESS
real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
# endif
# ifdef DIAGNOSTICS_UV
real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
# endif
real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
#endif
!
! Local variable declarations.
!
integer :: i, j, k
real(r8) :: fac, fac1, fac2, fac3
real(r8) :: cff1, cff2, cff3, cff4
#ifdef WJ_GRADP
real(r8) :: gamma
#endif
real(r8), dimension(IminS:ImaxS) :: phie
real(r8), dimension(IminS:ImaxS) :: phix
#include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Calculate pressure gradient in the XI-direction (m4/s2).
!-----------------------------------------------------------------------
!
! Compute surface baroclinic pressure gradient.
!
#ifdef ATM_PRESS
fac=100.0_r8/rho0
#endif
fac1=0.5_r8*g/rho0
fac2=1000.0_r8*g/rho0
fac3=0.25_r8*g/rho0
DO j=Jstr,Jend
DO i=IstrU,Iend
#ifdef TIDE_GENERATING_FORCES
cff1=z_w(i ,j,N(ng))-eq_tide(i ,j)-z_r(i ,j,N(ng))+ &
& z_w(i-1,j,N(ng))-eq_tide(i-1,j)-z_r(i-1,j,N(ng))
#else
cff1=z_w(i ,j,N(ng))-z_r(i ,j,N(ng))+ &
& z_w(i-1,j,N(ng))-z_r(i-1,j,N(ng))
#endif
phix(i)=fac1*(rho(i,j,N(ng))-rho(i-1,j,N(ng)))*cff1
#ifdef WEC_VF
phix(i)=phix(i)+zetat(i,j)-zetat(i-1,j)
#endif
#ifdef ATM_PRESS
phix(i)=phix(i)+fac*(Pair(i,j)-Pair(i-1,j))
#endif
#ifdef RHO_SURF
phix(i)=phix(i)+ &
& (fac2+fac1*(rho(i,j,N(ng))+rho(i-1,j,N(ng))))* &
& (z_w(i,j,N(ng))-z_w(i-1,j,N(ng)))
#endif
ru(i,j,N(ng),nrhs)=-0.5_r8*(Hz(i,j,N(ng))+Hz(i-1,j,N(ng)))* &
& phix(i)*on_u(i,j)
#ifdef WET_DRY
ru(i,j,N(ng),nrhs)=ru(i,j,N(ng),nrhs)*umask_wet(i,j)
#endif
#ifdef DIAGNOSTICS_UV
DiaRU(i,j,N(ng),nrhs,M3pgrd)=ru(i,j,N(ng),nrhs)
#endif
END DO
!
! Compute interior baroclinic pressure gradient. Differentiate and
! then vertically integrate.
!
DO k=N(ng)-1,1,-1
DO i=IstrU,Iend
#ifdef WJ_GRADP
cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* &
& (z_r(i-1,j,k+1)-z_r(i-1,j,k)))
cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ &
& z_r(i ,j,k+1)-z_r(i-1,j,k+1)
cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- &
& z_r(i-1,j,k+1)+z_r(i-1,j,k )
gamma=0.125_r8*cff1*cff2*cff3
cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i-1,j,k+1))+ &
& (1.0_r8-gamma)*(rho(i,j,k )-rho(i-1,j,k ))
cff2=rho(i,j,k+1)+rho(i-1,j,k+1)- &
& rho(i,j,k )-rho(i-1,j,k )
cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
& z_r(i,j,k )-z_r(i-1,j,k )
cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ &
& (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k ))
phix(i)=phix(i)+ &
& fac3*(cff1*cff3-cff2*cff4)
#else
cff1=rho(i,j,k+1)-rho(i-1,j,k+1)+ &
& rho(i,j,k )-rho(i-1,j,k )
cff2=rho(i,j,k+1)+rho(i-1,j,k+1)- &
& rho(i,j,k )-rho(i-1,j,k )
cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
& z_r(i,j,k )-z_r(i-1,j,k )
cff4=z_r(i,j,k+1)-z_r(i-1,j,k+1)+ &
& z_r(i,j,k )-z_r(i-1,j,k )
phix(i)=phix(i)+ &
& fac3*(cff1*cff3-cff2*cff4)
#endif
ru(i,j,k,nrhs)=-0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))* &
& phix(i)*on_u(i,j)
#ifdef WET_DRY
ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
#endif
#ifdef DIAGNOSTICS_UV
DiaRU(i,j,k,nrhs,M3pgrd)=ru(i,j,k,nrhs)
#endif
END DO
END DO
!
!-----------------------------------------------------------------------
! Calculate pressure gradient in the ETA-direction (m4/s2).
!-----------------------------------------------------------------------
!
! Compute surface baroclinic pressure gradient.
!
IF (j.ge.JstrV) THEN
DO i=Istr,Iend
#ifdef TIDE_GENERATING_FORCES
cff1=z_w(i,j ,N(ng))-eq_tide(i,j )-z_r(i,j ,N(ng))+ &
& z_w(i,j-1,N(ng))-eq_tide(i,j-1)-z_r(i,j-1,N(ng))
#else
cff1=z_w(i,j ,N(ng))-z_r(i,j ,N(ng))+ &
& z_w(i,j-1,N(ng))-z_r(i,j-1,N(ng))
#endif
phie(i)=fac1*(rho(i,j,N(ng))-rho(i,j-1,N(ng)))*cff1
#ifdef WEC_VF
phie(i)=phie(i)+zetat(i,j)-zetat(i,j-1)
#endif
#ifdef ATM_PRESS
phie(i)=phie(i)+fac*(Pair(i,j)-Pair(i,j-1))
#endif
#ifdef RHO_SURF
phie(i)=phie(i)+ &
& (fac2+fac1*(rho(i,j,N(ng))+rho(i,j-1,N(ng))))* &
& (z_w(i,j,N(ng))-z_w(i,j-1,N(ng)))
#endif
rv(i,j,N(ng),nrhs)=-0.5_r8*(Hz(i,j,N(ng))+Hz(i,j-1,N(ng)))* &
& phie(i)*om_v(i,j)
#ifdef WET_DRY
rv(i,j,N(ng),nrhs)=rv(i,j,N(ng),nrhs)*vmask_wet(i,j)
#endif
#ifdef DIAGNOSTICS_UV
DiaRV(i,j,N(ng),nrhs,M3pgrd)=rv(i,j,N(ng),nrhs)
#endif
END DO
!
! Compute interior baroclinic pressure gradient. Differentiate and
! then vertically integrate.
!
DO k=N(ng)-1,1,-1
DO i=Istr,Iend
#ifdef WJ_GRADP
cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* &
& (z_r(i,j-1,k+1)-z_r(i,j-1,k)))
cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ &
& z_r(i,j ,k+1)-z_r(i,j-1,k+1)
cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- &
& z_r(i,j-1,k+1)+z_r(i,j-1,k )
gamma=0.125_r8*cff1*cff2*cff3
cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i,j-1,k+1))+ &
& (1.0_r8-gamma)*(rho(i,j,k )-rho(i,j-1,k ))
cff2=rho(i,j,k+1)+rho(i,j-1,k+1)- &
& rho(i,j,k )-rho(i,j-1,k )
cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
& z_r(i,j,k )-z_r(i,j-1,k )
cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ &
& (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k ))
phie(i)=phie(i)+ &
& fac3*(cff1*cff3-cff2*cff4)
#else
cff1=rho(i,j,k+1)-rho(i,j-1,k+1)+ &
& rho(i,j,k )-rho(i,j-1,k )
cff2=rho(i,j,k+1)+rho(i,j-1,k+1)- &
& rho(i,j,k )-rho(i,j-1,k )
cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
& z_r(i,j,k )-z_r(i,j-1,k )
cff4=z_r(i,j,k+1)-z_r(i,j-1,k+1)+ &
& z_r(i,j,k )-z_r(i,j-1,k )
phie(i)=phie(i)+ &
& fac3*(cff1*cff3-cff2*cff4)
#endif
rv(i,j,k,nrhs)=-0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k))* &
& phie(i)*om_v(i,j)
#ifdef WET_DRY
rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
#endif
#ifdef DIAGNOSTICS_UV
DiaRV(i,j,k,nrhs,M3pgrd)=rv(i,j,k,nrhs)
#endif
END DO
END DO
END IF
END DO
!
RETURN
END SUBROUTINE prsgrd31_tile
END MODULE prsgrd_mod