model/src/calc_gt.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.4 1998/05/25 20:05:55 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 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

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

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	cnh	1.5	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.4 1998/05/25 20:05:55 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 dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89			_RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90	cnh	1.1
91			afFacT = 1. _d 0
92			dfFacT = 1. _d 0
93
94			C--- Calculate advective and diffusive fluxes between cells.
95
96			C-- Zonal flux (fZon is at west face of "theta" cell)
97			C Advective component of zonal flux
98			DO j=jMin,jMax
99			DO i=iMin,iMax
100			af(i,j) =
101			& uTrans(i,j)(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))0.5 _d 0
102			ENDDO
103			ENDDO
104	adcroft	1.3	C Zonal tracer gradient
105			DO j=jMin,jMax
106			DO i=iMin,iMax
107	cnh	1.5	dTdx(i,j) = _rdxC(i,j,bi,bj)*
108	adcroft	1.3	& (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
109			ENDDO
110			ENDDO
111	cnh	1.1	C Diffusive component of zonal flux
112			DO j=jMin,jMax
113			DO i=iMin,iMax
114	adcroft	1.3	df(i,j) = -(diffKhT+0.5(KapGM(i,j)+KapGM(i-1,j)))
115			& xA(i,j)*dTdx(i,j)
116	cnh	1.1	ENDDO
117			ENDDO
118			C Net zonal flux
119			DO j=jMin,jMax
120			DO i=iMin,iMax
121			fZon(i,j) = afFacTaf(i,j) + dfFacTdf(i,j)
122			ENDDO
123			ENDDO
124
125			C-- Meridional flux (fMer is at south face of "theta" cell)
126			C Advective component of meridional flux
127			DO j=jMin,jMax
128			DO i=iMin,iMax
129			C Advective component of meridional flux
130			af(i,j) =
131			& vTrans(i,j)(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))0.5 _d 0
132			ENDDO
133			ENDDO
134	adcroft	1.3	C Zonal tracer gradient
135			DO j=jMin,jMax
136			DO i=iMin,iMax
137			dTdy(i,j) = rdyC(i,j,bi,bj)*
138			& (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
139			ENDDO
140			ENDDO
141	cnh	1.1	C Diffusive component of meridional flux
142			DO j=jMin,jMax
143			DO i=iMin,iMax
144	adcroft	1.3	df(i,j) = -(diffKhT+0.5(KapGM(i,j)+KapGM(i,j-1)))
145			& yA(i,j)*dTdy(i,j)
146	cnh	1.1	ENDDO
147			ENDDO
148			C Net meridional flux
149			DO j=jMin,jMax
150			DO i=iMin,iMax
151			fMer(i,j) = afFacTaf(i,j) + dfFacTdf(i,j)
152			ENDDO
153			ENDDO
154
155	adcroft	1.3	C-- Interpolate terms for Redi/GM scheme
156			DO j=jMin,jMax
157			DO i=iMin,iMax
158			dTdx(i,j) = 0.5*(
159	cnh	1.5	& +0.5(maskW(i+1,j,k,bi,bj)_rdxC(i+1,j,bi,bj)*
160	adcroft	1.3	& (theta(i+1,j,k,bi,bj)-theta(i,j,k,bi,bj))
161	cnh	1.5	& +maskW(i,j,k,bi,bj)_rdxC(i,j,bi,bj)
162	adcroft	1.3	& (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj)))
163	cnh	1.5	& +0.5(maskW(i+1,j,km1,bi,bj)_rdxC(i+1,j,bi,bj)*
164	adcroft	1.3	& (theta(i+1,j,km1,bi,bj)-theta(i,j,km1,bi,bj))
165	cnh	1.5	& +maskW(i,j,km1,bi,bj)_rdxC(i,j,bi,bj)
166	adcroft	1.3	& (theta(i,j,km1,bi,bj)-theta(i-1,j,km1,bi,bj)))
167			& )
168			ENDDO
169			ENDDO
170			DO j=jMin,jMax
171			DO i=iMin,iMax
172			dTdy(i,j) = 0.5*(
173			& +0.5(maskS(i,j,k,bi,bj)rdyC(i,j,bi,bj)*
174			& (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
175			& +maskS(i,j+1,k,bi,bj)rdyC(i,j+1,bi,bj)
176			& (theta(i,j+1,k,bi,bj)-theta(i,j,k,bi,bj)))
177			& +0.5(maskS(i,j,km1,bi,bj)rdyC(i,j,bi,bj)*
178			& (theta(i,j,km1,bi,bj)-theta(i,j-1,km1,bi,bj))
179			& +maskS(i,j+1,km1,bi,bj)rdyC(i,j+1,bi,bj)
180			& (theta(i,j+1,km1,bi,bj)-theta(i,j,km1,bi,bj)))
181			& )
182			ENDDO
183			ENDDO
184
185	cnh	1.1	C-- Vertical flux (fVerT) above
186			C Advective component of vertical flux
187	adcroft	1.3	C Note: For K=1 then KM1=1 this gives a barZ(T) = T
188			C (this plays the role of the free-surface correction)
189	cnh	1.1	DO j=jMin,jMax
190			DO i=iMin,iMax
191			af(i,j) =
192			& wTrans(i,j)(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))0.5 _d 0
193			ENDDO
194			ENDDO
195			C Diffusive component of vertical flux
196	adcroft	1.3	C Note: For K=1 then KM1=1 this gives a dT/dz = 0 upper
197			C boundary condition.
198	cnh	1.1	DO j=jMin,jMax
199			DO i=iMin,iMax
200	adcroft	1.3	df(i,j) = zA(i,j,bi,bj)*(
201			& -(diffKzT+KapGM(i,j)K33(i,j,k))rdzC(k)
202	cnh	1.1	& *(theta(i,j,kM1,bi,bj)-theta(i,j,k,bi,bj))
203	adcroft	1.3	& -KapGM(i,j)K13(i,j,k)dTdx(i,j)
204			& -KapGM(i,j)K23(i,j,k)dTdy(i,j)
205			& )
206	cnh	1.1	ENDDO
207			ENDDO
208			C Net vertical flux
209			DO j=jMin,jMax
210			DO i=iMin,iMax
211			fVerT(i,j,kUp) = (afFacTaf(i,j) + dfFacTdf(i,j))*maskUp(i,j)
212			ENDDO
213			ENDDO
214
215			C-- Tendency is minus divergence of the fluxes.
216			C Note. Tendency terms will only be correct for range
217			C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
218			C will contain valid floating point numbers but
219			C they are not algorithmically correct. These points
220			C are not used.
221			DO j=jMin,jMax
222			DO i=iMin,iMax
223			gT(i,j,k,bi,bj)=
224			& -rHFacC(i,j,k,bi,bj)rdzF(k)rDxF(i,j,bi,bj)*rDyF(i,j,bi,bj)
225			& *(
226			& +( fZon(i+1,j)-fZon(i,j) )
227			& +( fMer(i,j+1)-fMer(i,j) )
228			& +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )
229			& )
230			ENDDO
231			ENDDO
232
233			C-- External thermal forcing term(s)
234
235			RETURN
236			END