model/src/timestep_wvel.F

C $Header: /u/gcmpack/MITgcm/model/src/timestep_wvel.F,v 1.5 2009/11/29 03:12:32 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: TIMESTEP_WVEL
C     !INTERFACE:
      SUBROUTINE TIMESTEP_WVEL(
     I                          bi,bj, myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R TIMESTEP_WVEL
C     | o Step model vertical velocity 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 "DYNVARS.h"
#include "NH_VARS.h"
c #include "SURFACE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
      INTEGER bi,bj
      _RL     myTime
      INTEGER myIter, myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER iMin,iMax,jMin,jMax
      INTEGER i, j, k, km1
      _RL     gWtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     tmpFac, nh_Fac, igwFac
CEOP

      iMin = 1
      iMax = sNx
      jMin = 1
      jMax = sNy

      igwFac = 0. _d 0
      IF ( implicitIntGravWave ) igwFac = 1. _d 0

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF ( nonHydrostatic ) THEN
        nh_Fac = 0.
        IF ( nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2

        DO k=1,Nr
         km1 = MAX( k-1, 1 )

         IF ( implicitNHPress.NE.1. _d 0 ) THEN
C--   add explicit part of NH pressure gradient:
          tmpFac = pfFacMom*(1. _d 0 - implicitNHPress)
     &           * wUnit2rVel(k)*wUnit2rVel(k)*recip_rhoFacF(k)
          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) )
            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) )
            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
            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)
          ENDDO
         ENDDO
C     Step forward vertical velocity
         tmpFac = nh_Fac + igwFac*wUnit2rVel(k)*wUnit2rVel(k)
     &                           *dBdrRef(k)*deltaTMom*dTtracerLev(k)
         IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
         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)
          ENDDO
         ENDDO

C-    End of k loop
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_OBCS
C--   This call is aesthetic: it makes the W field
C     consistent with the OBs but this has no algorithmic
C     impact. This is purely for diagnostic purposes.
        IF (useOBCS) THEN
         DO k=1,Nr
          CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS */

      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
            gW(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO

C-    End if nonHydrostatic / elseif implicitIntGravWave
      ENDIF

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/model/src/timestep_wvel.F,v 1.5 2009/11/29 03:12:32 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: TIMESTEP_WVEL
9			C !INTERFACE:
10			SUBROUTINE TIMESTEP_WVEL(
11	jmc	1.3	I bi,bj, myTime, myIter, myThid )
12	jmc	1.1	C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| S/R TIMESTEP_WVEL
15			C \| o Step model vertical velocity forward in time
16			C ==========================================================
17			C \ev
18
19			C !USES:
20			IMPLICIT NONE
21			C == Global variables ==
22			#include "SIZE.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25			#include "GRID.h"
26			#include "DYNVARS.h"
27			#include "NH_VARS.h"
28			c #include "SURFACE.h"
29
30			C !INPUT/OUTPUT PARAMETERS:
31			C == Routine Arguments ==
32	jmc	1.3	INTEGER bi,bj
33	jmc	1.1	_RL myTime
34			INTEGER myIter, myThid
35
36			#ifdef ALLOW_NONHYDROSTATIC
37			C !LOCAL VARIABLES:
38			C == Local variables ==
39			INTEGER iMin,iMax,jMin,jMax
40	jmc	1.6	INTEGER i, j, k, km1
41	jmc	1.1	_RL gWtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42	jmc	1.2	_RL tmpFac, nh_Fac, igwFac
43	jmc	1.1	CEOP
44
45			iMin = 1
46			iMax = sNx
47			jMin = 1
48			jMax = sNy
49
50	jmc	1.4	igwFac = 0. _d 0
51			IF ( implicitIntGravWave ) igwFac = 1. _d 0
52	jmc	1.2
53	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
54
55	jmc	1.3	IF ( nonHydrostatic ) THEN
56	jmc	1.2	nh_Fac = 0.
57			IF ( nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2
58
59	jmc	1.6	DO k=1,Nr
60			km1 = MAX( k-1, 1 )
61	jmc	1.1
62	jmc	1.5	IF ( implicitNHPress.NE.1. _d 0 ) THEN
63			C-- add explicit part of NH pressure gradient:
64			tmpFac = pfFacMom*(1. _d 0 - implicitNHPress)
65	jmc	1.6	& * wUnit2rVel(k)wUnit2rVel(k)recip_rhoFacF(k)
66			IF ( k.GE.2 ) THEN
67			DO j=jMin,jMax
68			DO i=iMin,iMax
69			gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
70			& - tmpFacrkSignrecip_drC(k)
71			& *( phi_nh(i,j,k,bi,bj) - phi_nh(i,j,k-1,bi,bj) )
72			ENDDO
73	jmc	1.5	ENDDO
74	jmc	1.6	ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
75			DO j=jMin,jMax
76			DO i=iMin,iMax
77			gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
78			& - tmpFacrkSignrecip_drC(k)
79			& *( phi_nh(i,j,k,bi,bj) - dPhiNH(i,j,bi,bj) )
80			ENDDO
81			ENDDO
82			ENDIF
83	jmc	1.5	ENDIF
84	jmc	1.2	C apply mask to gW and keep a copy of wVel in gW:
85			DO j=1-Oly,sNy+Oly
86			DO i=1-Olx,sNx+Olx
87			gWtmp(i,j) = gW(i,j,k,bi,bj)
88	jmc	1.6	& maskC(i,j,k,bi,bj)maskC(i,j,km1,bi,bj)
89	jmc	1.2	gW(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
90	jmc	1.1	ENDDO
91			ENDDO
92	jmc	1.2	C Step forward vertical velocity
93	jmc	1.4	tmpFac = nh_Fac + igwFacwUnit2rVel(k)wUnit2rVel(k)
94			& dBdrRef(k)deltaTMom*dTtracerLev(k)
95	jmc	1.2	IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
96	jmc	1.1	DO j=jMin,jMax
97			DO i=iMin,iMax
98			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
99	jmc	1.2	& + deltaTmomtmpFacgWtmp(i,j)
100	jmc	1.1	ENDDO
101			ENDDO
102
103			C- End of k loop
104			ENDDO
105
106			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
107
108			#ifdef ALLOW_OBCS
109			C-- This call is aesthetic: it makes the W field
110			C consistent with the OBs but this has no algorithmic
111			C impact. This is purely for diagnostic purposes.
112			IF (useOBCS) THEN
113			DO k=1,Nr
114			CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
115			ENDDO
116			ENDIF
117			#endif /* ALLOW_OBCS */
118
119	jmc	1.3	ELSEIF ( implicitIntGravWave ) THEN
120	jmc	1.2	C keep a copy of wVel in gW:
121			DO k=1,Nr
122			DO j=1-Oly,sNy+Oly
123			DO i=1-Olx,sNx+Olx
124			gW(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
125			ENDDO
126			ENDDO
127			ENDDO
128
129			C- End if nonHydrostatic / elseif implicitIntGravWave
130	jmc	1.3	ENDIF
131	jmc	1.1
132			#endif /* ALLOW_NONHYDROSTATIC */
133
134			RETURN
135			END