IMPLICIT REAL*8 (A-H,O-Z) LOGICAL*1 LOG(12,960),LAST,LBLK DIMENSION ADATA(12,960),BDATA(12,960),XS(12),YS(12),DATA(12), @X(12),Y(12),Z(12),ALAT(12),ALON(12),P(6,3),R(12,3),S(3),Q(6,2), @ SNAM(12), @KUSE(3),IER(4),KODE(2),IOS(12) DATA CONV,DENS,ROT/3.386389D0,1.2255D0,7.292115D-5/, @ KP,KY,KA,KBLK,KIND,KODE/'P','Y','A',' ','->','WN','WE'/, @ LAST,LBLK/'*','.'/, MS/12/, KF/'F'/ C C####################################################################### C RHOD=3.141592653589793D0/180.D0 READ (5,501) LOUT,IYRE,IYRS,IYRF,IPUN,IRED,IFOR,A,F,IOS,LIS 501 FORMAT (I2,3I5,2A1,I1,2F10.2,12I1,A1) DO 400 I=1,MS DO 400 J=1,960 LOG(I,J)=LBLK ADATA(I,J)=0.D0 400 CONTINUE JS=1 IF (IYRS.NE.0) JS=(IYRS-IYRE)*12+1 TRO=2.D0*DENS*ROT MJ=0 NS=0 102 READ (5,502,END=100) KOD,ANAM,IYR,DATA 502 FORMAT (A1,A6,I4,12F5.2) IF (KOD.EQ.KP.OR.KOD.EQ.KF) GO TO 101 NS=NS+1 IF (NS.GT.MS) GO TO 103 SNAM(NS)=ANAM ALAT(NS)=DATA(1)+40.D0 ALON(NS)=-DATA(2)-60.D0 IF (DATA(2).LT.2.D0) ALON(NS)=ALON(NS)-10.D0 GO TO 102 101 IOFF=(IYR-IYRE)*12 DO 401 I=1,12 IF (DATA(I).LE.0.D0) GO TO 401 J=IOFF+I IF (J.LT.1.OR.J.GT.960) GO TO 104 ADATA(NS,J)=DATA(I)*CONV MJ=MAX0(J,MJ) LOG(NS,J)=LAST 401 CONTINUE GO TO 102 103 WRITE (6,601) MS 601 FORMAT ('1*** ERROR *** TOO MANY STATIONS (MAX =',I3,') -- PROCES @SING ABANDONED') STOP 1111 104 WRITE (6,602) NS,ALAT(NS),ALON(NS),IYR,IYRE 602 FORMAT ('1*** ERROR *** DATA EXCEEDS ARRAY SIZE'/'0STN. NO.',I2, @' COORDS:',2F10.3,6X,'YEAR V',I5,' BEGIN YEAR V',I5) STOP 2222 100 IF (NS.GT.0) GO TO 115 WRITE (6,607) 607 FORMAT ('1*** ERROR *** NO INPUT DATA FOUND') STOP 3333 115 IF (LIS.EQ.KY) WRITE (6,619) ((ADATA(I,J),I=1,12),J=1,960) 619 FORMAT (1X,12F7.2) CALL PLANE (A,F,ALAT,ALON,NS,XS,YS,ALATO,ALONO) DO 402 I=1,NS XS(I)=XS(I)/1.D3 YS(I)=YS(I)/1.D3 402 CONTINUE WRITE (6,603) LOUT,IYRE,IYRS,IYRF,MJ,NS,IPUN,A,F,(I,SNAM(I), @ ALAT(I),ALON(I),XS(I),YS(I),I=1,NS) 603 FORMAT ('1',10X,17('*')/11X,'*',15X,'*'/11X,'* G E O W I N D *'/ @ 11X,'*',15X,'*'/11X,17('*')/////11X,'PARAMETERS:'/11X,11('=')/// @ 11X,'OUTPUT UNIT =',I5//11X,'BEGIN YEAR =',I5// @ 11X,'START/STOP =',2I5//11X,'MAX. MONTHS =',I5// @ 11X,'NO. STNS. =',I5//11X,'PUNCH O/P = ',A1// @ 11X,'ELLIPSOID = ',2F10.2////11X,'STATION COORDINATES:'/ @ 11X,20('=')///15X,53('-')/15X,'NO STN.',7X,'LAT',7X,'LON',9X, @ 'X (KM)',5X,'Y (KM)'/15X,53('-')/(15X,I2,1X,A6,2F10.3,2X,2F11.3)) WRITE (6,604) 604 FORMAT (15X,53('-')) IF (IPUN.EQ.KY) WRITE (7,703) (I,ALAT(I),ALON(I),XS(I),YS(I), @ I=1,NS) 703 FORMAT ('*',4X,I2,1X,4F12.3) DO 405 IB=1,8 IF (MOD(IB,3).EQ.1) WRITE (6,605) 605 FORMAT ('1AVAILABLE PRESSURE DATA'/1X,23('=')) IYS=IYRE+(IB-1)*10 IYF=IYS+9 JB=(IB-1)*120+1 JF=JB+119 WRITE (6,606) (IY,IY=IYS,IYF),(I,(LOG(I,J),J=JB,JF),I=1,MS) 606 FORMAT (///'0',I6,9I12/5X,'|',9(11X,'|')/(1X,I2,2X,120A1)) 405 CONTINUE IF (IYRF.NE.0) MJ=(IYRF-IYRE)*12+12 DO 403 J=JS,MJ IY=IYRE+(J-1)/12 MON=MOD(J,12) IF (MON.EQ.0) MON=12 IF (MOD(MON,6).EQ.1) WRITE (6,608) 608 FORMAT ('1') N=0 DO 404 I=1,NS IF (ADATA(I,J).LE.0.D0) GO TO 404 N=N+1 X(N)=XS(I) Y(N)=YS(I) Z(N)=ADATA(I,J) 404 CONTINUE IF (N.GT.2) GO TO 105 WRITE (6,609) IY,MON 609 FORMAT (1X,I4,I3,6X,'*** INSUFFICIENT DATA TO FIT A SURFACE ***' @ /////////) MIN=0 GO TO 111 105 IERR=IPRBLD (X,Y,Z,N,P,R,MS,S,Q) IF (IFOR.LE.0) GOTO 116 MIN=(IFOR-1)/1.6 GO TO 117 116 MIN=0 SM=1.D50 DO 406 I=1,3 NP=I*1.6+2 KUSE(I)=KBLK IT=IERR/10 IER(I)=IERR-IT*10 IERR=IT IF (IER(I).EQ.1.OR.S(I).GT.SM) GO TO 406 IF (IRED.EQ.KY.AND.NP.EQ.N) GO TO 406 SM=S(I) MIN=I 406 CONTINUE 117 WRITE (6,610) IY,MON,N 610 FORMAT (1X,I4,'/',I2,I4) IF (MIN.GT.0) KUSE(MIN)=KIND DO 407 I=1,3 IF (IER(I).EQ.1) GO TO 106 NP=I*1.6+2 WRITE (6,611) KUSE(I),(P(IP,I),IP=1,NP) 611 FORMAT (8X,A2,6D15.6) GO TO 407 106 WRITE (6,612) 612 FORMAT (14X,'*** NO SOLUTION ***') 407 CONTINUE WRITE (6,613) 613 FORMAT (1X) DO 408 I=1,3 IF (IER(I).EQ.1) GO TO 107 WRITE (6,614) KUSE(I),S(I),(R(IP,I),IP=1,N) 614 FORMAT (8X,A2,D15.6,(12F9.5)) GO TO 408 107 WRITE (6,612) 408 CONTINUE WRITE (6,615) 615 FORMAT ('0') NP=MIN*1.6+2 111 IF (IPUN.NE.KY) GO TO 108 IF (MIN.EQ.0) GO TO 109 WRITE (7,701) KP,IY,MON,NP,(P(IP,MIN),IP=1,NP) 701 FORMAT (A1,I4,I2,I1,1P6D12.5) IF (MIN.EQ.1) GO TO 110 WRITE (7,701) KA,IY,MON,NP,(Q(IP,MIN-1),IP=1,NP) GO TO 108 110 WRITE (7,701) KA,IY,MON,NP GO TO 108 109 WRITE (7,701) KP,IY,MON,MIN WRITE (7,701) KA,IY,MON,MIN 108 IF (MIN.NE.0) GO TO 112 DO 409 I=1,NS ADATA(I,J)=0.D0 BDATA(I,J)=0.D0 409 CONTINUE GO TO 403 112 DO 410 I=1,NS TROS=TRO*DSIN(ALAT(I)*RHOD) ADATA(I,J)=(P(2,MIN)+2.D0*P(4,MIN)*XS(I)+P(6,MIN)*YS(I))/TROS BDATA(I,J)=-(P(3,MIN)+2.D0*P(5,MIN)*YS(I)+P(6,MIN)*XS(I))/TROS 410 CONTINUE 403 CONTINUE IYRS=IYRE+(JS-1)/12 IYRF=IYRE+(MJ-1)/12 WRITE (6,616) IYRS,IYRF 616 FORMAT ('1GEOSTROPHIC WIND SPEED COMPONENTS'/1X,33('=')/'-PERIOD:' @,I5,' TO',I5/'-*** LISTING OF CARD IMAGES FOLLOWS ***') DO 412 I=1,NS IF (IOS(I).NE.0) GO TO 412 WRITE (6,617) I,ALAT(I),ALON(I) 617 FORMAT ('1STATION NO.',I2,5X,'(',F5.2,'/',F6.2,')'/1X,32('-')) DO 411 K=1,2 DO 411 IY=IYRS,IYRF JS=(IY-IYRE)*12+1 JF=JS+11 GO TO (113,114),K 113 WRITE (6,618) KODE(K),I,IY,(ADATA(I,J),J=JS,JF) 618 FORMAT (1X,A2,I2,I4,12F6.2) IF (LOUT.GE.7) WRITE (LOUT,702) KODE(K),I,IY,(ADATA(I,J),J=JS,JF) 702 FORMAT (A2,I2,I4,12F6.2) GO TO 411 114 WRITE (6,618) KODE(K),I,IY,(BDATA(I,J),J=JS,JF) IF (LOUT.GE.7) WRITE (LOUT,702) KODE(K),I,IY,(BDATA(I,J),J=JS,JF) 411 CONTINUE 412 CONTINUE STOP END FUNCTION IPRBLD (X,Y,Z,N,P,R,NR,S,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),Y(N),Z(N),R(NR,3),P(6,3),S(3),A(6,2),V(6,6),W(6) DO 400 I=1,3 S(I)=-1.D0 DO 400 J=1,6 P(J,I)=0.D0 400 CONTINUE DO 401 I=1,6 W(I)=0.D0 DO 401 J=1,6 V(I,J)=0.D0 401 CONTINUE IF (N.GE.3) GO TO 101 IPRBLD=-1 RETURN 101 DO 402 I=1,N V(2,2)=V(2,2)+X(I)**2 V(3,3)=V(3,3)+Y(I)**2 W(1)=W(1)+Z(I) 402 CONTINUE SFXY=DSQRT(V(2,2)+V(3,3)) IEXP=DLOG10(SFXY) SFXY=10.D0**IEXP V(2,2)=V(2,2)/SFXY**2 V(3,3)=V(3,3)/SFXY**2 IEXP=DLOG10(DABS(W(1))) SFZ=10.D0**IEXP W(1)=W(1)/SFZ V(1,1)=N DO 403 I=1,N XI=X(I)/SFXY YI=Y(I)/SFXY V(1,2)=V(1,2)+XI V(1,3)=V(1,3)+YI V(2,3)=V(2,3)+XI*YI ZI=Z(I)/SFZ W(2)=W(2)+XI*ZI W(3)=W(3)+YI*ZI IF (N.LE.4) GO TO 403 XX=XI**2 YY=YI**2 V(2,4)=V(2,4)+XX*XI V(2,5)=V(2,5)+YY*XI V(3,4)=V(3,4)+XX*YI V(3,5)=V(3,5)+YY*YI V(4,4)=V(4,4)+XX*XX V(4,5)=V(4,5)+XX*YY V(5,5)=V(5,5)+YY*YY W(4)=W(4)+XX*ZI W(5)=W(5)+YY*ZI IF (N.EQ.5) GO TO 403 V(4,6)=V(4,6)+XX*XI*YI V(5,6)=V(5,6)+YY*YI*XI W(6)=W(6)+XI*YI*ZI 403 CONTINUE NP=3 IPRBLD=110 IF (N.LE.4) GO TO 102 V(1,4)=V(2,2) V(1,5)=V(3,3) NP=5 IPRBLD=100 IF (N.EQ.5) GO TO 102 V(1,6)=V(2,3) V(2,6)=V(3,4) V(3,6)=V(2,5) V(6,6)=V(4,5) NP=6 IPRBLD=0 102 DO 404 I=2,NP LI=I-1 DO 404 J=1,LI V(I,J)=V(J,I) 404 CONTINUE DO 405 I=1,3 NO=I*1.6+2 IPRBLD=IPRBLD+ISOLN(X,Y,Z,V,W,SFXY,SFZ,N,NO,P(1,I),R(1,I),S(I))* @ 10**(I-1) GO TO (103,104,105),I 104 A(1,1)=P(4,2) A(2,1)=P(5,2) A(3,1)=-P(2,2)/(2.D0*A(1,1)) A(4,1)=-P(3,2)/(2.D0*A(2,1)) A(5,1)=P(1,2)-A(1,1)*A(3,1)**2-A(2,1)*A(4,1)**2 103 IF (NP.EQ.NO) RETURN GO TO 405 105 O=DATAN2(P(6,3),P(4,3)-P(5,3)) A(6,2)=O*90.D0/3.141592653589793D0 TA=P(4,3)+P(5,3) TB=P(6,3)/DSIN(O) A(1,2)=(TA+TB)/2.D0 A(2,2)=(TA-TB)/2.D0 TA=4.D0*P(4,3)*P(5,3)-P(6,3)**2 A(3,2)=(P(3,3)*P(6,3)-2.D0*P(2,3)*P(5,3))/TA A(4,2)=(P(2,3)*P(6,3)-2.D0*P(3,3)*P(4,3))/TA A(5,2)=P(1,3)-P(4,3)*A(3,2)**2-P(5,3)*A(4,2)**2-P(6,3)*A(3,2)* @ A(4,2) 405 CONTINUE RETURN END FUNCTION ISOLN(X,Y,Z,V,W,SFXY,SFZ,N,NP,P,R,S) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),Y(N),Z(N),R(N),P(NP),V(6,6),W(6),IWA(6),IWB(6), @ U(6,6) DO 400 I=1,6 DO 400 J=1,6 U(I,J)=V(I,J) 400 CONTINUE ISOLN=1 CALL MINVD(U,6,NP,DET,IWA,IWB) IF (DABS(DET).LT.1.D-10) RETURN ISOLN=0 SF=SFZ DO 401 I=1,NP DO 402 J=1,NP P(I)=P(I)+U(I,J)*W(J) 402 CONTINUE IF (MOD(I,2).EQ.0.AND.I.NE.6) SF=SF/SFXY P(I)=P(I)*SF 401 CONTINUE S=0.D0 DO 403 I=1,N R(I)=P(1)+P(2)*X(I)+P(3)*Y(I)-Z(I) IF (NP.GE.5) R(I)=R(I)+P(4)*X(I)**2+P(5)*Y(I)**2 IF (NP.EQ.6) R(I)=R(I)+P(6)*X(I)*Y(I) S=S+R(I)**2 403 CONTINUE IDF=N-NP IF (IDF.EQ.0) IDF=1 S=DSQRT(S/IDF) RETURN END SUBROUTINE PLANE (A,F,ALAT,ALON,N,X,Y,ALATO,ALONO) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALAT(N),ALON(N),X(N),Y(N) DATA PI/3.141592653589793D0/ PHIM=0.D0 ALAM=0.D0 RD=PI/180.D0 DO 401 I=1,N PHIM=PHIM+ALAT(I) ALAM=ALAM+ALON(I) 401 CONTINUE ALATO=PHIM/N PHIM=ALATO*RD ALONO=ALAM/N ALAM=ALONO*RD ESQ=(1.D0/F)**2 B=DSQRT(1.D0-ESQ*DSIN(PHIM)**2) RHO=A*(1.D0-ESQ)/B**3 ANU=A/B DO 402 I=1,N X(I)=(ALAT(I)*RD-PHIM)*RHO Y(I)=(ALON(I)*RD-ALAM)*ANU*DCOS(PHIM) 402 CONTINUE RETURN END