| 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. |
Parent Directory
| 
Revision Log
| 
Revision Graph