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

Annotation of /MITgcm_contrib/jscott/igsm/src/md2g04.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:31 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     ! MD2G04.F: Lots of utility functions (some of which are unknown).
7     !
8     ! ----------------------------------------------------------
9     !
10     ! Revision History:
11     !
12     ! When Who What
13     ! ---- ---------- -------
14     ! 073100 Chien Wang repack based on CliChem3 & M24x11,
15     ! and add cpp.
16     ! 082900 Andrei some subroutines are taken out
17     !
18     ! ==========================================================
19    
20     SUBROUTINE CHECKT (N) 8742.
21     C**** 8743.
22     C**** THIS SUBROUTINE CHECKS WHETHER THE TEMPERATURES ARE REASONABLE 8744.
23     C**** FOR DEBUGGING PURPOSES. IT IS TURNED ON BY SETTING IDACC(11) 8745.
24     C**** TO BE POSITIVE. REMEMBER TO SET IDACC(11) BACK TO ZERO AFTER 8746.
25     C**** THE ERRORS ARE CORRECTED. 8747.
26     C**** 8748.
27    
28     #include "BD2G04.COM" 8749.
29    
30     COMMON U,V,T,P,Q 8750.
31     IF(IDACC(11).LE.0) RETURN 8751.
32     C**** 8752.
33     C**** CHECK WHETHER GDATA ARE REASONABLE AND CONSISTENT OVER EARTH 8753.
34     C**** 8754.
35     X=1.001 8755.
36     DO 110 J=1,JM 8756.
37     IMAX=IM 8757.
38     IF((J.EQ.1).OR.(J.EQ.JM)) IMAX=1 8758.
39     DO 110 I=1,IMAX 8759.
40     PEARTH=FDATA(I,J,2)*(1.-FDATA(I,J,3)) 8760.
41     IF(PEARTH.LE.0.) GO TO 110 8761.
42     IF(GDATA(I,J,2).GE.0..AND.GDATA(I,J,2)*GDATA(I,J,4).LE.0.)GO TO 508762.
43     WRITE (6,901) N,I,J,L,TAU,(GDATA(I,J,K),K=2,10) 8763.
44     STOP 9021 8764.
45     50 DO 70 L=1,2 8765.
46     TGL=GDATA(I,J,4*L) 8766.
47     WTRL=GDATA(I,J,4*L+1) 8767.
48     ACEL=GDATA(I,J,4*L+2) 8768.
49     IF((TGL+60.)*(60.-TGL).GT.0.) GO TO 60 8769.
50     WRITE (6,901) N,I,J,L,TAU,(GDATA(I,J,K),K=2,10) 8770.
51     60 IF(WTRL.GE.0..AND.ACEL.GE.0..AND.TGL*WTRL.GE.0..AND. 8771.
52     * TGL*ACEL.LE.0..AND.(WTRL+ACEL).LE.X*VDATA(I,J,8+L)) GO TO 70 8772.
53     WRITE (6,901) N,I,J,L,TAU,(GDATA(I,J,K),K=2,10) 8773.
54     STOP 9021 8774.
55     70 CONTINUE 8775.
56     110 CONTINUE 8776.
57     RETURN 8777.
58     C**** 8778.
59     901 FORMAT ('0GDATA UNREASONABLE, N,I,J,L,TAU=',4I4,E14.5/ 8779.
60     * ' GDATA(2-10)=',9F12.4) 8780.
61     END 8781.
62     c SUBROUTINE COSZR (IM,JM,SINL,COSL,SIND,COSD,ROT1,ROT2,COSZ,KR) 4795.5
63     SUBROUTINE COSZR (IM,JM,SIND,COSD,ROT1,ROT2,COSZ)
64     DIMENSION COSZ(IM,JM),SINL(46),COSL(46) 4795.51
65     COMMON/SCJL/SINL,COSL
66     DATA TWOPI/6.283185/,ZERO/0./ 4795.52
67     KR=0
68     DO 200 J=1,JM 4795.53
69     IF(SINL(J)*SIND-COSL(J)*COSD.GE.-ZERO) 4795.54
70     * GO TO 230 4795.55
71     IF(SINL(J)*SIND+COSL(J)*COSD.LE.ZERO) GO TO 250 4795.56
72     HALFD=ABS(ACOS(-SINL(J)*SIND/(COSL(J)*COSD))) 4795.57
73     HALFS=HALFD*86400./TWOPI 4795.58
74     SINHD=SIN(HALFD) 4795.59
75     TEMP=TWOPI*HALFS/86400. 4795.6
76     COSZ(1,J)=SINL(J)*SIND+COSL(J)*COSD*SINHD/TEMP 4795.61
77     IF(KR.EQ.1) COSZ(1,J)=COSZ(1,J)*2.*HALFS/86400. 4795.62
78     GO TO 200 4795.63
79     230 COSZ(1,J)=SINL(J)*SIND 4795.64
80     GO TO 200 4795.65
81     250 COSZ(1,J)=ZERO 4795.66
82     200 CONTINUE 4795.67
83     RETURN 4795.68
84     END 4795.69
85     SUBROUTINE FILTER 8001.
86     C********* 8526.
87     C MFILTR=0 NO SHAPIRO FILTER 8527.
88     C MFILTR=1 SHAPIRO FILTER 1D 2TH ORDER ON SEA LEVEL PRESSURE 8528.
89     C MFILTR=2 SHAPIRO FILTER 1D 2TH ORDER ON SLP & THETA 8529.
90     C MFILTR=3 SHAPIRO FILTER 1D 2TH ORDER ON U-WIND 8530.
91     C********* 8531.
92    
93     #include "BD2G04.COM" 8532.
94    
95     COMMON U,V,T,P,Q 8533.
96     COMMON/WORK2/X1JI(72,46),X2JI(72,46),X3JI(72,46),X1(72),X2(72), 8534.
97     * X3(72),X4(72) 8535.
98     LOGICAL HPRNT
99     DATA IFIRST/1/ 8536.
100     IF ( IFIRST.NE.1 ) GO TO 100 8537.
101     IFIRST=0 8538.
102     C**** INITIALIZE CERTAIN QUANTITIES 8539.
103     CUTPLD=0.5 8540.
104     BETAG=.0065/GRAV 8541.
105     GBYRB=1./(RGAS*BETAG) 8542.
106     AKAP=KAPA-.205 8543.
107     IMM1=IM-1 8544.
108     IMBY2=IM/2 8545.
109     NORDER=16 8546.
110     c NORDER=8
111     XXX=4.**NORDER 8547.
112     WRITE(6,4000) MFILTR,NORDER,XXX 8548.
113     4000 FORMAT(' PERFORM SHAPIRO FILTER MFILTR=',I2,' NORDER=',I2, 8549.
114     * ' XXX=',F15.0) 8550.
115     print *,' filter for ps, t, and q'
116     c print *,' filter for q only'
117     100 CONTINUE 8551.
118     HPRNT=.FALSE.
119     C
120     c go to 250
121     C
122     C**** PERFORM FILTER ON SEA LEVEL PRESSURE 8552.
123     c IF (MFILTR.NE.1) GO TO 200 8553.
124     DO 110 J=1,JM 8554.
125     DO 110 I=1,IM 8555.
126     X2JI(I,J)=(1.+BETAG*FDATA(I,J,1)/BLDATA(I,J,2))**GBYRB 8556.
127     X1JI(I,J)=P(I,J)*X2JI(I,J) 8557.
128     X3JI(I,J)=X1JI(I,J) 8558.
129     110 CONTINUE 8559.
130     CALL SHAP2D (MFILTR,NORDER,XXX,IM,JM,1,1) 8560.
131     DO 44 I=1,IM 8561.
132     DO 44 J=1,JM 8562.
133     P(I,J)=X1JI(I,J)/X2JI(I,J) 8563.
134     44 CONTINUE 8564.
135     DOPK=1. 8565.
136     200 continue
137     c 200 IF(MFILTR.NE.2) GO TO 300 8566.
138     C**** PERFORM FILTER ON POTENTIAL TEMPERATURE 8567.
139     DO 240 L=1,LM 8568.
140     DO 210 J=1,JM 8569.
141     DO 210 I=1,IM 8570.
142     X2JI(I,J)=1. 8571.
143     X1JI(I,J)=T(I,J,L)*X2JI(I,J) 8572.
144     X3JI(I,J)=X1JI(I,J) 8573.
145     210 CONTINUE 8574.
146     CALL SHAP2D (MFILTR,NORDER,XXX,IM,JM,1,1) 8575.
147     DO 240 J=1,JM 8576.
148     DO 240 I=1,IM 8577.
149     240 T(I,J,L)=X1JI(I,J)/X2JI(I,J) 8578.
150     250 continue
151     C**** PERFORM FILTER ON SPECIFIC HUMIDITY 8578.01
152     DO 280 L=1,LM 8578.02
153     DO 260 J=1,JM 8578.03
154     DO 260 I=1,IM 8578.04
155     X2JI(I,J)=1. 8578.05
156     X1JI(I,J)=Q(I,J,L)*X2JI(I,J) 8578.06
157     X3JI(I,J)=X1JI(I,J) 8578.07
158     260 CONTINUE 8578.08
159     CALL SHAP2D (MFILTR,NORDER,XXX,IM,JM,1,1) 8578.09
160     DO 280 J=1,JM 8578.1
161     DO 280 I=1,IM 8578.11
162     QTEMP=X1JI(I,J)/X2JI(I,J) 8578.12
163     IF(QTEMP.LE.0.) QTEMP=0. 8578.13
164     280 Q(I,J,L)=QTEMP 8578.15
165     300 IF (MFILTR.NE.3) RETURN 8579.
166     print *,' WRONG FILTER'
167     C**** PERFORM FILTER ON U-WIND 8580.
168     DO 340 L=1,LM 8581.
169     DO 310 J=2,JM 8582.
170     X1JI(1,J)=U(1,J,L) 8583.
171     310 X3JI(1,J)=U(1,J,L) 8584.
172     CALL SHAP2D (MFILTR,NORDER,XXX,IM,JM,2,3) 8585.
173     DO 340 J=2,JM 8586.
174     340 U(1,J,L)=X1JI(1,J) 8587.
175     C**** PERFORM FILTER ON V-WIND
176     if(HPRNT)then
177     print *,' from filter 1'
178     print *,' V(J,3)'
179     print *,(V(1,J,3),J=1,JM)
180     endif
181     DO 440 L=1,LM
182     DO 410 J=2,JM
183     X1JI(1,J)=V(1,J,L)
184     410 X3JI(1,J)=V(1,J,L)
185     CALL SHAP2D (MFILTR,NORDER,XXX,IM,JM,2,3)
186     DO 440 J=2,JM
187     440 V(1,J,L)=X1JI(1,J)
188     if(HPRNT)then
189     print *,' from filter 2'
190     print *,' V(J,3)'
191     print *,(V(1,J,3),J=1,JM)
192     endif
193     RETURN 8588.
194     END 8589.
195     SUBROUTINE DTDX1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 9112.4
196     DIMENSION XNOW(5),XREF(5),XDT0(5),KFOR(5) 9112.5
197     C FCO2(X)=DLOG(1.0+1.3365*X+9.29E-03*X*X+2.1766E-06*X**3) 9112.6
198     C FFNC(X,Y)=1.54*DLOG(1.0+1.11*X**0.77) - 0.068*DLOG(1.0+X*Y) 9112.7
199     C + +0.231*Y**0.65/(1.0+0.113*Y**0.71) 9112.8
200     C 9112.9
201     FCO2(X)=DLOG(1.0+0.942*X/(1.0+6.2E-04*X)+8.8E-03*X*X 9113.
202     + +3.26E-06*X**3+0.156*X**1.3*EXP(-X/760.0)) 9113.1
203     C 9113.2
204     FFNC(X,Y)=1.556*DLOG(1.0+1.098*X**0.77*(1.0+0.032*X) 9113.3
205     + /(1.0+0.0014*X*X)) 9113.4
206     + +(0.394*Y**0.66+0.16*Y*EXP(-1.6*Y)) 9113.5
207     + /(1.0+0.169*Y**0.62) 9113.6
208     + -0.14*DLOG(1.0+0.636*(X*Y)**0.75+0.007*Y*(X*Y)**1.52) 9113.7
209     C 9113.8
210     SUM=0.0 9113.9
211     DO 200 K=1,NFOR 9114.
212     XDT0(K)=0.0 9114.1
213     XX=XNOW(K) 9114.2
214     XR=XREF(K) 9114.3
215     IF(KFOR(K).EQ.0) GO TO 200 9114.4
216     GO TO (110,120,130,140,150),K 9114.5
217     GO TO 200 9114.6
218     110 XDTX=FCO2(XX) 9114.7
219     XDTR=FCO2(XR) 9114.8
220     GO TO 190 9114.9
221     120 YY=XNOW(3) 9115.
222     YR=XREF(3) 9115.1
223     XDTX=FFNC(XX,YR) 9115.2
224     XDTR=FFNC(XR,YR) 9115.3
225     GO TO 190 9115.4
226     130 XR=XREF(2) 9115.5
227     YY=XNOW(3) 9115.6
228     YR=XREF(3) 9115.7
229     XDTX=FFNC(XR,YY) 9115.8
230     XDTR=FFNC(XR,YR) 9115.9
231     GO TO 190 9116.
232     140 XX=XNOW(4) 9116.1
233     XR=XREF(4) 9116.2
234     XDTX=0.066*XX 9116.3
235     XDTR=0.066*XR 9116.4
236     GO TO 190 9116.5
237     150 XX=XNOW(5) 9116.6
238     XR=XREF(5) 9116.7
239     XDTX=0.084*XX 9116.8
240     XDTR=0.084*XR 9116.9
241     GO TO 190 9117.
242     190 XDT0(K)=XDTX-XDTR 9117.1
243     SUM=SUM+XDT0(K) 9117.2
244     200 CONTINUE 9117.3
245     SDT0=SUM 9117.4
246     RETURN 9117.5
247     END 9117.6
248     SUBROUTINE DTDX3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 9126.5
249     DIMENSION XNOW(5),XREF(5),XDT0(5),KFOR(5) 9126.6
250     C FCO2(X)=DLOG(1.0 + 0.71*X**1.1 + 0.0084*X**1.9 + 2.8E-07*X**3.3) 9126.7
251     FCO2(X)=DLOG(1.0 + 1.2*X + 0.005*X**2 + 1.4E-06*X**3) 9126.8
252     FFNC(X,Y)=3.1515*X/(1.0+0.21448*X**0.87) - 0.056*X*Y 9126.9
253     + +0.3234*Y/(1.0+0.052*Y**0.84) 9127.
254     C 9127.1
255     SUM=0.0 9127.2
256     DO 200 K=1,NFOR 9127.3
257     XDT0(K)=0.0 9127.4
258     IF(KFOR(K).EQ.0) GO TO 200 9127.5
259     GO TO (110,120,130,140,150),K 9127.6
260     GO TO 200 9127.7
261     110 XX=XNOW(1) 9127.8
262     XR=XREF(1) 9127.9
263     XDTX=FCO2(XX) 9128.
264     XDTR=FCO2(XR) 9128.1
265     GO TO 190 9128.2
266     120 XX=XNOW(2) 9128.3
267     XR=XREF(2) 9128.4
268     YR=XREF(3) 9128.5
269     XDTX=FFNC(XX,YR) 9128.6
270     XDTR=FFNC(XR,YR) 9128.7
271     GO TO 190 9128.8
272     130 XR=XREF(2) 9128.9
273     YY=XNOW(3) 9129.
274     YR=XREF(3) 9129.1
275     XDTX=FFNC(XR,YY) 9129.2
276     XDTR=FFNC(XR,YR) 9129.3
277     GO TO 190 9129.4
278     140 XX=XNOW(4) 9129.5
279     XR=XREF(4) 9129.6
280     XDTX=0.066*XX 9129.7
281     XDTR=0.066*XR 9129.8
282     GO TO 190 9129.9
283     150 XX=XNOW(5) 9130.
284     XR=XREF(5) 9130.1
285     XDTX=0.084*XX 9130.2
286     XDTR=0.084*XR 9130.3
287     GO TO 190 9130.4
288     190 XDT0(K)=XDTX-XDTR 9130.5
289     SUM=SUM+XDT0(K) 9130.6
290     200 CONTINUE 9130.7
291     SDT0=SUM 9130.8
292     RETURN 9130.9
293     END 9131.
294     SUBROUTINE DXDT1D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 9117.7
295     DIMENSION XNOW(5),XREF(5),XDT0(5),KFOR(5) 9117.8
296     C FCO2(X)=DLOG(1.0+1.3365*X+9.29E-03*X*X+2.1766E-06*X**3) 9117.9
297     C FFNC(X,Y)=1.54*DLOG(1.0+1.11*X**0.77) - 0.068*DLOG(1.0+X*Y) 9118.
298     C + +0.231*Y**0.65/(1.0+0.113*Y**0.71) 9118.1
299     C 9118.2
300     FCO2(X)=DLOG(1.0+0.942*X/(1.0+6.2E-04*X)+8.8E-03*X*X 9118.3
301     + +3.26E-06*X**3+0.156*X**1.3*EXP(-X/760.0)) 9118.4
302     C 9118.5
303     FFNC(X,Y)=1.556*DLOG(1.0+1.098*X**0.77*(1.0+0.032*X) 9118.6
304     + /(1.0+0.0014*X*X)) 9118.7
305     + +(0.394*Y**0.66+0.16*Y*EXP(-1.6*Y)) 9118.8
306     + /(1.0+0.169*Y**0.62) 9118.9
307     + -0.14*DLOG(1.0+0.636*(X*Y)**0.75+0.007*Y*(X*Y)**1.52) 9119.
308     C 9119.1
309     DO 200 K=1,NFOR 9119.2
310     XN=XNOW(K) 9119.3
311     XR=XREF(K) 9119.4
312     ITER=0 9119.5
313     DTK=XDT0(K) 9119.6
314     IF(XN.EQ.XR) GO TO 200 9119.7
315     IF(KFOR(K).EQ.0) GO TO 200 9119.8
316     GO TO (110,120,130,140,150),K 9119.9
317     GO TO 200 9120.
318     110 TR=FCO2(XR) 9120.1
319     111 TN=FCO2(XN) 9120.2
320     DTNR=TN-TR 9120.3
321     ITER=ITER+1 9120.4
322     SLOPE=DTNR/(XN-XR) 9120.5
323     DT=DTK-DTNR 9120.6
324     DX=DT/SLOPE 9120.7
325     XN=XN+DX 9120.8
326     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 111 9120.9
327     GO TO 190 9121.
328     120 YR=XREF(3) 9121.1
329     TR=FFNC(XR,YR) 9121.2
330     121 TN=FFNC(XN,YR) 9121.3
331     DTNR=TN-TR 9121.4
332     ITER=ITER+1 9121.5
333     SLOPE=DTNR/(XN-XR) 9121.6
334     DT=DTK-DTNR 9121.7
335     DX=DT/SLOPE 9121.8
336     XN=XN+DX 9121.9
337     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 121 9122.
338     GO TO 190 9122.1
339     130 XR=XREF(2) 9122.2
340     YN=XNOW(3) 9122.3
341     YR=XREF(3) 9122.4
342     TR=FFNC(XR,YR) 9122.5
343     131 TN=FFNC(XR,YN) 9122.6
344     DTNR=TN-TR 9122.7
345     ITER=ITER+1 9122.8
346     SLOPE=DTNR/(YN-YR) 9122.9
347     DT=DTK-DTNR 9123.
348     DY=DT/SLOPE 9123.1
349     YN=YN+DY 9123.2
350     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 131 9123.3
351     XN=YN 9123.4
352     GO TO 190 9123.5
353     140 XN=XNOW(4) 9123.6
354     XR=XREF(4) 9123.7
355     141 TN=0.066*XN 9123.8
356     TR=0.066*XR 9123.9
357     DTNR=TN-TR 9124.
358     ITER=ITER+1 9124.1
359     SLOPE=DTNR/(XN-XR) 9124.2
360     DT=DTK-DTNR 9124.3
361     DX=DT/SLOPE 9124.4
362     XN=XN+DX 9124.5
363     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 141 9124.6
364     TR=0.066*XR 9124.7
365     GO TO 190 9124.8
366     150 XN=XNOW(5) 9124.9
367     XR=XREF(5) 9125.
368     TR=0.084*XR 9125.1
369     151 TN=0.084*XN 9125.2
370     DTNR=TN-TR 9125.3
371     ITER=ITER+1 9125.4
372     SLOPE=DTNR/(XN-XR) 9125.5
373     DT=DTK-DTNR 9125.6
374     DX=DT/SLOPE 9125.7
375     XN=XN+DX 9125.8
376     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 151 9125.9
377     GO TO 190 9126.
378     190 XNOW(K)=XN 9126.1
379     200 CONTINUE 9126.2
380     RETURN 9126.3
381     END 9126.4
382     SUBROUTINE DXDT3D(XNOW,XREF,XDT0,SDT0,KFOR,NFOR) 9131.1
383     DIMENSION XNOW(5),XREF(5),XDT0(5),KFOR(5) 9131.2
384     C FCO2(X)=DLOG(1.0 + 0.71*X**1.1 + 0.0084*X**1.9 + 2.8E-07*X**3.3) 9131.3
385     FCO2(X)=DLOG(1.0 + 1.2*X + 0.005*X**2 + 1.4E-06*X**3) 9131.4
386     FFNC(X,Y)=3.1515*X/(1.0+0.21448*X**0.87) - 0.056*X*Y 9131.5
387     + +0.3234*Y/(1.0+0.052*Y**0.84) 9131.6
388     C 9131.7
389     DO 200 K=1,NFOR 9131.8
390     XN=XNOW(K) 9131.9
391     XR=XREF(K) 9132.
392     ITER=0 9132.1
393     DTK=XDT0(K) 9132.2
394     IF(XN.EQ.XR) GO TO 200 9132.3
395     IF(KFOR(K).EQ.0) GO TO 200 9132.4
396     GO TO (110,120,130,140,150),K 9132.5
397     GO TO 200 9132.6
398     110 TR=FCO2(XR) 9132.7
399     111 TN=FCO2(XN) 9132.8
400     DTNR=TN-TR 9132.9
401     ITER=ITER+1 9133.
402     SLOPE=DTNR/(XN-XR) 9133.1
403     DT=DTK-DTNR 9133.2
404     DX=DT/SLOPE 9133.3
405     XN=XN+DX 9133.4
406     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 111 9133.5
407     GO TO 190 9133.6
408     120 YR=XREF(3) 9133.7
409     TR=FFNC(XR,YR) 9133.8
410     121 TN=FFNC(XN,YR) 9133.9
411     DTNR=TN-TR 9134.
412     ITER=ITER+1 9134.1
413     SLOPE=DTNR/(XN-XR) 9134.2
414     DT=DTK-DTNR 9134.3
415     DX=DT/SLOPE 9134.4
416     XN=XN+DX 9134.5
417     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 121 9134.6
418     GO TO 190 9134.7
419     130 XR=XREF(2) 9134.8
420     YN=XNOW(3) 9134.9
421     YR=XREF(3) 9135.
422     TR=FFNC(XR,YR) 9135.1
423     131 TN=FFNC(XR,YN) 9135.2
424     DTNR=TN-TR 9135.3
425     ITER=ITER+1 9135.4
426     SLOPE=DTNR/(YN-YR) 9135.5
427     DT=DTK-DTNR 9135.6
428     DY=DT/SLOPE 9135.7
429     YN=YN+DY 9135.8
430     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 131 9135.9
431     XN=YN 9136.
432     GO TO 190 9136.1
433     140 XN=XNOW(4) 9136.2
434     XR=XREF(4) 9136.3
435     141 TN=0.066*XN 9136.4
436     TR=0.066*XR 9136.5
437     DTNR=TN-TR 9136.6
438     ITER=ITER+1 9136.7
439     SLOPE=DTNR/(XN-XR) 9136.8
440     DT=DTK-DTNR 9136.9
441     DX=DT/SLOPE 9137.
442     XN=XN+DX 9137.1
443     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 141 9137.2
444     TR=0.066*XR 9137.3
445     GO TO 190 9137.4
446     150 XN=XNOW(5) 9137.5
447     XR=XREF(5) 9137.6
448     TR=0.084*XR 9137.7
449     151 TN=0.084*XN 9137.8
450     DTNR=TN-TR 9137.9
451     ITER=ITER+1 9138.
452     SLOPE=DTNR/(XN-XR) 9138.1
453     DT=DTK-DTNR 9138.2
454     DX=DT/SLOPE 9138.3
455     XN=XN+DX 9138.4
456     IF(ITER.LT.15.AND.ABS(DT).GT.1.E-05) GO TO 151 9138.5
457     GO TO 190 9138.6
458     190 XNOW(K)=XN 9138.7
459     200 CONTINUE 9138.8
460     RETURN 9138.9
461     END 9139.
462    
463     FUNCTION FUNC1(KX,KY,KZ) 2666.5
464     CONST=0.5 2667.
465     KA=KX-1 2667.5
466     KB=KY-1 2668.
467     KC=KZ-1 2668.5
468     K1=KA+KB 2669.
469     K2=KA+KC 2669.5
470     K3=KB+KC 2670.
471     FUNC1=0. 2670.5
472     IF((KA.EQ.K3).OR.(KB.EQ.K2).OR.(KC.EQ.K1)) FUNC1=CONST 2671.
473     RETURN 2671.5
474     C 2672.
475     ENTRY FUNC2(KX,KY,KZ) 2672.5
476     CONST=0.5 2673.
477     KA=KX-1 2673.5
478     KB=KY-1 2674.
479     KC=KZ-1 2674.5
480     K1=KA+KB 2675.
481     K2=KA+KC 2675.5
482     K3=KB+KC 2676.
483     FUNC2=0. 2676.5
484     IF((KA.GT.KB).AND.(KA.GT.KC).AND.(KA.EQ.K3)) FUNC2 =-CONST 2677.
485     IF((KA.LT.KB).AND.(KB.EQ.K2)) FUNC2=CONST 2677.5
486     IF((KA.LT.KC).AND.(KC.EQ.K1)) FUNC2=CONST 2678.
487     RETURN 2678.5
488     END 2679.
489    
490     c FUNCTION RANDU (X) 3801.
491     SUBROUTINE RANDUU(RX,X)
492     integer IX,IY
493     C**** THIS FUNCTION GENERATES RANDOM NUMBERS ON AN IBM 360 OR 370 3802.
494     10 IY=IX*65539 3803.
495     IF(IY) 20,40,30 3804.
496     20 IY=(IY+2147483647)+1 3805.
497     30 IX=IY 3806.
498     c RANDU=FLOAT(IY)*.4656613E-9 3807.
499     RX=FLOAT(IY)*.4656613E-9
500     RETURN 3808.
501     40 IX=1 3809.
502     GO TO 10 3810.
503     ENTRY RINIT (INIT) 3811.
504     IX=INIT 3812.
505     RETURN 3813.
506     ENTRY RFINAL (IFINAL) 3814.
507     IFINAL=IX 3815.
508     RETURN 3816.
509     END 3817.
510     SUBROUTINE SHAP2D (MFILTR,NORDER,XXX,IM,JM,J1,ITYPE) 8590.
511     COMMON/WORK2/X1JI(72,46),X2JI(72,46),X3JI(72,46),X1(72),X2(72), 8591.
512     * X3(72),X4(72),XM1(72),XJMP1(72) 8592.
513     C VARIABLE ITYPE DETERMINES TYPE OF BOUNDARY CONDITIONS
514     C ITYPE=1 FOR PS,T AND Q ( XM1=X2)
515     C ITYPE=2 FOR U (XM1=X1)
516     C ITYPE=3 FOR V (XM1=-X1)
517     JMM1=JM-1 8593.
518     J2=J1+1 8594.
519     IMBY2=1 8595.
520     DO 145 N=1,NORDER 8596.
521     DO 146 K=1,IM 8597.
522     X1(K)=X1JI(K,J1) 8598.
523     X2(K)=X1JI(K,J2) 8599.
524     X3(K)=X1JI(K,JMM1) 8600.
525     X4(K)=X1JI(K,JM) 8601.
526     IF(ITYPE.EQ.1)THEN
527     XM1(K)=X1JI(K,J2)
528     XJMP1(K)=X1JI(K,JMM1)
529     ELSEIF(ITYPE.EQ.2)THEN
530     XM1(K)=X1JI(K,J1)
531     XJMP1(K)=X1JI(K,JM)
532     ELSE
533     XM1(K)=-X1JI(K,J1)
534     XJMP1(K)=-X1JI(K,JM)
535     ENDIF
536     146 CONTINUE 8602.
537     DO 142 I=1,IM 8603.
538     X1IM1=X1JI(I,J1) 8604.
539     DO 142 J=J2,JMM1 8605.
540     X1I=X1JI(I,J) 8606.
541     X1JI(I,J)=X1IM1-X1I-X1I+X1JI(I,J+1) 8607.
542     X1IM1=X1I 8608.
543     142 CONTINUE 8609.
544     SUM1=0. 8610.
545     SUMJM=0. 8611.
546     DO 144 K=1,IMBY2 8612.
547     ccc SUM1 =SUM1 +X2(K)-X1(K)-X1(K)+X2(K) 8613.
548     SUM1 =SUM1 +XM1(K)-X1(K)-X1(K)+X2(K)
549     ccc SUMJM=SUMJM+X3(K)-X4(K)-X4(K)+X3(K) 8614.
550     SUMJM=SUMJM+X3(K)-X4(K)-X4(K)+XJMP1(K)
551     144 CONTINUE 8615.
552     X1SUM =SUM1 /IMBY2 8616.
553     XJMSUM =SUMJM/IMBY2 8617.
554     DO 147 K=1,IM 8618.
555     X1JI(K,JM)=XJMSUM 8619.
556     147 X1JI(K,J1)= X1SUM 8620.
557     145 CONTINUE 8621.
558     DO 160 I=1,IM 8622.
559     DO 160 J=J1,JM 8623.
560     X1JI(I,J)=(X3JI(I,J)-X1JI(I,J)/XXX) 8624.
561     160 CONTINUE 8625.
562     RETURN 8626.
563     END 8627.

  ViewVC Help
Powered by ViewVC 1.1.22