model/src/post_cg3d.F

C $Header: /u/gcmpack/MITgcm/model/src/post_cg3d.F,v 1.7 2016/05/04 22:10:37 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: POST_CG3D
C     !INTERFACE:
      SUBROUTINE POST_CG3D(
     I                      zeroPsNH, zeroMeanPnh,
     I                      myTime, myIter, myThid )

C     !DESCRIPTION:
C     Called from SOLVE_FOR_PRESSURE, after 3-D solver (cg3d):
C     Finish computation of Non-hydrostatic pressure from 3-D solver solution

C     !USES:
      IMPLICIT NONE
C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
c#include "FFIELDS.h"
#include "DYNVARS.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "NH_VARS.h"
#endif

C     === Functions ====
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     zeroPsNH    :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
C     zeroMeanPnh :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
C     myTime      :: Current time in simulation
C     myIter      :: Current iteration number in simulation
C     myThid      :: My Thread Id. number
      LOGICAL zeroPsNH, zeroMeanPnh
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j,k,bi,bj
      INTEGER ks
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL     locGamma
CEOP

C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
      IF ( nonHydrostatic .AND. exactConserv ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

         IF ( selectNHfreeSurf.GE.1 ) THEN
          DO j=1,sNy
           DO i=1,sNx
            locGamma = drC(1)*recip_Bo(i,j,bi,bj)
     &               /( deltaTMom*deltaTFreeSurf
     &                 *implicitNHPress*implicDiv2DFlow )
            ks = 1
c           ks = kSurfC(i,j,bi,bj)
c           IF ( ks.LE.Nr ) THEN
             dPhiNH(i,j,bi,bj) = ( phi_nh(i,j,ks,bi,bj)
     &           + locGamma*Bo_surf(i,j,bi,bj)
     &                     *implicDiv2DFlow*deltaTFreeSurf
c    &                     *( wVel(i,j,ks,bi,bj) - wSurfP2d(i,j) )
     &                     *( wVel(i,j,ks,bi,bj) - dPhiNH(i,j,bi,bj) )
     &                           )/(1. _d 0 + locGamma )
c           ENDIF
           ENDDO
          ENDDO
         ELSEIF ( uniformFreeSurfLev ) THEN
C-       Z coordinate: assume surface @ level k=1
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
           ENDDO
          ENDDO
         ELSE
C-       Other than Z coordinate: no assumption on surface level index
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            ks = kSurfC(i,j,bi,bj)
            IF ( ks.LE.Nr ) THEN
             dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
            ELSE
             dPhiNH(i,j,bi,bj) = 0.
            ENDIF
           ENDDO
          ENDDO
         ENDIF

        ENDDO
       ENDDO
       IF ( selectNHfreeSurf.GE.1 .AND.
     &  ( implicitNHPress.LT.oneRL .OR. selectP_inEOS_Zc.EQ.3 ) ) THEN
         CALL EXCH_XY_RL( dPhiNH, myThid )
       ENDIF
      ENDIF

C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
      IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
       IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
        CALL WRITE_FLD_XYZ_RL( 'cg3d_x','I10', phi_nh, myIter, myThid )
       ENDIF
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

         DO k=1,Nr
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            phi_nh(i,j,k,bi,bj) = ( phi_nh(i,j,k,bi,bj)
     &                            - dPhiNH(i,j,bi,bj)
     &                            )*maskC(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
            etaN(i,j,bi,bj) = etaN(i,j,bi,bj)
     &                      + recip_Bo(i,j,bi,bj)*dPhiNH(i,j,bi,bj)
            dPhiNH(i,j,bi,bj) = 0. _d 0
          ENDDO
         ENDDO

        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/model/src/post_cg3d.F,v 1.7 2016/05/04 22:10:37 jmc Exp $
2	C $Name: $
3
4	#include "CPP_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: POST_CG3D
8	C !INTERFACE:
9	SUBROUTINE POST_CG3D(
10	I zeroPsNH, zeroMeanPnh,
11	I myTime, myIter, myThid )
12
13	C !DESCRIPTION:
14	C Called from SOLVE_FOR_PRESSURE, after 3-D solver (cg3d):
15	C Finish computation of Non-hydrostatic pressure from 3-D solver solution
16
17	C !USES:
18	IMPLICIT NONE
19	C == Global variables
20	#include "SIZE.h"
21	#include "EEPARAMS.h"
22	#include "PARAMS.h"
23	#include "GRID.h"
24	#include "SURFACE.h"
25	c#include "FFIELDS.h"
26	#include "DYNVARS.h"
27	#ifdef ALLOW_NONHYDROSTATIC
28	#include "NH_VARS.h"
29	#endif
30
31	C === Functions ====
32	LOGICAL DIFFERENT_MULTIPLE
33	EXTERNAL DIFFERENT_MULTIPLE
34
35	C !INPUT/OUTPUT PARAMETERS:
36	C == Routine arguments ==
37	C zeroPsNH :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
38	C zeroMeanPnh :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
39	C myTime :: Current time in simulation
40	C myIter :: Current iteration number in simulation
41	C myThid :: My Thread Id. number
42	LOGICAL zeroPsNH, zeroMeanPnh
43	_RL myTime
44	INTEGER myIter
45	INTEGER myThid
46
47	#ifdef ALLOW_NONHYDROSTATIC
48	C !LOCAL VARIABLES:
49	C == Local variables ==
50	INTEGER i,j,k,bi,bj
51	INTEGER ks
52	c CHARACTER*(MAX_LEN_MBUF) msgBuf
53	_RL locGamma
54	CEOP
55
56	C-- Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
57	C from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
58	IF ( nonHydrostatic .AND. exactConserv ) THEN
59	DO bj=myByLo(myThid),myByHi(myThid)
60	DO bi=myBxLo(myThid),myBxHi(myThid)
61
62	IF ( selectNHfreeSurf.GE.1 ) THEN
63	DO j=1,sNy
64	DO i=1,sNx
65	locGamma = drC(1)*recip_Bo(i,j,bi,bj)
66	& /( deltaTMom*deltaTFreeSurf
67	& implicitNHPressimplicDiv2DFlow )
68	ks = 1
69	c ks = kSurfC(i,j,bi,bj)
70	c IF ( ks.LE.Nr ) THEN
71	dPhiNH(i,j,bi,bj) = ( phi_nh(i,j,ks,bi,bj)
72	& + locGamma*Bo_surf(i,j,bi,bj)
73	& implicDiv2DFlowdeltaTFreeSurf
74	c & *( wVel(i,j,ks,bi,bj) - wSurfP2d(i,j) )
75	& *( wVel(i,j,ks,bi,bj) - dPhiNH(i,j,bi,bj) )
76	& )/(1. _d 0 + locGamma )
77	c ENDIF
78	ENDDO
79	ENDDO
80	ELSEIF ( uniformFreeSurfLev ) THEN
81	C- Z coordinate: assume surface @ level k=1
82	DO j=1-OLy,sNy+OLy
83	DO i=1-OLx,sNx+OLx
84	dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
85	ENDDO
86	ENDDO
87	ELSE
88	C- Other than Z coordinate: no assumption on surface level index
89	DO j=1-OLy,sNy+OLy
90	DO i=1-OLx,sNx+OLx
91	ks = kSurfC(i,j,bi,bj)
92	IF ( ks.LE.Nr ) THEN
93	dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
94	ELSE
95	dPhiNH(i,j,bi,bj) = 0.
96	ENDIF
97	ENDDO
98	ENDDO
99	ENDIF
100
101	ENDDO
102	ENDDO
103	IF ( selectNHfreeSurf.GE.1 .AND.
104	& ( implicitNHPress.LT.oneRL .OR. selectP_inEOS_Zc.EQ.3 ) ) THEN
105	CALL EXCH_XY_RL( dPhiNH, myThid )
106	ENDIF
107	ENDIF
108
109	C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
110	IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
111	IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
112	CALL WRITE_FLD_XYZ_RL( 'cg3d_x','I10', phi_nh, myIter, myThid )
113	ENDIF
114	DO bj=myByLo(myThid),myByHi(myThid)
115	DO bi=myBxLo(myThid),myBxHi(myThid)
116
117	DO k=1,Nr
118	DO j=1-OLy,sNy+OLy
119	DO i=1-OLx,sNx+OLx
120	phi_nh(i,j,k,bi,bj) = ( phi_nh(i,j,k,bi,bj)
121	& - dPhiNH(i,j,bi,bj)
122	& )*maskC(i,j,k,bi,bj)
123	ENDDO
124	ENDDO
125	ENDDO
126	DO j=1-OLy,sNy+OLy
127	DO i=1-OLx,sNx+OLx
128	etaN(i,j,bi,bj) = etaN(i,j,bi,bj)
129	& + recip_Bo(i,j,bi,bj)*dPhiNH(i,j,bi,bj)
130	dPhiNH(i,j,bi,bj) = 0. _d 0
131	ENDDO
132	ENDDO
133
134	ENDDO
135	ENDDO
136	ENDIF
137	#endif /* ALLOW_NONHYDROSTATIC */
138
139	RETURN
140	END