/[MITgcm]/MITgcm_contrib/jscott/igsm/src/comp1.F
ViewVC logotype

Contents of /MITgcm_contrib/jscott/igsm/src/comp1.F

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


Revision 1.1 - (show annotations) (download)
Fri Aug 11 19:35:30 2006 UTC (18 years, 11 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
atm2d package

1
2 #include "ctrparam.h"
3
4 ! ==========================================================
5 !
6 ! COMP1.F: NUMERICALLY INTEGRATING THE NON-SOURCE TERMS OF
7 ! THE 2-DIMENSIONAL CLIMATE MODEL
8 !
9 ! ----------------------------------------------------------
10 !
11 ! Author of Chemistry Modules: Chien Wang
12 !
13 ! ----------------------------------------------------------
14 !
15 ! Revision History:
16 !
17 ! When Who What
18 ! ---- ---------- -------
19 ! 073100 Chien Wang repack based on CliChem3 and add cpp
20 !
21 ! ==========================================================
22
23
24 SUBROUTINE COMP1 (UT,VT,TT,PT,QT,U,V,T,P,Q,DT1,NS) 2001.
25 C 2001.5
26 C NUMERICALLY INTEGRATING THE NON-SOURCE TERMS OF 2002.
27 C THE 2-DIMENSIONAL CLIMATE MODEL 2002.5
28 C 2003.
29
30 #if ( defined CPL_CHEM )
31 !
32 #include "chem_para"
33 #include "chem_com"
34 !
35 #endif
36
37 #include "BD2G04.COM" 2003.5
38
39 COMMON/WORK1/PIT(IM0,JM0),SD(IM0,JM0,LM0-1) 2004.
40 COMMON/WORK3/PHI(IM0,JM0,LM0),SPA(IM0,JM0,LM0) 2004.5
41 COMMON/WORK4/FD(IM0,JM0),FLUXQ(72),DUMMYS(72),DUMMYN(72) 2005.
42 COMMON/WORK5/DUT(IO0,JM0,LM0),DVT(IO0,JM0,LM0) 2005.5
43 DIMENSION UT(IM0,JM0,LM0),VT(IM0,JM0,LM0),TT(IM0,JM0,LM0)
44 & ,PT(IM0,JM0), 2006.
45 * QT(IM0,JM0,LM0),PU(IM0,JM0),PV(IM0,JM0),CONV(IM0,JM0,LM0) 2006.5
46 DIMENSION SINL(JM0),COSL(JM0), 2007.
47 * DXU(JM0),DYU(JM0),DXYU(JM0),PHIS(IO0,JM0),UX(IO0,JM0,2*LM0) 2007.5
48 EQUIVALENCE (SINL(1),SINP(1)),(COSL(1),COSP(1)), 2008.
49 * (DXU(1),DXV(1)),(DYU(1),DYV(1)),(DXYU(1),DXYV(1)), 2008.5
50 * (PHIS,FDATA) 2009.
51 EQUIVALENCE (CONV,PIT) 2009.5
52 c REAL*4 KAPAP1 2010.
53 REAL KAPAP1
54 COMMON/EPARA/VTH(JM0,LM0),WTH(JM0,LM0),VU(JM0,LM0),VV(JM0,LM0)
55 & ,DQSDT(JM0,LM0) 2010.5
56 * ,DWV(JM0),PHIT(JM0,LM0),TPRIM2(JM0,LM0),WU(JM0,LM0),CKS,CKN 2011.
57 * ,WQ(JM0,LM0),VQ(JM0,LM0),MRCHT 2011.5
58 COMMON/INTA/COE1(01,01,01),COE2(01,01,01) 2012.
59 COMMON/SPEC1/ 2012.5
60 * XA(IM0,JM0),XB(IM0,JM0),YA(IM0,JM0),YB(IM0,JM0),ZA(IM0,JM0)
61 & ,ZB(IM0,JM0) 2013.
62 COMMON/SPEC2/KM,KINC,COEK,C3LAND(IO0,JM0),C3OICE(IO0,JM0) 2013.5
63 * ,C3LICE(IO0,JM0),WMGE(IO0,JM0),TSSFC(IO0,JM0,4) 2014.
64 LOGICAL SKIPSI,HPRNT,EDDFLX,FIRST,HPRNT1 2014.5
65 common/PRNT1/NCOMP
66 common/conprn/HPRNT1,JPR,LPR
67 data FIRST /.true./
68 C 2015.
69 HPRNT=TAU.gt.115.0.and.TAU.lt.125.0
70 HPRNT=TAU.lt.125.0
71 HPRNT=.false.
72 NCOMP=NCOMP+1
73 KBGN=KINC+1 2015.5
74 KM2=KM*2-1 2016.
75 KM3=KM2 2016.5
76 INEDDY=1
77 c ICHK=(1+(NDYN-1)*2)*5 2016.6
78 ICHK=(1+(NDYN-1)*2)*INEDDY
79 ICHK1=(1+(NDYN-1)*2)
80 IX=IM 2017.
81 IF(KM.EQ.1) IX=1 2017.5
82 IS=IX 2018.
83 FIS=IS 2018.5
84 CP=RGAS/KAPA 2019.
85 SHA=RGAS/KAPA 2019.5
86 KAPAP1=KAPA+1. 2020.
87 JMM2=JM-2 2020.5
88 NTRACE=0 2021.
89 NTLM=NTRACE*LM 2021.5
90 NTP1LM=NTLM+LM 2022.
91 NCOMP3=NDYN 2022.5
92 LMT2=LM+LM 2023.
93 LMM2=LM-2 2023.5
94 c DTE=NDYN*DT*5 2024.
95 DTE=NDYN*DT*INEDDY
96 MRCHT=MRCHT+IABS(MRCH) 2024.1
97 DTT2=DT1*2. 2024.5
98 DT2=DT1/2. 2025.
99 DT4=DT1/4. 2025.5
100 DT8=DT1/8. 2026.
101 DT12=DT1/12. 2026.5
102 DT24=DT1/24. 2027.
103 DTCP=DT1/CP 2027.5
104 DTDLN=DT1*DLON 2028.
105 DTDLN2=DT2*DLON 2028.5
106 SKIPSI=.TRUE. 2029.
107 FMU=10.E7 2029.5
108 FXCO=DT4/CP 2030.
109 FXCO1=DT2/CP 2030.5
110 HLAT=LHE 2031.
111 CLH=HLAT/CP 2031.5
112 if(FIRST)then
113 print *,' EDDY DIFFUSION is CALCULATED EVERY ',INEDDY,' HOURS'
114 print *,' DTE=',DTE,DTE/3600.
115 print *,' ICHK=',ICHK
116 FIRST=.false.
117 endif
118 DO 5 L=1,LMT2 2032.
119 DO 5 J=1,JM 2032.5
120 DO 5 I=1,IO 2033.
121 5 UX(I,J,L)=U(1,J,L) 2033.5
122 IF(MODD5K.LT.MRCH) CALL DIAG5F(UX) 2034.
123 C 2034.5
124 C SCALE THE PROGNOSTIC VARIABLES 2035.
125 C 2035.5
126 DO 50 J=1,JM 2036.
127 DO 50 I=1,IS 2036.5
128 50 FD(I,J)=PT(I,J)*DXYP(J) 2037.
129
130 #if ( defined CPL_CHEM )
131 !
132 ! --- get PAI for chemical advection:
133 !
134 if(mrch.ne.0)goto 6001
135 do i=1,n2dh
136 p00(i,1)=fd(i,1)
137 !
138 ! 092595
139 !
140 p4chem0(i,1) = pt(i,1)
141 enddo
142 6001 continue
143 !
144 #endif
145
146 c print *,'From comp1 TAU=',TAU,' MRCH=',MRCH
147 c print '(12f6.1/11f6.1)',(P(1,j),j=1,jm)
148 c print *,' U(jm-1),U(jm),v(jm)'
149 c print '(11f6.1)',(U(1,JM-1,l),l=1,lm)
150 c print '(11f6.1)',(U(1,JM,l),l=1,lm)
151 c print '(11f6.1)',(V(1,JM,l),l=1,lm)
152 c print *,' UT(jm-1),UT(jm),VT(jm)'
153 c print '(11f6.1)',(UT(1,JM-1,l),l=1,lm)
154 c print '(11f6.1)',(UT(1,JM,l),l=1,lm)
155 c print '(11f6.1)',(VT(1,JM,l),l=1,lm)
156 IF(MRCH.EQ.0) GO TO 58 2037.5
157 DO 57 L=1,NTP1LM 2038.
158 DO 57 J=1,JM 2038.5
159 DO 57 I=1,IS 2039.
160 57 QT (I,J,L)=QT (I,J,L)*FD(I,J) 2039.5
161 58 CONTINUE 2040.
162 DO 55 L=1,LM 2040.5
163 DO 55 J=1,JM 2041.
164 DO 55 I=1,IX 2041.5
165 55 TT(I,J,L)=TT(I,J,L)*FD(I,J) 2042.
166 DO 56 I=1,IX 2042.5
167 FD(I,1)=2.*FD(I,1) 2043.
168 56 FD(I,JM)=2.*FD(I,JM) 2043.5
169 DO 65 J=2,JM 2044.
170 DO 65 I=1,IX 2044.5
171 FDU=.5*(FD(I,J)+FD(I,J-1)) 2045.
172 DO 65 L=1,LMT2 2045.5
173 65 UT(I,J,L)=UT(I,J,L)*FDU 2046.
174 DO 66 L=1,LMT2 2046.5
175 DO 66 J=2,JM 2047.
176 DO 66 I=1,IO 2047.5
177 66 DUT(I,J,L)=0. 2048.
178 if(HPRNT)then
179 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH,' J=',JPR,' L=',LPR
180 print *,' after 66'
181 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
182 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
183 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
184 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
185 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
186 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
187 endif
188 C 2048.5
189 C BEGINNING OF LAYER LOOP 2049.
190 C 2049.5
191 ccc L=LM 2050.
192 c print *,' before 5934 MRCH=',MRCH
193 c RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
194 c print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
195 c RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
196 c print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
197 c print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
198 do 5934 l25=1,LM
199 L=LM+1-l25
200 C 2050.5
201 C COMPUTATION OF MASS FLUX P,T,PU PRIMARY GRID ROW 2051.
202 C PV,U,V SECONDARY GRID ROW 2051.5
203 C 2052.
204 2150 CONTINUE 2052.5
205 DO 2168 J=2,JM 2053.
206 TEMP=DXU(J)*0.5 2053.5
207 DO 2168 I=1,IX 2054.
208 PV(I,J)=(P(I,J)+P(I,J-1))*V(I,J,L)*TEMP 2054.5
209
210 #if ( defined CPL_CHEM )
211 !
212 ! --- get V for chemical advection:
213 !
214 if(mrch.eq.0)then
215 pvv(i,j,l)=v(i,j,l)
216 endif
217 !
218 #endif
219
220 2168 continue
221 if(HPRNT)then
222 print *,' PV L=',L
223 print *,PV(1,JPR),P(1,JPR),P(1,JPR-1),V(1,JPR,L),TEMP
224 endif
225 C 2055.
226 C HORIZONTAL ADVECTION OF MOISTURE 2055.5
227 C 2056.
228 IF(MRCH.EQ.0) GO TO 1400 2056.5
229 DO 1360 K=L,NTP1LM,LM 2057.
230 C**** SOUTH-NORTH ADVECTION OF MOISTURE AND TRACE COMPOUNDS 2057.5
231 DO 1330 I=1,IS 2058.
232 FLUXQ(I)=DT2*PV(I,2)*(Q (I,2,K)+Q (I,1,K)) 2058.5
233 IF(FLUXQ(I).GT. QT(I,1,K)) FLUXQ(I)= QT(I,1,K) 2059.
234 IF(FLUXQ(I).LT.-.5* QT(I,2,K)) FLUXQ(I)=-.5* QT(I,2,K) 2059.5
235 1330 QT(I,1,K)= QT(I,1,K)-FLUXQ(I) 2060.
236 DO 1340 J=2,JMM2 2060.5
237 DO 1340 I=1,IS 2061.
238 FLUX=DT2*PV(I,J+1)*( Q(I,J+1,K)+ Q(I,J,K)) 2061.5
239 IF(FLUX.GT..5* QT(I,J,K)) FLUX=.5* QT(I,J,K) 2062.
240 IF(FLUX.LT.-.5* QT(I,J+1,K)) FLUX=-.5* QT(I,J+1,K) 2062.5
241 QT(I,J,K)= QT(I,J,K)+(FLUXQ(I)-FLUX) 2063.
242 1340 FLUXQ(I)=FLUX 2063.5
243 DO 1350 I=1,IS 2064.
244 FLUX=DT2*PV(I,JM)*( Q(I,JM,K)+ Q(I,JMM1,K)) 2064.5
245 IF(FLUX.GT..5* QT(I,JMM1,K)) FLUX=.5* QT(I,JMM1,K) 2065.
246 IF(FLUX.LT.- QT(I,JM,K)) FLUX=- QT(I,JM,K) 2065.5
247 QT(I,JMM1,K)= QT(I,JMM1,K)+(FLUXQ(I)-FLUX) 2066.
248 1350 QT(I,JM,K)= QT(I,JM,K)+FLUX 2066.5
249 1360 CONTINUE 2067.
250 1400 CONTINUE 2067.5
251 C 2068.
252 C HORIZONTAL ADVECTION OF HEAT 2068.5
253 C 2069.
254 DO 2215 J=2,JM 2069.5
255 DO 2215 I=1,IX 2070.
256 FLUX=DT2*PV(I,J) 2070.5
257 FLUXT=FLUX*(T(I,J,L)+T(I,J-1,L)) 2071.
258 TT(I,J,L)=TT(I,J,L)+FLUXT 2071.5
259 2215 TT(I,J-1,L)=TT(I,J-1,L)-FLUXT 2072.
260 C 2072.5
261 C MERIDIONAL ADVECTION OF U AND V-MOMENTUM 2073.
262 C 2073.5
263 DO 2336 J=2,JMM1 2074.
264 DO 2336 I=1,IX 2074.5
265 FLUX=DT4*(PV(I,J)+PV(I,J+1)) 2075.
266 FLUXU=FLUX*(U(I,J,L)+U(I,J+1,L)) 2075.5
267 C*** CONTRIBUTION BY SYMMETRIC INSTABILITY 2076.
268 IF (MRCH.EQ.0.OR.SKIPSI) GO TO 2335 2076.5
269 DUDY=(U(I,J+1,L)*COSV(J+1)-U(I,J,L)*COSV(J))/ 2077.
270 * (DYP(J)*COSP(J)) 2077.5
271 FTEMP=F(J)/DXYP(J) 2078.
272 CRI=FTEMP*(FTEMP-DUDY) 2078.5
273 IF (CRI.GE.0.) GO TO 2335 2079.
274 FLUXDI=-DT1*FMU*P(I,J)*DUDY*DXP(J) 2079.5
275 FLUXU=FLUXU+FLUXDI 2080.
276 2335 CONTINUE 2080.5
277 UT(I,J+1,L)=UT(I,J+1,L)+FLUXU 2081.
278 UT(I,J,L)=UT(I,J,L)-FLUXU 2081.5
279 FLUXV=FLUX*(V(I,J,L)+V(I,J+1,L)) 2082.
280 VT(I,J+1,L)=VT(I,J+1,L)+FLUXV 2082.5
281 VT(I,J,L)=VT(I,J,L)-FLUXV 2083.
282 if(HPRNT.and.J.eq.JPR.and.L.eq.LPR)then
283 print *,' before 2336 c',' J=',J,' L=',L
284 print *,'U(I,J+1,L)=',U(I,J+1,L),' U(I,J,L)=',U(I,J,L)
285 print *,' DUDY=',DUDY,' F(J)=',F(J),' CRI=',CRI
286 print *,'V(I,J,L)=',V(I,J,L),' V(I,J+1,L)=',V(I,J+1,L)
287 print *,' PV(I,J)=',PV(I,J),' PV(I,J+1)=',PV(I,J+1)
288 print *,' FLUX=',FLUX,' FLUXV=',FLUXV
289 print *,' FLUXU=',FLUXU,' FLUXDI=',FLUXDI
290 endif
291 2336 continue
292 if(l25.eq.-LM)then
293 print *,' after 2336'
294 RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
295 print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
296 RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
297 print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
298 print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
299 endif
300 if(HPRNT)then
301 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
302 print *,' after 2336'
303 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
304 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
305 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
306 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
307 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
308 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
309 endif
310 C DO 2337 J=2,JMM1 2083.1
311 C2337 IF(L.EQ.8) WRITE(6,2339) J,MRCH,V(1,J,8),VT(1,J,8),U(1,J,8), 2083.2
312 C * UT(1,J,8) 2083.21
313 2339 FORMAT(1X,'AFTER M ADV---V VT U UT',2I3,4E11.3) 2083.3
314 C 2083.5
315 C COMPUTE MASS CONVERGENCE 2084.
316 C 2084.5
317 DSIGL=DSIG(L) 2085.
318 PVS=PV(1,2) 2085.5
319 PVN=PV(1,JM) 2086.
320 C
321 DO 2400 J=2,JMM1 2086.5
322 DO 2400 I=1,IX 2087.
323 2400 CONV(I,J,L)=(PV(I,J)-PV(I,J+1))*DSIGL 2087.5
324 CONV(1,1,L)=-PVS*DSIGL 2088.
325 CONV(1,JM,L)=PVN*DSIGL 2088.5
326 if(HPRNT)then
327 print *,' PV 3 L=',L,' JPR=',JPR
328 print *,PV(1,JPR),PV(1,JPR+1),DSIGL,CONV(1,JPR,L)
329 print *,DSIG(1),CONV(1,1,1)
330 endif
331 C
332 c2409 L=L-1 2089.
333 c IF(L.GE.1) GO TO 2150 2089.5
334 5934 continue
335 if(HPRNT)then
336 print *,' comp1 2 TAU=',TAU,' MRCH=',MRCH
337 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
338 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
339 print *,' PV(1,J) PV(1,J+1) PIT(1,J)'
340 print *,PV(1,JPR),PV(1,JPR+1),PIT(1,JPR)
341 print *,' CONV(1,JPR,L)'
342 print *,(CONV(1,JPR,L),L=1,LM)
343 endif
344 C 2090.
345 C END OF LAYER LOOP 2090.5
346 C 2091.
347 C COMPUTE PRESSURE TENDENCY AND SIGMA DOT 2091.5
348 C 2092.
349 C PIT(I,J)=CONV(I,J,1) 2092.5
350 DO 2420 LX=1,LMM1 2093.
351 L=LM+1-LX 2093.5
352 DO 2420 J=1,JM 2094.
353 DO 2420 I=1,IS 2094.5
354 2420 PIT(I,J)=PIT(I,J)+CONV(I,J,L) 2095.
355 DO 2430 J=1,JM 2095.5
356 DO 2430 I=1,IS 2096.
357 2430 SD(I,J,LMM1)=CONV(I,J,LM)-DSIG(LM)*PIT(I,J) 2096.5
358 DO 2440 LX=2,LMM1 2097.
359 L=LM-LX 2097.5
360 DO 2440 J=1,JM 2098.
361 DO 2440 I=1,IS 2098.5
362 2440 SD(I,J,L)=SD(I,J,L+1)+CONV(I,J,L+1)-DSIG(L+1)*PIT(I,J) 2099.
363
364 #if ( defined CPL_CHEM )
365 !
366 ! --- get Omiga for chemical advection:
367 !
368 if(mrch.ne.0)goto 6009
369 do 6008 i=1,nlon
370 do 6008 k=1,nlev
371 do 6008 j=1,nlat
372 pww(i,j,k)=sd(i,j,k)/p00(i,j)
373 6008 continue
374
375 6009 continue
376 !
377 #endif
378
379 if(HPRNT)then
380 j=jm
381 print *,' comp 21 J=',J,' L=',L
382 print *,' PIT(1,J)=',PIT(1,J)
383 print *,' SD ',(SD(1,J,L),l=1,lm-1)
384 print *,' CONV(1,J,L)'
385 print *,(CONV(1,J,L),L=1,LM)
386 endif
387 C 2099.5
388 C COMPUTE THE NEW SURFACE PRESSURE 2100.
389 C 2100.5
390 DO 2450 J=1,JM 2101.
391 DO 2450 I=1,IS 2101.5
392 PTNEW=PT(I,J)+DT1*PIT(I,J)/DXYP(J) 2102.
393 IF(PTNEW.LE.1150..AND.PTNEW.GT.400.) GO TO 2450 2102.5
394 4563 CONTINUE
395 print *,' TAU=',TAU
396 WRITE(6,901) I,J,MRCH,P(I,J),PTNEW,(FDATA(1,J,K),K=1,22), 2103.
397 * (V(1,J,L),L=1,LM),(V(1,J+1,L),L=1,LM), 2103.1
398 * (T(1,J,L),L=1,LM),(Q(1,J,L),L=1,LM) 2103.2
399 901 FORMAT(' PRESSURE DIAGNOSTIC I,J,MRCH,P,PT=',3I5,2F10.1/, 2103.5
400 * 2(1X,11F10.2/),4(1X,11E11.3/))
401 c * 2(1X,11F10.2/),4(10X,9E11.3/)) 2103.6
402 stop
403 2450 PT(I,J)=PTNEW 2104.
404 C 2104.5
405 C VERTICAL ADVECTION OF MOMENTUM 2105.
406 C 2105.5
407 c print *,' before 2480'
408 c RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
409 c print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
410 c RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
411 c print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
412 c print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
413 DO 2480 L=1,LMM1 2106.
414 LP1=L+1 2106.5
415 WTS=2. 2107.
416 WTN=1. 2107.5
417 DO 2480 J=2,JM 2108.
418 IF(J.EQ.JM) WTN=2. 2108.5
419 DO 2470 I=1,IX 2109.
420 SDU=DT4*(SD(I,J,L)*WTN+SD(I,J-1,L)*WTS) 2109.5
421 SDUDN=SDU/DSIG(L) 2110.
422 SDUUP=SDU/DSIG(LP1) 2110.5
423 TEMU=U(I,J,L)+U(I,J,LP1) 2111.
424 TEMV=V(I,J,L)+V(I,J,LP1) 2111.5
425 UT(I,J,L)=UT(I,J,L)+SDUDN*TEMU 2112.
426 UT(I,J,LP1)=UT(I,J,LP1)-SDUUP*TEMU 2112.5
427 VT(I,J,L)=VT(I,J,L)+SDUDN*TEMV 2113.
428 2470 VT(I,J,LP1)=VT(I,J,LP1)-SDUUP*TEMV 2113.5
429 2480 WTS=1. 2114.
430 c print *,' after 2480'
431 c RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
432 c print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
433 c RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
434 c print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
435 c print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
436 if(HPRNT)then
437 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
438 print *,' after 2480'
439 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
440 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
441 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
442 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
443 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
444 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
445 endif
446 C DO 2481 J=2,JMM1 2114.1
447 C2481 WRITE(6,2482) J,MRCH,V(1,J,8),VT(1,J,8),U(1,J,8),UT(1,J,8), 2114.2
448 C * SD(1,J,8) 2114.21
449 2482 FORMAT(1X,'AFTER V ADV---V VT U UT SD',2I3,5E11.3) 2114.3
450 IF(MODD5K.LT.MRCH) CALL DIAG5A(4,MRCH) 2114.5
451 C IF(MRCH.GT.0) CALL DIAG9D (1,DT1,U,V) 2115.
452 C 2115.5
453 C VERTICAL ADVECTION OF MOISTURE 2116.
454 C 2116.5
455 IF(MRCH.EQ.0) GO TO 1800 2117.
456 DO 1730 L=1,LMM1 2117.5
457 LP1=L+1 2118.
458 DO 1710 J=1,JM 2118.5
459 DO 1710 I=1,IS 2119.
460 FLUX=DT1*SD(I,J,L)*2.* Q(I,J,L)* Q(I,J,LP1)/( Q(I,J,L)+ 2119.5
461 * Q(I,J,LP1)+1.E-20) 2120.
462 IF(FLUX.GT..5* QT(I,J,LP1)*DSIG(LP1)) FLUX=.5* QT(I,J,LP1)* 2120.5
463 * DSIG(LP1) 2121.
464 IF(FLUX.LT.-.5* QT(I,J,L)*DSIG(L)) FLUX=-.5* QT(I,J,L)*DSIG(L) 2121.5
465 QT(I,J,L)= QT(I,J,L)+FLUX/DSIG(L) 2122.
466 1710 QT(I,J,LP1)= QT(I,J,LP1)-FLUX/DSIG(LP1) 2122.5
467 1730 CONTINUE 2123.
468 1800 CONTINUE 2123.5
469 C 2124.
470 C CORIOLIS FORCE 2124.5
471 c print *,' before 3130'
472 c RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
473 c print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
474 c RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
475 c print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
476 c print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
477 C 2125.
478 DO 3100 I=1,IX 2125.5
479 FD(I,1)=0. 2126.
480 3100 FD(I,JM)=0. 2126.5
481 DO 3130 L=1,LM 2127.
482 DO 3110 J=2,JMM1 2127.5
483 DO 3110 I=1,IX 2128.
484 3110 FD(I,J)=F(J) +0.5*(U(I,J,L)+U(I,J+1,L))*(DXU(J)-DXU(J+1)) 2128.5
485 DO 3115 J=2,JM 2129.
486 DO 3115 I=1,IX 2129.5
487 PU(I,J)=DT4*(P(I,J)+P(I,J-1))*(FD(I,J)+FD(I,J-1)) 2130.
488 3115 CONTINUE 2132.
489 DO 3120 J=2,JM 2132.5
490 DO 3120 I=1,IX 2133.
491 UT(I,J,L)=UT(I,J,L)+ V(I,J,L)*PU(I,J) 2133.5
492 3120 VT(I,J,L)=VT(I,J,L)- U(I,J,L)*PU(I,J) 2134.
493 3130 CONTINUE 2134.5
494 c print *,' after 3130'
495 c RFDP=2./(P(1,JM-2)*DXYP(JM-2)+P(1,JM-1)*DXYP(JM-1))
496 c print '(11f6.1)',(UT(1,JM-1,l1)*RFDP,l1=1,lm)
497 c RFDP=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
498 c print '(11f6.1)',(UT(1,JM,l1)*RFDP,l1=1,lm)
499 c print '(11f6.1)',(VT(1,JM,l1)*RFDP,l1=1,lm)
500 if(HPRNT)then
501 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
502 print *,' after 3130'
503 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
504 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
505 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
506 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
507 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
508 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
509 endif
510 C DO 3131 J=2,JMM1 2134.6
511 C3131 WRITE (6,3132) J,MRCH,V(1,J,8),VT(1,J,8),U(1,J,8),UT(1,J,8) 2134.7
512 3132 FORMAT(1X,'AFTER COR---V VT U UT',2I3,4E11.3) 2134.8
513 IF(MODD5K.LT.MRCH) CALL DIAG5A(5,MRCH) 2135.
514 C IF(MRCH.GT.0) CALL DIAG9D (2,DT1,U,V) 2135.5
515 C 2136.
516 C COMPUTE SPA, PHI AND VERTICAL ADVECTION OF POTENTIAL TEMPERATURE 2136.5
517 C 2137.
518 DO 3070 J=1,JM 2137.5
519 DO 3070 I=1,IX 2138.
520 SUM1=0. 2138.5
521 SUM2=0. 2139.
522 SP=P(I,J) 2139.5
523 PDN=SIG(1)*SP+PTOP 2140.
524 PKDN=EXPBYK(PDN) 2140.5
525 FLUXDN=0. 2141.
526 DO 3040 L=1,LMM1 2141.5
527 LP1=L+1 2142.
528 SPA(I,J,L)=RGAS*SP*SIG(L)*T(I,J,L)*PKDN/PDN 2142.5
529 SUM1=SUM1+SPA(I,J,L)*DSIG(L) 2143.
530 PUP=SIG(LP1)*SP+PTOP 2143.5
531 PKUP=EXPBYK(PUP) 2144.
532 THETA=0.5*(T(I,J,LP1)+T(I,J,L)) 2144.5
533 FLUX=DT1*SD(I,J,L)*THETA 2145.
534 TT(I,J,L)=TT(I,J,L)+(FLUX-FLUXDN)/DSIG(L) 2145.5
535 FLUXDN=FLUX 2146.
536 PHI(I,J,LP1)=CP*THETA*(PKDN-PKUP) 2146.5
537 PDN=PUP 2147.
538 PKDN=PKUP 2147.5
539 3040 SUM2=SUM2+SIGE(LP1)*PHI(I,J,LP1) 2148.
540 SPA(I,J,LM)=SIG(LM)*SP*RGAS*T(I,J,LM)*PKDN/PDN 2148.5
541 SUM1=SUM1+SPA(I,J,LM)*DSIG(LM) 2149.
542 TT(I,J,LM)=TT(I,J,LM)-FLUXDN/DSIG(LM) 2149.5
543 3050 PHI(I,J,1)=PHIS(I,J)+SUM1-SUM2 2150.
544 DO 3070 L=2,LM 2150.5
545 3070 PHI(I,J,L)=PHI(I,J,L)+PHI(I,J,L-1) 2151.
546 if(HPRNT)then
547 print *,' comp1 3 TAU=',TAU,' MRCH=',MRCH
548 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
549 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
550 endif
551 C
552 C 2151.5
553 DO 3340 L=1,LM 2152.
554 C 2152.5
555 C PRESSURE GRADIENT FORCE 2153.
556 C 2153.5
557 C NORTH-SOUTH DERIVATIVE AFFECTS THE V-MOMENTUM 2154.
558 C 2154.5
559 DO 3236 J=2,JM 2155.
560 DO 3236 I=1,IX 2155.5
561 FLUX=DT2*((P(I,J)+P(I,J-1))*(PHI(I,J,L)-PHI(I,J-1,L))+ 2156.
562 * (SPA(I,J,L)+SPA(I,J-1,L))*(P(I,J)-P(I,J-1)))*DXU(J) 2156.5
563 3236 VT(I,J,L)=VT(I,J,L)-FLUX 2157.
564 3340 CONTINUE 2157.5
565 if(HPRNT)then
566 print *,' after 3340'
567 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
568 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
569 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
570 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
571 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
572 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
573 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
574 endif
575 C DO 3341 J=2,JMM1 2157.6
576 C3341 WRITE(6,3342) J,MRCH,P(1,J),VT(1,J,8),UT(1,J,8),PHI(1,J,8), 2157.7
577 C * SPA(1,J,8) 2157.71
578 3342 FORMAT(1X,'AFTER P GRAD---P VT UT PHI SPA',2I3,5E11.3) 2157.8
579 IF(MODD5K.LT.MRCH) CALL DIAG5A(6,MRCH) 2158.
580 C IF(MRCH.GT.0) CALL DIAG9D (3,DT1,U,V) 2158.5
581 C 2159.
582 C END OF THE PRESSURE GRADIENT FORCE LAYER LOOP 2159.5
583 C 2160.
584 C CALL EDDY PARA. ROUTINE IF REQUIRED 2160.5
585 C 2161.
586 C IF (MRCH.NE.0) CALL EDDYPA (UT,VT,TT,PT,QT,U,V,T,P,Q,DT1) 2161.5
587 MODEDY=MOD(MRCHT,ICHK) 2161.6
588 MODTVR=MOD(MRCHT,ICHK1)
589
590 #if ( defined CPL_CHEM )
591 !
592 ! --- Get time parameter for eddy calculation:
593 !
594 meddy1=0
595
596 IF (MODEDY.EQ.0.AND.MRCH.NE.0) THEN
597
598 meddy1=1
599 !
600 #else
601 IF (MODEDY.EQ.0.AND.MRCH.NE.0) THEN
602 #endif
603
604 EDDFLX=.false.
605 IF (MODEDY.eq.0)EDDFLX=.true.
606 EDDFLX=.true.
607 if(HPRNT)then
608 print *,' before eddypa'
609 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
610 print *,' JPR=',jpr,' LPR=',lpr
611 print *,' T=',T(1,JPR,LPR),' Q=',Q(1,JPR,LPR)
612 print *,' TT=',TT(1,JPR,LPR),' QT=',QT(1,JPR,LPR)
613 print *,' V(J,L)=',V(1,JPR,LPR)/FD(1,J)
614 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
615 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
616 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
617 print *,' U(13,1),U(13,2)'
618 print *,' U(J,L)=',U(1,13,1),' U(J,L)=',U(1,13,2)
619 endif
620 ! print *,' call eddypa TAU=',TAU
621 ! print *,' comp1 MRCHT=',MRCHT,' MRCH=',MRCH
622 c RFDU=2./(2.*P(1,JM)*DXYP(JM)+P(1,JM-1)*DXYP(JM-1))
623 c print *,(VT(1,J,1),J=1,Jm)
624 CALL EDDYPA (UT,VT,TT,PT,QT,U,V,T,P,Q,DTE,EDDFLX)
625 c print *,' after eddypa'
626 c print *,(VT(1,J,1),J=1,Jm)
627 if(HPRNT)then
628 print *,' after eddypa'
629 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
630 print *,' T=',T(1,JPR,LPR),' Q=',Q(1,JPR,LPR)
631 print *,' TT=',TT(1,JPR,LPR),' QT=',QT(1,JPR,LPR)
632 print *,' V(J,L)=',V(1,JPR,LPR)/FD(1,J)
633 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
634 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
635 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
636 endif
637 END IF
638 C
639 C 2162.
640 C UNDO SCALING FOR MOISTURE AND TRACE COMPOUNDS 2162.5
641 C 2163.
642 DO 1910 J=1,JM 2163.5
643 DO 1910 I=1,IM 2164.
644 1910 FD(I,J)=PT(I,J)*DXYP(J) 2164.5
645
646 #if ( defined CPL_CHEM )
647 !
648 ! --- get PAI for chemical advection (mass remapping):
649 !
650 if(mrch.ne.2)goto 6003
651 do i=1,n2dh
652 p11(i,1)=fd(i,1)
653 !
654 ! 092595
655 !
656 p4chem1(i,1) = pt(i,1)
657 enddo
658 6003 continue
659 !
660 #endif
661
662 IF(MRCH.EQ.0) GO TO 2000 2165.
663 DO 1940 L=1,NTP1LM 2165.5
664 DO 1928 J=1,JM 2166.
665 DO 1928 I=1,IM 2166.5
666 1928 QT(I,J,L)= QT(I,J,L)/FD(I,J) 2167.
667 1940 CONTINUE 2167.5
668 2000 CONTINUE 2168.
669 C 2168.5
670 C UNDO SCALING FOR TEMPERATURE AND WINDS 2169.
671 C 2169.5
672 DO 3510 L=1,LM 2170.
673 DO 3512 J=1,JM 2170.5
674 DO 3512 I=1,IM 2171.
675 3512 TT(I,J,L)=TT(I,J,L)/FD(I,J) 2171.5
676 DO 3513 J=1,JM 2172.
677 DO 3513 I=1,IX 2172.5
678 c IF(TT(I,J,L).LT.20..OR.TT(I,J,L).GT.120.) 2173.
679 IF(TT(I,J,L).LT.20..OR.TT(I,J,L).GT.130.) then
680 WRITE(6,902) I,J,L,MRCH,T(I,J,L),TT(I,J,L),P(I,J),PT(I,J) 2173.5
681 stop
682 endif
683 902 FORMAT(' POTENTIAL TEMPERATURE DIAGNOSTIC I,J,L,MRCH,T,TT,P,PT2174.
684 * =',4I4,4F8.1) 2174.5
685 3513 CONTINUE 2175.
686 3510 CONTINUE 2175.5
687 DO 3520 I=1,IM 2176.
688 FD(I,1)=2.*FD(I,1) 2176.5
689 3520 FD(I,JM)=2.*FD(I,JM) 2177.
690 DO 3542 J=2,JM 2177.5
691 DO 3542 I=1,IM 2178.
692 RFDU=2./(FD(I,J)+FD(I,J-1)) 2178.5
693 DO 3542 L=1,LMT2 2179.
694 3542 UT(I,J,L)=UT(I,J,L)*RFDU 2179.5
695 if(HPRNT)then
696 print *,' after 3542'
697 print *,' comp1 1 TAU=',TAU,' MRCH=',MRCH
698 print *,' T(J,L)=',T(1,JPR,LPR),' Q(J,L)=',Q(1,JPR,LPR)
699 print *,' TT(J,L)=',TT(1,JPR,LPR),' QT(J,L)=',QT(1,JPR,LPR)
700 print *,' V(J,L)=',V(1,JPR,LPR),' V(J,L)=',V(1,JPR,LPR)
701 print *,' VT(J,L)=',VT(1,JPR,LPR),' VT(J,L)=',VT(1,JPR,LPR)
702 print *,' U(J,L)=',U(1,JPR,LPR),' U(J,L)=',U(1,JPR,LPR)
703 print *,' UT(J,L)=',UT(1,JPR,LPR),' UT(J,L)=',UT(1,JPR,LPR)
704 endif
705 c print '(12f6.1/11f6.1)',(PT(1,j),j=1,jm)
706 RETURN 2180.
707 END 2180.5

  ViewVC Help
Powered by ViewVC 1.1.22