Page MenuHomec4science
Files F122198140

timestep.F
No OneTemporary
Actions

File Metadata

Created
Wed, Jul 16, 13:18

timestep.F
View Options

C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.60 2014/12/10 22:22:31 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#ifdef ALLOW_CD_CODE
#include "CD_CODE_OPTIONS.h"
#endif
CBOP
C !ROUTINE: TIMESTEP
C !INTERFACE:
SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
I guDissip, gvDissip,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R TIMESTEP
C | o Step model fields forward in time
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "RESTART.h"
#include "DYNVARS.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "NH_VARS.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
C phiSurfY :: or geopotential (atmos) in X and Y direction
C guDissip :: dissipation tendency (all explicit terms), u component
C gvDissip :: dissipation tendency (all explicit terms), v component
INTEGER bi,bj,iMin,iMax,jMin,jMax
INTEGER k
_RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL myTime
INTEGER myIter, myThid
C !LOCAL VARIABLES:
C == Local variables ==
C guExt :: forcing tendency, u component
C gvExt :: forcing tendency, v component
C gu_AB :: tendency increment from Adams-Bashforth, u component
C gv_AB :: tendency increment from Adams-Bashforth, v component
INTEGER i,j
_RL phxFac,phyFac, psFac
_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)
_RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gu_AB(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL gv_AB(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-- factors for gradient (X & Y directions) of Hydrostatic Potential
phxFac = pfFacMom
phyFac = pfFacMom
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C- Initialize local arrays (not really necessary for all but safer)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
guExt(i,j) = 0. _d 0
gvExt(i,j) = 0. _d 0
gUtmp(i,j) = 0. _d 0
gVtmp(i,j) = 0. _d 0
#ifdef ALLOW_CD_CODE
guCor(i,j) = 0. _d 0
gvCor(i,j) = 0. _d 0
#endif
ENDDO
ENDDO
IF ( momForcing ) THEN
C-- Collect forcing term in local array guExt,gvExt:
CALL APPLY_FORCING_U(
U guExt,
I iMin,iMax,jMin,jMax, k, bi,bj,
I myTime, myIter, myThid )
CALL APPLY_FORCING_V(
U gvExt,
I iMin,iMax,jMin,jMax, k, bi,bj,
I myTime, myIter, myThid )
ENDIF
IF ( .NOT.staggerTimeStep .AND. .NOT. implicitIntGravWave ) THEN
C-- Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
DO j=jMin,jMax
DO i=iMin,iMax
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
ENDDO
ENDDO
phxFac = 0.
phyFac = 0.
c ELSE
C-- Stagger time step: grad Phi_Hyp will be added later
ENDIF
C-- Dissipation term inside the Adams-Bashforth:
IF ( momViscosity .AND. momDissip_In_AB) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
ENDDO
ENDDO
ENDIF
C-- Forcing term inside the Adams-Bashforth:
IF ( momForcing .AND. momForcingOutAB.NE.1 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guExt(i,j)
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvExt(i,j)
ENDDO
ENDDO
ENDIF
#ifdef CD_CODE_NO_AB_MOMENTUM
IF ( useCDscheme ) THEN
C- CD-scheme, before doing AB, store gU,Vtmp = gU,V^n (+dissip. +forcing)
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gU(i,j,k,bi,bj)
gVtmp(i,j) = gV(i,j,k,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* CD_CODE_NO_AB_MOMENTUM */
C- Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
C and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
#ifdef ALLOW_ADAMSBASHFORTH_3
CALL ADAMS_BASHFORTH3(
I bi, bj, k, Nr,
U gU(1-OLx,1-OLy,1,bi,bj), guNm,
O gu_AB,
I mom_StartAB, myIter, myThid )
CALL ADAMS_BASHFORTH3(
I bi, bj, k, Nr,
U gV(1-OLx,1-OLy,1,bi,bj), gvNm,
O gv_AB,
I mom_StartAB, myIter, myThid )
#else /* ALLOW_ADAMSBASHFORTH_3 */
CALL ADAMS_BASHFORTH2(
I bi, bj, k, Nr,
U gU(1-OLx,1-OLy,1,bi,bj),
U guNm1(1-OLx,1-OLy,1,bi,bj),
O gu_AB,
I mom_StartAB, myIter, myThid )
CALL ADAMS_BASHFORTH2(
I bi, bj, k, Nr,
U gV(1-OLx,1-OLy,1,bi,bj),
U gvNm1(1-OLx,1-OLy,1,bi,bj),
O gv_AB,
I mom_StartAB, myIter, myThid )
#endif /* ALLOW_ADAMSBASHFORTH_3 */
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(gu_AB,'AB_gU ',k,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(gv_AB,'AB_gV ',k,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C- Make a local copy in gU,Vtmp of gU,V^n+1/2 (+dissip. +forcing)
#ifdef CD_CODE_NO_AB_MOMENTUM
IF ( .NOT.useCDscheme ) THEN
#endif
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gU(i,j,k,bi,bj)
gVtmp(i,j) = gV(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef CD_CODE_NO_AB_MOMENTUM
ENDIF
#endif
C-- Forcing term outside the Adams-Bashforth:
IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j) + guExt(i,j)
gVtmp(i,j) = gVtmp(i,j) + gvExt(i,j)
ENDDO
ENDDO
ENDIF
C-- Dissipation term outside the Adams-Bashforth:
IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_CD_CODE
IF ( useCDscheme ) THEN
C- Step forward D-grid velocity using C-grid gU,Vtmp = gU,V +dissip +forcing
C and return coriolis terms on C-grid (guCor,gvCor)
CALL CD_CODE_SCHEME(
I bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
O guCor,gvCor,
I myTime, myIter, myThid)
#ifdef CD_CODE_NO_AB_MOMENTUM
IF ( momForcing .AND. momForcingOutAB.EQ.1 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = ( gU(i,j,k,bi,bj) + guExt(i,j) ) + guCor(i,j)
gVtmp(i,j) = ( gV(i,j,k,bi,bj) + gvExt(i,j) ) + gvCor(i,j)
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gU(i,j,k,bi,bj) + guCor(i,j)
gVtmp(i,j) = gV(i,j,k,bi,bj) + gvCor(i,j)
ENDDO
ENDDO
ENDIF
IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
ENDDO
ENDDO
ENDIF
#else /* CD_CODE_NO_AB_MOMENTUM */
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j) + guCor(i,j)
gVtmp(i,j) = gVtmp(i,j) + gvCor(i,j)
ENDDO
ENDDO
#endif /* CD_CODE_NO_AB_MOMENTUM */
ENDIF
#endif /* ALLOW_CD_CODE */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef NONLIN_FRSURF
IF ( .NOT. vectorInvariantMomentum
& .AND. nonlinFreeSurf.GT.1 ) THEN
IF ( select_rStar.GT.0 ) THEN
# ifndef DISABLE_RSTAR_CODE
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
ENDDO
ENDDO
# endif /* DISABLE_RSTAR_CODE */
ELSEIF ( selectSigmaCoord.NE.0 ) THEN
# ifndef DISABLE_SIGMA_CODE
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j)
& /( 1. _d 0 + dEtaWdt(i,j,bi,bj)*deltaTFreeSurf
& *dBHybSigF(k)*recip_drF(k)
& *recip_hFacW(i,j,k,bi,bj)
& )
gVtmp(i,j) = gVtmp(i,j)
& /( 1. _d 0 + dEtaSdt(i,j,bi,bj)*deltaTFreeSurf
& *dBHybSigF(k)*recip_drF(k)
& *recip_hFacS(i,j,k,bi,bj)
& )
ENDDO
ENDDO
# endif /* DISABLE_SIGMA_CODE */
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
IF ( k.EQ.kSurfW(i,j,bi,bj) ) THEN
gUtmp(i,j) = gUtmp(i,j)
& *_hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
ENDIF
IF ( k.EQ.kSurfS(i,j,bi,bj) ) THEN
gVtmp(i,j) = gVtmp(i,j)
& *_hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
ENDIF
ENDDO
ENDDO
ENDIF
ENDIF
#endif /* NONLIN_FRSURF */
#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
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(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) )
& )
gVtmp(i,j) = gVtmp(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) )
& )
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
gUtmp(i,j) = gUtmp(i,j)
& - nhFac*_recip_dxC(i,j,bi,bj)*
& (phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
gVtmp(i,j) = gVtmp(i,j)
& - nhFac*_recip_dyC(i,j,bi,bj)*
& (phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
ENDDO
ENDDO
ENDIF
ENDIF
#endif /* ALLOW_NONHYDROSTATIC */
C Step forward zonal velocity (store in Gu)
DO j=jMin,jMax
DO i=iMin,iMax
gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
& +deltaTMom*(
& gUtmp(i,j)
& - psFac*phiSurfX(i,j)
& - phxFac*dPhiHydX(i,j)
& )*_maskW(i,j,k,bi,bj)
ENDDO
ENDDO
C Step forward meridional velocity (store in Gv)
DO j=jMin,jMax
DO i=iMin,iMax
gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
& +deltaTMom*(
& gVtmp(i,j)
& - psFac*phiSurfY(i,j)
& - phyFac*dPhiHydY(i,j)
& )*_maskS(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_DIAGNOSTICS
IF ( momViscosity .AND. useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL( guDissip,'Um_Diss ',k,1,2,bi,bj,myThid )
CALL DIAGNOSTICS_FILL( gvDissip,'Vm_Diss ',k,1,2,bi,bj,myThid )
ENDIF
IF ( momForcing .AND. useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL( guExt,'Um_Ext ',k,1,2,bi,bj,myThid )
CALL DIAGNOSTICS_FILL( gvExt,'Vm_Ext ',k,1,2,bi,bj,myThid )
ENDIF
#ifdef ALLOW_CD_CODE
IF ( useCDscheme .AND. useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL( guCor,'Um_Cori ',k,1,2,bi,bj,myThid )
CALL DIAGNOSTICS_FILL( gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid )
ENDIF
#endif /* ALLOW_CD_CODE */
#endif /* ALLOW_DIAGNOSTICS */
RETURN
END

Event Timeline

c4science · Help