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

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

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


Revision 1.1 - (hide 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 jscott 1.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