/[MITgcm]/MITgcm/compare01/src/update_t.F
ViewVC logotype

Annotation of /MITgcm/compare01/src/update_t.F

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


Revision 1.2 - (hide annotations) (download)
Mon Jun 15 05:13:55 1998 UTC (26 years ago) by cnh
Branch: MAIN
CVS Tags: branch-atmos-merge-phase6, checkpoint24, checkpoint7, checkpoint26, branch-atmos-merge-start, checkpoint27, checkpoint9, checkpoint8, checkpoint11, checkpoint10, checkpoint13, checkpoint12, checkpoint15, checkpoint18, checkpoint17, checkpoint16, checkpoint19, checkpoint32, checkpoint31, branch-atmos-merge-zonalfilt, branch-atmos-merge-shapiro, branch-atmos-merge-freeze, branch-point-rdot, checkpoint14, checkpoint28, checkpoint29, branch-atmos-merge-phase5, branch-atmos-merge-phase4, branch-atmos-merge-phase7, checkpoint23, branch-atmos-merge-phase1, checkpoint25, branch-atmos-merge-phase3, branch-atmos-merge-phase2, checkpoint20, checkpoint21, checkpoint22
Branch point for: branch-atmos-merge, checkpoint7-4degree-ref, branch-rdot
Changes since 1.1: +9 -7 lines
Fairly coplete 4 degree global intercomparison
setup.
 Includes changes to make convective adjustment and hydrostatic
pressure correct as well as IO for climatological datasets

1 cnh 1.2 C $Id: update_t.F,v 1.1 1998/05/25 20:21:06 cnh Exp $
2 cnh 1.1 #include "CPP_OPTIONS.h"
3     #include "CPP_MACROS.h"
4    
5     #define _DB0 0.
6    
7     C===================================================================
8     C Procedure name: UPDATE_T |
9     C Function: Step forward the temperature field. |
10     C Comments: |
11     C===================================================================
12     CStartofinterface
13     SUBROUTINE UPDATE_T(
14     & U, V, W, K, T, GT )
15     IMPLICIT NONE
16     C============== Global data ==========================================
17     #include "SIZE.h"
18     #include "AJAINF.h"
19     #include "GRID.h"
20     #include "PARAMS.h"
21     #include "OPERATORS.h"
22     #include "OLDG.h"
23     #include "MASKS.h"
24     #include "FORCING.h"
25     C============= Routine arguments =====================================
26     C /--------------------------------------------------------------\
27     C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). |
28     C | K - Working level. |
29     C | T - Potential temperature (oC). |
30     C \--------------------------------------------------------------/
31     REAL U (Nx,Ny,Nz)
32     REAL V (Nx,Ny,Nz)
33     REAL W (Nx,Ny,Nz)
34     INTEGER K
35     REAL T (Nx,Ny,Nz)
36     REAL gt (Nx,Ny,Nz)
37     CEndofinterface
38     C============= Local variables ======================================
39     C /--------------------------------------------------------------\
40     C |Slice Notation |
41     C |phiK - XY slice of variable "phi" at level K. |
42     C |phiKP1 - XY slice of variable "phi" at level K+1. |
43     C |phiKM1 - XY slice of variable "phi" at level K-1. |
44     C |phiNM1 - Variable "phi" at time level N-1. |
45     C |phiU - Varaible "phi" interpolated to U grid. |
46     C |phiV - Varaible "phi" interpolated to V grid. |
47     C \--------------------------------------------------------------/
48     C /--------------------------------------------------------------\
49     C |Geometry terms |
50     C | o xa,ya,za - X, Y, Z face areas ( m.Pa, m.Pa, m**2 ). |
51     C | o rV - 1/Volume ( m**2.Pa ) |
52     C | o ?Grad?Op - Generic gradient operator A[XYZ][t].d/d[xyz] |
53     C | tGradxOp: AXt/DXt, tGradyOp: AYt/DYt etc... |
54     C \--------------------------------------------------------------/
55     REAL xaK (Nx+1,0:Ny)
56     REAL yaK (0:Nx,Ny+1)
57     REAL tGradxOp (Nx ,Ny )
58     REAL tGradyOp (Nx ,Ny )
59     REAL tGradzOpK (Nx ,Ny )
60     REAL tGradzOpKP1(Nx ,Ny )
61     REAL zaK (0:Nx,0:Ny)
62     REAL zaKP1 (Nx,Ny)
63     REAL rVk (Nx ,Ny )
64     C /--------------------------------------------------------------\
65     C |State variables |
66     C | o u, v, w - X,Y,Z velocity slices. |
67     C | o t - Potential temperature. |
68     C | o ph - Hydrostatic pressure. |
69     C \--------------------------------------------------------------/
70     REAL uK (Nx+1,0:Ny)
71     REAL wK (0:Nx,0:Ny)
72     REAL wkP1 (0:Nx,0:Ny)
73     REAL vK (0:Nx,Ny+1)
74     REAL tK (0:Nx+1,0:Ny+1)
75     REAL phK (60)
76     C /--------------------------------------------------------------\
77     C |Work space arrays |
78     C \--------------------------------------------------------------/
79     REAL tmp (0:Nx+1,0:Ny+1)
80     REAL tmp2 (0:Nx+1,0:Ny+1)
81     C /--------------------------------------------------------------\
82     C |Arrays for accumulating fluxes |
83     C | o fX, fY, fZ - Flux in X, Y and Z. |
84     C \--------------------------------------------------------------/
85     REAL fX (0:Nx+1,Ny)
86     REAL fY (Nx,0:Ny+1)
87     REAL fZK (Nx,Ny)
88     REAL fZKP1 (Nx,Ny)
89     C Loop counters:
90     INTEGER I, J
91     C Flags:
92     LOGICAL TOP_LAYER
93     LOGICAL BOTTOM_LAYER
94     C === Set up controlling flags ========================================
95     TOP_LAYER = K .EQ. 1
96     BOTTOM_LAYER = K .EQ. Nz
97     C _D(( ' K, Nz ', K, Nz ))
98     C _D(( ' BOTTOM_LAYER ', BOTTOM_LAYER ))
99     C === Set up the slices ===============================================
100     C /---------------------------------------------------------------\
101     C | xaK <- xa(K) |
102     C | yaK <- ya(K) |
103     C | za[K,KP1] <- za[(K),(K+1)] |
104     C \---------------------------------------------------------------/
105     DO J=1,Ny
106     DO I=1,Nx
107     xaK (I,J) = xa (_I3(K,I,J))
108     yaK (I,J) = ya (_I3(K,I,J))
109     zaK (I,J) = za (_I3(K,I,J))
110     ENDDO
111     ENDDO
112     IF ( .NOT. BOTTOM_LAYER ) THEN
113     DO J=1,Ny
114     DO I=1,Nx
115     zaKP1(I,J) = za(_I3(K+1,I,J))
116     ENDDO
117     ENDDO
118     ENDIF
119     fZK = 0.
120     fZKP1 = 0.
121     C /---------------------------------------------------------------\
122     C | tGradXOp <- 1/DX*xa |
123     C | tGradYOp <- 1/DY*ya |
124     C | tGradZOp[K,KP1] <- 1/DZ*ZA[K,K+1] |
125     C \---------------------------------------------------------------/
126     DO J = 1, Ny
127     DO I = 1, Nx
128     tGradxOp (I,J) = xaK(I,J)*pDdxOp(_I3(K,I,J))
129     tGradyOp (I,J) = yaK(I,J)*pDdyOp(_I3(K,I,J))
130     tGradzOpK(I,J) = zaK(I,J)*pDdzOp(_I3(K,I,J))
131     ENDDO
132     ENDDO
133     IF ( .NOT. BOTTOM_LAYER ) THEN
134     DO J = 1, Ny
135     DO I = 1, Nx
136     tGradzOpKP1(I,J) = zaKP1(I,J)*pDdzOp(I,J,K+1)
137     ENDDO
138     ENDDO
139     ENDIF
140     C /---------------------------------------------------------------\
141     C | w[K,KP1] <- W[(K),(K+1)] |
142     C \---------------------------------------------------------------/
143     wKP1 = 0.
144     DO J=1,Ny
145     DO I=1,Nx
146     wK(I,J) = W(_I3(K,I,J))
147     ENDDO
148     ENDDO
149     IF ( .NOT. BOTTOM_LAYER ) THEN
150     DO J=1,Ny
151     DO I=1,Nx
152     wKP1(I,J) = W(_I3(K+1,I,J))
153     ENDDO
154     ENDDO
155     ENDIF
156     C /---------------------------------------------------------------\
157     C | uK <- U(K) |
158     C \---------------------------------------------------------------/
159     DO J=1,Ny
160     DO I=1,Nx
161     uK(I,J) = U(_I3(K,I,J))
162     ENDDO
163     ENDDO
164     C /---------------------------------------------------------------\
165     C | vK <- V(K) |
166     C \---------------------------------------------------------------/
167     DO J=1,Ny
168     DO I=1,Nx
169     vK(I,J) = V(_I3(K,I,J))
170     ENDDO
171     ENDDO
172     C /---------------------------------------------------------------\
173     C | tK <- T(K) |
174     C \---------------------------------------------------------------/
175     DO J=1,Ny
176     DO I=1,Nx
177     tK(I,J) = T(_I3(K,I,J))
178     ENDDO
179     ENDDO
180     DO J=1,Ny
181     tK(0,J) = tK(Nx,J)
182     tK(Nx+1,J) = tK(1 ,J)
183     ENDDO
184     DO I=0,Nx+1
185     tK(I,0) = tK(I , Ny)
186     tK(I,Ny+1) = tK(I , 1)
187     ENDDO
188     C /---------------------------------------------------------------\
189     C | Gt Section |
190     C |===============================================================|
191     C | In this section Gt = GtDiffusion |
192     C | + GtAdvection |
193     C | + GtForcing |
194     C \---------------------------------------------------------------/
195     C /---------------------------------------------------------------\
196     C | Initialise local terms |
197     C | rVk <- 1/V |
198     C \---------------------------------------------------------------/
199     DO J=1,Ny
200     DO I=1,Nx
201     rVk(I,J) = rPvol(_I3(K,I,J))
202     ENDDO
203     ENDDO
204     fZK = 0.
205     fZKP1 = 0.
206     C /---------------------------------------------------------------\
207     C | XY GtDiffusion + XY GtAdvection |
208     C | Note: XY GtDiffusion = Laplacian + Biharmonic |
209     C | o Laplacian ( -a2TempXY.del^2(u) ) |
210     C | o Biharmonic ( +a4TempXY.del^4(u) ) |
211     C \---------------------------------------------------------------/
212     DO J = 1, Ny
213     DO I = 1, Nx
214     C /-------------------------------------------------------------\
215     C | fX <- 1/DX*XA*Ddx{T} |
216     C | fY <- 1/DY*YA*Ddy{T} |
217     C \-------------------------------------------------------------/
218     fX(I,J) = tGradxOp(I,J)*(tK(I,J)-tK(I-1,J))
219     fY(I,J) = tGradyOp(I,J)*(tK(I,J)-tK(I,J-1))
220     ENDDO
221     ENDDO
222     DO J=1,Ny
223     fX(Nx+1,J) = fX(1 ,J)
224     ENDDO
225     DO I=1,Nx
226     fY(I,Ny+1) = fY(I, 1)
227     ENDDO
228     DO J=1,Ny
229     DO I=1,Nx
230     C /-------------------------------------------------------------\
231     C | tmp <- 1/V*(Dx{fX}+Dy{fY}) |
232     C \-------------------------------------------------------------/
233     tmp(I,J)= rVk(I,J)*fX(I+1,J)-rVk(I,J)*fX(I,J)
234     & +rVk(I,J)*fY(I,J+1)-rVk(I,J)*fY(I,J)
235     ENDDO
236     ENDDO
237     tmp(0,1:Ny) = tmp(Nx,1:Ny)
238     tmp(1:Nx,0) = tmp(1:Nx,Ny)
239     DO J = 1,Ny
240     DO I = 1,Nx
241     C /-------------------------------------------------------------\
242     C | fX <- -k2*fX + ( Laplacian,<Line 1> ) |
243     C | k4/DX*XA*Ddx{tmp} + ( Biharmonic,<Line 2> ) |
244     C | XA*UBx{T} ( Advection,<Line 3> ) |
245     C \-------------------------------------------------------------/
246     fX(I,J) = 0.
247     _LPT(& -a2TempXY*fX(I,J) )
248     _BHT(& +a4TempXY*tGradxOp(I,J)*(tmp(I,J)-tmp(I-1,J)) )
249     _ADT(& +0.5*uK(I,J)*XAK(I,J)*( tK(I-1,J)+tK(I,J) ) )
250     C /-------------------------------------------------------------\
251     C | fY <- -k2*fY + ( Laplacian,<Line 1> ) |
252     C | k4/DY*YA*Ddy{tmp} + ( Biharmonic,<Line 2> ) |
253     C | YA*V*By{T} ( Advection,<Line 3> ) |
254     C \-------------------------------------------------------------/
255     fY(I,J) = 0.
256     _LPT(& -a2TempXY*fY(I,J) )
257     _BHT(& +a4TempXY*tGradyOp(I,J)*(tmp(I,J)-tmp(I,J-1)) )
258     _ADT(& +0.5*vK(I,J)*yaK(I,J)*( tK(I,J-1)+tK(I,J) ) )
259     ENDDO
260     ENDDO
261     C /---------------------------------------------------------------\
262     C | Z GtDiffusion + Z GtAdvection |
263     C \---------------------------------------------------------------/
264     IF ( .NOT. TOP_LAYER ) THEN
265     DO J = 1, Ny
266     DO I = 1, Nx
267     C /------------------------------------------------------------\
268     C | fZK <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) |
269     C | ZA*-w*Bz{T} ( Advection,<Line 2> ) |
270     C \------------------------------------------------------------/
271     fZK(I,J) = 0.
272     _LPT(& -a2TempZ*tGradzOpK(I,J)*(T(_I3(K-1,I,J))-tK(I,J) ) )
273     _ADT(& +0.5*(tK(I,J)+T(_I3(K-1,I,J)))*( -wK(I,J)*ZAK(I,J) ) )
274     ENDDO
275     ENDDO
276     ELSE
277     DO J = 1, Ny
278     DO I = 1, Nx
279     C /------------------------------------------------------------\
280     C | fZK <- ZA*-w*Bz{T} ( Advection,<Line 1> ) |
281     C \------------------------------------------------------------/
282     fZK(I,J) = 0.
283     _ADT(& +0.5*freeSurfFac*(tK(I,J)+tK(I,J))*(-wK(I,J)*ZAK(I,J) ) )
284     ENDDO
285     ENDDO
286     ENDIF
287     IF ( .NOT. BOTTOM_LAYER ) THEN
288     DO J = 1, Ny
289     DO I = 1, Nx
290     C /------------------------------------------------------------\
291     C | fZKP1 <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) |
292     C | ZA*-w*Bz{T} ( Advection,<Line 2> ) |
293     C \------------------------------------------------------------/
294     fZKP1(I,J) = 0.
295     _LPT(& -a2TempZ*tGradzOpKP1(I,J)*( tK(I,J)-T(_I3(K+1,I,J)) ) )
296     _ADT(& +0.5*(tK(I,J)+T(_I3(K+1,I,J)))*(-wKP1(I,J))*zaKP1(I,J) )
297     ENDDO
298     ENDDO
299     ENDIF
300     C /---------------------------------------------------------------\
301     C | gT <- GtDiffusion + GtAdvection = -div(fX,fY,fZ) |
302     C \---------------------------------------------------------------/
303     fX(Nx+1,1:Ny) = fX(1,1:Ny)
304     fY(:,Ny+1) = fY(:,1)
305     DO J = 1, Ny
306     DO I = 1, Nx
307     C /-------------------------------------------------------------\
308     C | gT <- -1/V*(Ddx{fX}+Ddy{fY}) |
309     C \-------------------------------------------------------------/
310     gt(_I3(K,I,J))=-fX(I+1,J)*rVk(I,J)+fX(I,J)*rVk(I,J)
311     & -fY(I,J+1)*rVk(I,J)+fY(I,J)*rVk(I,J)
312     ENDDO
313     ENDDO
314     DO J = 1, Ny
315     DO I = 1, Nx
316     C /------------------------------------------------------------\
317     C | gT <- gT-1/Vu*fZK |
318     C \------------------------------------------------------------/
319     gt(_I3(K,I,J))=gt(_I3(K,I,J))-fZK(I,J)*rVk(I,J)
320     ENDDO
321     ENDDO
322     IF ( .NOT. BOTTOM_LAYER ) THEN
323     DO J = 1, Ny
324     DO I = 1, Nx
325     C /------------------------------------------------------------\
326     C | gT <- gT+1/Vu*fZP1 |
327     C \------------------------------------------------------------/
328     gt(_I3(K,I,J))=gt(_I3(K,I,J))+fZKP1(I,J)*rVk(I,J)
329     ENDDO
330     ENDDO
331     ENDIF
332     C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MAXVAL(Gt )', MAXVAL(Gt)))
333     C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MINVAL(Gt )', MINVAL(Gt)))
334     C /---------------------------------------------------------------\
335     C | gT <- gT + GtForcing |
336     C \---------------------------------------------------------------/
337 cnh 1.2 C IF ( TOP_LAYER ) THEN
338     C DO J =1,Ny
339     C DO I=1, Nx
340 cnh 1.1 Cdbg gT(_I3(K,I,J)) = gT(_I3(K,I,J)) + HEAT(I,J)
341 cnh 1.2 C gT(_I3(K,I,J)) = gT(_I3(K,I,J)) -
342     C & 1./(30.*oneDay)*( T(_I3(K,I,J)) - Tsurf(I,J,K) )
343     C ENDDO
344     C ENDDO
345     C ENDIF
346 cnh 1.1
347     C Final time derivative from AB2.
348     gt(_I3(K,1:Nx,1:Ny)) = Gt(_I3(K,1:Nx,1:Ny))*pMask(_I3(K,:,:))
349     tmp = 0.
350     tmp2 = 0.
351     tmp (1:Nx,1:Ny) = gtNM1(_I3(K,:,:))
352     tmp2(1:Nx,1:Ny) = gt (_I3(K,:,:))
353    
354     C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(GtNM1)', MAXVAL(GtNm1) ))
355     C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(Gt)', MAXVAL(Gt) ))
356    
357     gtNM1(_I3(K,:,:)) = gt(_I3(K,:,:))
358     tmp2=(1.5+abEps)*tmp2-(0.5+abEps)*tmp
359     gt(_I3(K,1:Nx,1:Ny))=tmp2(1:Nx,1:Ny)
360     C
361     RETURN
362     END

  ViewVC Help
Powered by ViewVC 1.1.22