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

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

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


Revision 1.7 - (hide annotations) (download)
Mon May 25 16:17:36 1998 UTC (26 years ago) by cnh
Branch: MAIN
Changes since 1.6: +59 -58 lines
Added changes to support implicit free-surface.
 - included gBaro a "barotropic" gravity that can
   be set differently to the g.rhoprime gravity.
 - discovered and fixed coding error in dynamics
   loop. Per tile temporaries that needed correct
   initial values were not being reset for each tile.

1 cnh 1.7 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.6 1998/05/20 21:29:31 adcroft Exp $
2 cnh 1.1
3     #include "CPP_EEOPTIONS.h"
4    
5     SUBROUTINE DYNAMICS(myThid)
6     C /==========================================================\
7     C | SUBROUTINE DYNAMICS |
8     C | o Controlling routine for the explicit part of the model |
9     C | dynamics. |
10     C |==========================================================|
11     C | This routine evaluates the "dynamics" terms for each |
12     C | block of ocean in turn. Because the blocks of ocean have |
13     C | overlap regions they are independent of one another. |
14     C | If terms involving lateral integrals are needed in this |
15     C | routine care will be needed. Similarly finite-difference |
16     C | operations with stencils wider than the overlap region |
17     C | require special consideration. |
18     C | Notes |
19     C | ===== |
20     C | C*P* comments indicating place holders for which code is |
21     C | presently being developed. |
22     C \==========================================================/
23    
24     C == Global variables ===
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "CG2D.h"
28 adcroft 1.6 #include "PARAMS.h"
29 adcroft 1.3 #include "DYNVARS.h"
30 cnh 1.1
31     C == Routine arguments ==
32     C myThid - Thread number for this instance of the routine.
33     INTEGER myThid
34    
35     C == Local variables
36     C xA, yA - Per block temporaries holding face areas
37     C uTrans, vTrans, wTrans - Per block temporaries holding flow transport
38     C o uTrans: Zonal transport
39     C o vTrans: Meridional transport
40     C o wTrans: Vertical transport
41     C maskC,maskUp o maskC: land/water mask for tracer cells
42     C o maskUp: land/water mask for W points
43     C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
44     C mTerm, pTerm, tendency equations.
45     C fZon, fMer, fVer[STUV] o aTerm: Advection term
46     C o xTerm: Mixing term
47     C o cTerm: Coriolis term
48     C o mTerm: Metric term
49     C o pTerm: Pressure term
50     C o fZon: Zonal flux term
51     C o fMer: Meridional flux term
52     C o fVer: Vertical flux term - note fVer
53     C is "pipelined" in the vertical
54     C so we need an fVer for each
55     C variable.
56     C iMin, iMax - Ranges and sub-block indices on which calculations
57     C jMin, jMax are applied.
58     C bi, bj
59     C k, kUp, kDown, kM1 - Index for layer above and below. kUp and kDown
60     C are switched with layer to be the appropriate index
61     C into fVerTerm
62     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73     _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
77     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
78     _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
79     _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
80     _RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
81 adcroft 1.3 _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82     _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 adcroft 1.4 _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84     _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85 adcroft 1.6 _RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
86     _RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
87     _RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
88     _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 cnh 1.1 INTEGER iMin, iMax
90     INTEGER jMin, jMax
91     INTEGER bi, bj
92     INTEGER i, j
93     INTEGER k, kM1, kUp, kDown
94    
95     C-- Set up work arrays with valid (i.e. not NaN) values
96     C These inital values do not alter the numerical results. They
97     C just ensure that all memory references are to valid floating
98     C point numbers. This prevents spurious hardware signals due to
99     C uninitialised but inert locations.
100     DO j=1-OLy,sNy+OLy
101     DO i=1-OLx,sNx+OLx
102 adcroft 1.5 xA(i,j) = 0. _d 0
103     yA(i,j) = 0. _d 0
104     uTrans(i,j) = 0. _d 0
105     vTrans(i,j) = 0. _d 0
106     aTerm(i,j) = 0. _d 0
107     xTerm(i,j) = 0. _d 0
108     cTerm(i,j) = 0. _d 0
109     mTerm(i,j) = 0. _d 0
110     pTerm(i,j) = 0. _d 0
111     fZon(i,j) = 0. _d 0
112     fMer(i,j) = 0. _d 0
113 cnh 1.1 DO K=1,nZ
114 adcroft 1.5 pH (i,j,k) = 0. _d 0
115 adcroft 1.6 K13(i,j,k) = 0. _d 0
116     K23(i,j,k) = 0. _d 0
117     K33(i,j,k) = 0. _d 0
118 cnh 1.1 ENDDO
119 adcroft 1.5 rhokm1(i,j) = 0. _d 0
120     rhokp1(i,j) = 0. _d 0
121 cnh 1.1 ENDDO
122     ENDDO
123    
124     DO bj=myByLo(myThid),myByHi(myThid)
125     DO bi=myBxLo(myThid),myBxHi(myThid)
126    
127 cnh 1.7 C-- Boundary condition on hydrostatic pressure is pH(z=0)=0
128 adcroft 1.3 DO j=1-OLy,sNy+OLy
129     DO i=1-OLx,sNx+OLx
130     pH(i,j,1) = 0. _d 0
131 adcroft 1.6 K13(i,j,1) = 0. _d 0
132     K23(i,j,1) = 0. _d 0
133     K33(i,j,1) = 0. _d 0
134     KapGM(i,j) = 0. _d 0
135 adcroft 1.3 ENDDO
136     ENDDO
137    
138 cnh 1.7 C-- Set up work arrays that need valid initial values
139     DO j=1-OLy,sNy+OLy
140     DO i=1-OLx,sNx+OLx
141     wTrans(i,j) = 0. _d 0
142     fVerT(i,j,1) = 0. _d 0
143     fVerT(i,j,2) = 0. _d 0
144     fVerS(i,j,1) = 0. _d 0
145     fVerS(i,j,2) = 0. _d 0
146     fVerU(i,j,1) = 0. _d 0
147     fVerU(i,j,2) = 0. _d 0
148     fVerV(i,j,1) = 0. _d 0
149     fVerV(i,j,2) = 0. _d 0
150     ENDDO
151     ENDDO
152    
153 cnh 1.1 iMin = 1-OLx+1
154     iMax = sNx+OLx
155     jMin = 1-OLy+1
156     jMax = sNy+OLy
157    
158 adcroft 1.4 C-- Calculate gradient of surface pressure
159     CALL GRAD_PSURF(
160     I bi,bj,iMin,iMax,jMin,jMax,
161     O pSurfX,pSurfY,
162     I myThid)
163    
164     C-- Update fields in top level according to tendency terms
165 cnh 1.1 CALL TIMESTEP(
166 adcroft 1.4 I bi,bj,iMin,iMax,jMin,jMax,1,pSurfX,pSurfY,myThid)
167 cnh 1.1
168 cnh 1.7 C-- Density of 1st level (below W(1)) reference to level 1
169     CALL FIND_RHO(
170     I bi, bj, iMin, iMax, jMin, jMax, 1, 1, 'LINEAR',
171     O rhoKm1,
172     I myThid )
173     C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
174     CALL CALC_PH(
175     I bi,bj,iMin,iMax,jMin,jMax,1,rhoKm1,rhoKm1,
176     U pH,
177 adcroft 1.5 I myThid )
178    
179 adcroft 1.3 DO K=2,Nz
180 adcroft 1.4 C-- Update fields in Kth level according to tendency terms
181     CALL TIMESTEP(
182     I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)
183 cnh 1.7 C-- Density of K-1 level (above W(K)) reference to K level
184     CALL FIND_RHO(
185     I bi, bj, iMin, iMax, jMin, jMax, K-1, K, 'LINEAR',
186     O rhoKm1,
187     I myThid )
188     C-- Density of K level (below W(K)) reference to K level
189     CALL FIND_RHO(
190     I bi, bj, iMin, iMax, jMin, jMax, K, K, 'LINEAR',
191     O rhoKp1,
192     I myThid )
193 adcroft 1.6 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
194 cnh 1.7 CALL CALC_ISOSLOPES(
195     I bi, bj, iMin, iMax, jMin, jMax, K,
196     I rhoKm1, rhoKp1,
197     O K13, K23, K33, KapGM,
198     I myThid )
199 cnh 1.1 C-- Calculate static stability and mix where convectively unstable
200 cnh 1.7 CALL CONVECT(
201     I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,myThid)
202     C-- Density of K-1 level (above W(K)) reference to K-1 level
203     CALL FIND_RHO(
204     I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, 'LINEAR',
205     O rhoKm1,
206     I myThid )
207     C-- Density of K level (below W(K)) referenced to K level
208     CALL FIND_RHO(
209     I bi, bj, iMin, iMax, jMin, jMax, K, K, 'LINEAR',
210     O rhoKp1,
211     I myThid )
212     C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
213     CALL CALC_PH(
214     I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,
215     U pH,
216 adcroft 1.3 I myThid )
217 cnh 1.1
218 cnh 1.7 ENDDO
219 adcroft 1.5
220 cnh 1.1 DO K = Nz, 1, -1
221     kM1 =max(1,k-1) ! Points to level above k (=k-1)
222     kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above
223     kDown=1+MOD(k,2) ! Cycles through 2,1 to point to current layer
224     iMin = 1-OLx+2
225     iMax = sNx+OLx-1
226     jMin = 1-OLy+2
227     jMax = sNy+OLy-1
228    
229     C-- Get temporary terms used by tendency routines
230     CALL CALC_COMMON_FACTORS (
231     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
232     O xA,yA,uTrans,vTrans,wTrans,maskC,maskUp,
233     I myThid)
234    
235     C-- Calculate accelerations in the momentum equations
236     CALL CALC_MOM_RHS(
237     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
238     I xA,yA,uTrans,vTrans,wTrans,maskC,
239     I pH,
240     U aTerm,xTerm,cTerm,mTerm,pTerm,
241     U fZon, fMer, fVerU, fVerV,
242     I myThid)
243    
244     C-- Calculate active tracer tendencies
245     CALL CALC_GT(
246     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
247     I xA,yA,uTrans,vTrans,wTrans,maskUp,
248 adcroft 1.6 I K13,K23,K33,KapGM,
249 cnh 1.1 U aTerm,xTerm,fZon,fMer,fVerT,
250     I myThid)
251     Cdbg CALL CALC_GS(
252     Cdbg I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
253     Cdbg I xA,yA,uTrans,vTrans,wTrans,maskUp,
254 adcroft 1.6 Cdbg I K13,K23,K33,KapGM,
255 cnh 1.1 Cdbg U aTerm,xTerm,fZon,fMer,fVerS,
256     Cdbg I myThid)
257    
258 cnh 1.7 ENDDO
259 cnh 1.1
260     ENDDO
261     ENDDO
262 adcroft 1.6
263     !dbg write(0,*) 'dynamics: pS',minval(cg2d_x),maxval(cg2d_x)
264     !dbg write(0,*) 'dynamics: U',minval(uVel(1:sNx,1:sNy,:,:,:)),
265     !dbg & maxval(uVel(1:sNx,1:sNy,:,:,:))
266     !dbg write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)),
267     !dbg & maxval(vVel(1:sNx,1:sNy,:,:,:))
268     !dbg write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)),
269     !dbg & maxval(gT(1:sNx,1:sNy,:,:,:))
270     !dbg write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)),
271     !dbg & maxval(Theta(1:sNx,1:sNy,:,:,:))
272     !dbg write(0,*) 'dynamics: pH',minval(pH/(Gravity*Rhonil)),
273     !dbg & maxval(pH/(Gravity*Rhonil))
274 cnh 1.1
275     RETURN
276     END

  ViewVC Help
Powered by ViewVC 1.1.22