! ! IGDRVS - modified to reorder YZ data. ! 1st dimension of YZ() et al. set to 2 ! SUBROUTINE IGDRVS(YZ,DYZDT) PARAMETER (NPTMAX=160,NZMAX=401) DIMENSION YZ(2,NPTMAX),DYZDT(2,NPTMAX) COMMON /NGAM/GAM(NPTMAX),CGAM(NPTMAX),NPTS COMMON /KFLAGB/KUCR,KFLAG COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) X ,QZ(NZMAX),Q(NZMAX) X ,UCRZ(NZMAX),UCR(NZMAX) X ,UHWZ(NZMAX),UHW(NZMAX) & ,PRESS(NZMAX),THETA(NZMAX),BVFSQ(NZMAX),RHO(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /IZTOPB/ZLAST,QB,QB1,QB4 & ,IZTOPT,IZTOPQ,IZTOPU,IZTOPHW,NUTTOP COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3 DIMENSION Y(NPTMAX),Z(NPTMAX),DYDT(NPTMAX),DZDT(NPTMAX) SMALL = 1.0e-30 C CC PRINT 995, T,(YZ(1,I),I=1,6),(YZ(2,I),I=1,6) CC & ,(GAM(I),I=1,6),R2PI,NPTS CC 995 FORMAT('T,YZ = ',25E13.4/'GAM,R2PI,NPTS = ',14E13.4,I5) ! set up crosswind conditional so it's not in the loop below IF(KUCR.EQ.1 .AND. KFLAG.EQ.1) THEN FUCR=0.0 ELSE FUCR=1.0 ENDIF DO 10 I=1,NPTS Y(I)=YZ(1,I) Z(I)=YZ(2,I) 10 CONTINUE ! for the Jth vortex ... DO 30 J=1,NPTS SUMDY =0. SUMDZ =0. ! compute crosswind and headwind at vortex altitude UCRK = FK(Z(J),UCRZ,UCR,NZUCR,IZTOPU) UHWK = FK(Z(J),UHWZ,UHW,NZUHW,IZTOPU) ! compute velocity contribution from all vortices and their images DO 20 I=1,NPTS !!! IF(I.EQ.J) GO TO 20 DELY=Y(J)-Y(I) DELZ=Z(J)-Z(I) SUMZ=Z(J)+Z(I) DSQ=DELY*DELY+DELZ*DELZ+SMALL DSI=DELY*DELY+SUMZ*SUMZ CC PRINT 993, J,I,DELY,DELZ,DSQ CC 993 FORMAT(/' J,I,DELY,DELZ,DSQ = ',2I10,3E13.4) SUMDY=SUMDY+CGAM(I)*(SUMZ/DSI-DELZ/DSQ) SUMDZ=SUMDZ+CGAM(I)*(DELY/DSQ-DELY/DSI) CC PRINT 994, SUMDY,SUMDZ CC 994 FORMAT(/' SUMDY,SUMDZ = ',2E13.4) 20 CONTINUE ! total advected velocity is given by ... DYDT(J) = ( R2PI*SUMDY + UCRK*FUCR ) * ( 1 - UHWK*ACSPEEDI ) DZDT(J) = ( R2PI*SUMDZ - UHWK*TAN3 ) * ( 1 - UHWK*ACSPEEDI ) ! DZDT(J) = ( R2PI*SUMDZ - UHWK* 0 ) * ( 1 - UHWK*ACSPEEDI ) !!! !!! IF(J.EQ.1 .OR. J.EQ.2 .OR. J.EQ.5 .OR. J.EQ.6 .OR. !!! X J.EQ.9 .OR. J.EQ.10) THEN !!! UCRK = FK(Z(J),UCRZ,UCR,NZUCR,IZTOPU) !!! UHWK = FK(Z(J),UHWZ,UHW,NZUHW,IZTOPU) !!! IF(KUCR.EQ.1 .AND. KFLAG.EQ.1) UCRK = 0. !!!CC PRINT *, "J,Z(J),UCRK = ",J,Z(J),UCRK !!! DYDT(J)=-R2PI*SUMDY + UCRK !!! ELSEIF(J.EQ.3) THEN !!! DYDT(J) = DYDT(1) !!! ELSEIF(J.EQ.4) THEN !!! DYDT(J) = DYDT(2) !!! ELSEIF(J.EQ.7) THEN !!! DYDT(J) = DYDT(5) !!! ELSEIF(J.EQ.8) THEN !!! DYDT(J) = DYDT(6) !!! ELSEIF(J.EQ.11) THEN !!! DYDT(J) = DYDT(9) !!! ELSEIF(J.EQ.12) THEN !!! DYDT(J) = DYDT(10) !!! ENDIF !!! DZDT(J)=R2PI*SUMDZ !!! !wj write(2,*)"hi joe: J,DYDT,DZDT =",J,DYDT(J),DZDT(J) 30 CONTINUE C DO 90 J=1,NPTS DYZDT(1,J)=DYDT(J) DYZDT(2,J)=DZDT(J) ! write(6,*)"IGDRVS - J,Y,Z,CGAM =",J,Y(J),Z(J),CGAM(J) 90 CONTINUE RETURN END SUBROUTINE IGDRVS ! ! RKQC - this is a general-use routine with no common blocks; there ! is no need to reveal details of the variable dimension structure. ! ! added variable GAM to disregard zero circulation vortices from ! error calculation ! SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT & ,GAM,GAMSMALL,ZMIN) PARAMETER (NPTMAX=160,NMAX=2*NPTMAX,FCOR=.0666666667, & ONE=1.,SAFETY=0.9) DIMENSION Y(N),DYDX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),DYSAV(NMAX) DIMENSION GAM(N/2),ZMIN(N/2) PGROW=-0.10 PSHRNK=-0.125 ERRCON = (4./SAFETY)**(1./PGROW) HMIN = 0.01 HMAX = 0.4 HTRY = AMAX1(HMIN,HTRY) HTRY = AMIN1(HMAX,HTRY) XSAV=X do i=1,n/2 y(i*2) = amax1( y(i*2) , zmin(i) ) yscal(i*2) = amax1( yscal(i*2) , zmin(i) ) enddo DO 11 I=1,N YSAV( I)=Y( I) DYSAV(I)=DYDX(I) 11 CONTINUE H=HTRY 1 HH=0.5*H CC print 96 CC 96 format(' rk4old try1') CALL RK4OLD(YSAV,DYSAV,N,XSAV,HH,YTEMP,ZMIN) do i=1,n/2 ytemp(i*2) = amax1( ytemp(i*2) , zmin(i) ) enddo CC print 99 CC 99 format(' rk4old ok1') X=XSAV+HH CALL IGDRVS(YTEMP,DYDX) CALL RK4OLD(YTEMP,DYDX,N,X,HH,Y,ZMIN) do i=1,n/2 y(i*2) = amax1( y(i*2) , zmin(i) ) enddo CC print 98 CC 98 format(' rk4old ok2') X=XSAV+H C CCC IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.' IF(X.EQ.XSAV)THEN WRITE(2,991) 991 FORMAT(/'Stepsize not significant in RKQC.') CLOSE(2) CCC CLOSE(3) STOP ENDIF C CALL RK4OLD(YSAV,DYSAV,N,XSAV,H,YTEMP,ZMIN) do i=1,n/2 ytemp(i*2) = amax1( ytemp(i*2) , zmin(i) ) enddo CC print 97 CC 97 format(' rk4old ok3') ERRMAX=0. DO 12 I=1,N YTEMP(I)=Y(I)-YTEMP(I) I2 = (I+1)/2 if(abs(GAM(I2)).ge.GAMSMALL) then ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I))) endif 12 CONTINUE ERRMAX=ERRMAX/EPS CC PRINT 997, SAFETY,H,ERRMAX,PSHRNK CC 997 FORMAT('SAFETY,H,ERRMAX,PSHRNK = ',4E15.5) CC PRINT *, Y !JOE IF(ERRMAX.GT.ONE) THEN c IF(ERRMAX.GT.ONE .AND. H.GT.HMIN) THEN c H=SAFETY*H*(ERRMAX**PSHRNK) c GOTO 1 c ELSE HDID=H c IF(ERRMAX.GT.ERRCON)THEN c HNEXT=SAFETY*H*(ERRMAX**PGROW) c ELSE c HNEXT=amax1(4.*H,HMAX) c ENDIF hnext = h c ENDIF DO 13 I=1,N Y(I)=Y(I)+YTEMP(I)*FCOR 13 CONTINUE do i=1,n/2 y(i*2) = amax1( y(i*2) , zmin(i) ) enddo RETURN END SUBROUTINE RKQC C ! RK4OLD - this is a general-use routine with no common blocks; there ! is no need to reveal details of the variable dimension structure. ! SUBROUTINE RK4OLD(Y,DYDX,N,X,H,YOUT,ZMIN) PARAMETER (NPTMAX=160,NMAX=2*NPTMAX) DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX) & ,ZMIN(N/2) HH=H*0.5 H6=H/6. XH=X+HH DO 11 I=1,N YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE do i=1,n/2 yt(2*i) = amax1(yt(2*i),zmin(i)) enddo CC print 91 CC 91 format(' igdrvs try1') CALL IGDRVS(YT,DYT) CC print 92 CC 92 format(' igdrvs ok1') DO 12 I=1,N YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE do i=1,n/2 yt(2*i) = amax1(yt(2*i),zmin(i)) enddo CC print 93 CC 93 format(' igdrvs try2') CALL IGDRVS(YT,DYM) CC print 94 CC 94 format(' igdrvs ok2') DO 13 I=1,N YT(I)=Y(I)+H*DYM(I) 13 CONTINUE do i=1,n/2 yt(2*i) = amax1(yt(2*i),zmin(i)) enddo do I=1,N DYM(I)=DYT(I)+DYM(I) enddo CC print 95 CC 95 format(' igdrvs try3') CALL IGDRVS(YT,DYT) CC print 96 CC 96 format(' igdrvs ok3') DO 14 I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE do i=1,n/2 yout(2*i) = amax1(yout(2*i),zmin(i)) enddo RETURN END SUBROUTINE RK4OLD C ! GTOGL - this is a general-use routine with no common blocks; there ! is no need to reveal details of the variable dimension structure. ! 1D interpolation, that's all that's going on here.. ! SUBROUTINE GTOGL(NX,X,Y,NXI,XI,YI) DIMENSION X(*),Y(*),XI(*),YI(*) INTEGER IPT,IXI,IX1,IX,IXP1 REAL FPT,PT C IPT=1 DO 20 IXI=1,NXI FPT=0. PT=XI(IXI) IF(PT.LT.X(1) .OR. PT.GT.X(NX)) GO TO 19 C IX1=IPT DO 10 IX=IX1,NX IXP1=IX+1 IF(PT.LT.X(IX) .OR. PT.GT.X(IXP1)) GO TO 10 FPT=Y(IX) + (Y(IXP1)-Y(IX))*(PT-X(IX))/(X(IXP1)-X(IX)) IPT=IX GO TO 19 10 CONTINUE C 19 CONTINUE YI(IXI)=FPT 20 CONTINUE RETURN END SUBROUTINE GTOGL C REAL FUNCTION ANGLE(DY,DZ) C COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC C HYP = SQRT(DY*DY+DZ*DZ) C IF (DY .GE. 0. .AND. DZ .GE. 0.) THEN ANGLE = ASIN(DZ/HYP) ELSEIF (DY .LT. 0. .AND. DZ .GE. 0.) THEN ANGLE = PI - ASIN(DZ/HYP) ELSEIF (DY .LE. 0. .AND. DZ .LT. 0.) THEN ANGLE = PI + ASIN(-DZ/HYP) ELSEIF (DY .GT. 0. .AND. DZ .LT. 0.) THEN ANGLE = TWOPI - ASIN(-DZ/HYP) ENDIF C RETURN END FUNCTION ANGLE C REAL FUNCTION THFACF(THIN) C COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC C TH = THIN IF (TH.LT.THBOT) TH = TH + TWOPI IF (TH.GT.THTOP) TH = TH - TWOPI C IF(THMAX1.LE.TH .AND. TH.LE.THMAX2) THEN THFACF = FACMAX ELSEIF(THMAX2.LT.TH .AND. TH.LT.THMIN1) THEN THFACF = FACMIN + DELFAC*(THMIN1-TH)/DELTH1 ELSEIF(THMIN1.LE.TH .AND. TH.LE.THMIN2) THEN THFACF = FACMIN ELSE THFACF = FACMIN + DELFAC*(TH-THMIN2)/DELTH2 ENDIF C RETURN END FUNCTION THFACF