From 2c28ab55566bba1421bb23e2699f90e7ff03fa6e Mon Sep 17 00:00:00 2001 From: Jean-Michel Campin Date: Fri, 3 Jan 2025 12:46:41 -0500 Subject: [PATCH] implement solveForDeltaP=T option --- model/src/dynamics.F | 2 +- model/src/timestep.F | 41 ++++++++++++++++++++------------------- model/src/timestep_wvel.F | 31 +++++++++++++++++------------ 3 files changed, 41 insertions(+), 33 deletions(-) diff --git a/model/src/dynamics.F b/model/src/dynamics.F index 6847ed73bc..693980c610 100644 --- a/model/src/dynamics.F +++ b/model/src/dynamics.F @@ -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, diff --git a/model/src/timestep.F b/model/src/timestep.F index 024db1f536..acf7ab80ea 100644 --- a/model/src/timestep.F +++ b/model/src/timestep.F @@ -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) @@ -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 @@ -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) @@ -336,22 +339,16 @@ 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 @@ -359,10 +356,14 @@ SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k, 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 diff --git a/model/src/timestep_wvel.F b/model/src/timestep_wvel.F index a1f6edb3bd..021cf9102e 100644 --- a/model/src/timestep_wvel.F +++ b/model/src/timestep_wvel.F @@ -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) @@ -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 @@ -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