/[MITgcm]/MITgcm/model/src/calc_gs.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_gs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.7 - (hide annotations) (download)
Wed Jun 10 16:05:39 1998 UTC (25 years, 11 months ago) by adcroft
Branch: MAIN
Changes since 1.6: +98 -25 lines
Added code to bring "salt" up-to-date with "theta".
One caveat is that implicit diffusion of salt is done with the
diffusivity of theta. We'll sort this out later. In explicit
mode, diffKzS is used.

1 adcroft 1.7 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.10 1998/06/08 21:43:00 cnh Exp $
2 cnh 1.1
3     #include "CPP_EEOPTIONS.h"
4    
5     CStartOfInterFace
6     SUBROUTINE CALC_GS(
7     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
8     I xA,yA,uTrans,vTrans,wTrans,maskup,
9 adcroft 1.7 I K13,K23,KappaZS,KapGM,
10     U af,df,fZon,fMer,fVerS,
11 cnh 1.1 I myThid )
12     C /==========================================================\
13     C | SUBROUTINE CALC_GS |
14 adcroft 1.7 C | o Calculate the salt tendency terms. |
15 cnh 1.1 C |==========================================================|
16     C | A procedure called EXTERNAL_FORCING_S is called from |
17     C | here. These procedures can be used to add per problem |
18 adcroft 1.7 C | E-P flux source terms. |
19 cnh 1.1 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 adcroft 1.7 C fVerS - Flux of salt (S) in the vertical
52 cnh 1.1 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 adcroft 1.7 C myThid - Instance number for this innvocation of CALC_GT
64 cnh 1.1 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL fVerS (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.7 _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 KappaZS(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.7 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.7 INTEGER i,j
87     LOGICAL TOP_LAYER
88     _RL afFacS, dfFacS
89     _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90     _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 cnh 1.1
92     afFacS = 1. _d 0
93     dfFacS = 1. _d 0
94 adcroft 1.7 TOP_LAYER = K .EQ. 1
95 cnh 1.1
96     C--- Calculate advective and diffusive fluxes between cells.
97    
98     C-- Zonal flux (fZon is at west face of "salt" cell)
99     C Advective component of zonal flux
100     DO j=jMin,jMax
101     DO i=iMin,iMax
102     af(i,j) =
103     & uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0
104     ENDDO
105     ENDDO
106 adcroft 1.7 C Zonal tracer gradient
107     DO j=jMin,jMax
108     DO i=iMin,iMax
109     dSdx(i,j) = _rdxC(i,j,bi,bj)*
110     & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
111     ENDDO
112     ENDDO
113 cnh 1.1 C Diffusive component of zonal flux
114     DO j=jMin,jMax
115     DO i=iMin,iMax
116 adcroft 1.7 df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*
117     & xA(i,j)*dSdx(i,j)
118 cnh 1.1 ENDDO
119     ENDDO
120     C Net zonal flux
121     DO j=jMin,jMax
122     DO i=iMin,iMax
123     fZon(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
124     ENDDO
125     ENDDO
126    
127     C-- Meridional flux (fMer is at south face of "salt" cell)
128     C Advective component of meridional flux
129     DO j=jMin,jMax
130     DO i=iMin,iMax
131     C Advective component of meridional flux
132     af(i,j) =
133     & vTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*0.5 _d 0
134     ENDDO
135     ENDDO
136 adcroft 1.7 C Zonal tracer gradient
137     DO j=jMin,jMax
138     DO i=iMin,iMax
139     dSdy(i,j) = _rdyC(i,j,bi,bj)*
140     & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
141     ENDDO
142     ENDDO
143 cnh 1.1 C Diffusive component of meridional flux
144     DO j=jMin,jMax
145     DO i=iMin,iMax
146 adcroft 1.7 df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i,j-1)))*
147     & yA(i,j)*dSdy(i,j)
148 cnh 1.1 ENDDO
149     ENDDO
150     C Net meridional flux
151     DO j=jMin,jMax
152     DO i=iMin,iMax
153     fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
154     ENDDO
155     ENDDO
156    
157 adcroft 1.7 C-- Interpolate terms for Redi/GM scheme
158     DO j=jMin,jMax
159     DO i=iMin,iMax
160     dSdx(i,j) = 0.5*(
161     & +0.5*(_maskW(i+1,j,k,bi,bj)*_rdxC(i+1,j,bi,bj)*
162     & (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))
163     & +_maskW(i,j,k,bi,bj)*_rdxC(i,j,bi,bj)*
164     & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))
165     & +0.5*(_maskW(i+1,j,km1,bi,bj)*_rdxC(i+1,j,bi,bj)*
166     & (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))
167     & +_maskW(i,j,km1,bi,bj)*_rdxC(i,j,bi,bj)*
168     & (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))
169     & )
170     ENDDO
171     ENDDO
172     DO j=jMin,jMax
173     DO i=iMin,iMax
174     dSdy(i,j) = 0.5*(
175     & +0.5*(_maskS(i,j,k,bi,bj)*_rdyC(i,j,bi,bj)*
176     & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
177     & +_maskS(i,j+1,k,bi,bj)*_rdyC(i,j+1,bi,bj)*
178     & (salt(i,j+1,k,bi,bj)-salt(i,j,k,bi,bj)))
179     & +0.5*(_maskS(i,j,km1,bi,bj)*_rdyC(i,j,bi,bj)*
180     & (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))
181     & +_maskS(i,j+1,km1,bi,bj)*_rdyC(i,j+1,bi,bj)*
182     & (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))
183     & )
184     ENDDO
185     ENDDO
186    
187 cnh 1.1 C-- Vertical flux (fVerS) above
188     C Advective component of vertical flux
189 adcroft 1.7 C Note: For K=1 then KM1=1 this gives a barZ(T) = T
190     C (this plays the role of the free-surface correction)
191 cnh 1.1 DO j=jMin,jMax
192     DO i=iMin,iMax
193     af(i,j) =
194     & wTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
195     ENDDO
196     ENDDO
197     C Diffusive component of vertical flux
198 adcroft 1.7 C Note: For K=1 then KM1=1 this gives a dS/dz = 0 upper
199     C boundary condition.
200 cnh 1.1 DO j=jMin,jMax
201     DO i=iMin,iMax
202 adcroft 1.7 df(i,j) = _zA(i,j,bi,bj)*(
203     & -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)
204     & -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)
205     & )
206 cnh 1.1 ENDDO
207     ENDDO
208 adcroft 1.7 IF (.NOT.implicitDiffusion) THEN
209     DO j=jMin,jMax
210     DO i=iMin,iMax
211     df(i,j) = df(i,j) + _zA(i,j,bi,bj)*(
212     & -KappaZS(i,j,k)*rdzC(k)
213     & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))
214     & )
215     ENDDO
216     ENDDO
217     ENDIF
218 cnh 1.1 C Net vertical flux
219     DO j=jMin,jMax
220     DO i=iMin,iMax
221 adcroft 1.7 fVerS(i,j,kUp) = ( afFacS*af(i,j)+ dfFacS*df(i,j) )*maskUp(i,j)
222 cnh 1.1 ENDDO
223     ENDDO
224 adcroft 1.7 IF ( TOP_LAYER ) THEN
225     DO j=jMin,jMax
226     DO i=iMin,iMax
227     fVerS(i,j,kUp) = afFacS*af(i,j)*freeSurfFac
228     ENDDO
229     ENDDO
230     ENDIF
231 cnh 1.1
232     C-- Tendency is minus divergence of the fluxes.
233     C Note. Tendency terms will only be correct for range
234     C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
235     C will contain valid floating point numbers but
236     C they are not algorithmically correct. These points
237     C are not used.
238     DO j=jMin,jMax
239     DO i=iMin,iMax
240 adcroft 1.7 C & -_rhFacC(i,j,k,bi,bj)*rdzF(k)*_rdxF(i,j,bi,bj)*_rdyF(i,j,bi,bj)
241     C & -_rhFacC(i,j,k,bi,bj)*rdzF(k)/_zA(i,j,bi,bj)
242     C #define _rVolS(i,j,k,bi,bj) _rhFacC(i,j,k,bi,bj)*rdzF(k)*_rdxF(i,j,bi,bj)*_rdyF(i,j,bi,bj)
243     #define _rVolS(i,j,k,bi,bj) _rhFacC(i,j,k,bi,bj)*rdzF(k)/_zA(i,j,bi,bj)
244 cnh 1.1 gS(i,j,k,bi,bj)=
245 adcroft 1.7 & -_rVolS(i,j,k,bi,bj)
246 cnh 1.1 & *(
247     & +( fZon(i+1,j)-fZon(i,j) )
248     & +( fMer(i,j+1)-fMer(i,j) )
249     & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )
250     & )
251     ENDDO
252     ENDDO
253    
254 adcroft 1.7 C-- External P-E forcing term(s)
255 cnh 1.1
256     RETURN
257     END

  ViewVC Help
Powered by ViewVC 1.1.22