#include "cppdefs.h"
#undef AC_U3HADVECTION
#define AC_HSIMT
MODULE prestep_inw_mod
#if defined INWAVE_MODEL
!
!=======================================================================
! !
! This subroutine initialize computations for new time step of the !
! inwave model. !
! !
!=======================================================================
!
implicit none
PRIVATE
PUBLIC :: prestep_inw
CONTAINS
!
!***********************************************************************
SUBROUTINE prestep_inw (ng, tile)
!***********************************************************************
!
USE mod_param
USE mod_grid
USE mod_stepping
USE mod_ocean
USE mod_inwave_vars
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
!
! Local variable declarations.
!
# include "tile.h"
!
# ifdef PROFILE
CALL wclock_on (ng, iNLM, 22)
# endif
CALL prestep_inw_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs(ng), nstp(ng), nnew(ng), &
# ifdef MASKING
& GRID(ng) % rmask, &
& GRID(ng) % umask, &
& GRID(ng) % vmask, &
# endif
& GRID(ng) % pm, &
& GRID(ng) % pn, &
& GRID(ng) % on_u, &
& GRID(ng) % om_v, &
& OCEAN(ng) % u, &
& OCEAN(ng) % v, &
& WAVEP(ng) % AC, &
& WAVEP(ng) % cx, &
& WAVEP(ng) % cy, &
& WAVEP(ng) % ct, &
& WAVEP(ng) % Tr, &
& WAVEP(ng) % kwc, &
& WAVEG(ng) % pd)
# ifdef PROFILE
CALL wclock_off (ng, iNLM, 22)
# endif
RETURN
END SUBROUTINE prestep_inw
!
!***********************************************************************
SUBROUTINE prestep_inw_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs, nstp, nnew, &
# ifdef MASKING
& rmask, umask, vmask, &
# endif
& pm, pn, on_u, om_v, &
& u, v, &
& AC, cx, cy, ct, Tr, kwc, pd)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_coupling
USE mod_inwave_params
USE mod_inwave_vars
USE exchange_3d_mod, ONLY : exchange_AC3d_tile
# ifdef DISTRIBUTE
USE mp_exchange_mod, ONLY : mp_exchange3d
# endif
USE AC3dbc_mod, ONLY : AC3dbc_tile
!
! 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, nstp, nnew
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(in) :: on_u(LBi:,LBj:)
real(r8), intent(in) :: om_v(LBi:,LBj:)
real(r8), intent(in) :: u(LBi:,LBj:,:,:)
real(r8), intent(in) :: v(LBi:,LBj:,:,:)
real(r8), intent(inout) :: AC(LBi:,LBj:,:,:)
real(r8), intent(in) :: cx(LBi:,LBj:,:)
real(r8), intent(in) :: cy(LBi:,LBj:,:)
real(r8), intent(in) :: ct(LBi:,LBj:,:)
real(r8), intent(in) :: Tr(LBi:,LBj:,:)
real(r8), intent(in) :: kwc(LBi:,LBj:,:)
real(r8), intent(in) :: pd
# else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
real(r8), intent(inout) :: AC(LBi:UBi,LBj:UBj,ND,3)
real(r8), intent(in) :: cx(LBi:UBi,LBj:UBj,ND)
real(r8), intent(in) :: cy(LBi:UBi,LBj:UBj,ND)
real(r8), intent(in) :: ct(LBi:UBi,LBj:UBj,ND+1)
real(r8), intent(in) :: Tr(LBi:UBi,LBj:UBj,ND)
real(r8), intent(in) :: kwc(LBi:UBi,LBj:UBj,ND)
real(r8), intent(in) :: pd
# endif
!
! Local variable declarations.
!
integer :: i, indx, is, itrc, j, d, ltrc
# if defined AC_MPDATA || defined AC_HSIMT
real(r8), parameter :: Gamma = 0.5_r8
# else
real(r8), parameter :: Gamma = 1.0_r8/6.0_r8
# endif
real(r8), parameter :: eps = 1.0E-16_r8
real(r8) :: cff, cff1, cff2, cff3, cff4, opd
real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
real(r8), dimension(IminS:ImaxS,0:ND+2) :: FD
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curv
real(r8), dimension(IminS:ImaxS,0:ND+1) :: curvd
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute intermediate action density at n+1/2 time-step, AC(i,j,dir,3)
!-----------------------------------------------------------------------
!
! Compute time rate of change of intermediate AC due to
! horizontal advection.
!
D_LOOP: DO d=1,ND
# ifdef AC_U3HADVECTION
! Third-order, uptream-biased horizontal advective fluxes.
DO j=Jstr,Jend
DO i=Istrm1,Iendp2
FX(i,j)=AC(i ,j,d,nstp)- &
& AC(i-1,j,d,nstp)
# ifdef MASKING
FX(i,j)=FX(i,j)*umask(i,j)
# endif
END DO
END DO
IF (.not.(CompositeGrid(iwest,ng).or.EWperiodic(ng))) THEN
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=Jstr,Jend
FX(Istr-1,j)=FX(Istr,j)
END DO
END IF
END IF
IF (.not.(CompositeGrid(ieast,ng).or.EWperiodic(ng))) THEN
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=Jstr,Jend
FX(Iend+2,j)=FX(Iend+1,j)
END DO
END IF
END IF
!
DO j=Jstr,Jend
DO i=Istr-1,Iend+1
curv(i,j)=FX(i+1,j)-FX(i,j)
END DO
END DO
!
cff1=1.0_r8/6.0_r8
cff2=1.0_r8/3.0_r8
DO j=Jstr,Jend
DO i=Istr,Iend+1
cff=cx(i,j,d)*on_u(i,j)
FX(i,j)=cff*0.5_r8* &
& (AC(i-1,j,d,nstp)+ &
& AC(i ,j,d,nstp))- &
& cff1*(curv(i-1,j)*MAX(cff,0.0_r8)+ &
& curv(i ,j)*MIN(cff,0.0_r8))
END DO
END DO
!
DO j=Jstrm1,Jendp2
DO i=Istr,Iend
FE(i,j)=AC(i,j ,d,nstp)- &
& AC(i,j-1,d,nstp)
# ifdef MASKING
FE(i,j)=FE(i,j)*vmask(i,j)
# endif
END DO
END DO
IF (.not.(CompositeGrid(isouth,ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,Iend
FE(i,Jstr-1)=FE(i,Jstr)
END DO
END IF
END IF
IF (.not.(CompositeGrid(inorth,ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,Iend
FE(i,Jend+2)=FE(i,Jend+1)
END DO
END IF
END IF
!
DO j=Jstr-1,Jend+1
DO i=Istr,Iend
curv(i,j)=FE(i,j+1)-FE(i,j)
END DO
END DO
!
cff1=1.0_r8/6.0_r8
cff2=1.0_r8/3.0_r8
DO j=Jstr,Jend+1
DO i=Istr,Iend
cff=cy(i,j,d)*om_v(i,j)
FE(i,j)=cff*0.5_r8* &
& (AC(i,j-1,d,nstp)+ &
& AC(i,j ,d,nstp))- &
& cff1*(curv(i,j-1)*MAX(cff,0.0_r8)+ &
& curv(i,j )*MIN(cff,0.0_r8))
END DO
END DO
# elif defined AC_MPDATA || defined AC_HSIMT
!
! First-order, upstream differences horizontal advective fluxes.
!
DO j=Jstr,Jend
DO i=Istr,Iend+1
cff=cx(i,j,d)*on_u(i,j)
cff1=MAX(cff,0.0_r8)
cff2=MIN(cff,0.0_r8)
FX(i,j)=cff1*AC(i-1,j,d,nstp)+ &
& cff2*AC(i ,j,d,nstp)
END DO
END DO
DO j=Jstr,Jend+1
DO i=Istr,Iend
cff=cy(i,j,d)*om_v(i,j)
cff1=MAX(cff,0.0_r8)
cff2=MIN(cff,0.0_r8)
FE(i,j)=cff1*AC(i,j-1,d,nstp)+ &
& cff2*AC(i,j ,d,nstp)
END DO
END DO
# endif
!
! Time-step horizontal advection.
!
IF (iic(ng).eq.ntfirst(ng)) THEN
cff=0.5_r8*dt(ng)
cff1=1.0_r8
cff2=0.0_r8
ELSE
cff=(1.0_r8-Gamma)*dt(ng)
cff1=0.5_r8+Gamma
cff2=0.5_r8-Gamma
END IF
DO j=Jstr,Jend
DO i=Istr,Iend
AC(i,j,d,3)=(cff1*AC(i,j,d,nstp)+ &
& cff2*AC(i,j,d,nnew))- &
& cff*pm(i,j)*pn(i,j)* &
& (FX(i+1,j)-FX(i,j)+ &
& FE(i,j+1)-FE(i,j))
END DO
END DO
END DO D_LOOP
!
! Advection in theta space.
! Need to wrap around in theta dir. NOt always
!
opd=1.0_r8/pd
J_LOOP: DO j=Jstr,Jend
DO i=Istr,Iend
# if defined THETA_AC_PERIODIC
FD(i,0)=AC(i,j,ND ,nstp)- &
& AC(i,j,ND-1,nstp)
FD(i,1)=AC(i,j,1 ,nstp)- &
& AC(i,j,ND ,nstp)
# else
!!IN THIS POINT IT DOESNT MATTER THE BOUNDARY CONDITION,
!!WE JUST PUT IT AS IF IT WAS A NO GRADIENT
!!THE WALL BOUNDARY CONDITION WILL BE STABLISHED LATER
FD(i,0)=0.0_r8
FD(i,1)=0.0_r8
# endif
DO d=2,ND
FD(i,d)=AC(i,j,d ,nstp)- &
& AC(i,j,d-1,nstp)
END DO
# if defined THETA_AC_PERIODIC
FD(i,ND+1)=FD(i,1)
FD(i,ND+2)=FD(i,2)
# else
FD(i,ND+1)=0.0_r8
FD(i,ND+2)=0.0_r8
# endif
END DO
!
DO i=Istr,Iend
DO d=0,ND+1
curvd(i,d)=FD(i,d+1)-FD(i,d)
END DO
END DO
!
cff1=1.0_r8/6.0_r8
cff2=1.0_r8/3.0_r8
DO i=Istr,Iend
DO d=1,1
# if defined THETA_AC_PERIODIC
cff=ct(i,j,d)*opd
# else
# if defined THETA_AC_WALL
cff=0.0_r8
# else
cff=ct(i,j,d)*opd
# endif
# endif
FD(i,d)=cff*0.5_r8* &
# if defined THETA_AC_PERIODIC
& (AC(i,j,ND,nstp)+ &
& AC(i,j,d ,nstp))- &
# else
& (AC(i,j,d ,nstp)+ &
& AC(i,j,d ,nstp))- &
# endif
& cff1*(curvd(i,d-1)*MAX(cff,0.0_r8)+ &
& curvd(i,d )*MIN(cff,0.0_r8))
END DO
DO d=2,ND
cff=ct(i,j,d)*opd
FD(i,d)=cff*0.5_r8* &
& (AC(i,j,d-1,nstp)+ &
& AC(i,j,d ,nstp))- &
& cff1*(curvd(i,d-1)*MAX(cff,0.0_r8)+ &
& curvd(i,d )*MIN(cff,0.0_r8))
END DO
DO d=ND+1,ND+1
# if defined THETA_AC_PERIODIC
cff=ct(i,j,d)*opd
# else
# if defined THETA_AC_WALL
cff=0.0_r8
# else
cff=ct(i,j,d)*opd
# endif
# endif
FD(i,d)=cff*0.5_r8* &
# if defined THETA_AC_PERIODIC
& (AC(i,j,ND,nstp)+ &
& AC(i,j,1 ,nstp))- &
# else
& (AC(i,j,ND,nstp)+ &
& AC(i,j,ND,nstp))- &
# endif
& cff1*(curvd(i,d-1)*MAX(cff,0.0_r8)+ &
& curvd(i,d )*MIN(cff,0.0_r8))
END DO
END DO
!
! Time-step directional advection.
!
IF (iic(ng).eq.ntfirst(ng)) THEN
cff=0.5_r8*dt(ng)
ELSE
cff=(1.0_r8-Gamma)*dt(ng)
END IF
DO d=1,ND
DO i=Istr,Iend
AC(i,j,d,3)=AC(i,j,d,3)- &
& cff*pd* &
& (FD(i,d+1)-FD(i,d))
AC(i,j,d,3)=MAX(0.0_r8,AC(i,j,d,3))
END DO
END DO
END DO J_LOOP
!
!=======================================================================
! Apply lateral boundary conditions.
!=======================================================================
!
! Apply no periodic boundary conditions.
CALL AC3dbc_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nstp, 3, &
& AC)
!
! Apply periodic boundary conditions.
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL exchange_AC3d_tile (ng, tile, &
& LBi, UBi, LBj, UBj, 1, ND, &
& AC(:,:,:,3))
END IF
# ifdef DISTRIBUTE
CALL mp_exchange3d (ng, tile, iNLM, 1, &
& LBi, UBi, LBj, UBj, 1, ND, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& AC(:,:,:,3))
# endif
#endif
RETURN
END SUBROUTINE prestep_inw_tile
END MODULE prestep_inw_mod