model/src/calc_gt.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.2 1998/04/24 02:05:40 cnh Exp $

#include "CPP_EEOPTIONS.h"

CStartOfInterFace
      SUBROUTINE CALC_GT( 
     I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
     I           xA,yA,uTrans,vTrans,wTrans,maskup,
     I           K13,K23,K33,KapGM,
     U           af,df,fZon,fMer,fVerT,
     I           myThid )
C     /==========================================================\
C     | SUBROUTINE CALC_GT                                       |
C     | o Calculate the temperature tendency terms.              |
C     |==========================================================|
C     | A procedure called EXTERNAL_FORCING_T is called from     |
C     | here. These procedures can be used to add per problem    |
C     | heat flux source terms.                                  |
C     | Note: Although it is slightly counter-intuitive the      |
C     |       EXTERNAL_FORCING routine is not the place to put   |
C     |       file I/O. Instead files that are required to       |
C     |       calculate the external source terms are generally  |
C     |       read during the model main loop. This makes the    |
C     |       logisitics of multi-processing simpler and also    |
C     |       makes the adjoint generation simpler. It also      |
C     |       allows for I/O to overlap computation where that   |
C     |       is supported by hardware.                          |
C     | Aside from the problem specific term the code here       |
C     | forms the tendency terms due to advection and mixing     |
C     | The baseline implementation here uses a centered         |
C     | difference form for the advection term and a tensorial   |
C     | divergence of a flux form for the diffusive term. The    |
C     | diffusive term is formulated so that isopycnal mixing and|
C     | GM-style subgrid-scale terms can be incorporated b simply|
C     | setting the diffusion tensor terms appropriately.        |
C     \==========================================================/
      IMPLICIT NONE

C     == GLobal variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     == Routine arguments ==
C     fZon    - Work array for flux of temperature in the east-west 
C               direction at the west face of a cell.
C     fMer    - Work array for flux of temperature in the north-south 
C               direction at the south face of a cell.
C     fVerT   - Flux of temperature (T) in the vertical 
C               direction at the upper(U) and lower(D) faces of a cell.
C     maskUp  - Land mask used to denote base of the domain.
C     xA      - Tracer cell face area normal to X
C     yA      - Tracer cell face area normal to X
C     uTrans  - Zonal volume transport through cell face
C     vTrans  - Meridional volume transport through cell face
C     wTrans  - Vertical volume transport through cell face
C     af      - Advective flux component work array
C     df      - Diffusive flux component work array
C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
C                                      results will be set.
C     myThid - Instance number for this innvocation of CALC_GT
      _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
      _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
      _RL K33   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
      _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER k,kUp,kDown,kM1
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER myThid
CEndOfInterface

C     == Local variables ==
C     I, J, K - Loop counters
      INTEGER i,j
      _RL afFacT, dfFacT
      _RL dutdxFac
      _RL dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

      afFacT = 1. _d 0
      dfFacT = 1. _d 0
      dutdxFac = afFacT

C---
C---  Calculate advective and diffusive fluxes between cells.
C---

C--   Zonal flux (fZon is at west face of "theta" cell)
C     Advective component of zonal flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        af(i,j) = 
     &   uTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*0.5 _d 0
       ENDDO
      ENDDO
C     Zonal tracer gradient
      DO j=jMin,jMax
       DO i=iMin,iMax
        dTdx(i,j) = rdxC(i,j,bi,bj)*
     &  (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
       ENDDO
      ENDDO
C     Diffusive component of zonal flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i-1,j)))*
     &            xA(i,j)*dTdx(i,j)
       ENDDO
      ENDDO
C     Net zonal flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        fZon(i,j) = afFacT*af(i,j) + dfFacT*df(i,j)
       ENDDO
      ENDDO

C--   Meridional flux (fMer is at south face of "theta" cell)
C     Advective component of meridional flux
      DO j=jMin,jMax
       DO i=iMin,iMax
C       Advective component of meridional flux
        af(i,j) =
     &   vTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*0.5 _d 0
       ENDDO
      ENDDO
C     Zonal tracer gradient
      DO j=jMin,jMax
       DO i=iMin,iMax
        dTdy(i,j) = rdyC(i,j,bi,bj)*
     &  (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
       ENDDO
      ENDDO
C     Diffusive component of meridional flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i,j-1)))*
     &            yA(i,j)*dTdy(i,j)
       ENDDO
      ENDDO
C     Net meridional flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        fMer(i,j) = afFacT*af(i,j) + dfFacT*df(i,j)
       ENDDO
      ENDDO

C--   Interpolate terms for Redi/GM scheme
      DO j=jMin,jMax
       DO i=iMin,iMax
        dTdx(i,j) = 0.5*( 
     &   +0.5*(maskW(i+1,j,k,bi,bj)*rdxC(i+1,j,bi,bj)*
     &           (theta(i+1,j,k,bi,bj)-theta(i,j,k,bi,bj))
     &        +maskW(i,j,k,bi,bj)*rdxC(i,j,bi,bj)*
     &           (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj)))
     &   +0.5*(maskW(i+1,j,km1,bi,bj)*rdxC(i+1,j,bi,bj)*
     &           (theta(i+1,j,km1,bi,bj)-theta(i,j,km1,bi,bj))
     &        +maskW(i,j,km1,bi,bj)*rdxC(i,j,bi,bj)*
     &           (theta(i,j,km1,bi,bj)-theta(i-1,j,km1,bi,bj)))
     &       )
       ENDDO
      ENDDO
      DO j=jMin,jMax
       DO i=iMin,iMax
        dTdy(i,j) = 0.5*(
     &   +0.5*(maskS(i,j,k,bi,bj)*rdyC(i,j,bi,bj)*
     &           (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
     &        +maskS(i,j+1,k,bi,bj)*rdyC(i,j+1,bi,bj)*
     &           (theta(i,j+1,k,bi,bj)-theta(i,j,k,bi,bj)))
     &   +0.5*(maskS(i,j,km1,bi,bj)*rdyC(i,j,bi,bj)*
     &           (theta(i,j,km1,bi,bj)-theta(i,j-1,km1,bi,bj))
     &        +maskS(i,j+1,km1,bi,bj)*rdyC(i,j+1,bi,bj)*
     &           (theta(i,j+1,km1,bi,bj)-theta(i,j,km1,bi,bj)))
     &       )
       ENDDO
      ENDDO

C--   Vertical flux (fVerT) above
C     Advective component of vertical flux
C     Note: For K=1 then KM1=1 this gives a barZ(T) = T
C     (this plays the role of the free-surface correction)
      DO j=jMin,jMax
       DO i=iMin,iMax
        af(i,j) =
     &   wTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0
       ENDDO
      ENDDO
C     Diffusive component of vertical flux
C     Note: For K=1 then KM1=1 this gives a dT/dz = 0 upper
C           boundary condition.
      DO j=jMin,jMax
       DO i=iMin,iMax
        df(i,j) = zA(i,j,bi,bj)*(
     &   -(diffKzT+KapGM(i,j)*K33(i,j,k))*rdzC(k)
     &   *(theta(i,j,kM1,bi,bj)-theta(i,j,k,bi,bj))
     &   -KapGM(i,j)*K13(i,j,k)*dTdx(i,j)
     &   -KapGM(i,j)*K23(i,j,k)*dTdy(i,j)
     &   )
       ENDDO
      ENDDO
C     Net vertical flux
      DO j=jMin,jMax
       DO i=iMin,iMax
        fVerT(i,j,kUp) = (afFacT*af(i,j) + dfFacT*df(i,j))*maskUp(i,j)
       ENDDO
      ENDDO

C--   Tendency is minus divergence of the fluxes.
C     Note. Tendency terms will only be correct for range
C           i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
C           will contain valid floating point numbers but 
C           they are not algorithmically correct. These points
C           are not used.
      DO j=jMin,jMax
       DO i=iMin,iMax
        gT(i,j,k,bi,bj)=
     &   -rHFacC(i,j,k,bi,bj)*rdzF(k)*rDxF(i,j,bi,bj)*rDyF(i,j,bi,bj)
     &   *(
     &    +( fZon(i+1,j)-fZon(i,j) )
     &    +( fMer(i,j+1)-fMer(i,j) )
     &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )
     &    )
       ENDDO
      ENDDO

C--   External thermal forcing term(s)

      RETURN
      END
1	adcroft	1.3	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.2 1998/04/24 02:05:40 cnh Exp $
2	cnh	1.1
3			#include "CPP_EEOPTIONS.h"
4
5			CStartOfInterFace
6			SUBROUTINE CALC_GT(
7			I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
8			I xA,yA,uTrans,vTrans,wTrans,maskup,
9	adcroft	1.3	I K13,K23,K33,KapGM,
10	cnh	1.1	U af,df,fZon,fMer,fVerT,
11			I myThid )
12			C /==========================================================\
13			C \| SUBROUTINE CALC_GT \|
14			C \| o Calculate the temperature tendency terms. \|
15			C \|==========================================================\|
16			C \| A procedure called EXTERNAL_FORCING_T is called from \|
17			C \| here. These procedures can be used to add per problem \|
18			C \| heat flux source terms. \|
19			C \| Note: Although it is slightly counter-intuitive the \|
20			C \| EXTERNAL_FORCING routine is not the place to put \|
21			C \| file I/O. Instead files that are required to \|
22			C \| calculate the external source terms are generally \|
23			C \| read during the model main loop. This makes the \|
24			C \| logisitics of multi-processing simpler and also \|
25			C \| makes the adjoint generation simpler. It also \|
26			C \| allows for I/O to overlap computation where that \|
27			C \| is supported by hardware. \|
28			C \| Aside from the problem specific term the code here \|
29			C \| forms the tendency terms due to advection and mixing \|
30			C \| The baseline implementation here uses a centered \|
31			C \| difference form for the advection term and a tensorial \|
32			C \| divergence of a flux form for the diffusive term. The \|
33			C \| diffusive term is formulated so that isopycnal mixing and\|
34			C \| GM-style subgrid-scale terms can be incorporated b simply\|
35			C \| setting the diffusion tensor terms appropriately. \|
36			C \==========================================================/
37			IMPLICIT NONE
38
39			C == GLobal variables ==
40			#include "SIZE.h"
41			#include "DYNVARS.h"
42			#include "EEPARAMS.h"
43			#include "PARAMS.h"
44			#include "GRID.h"
45
46			C == Routine arguments ==
47			C fZon - Work array for flux of temperature in the east-west
48			C direction at the west face of a cell.
49			C fMer - Work array for flux of temperature in the north-south
50			C direction at the south face of a cell.
51			C fVerT - Flux of temperature (T) in the vertical
52			C direction at the upper(U) and lower(D) faces of a cell.
53			C maskUp - Land mask used to denote base of the domain.
54			C xA - Tracer cell face area normal to X
55			C yA - Tracer cell face area normal to X
56			C uTrans - Zonal volume transport through cell face
57			C vTrans - Meridional volume transport through cell face
58			C wTrans - Vertical volume transport through cell face
59			C af - Advective flux component work array
60			C df - Diffusive flux component work array
61			C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
62			C results will be set.
63			C myThid - Instance number for this innvocation of CALC_GT
64			_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65			_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66			_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
67			_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69			_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70			_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71			_RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72			_RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73	adcroft	1.3	_RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
74			_RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
75			_RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
76			_RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77	cnh	1.1	_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78			_RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79	adcroft	1.3	INTEGER k,kUp,kDown,kM1
80	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
81			INTEGER myThid
82			CEndOfInterface
83
84			C == Local variables ==
85			C I, J, K - Loop counters
86	adcroft	1.3	INTEGER i,j
87			_RL afFacT, dfFacT
88			_RL dutdxFac
89			_RL dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90			_RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91	cnh	1.1
92			afFacT = 1. _d 0
93			dfFacT = 1. _d 0
94			dutdxFac = afFacT
95
96			C---
97			C--- Calculate advective and diffusive fluxes between cells.
98			C---
99
100			C-- Zonal flux (fZon is at west face of "theta" cell)
101			C Advective component of zonal flux
102			DO j=jMin,jMax
103			DO i=iMin,iMax
104			af(i,j) =
105			& uTrans(i,j)(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))0.5 _d 0
106			ENDDO
107			ENDDO
108	adcroft	1.3	C Zonal tracer gradient
109			DO j=jMin,jMax
110			DO i=iMin,iMax
111			dTdx(i,j) = rdxC(i,j,bi,bj)*
112			& (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
113			ENDDO
114			ENDDO
115	cnh	1.1	C Diffusive component of zonal flux
116			DO j=jMin,jMax
117			DO i=iMin,iMax
118	adcroft	1.3	df(i,j) = -(diffKhT+0.5(KapGM(i,j)+KapGM(i-1,j)))
119			& xA(i,j)*dTdx(i,j)
120	cnh	1.1	ENDDO
121			ENDDO
122			C Net zonal flux
123			DO j=jMin,jMax
124			DO i=iMin,iMax
125			fZon(i,j) = afFacTaf(i,j) + dfFacTdf(i,j)
126			ENDDO
127			ENDDO
128
129			C-- Meridional flux (fMer is at south face of "theta" cell)
130			C Advective component of meridional flux
131			DO j=jMin,jMax
132			DO i=iMin,iMax
133			C Advective component of meridional flux
134			af(i,j) =
135			& vTrans(i,j)(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))0.5 _d 0
136			ENDDO
137			ENDDO
138	adcroft	1.3	C Zonal tracer gradient
139			DO j=jMin,jMax
140			DO i=iMin,iMax
141			dTdy(i,j) = rdyC(i,j,bi,bj)*
142			& (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
143			ENDDO
144			ENDDO
145	cnh	1.1	C Diffusive component of meridional flux
146			DO j=jMin,jMax
147			DO i=iMin,iMax
148	adcroft	1.3	df(i,j) = -(diffKhT+0.5(KapGM(i,j)+KapGM(i,j-1)))
149			& yA(i,j)*dTdy(i,j)
150	cnh	1.1	ENDDO
151			ENDDO
152			C Net meridional flux
153			DO j=jMin,jMax
154			DO i=iMin,iMax
155			fMer(i,j) = afFacTaf(i,j) + dfFacTdf(i,j)
156			ENDDO
157			ENDDO
158
159	adcroft	1.3	C-- Interpolate terms for Redi/GM scheme
160			DO j=jMin,jMax
161			DO i=iMin,iMax
162			dTdx(i,j) = 0.5*(
163			& +0.5(maskW(i+1,j,k,bi,bj)rdxC(i+1,j,bi,bj)*
164			& (theta(i+1,j,k,bi,bj)-theta(i,j,k,bi,bj))
165			& +maskW(i,j,k,bi,bj)rdxC(i,j,bi,bj)
166			& (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj)))
167			& +0.5(maskW(i+1,j,km1,bi,bj)rdxC(i+1,j,bi,bj)*
168			& (theta(i+1,j,km1,bi,bj)-theta(i,j,km1,bi,bj))
169			& +maskW(i,j,km1,bi,bj)rdxC(i,j,bi,bj)
170			& (theta(i,j,km1,bi,bj)-theta(i-1,j,km1,bi,bj)))
171			& )
172			ENDDO
173			ENDDO
174			DO j=jMin,jMax
175			DO i=iMin,iMax
176			dTdy(i,j) = 0.5*(
177			& +0.5(maskS(i,j,k,bi,bj)rdyC(i,j,bi,bj)*
178			& (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
179			& +maskS(i,j+1,k,bi,bj)rdyC(i,j+1,bi,bj)
180			& (theta(i,j+1,k,bi,bj)-theta(i,j,k,bi,bj)))
181			& +0.5(maskS(i,j,km1,bi,bj)rdyC(i,j,bi,bj)*
182			& (theta(i,j,km1,bi,bj)-theta(i,j-1,km1,bi,bj))
183			& +maskS(i,j+1,km1,bi,bj)rdyC(i,j+1,bi,bj)
184			& (theta(i,j+1,km1,bi,bj)-theta(i,j,km1,bi,bj)))
185			& )
186			ENDDO
187			ENDDO
188
189	cnh	1.1	C-- Vertical flux (fVerT) above
190			C Advective component of vertical flux
191	adcroft	1.3	C Note: For K=1 then KM1=1 this gives a barZ(T) = T
192			C (this plays the role of the free-surface correction)
193	cnh	1.1	DO j=jMin,jMax
194			DO i=iMin,iMax
195			af(i,j) =
196			& wTrans(i,j)(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))0.5 _d 0
197			ENDDO
198			ENDDO
199			C Diffusive component of vertical flux
200	adcroft	1.3	C Note: For K=1 then KM1=1 this gives a dT/dz = 0 upper
201			C boundary condition.
202	cnh	1.1	DO j=jMin,jMax
203			DO i=iMin,iMax
204	adcroft	1.3	df(i,j) = zA(i,j,bi,bj)*(
205			& -(diffKzT+KapGM(i,j)K33(i,j,k))rdzC(k)
206	cnh	1.1	& *(theta(i,j,kM1,bi,bj)-theta(i,j,k,bi,bj))
207	adcroft	1.3	& -KapGM(i,j)K13(i,j,k)dTdx(i,j)
208			& -KapGM(i,j)K23(i,j,k)dTdy(i,j)
209			& )
210	cnh	1.1	ENDDO
211			ENDDO
212			C Net vertical flux
213			DO j=jMin,jMax
214			DO i=iMin,iMax
215			fVerT(i,j,kUp) = (afFacTaf(i,j) + dfFacTdf(i,j))*maskUp(i,j)
216			ENDDO
217			ENDDO
218
219			C-- Tendency is minus divergence of the fluxes.
220			C Note. Tendency terms will only be correct for range
221			C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
222			C will contain valid floating point numbers but
223			C they are not algorithmically correct. These points
224			C are not used.
225			DO j=jMin,jMax
226			DO i=iMin,iMax
227			gT(i,j,k,bi,bj)=
228			& -rHFacC(i,j,k,bi,bj)rdzF(k)rDxF(i,j,bi,bj)*rDyF(i,j,bi,bj)
229			& *(
230			& +( fZon(i+1,j)-fZon(i,j) )
231			& +( fMer(i,j+1)-fMer(i,j) )
232			& +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )
233			& )
234			ENDDO
235			ENDDO
236
237			C-- External thermal forcing term(s)
238
239			RETURN
240			END