model/src/update_etah.F

C $Header: /u/gcmpack/MITgcm/model/src/update_etah.F,v 1.11 2011/04/27 22:20:22 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: UPDATE_ETAH
C     !INTERFACE:
      SUBROUTINE UPDATE_ETAH( myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE UPDATE_ETAH
C     | o Update etaH after mom-correction-step/integr_continuity
C     |  (required with NLFS to derive surface layer thickness)
C     | o Also derive SSH tendency at grid-cell Western and
C     |   Southern edges (for hybrid sigma-coordinate)
C     *==========================================================*
C     \ev

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

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime  :: Current time in simulation
C     myIter  :: Current iteration number in simulation
C     myThid  :: Thread number for this instance of the routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
#ifdef EXACT_CONSERV
C     Local variables in common block

C     Local variables
C     i,j,bi,bj  :: Loop counters
      INTEGER i,j,bi,bj
CEOP


      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

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

C--   before updating etaH, save current etaH field in etaHnm1
         DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
             etaHnm1(i,j,bi,bj) = etaH(i,j,bi,bj)
           ENDDO
         ENDDO

C--   Update etaH at the end of the time step :
C     Incorporate the Explicit part of -Divergence(Barotropic_Flow)

        IF (implicDiv2Dflow.EQ. 1. _d 0) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
          ENDDO
         ENDDO

        ELSE
         DO j=1,sNy
          DO i=1,sNx
            etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
     &       + (1. - implicDiv2Dflow)*dEtaHdt(i,j,bi,bj)
     &                               *deltaTfreesurf
          ENDDO
         ENDDO
        ENDIF

#ifdef ALLOW_OBCS
C--    Apply OBC to etaH (NonLin-FreeSurf): needed since viscous terms
C       depend on hFacZ which is not only function of boundary hFac values.
       IF ( useOBCS.AND.nonlinFreeSurf.GT.0 )
     &    CALL OBCS_APPLY_ETA( bi, bj, etaH, myThid )
#endif /* ALLOW_OBCS */

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

C- end bi,bj loop.
       ENDDO
      ENDDO

      IF ( implicDiv2Dflow .NE. 1. _d 0 .OR.
     &    ( useOBCS.AND.nonlinFreeSurf.GT.0 ) )
     &    CALL EXCH_XY_RL( etaH, myThid )

c     IF (useRealFreshWaterFlux .AND. myTime.EQ.startTime)
c    &    _EXCH_XY_RL( PmEpR, myThid )

#ifdef NONLIN_FRSURF
# ifndef DISABLE_SIGMA_CODE
      IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
C-     2nd bi,bj loop :

C-- copy etaHX -> dEtaXdt
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            dEtaWdt(i,j,bi,bj) = etaHw(i,j,bi,bj)
            dEtaSdt(i,j,bi,bj) = etaHs(i,j,bi,bj)
          ENDDO
         ENDDO

         DO j=1,sNy+1
          DO i=1,sNx+1
            etaHw(i,j,bi,bj)   = ( etaH (i-1,j,bi,bj)
     &                           + etaH ( i ,j,bi,bj) )*0.5 _d 0
            etaHs(i,j,bi,bj)   = ( etaH (i,j-1,bi,bj)
     &                           + etaH (i, j ,bi,bj) )*0.5 _d 0
c           etaHw(i,j,bi,bj)   = 0.5 _d 0
c    &                         *(   etaH (i-1,j,bi,bj)*rA(i-1,j,bi,bj)
c    &                            + etaH ( i ,j,bi,bj)*rA( i ,j,bi,bj)
c    &                          )*recip_rAw(i,j,bi,bj)
c           etaHs(i,j,bi,bj)   = 0.5 _d 0
c    &                         *(   etaH (i,j-1,bi,bj)*rA(i,j-1,bi,bj)
c    &                            + etaH (i, j ,bi,bj)*rA(i, j ,bi,bj)
c    &                          )*recip_rAs(i,j,bi,bj)
          ENDDO
         ENDDO

C- end 2nd bi,bj loop.
        ENDDO
       ENDDO

       CALL EXCH_UV_XY_RL( etaHw, etaHs, .FALSE., myThid )
       CALL EXCH_XY_RL( dEtaHdt, myThid )

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
C-     3rd bi,bj loop :

         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            dEtaWdt(i,j,bi,bj) = ( etaHw(i,j,bi,bj)
     &                           - dEtaWdt(i,j,bi,bj) )/deltaTfreesurf
            dEtaSdt(i,j,bi,bj) = ( etaHs(i,j,bi,bj)
     &                           - dEtaSdt(i,j,bi,bj) )/deltaTfreesurf
          ENDDO
         ENDDO

C- end 3rd bi,bj loop.
        ENDDO
       ENDDO

      ENDIF
# endif /* DISABLE_SIGMA_CODE */
#endif /* NONLIN_FRSURF */

#endif /* EXACT_CONSERV */

      RETURN
      END
1	jmc	1.12	C $Header: /u/gcmpack/MITgcm/model/src/update_etah.F,v 1.11 2011/04/27 22:20:22 jmc Exp $
2	jmc	1.1	C $Name: $
3
4	edhill	1.3	#include "PACKAGES_CONFIG.h"
5	jmc	1.1	#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: UPDATE_ETAH
9			C !INTERFACE:
10			SUBROUTINE UPDATE_ETAH( myTime, myIter, myThid )
11			C !DESCRIPTION: \bv
12			C ==========================================================
13	jmc	1.10	C \| SUBROUTINE UPDATE_ETAH
14	jmc	1.12	C \| o Update etaH after mom-correction-step/integr_continuity
15	jmc	1.1	C \| (required with NLFS to derive surface layer thickness)
16	jmc	1.12	C \| o Also derive SSH tendency at grid-cell Western and
17			C \| Southern edges (for hybrid sigma-coordinate)
18	jmc	1.1	C ==========================================================
19			C \ev
20
21			C !USES:
22			IMPLICIT NONE
23			C == Global variables
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "DYNVARS.h"
28			#include "GRID.h"
29			#include "SURFACE.h"
30			#include "FFIELDS.h"
31
32			C !INPUT/OUTPUT PARAMETERS:
33			C == Routine arguments ==
34			C myTime :: Current time in simulation
35			C myIter :: Current iteration number in simulation
36			C myThid :: Thread number for this instance of the routine.
37			_RL myTime
38			INTEGER myIter
39			INTEGER myThid
40
41			C !LOCAL VARIABLES:
42			#ifdef EXACT_CONSERV
43			C Local variables in common block
44
45			C Local variables
46	jmc	1.7	C i,j,bi,bj :: Loop counters
47			INTEGER i,j,bi,bj
48	jmc	1.1	CEOP
49
50
51			DO bj=myByLo(myThid),myByHi(myThid)
52	jmc	1.10	DO bi=myBxLo(myThid),myBxHi(myThid)
53	jmc	1.1
54			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
55
56	jmc	1.5	C-- before updating etaH, save current etaH field in etaHnm1
57			DO j=1-Oly,sNy+Oly
58			DO i=1-Olx,sNx+Olx
59	jmc	1.10	etaHnm1(i,j,bi,bj) = etaH(i,j,bi,bj)
60	jmc	1.5	ENDDO
61			ENDDO
62
63	jmc	1.10	C-- Update etaH at the end of the time step :
64	jmc	1.1	C Incorporate the Explicit part of -Divergence(Barotropic_Flow)
65
66	jmc	1.5	IF (implicDiv2Dflow.EQ. 1. _d 0) THEN
67			DO j=1-Oly,sNy+Oly
68			DO i=1-Olx,sNx+Olx
69	jmc	1.10	etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
70	jmc	1.5	ENDDO
71	jmc	1.1	ENDDO
72
73	jmc	1.5	ELSE
74			DO j=1,sNy
75			DO i=1,sNx
76	jmc	1.10	etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
77	jmc	1.5	& + (1. - implicDiv2Dflow)*dEtaHdt(i,j,bi,bj)
78			& *deltaTfreesurf
79			ENDDO
80	jmc	1.1	ENDDO
81	jmc	1.5	ENDIF
82	jmc	1.1
83			#ifdef ALLOW_OBCS
84	jmc	1.12	C-- Apply OBC to etaH (NonLin-FreeSurf): needed since viscous terms
85			C depend on hFacZ which is not only function of boundary hFac values.
86			IF ( useOBCS.AND.nonlinFreeSurf.GT.0 )
87			& CALL OBCS_APPLY_ETA( bi, bj, etaH, myThid )
88	jmc	1.1	#endif /* ALLOW_OBCS */
89
90			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
91
92			C- end bi,bj loop.
93			ENDDO
94	jmc	1.10	ENDDO
95	jmc	1.1
96	jmc	1.12	IF ( implicDiv2Dflow .NE. 1. _d 0 .OR.
97			& ( useOBCS.AND.nonlinFreeSurf.GT.0 ) )
98	jmc	1.10	& CALL EXCH_XY_RL( etaH, myThid )
99	jmc	1.1
100	jmc	1.6	c IF (useRealFreshWaterFlux .AND. myTime.EQ.startTime)
101	jmc	1.11	c & _EXCH_XY_RL( PmEpR, myThid )
102	jmc	1.2
103	jmc	1.10	#ifdef NONLIN_FRSURF
104			# ifndef DISABLE_SIGMA_CODE
105			IF ( nonlinFreeSurf.GT.0 .AND. selectSigmaCoord.NE.0 ) THEN
106
107			DO bj=myByLo(myThid),myByHi(myThid)
108			DO bi=myBxLo(myThid),myBxHi(myThid)
109			C- 2nd bi,bj loop :
110
111			C-- copy etaHX -> dEtaXdt
112			DO j=1-Oly,sNy+Oly
113			DO i=1-Olx,sNx+Olx
114			dEtaWdt(i,j,bi,bj) = etaHw(i,j,bi,bj)
115			dEtaSdt(i,j,bi,bj) = etaHs(i,j,bi,bj)
116			ENDDO
117			ENDDO
118
119			DO j=1,sNy+1
120			DO i=1,sNx+1
121			etaHw(i,j,bi,bj) = ( etaH (i-1,j,bi,bj)
122			& + etaH ( i ,j,bi,bj) )*0.5 _d 0
123			etaHs(i,j,bi,bj) = ( etaH (i,j-1,bi,bj)
124			& + etaH (i, j ,bi,bj) )*0.5 _d 0
125			c etaHw(i,j,bi,bj) = 0.5 _d 0
126			c & ( etaH (i-1,j,bi,bj)rA(i-1,j,bi,bj)
127			c & + etaH ( i ,j,bi,bj)*rA( i ,j,bi,bj)
128			c & )*recip_rAw(i,j,bi,bj)
129			c etaHs(i,j,bi,bj) = 0.5 _d 0
130			c & ( etaH (i,j-1,bi,bj)rA(i,j-1,bi,bj)
131			c & + etaH (i, j ,bi,bj)*rA(i, j ,bi,bj)
132			c & )*recip_rAs(i,j,bi,bj)
133			ENDDO
134			ENDDO
135
136			C- end 2nd bi,bj loop.
137			ENDDO
138			ENDDO
139
140			CALL EXCH_UV_XY_RL( etaHw, etaHs, .FALSE., myThid )
141			CALL EXCH_XY_RL( dEtaHdt, myThid )
142
143			DO bj=myByLo(myThid),myByHi(myThid)
144			DO bi=myBxLo(myThid),myBxHi(myThid)
145			C- 3rd bi,bj loop :
146
147			DO j=1-Oly,sNy+Oly
148			DO i=1-Olx,sNx+Olx
149			dEtaWdt(i,j,bi,bj) = ( etaHw(i,j,bi,bj)
150			& - dEtaWdt(i,j,bi,bj) )/deltaTfreesurf
151			dEtaSdt(i,j,bi,bj) = ( etaHs(i,j,bi,bj)
152			& - dEtaSdt(i,j,bi,bj) )/deltaTfreesurf
153			ENDDO
154			ENDDO
155
156			C- end 3rd bi,bj loop.
157			ENDDO
158			ENDDO
159
160			ENDIF
161			# endif /* DISABLE_SIGMA_CODE */
162			#endif /* NONLIN_FRSURF */
163
164	jmc	1.4	#endif /* EXACT_CONSERV */
165	jmc	1.1
166			RETURN
167			END