model/src/correction_step.F

C $Header: /u/gcmpack/MITgcm/model/src/correction_step.F,v 1.28 2016/11/29 22:43:58 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: CORRECTION_STEP
C     !INTERFACE:
      SUBROUTINE CORRECTION_STEP( bi, bj, iMin, iMax, jMin, jMax,
     I                      phiSurfX, phiSurfY,
     I                      myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R CORRECTION_STEP
C     | o Corrects the horizontal flow fields with the surface
C     |   pressure (and Non-Hydrostatic pressure).
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"
#ifdef ALLOW_NONHYDROSTATIC
#include "NH_VARS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     bi, bj              :: Tile indices
C     iMin,iMax,jMin,jMax :: Loop counters range
C     phiSurfX, phiSurfY  :: Surface Potential gradient
C     myTime              :: Current time in simulation
C     myIter              :: Current iteration number in simulation
C     myThid              :: my Thread Id number
      _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER bi, bj
      INTEGER iMin, iMax, jMin, jMax
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     i, j         :: Loop counters
C     k            :: Level index
C     psFac, nhFac :: Scaling parameters for supressing gradients
C     gU_dpx       :: implicit part of pressure gradient tendency
C     gV_dpy       :: implicit part of pressure gradient tendency
      INTEGER i,j
      INTEGER k
      _RL     psFac, nhFac
      _RL     gU_dpx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gV_dpy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

C--   Loop over all layers, top to bottom
      DO k=1,Nr

#ifdef ALLOW_SOLVE4_PS_AND_DRAG
       IF ( selectImplicitDrag.EQ.2 ) THEN

C     On/off scaling parameter
        psFac = pfFacMom*implicSurfPress

C     Pressure gradient tendency (zonal mom): Implicit part
        DO j=jMin,jMax
         DO i=iMin,iMax
          gU_dpx(i,j) =
     &         -psFac*dU_psFacX(i,j,k,bi,bj)*phiSurfX(i,j)
c    &                  *_maskW(i,j,k,bi,bj)
         ENDDO
        ENDDO

C     Pressure gradient tendency (merid mom): Implicit part
        DO j=jMin,jMax
         DO i=iMin,iMax
          gV_dpy(i,j) =
     &         -psFac*dV_psFacY(i,j,k,bi,bj)*phiSurfY(i,j)
c    &                  *_maskS(i,j,k,bi,bj)
         ENDDO
        ENDDO

       ELSE
#endif /* ALLOW_SOLVE4_PS_AND_DRAG */

C     On/off scaling parameters (including anelastic & deep-model factors)
        psFac = pfFacMom*implicSurfPress
     &         *recip_deepFacC(k)*recip_rhoFacC(k)
        IF ( use3Dsolver ) THEN
         nhFac = pfFacMom*implicitNHPress
     &          *recip_deepFacC(k)*recip_rhoFacC(k)
        ELSE
         nhFac = 0.
        ENDIF

C     Pressure gradient tendency (zonal mom): Implicit part
        DO j=jMin,jMax
         DO i=iMin,iMax
          gU_dpx(i,j) = -(
     &          psFac*phiSurfX(i,j)
#ifdef ALLOW_NONHYDROSTATIC
     &        + nhFac*_recip_dxC(i,j,bi,bj)
     &           *(phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
#endif
     &                   )*_maskW(i,j,k,bi,bj)
         ENDDO
        ENDDO

C     Pressure gradient tendency (merid mom): Implicit part
        DO j=jMin,jMax
         DO i=iMin,iMax
          gV_dpy(i,j) = -(
     &          psFac*phiSurfY(i,j)
#ifdef ALLOW_NONHYDROSTATIC
     &        + nhFac*_recip_dyC(i,j,bi,bj)
     &           *(phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
#endif
     &                   )*_maskS(i,j,k,bi,bj)
         ENDDO
        ENDDO

#ifdef ALLOW_SOLVE4_PS_AND_DRAG
       ENDIF
#endif /* ALLOW_SOLVE4_PS_AND_DRAG */

C     Update zonal velocity: add implicit pressure gradient tendency
       DO j=jMin,jMax
        DO i=iMin,iMax
          uVel(i,j,k,bi,bj)=( gU(i,j,k,bi,bj)
     &                      + deltaTMom*gU_dpx(i,j)
     &                      )*_maskW(i,j,k,bi,bj)
#ifdef ALLOW_OBCS
     &                       *maskInW(i,j,bi,bj)
#endif
        ENDDO
       ENDDO

C     Update merid. velocity: add implicit pressure gradient tendency
       DO j=jMin,jMax
        DO i=iMin,iMax
          vVel(i,j,k,bi,bj)=( gV(i,j,k,bi,bj)
     &                      + deltaTMom*gV_dpy(i,j)
     &                      )*_maskS(i,j,k,bi,bj)
#ifdef ALLOW_OBCS
     &                       *maskInS(i,j,bi,bj)
#endif
        ENDDO
       ENDDO

#ifdef ALLOW_DIAGNOSTICS
c      IF ( useDiagnostics ) THEN
c       CALL DIAGNOSTICS_FILL( gU_dpx,
c    &                         'UDIAG7  ', k, 1, 2, bi, bj, myThid )
c       CALL DIAGNOSTICS_FILL( gV_dpy,
c    &                         'UDIAG8  ', k, 1, 2, bi, bj, myThid )
c      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C-    end of k loop
      ENDDO

      RETURN
      END
1	jmc	1.29	C $Header: /u/gcmpack/MITgcm/model/src/correction_step.F,v 1.28 2016/11/29 22:43:58 jmc Exp $
2	jmc	1.15	C $Name: $
3	adcroft	1.1
4	jmc	1.26	#include "PACKAGES_CONFIG.h"
5	cnh	1.10	#include "CPP_OPTIONS.h"
6	adcroft	1.1
7	cnh	1.19	CBOP
8			C !ROUTINE: CORRECTION_STEP
9			C !INTERFACE:
10	adcroft	1.1	SUBROUTINE CORRECTION_STEP( bi, bj, iMin, iMax, jMin, jMax,
11	jmc	1.29	I phiSurfX, phiSurfY,
12			I myTime, myIter, myThid )
13	cnh	1.19	C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.23	C \| S/R CORRECTION_STEP
16			C \| o Corrects the horizontal flow fields with the surface
17			C \| pressure (and Non-Hydrostatic pressure).
18	cnh	1.19	C ==========================================================
19			C \ev
20
21			C !USES:
22	cnh	1.7	IMPLICIT NONE
23	heimbach	1.12	C == Global variables ==
24	adcroft	1.1	#include "SIZE.h"
25	cnh	1.3	#include "EEPARAMS.h"
26	adcroft	1.1	#include "PARAMS.h"
27			#include "GRID.h"
28	jmc	1.22	#include "DYNVARS.h"
29	adcroft	1.11	#ifdef ALLOW_NONHYDROSTATIC
30	jmc	1.22	#include "NH_VARS.h"
31	adcroft	1.11	#endif
32	cnh	1.19
33			C !INPUT/OUTPUT PARAMETERS:
34	adcroft	1.1	C == Routine Arguments ==
35	jmc	1.25	C bi, bj :: Tile indices
36			C iMin,iMax,jMin,jMax :: Loop counters range
37			C phiSurfX, phiSurfY :: Surface Potential gradient
38			C myTime :: Current time in simulation
39	jmc	1.29	C myIter :: Current iteration number in simulation
40	jmc	1.25	C myThid :: my Thread Id number
41			_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43	jmc	1.29	INTEGER bi, bj
44			INTEGER iMin, iMax, jMin, jMax
45	jmc	1.25	_RL myTime
46	jmc	1.29	INTEGER myIter
47	adcroft	1.1	INTEGER myThid
48	cnh	1.7
49	cnh	1.19	C !LOCAL VARIABLES:
50	adcroft	1.1	C == Local variables ==
51	jmc	1.29	C i, j :: Loop counters
52			C k :: Level index
53	jmc	1.25	C psFac, nhFac :: Scaling parameters for supressing gradients
54	jmc	1.29	C gU_dpx :: implicit part of pressure gradient tendency
55			C gV_dpy :: implicit part of pressure gradient tendency
56	adcroft	1.1	INTEGER i,j
57	jmc	1.29	INTEGER k
58	jmc	1.25	_RL psFac, nhFac
59	jmc	1.29	_RL gU_dpx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60			_RL gV_dpy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61	cnh	1.19	CEOP
62	adcroft	1.1
63	jmc	1.29	C-- Loop over all layers, top to bottom
64			DO k=1,Nr
65
66	jmc	1.27	#ifdef ALLOW_SOLVE4_PS_AND_DRAG
67	jmc	1.29	IF ( selectImplicitDrag.EQ.2 ) THEN
68	jmc	1.27
69	jmc	1.29	C On/off scaling parameter
70			psFac = pfFacMom*implicSurfPress
71	jmc	1.27
72	jmc	1.29	C Pressure gradient tendency (zonal mom): Implicit part
73			DO j=jMin,jMax
74			DO i=iMin,iMax
75			gU_dpx(i,j) =
76			& -psFacdU_psFacX(i,j,k,bi,bj)phiSurfX(i,j)
77			c & *_maskW(i,j,k,bi,bj)
78			ENDDO
79	jmc	1.27	ENDDO
80
81	jmc	1.29	C Pressure gradient tendency (merid mom): Implicit part
82			DO j=jMin,jMax
83			DO i=iMin,iMax
84			gV_dpy(i,j) =
85			& -psFacdV_psFacY(i,j,k,bi,bj)phiSurfY(i,j)
86			c & *_maskS(i,j,k,bi,bj)
87			ENDDO
88	jmc	1.27	ENDDO
89
90	jmc	1.29	ELSE
91	jmc	1.27	#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
92
93	jmc	1.24	C On/off scaling parameters (including anelastic & deep-model factors)
94	jmc	1.29	psFac = pfFacMom*implicSurfPress
95	jmc	1.25	& recip_deepFacC(k)recip_rhoFacC(k)
96	jmc	1.29	IF ( use3Dsolver ) THEN
97			nhFac = pfFacMom*implicitNHPress
98			& recip_deepFacC(k)recip_rhoFacC(k)
99			ELSE
100			nhFac = 0.
101			ENDIF
102
103			C Pressure gradient tendency (zonal mom): Implicit part
104			DO j=jMin,jMax
105			DO i=iMin,iMax
106			gU_dpx(i,j) = -(
107			& psFac*phiSurfX(i,j)
108			#ifdef ALLOW_NONHYDROSTATIC
109			& + nhFac*_recip_dxC(i,j,bi,bj)
110			& *(phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
111			#endif
112			& )*_maskW(i,j,k,bi,bj)
113			ENDDO
114			ENDDO
115
116			C Pressure gradient tendency (merid mom): Implicit part
117			DO j=jMin,jMax
118			DO i=iMin,iMax
119			gV_dpy(i,j) = -(
120			& psFac*phiSurfY(i,j)
121			#ifdef ALLOW_NONHYDROSTATIC
122			& + nhFac*_recip_dyC(i,j,bi,bj)
123			& *(phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
124			#endif
125			& )*_maskS(i,j,k,bi,bj)
126			ENDDO
127			ENDDO
128
129			#ifdef ALLOW_SOLVE4_PS_AND_DRAG
130	jmc	1.27	ENDIF
131	jmc	1.29	#endif /* ALLOW_SOLVE4_PS_AND_DRAG */
132	adcroft	1.1
133	jmc	1.29	C Update zonal velocity: add implicit pressure gradient tendency
134	jmc	1.27	DO j=jMin,jMax
135			DO i=iMin,iMax
136	jmc	1.29	uVel(i,j,k,bi,bj)=( gU(i,j,k,bi,bj)
137			& + deltaTMom*gU_dpx(i,j)
138			& )*_maskW(i,j,k,bi,bj)
139	jmc	1.26	#ifdef ALLOW_OBCS
140	jmc	1.29	& *maskInW(i,j,bi,bj)
141	jmc	1.26	#endif
142	jmc	1.27	ENDDO
143	adcroft	1.1	ENDDO
144
145	jmc	1.29	C Update merid. velocity: add implicit pressure gradient tendency
146	jmc	1.27	DO j=jMin,jMax
147			DO i=iMin,iMax
148	jmc	1.29	vVel(i,j,k,bi,bj)=( gV(i,j,k,bi,bj)
149			& + deltaTMom*gV_dpy(i,j)
150			& )*_maskS(i,j,k,bi,bj)
151	jmc	1.26	#ifdef ALLOW_OBCS
152	jmc	1.29	& *maskInS(i,j,bi,bj)
153	jmc	1.26	#endif
154	jmc	1.27	ENDDO
155	adcroft	1.4	ENDDO
156	jmc	1.27
157	jmc	1.29	#ifdef ALLOW_DIAGNOSTICS
158			c IF ( useDiagnostics ) THEN
159			c CALL DIAGNOSTICS_FILL( gU_dpx,
160			c & 'UDIAG7 ', k, 1, 2, bi, bj, myThid )
161			c CALL DIAGNOSTICS_FILL( gV_dpy,
162			c & 'UDIAG8 ', k, 1, 2, bi, bj, myThid )
163			c ENDIF
164			#endif /* ALLOW_DIAGNOSTICS */
165
166			C- end of k loop
167			ENDDO
168	adcroft	1.1
169			RETURN
170			END