*DECK MAIN PROGRAM MAIN C------------------------------------------------------------------------ C NOTE: for a more complete code that also interpolate data see GRID.f -- C------------------------------------------------------------------------ C PROGRAM TO GENERATE A CURVILINEAR ORTHOGONAL GRID. HERE, THE C BORDERS OF THE DOMAIN ARE SPECIFIED IN DATA STATEMENTS. C The program is conceptually simple, but is not without its C its pitfalls. For example, a bad guess of the border points C will lead to focusing of adjacent coordinates and failure. C We also use another scheme but it is long and not yet C documented. C------------------------------------------------------------------------ PARAMETER (IM=120,JM=39,IMJM=IM*JM,II=IM/3,JJ=JM/3) PARAMETER (NT=7,NB=7,NL=2,NR=2) REAL LFTX,LFTY DIMENSION ITOP(NT), IBOT(NT), & TOPX(NT),TOPY(NT),BOTX(NT),BOTY(NT), & JRHT(NR), JLFT(NL), & RHTX(NR),RHTY(NR),LFTX(NL),LFTY(NL) DIMENSION RITOP(NT),RIBOT(NB),RJLFT(NL),RJRHT(NR), & RI(IM),RJ(IM) DIMENSION XXX(IM),YYY(IM), & FM(IM,JM),DNI(IM,JM), & A(IM,JM),P(IM,JM),Q(IM,JM),R(IM,JM),S(IM,JM), & AXP(IM,JM),AXM(IM,JM),AYP(IM,JM),AYM(IM,JM) REAL X(IM,JM),Y(IM,JM),XX(IM,JM),YY(IM,JM),DX(IM,JM),DY(IM,JM) REAL X2(II,JJ),Y2(II,JJ),DX2(II,JJ),DY2(II,JJ) C DATA ZER0/1.E-38/ DATA PI/3.14159/ C DATA ITOP/ 1, 20, 40, 60, 80, 100,120/ DATA TOPX/-8.,-1., 6., 12.5,22.,31.5, 38./ DATA TOPY/ 37.5,37.5,44.,46.1,42., 39., 37./ C DATA IBOT/ 1 , 20, 40, 60, 80,100, 120 / DATA BOTX/-8., -1.,6.,12.5,22., 31.5, 38./ DATA BOTY/ 34.,34.,34.,32.,30., 30.5,31./ C DATA JLFT/ 1,39/ DATA LFTX/-8.,-8./ DATA LFTY/34. ,37.5/ C C DATA JRHT/ 1, 39/ C DATA RHTX/38.,38./ C DATA RHTY/31.,37./ C DO 1000 J=1,JM DO 1000 I=1,IM X(I,J)=0. 1000 Y(I,J)=0. DO 1001 J=1,JJ DO 1001 I=1,II DX2(I,J)=0. 1001 DY2(I,J)=0. DO 1002 I=1,IM XXX(I)=0. 1002 YYY(I)=0. C C------ FILL IN BORDER POINTS ------------------------------ DO 10 I=1,IM 10 RI(I)=I DO 20 I=1,NT 20 RITOP(I)=ITOP(I) CALL SPLINE (RITOP,TOPX,NT,RI,XXX,IM) CALL SPLINE (RITOP,TOPY,NT,RI,YYY,IM) DO 30 I=1,IM X(I,JM)=XXX(I) Y(I,JM)=YYY(I) 30 CONTINUE C DO 40 I=1,NB 40 RIBOT(I)=IBOT(I) CALL SPLINE (RIBOT,BOTX,NB,RI,XXX,IM) CALL SPLINE (RIBOT,BOTY,NB,RI,YYY,IM) DO 50 I=1,IM X(I,1)=XXX(I) Y(I,1)=YYY(I) 50 CONTINUE C DO 60 J=1,JM 60 RJ(J)=FLOAT(J) DO 70 J=1,5 70 RJLFT(J)=JLFT(J) CALL SPLINE (RJLFT,LFTX,NL,RJ,XXX,JM) CALL SPLINE (RJLFT,LFTY,NL,RJ,YYY,JM) DO 80 I=1,JM X(1,I)=XXX(I) Y(1,I)=YYY(I) 80 CONTINUE C C DO 90 J=1,NR C 90 RJRHT(J)=JRHT(J) C CALL SPLINE (RJRHT,RHTX,NR,RJ,XXX,JM) C CALL SPLINE (RJRHT,RHTY,NR,RJ,YYY,JM) C DO 100 I=1,JM C X(IM,I)=XXX(I) C Y(IM,I)=YYY(I) C 100 CONTINUE C CALL PRXY('X BEFORE FILLING',0.,X,IM,1,JM,1,0.) CALL PRXY('Y BEFORE FILLING',0.,Y,IM,1,JM,1,0.) C C --------------- SET FIRST GUESS ------------------ DO 15 J=2,JM-1 DO 15 I=2,IM X(I,J)=X(I,1)+(X(I,JM)-X(I,1))*(Y(1,J)-Y(1,1))/(Y(1,JM)-Y(1,1)) Y(I,J)=Y(I,1)+(Y(I,JM)-Y(I,1))*(Y(1,J)-Y(1,1))/(Y(1,JM)-Y(1,1)) 15 CONTINUE C CALL PRXY ('X BEFORE ORTHOG',0.,X,IM,1,JM,1,0.) CALL PRXY ('Y BEFORE ORTHOG',0.,Y,IM,1,JM,1,0.) C C CALL POISON(Y,FM,P,Q,R,S,A,DNI,AXP,AXM,AYP,AYM,IM,JM) C CALL ORTHOG (X,Y,IM,JM) CALL ORTHOG (X,Y,IM,JM) C CALL PRXY ('X AFTER ORTHOG',0.,X,IM,1,JM,1,0.) CALL PRXY ('Y AFTER ORTHOG',0.,Y,IM,1,JM,1,0.) C------ DECIMATE FOR PLOTTING PURPOSES --------------------- DO 35 J=1,JJ DO 35 I=1,II X2(I,J)=X(3*I-1,3*J-1) 35 Y2(I,J)=Y(3*I-1,3*J-1) CALL PRXY ('X2',0.,X2,II,1,JJ,1,0.) CALL PRXY ('Y2',0.,Y2,II,1,JJ,1,0.) WRITE(78,788) II,JJ 788 FORMAT(2I5) WRITE(78,787) X2,Y2 787 FORMAT(6E12.5) C STOP END *DECK ORTHOG SUBROUTINE ORTHOG (X,Y,IM,JM) PARAMETER(IMM=120,JMM=39,IJMM=IMM*JMM) DIMENSION X(IM,JM),Y(IM,JM),DSI(IMM,JMM),DSJ(IMM,JMM),CS(IMM,JMM) DATA PI/3.141593/ DATA DSI/IJMM*0./,DSJ/IJMM*0./,CS/IJMM*0./ C DO 30 J=2,JM DO 30 I=1,IM 30 CS(I,J)=COS(PI*.5*(Y(I,J)+Y(I,J-1))/180.) DO 40 J=2,JM I=1 DSJ(I,J)=SQRT(((X(I,J)-X(I,J-1))*CS(I,J))**2+(Y(I,J)-Y(I,J-1))**2) I=IM DSJ(I,J)=SQRT(((X(I,J)-X(I,J-1))*CS(I,J))**2+(Y(I,J)-Y(I,J-1))**2) DO 40 I=2,IM-1 DSJ(I,J)=SQRT(((X(I,J)-X(I,J-1))*CS(I,J))**2+(Y(I,J)-Y(I,J-1))**2) 40 CONTINUE C DO 400 J=2,JM C DO 300 LOOP=1,100 I=1 DSI(I,J)=.5*(SQRT(((X(I,J)-X(I+1,J)+X(I,J-1)-X(I+1,J-1)) 1 *CS(I,J))**2 +(Y(I,J)-Y(I+1,J)+Y(I,J-1)-Y(I+1,J-1))**2)) I=IM DSI(I,J)=.5*(SQRT(((X(I,J)-X(I-1,J)+X(I,J-1)-X(I-1,J-1)) 1 *CS(I,J))**2 +(Y(I,J)-Y(I-1,J)+Y(I,J-1)-Y(I-1,J-1))**2)) DO 150 I=2,IM-1 DSI(I,J)=.25*(SQRT(((X(I+1,J)-X(I-1,J)+X(I+1,J-1)-X(I-1,J-1)) 1 *CS(I,J))**2 +(Y(I+1,J)-Y(I-1,J)+Y(I+1,J-1)-Y(I-1,J-1))**2)) 150 CONTINUE X(1,J)=X(1,J-1)-(DSJ(1,J)/DSI(1,J)/CS(1,J)) 1 *.5*(Y(2,J)-Y(1,J)+Y(2,J-1)-Y(1,J-1)) Y(1,J)=Y(1,J-1)+(DSJ(1,J)/DSI(1,J))*CS(1,J) 1 *.5*(X(2,J)-X(1,J)+X(2,J-1)-X(1,J-1)) X(IM,J)=X(IM,J-1)-(DSJ(IM,J)/DSI(IM,J)/CS(IM,J)) 1 *.5*(Y(IM,J)-Y(IM-1,J)+Y(IM,J-1)-Y(IM-1,J-1)) Y(IM,J)=Y(IM,J-1)+(DSJ(IM,J)/DSI(IM,J))*CS(IM,J) 1 *.5*(X(IM,J)-X(IM-1,J)+X(IM,J-1)-X(IM-1,J-1)) DO 100 I=2,IM-1 X(I,J)=X(I,J-1)-(DSJ(I,J)/DSI(I,J)/CS(I,J)) 1 *.25*(Y(I+1,J)-Y(I-1,J)+Y(I+1,J-1)-Y(I-1,J-1)) Y(I,J)=Y(I,J-1)+(DSJ(I,J)/DSI(I,J))*CS(I,J) 1 *.25*(X(I+1,J)-X(I-1,J)+X(I+1,J-1)-X(I-1,J-1)) C 100 CONTINUE DO 200 I=IM-1,2,-1 X(I,J)=X(I,J-1)-(DSJ(I,J)/DSI(I,J)/CS(I,J)) 1 *.25*(Y(I+1,J)-Y(I-1,J)+Y(I+1,J-1)-Y(I-1,J-1)) Y(I,J)=Y(I,J-1)+(DSJ(I,J)/DSI(I,J))*CS(I,J) 1 *.25*(X(I+1,J)-X(I-1,J)+X(I+1,J-1)-X(I-1,J-1)) 200 CONTINUE DO 250 I=1,IM 250 CS(I,J)=COS(PI*.5*(Y(I,J)+Y(I,J-1))/180.) 300 CONTINUE C 400 CONTINUE C CALL PRXY('DSJ ',0.,DSJ,IM,1,JM,1,0.) C CALL PRXY('DSI ',0.,DSI,IM,1,JM,1,0.) RETURN END *DECK POISON SUBROUTINE POISON (F, FM, P, Q, R, S, A, DNI, AXP, AXM, AYP, AYM, 1 IX, JY) C ****************************************************************** C * SOLUTION OF: * C * APX(I,J)*F(I+1,J)+AMX(I,J)*F(I-1,J) * C * +APY(I,J)*F(I,J+1)+AMY(I,J)*F(I,J-1) * C * -(APX)I,J)+AMX(I,J)+APY(I,J)+AMY(I,J))*F(I,J) * C * =S(I,J) * C * * C * H J HERRING * C * DYNALYSIS OF PRINCETON * C * FEBRUARY 1980 * C * * C ****************************************************************** C * * C * NOTATION * C * P,Q,R,S DYNAMIC STORAGE * C * A,DNI DYNAMIC STORAGE * C * AXP,AXM DYNAMIC STORAGE * C * AYP,AYM DYNAMIC STORAGE * C * IX,JY DIMENSIONS OF I AND J INDICES OF MATRICES * C * * C * VERSION (11/4/82) * C ****************************************************************** C VERSION (11/4/82) DIMENSION F(IX,JY), FM(IX,JY), A(IX,JY), DNI(IX,JY) DIMENSION P(IX,JY), Q(IX,JY), R(IX,JY), S(IX,JY) DIMENSION AXP(IX,JY), AXM(IX,JY), AYP(IX,JY), AYM(IX,JY) DO 1010 J=1,JM DO 1010 I=1,IM AXP(I,J)=1. AXM(I,J)=1. AYP(I,J)=1. AYM(I,J)=1. S(I,J)=0. FM(I,J)=1. 1010 CONTINUE DO 220 J=1, JY DO 210 I=1, IX A(I,J)=AXP(I,J)+AXM(I,J)+AYP(I,J)+AYM(I,J) 210 CONTINUE 220 CONTINUE DNI(1,1)=1./A(1,1) C P(1,1)=AXP(1,1)*DNI(1,1) C Q(1,1)=AYP(1,1)*DNI(1,1) P(1,1)=0. Q(1,1)=0. DO 300 I=2, IX DNI(I,1)=1./(A(I,1)-AYM(I,1)) C P(I,1)=AXP(I,1)*DNI(I,1) C Q(I,1)=AYP(I,1)*DNI(I,1) P(I,1)=0. Q(I,1)=0. 300 CONTINUE DO 310 J=2, JY DNI(1,J)=1./(A(1,J)-AXM(1,J)) C P(1,J)=AXP(1,J)*DNI(1,J) C Q(1,J)=AYP(1,J)*DNI(1,J) P(1,J)=0. Q(1,J)=0. 310 CONTINUE DO 340 J=2, JY DO 320 I=2, IX DNI(I,J)=1./(A(I,J)-AXM(I,J)*P(I-1,J)-AYM(I,J)*Q(I,J-1)) P(I,J)=AXP(I,J)*DNI(I,J) Q(I,J)=AYP(I,J)*DNI(I,J) 320 CONTINUE 340 CONTINUE IXM=IX-1 JYM=JY-1 C------------------------------------------------------------------ C ITERATE TOWARDS SOLUTION C------------------------------------------------------------------ DO 900 L=1, 100 C R(1,1)=S(1,1)*DNI(1,1) R(1,1)=F(1,1) DO 200 I=2, IX C R(I,1)=(S(I,1)+AXM(I,1)*(Q(I-1,1)*F(I-1,2)+R(I-1,1)))*DNI(I,1) R(I,1)=F(I,1) 200 CONTINUE DO 250 J=2, JYM C R(1,J)=(S(1,J)+AYM(1,J)*(P(1,J-1)*F(2,J-1)+R(1,J-1)))*DNI(1,J) R(1,J)=F(1,J) 250 CONTINUE DO 400 J=2, JYM DO 350 I=2, IXM R(I,J)=(S(I,J)+AXM(I,J)*(Q(I-1,J)*F(I-1,J+1)+R(I-1,J)) 1 +AYM(I,J)*(P(I,J-1)*F(I+1,J-1)+R(I,J-1)))*DNI(I,J) 350 CONTINUE C R(IX,J)=(S(IX,J)+AXM(IX,J)*(Q(IXM,J)*F(IXM,J+1)+R(IXM,J)) C 1 +AYM(IX,J)*(P(IX,J-1)*F(IX,J-1)+R(IX,J-1)))*DNI(IX,J) 400 CONTINUE C R(1,JY)=(S(1,JY)+AYM(1,JY)*(P(1,JYM)*F(2,JYM)+R(1,JYM)))*DNI(1,JY) C DO 500 I=2, IXM C R(I,JY)=(S(I,JY)+AXM(I,JY)*(Q(I-1,JY)*F(I-1,JY)+R(I-1,JY)) C 1 +AYM(I,JY)*(P(I,JYM)*F(I+1,JYM)+R(I,JYM)))*DNI(I,JY) C 500 CONTINUE C R(IX,JY)=(S(IX,JY)+AXM(IX,JY)*(Q(IXM,JY)*F(IXM,JY)+R(IXM,JY)) C 1 +AYM(IX,JY)*(P(IX,JYM)*F(IX,JYM)+R(IX,JYM)))*DNI(IX,JY) C F(IX,JY)=FM(IX,JY)*R(IX,JY)/(1.-P(IX,JY)-Q(IX,JY)) C DO 600 II=2, IX C I=IX-II+1 C F(I,JY)=FM(I,JY)*(P(I,JY)*F(I+1,JY)+R(I,JY))/(1.-Q(I,JY)) C 600 CONTINUE DO 800 JJ=2, JY J=JY-JJ+1 C F(IX,J)=FM(IX,J)*(Q(IX,J)*F(IX,J+1)+R(IX,J))/(1.-P(IX,J)) DO 700 II=2, IX I=IX-II+1 F(I,J)=FM(I,J)*(P(I,J)*F(I+1,J)+Q(I,J)*F(I,J+1)+R(I,J)) 700 CONTINUE 800 CONTINUE C WRITE(6,'(/,'' L ='',I4)') L C WRITE(6,'(1X,11F8.2)') ((F(I,J),I=1,IX),J=1,JY) 900 CONTINUE RETURN END *DECK PRXY SUBROUTINE PRXY (LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C DIMENSION A(IM,JM),NUM(500),LINE(500) CHARACTER LABEL*(*) DATA ZERO /1.E-10/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE SCALE=10.**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.0/SCALE WRITE(6,170) LABEL 170 FORMAT(1X,A40) WRITE(6,180) TIME,SCALE 180 FORMAT(' TIME =',F9.4,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) IB=1 DO 190 I=1,IM 190 NUM(I)=I C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(6,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,2X,24I5,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+24*ISKP GO TO 200 END *DECK ROTATE SUBROUTINE ROTATE(X,Y,N) REAL X(N),Y(N) XPI=4.*ATAN(1.) C DO 2 I=1,N B=(90.-Y(I))*XPI/180. A=X(I)*XPI/180. X(I)= ATAN(-TAN(B) * COS(A) ) Y(I)= ACOS(-SIN(B) * SIN(A) ) 2 CONTINUE C DO 3 I=1,N X(I)= 180.*X(I)/XPI 3 Y(I)= 90.-180.*Y(I)/XPI C RETURN END *DECK ROTBACK SUBROUTINE ROTBACK(X,Y,IM,JM) REAL X(IM,JM),Y(IM,JM) XPI=4.*ATAN(1.) C DO 2 J=1,JM DO 2 I=1,IM A = XPI*X(I,J)/180. B = XPI*(90.-Y(I,J))/180. Y(I,J)= ACOS( SIN(B) * COS(A) ) X(I,J)= ATAN( COS(B) /( SIN(B) * SIN(A) ) ) IF(A.GT.0..AND.B.GT.0.) X(I,J)=X(I,J)+XPI IF(A.GT.0..AND.B.LT.0.) X(I,J)=X(I,J)+XPI/2. 2 CONTINUE C DO 3 J=1,JM DO 3 I=1,IM Y(I,J)=90.-Y(I,J)*180./XPI 3 X(I,J)=X(I,J)*180./XPI C RETURN END *DECK SPLINE SUBROUTINE SPLINE(X,Y,N,XNEW,YNEW,M) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(NMAX),U(NMAX),XNEW(M),YNEW(M) YP1=1.E30 YPN=1.E30 IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE C DO 20 I =1,M CALL SPLINT(X,Y,Y2,N,XNEW(I),YNEW(I)) 20 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) C IF (H.EQ.0.) PAUSE 'B^A^D XA ^I^N^P^U^T.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END