Skip to content

Commit

Permalink
implement solveForDeltaP=T option
Browse files Browse the repository at this point in the history
  • Loading branch information
@jm-c
jm-c committed Jan 3, 2025
1 parent b4d96b0 commit 2c28ab5
Show file tree
Hide file tree
Showing 3 changed files with 41 additions and 33 deletions.
2 changes: 1 addition & 1 deletion model/src/dynamics.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,7 @@ SUBROUTINE DYNAMICS(myTime, myIter, myThid)

C-- Explicit part of the Surface Potential Gradient (add in TIMESTEP)
C (note: this loop will be replaced by CALL CALC_GRAD_ETA)
IF (implicSurfPress.NE.1.) THEN
IF ( solveForDeltaP .OR. implicSurfPress.NE.oneRL ) THEN
CALL CALC_GRAD_PHI_SURF(
I bi,bj,iMin,iMax,jMin,jMax,
I etaN,
Expand Down
41 changes: 21 additions & 20 deletions model/src/timestep.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
C gu_AB :: tendency increment from Adams-Bashforth, u component
C gv_AB :: tendency increment from Adams-Bashforth, v component
INTEGER i,j
_RL phFac, psFac
_RL phFac, psFac, nhFac
_RL guExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gvExt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
Expand All @@ -67,18 +67,21 @@ SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
_RL gv_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gUdPx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gVdPy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#ifdef ALLOW_NONHYDROSTATIC
_RL nhFac
#endif
#ifdef ALLOW_CD_CODE
_RL guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif
CEOP

C-- explicit part of the surface potential gradient is added in this S/R
psFac = pfFacMom*(1. _d 0 - implicSurfPress)
& *recip_deepFacC(k)*recip_rhoFacC(k)
C-- solveForDeltaP=T : explicit surface potential gradient is added in this S/R
C solveForDeltaP=F : only the explicit part (1-implicSurfPress) is added
C and same for Non-Hydrostatic pressure/potential horizontal gradient.
psFac = pfFacMom*recip_deepFacC(k)*recip_rhoFacC(k)
nhFac = psFac
IF ( .NOT.solveForDeltaP ) THEN
psFac = psFac*(1. _d 0 - implicSurfPress)
nhFac = nhFac*(1. _d 0 - implicitNHPress)
ENDIF

C-- factors for gradient (X & Y directions) of Hydrostatic Potential
phFac = pfFacMom
Expand Down Expand Up @@ -325,7 +328,7 @@ SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
gVdPy(i,j) = -phFac*dPhiHydY(i,j) - psFac*phiSurfY(i,j)
ENDDO
ENDDO
ELSEIF ( implicSurfPress.NE.oneRL ) THEN
ELSEIF ( psFac.NE.zeroRL ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUdPx(i,j) = -psFac*phiSurfX(i,j)
Expand All @@ -336,33 +339,31 @@ SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,

#ifdef ALLOW_NONHYDROSTATIC
C-- explicit part of the NH pressure gradient is added here
IF ( use3Dsolver .AND. implicitNHPress.NE.1. _d 0 ) THEN
nhFac = pfFacMom*(1. _d 0 - implicitNHPress)
& *recip_deepFacC(k)*recip_rhoFacC(k)
IF ( exactConserv ) THEN
IF ( use3Dsolver .AND. nhFac.NE.zeroRL ) THEN
IF ( solveForDeltaP .OR. .NOT.exactConserv ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUdPx(i,j) = gUdPx(i,j)
& - nhFac*_recip_dxC(i,j,bi,bj)
& *( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
& -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
& )
& *( phi_nh(i,j,k,bi,bj) - phi_nh(i-1,j,k,bi,bj) )
gVdPy(i,j) = gVdPy(i,j)
& - nhFac*_recip_dyC(i,j,bi,bj)
& *( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
& -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
& )
& *( phi_nh(i,j,k,bi,bj) - phi_nh(i,j-1,k,bi,bj) )
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
gUdPx(i,j) = gUdPx(i,j)
& - nhFac*_recip_dxC(i,j,bi,bj)
& *( phi_nh(i,j,k,bi,bj) - phi_nh(i-1,j,k,bi,bj) )
& *( (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
& -( dPhiNH(i,j,bi,bj) - dPhiNH(i-1,j,bi,bj) )
& )
gVdPy(i,j) = gVdPy(i,j)
& - nhFac*_recip_dyC(i,j,bi,bj)
& *( phi_nh(i,j,k,bi,bj) - phi_nh(i,j-1,k,bi,bj) )
& *( (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
& -( dPhiNH(i,j,bi,bj) - dPhiNH(i,j-1,bi,bj) )
& )
ENDDO
ENDDO
ENDIF
Expand Down
31 changes: 19 additions & 12 deletions model/src/timestep_wvel.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,31 +60,38 @@ SUBROUTINE TIMESTEP_WVEL(
DO k=1,Nr
km1 = MAX( k-1, 1 )

IF ( implicitNHPress.NE.1. _d 0 ) THEN
IF ( solveForDeltaP .OR. implicitNHPress.NE.oneRL ) THEN
C-- add explicit part of NH pressure gradient:
tmpFac = pfFacMom*(1. _d 0 - implicitNHPress)
& * wUnit2rVel(k)*wUnit2rVel(k)*recip_rhoFacF(k)
tmpFac = pfFacMom*wUnit2rVel(k)*wUnit2rVel(k)
& *recip_rhoFacF(k)*recip_drC(k)*rkSign
IF ( .NOT.solveForDeltaP )
& tmpFac = tmpFac*(1. _d 0 - implicitNHPress)
IF ( k.GE.2 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
& - tmpFac*rkSign*recip_drC(k)
& *( phi_nh(i,j,k,bi,bj) - phi_nh(i,j,k-1,bi,bj) )
& - tmpFac*( phi_nh(i,j,k,bi,bj) - phi_nh(i,j,k-1,bi,bj) )
ENDDO
ENDDO
ELSEIF ( selectNHfreeSurf.GE.1 .AND. solveForDeltaP ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
& - tmpFac * phi_nh(i,j,k,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
& - tmpFac*rkSign*recip_drC(k)
& *( phi_nh(i,j,k,bi,bj) - dPhiNH(i,j,bi,bj) )
& - tmpFac*( phi_nh(i,j,k,bi,bj) - dPhiNH(i,j,bi,bj) )
ENDDO
ENDDO
ENDIF
ENDIF
C apply mask to gW and keep a copy of wVel in gW:
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gWtmp(i,j) = gW(i,j,k,bi,bj)
& *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
gW(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
Expand All @@ -97,7 +104,7 @@ SUBROUTINE TIMESTEP_WVEL(
DO j=jMin,jMax
DO i=iMin,iMax
wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
& + deltaTmom*tmpFac*gWtmp(i,j)
& + deltaTMom*tmpFac*gWtmp(i,j)
ENDDO
ENDDO

Expand All @@ -109,8 +116,8 @@ SUBROUTINE TIMESTEP_WVEL(
ELSEIF ( implicitIntGravWave ) THEN
C keep a copy of wVel in gW:
DO k=1,Nr
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gW(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
ENDDO
ENDDO
Expand Down

0 comments on commit 2c28ab5

Please sign in to comment.