C kap command line: /usr/csrd/kap SRS.f -ct=30 -f=SRA.f C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 c loop bicong_do7 cv stripmineing turnt off for a Kap bug c C************************************************************************* C************************************************************************* C************************************************************************* C ***************** C ARC2D : EULER VERSION FOR BENCHMARK USE. ***************** C 11/05/86 THOMAS H. PULLIAM ***************** C NASA AMES RESEARCH CENTER ***************** C ***************** C************************************************************************* C ARC2D SOLVES THE EULER EQUATIONS IN GENERALIZED CURVILINEAR COORDINATES C USING AN IMPLICIT FINITE DIFFERENCE ALGORITHM WITH APPROXIMATED C FACTORIZATION AND A DIAGONALIZATION OF THE IMPLICIT OPERATORS. C C C EULER EQUATIONS : Q_T + E_XI + F_ETA = 0 C C -1 | RHO | -1| RHO U | -1| RHO V | C Q = J |RHO U|, E = J |RHO UU + XI_X P|, F = J |RHO U V + ETA_X P| C |RHO V| |RHO VU + XI_Y P| |RHO V V + ETA_Y P| C | E | |U(E+P) - XI_T P| |V(E+P) - ETA_T P| C C U = XI_T + XI_X U + XI_Y V, V = ETA_T + ETA_X U + ETA_Y V C C -1 C J = (X_XI Y_ETA - X_ETA Y_XI) C C XI_X = J Y_ETA, XI_Y = -J X_ETA, XI_T = -X_T XI_X - Y_T XI_Y C ETA_X = -J Y_XI, ETA_Y = J X_XI, ETA_T = -X_T ETA_X - Y_T ETA_Y C C P = (GAMMA-1) (E - 0.5 RHO (U^2 + V^2)) PRESSURE C U IS X COMPONENT OF VELOCITY, V IS Y COMPONENT OF VELOCITY C E IS TOTAL ENERGY, RHO IS DENSITY C GAMMA IS THE RATIO OF SPECIFIC HEATS (USUALLY 1.4) C C NONDIMENSIONALIZATION : RHO_INFINITY FOR RHO C A_INFINITY (SPEED OF SOUND) FOR U AND V C RHO_INFINITY A_INFINITY^2 FOR E C REFERENCE LENGTH IS L C (USUALLY CHORD FOR AN AIRFOIL) C NOTE THIS GIVES P AT INFINITY = 1./GAMMA C C C THE NUMERICAL ALGORITHM IS THE CLASS OF TWO TIME LEVEL 1ST OR 2ND C ORDER ACCURATE TIME IMPLICIT SCHEMES WITH SECOND ORDER CENTRAL C DIFFERENCES IN COMPUTATIONAL SPACE. C THE IMPLICIT SIDE IS APPROXIMATELY FACTORED AND DIAGONALIZED. C SEE: PULLIAM T. H. AND CHAUSSEE D. S. , C "A DIAGONAL FORM OF AN IMPLICIT APPROXIMATE FACTORIZATION C ALGORITHM", JCP 39, 1981, P347. C C THE BENCHMARK TEST CASE SET UP HERE IS FOR THE UPPER HALF OF A C BICONVEX AIRFOIL SITTING IN A STRETCH RECTANGULAR DOMAIN. C THE GRID IS GENERATED INTERNALLY, SEE BICONG. C THE BOUNDARY CONDITIONS ARE THIN AIRFOIL TANGENCY CONDITIONS ON THE C AIRFOIL SURFACE, QUASI-ONE-DIMENSIONAL RIEMANN INVARIANT CONDITIONS C AT INFLOW, OUTFLOW AND THE TOP OF THE DOMAIN, SEE BC. C C INPUT IS HANDLED IN INPUT AND INITIA. C OUPUT IS HANDLED IN IOALL. C CONDITIONS FOR THE SUPPLIED TEST CASES ARE VARIOUS GRID SIZES AT C A MACH NUMBER = 0.8. THIS GIVES TRANSONIC FLOW AND A SHOCK FOR C A THICKNESS OF THE BICONVEX = 0.1 C C THE CODE IS SELF DOCUMENTED, SEE THE DEFINITIONS GIVEN BELOW FOR MORE C DETAIL. CONTACT THOMAS H. PULLIAM, MS 202A-1 NASA AMES RESEARCH CENTER, C MOFFETT FIELD, CA., (415-694-6417) (EMAIL: PULLIAM@PRANDTL.NAS.NASA.GOV) C TO REPORT PROBLEMS OR FOR QUESTIONS. C C************************************************************************* C************************************************************************* C C COMMON/BASE/ C C------------------------------------------------------------------------ C INDICES J,K C------------------------------------------------------------------------ C JMAX: TOTAL NUMBER OF POINTS IN XI DIRECTION (DEFAULT=193) C KMAX: TOTAL NUMBER OF POINTS IN ETA DIRECTION (DEFAULT=33 ) C JM: JMAX - 1 C KM: KMAX - 1 C JBEGIN, JEND : ABSOLUTE LOWER AND UPPER LIMITS ON J C TYPICALLY 1 AND JMAX C KBEGIN, KEND : ABSOLUTE LOWER AND UPPER LIMITS ON K C TYPICALLY 1 AND KMAX C------------------------------------------------------------------------ C AUTOMATIC PERIODIC AND NONPERIODIC LOGIC IS AVAILABLE. C C PERIODIC SWITCH FOR NONPERIODIC VRS PERIODIC CODING C PERIODIC = F NONPERIODIC OPTION C PERIODIC = T PERIODIC OPTION (DEFAULT=FALSE) C JPLUS INDEX ARRAY FOR J DIFFERENCING C JPLUS(J) = J+1 C JPLUS(JMAX) = JMAX FOR PERIODIC = F C JPLUS(JMAX) = 1 FOR PERIODIC = T C JMINU INDEX ARRAY FOR J DIFFERENCING C JMINU (J) = J+1 C JMINU (1) = 1 FOR PERIODIC = F C JMINU (1) = JMAX FOR PERIODIC = T C JLOW LOWER DO LOOP LIMIT IN J DIRECTION C JLOW = 2 FOR PERIODIC = F C JLOW = 1 FOR PERIODIC = T C JUP UPPER DO LOOP LIMIT IN J DIRECTION C JUP = JMAX-1 FOR PERIODIC = F C JUP = JMAX FOR PERIODIC = T C KLOW LOWER DO LOOP LIMIT IN K DIRECTION C KLOW = 2 TYPICALLY C KUP UPPER DO LOOP LIMIT IN K DIRECTION C KUP = KMAX-1 TYPICALLY C------------------------------------------------------------------------ C C NP: CALLS IOALL(15) TO SHOW RHO, U, V, E, P, CP, MACH, X, Y C (DEFAULT=10000) C C DT: TIME STEP (DEFAULT=5.0) C EITHER A TRUE DT FOR TIME ACCURATE COMPUTATION C OR A SCALING FACTOR FOR THE SPATIALLY VARIABLE TIME STEP C [ IN STEADY STATE CALCULATION DT IS ORDER(1.0) ] C C CPUSEC: ELAPSED CPU TIME OF THE RUN. C NOTE CHECK FUNCTION SECOND FOR INSTRUCTIONS FOR C MACHINE DEPENDENT CHANGES. C------------------------------------------------------------------------ C FSMACH: FREE STREAM MACH NUMBER (DEFAULT=0.397) C READ IN IN IOALL(1) C ALPHA: ANGLE OF ATTACK OF AIRFOIL (DEFAULT=0.0 ) C READ IN IN IOALL(1) C GAMMA: GAS CONSTANT ( 1.4 ) C GAMI: GAMMA - 1 C PI : 4.*ATAN(1.0) C------------------------------------------------------------------------ C DISSIPATION PARAMETERS C------------------------------------------------------------------------ C C DIS2X: SECOND ORDER X COEFFICIENT FOR FILTER (DEFAULT=0.0) C DIS2Y: SECOND ORDER Y COEFFICIENT FOR FILTER (DEFAULT=0.0) C READ IN IN IOALL(1) ONLY USED FOR METH = 2,3 C USED IN FILERX, FILERY, STEPFX, STEPFY C DIS4X: FOURTH ORDER X COEFFICIENT FOR FILTER (DEFAULT=0.64) C DIS4Y: FOURTH ORDER Y COEFFICIENT FOR FILTER (DEFAULT=0.64) C READ IN IN IOALL(1) ONLY USED FOR METH = 2,3 C USED IN FILERX, FILERY, STEPFX, STEPFY C------------------------------------------------------------------------ C JACDT: SWITCH TO TURN ON VARIABLE TIME STEP. (DEFAULT = 1) C FOR UNSTEADY PROBLEMS USE CONSTANT TIME STEP. C FOR STEADY - STATE USE VARIABLE TIME STEP. C READ IN IN IOALL(1) C =0 : CONSTANT TIME STEP C =1 : TIME STEP SCALED BY JACOBIAN C =2 : TIME STEP BASED ON MINIMUM 1/EIGENVALUES C =3 : TIME STEP BASED CONSTANT CFL AT EACH POINT C C------------------------------------------------------------------------ C------------------------------------------------------------------------ C RESID: ASSIGNED TO IN RESIL2 OUTPUT IN IOALL(8) ONTO TAPE7 C OUTPUT IN IOALL(5) TO TAPE6 C ZEROED IN INITIA C 2-NORM OF RESIDUAL C C IPRINT: OUTPUT SWITCH FOR FREQUENCY OF RESIDUAL HISTORY OUTPUT C ON UNIT 6 ( MAIN ) (DEFAULT=1) C C------------------------------------------------------------------------ C GRID PARAMETERS C------------------------------------------------------------------------ C DYM: MINIMUM SPACING AT BODY FOR BICONVEX AIRFOIL C YMAX: UPPER DIMENSION IN CORD LENGTHS FOR BICONVEX GRID C XMIN: LEFT DIMENSION IN CORD LENGTHS FOR BICONVEX GRID C XMAX: RIGHT DIMENSION IN CORD LENGTHS FOR BICONVEX GRID C THICK: THICKNESS RATIO FOR BICONVEX AIRFOIL C JTAIL1: FIRST SOLID BODY POINT C JTAIL2: LAST SOLID BODY POINT C SLOPE: SLOPE OF THE AIRFOIL, USED FOR BICONVEX CASE C SET IN INITIA, USED IN BC C C------------------------------------------------------------------------ C TIME ACCURACY PARAMETERS C------------------------------------------------------------------------ C A CLASS OF THREE STEP FIRST OR SECOND ORDER TIME ACCURACTE SCHEMES C CAN BE EMPLOYED. FOR STEADY STATE EULER IMPLICIT PERFORMS THE BEST. C FOR UNSTEADY, 3PT IMPLICIT IS A GOOD CHOICE. C C PHIDT TIME INTEGRATION PARAMETERS : C THETAD PHIDT = 0 THETAD = 1 EULER IMPLICIT (DEFAULT) C PHIDT = 0 THETAD = 1/2 TRAP. IMPLICIT C PHIDT = 1/2 THETAD = 1 3PT IMPLICIT C C------------------------------------------------------------------------ C NUMITE : ITERATION COUNTER FOR MAIN STEPS C ISTART: THE STARTING ITERATION, C READ IN FROM THE RESTART C NSTEPS: TOTAL NUMBER OF ITERATION STEPS C---------------------------------------------------------------------- C----------------------- VARIABLES -------------------------------- C---------------------------------------------------------------------- C Q: CONSERVATIVE VARIABLES C Q(I,J, 1) DENSITY / DETERMINANT OF METRIC JACOBIAN C Q(I,J, 2) DENSITY * U / DETERMINANT OF METRIC JACOBIAN C Q(I,J, 3) DENSITY * V / DETERMINANT OF METRIC JACOBIAN C Q(I,J, 4) ENERGY / DETERMINANT OF METRIC JACOBIAN C C PRESS : PRESSURE = GAMI*(Q(4) - 0.5*(Q(2)**2+Q(3)**2)/Q(1)) C NOTE: HAS JACOBIAN IN IT C SNDSP : A = SQRT(GAMMA*PRESS/Q(1)) C NOTE: DOESN'T HAVE JACOBIAN IN IT C S: C USED FOR EXPLICIT RIGHT HAND SIDE IN RHS, C ALTERED IN VISRHS, C SMOOTHED IN SMOOTH,FILTER C USED FOR INTERMEDIATE QUANTITIES IN THE AF SCHEME IN STEP, ETC. C AND ULTIMATELY HOLDS THE CHANGE IN Q. C XY: METRIC TRANSFORMATIONS C XY(I,J, 1) D XI / DX C XY(I,J, 2) D XI / DY C XY(I,J, 3) D ETA / DX C XY(I,J, 4) D ETA / DY C XYJ: JACOBIANS OF METRIC TRANSFORMATIONS C = ( DXI/DX * DETA/DY ) - ( DXI/DY * DETA/DX ) C XIT: D XI / DT C ETT: D ETA / DT C DS: SPATIALLY VARIABLE TIME STEP, SCALED ON METRIC JACOBIAN, C SET IN INITIA, USED IN STEP, STEPF, C EIGEN, FLMATX2, FLMATY2 C X: CARTESIAN X COORDINATES OF THE GRID. C Y: CARTESIAN Y COORDINATES OF THE GRID. C C C TO CHANGE DIMENSIONS::: MODIFY ALL THE PARAMETRER STATEMENTS C MAXJ = MAXIMUM J DIMENSION C MAXK = MAXIMUM K DIMENSION C MAXJK = PRODUCT MAXJ*MAXK C PROGRAM ARC2D IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(MAXJK*4) DIMENSION PRESS(MAXJK),SNDSP(MAXJK) DIMENSION S(MAXJK*4),XY(MAXJK*4),XYJ(MAXJK) DIMENSION XIT(MAXJK),ETT(MAXJK),DS(MAXJK) DIMENSION X(MAXJK),Y(MAXJK) DIMENSION DELTQN(MAXJK*4),SLOPE(MAXJ) DIMENSION UU(MAXJK),VV(MAXJK),CCX(MAXJK),CCY(MAXJK) DIMENSION COEF4X(MAXJK),COEF4Y(MAXJK) DIMENSION COEF2X(MAXJK),COEF2Y(MAXJK) C C COMMON/WORKSP/WORK2(MAXJ,92) COMMON/WORK/WORK1(MAXJK*10) C REAL T1,T2,T3,T4, T5, T6 COMMON/IO/IREAD,IWRIT,IVER C C EXTERNAL SECOND C SECOND REPLACED BY CPUTIME AT CSRD JULY, 1989 BY L. POINTER C C D CALL INIT_TRACE(1) IREAD = 5 IWRIT = 6 IVER = 7 C OPEN(UNIT=IREAD,FILE='SRI5',FORM='FORMATTED',STATUS='OLD') OPEN(UNIT=IWRIT,FILE='SRO6',FORM='FORMATTED',STATUS='NEW') OPEN(UNIT=IVER, FILE='SRV', FORM='FORMATTED',STATUS='NEW') CALL CPUTIM(T1) CALL ELAPSE(T2) C INPUT DATA C CALL INPUT(MAXJ,MAXJK) C C INITIALIZE CPU TIME COUNTER C C CPUSC0 = CPUTIM(DUMMY) CALL CPUTIM(T5) C WRITE(IWRIT,*) '--------------------------------------------' C WRITE(IWRIT,*)' CPUSEC AT START = ',CPUSC0 WRITE(IWRIT,*)' CPUSEC AT START = ',T5 WRITE(IWRIT,*) '--------------------------------------------' C C C INITIALIZE VARIABLES C C SET JDIMEN,KDIMEN C JDIMEN = JMAX KDIMEN = KMAX C C CALL INITIA(JDIMEN,KDIMEN,Q,PRESS,SNDSP,S,XY,XYJ, * XIT,ETT,DS,X,Y,SLOPE) C C C PERFORM INTEGRATION C ISTART = 1 IEND = NSTEPS d call trace_enter (2001) DO 10 NUMITE = ISTART,IEND CALL INTEGR(JDIMEN,KDIMEN,Q,S,PRESS,SNDSP,SLOPE, > X,Y,XY,XYJ,XIT,ETT,DS,DELTQN, > UU,VV,CCX,CCY,COEF2X,COEF2Y,COEF4X,COEF4Y) C C COMPUTE CPUSEC C C CPUSEC = CPUTIM(DUMMY) - CPUSC0 CALL CPUTIM(T6) CPUSEC = T6 -T5 C C********************************************************************* C*** RESIDUAL, SS PTS, *** C********************************************************************* IP = 5 IF( MOD( NUMITE - ISTART + 1, IPRINT) .EQ. 0 ) * CALL IOALL(IP,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y) C********************************************************************* C*** WRITE OUT RHO, U, V, E, P, CP, MACH, X, Y *** C*** ONTO TAPE6 EVERY NP ITERATIONS *** C********************************************************************* IP = 15 IF ( MOD( NUMITE - ISTART + 1, NP) .EQ. 0 ) * CALL IOALL(IP,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y) C C********************************************************************* 10 CONTINUE C********************************************************************* C C C FINAL OUTPUT C C Q SOLUTION FILE FOR PLOTTING d call trace_exit (2002) CALL IOALL(10,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y) C XY GRID FILE FOR PLOTTING CALL IOALL(16,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y) C CP VRS X CALL IOALL(6,JDIMEN,KDIMEN,Q,PRESS,SNDSP,XY,XYJ,X,Y) C*********************************************************************** CALL CPUTIM(T3) CALL ELAPSE(T4) WRITE(IVER,*) 'ELAPSED CPU TIME IN SECONDS: ',T3-T1 WRITE(IVER,*) 'MLOPS RATE: ', 1842.113398/(T3-T1) WRITE(IVER,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ', * MOD(T4-T2+1000000.0,1000000.0) CLOSE(IREAD) CLOSE(IWRIT) CLOSE(IVER) d CALL EXIT_TRACE() STOP END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C****************** BICONVEX AIRFOIL BC ******************** C*********************************************************************** SUBROUTINE BC(JDIM,KDIM,Q,PRESS,SNDSP,XY,XIT,ETT,XYJ,X,Y,SLOPE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4) DIMENSION PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM),SLOPE(JDIM) COMMON/INFNTY/RHOINF,UINF,VINF,EINF,PINF C COMMON/WORKSP/F(MAXJ),WORK0(MAXJ,91) INTEGER J1, J2, K1, K2, K3, K4, J3, J4, J5, J6, J7, J10, K5, K6, X K7, K8, K11, K14, K17, J13 C C BICONVEX AIFOIL BC C STRTIT = 12. T = ( NUMITE -1.)/ STRTIT T = MIN (T, 1.D0) SCAL = (10. -15.*T + 6.*T*T) *T**3 C C U AND V ON LOWER SURFACE C U IS EXTRAPOLATED C V IS SET TO SATISFY TANGENCY C A RECTILINEAR MESH IS ASSUMED C J4 = (JM + 30) / 32 J7 = MAX0 (MIN0 (31, JM - 2) + 1, 1) d call trace_enter (2003) CDOALL 25 J1=0,32*J4-32,32 DOUBLE PRECISION U21(J7) INTEGER J8, J9 J9 = J1 + 2 J8 = MIN0 (J1 + 2 + 32 - 1, JM) - J1 - 1 C C FIRST ORDER EXTRAPOLATION OF U ( DIVIDE BY RHO/J ) C U2 = 2. * Q(J,2,2)/Q(J,2,1) - Q(J,3,2)/Q(J,3,1) C C ZERO'TH ORDER EXTRAPOLATION OF U ( DIVIDE BY RHO/J ) U21(1:$J8) = Q(J9:$J8,2,2) / Q(J9:$J8,2,1) C C ZERO'TH ORDER EXTRAPOLATION OF U Q(J9:$J8,1,2) = U21(1:$J8) C C TANGENCY Q(J9:$J8,1,3) = U21(1:$J8) * SLOPE(J9:$J8) * SCAL C C 25 CONTINUE C C SATISFY MOMENTUM RELATION FOR NORMAL PRESSURE DERIVATION C d call trace_exit (2004) K = 1 J10 = (JM + 30) / 32 J7 = MAX0 (MIN0 (31, JM - 2) + 1, 1) d call trace_enter (2005) CDOALL 30 J2=0,32*J10-32,32 DOUBLE PRECISION P21(J7), P31(J7) INTEGER J11, J12 J12 = J2 + 2 J11 = MIN0 (J2 + 2 + 32 - 1, JM) - J2 - 1 C C PRESSURE IS EXTRAPOLATED TO FIRST ORDER INTO F C P21(1:$J11) = GAMI * (Q(J12:$J11,2,4) - 0.5 * (Q(J12:$J11,2,2) X ** 2 + Q(J12:$J11,2,3) ** 2) / Q(J12:$J11,2,1)) P31(1:$J11) = GAMI * (Q(J12:$J11,3,4) - 0.5 * (Q(J12:$J11,3,2) X ** 2 + Q(J12:$J11,3,3) ** 2) / Q(J12:$J11,3,1)) F(J12:$J11) = (2. * P21(1:$J11) * XYJ(J12:$J11,2) - P31(1:$J11) X * XYJ(J12:$J11,3)) / XYJ(J12:$J11,1) C C RESCALE Q2 AND Q3 TO NEW DENSITY SUCH THAT WHEN REFORMING PRESSURE C FROM Q1 TO Q4, THE NORMAL DERIVATIVE OF PRESSURE IS SATIFIED C C ZERO'TH ORDER EXTRAPOLATION OF RHO Q(J12:$J11,K,1) = Q(J12:$J11,K+1,1) * XYJ(J12:$J11,K+1) / XYJ( X J12:$J11,K) C C U EXTRAP ABOVE, THIS PUTS IN RHO ( WITH METRIC JACOBIAN ) Q(J12:$J11,K,2) = Q(J12:$J11,K,2) * Q(J12:$J11,K,1) C C V SET ABOVE, THIS PUTS IN RHO ( WITH METRIC JACOBIAN ) Q(J12:$J11,K,3) = Q(J12:$J11,K,3) * Q(J12:$J11,K,1) C C ENERGY C RECALCULATE THE ENERGY USING THE NEW RHO, U, V C AND EXTRAPOLATED PRESSURE. Q(J12:$J11,K,4) = F(J12:$J11) / GAMI + .5 * (Q(J12:$J11,K,2) X ** 2 + Q(J12:$J11,K,3) ** 2) / Q(J12:$J11,K,1) C 30 CONTINUE C..................................................................... C INFLOW/OUTFLOW BOUNDARIES C CHECK FOR SUPERSONIC FREE STREAM C..................................................................... C d call trace_exit (2006) IF (FSMACH .LE. 1.) THEN C C C INFLOW BOUNDARY C GI = 1./GAMMA GM1I = 1./GAMI C AINF = SQRT(GAMMA*PINF/RHOINF) HSTFS = 1./GAMI + 0.5*FSMACH**2 C J = 1 K5 = (KM + 30) / 32 K8 = MAX0 (MIN0 (31, KM - 2) + 1, 1) d call trace_enter (2007) CDOALL 50 K1=0,32*K5-32,32 DOUBLE PRECISION SNORM1(K8), XY1H1(K8), XY2H1(K8), XYJM11(K8), X RHOEXT1(K8), RJM11(K8), UEXT1(K8), VEXT1(K8), EEXT1(K8), PEXT1( X K8), RFIX1(K8), REXT1(K8), QN1(K8), CSPE1(K8), C21(K8), QT1(K8) X , ENTRO1(K8), U1(K8), V1(K8), PRES1(K8), RJ1(K8) INTEGER K9, K10 K10 = K1 + 2 K9 = MIN0 (K1 + 2 + 32 - 1, KM) - K1 - 1 SNORM1(1:$K9) = 1. / SQRT (XY(J,K10:$K9,1) ** 2 + XY(J,K10:$K9,2) X ** 2) XY1H1(1:$K9) = XY(J,K10:$K9,1) * SNORM1(1:$K9) XY2H1(1:$K9) = XY(J,K10:$K9,2) * SNORM1(1:$K9) C C EXTRAPOLATE TEMPORARY VARIABLES C XYJM11(1:$K9) = XYJ(J+1,K10:$K9) RHOEXT1(1:$K9) = Q(J+1,K10:$K9,1) * XYJM11(1:$K9) RJM11(1:$K9) = 1. / Q(J+1,K10:$K9,1) UEXT1(1:$K9) = Q(J+1,K10:$K9,2) * RJM11(1:$K9) VEXT1(1:$K9) = Q(J+1,K10:$K9,3) * RJM11(1:$K9) EEXT1(1:$K9) = Q(J+1,K10:$K9,4) * XYJM11(1:$K9) PEXT1(1:$K9) = GAMI * (EEXT1(1:$K9) - 0.5 * RHOEXT1(1:$K9) * ( X UEXT1(1:$K9) ** 2 + VEXT1(1:$K9) ** 2)) C C FIX INCOMING RIEMANN INVARIANT RFIX1(1:$K9) = XY1H1(1:$K9) * UINF + XY2H1(1:$K9) * VINF + 2. * X AINF * GM1I REXT1(1:$K9) = XY1H1(1:$K9) * UEXT1(1:$K9) + XY2H1(1:$K9) * VEXT1( X 1:$K9) - 2. * SQRT (GAMMA * PEXT1(1:$K9) / RHOEXT1(1:$K9)) * X GM1I C C COMPUTE FLOW VARIABLES BASED ON ABOVE CALC C QN1(1:$K9) = (RFIX1(1:$K9) + REXT1(1:$K9)) * 0.5 CSPE1(1:$K9) = (RFIX1(1:$K9) - REXT1(1:$K9)) * GAMI * 0.25 C21(1:$K9) = CSPE1(1:$K9) ** 2 C QT1(1:$K9) = (-XY2H1(1:$K9) * UINF + XY1H1(1:$K9) * VINF) ENTRO1(1:$K9) = GAMMA C C COMPUTE FLOW VARIABLES C U1(1:$K9) = XY1H1(1:$K9) * QN1(1:$K9) - XY2H1(1:$K9) * QT1(1:$K9) V1(1:$K9) = XY2H1(1:$K9) * QN1(1:$K9) + XY1H1(1:$K9) * QT1(1:$K9) Q(J,K10:$K9,1) = (C21(1:$K9) * ENTRO1(1:$K9) * GI) ** GM1I PRES1(1:$K9) = C21(1:$K9) * Q(J,K10:$K9,1) * GI C C ADD JACOBIAN C RJ1(1:$K9) = 1. / XYJ(J,K10:$K9) Q(J,K10:$K9,1) = Q(J,K10:$K9,1) * RJ1(1:$K9) Q(J,K10:$K9,2) = Q(J,K10:$K9,1) * U1(1:$K9) Q(J,K10:$K9,3) = Q(J,K10:$K9,1) * V1(1:$K9) Q(J,K10:$K9,4) = PRES1(1:$K9) * GM1I * RJ1(1:$K9) + 0.5 * Q(J, X K10:$K9,1) * (U1(1:$K9) ** 2 + V1(1:$K9) ** 2) 50 CONTINUE C C C OUTFLOW BOUNDARY C C d call trace_exit (2008) J = JMAX K11 = (KM + 30) / 32 K8 = MAX0 (MIN0 (31, KM - 2) + 1, 1) d call trace_enter (2009) CDOALL 55 K2=0,32*K11-32,32 DOUBLE PRECISION SNORM2(K8), XY1H2(K8), XY2H2(K8), XYJM12(K8), X RHOEXT2(K8), RJM12(K8), UEXT2(K8), VEXT2(K8), EEXT2(K8), PEXT2( X K8), RFIX2(K8), REXT2(K8), QN2(K8), CSPE2(K8), C22(K8), QT2(K8) X , ENTRO2(K8), U3(K8), V2(K8), PRES2(K8), RJ2(K8) INTEGER K12, K13 K13 = K2 + 2 K12 = MIN0 (K2 + 2 + 32 - 1, KM) - K2 - 1 SNORM2(1:$K12) = 1. / SQRT (XY(J,K13:$K12,1) ** 2 + XY(J,K13:$K12, X 2) ** 2) XY1H2(1:$K12) = XY(J,K13:$K12,1) * SNORM2(1:$K12) XY2H2(1:$K12) = XY(J,K13:$K12,2) * SNORM2(1:$K12) C C EXTRAPOLATE TEMPORARY VARIABLES C XYJM12(1:$K12) = XYJ(J-1,K13:$K12) RHOEXT2(1:$K12) = Q(J-1,K13:$K12,1) * XYJM12(1:$K12) RJM12(1:$K12) = 1. / Q(J-1,K13:$K12,1) UEXT2(1:$K12) = Q(J-1,K13:$K12,2) * RJM12(1:$K12) VEXT2(1:$K12) = Q(J-1,K13:$K12,3) * RJM12(1:$K12) EEXT2(1:$K12) = Q(J-1,K13:$K12,4) * XYJM12(1:$K12) PEXT2(1:$K12) = GAMI * (EEXT2(1:$K12) - 0.5 * RHOEXT2(1:$K12) * ( X UEXT2(1:$K12) ** 2 + VEXT2(1:$K12) ** 2)) C C FIX INCOMING RIEMANN INVARIANT RFIX2(1:$K12) = XY1H2(1:$K12) * UINF + XY2H2(1:$K12) * VINF - 2. * X AINF * GM1I REXT2(1:$K12) = XY1H2(1:$K12) * UEXT2(1:$K12) + XY2H2(1:$K12) * X VEXT2(1:$K12) + 2. * SQRT (GAMMA * PEXT2(1:$K12) / RHOEXT2(1:$ X K12)) * GM1I C C COMPUTE FLOW VARIABLES BASED ON ABOVE CALC C QN2(1:$K12) = (RFIX2(1:$K12) + REXT2(1:$K12)) * 0.5 CSPE2(1:$K12) = (REXT2(1:$K12) - RFIX2(1:$K12)) * GAMI * 0.25 C22(1:$K12) = CSPE2(1:$K12) ** 2 C QT2(1:$K12) = (-XY2H2(1:$K12) * UEXT2(1:$K12) + XY1H2(1:$K12) X * VEXT2(1:$K12)) ENTRO2(1:$K12) = RHOEXT2(1:$K12) ** GAMMA / PEXT2(1:$K12) C C COMPUTE FLOW VARIABLES C U3(1:$K12) = XY1H2(1:$K12) * QN2(1:$K12) - XY2H2(1:$K12) * QT2(1:$ X K12) V2(1:$K12) = XY2H2(1:$K12) * QN2(1:$K12) + XY1H2(1:$K12) * QT2(1:$ X K12) Q(J,K13:$K12,1) = (C22(1:$K12) * ENTRO2(1:$K12) * GI) ** GM1I PRES2(1:$K12) = C22(1:$K12) * Q(J,K13:$K12,1) * GI C C ADD JACOBIAN C RJ2(1:$K12) = 1. / XYJ(J,K13:$K12) Q(J,K13:$K12,1) = Q(J,K13:$K12,1) * RJ2(1:$K12) Q(J,K13:$K12,2) = Q(J,K13:$K12,1) * U3(1:$K12) Q(J,K13:$K12,3) = Q(J,K13:$K12,1) * V2(1:$K12) Q(J,K13:$K12,4) = PRES2(1:$K12) * GM1I * RJ2(1:$K12) + 0.5 * Q X (J,K13:$K12,1) * (U3(1:$K12) ** 2 + V2(1:$K12) ** 2) 55 CONTINUE d call trace_exit (2010) ELSE C 75 CONTINUE C C FOR SUPERSONIC INFLOW ALL VARIABLES FIXED TO FREE STREAM C J = 1 K14 = (KMAX + 31) / 32 K8 = MAX0 (MIN0 (31, KMAX - 1) + 1, 1) d call trace_enter (2011) CDOALL 80 K3=0,32*K14-32,32 DOUBLE PRECISION RJJ1(K8) INTEGER K15, K16 K16 = K3 + 1 K15 = MIN0 (K3 + 1 + 32 - 1, KMAX) - K3 RJJ1(1:$K15) = 1. / XYJ(J,K16:$K15) Q(J,K16:$K15,1) = RHOINF * RJJ1(1:$K15) Q(J,K16:$K15,2) = UINF * Q(J,K16:$K15,1) Q(J,K16:$K15,3) = VINF * Q(J,K16:$K15,1) Q(J,K16:$K15,4) = EINF * RJJ1(1:$K15) 80 CONTINUE C C FOR SUPERSONIC OUTFLOW EXTRAPOLATE ALL VARIABLES C d call trace_exit (2012) J = JMAX K17 = (KMAX + 31) / 32 K8 = MAX0 (MIN0 (31, KMAX - 1) + 1, 1) d call trace_enter (2013) CDOALL 85 K4=0,32*K17-32,32 DOUBLE PRECISION RJJ2(K8) INTEGER K18, K19 K19 = K4 + 1 K18 = MIN0 (K4 + 1 + 32 - 1, KMAX) - K4 RJJ2(1:$K18) = XYJ(J-1,K19:$K18) / XYJ(J,K19:$K18) Q(J,K19:$K18,1) = Q(J-1,K19:$K18,1) * RJJ2(1:$K18) Q(J,K19:$K18,2) = Q(J-1,K19:$K18,2) * RJJ2(1:$K18) Q(J,K19:$K18,3) = Q(J-1,K19:$K18,3) * RJJ2(1:$K18) Q(J,K19:$K18,4) = Q(J-1,K19:$K18,4) * RJJ2(1:$K18) 85 CONTINUE d call trace_exit (2014) END IF C 100 CONTINUE C C TOP BOUNDARY C C ***** FAR-FIELD STUFF ***** C C TOP BOUNDARY C C AINF = SQRT(GAMMA*PINF/RHOINF) HSTFS = 1./GAMI + 0.5*FSMACH**2 GM1I = 1./GAMI GI = 1./GAMMA K = KMAX IF(FSMACH.LT.1.0)THEN C C SUBSONIC FREESTREAM C J13 = (JM + 30) / 32 J7 = MAX0 (MIN0 (31, JM - 2) + 1, 1) d call trace_enter (2015) CDOALL 60 J3=0,32*J13-32,32 DOUBLE PRECISION UF1(J7), VF1(J7), AF21(J7), AF1(J7), SNORM3(J7), X XY3H1(J7), XY4H1(J7), RHOEXT3(J7), RJM13(J7), UEXT3(J7), VEXT3( X J7), EEXT3(J7), PEXT3(J7), R11(J7), R21(J7), QN3(J7), CSPE3(J7) X , C23(J7), QT3(J7), ENTRO3(J7), U4(J7), V3(J7), PRES3(J7), RJJ3( X J7) INTEGER J14, J15 J15 = J3 + 2 J14 = MIN0 (J3 + 2 + 32 - 1, JM) - J3 - 1 UF1(1:$J14) = UINF VF1(1:$J14) = VINF AF21(1:$J14) = GAMI * (HSTFS - .5 * (UF1(1:$J14) ** 2 + VF1(1:$ X J14) ** 2)) AF1(1:$J14) = SQRT (AF21(1:$J14)) C C CHOOSE A REFERENCE FRAME IN TERMS OF NORMAL AND C TANGENTIAL COMPONENTS C C METRIC TERMS C SNORM3(1:$J14) = 1. / SQRT (XY(J15:$J14,K,3) ** 2 + XY(J15:$J14,K, X 4) ** 2) XY3H1(1:$J14) = XY(J15:$J14,K,3) * SNORM3(1:$J14) XY4H1(1:$J14) = XY(J15:$J14,K,4) * SNORM3(1:$J14) C C CHECK FOR INFLOW OR OUTFLOW C FOR INFLOW : THREE VARIABLES ARE SPECIFIED R1 AND C QT ( Q_TANGENTIAL) AND S~ R**GAMMA/P (~ENTROPY) C WITH ONE VAR. COMPUTED. R2 C FOR OUTFLOW : ONE VARIABLE IS FIXED R1 WITH THREE COMPUTED C R2, QT AND S C C COMPUTE RIEMANN INVARIANTS C FIX R1 = QN - 2.*A/(GAMMA-1) AT FREE STREAM C COMPUTE R2 = QN + 2.*A/(GAMMA-1) FROM INTERIOR C EXTRAPOLATION OF FLOW VARIABLES C C GET EXTRAPOLATED VARIABLES C RHOEXT3(1:$J14) = Q(J15:$J14,K-1,1) * XYJ(J15:$J14,K-1) RJM13(1:$J14) = 1. / Q(J15:$J14,K-1,1) UEXT3(1:$J14) = Q(J15:$J14,K-1,2) * RJM13(1:$J14) VEXT3(1:$J14) = Q(J15:$J14,K-1,3) * RJM13(1:$J14) EEXT3(1:$J14) = Q(J15:$J14,K-1,4) * XYJ(J15:$J14,K-1) PEXT3(1:$J14) = GAMI * (EEXT3(1:$J14) - 0.5 * RHOEXT3(1:$J14) * ( X UEXT3(1:$J14) ** 2 + VEXT3(1:$J14) ** 2)) C C SET RIEMANN INVARIANTS C R11(1:$J14) = XY3H1(1:$J14) * UF1(1:$J14) + XY4H1(1:$J14) * VF1(1 X :$J14) - 2. * AF1(1:$J14) * GM1I R21(1:$J14) = XY3H1(1:$J14) * UEXT3(1:$J14) + XY4H1(1:$J14) * X VEXT3(1:$J14) + 2. * SQRT (GAMMA * PEXT3(1:$J14) / RHOEXT3(1:$ X J14)) * GM1I C QN3(1:$J14) = (R11(1:$J14) + R21(1:$J14)) * 0.5 CSPE3(1:$J14) = (R21(1:$J14) - R11(1:$J14)) * GAMI * 0.25 C23(1:$J14) = CSPE3(1:$J14) ** 2 C C SET OTHER FIXED OR EXTRAPOLATED VARIABLES C WHERE (QN3(1:$J14) .LE. 0.0) QT3(1:$J14) = XY4H1(1:$J14) * UF1(1:$J14) - XY3H1(1:$J14) * X VF1(1:$J14) ENTRO3(1:$J14) = GAMMA OTHERWISE QT3(1:$J14) = XY4H1(1:$J14) * UEXT3(1:$J14) - XY3H1(1:$J14) * X VEXT3(1:$J14) ENTRO3(1:$J14) = RHOEXT3(1:$J14) ** GAMMA / PEXT3(1:$J14) END WHERE C C COMPUTE FLOW VARIABLES C U4(1:$J14) = XY3H1(1:$J14) * QN3(1:$J14) + XY4H1(1:$J14) * QT3(1:$ X J14) V3(1:$J14) = XY4H1(1:$J14) * QN3(1:$J14) - XY3H1(1:$J14) * QT3(1:$ X J14) C Q(J15:$J14,K,1) = (C23(1:$J14) * ENTRO3(1:$J14) * GI) ** GM1I PRES3(1:$J14) = C23(1:$J14) * Q(J15:$J14,K,1) * GI C C ADD JACOBIAN C RJJ3(1:$J14) = 1. / XYJ(J15:$J14,K) Q(J15:$J14,K,1) = Q(J15:$J14,K,1) * RJJ3(1:$J14) Q(J15:$J14,K,2) = Q(J15:$J14,K,1) * U4(1:$J14) Q(J15:$J14,K,3) = Q(J15:$J14,K,1) * V3(1:$J14) Q(J15:$J14,K,4) = PRES3(1:$J14) * GM1I * RJJ3(1:$J14) + 0.5 * X Q(J15:$J14,K,1) * (U4(1:$J14) ** 2 + V3(1:$J14) ** 2) 60 CONTINUE d call trace_exit (2016) ELSE C C SUPERSONIC FREESTREAM C K = KMAX d call trace_enter (2017) DO 70 J = 2,JM C C CHOOSE A REFERENCE FRAME IN TERMS OF NORMAL AND C TANGENTIAL COMPONENTS C C C METRIC TERMS C C USE EXTRAPOLATED QN TO DETERMINE INFLOW/OUTFLOW C RHOINV = 1./Q(J,K-1,1) U = Q(J,K-1,2)*RHOINV V = Q(J,K-1,3)*RHOINV SNORM = 1./SQRT(XY(J,K-1,3)**2+XY(J,K-1,4)**2) XY3H = XY(J,K-1,3)*SNORM XY4H = XY(J,K-1,4)*SNORM QN1EXT = XY3H*U + XY4H*V RHOINV = 1./Q(J,K-2,1) U = Q(J,K-2,2)*RHOINV V = Q(J,K-2,3)*RHOINV SNORM = 1./SQRT(XY(J,K-2,3)**2+XY(J,K-2,4)**2) XY3H = XY(J,K-2,3)*SNORM XY4H = XY(J,K-2,4)*SNORM QN2EXT = XY3H*U + XY4H*V QNEXT = QN1EXT C IF(QNEXT.GT.0.)THEN RMET1 = XYJ(J,K-1)/XYJ(J,K) RMET2 = XYJ(J,K-2) / XYJ(J,K) CDOALL 65 N=1,4 Q(J,K,N) = Q(J,K-1,N)*RMET1 65 CONTINUE ENDIF 70 CONTINUE d call trace_exit (2018) ENDIF C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C******************** DEFAULT GRID FOR ********************************* C**************** BICONVEX AIFOIL CALCULATIONS ************************* C*********************************************************************** SUBROUTINE BICONG(JDIM,KDIM,NX,NY,J1,J2,DY,THICK1,TOP,RIGHT,ALEFT > ,X,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM) COMMON/WORKSP/CLUS(MAXJ), WORK2(MAXJ,91) INTEGER J3, I2, J4, K1, J5, J6, J7, J8, I3, I4, I5, I6, J11, K2, X K3, K4, K5 DOUBLE PRECISION XJ1 (:) C C C J REFERS TO THE X ( STREAMWISE ) DIRECTION C K REFERS TO THE Y DIRECTION C C C POINTS IN X DIRECTION NJ = NX C POINTS IN Y DIRECTION NK = NY C C THICKNESS OF AIRFOIL ( 10% PARABOLIC BICONVEX ) THICK = THICK1 C C BELOW, REGION IS [ -X1, X2 ] X [ 0, Y1 ] C AND AIRFOIL IS [ 0, 1] ON THE X AXIS C -X1 IS LEFT BOUNDARY OF REGION X1 = -ALEFT C X2 IS LENGTH OF REGION TO THE RIGHT OF THE AIRFOIL X2 = RIGHT - 1. C TOP OF GRID, NUMBER OF CHORD LENGTHS PERP. TO C STREAMWISE DIRECTION Y1 = TOP C C DEFINES THE FINE MESH C NJM1=NJ-1 NKM1=NK-1 NJM2=NJ-2 NKM2=NK-2 C C NUMBER OF POINTS ON AIRFOIL C JBODY=J2-J1+1 PI=3.141592654 C C PUT THE POINTS ON THE AIRFOIL C DANG=PI/JBODY C TWOTH=2.*THICK J5 = (J2 - J1 + 32) / 32 J8 = MAX0 (MIN0 (J1 + 31, J2) - J1 + 1, 1) d call trace_enter (2019) CDOALL 1 J3=0,32*J5-32,32 DOUBLE PRECISION XI1(J8), DEV1(J8) INTEGER J9, J10 J10 = J1 + J3 J9 = MIN0 (J1 + J3 + 32 - 1, J2) - (J1 + J3) + 1 XI1(1:$J9) = [J10 - J1 :$ J9] + .5 C THIS USES COS ON [ 0, PI ] DEV1(1:$J9) = COS (XI1(1:$J9) * DANG) C X(J,1)=.5*(1.-DEV) CJB CJB THIS ONE GIVES A LOT OF STRETCHING CJB X(1,J)=.5*(1.-DEV) CJB C THIS GIVES LESS STRETCHING X(J10:$J9,1) = .5 * (.5 * (1. - DEV1(1:$J9)) + XI1(1:$J9) X / JBODY) C 1 CONTINUE CJB CJB PUT THE POINTS BEFORE THE AIRFOIL CJB d call trace_exit (2020) DX = X(J1,1) CALL CLUSTR(JDIM,ALEFT,0.,DX,J1,CLUS) I3 = (J1 + 30) / 32 d call trace_enter (2021) CDOALL 2 I2=0,32*I3-32,32 INTEGER I7, I8, JJ I8 = I2 + 2 I7 = MIN0 (I2 + 2 + 32 - 1, J1) - I2 - 1 X(J1-I8+1:$I7:(-1),1) = -CLUS(I8:$I7) 2 CONTINUE CJB CJB PUT THE POINTS AFTER THE AIRFOIL CJB d call trace_exit (2022) DX = 1. - X(J2,1) NI = NJ-J2+1 XMM = RIGHT - 1.0 CALL CLUSTR(JDIM,XMM,0.,DX,NI,CLUS) J11 = (NX - J2 + 31) / 32 d call trace_enter (2023) CDOALL 3 J4=0,32*J11-32,32 INTEGER J12, J13, I J13 = J2 + J4 + 1 J12 = MIN0 (J2 + J4 + 1 + 32 - 1, NX) - (J2 + J4) X(J13:$J12,1) = 1.0 + CLUS(J13-J2+1:$J12) 3 CONTINUE C c*$* cv stripmine off d call trace_exit (2024) K2 = NY - 1 ALLOCATE (XJ1(NX)) CDOALL J=1,NX XJ1(J) = X(J,1) END CDOALL d call trace_enter (2025) CDOALL 7 K=2,K2+1 X(1:$NX,K) = XJ1(1:$NX) 7 CONTINUE d call trace_exit (2026) DEALLOCATE (XJ1) CJB CJB PUT THE POINTS ABOVE THE AIRFOIL CJB CALL CLUSTR(JDIM,TOP,0.,DY,NY,CLUS) K3 = (NY + 31) / 32 d call trace_enter (2027) CDOALL 5 K1=0,32*K3-32,32 INTEGER K6, K7 K7 = K1 + 1 K6 = MIN0 (K1 + 1 + 32 - 1, NY) - K1 Y(1,K7:$K6) = CLUS(K7:$K6) 5 CONTINUE d call trace_exit (2028) J6 = NX - 1 d call trace_enter (2029) CDOALL 6 K=1,NY Y(2:$J6,K) = Y(1,K) 6 CONTINUE C d call trace_exit (2030) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C********************** METRIC CALCULATIONS **************************** C*********************************************************************** SUBROUTINE CALCME(JDIM,KDIM,XY,XIT,ETT,XYJ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK COMMON/IO/IREAD,IWRIT,IVER C DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM) INTEGER K1, J1 C C FORM JACOBIAN C d call trace_enter (2031) DO 30 K=KBEGIN,KEND DO 30 J=JBEGIN,JEND DINV = 1./ ( XY(J,K,4) * XY(J,K,1) - XY(J,K,3) * XY(J,K,2) ) IF(DINV .LE. 0.0)WRITE(IWRIT,*)'JACOBIAN( ',J,K,' )= ',DINV XYJ(J,K) = DINV 30 CONTINUE d call trace_exit (2032) K1 = KEND - KBEGIN + 1 J1 = JEND - JBEGIN + 1 C C CALCULATE METRICS C d call trace_enter (2033) CDOALL 32 K=KBEGIN,KBEGIN+K1-1 DOUBLE PRECISION DINV2(J1) DINV2(1:$J1) = XYJ(JBEGIN:$J1,K) C C XIX = XY1, XIY=XY2,ETAX=XY3, ETAY = XY4 C XY(JBEGIN:$J1,K,1) = XY(JBEGIN:$J1,K,1) * DINV2(1:$J1) XY(JBEGIN:$J1,K,2) = -XY(JBEGIN:$J1,K,2) * DINV2(1:$J1) XY(JBEGIN:$J1,K,3) = -XY(JBEGIN:$J1,K,3) * DINV2(1:$J1) XY(JBEGIN:$J1,K,4) = XY(JBEGIN:$J1,K,4) * DINV2(1:$J1) 32 CONTINUE C d call trace_exit (2034) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************* CALCULATE PRESSURE AND SOUND SPEED ********************* C*********************************************************************** SUBROUTINE CALCPS(JDIM,KDIM,Q,PRESS,SNDSP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) INTEGER K1, J1 K1 = KEND - KBEGIN + 1 J1 = JEND - JBEGIN + 1 C C CALCULATE PRESSURE AND SOUND SPEED AND STORE C C NOTE PRESSURE HAS JACOBIAN C SOUND SPEED DOES NOT HAVE JACOBIAN C CALL AFTER BC ROUTINE SO THAT BOUNDARY VALUES ARE CORRECT C d call trace_enter (2035) CDOALL 1 K=KBEGIN,KBEGIN+K1-1 C PRESS(JBEGIN:$J1,K) = GAMI * (Q(JBEGIN:$J1,K,4) - 0.5 * (Q( X JBEGIN:$J1,K,2) ** 2 + Q(JBEGIN:$J1,K,3) ** 2) / Q(JBEGIN:$J1,K, X 1)) SNDSP(JBEGIN:$J1,K) = SQRT (GAMMA * PRESS(JBEGIN:$J1,K) / Q( X JBEGIN:$J1,K,1)) C 1 CONTINUE C d call trace_exit (2036) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C !********** SUBROUTINE CLUSTR SUBROUTINE CLUSTR(JDIM,YMAX,YMIN,DYM,LMAX,R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(JDIM) C C C EPS = EPSIL(YMAX,YMIN,DYM,LMAX,0.001,100,1,2) K = -1 R(1) = YMIN d call trace_enter (2037) DO 5 L = 2,LMAX K = K+1 GN = (1.0 +EPS/SQRT(FLOAT(K+4)))**K R(L) = R(L-1) + DYM * GN 5 CONTINUE C d call trace_exit (2038) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************** COMPUTE NONLINEAR 2ND AND 4TH ORDER COEFFICIENTS ***** C*********************************************************************** SUBROUTINE COEF24(JDIM,KDIM,COEF2X,COEF4X,COEF2Y,COEF4Y, > PCOEF,SPECT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION COEF2X(JDIM,KDIM),COEF4X(JDIM,KDIM) DIMENSION COEF2Y(JDIM,KDIM),COEF4Y(JDIM,KDIM) DIMENSION PCOEF(JDIM,KDIM),SPECT(JDIM,KDIM) INTEGER J1, K1, J2, J3, J4, J5, J6 K1 = KEND - KBEGIN + 1 J2 = JEND - JBEGIN + 1 C C COEF2X,Y COMES IN AS COEFX,Y (PRESSURE GRADIANT COEFFICIENT) C COEF4X,Y COMES IN AS SPECX,Y (SPECTRAL RADIUS) C d call trace_enter (2039) CDOALL 1 K=KBEGIN,KBEGIN+K1-1 PCOEF(JBEGIN:$J2,K) = COEF2Y(JBEGIN:$J2,K) SPECT(JBEGIN:$J2,K) = COEF4Y(JBEGIN:$J2,K) 1 CONTINUE C d call trace_exit (2040) EPS4Y = DIS4Y/64. EPS2Y = DIS2Y K1 = KUP - KBEGIN + 1 J2 = JUP - JLOW + 1 C C FORM 2ND AND 4TH ORDER DISSIPATION COEFFICIENTS IN Y C d call trace_enter (2041) CDOALL 10 K=KBEGIN,KBEGIN+K1-1 DOUBLE PRECISION FIL1(J2), C21(J2), C41(J2) C FIL1(1:$J2) = SPECT(JLOW:$J2,K+1) + SPECT(JLOW:$J2,K) C21(1:$J2) = PCOEF(JLOW:$J2,K) * FIL1(1:$J2) * EPS2Y C41(1:$J2) = EPS4Y * FIL1(1:$J2) C41(1:$J2) = C41(1:$J2) - MIN (C41(1:$J2), C21(1:$J2)) COEF2Y(JLOW:$J2,K) = C21(1:$J2) COEF4Y(JLOW:$J2,K) = C41(1:$J2) 10 CONTINUE d call trace_exit (2042) K1 = KEND - KBEGIN + 1 J2 = JEND - JBEGIN + 1 C C d call trace_enter (2043) CDOALL 2 K=KBEGIN,KBEGIN+K1-1 PCOEF(JBEGIN:$J2,K) = COEF2X(JBEGIN:$J2,K) SPECT(JBEGIN:$J2,K) = COEF4X(JBEGIN:$J2,K) 2 CONTINUE C d call trace_exit (2044) EPS4X = DIS4X/64. EPS2X = DIS2X K1 = KUP - KLOW + 1 C C FORM 2ND AND 4TH ORDER DISSIPATION COEFFICIENTS IN Y C J3 = (JUP - JBEGIN + 32) / 32 J2 = JUP - JBEGIN + 1 CDOALL J1=0,32*J3-32,32 INTEGER J7, J8, JPL J8 = JBEGIN + J1 J7 = MIN0 (JBEGIN + J1 + 32 - 1, JUP) - (JBEGIN + J1) + 1 END CDOALL d call trace_enter (2045) CDOALL 20 K=KLOW,KLOW+K1-1 DOUBLE PRECISION FIL2(J2), C22(J2), C42(J2) C FIL2(1:$J2) = SPECT(JPLUS(JBEGIN:$J2),K) + SPECT(JBEGIN:$J2,K) C22(1:$J2) = PCOEF(JBEGIN:$J2,K) * FIL2(1:$J2) * EPS2X C42(1:$J2) = EPS4X * FIL2(1:$J2) C42(1:$J2) = C42(1:$J2) - MIN (C42(1:$J2), C22(1:$J2)) COEF2X(JBEGIN:$J2,K) = C22(1:$J2) COEF4X(JBEGIN:$J2,K) = C42(1:$J2) 20 CONTINUE d call trace_exit (2046) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************** COMPUTE EIGENVALUES *********************************** C*********************************************************************** SUBROUTINE EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ,XY,TMET, * UU,CC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C IX = 1 IY = 2 FOR XI DIRECTION METRICS C IX = 3 IY = 4 FOR ETA DIRECTION METRICS C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION XY(JDIM,KDIM,4) DIMENSION TMET(JDIM,KDIM),XYJ(JDIM,KDIM) C DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM) INTEGER K1, J1 K1 = KEND - KBEGIN + 1 J1 = JEND - JBEGIN + 1 C C COMPUTE SPECTRAL RADIUS SCALING C d call trace_enter (2047) CDOALL 100 K=KBEGIN,KBEGIN+K1-1 DOUBLE PRECISION RI1(J1), U1(J1), V1(J1) C RI IS J/RHO RI1(1:$J1) = 1. / Q(JBEGIN:$J1,K,1) U1(1:$J1) = Q(JBEGIN:$J1,K,2) * RI1(1:$J1) V1(1:$J1) = Q(JBEGIN:$J1,K,3) * RI1(1:$J1) C C NOTE: TIME METRICS PUT IN (SLOPPY BUT THEY WORK) C UU(JBEGIN:$J1,K) = TMET(JBEGIN:$J1,K) + U1(1:$J1) * XY(JBEGIN:$ X J1,K,IX) + V1(1:$J1) * XY(JBEGIN:$J1,K,IY) CC(JBEGIN:$J1,K) = SNDSP(JBEGIN:$J1,K) * SQRT (XY(JBEGIN:$J1,K, X IX) ** 2 + XY(JBEGIN:$J1,K,IY) ** 2) C 100 CONTINUE C d call trace_exit (2048) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C DOUBLE PRECISION FUNCTION EPSIL * (FMX,FMIN,DFM,NPT,FPCC,ICC,NCALL,KEY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/IO/IREAD,IWRIT,IVER INTEGER M1 C C FMXL = FMX FMINL = FMIN DFML = DFM FPCCL = FPCC ICCL = ICC C C SLOWER STRETCHING IF KEY = 2 C IF (KEY .NE. 2) THEN C FNPTM2 = NPT-2 IF(NCALL.EQ.1) EPS = (FMXL/DFML)**(1./FNPTM2)-1.0 C d call trace_enter (2049) DO 3 NIT = 1,ICCL EP1 = EPS+1. EP1TN = EP1**FNPTM2 DFMOE = DFML/EPS F = FMXL-FMINL-DFMOE*(EP1TN*EP1-1.0) IF(ABS(F).LT.FPCCL) GO TO 4 FPN = DFMOE*FNPTM2*EP1TN EPS = EPS+F/FPN 3 CONTINUE d call trace_exit (2050) ELSE 50 CONTINUE IF(NCALL.EQ.1) EPS=20. d call trace_enter (2051) DO 11 NIT=1,ICCL SUMG = 0.0 SUMGP = 0.0 M1 = NPT - 1 CDOACROSS M=2,M1+1 DOUBLE PRECISION GPN, GN INTEGER N GN = (1.0 + EPS / SQRT (FLOAT (M - 2 + 4))) ** (M - 2) GPN = (M - 2) * (1.0 + EPS / SQRT (FLOAT (M - 2 + 4))) ** (M - 2 - X 1) * (1.0 / SQRT (FLOAT (M - 2 + 4))) CALL AWAIT (1, 1, SUMG) SUMG = SUMG + GN CALL ADVANCE (1, SUMG) 12 CALL AWAIT (2, 1, SUMGP) SUMGP = SUMGP + GPN CALL ADVANCE (2, SUMGP) END CDOACROSS F = FMXL-FMINL-DFML*SUMG IF(ABS(F).LT.FPCCL) GO TO 4 FP=-DFML*SUMGP 11 EPS = EPS-F/FP d call trace_exit (2052) END IF C 5 CONTINUE EPSIL = EPS WRITE(IWRIT,100) RETURN 4 EPSIL = EPS WRITE(IWRIT,101) EPSIL,F,NIT RETURN C 100 FORMAT(/40HEXCEEDED MAX NO. OF ITERATIONS IN EPSIL.) 101 FORMAT(/6HEPSIL=,F12.5,5X,7H AND F=,F12.5,5X,6HAFTER ,I3, * 12H ITERATIONS.) C END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C********************** DIFFERENCE X, Y IN XI ************************* C*********************************************************************** SUBROUTINE ETADIF(JDIM,KDIM,X,Y,XY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION XY(JDIM,KDIM,4) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM) INTEGER K1, K2, K3, K4, K5 C C XY2 = XETA, XY1 = YETA C d call trace_enter (2053) DO 21 J=JBEGIN,JEND K2 = (KUP - KLOW + 32) / 32 CDOALL 20 K1=0,32*K2-32,32 INTEGER K6, K7 K7 = KLOW + K1 K6 = MIN0 (KLOW + K1 + 32 - 1, KUP) - (KLOW + K1) + 1 XY(J,K7:$K6,2) = (X(J,K7+1:$K6) - X(J,K7-1:$K6)) * .5 XY(J,K7:$K6,1) = (Y(J,K7+1:$K6) - Y(J,K7-1:$K6)) * .5 20 CONTINUE K = KBEGIN XY(J,K,2) = (-3. * X(J,K) + 4. * X(J,K+1) - X(J,K+2)) * .5 XY(J,K,1) = (-3. * Y(J,K) + 4. * Y(J,K+1) - Y(J,K+2)) * .5 K = KEND XY(J,K,2) = (3. * X(J,K) - 4. * X(J,K-1) + X(J,K-2)) * .5 XY(J,K,1) = (3. * Y(J,K) - 4. * Y(J,K-1) + Y(J,K-2)) * .5 21 CONTINUE C C d call trace_exit (2054) RETURN END C*********************************************************************** C*** SECOND - FOURTH : ORDER EXPLICIT FILTER ************************* C*********************************************************************** SUBROUTINE FILERX(JDIM,KDIM,Q,S,XYJ,COEF2,COEF4,WORK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM) C DIMENSION WORK(JDIM,KDIM,3) C C FOURTH ORDER SMOOTHING, ADDED EXPLICITLY TO RHS C SECOND ORDER NEAR SHOCKS WITH PRESSURE GRD COEFF. C C XI DIRECTION C C START DIFFERENCES EACH VARIABLE SEPARATELY C call trace_enter (2055) DTD = DT / (1. + PHIDT) CDOALL 16 K = KLOW,KUP double precision WORK(JDIM,KDIM,3),N DO 15 N = 1,4 WORK(JLOW:JUP,K,1) = Q(JLOW+1:JUP+1,K,N)*XYJ(JLOW+1:JUP+1,K) * - Q(JLOW:JUP,K,N)*XYJ(JLOW:JUP,K) C C C MESH BC C IF(.NOT.PERIDC )THEN WORK(JBEGIN,K,1) = Q(JBEGIN+1,K,N)*XYJ(JBEGIN+1,K) - * Q(JBEGIN,K,N)*XYJ(JBEGIN,K) WORK(JEND,K,1) = WORK(JEND-1,K,1) ENDIF WORK(JLOW:JUP,K,2) = WORK(JLOW+1:JUP+1,K,1) - * 2.* WORK(JLOW:JUP,K,1) + WORK(JLOW-1:JUP-1,K,1) IF(.NOT.PERIDC )THEN WORK(JBEGIN,K,2) = Q(JBEGIN+2,K,N)*XYJ(JBEGIN+2,K) - * 2.*Q(JBEGIN+1,K,N)*XYJ(JBEGIN+1,K) + * Q(JBEGIN,K,N)*XYJ(JBEGIN,K) WORK(JEND,K,2) = 0. ENDIF WORK(JLOW:JUP,K,3) = (COEF2(JLOW:JUP,K)*WORK(JLOW:JUP,K,1) - * COEF4(JLOW:JUP,K)*WORK(JLOW:JUP,K,2)) IF(.NOT.PERIDC )THEN WORK(JBEGIN,K,3) = (COEF2(JBEGIN,K)*WORK(JBEGIN,K,1) - > COEF4(JBEGIN,K)*WORK(JBEGIN,K,2)) WORK(JEND,K,3) = 0. ENDIF 98 S(JLOW:JUP,K,N) = S(JLOW:JUP,K,N) + (WORK(JLOW:JUP,K,3) - * WORK(JLOW-1:JUP-1,K,3))*DTD 15 CONTINUE 16 CONTINUE call trace_exit (2056) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C*** SECOND - FOURTH : ORDER EXPLICIT FILTER ************************* C*********************************************************************** SUBROUTINE FILERY(JDIM,KDIM,Q,S,XYJ,COEF2,COEF4,WORK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM) C DIMENSION WORK(JDIM,KDIM,3) INTEGER J1, K1, J2, J3, J4, J5 C C FOURTH ORDER SMOOTHING, ADDED EXPLICITLY TO RHS C SECOND ORDER NEAR SHOCKS WITH PRESSURE GRD COEFF. C C ETA DIRECTION C C START DIFFERENCES EACH VARIABLE SEPARATLY C d call trace_enter (2057) DTD = DT / (1. + PHIDT) DO 39 N = 1,4 C CDOALL 35 K=KBEGIN,KUP INTEGER KPL WORK(JLOW:JUP,K,1) = Q(JLOW:JUP,K+1,N) * XYJ(JLOW:JUP,K+1) - Q( X JLOW:JUP,K,N) * XYJ(JLOW:JUP,K) 35 CONTINUE C C FOR FOURTH ORDER C CDOALL 36 K=KLOW,KUP INTEGER KMI, KPL WORK(JLOW:JUP,K,2) = WORK(JLOW:JUP,K+1,1) - 2. * WORK(JLOW:JUP,K,1 X ) + WORK(JLOW:JUP,K-1,1) 36 CONTINUE C C BOUNDARY C J3 = (JUP - JLOW + 32) / 32 CDOALL 37 J1=0,32*J3-32,32 INTEGER J6, J7 J7 = JLOW + J1 J6 = MIN0 (JLOW + J1 + 32 - 1, JUP) - (JLOW + J1) + 1 WORK(J7:$J6,KBEGIN,2) = 0. WORK(J7:$J6,KUP,2) = WORK(J7:$J6,KUP-1,1) - WORK(J7:$J6,KUP,1) 37 CONTINUE C C C FORM DISSIPATION TERM C CDOALL 38 K=KBEGIN,KUP WORK(JLOW:JUP,K,3) = (COEF2(JLOW:JUP,K) * WORK(JLOW:JUP,K,1) - X COEF4(JLOW:JUP,K) * WORK(JLOW:JUP,K,2)) 38 CONTINUE C C ADD IN DISSIPATION C C C CDOALL 40 K=KLOW,KUP S(JLOW:JUP,K,N) = S(JLOW:JUP,K,N) + (WORK(JLOW:JUP,K,3) - WORK( X JLOW:JUP,K-1,3)) * DTD 40 CONTINUE C 39 CONTINUE C d call trace_exit (2058) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C*********************************************************************** C*** GRADIENT COEFFICENT FOR ARTIFICIAL DISSIPATION ******************** C*********************************************************************** SUBROUTINE GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF,PRSS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C IXY = 1 FOR XI DIRECTION C IXY = 2 FOR ETA DIRECTION COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION PRESS(JDIM,KDIM),XYJ(JDIM,KDIM) DIMENSION COEF(JDIM,KDIM),PRSS(JDIM,KDIM) INTEGER J1, K1, J2, J3, K2, J4, K3, J5, J6, J7, J8, J9, K4, K5, K6 X , J12, J15, K9, J18 C IF(IXY.EQ.1)THEN C JEDGE1 = JLOW JEDGE2 = JUP K3 = KEND - KBEGIN + 1 J5 = JEND - JBEGIN + 1 C C COMPUTE PRESSURE GRAD C d call trace_enter (2059) CDOALL 1 K=KBEGIN,KBEGIN+K3-1 COEF(JBEGIN:$J5,K) = 0. PRSS(JBEGIN:$J5,K) = PRESS(JBEGIN:$J5,K) * XYJ(JBEGIN:$J5,K) 1 CONTINUE C C BYPASS LOGIC IF DIS2X = 0.0 C d call trace_exit (2060) IF(DIS2X.EQ.0.0)RETURN K3 = KUP - KLOW + 1 C J6 = (JUP - JLOW + 32) / 32 J5 = JUP - JLOW + 1 CDOALL J1=0,32*J6-32,32 INTEGER J10, J11, JMI, JPL J11 = JLOW + J1 J10 = MIN0 (JLOW + J1 + 32 - 1, JUP) - (JLOW + J1) + 1 END CDOALL d call trace_enter (2061) CDOALL 2 K=KLOW,KLOW+K3-1 DOUBLE PRECISION R11(J5), R21(J5) R11(1:$J5) = PRSS(JPLUS(JLOW:$J5),K) - 2. * PRSS(JLOW:$J5,K) + X PRSS(JMINU(JLOW:$J5),K) R21(1:$J5) = PRSS(JPLUS(JLOW:$J5),K) + 2. * PRSS(JLOW:$J5,K) + X PRSS(JMINU(JLOW:$J5),K) COEF(JLOW:$J5,K) = ABS (R11(1:$J5) / R21(1:$J5)) 2 CONTINUE C C C MESH BC C d call trace_exit (2062) IF(.NOT.PERIDC )THEN K4 = (KUP - KLOW + 32) / 32 d call trace_enter (2063) CDOALL 3 K1=0,32*K4-32,32 INTEGER K7, K8 K8 = KLOW + K1 K7 = MIN0 (KLOW + K1 + 32 - 1, KUP) - (KLOW + K1) + 1 COEF(JBEGIN,K8:$K7) = COEF(JLOW,K8:$K7) COEF(JEND,K8:$K7) = COEF(JUP,K8:$K7) PRSS(JBEGIN,K8:$K7) = COEF(JLOW,K8:$K7) PRSS(JEND,K8:$K7) = COEF(JUP,K8:$K7) 3 CONTINUE d call trace_exit (2064) ENDIF K3 = KUP - KLOW + 1 C C SMOOTH OUT PRESSURE COEFFICIENT C J12 = (JUP - JLOW + 32) / 32 J5 = JUP - JLOW + 1 CDOALL J2=0,32*J12-32,32 INTEGER J13, J14, JPLPL, JMI, JPL J14 = JLOW + J2 J13 = MIN0 (JLOW + J2 + 32 - 1, JUP) - (JLOW + J2) + 1 END CDOALL d call trace_enter (2065) CDOALL 4 K=KLOW,KLOW+K3-1 PRSS(JLOW:$J5,K) = MAX (COEF(JMINU(JLOW:$J5),K), COEF(JLOW:$J5,K) X , COEF(JPLUS(JLOW:$J5),K), COEF(JPLUS(JPLUS(JLOW:$J5))+((-1)+K)* X JDIM,1)) 4 CONTINUE d call trace_exit (2066) K5 = KUP - KLOW + 1 C J15 = (JUP - JLOW + 32) / 32 J5 = JUP - JLOW + 1 CDOALL J3=0,32*J15-32,32 INTEGER J16, J17, JMI, JPL J17 = JLOW + J3 J16 = MIN0 (JLOW + J3 + 32 - 1, JUP) - (JLOW + J3) + 1 END CDOALL d call trace_enter (2067) CDOALL 5 K=KLOW,KLOW+K5-1 COEF(JLOW:$J5,K) = 0.5 * PRSS(JLOW:$J5,K) + 0.25 * (PRSS(JPLUS( X JLOW:$J5),K) + PRSS(JMINU(JLOW:$J5),K)) 5 CONTINUE C d call trace_exit (2068) IF(.NOT.PERIDC )THEN C TURN OFF SECOND ORDER NEAR TRAILING EDGE C C C MESH BC C K9 = (KUP - KLOW + 32) / 32 d call trace_enter (2069) CDOALL 12 K2=0,32*K9-32,32 INTEGER K10, K11 K11 = KLOW + K2 K10 = MIN0 (KLOW + K2 + 32 - 1, KUP) - (KLOW + K2) + 1 COEF(JBEGIN,K11:$K10) = COEF(JLOW,K11:$K10) COEF(JEND,K11:$K10) = COEF(JUP,K11:$K10) 12 CONTINUE C d call trace_exit (2070) ENDIF C ELSEIF(IXY .EQ. 2)THEN K6 = KEND - KBEGIN + 1 J5 = JEND - JBEGIN + 1 C C C GRADIANT COEFFICENT IN ETA DIRECTION C d call trace_enter (2071) CDOALL 21 K=KBEGIN,KBEGIN+K6-1 COEF(JBEGIN:$J5,K) = 0. PRSS(JBEGIN:$J5,K) = PRESS(JBEGIN:$J5,K) * XYJ(JBEGIN:$J5,K) 21 CONTINUE C C BYPASS LOGIC IF DIS2Y = 0.0 C d call trace_exit (2072) IF(DIS2Y.EQ.0.0)RETURN K3 = KUP - KLOW + 1 J7 = JUP - JLOW + 1 C d call trace_enter (2073) CDOALL 22 K=KLOW,KLOW+K3-1 DOUBLE PRECISION R12(J7), R22(J7) INTEGER KMI, KPL R12(1:$J7) = PRSS(JLOW:$J7,K+1) - 2. * PRSS(JLOW:$J7,K) + PRSS( X JLOW:$J7,K-1) R22(1:$J7) = PRSS(JLOW:$J7,K+1) + 2. * PRSS(JLOW:$J7,K) + PRSS( X JLOW:$J7,K-1) COEF(JLOW:$J7,K) = ABS (R12(1:$J7) / R22(1:$J7)) 22 CONTINUE C C BC C d call trace_exit (2074) J18 = (JUP - JLOW + 32) / 32 d call trace_enter (2075) CDOALL 24 J4=0,32*J18-32,32 INTEGER J19, J20 J20 = JLOW + J4 J19 = MIN0 (JLOW + J4 + 32 - 1, JUP) - (JLOW + J4) + 1 COEF(J20:$J19,KBEGIN) = COEF(J20:$J19,KLOW) COEF(J20:$J19,KEND) = COEF(J20:$J19,KUP) PRSS(J20:$J19,KBEGIN) = COEF(J20:$J19,KLOW) PRSS(J20:$J19,KEND) = COEF(J20:$J19,KUP) 24 CONTINUE d call trace_exit (2076) K5 = KUP - KLOW + 1 J7 = JUP - JLOW + 1 C C SMOOTH OUT PRESSURE COEFFICIENT C d call trace_enter (2077) CDOALL 26 K=KLOW,KLOW+K5-1 INTEGER KMI, KPL PRSS(JLOW:$J7,K) = MAX (COEF(JLOW:$J7,K-1), COEF(JLOW:$J7,K), COEF X (JLOW:$J7,K+1)) 26 CONTINUE d call trace_exit (2078) K6 = KUP - KLOW + 1 J8 = JUP - JLOW + 1 C d call trace_enter (2079) CDOALL 28 K=KLOW,KLOW+K6-1 INTEGER KMI, KPL COEF(JLOW:$J8,K) = 0.5 * PRSS(JLOW:$J8,K) + 0.25 * (PRSS(JLOW:$J8, X K+1) + PRSS(JLOW:$J8,K-1)) 28 CONTINUE C d call trace_exit (2080) ENDIF C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C****************** INITIALIZE VARIABLES AND DATA ********************** C*********************************************************************** SUBROUTINE INITIA(JDIM,KDIM,Q,PRESS,SNDSP,S,XY,XYJ, * XIT,ETT,DS,X,Y,SLOPE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4) DIMENSION PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION XIT(JDIM,KDIM),ETT(JDIM,KDIM),DS(JDIM,KDIM) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM),SLOPE(JDIM) COMMON/INFNTY/RHOINF,UINF,VINF,EINF,PINF INTEGER J1, J2, J3, J4, J5, J6, J7, J10, J13 C C C COMPUTE VARIOUS CONSTANTS C GAMMA = 1.4 GAMI = GAMMA -1. PI = 4.*ATAN(1.0) C COSANG = COS( PI*ALPHA/180.) SINANG = SIN(PI*ALPHA/180.) RHOINF = 1. UINF = FSMACH*COSANG VINF = FSMACH*SINANG EINF = 1./(GAMMA*GAMI) + .5*FSMACH**2 PINF = 1./GAMMA C C GRID CALCULATION C CALL BICONG(JDIM,KDIM,JMAX,KMAX,JTAIL1,JTAIL2,DYM,THICK, 1 YMAX,XMAX,XMIN,X,Y) C C INITIALIZE TO FREE STREAM C C LOAD FREE STREAM AT ANGLE OF ATTACK ALPHA C d call trace_enter (2081) CDOALL 32 K=1,KMAX Q(1:$JMAX,K,1) = RHOINF Q(1:$JMAX,K,2) = RHOINF * UINF Q(1:$JMAX,K,3) = RHOINF * VINF Q(1:$JMAX,K,4) = EINF 32 CONTINUE C C C SET UP VARIABLES, INDECIES, AND COMPUTE METRICS C C C JTAIL2 = JMAX FOR PERIDC MESHES C d call trace_exit (2082) IF(PERIDC )JTAIL2 = JMAX C C SET PERIDC INDEX C JM = JMAX-1 KM = KMAX-1 JBEGIN = 1 JEND = JMAX KBEGIN = 1 KEND = KMAX KLOW = 2 KUP = KM J4 = (JMAX + 31) / 32 d call trace_enter (2083) CDOALL 10 J1=0,32*J4-32,32 INTEGER J8, J9 J9 = J1 + 1 J8 = MIN0 (J1 + 1 + 32 - 1, JMAX) - J1 JPLUS(J9:$J8) = [J9 + 1 :$ J8] JMINU(J9:$J8) = [J9 - 1 :$ J8] 10 CONTINUE C d call trace_exit (2084) IF ( .NOT.PERIDC ) THEN JPLUS(JMAX) = JMAX JMINU (1) = 1 JLOW = 2 JUP = JMAX-1 ELSE JPLUS(JMAX) = 1 JMINU (1) = JMAX JLOW = 1 JUP = JMAX JTAIL1 = 1 JTAIL2 = JMAX ENDIF C d call trace_enter (2085) CDOALL 8 K=1,KMAX XIT(1:$JMAX,K) = 0. ETT(1:$JMAX,K) = 0. 8 CONTINUE C C C METRIC CALCULATION C S USED FOR WORK SPACE HERE C d call trace_exit (2086) CALL XYMETS(JDIM,KDIM,X,Y,XY,XIT,ETT,XYJ,S) C d call trace_enter (2087) CDOALL 9 N=1,4 INTEGER K DO 9 K=1,KMAX S(1:$JMAX,K,N) = 0. 9 CONTINUE C C INITIALIZE VARIABLE TIME STEP IF JACDT=1 C d call trace_exit (2088) FACTOR=1.0 IF(JACDT.EQ.0) FACTOR=0.0 d call trace_enter (2089) CDOALL 65 K=1,KMAX DS(1:$JMAX,K) = 1.0 / (1.0 + FACTOR * SQRT (XYJ(1:$JMAX,K))) 65 CONTINUE C C C SCALE VARIABLES BY JACOBIAN STORED BY XYMETS IN XYJ(J,K) C d call trace_exit (2090) d call trace_enter (2091) CDOALL 50 N=1,4 INTEGER K DO 50 K=1,KMAX Q(1:$JMAX,K,N) = Q(1:$JMAX,K,N) / XYJ(1:$JMAX,K) 50 CONTINUE C C C AIRFOIL BODY SLOPE C C USED ONLY IF RUNNING THE BICONVEX AIRFOIL C d call trace_exit (2092) J10 = (JMAX + 31) / 32 d call trace_enter (2093) CDOALL 22 J2=0,32*J10-32,32 INTEGER J11, J12 J12 = J2 + 1 J11 = MIN0 (J2 + 1 + 32 - 1, JMAX) - J2 SLOPE(J12:$J11) = 0.0 22 CONTINUE C d call trace_exit (2094) J13 = (JTAIL2 - JTAIL1 + 32) / 32 J7 = MAX0 (MIN0 (JTAIL1 + 31, JTAIL2) - JTAIL1 + 1, 1) d call trace_enter (2095) CDOALL 24 J3=0,32*J13-32,32 DOUBLE PRECISION XPT1(J7), DGDX1(J7) INTEGER J14, J15 J15 = JTAIL1 + J3 J14 = MIN0 (JTAIL1 + J3 + 32 - 1, JTAIL2) - (JTAIL1 + J3) + 1 XPT1(1:$J14) = X(J15:$J14,1) DGDX1(1:$J14) = 2. * THICK * (1. - 2. * XPT1(1:$J14)) SLOPE(J15:$J14) = ATAN (DGDX1(1:$J14)) 24 CONTINUE d call trace_exit (2096) 25 FORMAT(1H ,'****J,SLOPE',I5,2X,E20.8) C C INITIALIZE PRESS,SNDSP C CALL CALCPS(JDIM,KDIM,Q,PRESS,SNDSP) C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C****************** READ IN DATA AND SET PARAMETERS ******************** C*********************************************************************** SUBROUTINE INPUT(JLARGE,JKLARG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK COMMON/IO/IREAD,IWRIT,IVER C C C SET PARAMETER DEFAULTS C JMAX = 65 KMAX = 21 NP = 10000 NPCP = 1 IPRINT = 1 JTAIL1 = 15 JTAIL2 = 51 JACDT = 1 THETAD = 1. PHIDT = 0. PERIDC = .FALSE. FSMACH = 0.8 ALPHA = 0. DIS2X = 0.0 DIS4X = 0.64 DIS2Y = 0.0 DIS4Y = 0.64 DYM = 0.02 XMIN = 2.0 XMAX = 3.0 YMAX = 3.0 THICK = 0.1 NSTEPS = 100 DT = 1.0 C WRITE(IWRIT,*)' ', 1 ' * ****** **** ***** ***** ' WRITE(IWRIT,*)' ', 1 ' * * * * * * * * * ' WRITE(IWRIT,*)' ', 1 ' * * * * * * * *' WRITE(IWRIT,*)' ', 1 ' ******* ****** * * * *' WRITE(IWRIT,*)' ', 1 ' * * * * * * * * ' WRITE(IWRIT,*)' ', 1 ' * * * * **** ****** ***** ' WRITE(IWRIT,*) ' ' C WRITE(IWRIT,*)'ARC2D : AN IMPLICIT FINITE ', 1 'DIFFERENCE CODE FOR AIRFOILS' WRITE(IWRIT,*)'FOR FURTHER DETAILS SEE SELF ', 1 'DOCUMENTATION OR CONTACT ' WRITE(IWRIT,*)'THOMAS H. PULLIAM' WRITE(IWRIT,*)'MS 202A-1, NASA AMES RESEARCH CENTER' WRITE(IWRIT,*)'MOFFETT FIELD, CA., 94035' WRITE(IWRIT,*)'PHONE : 415 - 694 - 6417 ' C C READ IN THE DATA FOR THE RUN FROM STANDARD INPUT C *1* WRITE(IWRIT,*) 'ENTER NSTEPS, DT DEFAULTS ARE ', NSTEPS,DT READ(IREAD,*) NSTEPS,DT C *2* WRITE(IWRIT,*) 'ENTER JMAX, KMAX DEFAULTS ARE ', JMAX, KMAX READ(IREAD,*) JMAX, KMAX C *3* WRITE(IWRIT,*) 'ENTER FSMACH, ALPHA DEFAULTS ARE ', 1 FSMACH, ALPHA READ(IREAD,*) FSMACH, ALPHA C *4* WRITE(IWRIT,*) 'ENTER DIS2X, DIS2X DEFAULTS ARE ', * DIS2X, DIS4X READ(IREAD,*) DIS2X, DIS4X WRITE(IWRIT,*) 'ENTER DIS2Y, DIS4Y DEFAULTS ARE ', * DIS2Y, DIS4Y READ(IREAD,*) DIS2Y, DIS4Y C *5* WRITE(IWRIT,*) 'ENTER PERIDC DEFAULT IS ', PERIDC READ(IREAD,*) PERIDC IF(.NOT.PERIDC )THEN WRITE(IWRIT,*) 'ENTER JTAIL1, JTAIL2 DEFAULTS ARE ', > JTAIL1, JTAIL2 READ(IREAD,*) JTAIL1, JTAIL2 ELSE JTAIL1 = 1 JTAIL2 = JMAX ENDIF C *6* WRITE(IWRIT,*) 'ENTER GRID SPECIFICATIONS: ' WRITE(IWRIT,*) 'ENTER XMIN, XMAX, YMAX DEFAULTS ARE ', * XMIN,XMAX,YMAX READ(IREAD,*) XMIN, XMAX, YMAX WRITE(IWRIT,*)'ENTER DYM, THICK DEFAULTS ARE ',DYM, THICK READ(IREAD,*) DYM, THICK C *7* WRITE(IWRIT,*) 'ENTER PHIDT, THETAD DEFAULTS ARE ', * PHIDT, THETAD READ(IREAD,*) PHIDT, THETAD C *8* WRITE(IWRIT,*) 'ENTER JACDT DEFAULT ', JACDT READ(IREAD,*) JACDT C *9* WRITE(IWRIT,*) 'ENTER IPRINT DEFAULT ', IPRINT READ(IREAD,*) IPRINT C *10* WRITE(IWRIT,*) 'ENTER NP, NPCP DEFAULTS ARE ',NP,' ',NPCP READ(IREAD,*) NP, NPCP C *11* C C C CHECK MESH SIZES C IF ( JMAX .GT. JLARGE ) THEN WRITE(IWRIT, 36) JMAX, JLARGE 36 FORMAT( 1X, 'JMAX IS ', I5, 5X, 'ONLY DIMENSIONED FOR ', I5) STOP 'ERROR 1 IN IOALL' ENDIF IF ( JMAX*KMAX .GT. JKLARG ) THEN WRITE(IWRIT, 37) JMAX*KMAX, JKLARG 37 FORMAT( 1X,'JMAX*KMAX IS ',I5,5X,'ONLY DIMENSIONED FOR ',I5) STOP 'ERROR 1 IN IOALL' ENDIF IF ( JMAX .GT. 999 .OR. KMAX .GT. 999 ) THEN WRITE(IWRIT, 38) JMAX,KMAX 38 FORMAT( 1X, 'JMAX,KMAX ARE ', 2I5, 5X, * 'BIGGER THAN 999 ', I5) STOP 'ERROR 1 IN IOALL' ENDIF C WRITE(IWRIT,*) ' ' C WRITE(IWRIT,*)' ' WRITE(IWRIT,*)'...........FLOW FIELD CONDITIONS...............' WRITE(IWRIT,*)' FSMACH = ',FSMACH,' ALPHA = ',ALPHA C WRITE(IWRIT,*)' ' WRITE(IWRIT,*)'.............GRID SPECIFICATIONS...............' WRITE(IWRIT,*)' BICONVEX AIRFOIL GRID ', 1 'IS BEING INTERNALLY GENERATED' WRITE(IWRIT,*)'...............................................' WRITE(IWRIT,*)' MINIMUM SPACING AT BODY' WRITE(IWRIT,*)' DYM = ',DYM WRITE(IWRIT,*)' UPPER DIMENSION IN CHORDS (C=1) ' WRITE(IWRIT,*)' YMAX = ',YMAX WRITE(IWRIT,*)' LEFT XMIN = ',XMIN,' RIGHT XMAX = ',XMAX WRITE(IWRIT,*)' THICKNESS OF BICONVEX THICK = ',THICK WRITE(IWRIT,*)'...............................................' IF(PERIDC )THEN WRITE(IWRIT,*)' THIS IS A PERIDC GRID', 1 ' JMAX = ',JMAX,' KMAX = ',KMAX ELSE WRITE(IWRIT,*)' THIS IS A NONPERIDC GRID', 1 ' JMAX = ',JMAX,' KMAX = ',KMAX, 1 ' JTAIL1 = ',JTAIL1,' JTAIL2 = ',JTAIL2 ENDIF C C WRITE(IWRIT,*)' ' WRITE(IWRIT,*)'.............I/O SPECIFICATIONS.............' WRITE(IWRIT,*)' NP = ',NP,' IPRINT = ',IPRINT WRITE(IWRIT,*)'.............ALGORITHM PARAMETERS.............' WRITE(IWRIT,*)' THETAD = ',THETAD ,' PHIDT = ',PHIDT WRITE(IWRIT,*)' JACDT = ',JACDT C C IF(JACDT.EQ.0)WRITE(IWRIT,*) 1 ' <<<<<<<<<<<< TIME ACCURATE >>>>>>>>>>>>>>>>>>' IF(JACDT.EQ.1)WRITE(IWRIT,*) 1 ' <<<< VARIABLE DT = 1./(1+SQRT(J)) >>>>>>>>>>>>>' IF(JACDT.EQ.2)WRITE(IWRIT,*) 1 ' <<<< VARIABLE DT : MIN. EIGEN. >>>>>>>>>>>>>>>' IF(JACDT.EQ.3)WRITE(IWRIT,*) 1 ' <<<< VARIABLE DT : CONST. CFL >>>>>>>>>>>>>>>' WRITE(IWRIT,*)' ******** NUMBER OF TIME STEPS == ',NSTEPS WRITE(IWRIT,*)' ******** TIME STEP == ',DT C WRITE(IWRIT,*)'*************************************************' WRITE(IWRIT,*)'...... PENTADIAGONAL ALGORITHM ...........' WRITE(IWRIT,*)'*************************************************' WRITE(IWRIT,*)' DIS2X = ',DIS2X,' DIS4X = ',DIS4X WRITE(IWRIT,*)' DIS2Y = ',DIS2Y,' DIS4Y = ',DIS4Y C C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C********************************************************************* C********************** MAIN INTEGRATION ROUTINE ********************* C********************************************************************* C SUBROUTINE INTEGR(JDIM,KDIM,Q,S,PRESS,SNDSP,SLOPE,X,Y,XY,XYJ, > XIT,ETT,DS,DELTQN, > UU,VV,CCX,CCY,COEF2X,COEF2Y,COEF4X,COEF4Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C MAIN INTEGRATION CALLS C PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM*KDIM*4) DIMENSION PRESS(JDIM*KDIM),SNDSP(JDIM*KDIM) DIMENSION S(JDIM*KDIM*4),XY(JDIM*KDIM*4),XYJ(JDIM*KDIM) DIMENSION XIT(JDIM*KDIM),ETT(JDIM*KDIM),DS(JDIM*KDIM) DIMENSION X(JDIM*KDIM),Y(JDIM*KDIM),SLOPE(JDIM) C DIMENSION DELTQN(JDIM*KDIM*4) DIMENSION UU(JDIM*KDIM),CCX(JDIM*KDIM) DIMENSION VV(JDIM*KDIM),CCY(JDIM*KDIM) DIMENSION COEF2X(JDIM*KDIM),COEF2Y(JDIM*KDIM) DIMENSION COEF4X(JDIM*KDIM),COEF4Y(JDIM*KDIM) C COMMON/WORK/WORK1(MAXJK,10) C C C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C <<< EXPLICIT OPERATORS >>> C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C C********************************************************************* C*** 0. UPDATE PRESSURE AND SOUND SPEED *** C********************************************************************* C CALL CALCPS(JDIM,KDIM,Q,PRESS,SNDSP) C C********************************************************************* C*** 1. EXPLICIT BOUNDARY CONDITIONS *** C********************************************************************* C CALL BC(JDIM,KDIM,Q,PRESS,SNDSP,XY,XIT,ETT,XYJ,X,Y,SLOPE) C C********************************************************************* C*** 2. EXPLICIT RHS *** C********************************************************************* C--------------------------------------------------------------------- C--- EXPLICIT OPERATOR IN XI --- C--------------------------------------------------------------------- CALL RHSX(JDIM,KDIM,Q,S,PRESS,XY,XYJ,XIT,WORK1) C C--------------------------------------------------------------------- C EXPLICIT OPERATORS IN ETA C--------------------------------------------------------------------- CALL RHSY(JDIM,KDIM,Q,S,PRESS,XY,XYJ,ETT,WORK1) C C C--------------------------------------------------------------------- C--- 3-A-1. CALCULATE EIGENVALUES FOR XI AND ETA --- C--------------------------------------------------------------------- C IX = 1 IY = 2 CALL EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ, * XY,XIT,UU,CCX) IX = 3 IY = 4 CALL EIGVAL(IX,IY,JDIM,KDIM,Q,PRESS,SNDSP,XYJ, * XY,ETT,VV,CCY) C C C--------------------------------------------------------------------- C--- 3-A-2. FORM SPECTRAL RADII FOR DISSIPATION MODEL --- C--------------------------------------------------------------------- C CALL SPECT(JDIM,KDIM,UU,VV,CCX,CCY,XYJ,XY,COEF4X,COEF4Y) C C********************************************************************* C*** 3-B. DISSIPATION TERMS *** C********************************************************************* C C--------------------------------------------------------------------- C--- 3-B-1-A. FORM PRESSURE COEFFICIENTS --- C--------------------------------------------------------------------- IXY = 1 CALL GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF2X,WORK1) IXY = 2 CALL GRADCO(IXY,JDIM,KDIM,PRESS,XYJ,COEF2Y,WORK1) C C--------------------------------------------------------------------- C--- 3-B-1-B. ARTIFICIAL DISSIPATION COEFFICIENTS --- C--------------------------------------------------------------------- C CALL COEF24(JDIM,KDIM,COEF2X,COEF4X,COEF2Y,COEF4Y, > WORK1(1,1),WORK1(1,2)) C C--------------------------------------------------------------------- C--- 3-B-2. ARTIFICIAL DISSIPATION --- C--------------------------------------------------------------------- C CALL FILERX(JDIM,KDIM,Q,S,XYJ,COEF2X,COEF4X,WORK1) CALL FILERY(JDIM,KDIM,Q,S,XYJ,COEF2Y,COEF4Y,WORK1) C C--------------------------------------------------------------------- C C********************************************************************* C 4. COMPUTE RESIDUAL AND SCALE BY VARIABLE DT C********************************************************************* C--------------------------------------------------------------------- C--- 4-A. COMPUTE L2 AND MAX NORMS OF RESIDUAL --- C--------------------------------------------------------------------- C C DENSITY RESIDUAL C C C NRES = 1 CALL RESIL2(JDIM,KDIM,S,NRES) C********************************************************************* C********************************************************************* C C--------------------------------------------------------------------- C--- 4-B. SCALE WITH VARIABLE DT --- C--------------------------------------------------------------------- C CALL SCALDT(JDIM,KDIM,S,DS) C C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<< IMPLICIT OPERATORS >>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C--------------------------------------------------------------------- C 5. SCALAR PENTADIAGONAL ALGORITHM C--------------------------------------------------------------------- C C--------------------------------------------------------------------- C--- 5-A. INTEGRATION XI - ETA ORDER --- C--------------------------------------------------------------------- C C..................................................................... C... -1 ... C... 5-A-1. T S EIGENVECTOR INVERSE MULTIPLY BY RHS C... XI CALL WITH 1,2 FOR XI METRICS ... C..................................................................... C IX = 1 IY = 2 CALL TKINV(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY) C C..................................................................... C... 5-A-2. XI DIAGONAL INVERSIONS ... C..................................................................... C CALL STEPFX(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS, * COEF2X,COEF4X,UU,CCX,WORK1) C C..................................................................... C... N INVERSE MULTIPLY BETWEEN DIAGONAL OPERATORS ... C..................................................................... C IXX = 1 IXY = 2 IYX = 3 IYY = 4 CALL NINVER(JDIM,KDIM,S,XY,IXX,IXY,IYX,IYY) C..................................................................... C... ETA DIAGONAL INTEGRATION ... C..................................................................... C CALL STEPFY(JDIM,KDIM,Q,S,PRESS,SNDSP, * XY,XYJ,DS,COEF2Y,COEF4Y,VV,CCY,WORK1) C..................................................................... C... T S MATRIX MULTIPLY TO CLOSE DIAGONALIZATION ... C... ETA 3, 4 FOR ETA METRICS ... C..................................................................... C IX = 3 IY = 4 CALL TK(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY) C C..................................................................... C--------------------------------------------------------------------- C--- END 5-A. --- C--------------------------------------------------------------------- C********************************************************************* C*** END 5. *** C********************************************************************* C C********************************************************************* C C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<< FINISH INTEGRATION >>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C C********************************************************************* C*** 7. UPDATE *** C********************************************************************* C CALL UPDATE(JDIM,KDIM,Q,S,DELTQN) C C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** INPUT/OUTPUT ROUTINE ******************************** C*********************************************************************** SUBROUTINE IOALL(IWHAT,JDIM,KDIM,Q,PRESS,SNDSP,XY,XYJ,X,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4) DIMENSION PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM) C COMMON/WORKSP/CPBOD(MAXJ),XBOD(MAXJ), WORK2(MAXJ,90) C COMMON/IO/IREAD,IWRIT,IVER INTEGER J1 C GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17) IWHAT C 1 CONTINUE C GO TO 10000 C C C*********************************************************************** C 2 CONTINUE GO TO 10000 C C********************************************************************** C C 3 CONTINUE C C GO TO 10000 C C C*********************************************************************** C READ THE RESTART C 4 CONTINUE C C READ IN A UNFORMATTED RESTART FILE FROM TAPE3. C THE Q READ IN IS NOT YET SCALED BY THE METRIC JACOBIAN XYJ, C THE RESCALING IS DONE IN INITIA. C C CHANGED TO BE COMPATIBLE WITH NEW PLOT3D 4/2/86 C OPEN(UNIT=3, FILE='SRI3',STATUS='OLD', FORM='UNFORMATTED') READ(3) JMAX2,KMAX2 READ(3) FSMTEM,ALPTEM,RETEM,TOTIME C WRITE(IWRIT, 433) TOTIME, ISTART 433 FORMAT(1H0,12HTOTAL TIME =,F12.6, 1 10X, 20HSTARTING ITERATION =, I5) C READ(3) ( ( ( Q(J,K,N),J=1,JMAX2) ,K=1,KMAX2) ,N=1,4) READ(3) JTITEM, ISTART C GO TO 10000 C C C*********************************************************************** C ITERATION NUMBER AND RESIDUAL WRITTEN TO OUTPUT C 5 CONTINUE C C C COMPUTE NUMBER OF SUPERSONIC POINTS C NSUPER = 0 J1 = JUP - JLOW + 1 d call trace_enter (2097) CDOACROSS 501 K=1,KMAX DOUBLE PRECISION RI1(J1), UV21(J1), CSPED21(J1), ZS1(J1) RI1(1:$J1) = 1. / Q(JLOW:$J1,K,1) UV21(1:$J1) = (Q(JLOW:$J1,K,2) ** 2 + Q(JLOW:$J1,K,3) ** 2) * X RI1(1:$J1) ** 2 C CHECK MACH**2 ( SAME AS MACH) CSPED21(1:$J1) = SNDSP(JLOW:$J1,K) ** 2 ZS1(1:$J1) = UV21(1:$J1) / CSPED21(1:$J1) CALL AWAIT (1, 1, NSUPER) NSUPER = NSUPER + COUNT (MASK=ZS1(1:$J1) .GE. 1.) CALL ADVANCE (1, NSUPER) 501 CONTINUE C C LIFT AND DRAG CALCULATIONS AND OUTPUT ON TAPE6 C d call trace_exit (2098) RSD = RESID C C OUTPUT C C IF (NUMITE.EQ.100) THEN CALL VERIFY(NUMITE,RSD) ENDIF WRITE(IWRIT, 505)NUMITE , RSD, NSUPER, CPUSEC 505 FORMAT(1H ,' N = ',I6, ' RESID = ',E12.6, > ' # OF SS PTS = ',I5,' CPUSEC = ',F10.2) C GO TO 10000 C C C*********************************************************************** C COEFFICIENT OF PRESSURE DATA WRITTEN TO OUTPUT C 6 CONTINUE C C WRITE(IWRIT, 666) 666 FORMAT( 1X, 23HCOEFFICIENT OF PRESSURE ) CPC = 2./(GAMMA*FSMACH**2) C COMPUTE CP AT GRID POINTS AND STORE IN Z ARRAY d call trace_enter (2099) DO 665 J=JTAIL1,JTAIL2 PP = PRESS(J,1)*XYJ(J,1) CP = (PP*GAMMA -1.)*CPC WRITE(IWRIT, 667) J, X(J,1), CP 665 CONTINUE C d call trace_exit (2100) 667 FORMAT( 1X, I5, 5X, F14.7, 5X, F14.7) C C GO TO 10000 C C*********************************************************************** 7 CONTINUE GO TO 10000 C C*********************************************************************** C 8 CONTINUE C GO TO 10000 C C C*********************************************************************** C 9 CONTINUE C GO TO 10000 C C C*********************************************************************** C STORE A UNFORMATTED RESTART FILE C 10 CONTINUE C C NOTE THAT NUMITE IS ACTUALLY ONE GREATER THAN THE LAST USED. C IT GETS STORED, AND USED AS THE FIRST ITERATION NEXT TIME. C C REWIND BECAUSE THIS MAY BE CALLED PERIDC ALLY OPEN(UNIT=4, FILE='SRO4',STATUS='NEW', FORM='UNFORMATTED') WRITE(4) JMAX,KMAX WRITE(4) FSMACH,ALPHA,FSMACH,ALPHA WRITE(4) ( ( ( Q(J,K,N)*XYJ(J,K),J=1,JMAX) ,K=1,KMAX) ,N=1,4) WRITE(4) JTAIL1, NUMITE REWIND 4 C GO TO 10000 C C*********************************************************************** C 11 CONTINUE GOTO 10000 C*********************************************************************** C 12 CONTINUE GOTO 10000 C*********************************************************************** 13 CONTINUE GOTO 10000 C*********************************************************************** C*********************************************************************** 14 CONTINUE GOTO 10000 C*********************************************************************** C*********************************************************************** C 15 CONTINUE C C C OUTPUT ON TAPE6 9 VARIABLES ALONG THE BODY C K = 1 WRITE(IWRIT, 1510) K 1510 FORMAT( 1X, 4HK = , I4 ) WRITE(IWRIT, 1515) 1515 FORMAT( 1X, 1 51H J K RHO U V , 2 25HENERGY PRESSURE , 3 42H COEFF P LOCAL MACH X Y ) C d call trace_enter (2101) DO 1590 J = 1, JMAX RHO = Q( J, K, 1) * XYJ( J, K) U = Q( J, K, 2) / Q( J, K, 1) V = Q( J, K, 3) / Q( J, K, 1) UVSQ = U**2 + V**2 ENERGY = Q( J, K, 4) * XYJ( J, K) PPP = PRESS(J,K)*XYJ(J,K) CPC = 2./(GAMMA*FSMACH**2) CP = (PPP * GAMMA - 1. ) * CPC CSQ = SNDSP(J,K) EMACH = SQRT( UVSQ )/ CSQ XCOORD = X( J, K) YCOORD = Y( J, K) WRITE(IWRIT, 1520) J, K, RHO, U, V, ENERGY, PPP, 1 CP, EMACH, XCOORD, YCOORD 1520 FORMAT( 1X, 2( I3, 1X), 9( F11.6, 2X ) ) 1590 CONTINUE C d call trace_exit (2102) GOTO 10000 C C*********************************************************************** C*********************************************************************** C GRID OUTPUT C 16 CONTINUE C OPEN(UNIT=11, FILE='SRO11',STATUS='NEW', FORM='UNFORMATTED') WRITE( 11) JMAX,KMAX WRITE( 11) ( ( X(J,K),J=1,JMAX) ,K=1,KMAX) , & ( ( Y(J,K),J=1,JMAX) ,K=1,KMAX) C C GO TO 10000 C C C*********************************************************************** C*********************************************************************** C COMPUTE CIRCULATION, AND PRINT IT OUT ( NOT IN COMMON ) C 17 CONTINUE C C GO TO 10000 C C C*********************************************************************** C C 10000 RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** N INVERSE MULTIPLY ON S **************************** C*********************************************************************** SUBROUTINE NINVER(JDIM,KDIM,S,XY,IXX,IXY,IYX,IYY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C IXX = 1, IXY = 2, IYX = 3, IYY = 4 C -1 C N = T T C XI ETA C C IXX = 3, IXY = 4, IYX = 1, IYY = 2 C -1 C N = T T C ETA XI C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4) C DIMENSION RG(4,4),SR(4) INTEGER K1, J1 C C PARAMETERS FOR DIAGONALIZATION C BT = 1./SQRT(2.) BT2 = 0.5 K1 = KUP - KLOW + 1 J1 = JUP - JLOW + 1 C d call trace_enter (2103) CDOALL 1 K=KLOW,KLOW+K1-1 DOUBLE PRECISION RR11(J1), RR21(J1), RM11(J1), RM21(J1), S34M1( X J1), S34P1(J1), S21(J1), S31(J1), S41(J1) RR11(1:$J1) = XY(JLOW:$J1,K,IXX) ** 2 + XY(JLOW:$J1,K,IXY) ** X 2 RR21(1:$J1) = XY(JLOW:$J1,K,IYX) ** 2 + XY(JLOW:$J1,K,IYY) ** X 2 RR11(1:$J1) = 1. / SQRT (RR11(1:$J1)*RR21(1:$J1)) c RR21(1:$J1) = 1. / SQRT (RR21(1:$J1)) c squareroots rr11 and rr21 combined RM11(1:$J1) = (XY(JLOW:$J1,K,IXX) * XY(JLOW:$J1,K,IYX) + XY( X JLOW:$J1,K,IXY) * XY(JLOW:$J1,K,IYY)) * RR11(1:$J1) RM21(1:$J1) = (XY(JLOW:$J1,K,IXX) * XY(JLOW:$J1,K,IYY) - XY( X JLOW:$J1,K,IXY) * XY(JLOW:$J1,K,IYX)) * RR11(1:$J1) S34M1(1:$J1) = S(JLOW:$J1,K,3) - S(JLOW:$J1,K,4) S34P1(1:$J1) = S(JLOW:$J1,K,3) + S(JLOW:$J1,K,4) S21(1:$J1) = RM11(1:$J1) * S(JLOW:$J1,K,2) + BT * S34M1(1 X :$J1) * RM21(1:$J1) S31(1:$J1) = -BT * S(JLOW:$J1,K,2) * RM21(1:$J1) + BT2 * X (S34P1(1:$J1) + RM11(1:$J1) * S34M1(1:$J1)) S41(1:$J1) = BT * RM21(1:$J1) * S(JLOW:$J1,K,2) + BT2 * ( X S34P1(1:$J1) - RM11(1:$J1) * S34M1(1:$J1)) S(JLOW:$J1,K,2) = S21(1:$J1) S(JLOW:$J1,K,3) = S31(1:$J1) S(JLOW:$J1,K,4) = S41(1:$J1) 1 CONTINUE C d call trace_exit (2104) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************ COMPUTE L 2 NORM OF S , THE RESIDUAL ****************** C*********************************************************************** SUBROUTINE RESIL2(JDIM,KDIM,S,NRES) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION S(JDIM,KDIM,4) INTEGER J1 C C COMPUTE L2 NORM OF DENSITY RESIDUAL ( RIGHT HAND SIDE ) C RESID = 0. J1 = JUP - JLOW + 1 d call trace_enter (2105) DO 5 K=KLOW,KUP RESID = RESID + SUM (S(JLOW:$J1,K,NRES) ** 2) 5 CONTINUE d call trace_exit (2106) RESID = RESID/( (JUP-JLOW+1)*(KUP-KLOW+1)) RESID = SQRT( RESID ) / DT 6 CONTINUE RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************* EVALUATION OF RIGHT HAND SIDE ( EXPLICIT PART ) ********* C*********************************************************************** SUBROUTINE RHSX(JDIM,KDIM,Q,S,PRESS,XY,XYJ,XIT,FLUX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),S(JDIM,KDIM,4) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM),XIT(JDIM,KDIM) DIMENSION FLUX(JDIM,KDIM,4) INTEGER K1, J1 K1 = KUP - KLOW + 1 J1 = JEND - JBEGIN + 1 C C CENTRAL DIFFERENCING USED IN XI C COMPUTE FLUX VECTORS C d call trace_enter (2107) CDOALL 200 K=KLOW,KLOW+K1-1 DOUBLE PRECISION R11(J1), R21(J1), R01(J1), RR1(J1), U1(J1), V1(J1 X ), QS1(J1), PRESSJ1(J1) R11(1:$J1) = XY(JBEGIN:$J1,K,1) R21(1:$J1) = XY(JBEGIN:$J1,K,2) R01(1:$J1) = XIT(JBEGIN:$J1,K) C RR IS J/RHO RR1(1:$J1) = 1. / Q(JBEGIN:$J1,K,1) U1(1:$J1) = Q(JBEGIN:$J1,K,2) * RR1(1:$J1) V1(1:$J1) = Q(JBEGIN:$J1,K,3) * RR1(1:$J1) C C R1, R2 ARE EITHER D XI /DX, D ETA/DX, OR D XI /DY, D ETA/DY C QS IS EITHER CAP-U OR CAP-V ( CONTRAVARIANT VELOCITIES ) C QS1(1:$J1) = R01(1:$J1) + R11(1:$J1) * U1(1:$J1) + R21(1:$J1) * X V1(1:$J1) C C PRESSJ IS PRESSURE / J C PRESSJ1(1:$J1) = PRESS(JBEGIN:$J1,K) C FLUX(JBEGIN:$J1,K,1) = Q(JBEGIN:$J1,K,1) * QS1(1:$J1) FLUX(JBEGIN:$J1,K,2) = Q(JBEGIN:$J1,K,2) * QS1(1:$J1) + R11(1:$ X J1) * PRESSJ1(1:$J1) FLUX(JBEGIN:$J1,K,3) = Q(JBEGIN:$J1,K,3) * QS1(1:$J1) + R21(1:$ X J1) * PRESSJ1(1:$J1) FLUX(JBEGIN:$J1,K,4) = QS1(1:$J1) * (Q(JBEGIN:$J1,K,4) + X PRESSJ1(1:$J1)) - PRESSJ1(1:$J1) * R01(1:$J1) C 200 CONTINUE C C d call trace_exit (2108) RT0 = - .5*DT / (1. + PHIDT) J1 = JUP - JLOW + 1 K1 = KUP - KLOW + 1 d call trace_enter (2109) CDOALL 400 J=JLOW,JLOW+J1-1 INTEGER N, JM1, JP1 DO 400 N=1,4 S(J,KLOW:$K1,N) = S(J,KLOW:$K1,N) + RT0 * (FLUX(JPLUS(J), X KLOW:$K1,N) - FLUX(JMINU(J),KLOW:$K1,N)) 400 CONTINUE C C d call trace_exit (2110) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************* EVALUATION OF RIGHT HAND SIDE ( EXPLICIT PART ) ********* C*********************************************************************** SUBROUTINE RHSY(JDIM,KDIM,Q,S,PRESS,XY,XYJ,ETT,FLUX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),S(JDIM,KDIM,4) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM),ETT(JDIM,KDIM) DIMENSION FLUX(JDIM,KDIM,4) INTEGER K1, J1 K1 = KEND - KBEGIN + 1 J1 = JUP - JLOW + 1 C C CENTRAL DIFFERENCING USED IN ETA DIRECTION d call trace_enter (2111) CDOALL 20 K=KBEGIN,KBEGIN+K1-1 DOUBLE PRECISION R11(J1), R21(J1), R01(J1), RR1(J1), U1(J1), V1(J1 X ), QS1(J1), PRESSJ1(J1) R11(1:$J1) = XY(JLOW:$J1,K,3) R21(1:$J1) = XY(JLOW:$J1,K,4) R01(1:$J1) = ETT(JLOW:$J1,K) C RR IS J/RHO RR1(1:$J1) = 1. / Q(JLOW:$J1,K,1) U1(1:$J1) = Q(JLOW:$J1,K,2) * RR1(1:$J1) V1(1:$J1) = Q(JLOW:$J1,K,3) * RR1(1:$J1) C C R1, R2 ARE EITHER D XI /DX, D ETA/DX, OR D XI /DY, D ETA/DY C QS IS EITHER CAP-U OR CAP-V ( CONTRAVARIANT VELOCITIES ) QS1(1:$J1) = R01(1:$J1) + R11(1:$J1) * U1(1:$J1) + R21(1:$J1) * X V1(1:$J1) C C PRESSJ IS PRESSURE / J PRESSJ1(1:$J1) = PRESS(JLOW:$J1,K) C FLUX(JLOW:$J1,K,1) = Q(JLOW:$J1,K,1) * QS1(1:$J1) FLUX(JLOW:$J1,K,2) = Q(JLOW:$J1,K,2) * QS1(1:$J1) + R11(1:$J1) X * PRESSJ1(1:$J1) FLUX(JLOW:$J1,K,3) = Q(JLOW:$J1,K,3) * QS1(1:$J1) + R21(1:$J1) X * PRESSJ1(1:$J1) FLUX(JLOW:$J1,K,4) = QS1(1:$J1) * (Q(JLOW:$J1,K,4) + PRESSJ1(1 X :$J1)) - PRESSJ1(1:$J1) * R01(1:$J1) C 20 CONTINUE C d call trace_exit (2112) R0 = - .5*DT / (1. + PHIDT) J1 = JUP - JLOW + 1 C CENTRAL DIFFERENCE F d call trace_enter (2113) CDOALL 30 N=1,4 INTEGER K DO 30 K=KLOW,KUP S(JLOW:$J1,K,N) = R0 * (FLUX(JLOW:$J1,K+1,N) - FLUX(JLOW:$J1, X K-1,N)) + S(JLOW:$J1,K,N) 30 CONTINUE C d call trace_exit (2114) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C************* SCALE RHS WITH VARIABLE DT ************************* C*********************************************************************** SUBROUTINE SCALDT(JDIM,KDIM,S,DS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION S(JDIM,KDIM,4),DS(JDIM,KDIM) INTEGER J1 J1 = JEND - JBEGIN + 1 C C C SCALE WITH SPATIALLY VARIABLE TIME STEP C d call trace_enter (2115) CDOALL 10 M=1,4 INTEGER K DO 10 K=KBEGIN,KEND S(JBEGIN:$J1,K,M) = S(JBEGIN:$J1,K,M) * DS(JBEGIN:$J1,K) 10 CONTINUE C d call trace_exit (2116) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C C*********************************************************************** C************** COMPUTE EIGENVALUES *********************************** C*********************************************************************** SUBROUTINE SPECT(JDIM,KDIM,UU,VV,CCX,CCY,XYJ,XY,SPECTX,SPECTY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C C FORM SPECTX AND SPECTY INDEPENDENTLY C C SPECTX = (ABS(UU) + CCX)/XYJ C SPECTY = (ABS(VV) + CCY)/XYJ C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION XYJ(JDIM,KDIM),XY(JDIM,KDIM,4) DIMENSION SPECTX(JDIM,KDIM),UU(JDIM,KDIM),CCX(JDIM,KDIM) DIMENSION SPECTY(JDIM,KDIM),VV(JDIM,KDIM),CCY(JDIM,KDIM) INTEGER K1, J1 K1 = KEND - KBEGIN + 1 J1 = JEND - JBEGIN + 1 C C COMPUTE DISSIPATION SCALING C C C SPECTRAL RADIUS C d call trace_enter (2117) CDOALL 100 K=KBEGIN,KBEGIN+K1-1 DOUBLE PRECISION SIGX1(J1), SIGY1(J1), RJ1(J1) C SIGX1(1:$J1) = ABS (UU(JBEGIN:$J1,K)) + CCX(JBEGIN:$J1,K) SIGY1(1:$J1) = ABS (VV(JBEGIN:$J1,K)) + CCY(JBEGIN:$J1,K) C RJ1(1:$J1) = 1. / XYJ(JBEGIN:$J1,K) SPECTX(JBEGIN:$J1,K) = SIGX1(1:$J1) * RJ1(1:$J1) SPECTY(JBEGIN:$J1,K) = SIGY1(1:$J1) * RJ1(1:$J1) C 100 CONTINUE C d call trace_exit (2118) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** MAIN INTEGRATION ROUTINE **************************** C********************* DIAGONAL VERSION ******************************** C******** WITH IMPLICIT FOURTH ORDER DISSIPATION ******************** C ROUTINE VECTORIZED THP 10/18/83 C*********************************************************************** SUBROUTINE STEPFX(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS, * COEF2,COEF4,UU,CC,WORK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION DS(JDIM,KDIM) C DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM) DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM) C DIMENSION WORK(JDIM,KDIM,10) INTEGER J1, K1, K2, J2, K3, J3, J4, J5, J6, J7, K4, K5, K6, K9, X J10 C C NOTE : D USED FOR VECTORIZATION C DTD = DT * THETAD /(1. + PHIDT) HD = 0.5*DTD C C XI DIRECTION C C FILL MATRIX FOR INVERSION FOR DIAGONAL IN XI C C SET SWITCHES FOR EIGENVALUES C C DO FIRST TWO AT THE SAME TIME SINCE THEY HAVE THE SAME COEFFICIENTS C d call trace_enter (2119) DO 300 N = 2,4 IF ( N .EQ. 2) SN = 0. IF ( N .EQ. 3) SN = 1. IF ( N .EQ. 4) SN = -1. K3 = KUP - KLOW + 1 C C J3 = (JUP - JLOW + 32) / 32 J4 = JUP - JLOW + 1 C CDOALL 210 J1=0,32*J3-32,32 C INTEGER J8, J9, JM2, JP2, JM1, JP1 C J9 = JLOW + J1 C J8 = MIN0 (JLOW + J1 + 32 - 1, JUP) - (JLOW + J1) + 1 C210 CONTINUE CDOALL 212 K=KLOW,KLOW+K3-1 DOUBLE PRECISION C2M1(J4), C4M1(J4), C21(J4), C41(J4) C2M1(1:$J4) = COEF2(JMINU(JLOW:$J4),K) * DTD C4M1(1:$J4) = COEF4(JMINU(JLOW:$J4),K) * DTD C21(1:$J4) = COEF2(JLOW:$J4,K) * DTD C41(1:$J4) = COEF4(JLOW:$J4,K) * DTD WORK(JLOW:$J4,K,1) = XYJ(JMINU(JMINU(JLOW:$J4))+((-1)+K)* X JDIM,1) * C4M1(1:$J4) WORK(JLOW:$J4,K,2) = -(C2M1(1:$J4) + 3. * C4M1(1:$J4) + C41( X 1:$J4)) * XYJ(JMINU(JLOW:$J4),K) WORK(JLOW:$J4,K,3) = XYJ(JLOW:$J4,K) * (C2M1(1:$J4) + 3. * X C4M1(1:$J4) + C21(1:$J4) + 3. * C41(1:$J4)) WORK(JLOW:$J4,K,4) = -(C21(1:$J4) + 3. * C41(1:$J4) + C4M1(1 X :$J4)) * XYJ(JPLUS(JLOW:$J4),K) WORK(JLOW:$J4,K,5) = XYJ(JPLUS(JPLUS(JLOW:$J4))+((-1)+K)* X JDIM,1) * C41(1:$J4) 212 CONTINUE C IF(.NOT.PERIDC )THEN J = JLOW JP1 = JPLUS(J) JM1 = JMINU (J) JP2 = JPLUS(JP1) K4 = (KUP - KLOW + 32) / 32 K6 = MAX0 (MIN0 (KLOW + 31, KUP) - KLOW + 1, 1) CDOALL 220 K1=0,32*K4-32,32 DOUBLE PRECISION C2M2(K6), C4M2(K6), C22(K6), C42(K6) INTEGER K7, K8 K8 = KLOW + K1 K7 = MIN0 (KLOW + K1 + 32 - 1, KUP) - (KLOW + K1) + 1 C2M2(1:$K7) = COEF2(JM1,K8:$K7) * DTD C4M2(1:$K7) = COEF4(JM1,K8:$K7) * DTD C22(1:$K7) = COEF2(J,K8:$K7) * DTD C42(1:$K7) = COEF4(J,K8:$K7) * DTD WORK(J,K8:$K7,1) = 0. WORK(J,K8:$K7,2) = -(C2M2(1:$K7) + C4M2(1:$K7) + C42(1:$K7)) X * XYJ(JM1,K8:$K7) WORK(J,K8:$K7,3) = XYJ(J,K8:$K7) * (C2M2(1:$K7) + 2. * C4M2( X 1:$K7) + C22(1:$K7) + 3. * C42(1:$K7)) WORK(J,K8:$K7,4) = -(C22(1:$K7) + 3. * C42(1:$K7) + C4M2(1:$ X K7)) * XYJ(JP1,K8:$K7) WORK(J,K8:$K7,5) = XYJ(JP2,K8:$K7) * C42(1:$K7) 220 CONTINUE C J = JUP JP1 = JPLUS(J) JM1 = JMINU (J) JM2 = JMINU (JM1) K9 = (KUP - KLOW + 32) / 32 K6 = MAX0 (MIN0 (KLOW + 31, KUP) - KLOW + 1, 1) CDOALL 222 K2=0,32*K9-32,32 DOUBLE PRECISION C2M3(K6), C4M3(K6), C23(K6), C43(K6) INTEGER K10, K11 K11 = KLOW + K2 K10 = MIN0 (KLOW + K2 + 32 - 1, KUP) - (KLOW + K2) + 1 C2M3(1:$K10) = COEF2(JM1,K11:$K10) * DTD C4M3(1:$K10) = COEF4(JM1,K11:$K10) * DTD C23(1:$K10) = COEF2(J,K11:$K10) * DTD C43(1:$K10) = COEF4(J,K11:$K10) * DTD WORK(J,K11:$K10,1) = XYJ(JM2,K11:$K10) * C4M3(1:$K10) WORK(J,K11:$K10,2) = -(C2M3(1:$K10) + 3. * C4M3(1:$K10) + X C43(1:$K10)) * XYJ(JM1,K11:$K10) WORK(J,K11:$K10,3) = XYJ(J,K11:$K10) * (C2M3(1:$K10) + 3. * X C4M3(1:$K10) + C23(1:$K10) + 2. * C43(1:$K10)) WORK(J,K11:$K10,4) = -(C23(1:$K10) + C43(1:$K10) + C4M3(1:$ X K10)) * XYJ(JP1,K11:$K10) WORK(J,K11:$K10,5) = 0. 222 CONTINUE C ENDIF K3 = KUP - KLOW + 1 C C ADD IN FLUX JACOBIANS AND VARIABLE DT AND IDENTITY C C J10 = (JUP - JLOW + 32) / 32 J4 = JUP - JLOW + 1 C CDOALL 230 J2=0,32*J10-32,32 C INTEGER J11, J12, JM1, JP1 C J12 = JLOW + J2 C J11 = MIN0 (JLOW + J2 + 32 - 1, JUP) - (JLOW + J2) + 1 C230 CONTINUE CDOALL 232 K=KLOW,KLOW+K3-1 WORK(JLOW:$J4,K,2) = WORK(JLOW:$J4,K,2) - (UU(JMINU(JLOW:$J4 X ),K) + SN * CC(JMINU(JLOW:$J4),K)) * HD WORK(JLOW:$J4,K,4) = WORK(JLOW:$J4,K,4) + (UU(JPLUS(JLOW:$J4 X ),K) + SN * CC(JPLUS(JLOW:$J4),K)) * HD WORK(JLOW:$J4,K,1) = WORK(JLOW:$J4,K,1) * DS(JLOW:$J4,K) WORK(JLOW:$J4,K,2) = WORK(JLOW:$J4,K,2) * DS(JLOW:$J4,K) WORK(JLOW:$J4,K,3) = WORK(JLOW:$J4,K,3) * DS(JLOW:$J4,K) WORK(JLOW:$J4,K,4) = WORK(JLOW:$J4,K,4) * DS(JLOW:$J4,K) WORK(JLOW:$J4,K,5) = WORK(JLOW:$J4,K,5) * DS(JLOW:$J4,K) WORK(JLOW:$J4,K,3) = WORK(JLOW:$J4,K,3) + 1. 232 CONTINUE C C FOR N = 1 AND 2 DO BOTH TOGETHER C INVERT ON FIRST TWO S ELEMENTS C IF(N.EQ.2)THEN C IF(.NOT.PERIDC )THEN CALL XPENT2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * S(1,1,1),JLOW,JUP,KLOW,KUP) C ELSE C PERIDC CALL XPENP2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * WORK(1,1,9),WORK(1,1,10),S(1,1,1),JLOW,JUP,KLOW,KUP) ENDIF C ELSE C C N = 3 OR 4 C IF(.NOT.PERIDC )THEN CALL XPENTA(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * S(1,1,N),JLOW,JUP,KLOW,KUP) ELSE C PERIDC CALL XPENTP(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * WORK(1,1,9),WORK(1,1,10),S(1,1,N),JLOW,JUP,KLOW,KUP) ENDIF C ENDIF C 300 CONTINUE C C END XI DIRECTION C d call trace_exit (2120) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** MAIN INTEGRATION ROUTINE **************************** C********************* DIAGONAL VERSION ******************************** C******** WITH IMPLICIT FOURTH ORDER DISSIPATION ******************** C ROUTINE VECTORIZED THP 10/18/83 C*********************************************************************** SUBROUTINE STEPFY(JDIM,KDIM,Q,S,PRESS,SNDSP,XY,XYJ,DS, * COEF2,COEF4,UU,CC,WORK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),PRESS(JDIM,KDIM),SNDSP(JDIM,KDIM) DIMENSION S(JDIM,KDIM,4),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) DIMENSION DS(JDIM,KDIM) C DIMENSION COEF2(JDIM,KDIM),COEF4(JDIM,KDIM) DIMENSION UU(JDIM,KDIM),CC(JDIM,KDIM) C DIMENSION WORK(JDIM,KDIM,10) INTEGER J1, J2, J3, J4, J5, J6, J9, K1, K2 C C DTD = DT * THETAD /(1. + PHIDT) HD = 0.5*DTD C C ETA INVERSION PART OF DIAGONAL SCHEME C C SET EIGENVALUE SWTICHES C C DO N = 1 AND 2 TOGETHER C d call trace_enter (2121) DO 420 N = 2,4 IF ( N .EQ. 2) SN = 0. IF ( N .EQ. 3) SN = 1. IF ( N .EQ. 4) SN = -1. C C C FILL TRIDIAGONALS C K = KLOW KM1 = K-1 J3 = (JUP - JLOW + 32) / 32 J6 = MAX0 (MIN0 (JLOW + 31, JUP) - JLOW + 1, 1) CDOALL 422 J1=0,32*J3-32,32 DOUBLE PRECISION C2M1(J6), C21(J6), C41(J6) INTEGER J7, J8 J8 = JLOW + J1 J7 = MIN0 (JLOW + J1 + 32 - 1, JUP) - (JLOW + J1) + 1 C2M1(1:$J7) = COEF2(J8:$J7,KM1) * DTD C21(1:$J7) = COEF2(J8:$J7,K) * DTD C41(1:$J7) = COEF4(J8:$J7,K) * DTD WORK(J8:$J7,K,1) = 0. WORK(J8:$J7,K,2) = -(C2M1(1:$J7) + C41(1:$J7)) * XYJ(J8:$J7, X K-1) WORK(J8:$J7,K,3) = XYJ(J8:$J7,K) * (C2M1(1:$J7) + C21(1:$J7) X + 3. * C41(1:$J7)) WORK(J8:$J7,K,4) = -(C21(1:$J7) + 3. * C41(1:$J7)) * XYJ(J8 X :$J7,K+1) WORK(J8:$J7,K,5) = XYJ(J8:$J7,K+2) * C41(1:$J7) 422 CONTINUE C K = KUP KM1 = K-1 J9 = (JUP - JLOW + 32) / 32 J6 = MAX0 (MIN0 (JLOW + 31, JUP) - JLOW + 1, 1) CDOALL 426 J2=0,32*J9-32,32 DOUBLE PRECISION C2M2(J6), C4M1(J6), C22(J6), C42(J6) INTEGER J10, J11 J11 = JLOW + J2 J10 = MIN0 (JLOW + J2 + 32 - 1, JUP) - (JLOW + J2) + 1 C2M2(1:$J10) = COEF2(J11:$J10,KM1) * DTD C4M1(1:$J10) = COEF4(J11:$J10,KM1) * DTD C22(1:$J10) = COEF2(J11:$J10,K) * DTD C42(1:$J10) = COEF4(J11:$J10,K) * DTD WORK(J11:$J10,K,1) = XYJ(J11:$J10,K-2) * C4M1(1:$J10) WORK(J11:$J10,K,2) = -(C2M2(1:$J10) + 3. * C4M1(1:$J10) + X C42(1:$J10)) * XYJ(J11:$J10,K-1) WORK(J11:$J10,K,3) = XYJ(J11:$J10,K) * (C2M2(1:$J10) + 3. * X C4M1(1:$J10) + C22(1:$J10) + 2. * C42(1:$J10)) WORK(J11:$J10,K,4) = -(C22(1:$J10) + C42(1:$J10) + C4M1(1:$ X J10)) * XYJ(J11:$J10,K+1) WORK(J11:$J10,K,5) = 0. 426 CONTINUE K1 = KLOW + 1 K2 = KUP - KLOW - 1 J4 = JUP - JLOW + 1 C C FOURTH ORDER IN INTERIOR C CDOALL 428 K=K1,K1+K2-1 DOUBLE PRECISION C2M3(J4), C4M2(J4), C23(J4), C43(J4) INTEGER KM1 C2M3(1:$J4) = COEF2(JLOW:$J4,K-1) * DTD C4M2(1:$J4) = COEF4(JLOW:$J4,K-1) * DTD C23(1:$J4) = COEF2(JLOW:$J4,K) * DTD C43(1:$J4) = COEF4(JLOW:$J4,K) * DTD WORK(JLOW:$J4,K,1) = XYJ(JLOW:$J4,K-2) * C4M2(1:$J4) WORK(JLOW:$J4,K,2) = -(C2M3(1:$J4) + 3. * C4M2(1:$J4) + C43( X 1:$J4)) * XYJ(JLOW:$J4,K-1) WORK(JLOW:$J4,K,3) = XYJ(JLOW:$J4,K) * (C2M3(1:$J4) + 3. * X C4M2(1:$J4) + C23(1:$J4) + 3. * C43(1:$J4)) WORK(JLOW:$J4,K,4) = -(C23(1:$J4) + 3. * C43(1:$J4) + C4M2(1 X :$J4)) * XYJ(JLOW:$J4,K+1) WORK(JLOW:$J4,K,5) = XYJ(JLOW:$J4,K+2) * C43(1:$J4) 428 CONTINUE K1 = KUP - KLOW + 1 J5 = JUP - JLOW + 1 C C C ADD IN FLUX JACOBIANS, VARIABLE DT AND IDENTITY C CDOALL 435 K=KLOW,KLOW+K1-1 WORK(JLOW:$J5,K,2) = WORK(JLOW:$J5,K,2) - (UU(JLOW:$J5,K-1) + X SN * CC(JLOW:$J5,K-1)) * HD WORK(JLOW:$J5,K,4) = WORK(JLOW:$J5,K,4) + (UU(JLOW:$J5,K+1) + X SN * CC(JLOW:$J5,K+1)) * HD WORK(JLOW:$J5,K,1) = WORK(JLOW:$J5,K,1) * DS(JLOW:$J5,K) WORK(JLOW:$J5,K,2) = WORK(JLOW:$J5,K,2) * DS(JLOW:$J5,K) WORK(JLOW:$J5,K,3) = WORK(JLOW:$J5,K,3) * DS(JLOW:$J5,K) WORK(JLOW:$J5,K,4) = WORK(JLOW:$J5,K,4) * DS(JLOW:$J5,K) WORK(JLOW:$J5,K,5) = WORK(JLOW:$J5,K,5) * DS(JLOW:$J5,K) WORK(JLOW:$J5,K,3) = WORK(JLOW:$J5,K,3) + 1.0 435 CONTINUE C C COMBINE N = 1 AND 2 C INVERT ON FIRST TWO S ELEMENTS C IF(N.EQ.2)THEN C CALL YPENT2(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * S(1,1,1),KLOW,KUP,JLOW,JUP) C ELSE C DO N = 3 AND 4 C CALL YPENTA(JDIM,KDIM,WORK(1,1,1),WORK(1,1,2),WORK(1,1,3), * WORK(1,1,4),WORK(1,1,5),WORK(1,1,7),WORK(1,1,8), * S(1,1,N),KLOW,KUP,JLOW,JUP) C ENDIF C 420 CONTINUE C d call trace_exit (2122) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** TK MULTIPLY ON S ************************************ C*********************************************************************** SUBROUTINE TK(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C IX = 1, IY = 2 XI METRICS C IX = 3, IY = 4 ETA METRICS C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),SNDSP(JDIM,KDIM) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) INTEGER K1, J1 C BT = 1./SQRT(2.) RGAMI = 1./GAMI K1 = KUP - KLOW + 1 J1 = JUP - JLOW + 1 C C T MATRIX MULTIPLY FOR DIAGONAL ALGORITHM C K C d call trace_enter (2123) CDOALL 1 K=KLOW,KLOW+K1-1 DOUBLE PRECISION RX1(J1), RY1(J1), RRXY1(J1), RHO1(J1), RHOINV1(J1 X ), U1(J1), V1(J1), ALP01(J1), S11(J1), S21(J1), SA1(J1), SM1(J1) X , UV21(J1), RRX1(J1), RRY1(J1) RX1(1:$J1) = XY(JLOW:$J1,K,IX) RY1(1:$J1) = XY(JLOW:$J1,K,IY) RRXY1(1:$J1) = 1. / SQRT (RX1(1:$J1) ** 2 + RY1(1:$J1) ** 2) RX1(1:$J1) = RX1(1:$J1) * RRXY1(1:$J1) RY1(1:$J1) = RY1(1:$J1) * RRXY1(1:$J1) RHO1(1:$J1) = Q(JLOW:$J1,K,1) * XYJ(JLOW:$J1,K) RHOINV1(1:$J1) = 1. / Q(JLOW:$J1,K,1) U1(1:$J1) = Q(JLOW:$J1,K,2) * RHOINV1(1:$J1) V1(1:$J1) = Q(JLOW:$J1,K,3) * RHOINV1(1:$J1) ALP01(1:$J1) = RHO1(1:$J1) * BT / SNDSP(JLOW:$J1,K) S11(1:$J1) = S(JLOW:$J1,K,1) S21(1:$J1) = S(JLOW:$J1,K,2) SA1(1:$J1) = (S(JLOW:$J1,K,3) + S(JLOW:$J1,K,4)) * X ALP01(1:$J1) SM1(1:$J1) = (S(JLOW:$J1,K,3) - S(JLOW:$J1,K,4)) * RHO1 X (1:$J1) * BT UV21(1:$J1) = 0.5 * (U1(1:$J1) ** 2 + V1(1:$J1) ** 2) RRX1(1:$J1) = RHO1(1:$J1) * RX1(1:$J1) RRY1(1:$J1) = RHO1(1:$J1) * RY1(1:$J1) C S(JLOW:$J1,K,1) = S11(1:$J1) + SA1(1:$J1) S(JLOW:$J1,K,2) = U1(1:$J1) * S11(1:$J1) + RRY1(1:$J1) * S21(1:$J1 X ) + U1(1:$J1) * SA1(1:$J1) + RX1(1:$J1) * SM1(1:$J1) S(JLOW:$J1,K,3) = V1(1:$J1) * S11(1:$J1) - RRX1(1:$J1) * S21(1:$J1 X ) + V1(1:$J1) * SA1(1:$J1) + RY1(1:$J1) * SM1(1:$J1) S(JLOW:$J1,K,4) = UV21(1:$J1) * S11(1:$J1) + (RRY1(1:$J1) * U1(1:$ X J1) - RRX1(1:$J1) * V1(1:$J1)) * S21(1:$J1) + (UV21(1:$J1) + X SNDSP(JLOW:$J1,K) ** 2 * RGAMI) * SA1(1:$J1) + (RX1(1:$J1) * U1( X 1:$J1) + RY1(1:$J1) * V1(1:$J1)) * SM1(1:$J1) C 1 CONTINUE C d call trace_exit (2124) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** TK INVERSE MULTIPLY ON S **************************** C*********************************************************************** SUBROUTINE TKINV(JDIM,KDIM,Q,SNDSP,S,XYJ,XY,IX,IY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C IX = 1, IY = 2 XI METRICS C IX = 3, IY = 4 ETA METRICS C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),SNDSP(JDIM,KDIM) DIMENSION XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) INTEGER K1, J1 C BT = 1./SQRT(2.) K1 = KUP - KLOW + 1 J1 = JUP - JLOW + 1 C C T INVERSE MATRIX MULTIPLY FOR DIAGONAL ALGORITHM C K C d call trace_enter (2125) CDOALL 1 K=KLOW,KLOW+K1-1 DOUBLE PRECISION RX1(J1), RY1(J1), RRXY1(J1), RHO1(J1), RHOINV1(J1 X ), RHOINR1(J1), U1(J1), V1(J1), SNR1(J1), BET01(J1), TERM11(J1) X , S11(J1), S21(J1), TERM21(J1) RX1(1:$J1) = XY(JLOW:$J1,K,IX) RY1(1:$J1) = XY(JLOW:$J1,K,IY) RRXY1(1:$J1) = 1. / SQRT (RX1(1:$J1) ** 2 + RY1(1:$J1) ** 2) RX1(1:$J1) = RX1(1:$J1) * RRXY1(1:$J1) RY1(1:$J1) = RY1(1:$J1) * RRXY1(1:$J1) RHO1(1:$J1) = Q(JLOW:$J1,K,1) * XYJ(JLOW:$J1,K) RHOINV1(1:$J1) = 1. / RHO1(1:$J1) RHOINR1(1:$J1) = RHOINV1(1:$J1) * XYJ(JLOW:$J1,K) U1(1:$J1) = Q(JLOW:$J1,K,2) * RHOINR1(1:$J1) V1(1:$J1) = Q(JLOW:$J1,K,3) * RHOINR1(1:$J1) SNR1(1:$J1) = 1. / SNDSP(JLOW:$J1,K) BET01(1:$J1) = BT * RHOINV1(1:$J1) * SNR1(1:$J1) C TERM11(1:$J1) = GAMI * (U1(1:$J1) * S(JLOW:$J1,K,2) + V1( X 1:$J1) * S(JLOW:$J1,K,3) - S(JLOW:$J1,K,4) - 0.5 * (U1(1:$J1) ** X 2 + V1(1:$J1) ** 2) * S(JLOW:$J1,K,1)) C S11(1:$J1) = S(JLOW:$J1,K,1) + TERM11(1:$J1) * SNR1(1:$J1 X ) ** 2 S21(1:$J1) = ((-RY1(1:$J1) * U1(1:$J1) + RX1(1:$J1) * V1( X 1:$J1)) * S(JLOW:$J1,K,1) + RY1(1:$J1) * S(JLOW:$J1,K,2) - RX1(1 X :$J1) * S(JLOW:$J1,K,3)) * RHOINV1(1:$J1) C TERM11(1:$J1) = TERM11(1:$J1) * BET01(1:$J1) TERM21(1:$J1) = ((RX1(1:$J1) * U1(1:$J1) + RY1(1:$J1) * X V1(1:$J1)) * S(JLOW:$J1,K,1) - RX1(1:$J1) * S(JLOW:$J1,K,2) - X RY1(1:$J1) * S(JLOW:$J1,K,3)) * BT * RHOINV1(1:$J1) C S(JLOW:$J1,K,1) = S11(1:$J1) S(JLOW:$J1,K,2) = S21(1:$J1) S(JLOW:$J1,K,3) = -(TERM11(1:$J1) + TERM21(1:$J1)) S(JLOW:$J1,K,4) = -(TERM11(1:$J1) - TERM21(1:$J1)) C 1 CONTINUE C d call trace_exit (2126) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C*************** SCALAR TRIDIAGONAL ************************************ C*********************************************************************** SUBROUTINE TRIB(JDIM,A,B,C,X,F,NL,NU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C DIMENSION A(JDIM),B(JDIM),C(JDIM),X(JDIM),F(JDIM) C X(NL)=C(NL)/B(NL) F(NL)=F(NL)/B(NL) NLP1 = NL +1 d call trace_enter (2127) DO 1 J=NLP1,NU Z = 1. / (B(J) - A(J) * X(J-1)) X(J) = C(J) * Z F(J)=(F(J)-A(J)*F(J-1))*Z 1 CONTINUE C d call trace_exit (2128) NUPNL=NU+NL d call trace_enter (2129) DO 2 J1=NLP1,NU J=NUPNL-J1 F(J) = F(J) - X(J) * F(J+1) 2 CONTINUE C d call trace_exit (2130) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C*************** SCALAR PERIDC TRIDIAGONAL *************************** C*********************************************************************** SUBROUTINE TRIP(JDIM,A,B,C,F,Q,S,J1,J2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(JDIM),B(JDIM),C(JDIM),F(JDIM),Q(JDIM),S(JDIM) INTEGER I1, I2, I3, I4, I5 JA = J1 + 1 FN = F(J2) C C FORWARD ELIMINATION SWEEP C Q(J1) = -C(J1)/B(J1) F(J1) = F(J1)/B(J1) S(J1) = - A(J1)/B(J1) d call trace_enter (2131) DO 10 J=JA,J2 P = 1. / (B(J) + A(J) * Q(J-1)) Q(J) = -C(J) * P F(J) = ( F(J) - A(J)*F(J-1))*P S(J) = - A(J)*S(J-1)*P 10 CONTINUE C C BACKWARD PASS C d call trace_exit (2132) JJ = J1 + J2 Q(J2) = 0. S(J2) = 1. d call trace_enter (2133) DO 11 I=JA,J2 J = JJ - I S(J) = S(J) + Q(J) * S(J+1) Q(J) = F(J) + Q(J)*Q(J+1) 11 CONTINUE d call trace_exit (2134) F(J2) = ( FN - C(J2)*Q(J1) - A(J2)*Q(J2-1))/(C(J2)*S(J1) + 1 A(J2)*S(J2-1) +B(J2)) C C BACKWARD ELIMINATION PASS C I2 = (J2 - J1 + 31) / 32 d call trace_enter (2135) CDOALL 12 I1=0,32*I2-32,32 INTEGER I6, I7, J I7 = J1 + I1 + 1 I6 = MIN0 (J1 + I1 + 1 + 32 - 1, J2) - (J1 + I1) F(JJ-I7:$I6:(-1)) = F(J2) * S(JJ-I7:$I6:(-1)) + Q(JJ-I7:$I6:(-1 X )) 12 CONTINUE d call trace_exit (2136) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C***************** SOLUTION UPDATE ROUTINE **************************** C*********************************************************************** SUBROUTINE UPDATE(JDIM,KDIM,Q,S,DELTQN) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C C DIMENSION Q(JDIM,KDIM,4),S(JDIM,KDIM,4),DELTQN(JDIM,KDIM,4) INTEGER J1 C C PHIC = PHIDT/(1.+PHIDT) J1 = JUP - JLOW + 1 C C UPDATE SOLUTION C d call trace_enter (2137) CDOALL 600 N=1,4 INTEGER K DO 600 K=KLOW,KUP C DELTA Q N DELTQN(JLOW:$J1,K,N) = S(JLOW:$J1,K,N) C UPDATE Q(JLOW:$J1,K,N) = Q(JLOW:$J1,K,N) + S(JLOW:$J1,K,N) C START NEXT S(JLOW:$J1,K,N) = PHIC * S(JLOW:$J1,K,N) 600 CONTINUE C C d call trace_exit (2138) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 SUBROUTINE VERIFY(NUMITE,RSD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C VERIFY THE CORRECTNESS OF THE OUTPUT USING THE 100TH ITERATION. C THE VERIFICATION PARAMETER IS THE 2-NORM OF THE RESDIUAL. COMMON/IO/IREAD,IWRIT,IVER IF (NUMITE.EQ.100) THEN C WRITE(IVER,1) 1 FORMAT(1X,'OUTPUT FOR PERFECT CLUB BENCHMARK: ARC2D'/ * 1X,'$Revision: 1.1 $ $Author: sharma $'/) WRITE(IVER,*) WRITE(IVER,2) 2 FORMAT(1X,'VERIFICATION PARAMETERS:'/) C WRITE(IVER,3) NUMITE,RSD 3 FORMAT(1X,'ITERATION = ',I5,3X,'L2 RESIDUAL = ',E12.6) C IF (ABS((RSD-0.142979E-04)/.142979E-04).LT.0.1E-03) THEN WRITE(IVER,4) 4 FORMAT(//1X,'RESULTS FOR THIS RUN ARE: VALID') ELSE WRITE(IVER,5) 5 FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID') ENDIF ENDIF RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C********************************************************************** C********************** DIFFERENCE X, Y IN XI ************************* C*********************************************************************** SUBROUTINE XIDIF(JDIM,KDIM,X,Y,XY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, 1 JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, 1 JUP, 1 KLOW, KUP, PERIDC , NP, DT, 1 CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION XY(JDIM,KDIM,4) DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM) INTEGER J1, J2, J3, J4, J5 C C C XY4 =XXI,XY3 = YXI C d call trace_enter (2139) DO 11 K=KBEGIN,KEND J2 = (JUP - JLOW + 32) / 32 CDOALL 10 J1=0,32*J2-32,32 INTEGER J6, J7, JM1, JP1 J7 = JLOW + J1 J6 = MIN0 (JLOW + J1 + 32 - 1, JUP) - (JLOW + J1) + 1 XY(J7:$J6,K,4) = (X(JPLUS(J7:$J6),K) - X(JMINU(J7:$J6),K)) * X .5 XY(J7:$J6,K,3) = (Y(JPLUS(J7:$J6),K) - Y(JMINU(J7:$J6),K)) * X .5 10 CONTINUE C IF(.NOT.PERIDC )THEN J = JBEGIN XY(J,K,4) = .5 * (-3. * X(J,K) + 4. * X(J+1,K) - X(J+2,K)) XY(J,K,3) = .5 * (-3. * Y(J,K) + 4. * Y(J+1,K) - Y(J+2,K)) J = JEND XY(J,K,4) = (3. * X(J,K) - 4. * X(J-1,K) + X(J-2,K)) * .5 XY(J,K,3) = (3. * Y(J,K) - 4. * Y(J-1,K) + Y(J-2,K)) * .5 ENDIF C 11 CONTINUE C d call trace_exit (2140) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C********************************************************************* C********************************************************************* C************ VECTORIZED IN K SCALAR PERIDC PENTA SOLVER ********* C DOES TWO F'S AT A TIME ********* C********************************************************************* C********************************************************************* SUBROUTINE XPENP2(JDIM,KDIM,A,B,C,D,E,U1,U2,U3,U4,F, * JL,JU,KL,KU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) DOUBLE PRECISION LD,LD1,LD2,LDI,LDD,LU DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM) DIMENSION E(JDIM,KDIM),F(JDIM,KDIM,2),U4(JDIM,KDIM) DIMENSION U1(JDIM,KDIM),U2(JDIM,KDIM),U3(JDIM,KDIM) C COMMON/WORKSP/LDA(MAXJ),LDB(MAXJ),LUS(MAXJ), * LDS(MAXJ),WORK(MAXJ,88) DOUBLE PRECISION LDA,LDB,LUS,LDS INTEGER K1, K2, K3, K4, K5, K6, K7, K8, K9, K10, K11, K12, K13, X K14, K15, K16, K17, K20, K23, K26, K29, K32, K35, K38, K41, X K44, K47, K50, K53 C C ! START FORWARD GENERATION PROCESS AND SWEEP C I = JL K14 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2141) CDOALL 10 K1=0,32*K14-32,32 DOUBLE PRECISION LD3(K17), LDI1(K17) INTEGER K18, K19 K19 = KL + K1 K18 = MIN0 (KL + K1 + 32 - 1, KU) - (KL + K1) + 1 LD3(1:$K18) = C(I,K19:$K18) LDI1(1:$K18) = 1. / LD3(1:$K18) U1(I,K19:$K18) = D(I,K19:$K18) * LDI1(1:$K18) U2(I,K19:$K18) = E(I,K19:$K18) * LDI1(1:$K18) U3(I,K19:$K18) = A(I,K19:$K18) * LDI1(1:$K18) U4(I,K19:$K18) = B(I,K19:$K18) * LDI1(1:$K18) F(I,K19:$K18,1) = F(I,K19:$K18,1) * LDI1(1:$K18) F(I,K19:$K18,2) = F(I,K19:$K18,2) * LDI1(1:$K18) 10 CONTINUE C d call trace_exit (2142) I = JL+1 K20 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2143) CDOALL 20 K2=0,32*K20-32,32 DOUBLE PRECISION LD11(K17), LD4(K17), LDI2(K17) INTEGER K21, K22 K22 = KL + K2 K21 = MIN0 (KL + K2 + 32 - 1, KU) - (KL + K2) + 1 LD11(1:$K21) = B(I,K22:$K21) LD4(1:$K21) = C(I,K22:$K21) - LD11(1:$K21) * U1(I-1,K22:$K21) LDI2(1:$K21) = 1. / LD4(1:$K21) U1(I,K22:$K21) = (D(I,K22:$K21) - LD11(1:$K21) * U2(I-1,K22:$K21 X )) * LDI2(1:$K21) U2(I,K22:$K21) = E(I,K22:$K21) * LDI2(1:$K21) U3(I,K22:$K21) = -LD11(1:$K21) * U3(I-1,K22:$K21) * LDI2(1:$K21) U4(I,K22:$K21) = (A(I,K22:$K21) - LD11(1:$K21) * U4(I-1,K22:$K21 X )) * LDI2(1:$K21) F(I,K22:$K21,1) = (F(I,K22:$K21,1) - LD11(1:$K21) * F(I-1,K22:$ X K21,1)) * LDI2(1:$K21) F(I,K22:$K21,2) = (F(I,K22:$K21,2) - LD11(1:$K21) * F(I-1,K22:$ X K21,2)) * LDI2(1:$K21) 20 CONTINUE d call trace_exit (2144) K15 = KU - KL + 1 d call trace_enter (2145) CDOALL 11 K=KL,KL+K15-1 INTEGER I DOUBLE PRECISION LDI, LD2, LD1, LD C DO 1 I=JL+2,JU-4 LD2 = A(I,K) LD1 = B(I,K) - LD2 * U1(I-2,K) LD = C(I,K) - (LD2 * U2(I-2,K) + LD1 * U1(I-1,K)) LDI = 1. / LD U1(I,K) = (D(I,K) - LD1 * U2(I-1,K)) * LDI U2(I,K) = E(I,K) * LDI U3(I,K) = (-LD2 * U3(I-2,K) - LD1 * U3(I-1,K)) * LDI U4(I,K) = (-LD2 * U4(I-2,K) - LD1 * U4(I-1,K)) * LDI F(I,K,1) = (F(I,K,1) - LD2 * F(I-2,K,1) - LD1 * F(I-1,K, X 1)) * LDI F(I,K,2) = (F(I,K,2) - LD2 * F(I-2,K,2) - LD1 * F(I-1,K, X 2)) * LDI 1 CONTINUE 11 CONTINUE C d call trace_exit (2146) I = JU-3 K23 = (KU - KL + 32) / 32 K15 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2147) CDOALL 12 K3=0,32*K23-32,32 DOUBLE PRECISION LD22(K15), LD13(K15), LD6(K15), LDI4( X K15) INTEGER K24, K25 K25 = KL + K3 K24 = MIN0 (KL + K3 + 32 - 1, KU) - (KL + K3) + 1 LD22(1:$K24) = A(I,K25:$K24) LD13(1:$K24) = B(I,K25:$K24) - LD22(1:$K24) * U1(I-2,K25 X :$K24) LD6(1:$K24) = C(I,K25:$K24) - (LD22(1:$K24) * U2(I-2,K25 X :$K24) + LD13(1:$K24) * U1(I-1,K25:$K24)) LDI4(1:$K24) = 1. / LD6(1:$K24) U1(I,K25:$K24) = (D(I,K25:$K24) - LD13(1:$K24) * U2(I-1, X K25:$K24)) * LDI4(1:$K24) U3(I,K25:$K24) = (E(I,K25:$K24) - LD22(1:$K24) * U3(I-2, X K25:$K24) - LD13(1:$K24) * U3(I-1,K25:$K24)) * LDI4(1 X :$K24) U4(I,K25:$K24) = (-LD22(1:$K24) * U4(I-2,K25:$K24) - X LD13(1:$K24) * U4(I-1,K25:$K24)) * LDI4(1:$K24) F(I,K25:$K24,1) = (F(I,K25:$K24,1) - LD22(1:$K24) * F(I- X 2,K25:$K24,1) - LD13(1:$K24) * F(I-1,K25:$K24,1)) * LDI4(1:$K24) F(I,K25:$K24,2) = (F(I,K25:$K24,2) - LD22(1:$K24) * F(I- X 2,K25:$K24,2) - LD13(1:$K24) * F(I-1,K25:$K24,2)) * LDI4(1:$K24) 12 CONTINUE C d call trace_exit (2148) I = JU-2 K26 = (KU - KL + 32) / 32 K15 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2149) CDOALL 22 K4=0,32*K26-32,32 DOUBLE PRECISION LD23(K15), LD14(K15), LD7(K15), LDI5(K15) INTEGER K27, K28 K28 = KL + K4 K27 = MIN0 (KL + K4 + 32 - 1, KU) - (KL + K4) + 1 LD23(1:$K27) = A(I,K28:$K27) LD14(1:$K27) = B(I,K28:$K27) - LD23(1:$K27) * U1(I-2,K28:$K27) LD7(1:$K27) = C(I,K28:$K27) - (LD23(1:$K27) * U2(I-2,K28:$K27) + X LD14(1:$K27) * U1(I-1,K28:$K27)) LDI5(1:$K27) = 1. / LD7(1:$K27) U3(I,K28:$K27) = (D(I,K28:$K27) - LD23(1:$K27) * U3(I-2,K28:$K27 X ) - LD14(1:$K27) * U3(I-1,K28:$K27)) * LDI5(1:$K27) U4(I,K28:$K27) = (E(I,K28:$K27) - LD23(1:$K27) * U4(I-2,K28:$K27 X ) - LD14(1:$K27) * U4(I-1,K28:$K27)) * LDI5(1:$K27) F(I,K28:$K27,1) = (F(I,K28:$K27,1) - LD23(1:$K27) * F(I-2,K28:$ X K27,1) - LD14(1:$K27) * F(I-1,K28:$K27,1)) * LDI5(1:$K27) F(I,K28:$K27,2) = (F(I,K28:$K27,2) - LD23(1:$K27) * F(I-2,K28:$ X K27,2) - LD14(1:$K27) * F(I-1,K28:$K27,2)) * LDI5(1:$K27) 22 CONTINUE C d call trace_exit (2150) I = JU-1 K29 = (KU - KL + 32) / 32 d call trace_enter (2151) CDOALL 32 K5=0,32*K29-32,32 INTEGER K30, K31 K31 = KL + K5 K30 = MIN0 (KL + K5 + 32 - 1, KU) - (KL + K5) + 1 LDA(K31:$K30) = E(I,K31:$K30) LDB(K31:$K30) = -LDA(K31:$K30) * U1(JL,K31:$K30) F(I,K31:$K30,1) = F(I,K31:$K30,1) - LDA(K31:$K30) * F(JL,K31:$ X K30,1) - LDB(K31:$K30) * F(JL+1,K31:$K30,1) F(I,K31:$K30,2) = F(I,K31:$K30,2) - LDA(K31:$K30) * F(JL,K31:$ X K30,2) - LDB(K31:$K30) * F(JL+1,K31:$K30,2) LDS(K31:$K30) = LDA(K31:$K30) * U3(JL,K31:$K30) + LDB(K31:$K30) X * U3(JL+1,K31:$K30) U4(I,K31:$K30) = LDA(K31:$K30) * U4(JL,K31:$K30) + LDB(K31:$K30) X * U4(JL+1,K31:$K30) 32 CONTINUE d call trace_exit (2152) K16 = KU - KL + 1 d call trace_enter (2153) CDOALL 14 K=KL,KL+K16-1 INTEGER M DOUBLE PRECISION LDD DO 3 M=JL+2,JU-4 LDD = -LDA(K) * U2(M-2,K) - LDB(K) * U1(M-1,K) LDA(K) = LDB(K) LDB(K) = LDD LDS(K) = LDS(K) + LDD * U3(M,K) U4(I,K) = U4(I,K) + LDD * U4(M,K) F(I,K,1) = F(I,K,1) - LDD * F(M,K,1) F(I,K,2) = F(I,K,2) - LDD * F(M,K,2) 3 CONTINUE 14 CONTINUE C d call trace_exit (2154) M = JU-3 K32 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2155) CDOALL 15 K6=0,32*K32-32,32 DOUBLE PRECISION LDD2(K16) INTEGER K33, K34 K34 = KL + K6 K33 = MIN0 (KL + K6 + 32 - 1, KU) - (KL + K6) + 1 LDD2(1:$K33) = A(I,K34:$K33) - LDA(K34:$K33) * U2(M-2,K34:$K33) X - LDB(K34:$K33) * U1(M-1,K34:$K33) LDA(K34:$K33) = LDB(K34:$K33) LDB(K34:$K33) = LDD2(1:$K33) LDS(K34:$K33) = LDS(K34:$K33) + LDD2(1:$K33) * U3(M,K34:$K33) U4(I,K34:$K33) = U4(I,K34:$K33) + LDD2(1:$K33) * U4(M,K34:$K33) F(I,K34:$K33,1) = F(I,K34:$K33,1) - LDD2(1:$K33) * F(M,K34:$K33, X 1) F(I,K34:$K33,2) = F(I,K34:$K33,2) - LDD2(1:$K33) * F(M,K34:$K33, X 2) 15 CONTINUE d call trace_exit (2156) M = JU-2 K35 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2157) CDOALL 25 K7=0,32*K35-32,32 DOUBLE PRECISION LDD3(K16) INTEGER K36, K37 K37 = KL + K7 K36 = MIN0 (KL + K7 + 32 - 1, KU) - (KL + K7) + 1 LDD3(1:$K36) = B(I,K37:$K36) - LDA(K37:$K36) * U2(M-2,K37:$K36) X - LDB(K37:$K36) * U1(M-1,K37:$K36) LDA(K37:$K36) = LDB(K37:$K36) LDB(K37:$K36) = LDD3(1:$K36) LDS(K37:$K36) = LDS(K37:$K36) + LDD3(1:$K36) * U3(M,K37:$K36) U4(I,K37:$K36) = U4(I,K37:$K36) + LDD3(1:$K36) * U4(M,K37:$K36) F(I,K37:$K36,1) = F(I,K37:$K36,1) - LDD3(1:$K36) * F(M,K37:$K36, X 1) F(I,K37:$K36,2) = F(I,K37:$K36,2) - LDD3(1:$K36) * F(M,K37:$K36, X 2) 25 CONTINUE C d call trace_exit (2158) K38 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2159) CDOALL 35 K8=0,32*K38-32,32 DOUBLE PRECISION LDI6(K16) INTEGER K39, K40 K40 = KL + K8 K39 = MIN0 (KL + K8 + 32 - 1, KU) - (KL + K8) + 1 LDS(K40:$K39) = C(I,K40:$K39) - LDS(K40:$K39) LDI6(1:$K39) = 1. / LDS(K40:$K39) U4(I,K40:$K39) = (D(I,K40:$K39) - U4(I,K40:$K39)) * LDI6(1:$K39) F(I,K40:$K39,1) = F(I,K40:$K39,1) * LDI6(1:$K39) F(I,K40:$K39,2) = F(I,K40:$K39,2) * LDI6(1:$K39) 35 CONTINUE C C THIS BREAKUP OF K LOOPS IS VERY IMPORTANT C d call trace_exit (2160) I = JU K41 = (KU - KL + 32) / 32 d call trace_enter (2161) CDOALL 16 K9=0,32*K41-32,32 INTEGER K42, K43 K43 = KL + K9 K42 = MIN0 (KL + K9 + 32 - 1, KU) - (KL + K9) + 1 LDA(K43:$K42) = D(I,K43:$K42) LDB(K43:$K42) = E(I,K43:$K42) - LDA(K43:$K42) * U1(JL,K43:$K42) LUS(K43:$K42) = LDA(K43:$K42) * U3(JL,K43:$K42) + LDB(K43:$K42) X * U3(JL+1,K43:$K42) LDS(K43:$K42) = LDA(K43:$K42) * U4(JL,K43:$K42) + LDB(K43:$K42) X * U4(JL+1,K43:$K42) F(I,K43:$K42,1) = F(I,K43:$K42,1) - LDA(K43:$K42) * F(JL,K43:$ X K42,1) - LDB(K43:$K42) * F(JL+1,K43:$K42,1) F(I,K43:$K42,2) = F(I,K43:$K42,2) - LDA(K43:$K42) * F(JL,K43:$ X K42,2) - LDB(K43:$K42) * F(JL+1,K43:$K42,2) C 16 CONTINUE d call trace_exit (2162) K17 = KU - KL + 1 d call trace_enter (2163) CDOALL 17 K=KL,KL+K17-1 INTEGER M DOUBLE PRECISION LDD C DO 4 M=JL+2,JU-3 LDD = -LDA(K) * U2(M-2,K) - LDB(K) * U1(M-1,K) LDA(K) = LDB(K) LDB(K) = LDD LUS(K) = LUS(K) + LDD * U3(M,K) LDS(K) = LDS(K) + LDD * U4(M,K) F(I,K,1) = F(I,K,1) - LDD * F(M,K,1) F(I,K,2) = F(I,K,2) - LDD * F(M,K,2) 4 CONTINUE 17 CONTINUE C d call trace_exit (2164) M = JU-2 K44 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2165) CDOALL 18 K10=0,32*K44-32,32 DOUBLE PRECISION LDD5(K17) INTEGER K45, K46 K46 = KL + K10 K45 = MIN0 (KL + K10 + 32 - 1, KU) - (KL + K10) + 1 LDD5(1:$K45) = A(I,K46:$K45) - LDA(K46:$K45) * U2(M-2,K46:$K45) X - LDB(K46:$K45) * U1(M-1,K46:$K45) LDA(K46:$K45) = LDB(K46:$K45) LDB(K46:$K45) = LDD5(1:$K45) LUS(K46:$K45) = LUS(K46:$K45) + LDD5(1:$K45) * U3(M,K46:$K45) LDS(K46:$K45) = LDS(K46:$K45) + LDD5(1:$K45) * U4(M,K46:$K45) F(I,K46:$K45,1) = F(I,K46:$K45,1) - LDD5(1:$K45) * F(M,K46:$K45, X 1) F(I,K46:$K45,2) = F(I,K46:$K45,2) - LDD5(1:$K45) * F(M,K46:$K45, X 2) 18 CONTINUE d call trace_exit (2166) M = JU-1 K47 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2167) CDOALL 28 K11=0,32*K47-32,32 DOUBLE PRECISION LDD6(K17), LDI7(K17) INTEGER K48, K49 K49 = KL + K11 K48 = MIN0 (KL + K11 + 32 - 1, KU) - (KL + K11) + 1 LDD6(1:$K48) = B(I,K49:$K48) - LUS(K49:$K48) F(I,K49:$K48,1) = F(I,K49:$K48,1) - LDD6(1:$K48) * F(M,K49:$K48, X 1) F(I,K49:$K48,2) = F(I,K49:$K48,2) - LDD6(1:$K48) * F(M,K49:$K48, X 2) LDS(K49:$K48) = C(I,K49:$K48) - LDD6(1:$K48) * U4(M,K49:$K48) - X LDS(K49:$K48) LDI7(1:$K48) = 1. / LDS(K49:$K48) F(I,K49:$K48,1) = F(I,K49:$K48,1) * LDI7(1:$K48) F(I,K49:$K48,2) = F(I,K49:$K48,2) * LDI7(1:$K48) 28 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2168) K50 = (KU - KL + 32) / 32 d call trace_enter (2169) CDOALL 19 K12=0,32*K50-32,32 INTEGER K51, K52 K52 = KL + K12 K51 = MIN0 (KL + K12 + 32 - 1, KU) - (KL + K12) + 1 F(JU,K52:$K51,1) = F(JU,K52:$K51,1) F(JU,K52:$K51,2) = F(JU,K52:$K51,2) F(JU-1,K52:$K51,1) = F(JU-1,K52:$K51,1) - U4(JU-1,K52:$K51) * F(JU X ,K52:$K51,1) F(JU-1,K52:$K51,2) = F(JU-1,K52:$K51,2) - U4(JU-1,K52:$K51) * F(JU X ,K52:$K51,2) F(JU-2,K52:$K51,1) = F(JU-2,K52:$K51,1) - U4(JU-2,K52:$K51) * F(JU X ,K52:$K51,1) - U3(JU-2,K52:$K51) * F(JU-1,K52:$K51,1) F(JU-2,K52:$K51,2) = F(JU-2,K52:$K51,2) - U4(JU-2,K52:$K51) * F(JU X ,K52:$K51,2) - U3(JU-2,K52:$K51) * F(JU-1,K52:$K51,2) F(JU-3,K52:$K51,1) = F(JU-3,K52:$K51,1) - U4(JU-3,K52:$K51) * F(JU X ,K52:$K51,1) - U3(JU-3,K52:$K51) * F(JU-1,K52:$K51,1) - U1(JU-3, X K52:$K51) * F(JU-2,K52:$K51,1) F(JU-3,K52:$K51,2) = F(JU-3,K52:$K51,2) - U4(JU-3,K52:$K51) * F(JU X ,K52:$K51,2) - U3(JU-3,K52:$K51) * F(JU-1,K52:$K51,2) - U1(JU-3, X K52:$K51) * F(JU-2,K52:$K51,2) 19 CONTINUE C d call trace_exit (2170) d call trace_enter (2171) DO 2 I = 4,JU-JL IX = JU-I K53 = (KU - KL + 32) / 32 CDOALL 50 K13=0,32*K53-32,32 INTEGER K54, K55 K55 = KL + K13 K54 = MIN0 (KL + K13 + 32 - 1, KU) - (KL + K13) + 1 F(JU-I,K55:$K54,1) = F(JU-I,K55:$K54,1) - U4(JU-I,K55:$K54) * X F(JU,K55:$K54,1) - U3(JU-I,K55:$K54) * F(JU-1,K55:$K54,1) - U1( X JU-I,K55:$K54) * F(JU-I+1,K55:$K54,1) - U2(JU-I,K55:$K54) * F(JU X -I+2,K55:$K54,1) F(JU-I,K55:$K54,2) = F(JU-I,K55:$K54,2) - U4(JU-I,K55:$K54) * X F(JU,K55:$K54,2) - U3(JU-I,K55:$K54) * F(JU-1,K55:$K54,2) - U1( X JU-I,K55:$K54) * F(JU-I+1,K55:$K54,2) - U2(JU-I,K55:$K54) * F(JU X -I+2,K55:$K54,2) 50 CONTINUE 2 CONTINUE C d call trace_exit (2172) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C***************************************************************** C***************************************************************** C ******** VECTORIZED IN K SCALAR PENTA SOLVER *************** C SOLVES FOR TWO F'S AT A TIME *************** C***************************************************************** C***************************************************************** SUBROUTINE XPENT2(JDIM,KDIM,A,B,C,D,E,X,Y,F,JL,JU,KL,KU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LD,LD1,LD2,LDI DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM), * E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM,2) INTEGER K1, K2, K3, K4, K5, K6, K7, K8, K9, K10, K13, K16, K19, X K22, K25 C C ! START FORWARD GENERATION PROCESS AND SWEEP C J = JL K7 = (KU - KL + 32) / 32 K10 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2173) CDOALL 1 K1=0,32*K7-32,32 DOUBLE PRECISION LD3(K10), LDI1(K10) INTEGER K11, K12 K12 = KL + K1 K11 = MIN0 (KL + K1 + 32 - 1, KU) - (KL + K1) + 1 LD3(1:$K11) = C(J,K12:$K11) LDI1(1:$K11) = 1. / LD3(1:$K11) F(J,K12:$K11,1) = F(J,K12:$K11,1) * LDI1(1:$K11) F(J,K12:$K11,2) = F(J,K12:$K11,2) * LDI1(1:$K11) X(J,K12:$K11) = D(J,K12:$K11) * LDI1(1:$K11) Y(J,K12:$K11) = E(J,K12:$K11) * LDI1(1:$K11) 1 CONTINUE C d call trace_exit (2174) J = JL+1 K13 = (KU - KL + 32) / 32 K10 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2175) CDOALL 2 K2=0,32*K13-32,32 DOUBLE PRECISION LD11(K10), LD4(K10), LDI2(K10) INTEGER K14, K15 K15 = KL + K2 K14 = MIN0 (KL + K2 + 32 - 1, KU) - (KL + K2) + 1 LD11(1:$K14) = B(J,K15:$K14) LD4(1:$K14) = C(J,K15:$K14) - LD11(1:$K14) * X(J-1,K15:$K14) LDI2(1:$K14) = 1. / LD4(1:$K14) F(J,K15:$K14,1) = (F(J,K15:$K14,1) - LD11(1:$K14) * F(J-1,K15:$ X K14,1)) * LDI2(1:$K14) F(J,K15:$K14,2) = (F(J,K15:$K14,2) - LD11(1:$K14) * F(J-1,K15:$ X K14,2)) * LDI2(1:$K14) X(J,K15:$K14) = (D(J,K15:$K14) - LD11(1:$K14) * Y(J-1,K15:$K14)) X * LDI2(1:$K14) Y(J,K15:$K14) = E(J,K15:$K14) * LDI2(1:$K14) 2 CONTINUE d call trace_exit (2176) K8 = KU - KL + 1 d call trace_enter (2177) CDOALL 11 K=KL,KL+K8-1 INTEGER J DOUBLE PRECISION LDI, LD2, LD1, LD C DO 3 J=JL+2,JU-2 LD2 = A(J,K) LD1 = B(J,K) - LD2 * X(J-2,K) LD = C(J,K) - (LD2 * Y(J-2,K) + LD1 * X(J-1,K)) LDI = 1. / LD F(J,K,1) = (F(J,K,1) - LD2 * F(J-2,K,1) - LD1 * F(J-1,K, X 1)) * LDI F(J,K,2) = (F(J,K,2) - LD2 * F(J-2,K,2) - LD1 * F(J-1,K, X 2)) * LDI X(J,K) = (D(J,K) - LD1 * Y(J-1,K)) * LDI Y(J,K) = E(J,K) * LDI 3 CONTINUE 11 CONTINUE C d call trace_exit (2178) J = JU-1 K16 = (KU - KL + 32) / 32 K8 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2179) CDOALL 12 K3=0,32*K16-32,32 DOUBLE PRECISION LD22(K8), LD13(K8), LD6(K8), LDI4(K8) INTEGER K17, K18 K18 = KL + K3 K17 = MIN0 (KL + K3 + 32 - 1, KU) - (KL + K3) + 1 LD22(1:$K17) = A(J,K18:$K17) LD13(1:$K17) = B(J,K18:$K17) - LD22(1:$K17) * X(J-2,K18:$K17) LD6(1:$K17) = C(J,K18:$K17) - (LD22(1:$K17) * Y(J-2,K18:$K17) + X LD13(1:$K17) * X(J-1,K18:$K17)) LDI4(1:$K17) = 1. / LD6(1:$K17) F(J,K18:$K17,1) = (F(J,K18:$K17,1) - LD22(1:$K17) * F(J-2,K18:$ X K17,1) - LD13(1:$K17) * F(J-1,K18:$K17,1)) * LDI4(1:$K17) F(J,K18:$K17,2) = (F(J,K18:$K17,2) - LD22(1:$K17) * F(J-2,K18:$ X K17,2) - LD13(1:$K17) * F(J-1,K18:$K17,2)) * LDI4(1:$K17) X(J,K18:$K17) = (D(J,K18:$K17) - LD13(1:$K17) * Y(J-1,K18:$K17)) X * LDI4(1:$K17) 12 CONTINUE C d call trace_exit (2180) J = JU K19 = (KU - KL + 32) / 32 K8 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2181) CDOALL 13 K4=0,32*K19-32,32 DOUBLE PRECISION LD23(K8), LD14(K8), LD7(K8), LDI5(K8) INTEGER K20, K21 K21 = KL + K4 K20 = MIN0 (KL + K4 + 32 - 1, KU) - (KL + K4) + 1 LD23(1:$K20) = A(J,K21:$K20) LD14(1:$K20) = B(J,K21:$K20) - LD23(1:$K20) * X(J-2,K21 X :$K20) LD7(1:$K20) = C(J,K21:$K20) - (LD23(1:$K20) * Y(J-2,K21 X :$K20) + LD14(1:$K20) * X(J-1,K21:$K20)) LDI5(1:$K20) = 1. / LD7(1:$K20) F(J,K21:$K20,1) = (F(J,K21:$K20,1) - LD23(1:$K20) * F(J- X 2,K21:$K20,1) - LD14(1:$K20) * F(J-1,K21:$K20,1)) * LDI5(1:$K20) F(J,K21:$K20,2) = (F(J,K21:$K20,2) - LD23(1:$K20) * F(J- X 2,K21:$K20,2) - LD14(1:$K20) * F(J-1,K21:$K20,2)) * LDI5(1:$K20) 13 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2182) K22 = (KU - KL + 32) / 32 d call trace_enter (2183) CDOALL 14 K5=0,32*K22-32,32 INTEGER K23, K24 K24 = KL + K5 K23 = MIN0 (KL + K5 + 32 - 1, KU) - (KL + K5) + 1 F(JU,K24:$K23,1) = F(JU,K24:$K23,1) F(JU,K24:$K23,2) = F(JU,K24:$K23,2) F(JU-1,K24:$K23,1) = F(JU-1,K24:$K23,1) - X(JU-1,K24:$K23) * F( X JU,K24:$K23,1) F(JU-1,K24:$K23,2) = F(JU-1,K24:$K23,2) - X(JU-1,K24:$K23) * F( X JU,K24:$K23,2) 14 CONTINUE C d call trace_exit (2184) d call trace_enter (2185) DO 4 J = 2,JU-JL JX = JU-J K25 = (KU - KL + 32) / 32 CDOALL 15 K6=0,32*K25-32,32 INTEGER K26, K27 K27 = KL + K6 K26 = MIN0 (KL + K6 + 32 - 1, KU) - (KL + K6) + 1 F(JU-J,K27:$K26,1) = F(JU-J,K27:$K26,1) - X(JU-J,K27:$K26) * X F(JU-J+1,K27:$K26,1) - Y(JU-J,K27:$K26) * F(JU-J+2,K27:$K26,1) F(JU-J,K27:$K26,2) = F(JU-J,K27:$K26,2) - X(JU-J,K27:$K26) * X F(JU-J+1,K27:$K26,2) - Y(JU-J,K27:$K26) * F(JU-J+2,K27:$K26,2) 15 CONTINUE 4 CONTINUE C d call trace_exit (2186) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C***************************************************************** C***************************************************************** C ******** VECTORIZED IN K SCALAR PENTA SOLVER *************** C***************************************************************** C***************************************************************** SUBROUTINE XPENTA(JDIM,KDIM,A,B,C,D,E,X,Y,F,JL,JU,KL,KU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LD,LD1,LD2,LDI DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM), * E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM) INTEGER K1, K2, K3, K4, K5, K6, K7, K8, K9, K10, K13, K16, K19, X K22, K25 C C ! START FORWARD GENERATION PROCESS AND SWEEP C J = JL K7 = (KU - KL + 32) / 32 K10 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2187) CDOALL 1 K1=0,32*K7-32,32 DOUBLE PRECISION LD3(K10), LDI1(K10) INTEGER K11, K12 K12 = KL + K1 K11 = MIN0 (KL + K1 + 32 - 1, KU) - (KL + K1) + 1 LD3(1:$K11) = C(J,K12:$K11) LDI1(1:$K11) = 1. / LD3(1:$K11) F(J,K12:$K11) = F(J,K12:$K11) * LDI1(1:$K11) X(J,K12:$K11) = D(J,K12:$K11) * LDI1(1:$K11) Y(J,K12:$K11) = E(J,K12:$K11) * LDI1(1:$K11) 1 CONTINUE C d call trace_exit (2188) J = JL+1 K13 = (KU - KL + 32) / 32 K10 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2189) CDOALL 2 K2=0,32*K13-32,32 DOUBLE PRECISION LD11(K10), LD4(K10), LDI2(K10) INTEGER K14, K15 K15 = KL + K2 K14 = MIN0 (KL + K2 + 32 - 1, KU) - (KL + K2) + 1 LD11(1:$K14) = B(J,K15:$K14) LD4(1:$K14) = C(J,K15:$K14) - LD11(1:$K14) * X(J-1,K15:$K14) LDI2(1:$K14) = 1. / LD4(1:$K14) F(J,K15:$K14) = (F(J,K15:$K14) - LD11(1:$K14) * F(J-1,K15:$K14)) X * LDI2(1:$K14) X(J,K15:$K14) = (D(J,K15:$K14) - LD11(1:$K14) * Y(J-1,K15:$K14)) X * LDI2(1:$K14) Y(J,K15:$K14) = E(J,K15:$K14) * LDI2(1:$K14) 2 CONTINUE d call trace_exit (2190) K8 = KU - KL + 1 d call trace_enter (2191) CDOALL 11 K=KL,KL+K8-1 INTEGER J DOUBLE PRECISION LDI, LD2, LD1, LD C DO 3 J=JL+2,JU-2 LD2 = A(J,K) LD1 = B(J,K) - LD2 * X(J-2,K) LD = C(J,K) - (LD2 * Y(J-2,K) + LD1 * X(J-1,K)) LDI = 1. / LD F(J,K) = (F(J,K) - LD2 * F(J-2,K) - LD1 * F(J-1,K)) * X LDI X(J,K) = (D(J,K) - LD1 * Y(J-1,K)) * LDI Y(J,K) = E(J,K) * LDI 3 CONTINUE 11 CONTINUE C d call trace_exit (2192) J = JU-1 K16 = (KU - KL + 32) / 32 K8 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2193) CDOALL 12 K3=0,32*K16-32,32 DOUBLE PRECISION LD22(K8), LD13(K8), LD6(K8), LDI4(K8) INTEGER K17, K18 K18 = KL + K3 K17 = MIN0 (KL + K3 + 32 - 1, KU) - (KL + K3) + 1 LD22(1:$K17) = A(J,K18:$K17) LD13(1:$K17) = B(J,K18:$K17) - LD22(1:$K17) * X(J-2,K18:$K17) LD6(1:$K17) = C(J,K18:$K17) - (LD22(1:$K17) * Y(J-2,K18:$K17) + X LD13(1:$K17) * X(J-1,K18:$K17)) LDI4(1:$K17) = 1. / LD6(1:$K17) F(J,K18:$K17) = (F(J,K18:$K17) - LD22(1:$K17) * F(J-2,K18:$K17) X - LD13(1:$K17) * F(J-1,K18:$K17)) * LDI4(1:$K17) X(J,K18:$K17) = (D(J,K18:$K17) - LD13(1:$K17) * Y(J-1,K18:$K17)) X * LDI4(1:$K17) 12 CONTINUE C d call trace_exit (2194) J = JU K19 = (KU - KL + 32) / 32 K8 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2195) CDOALL 13 K4=0,32*K19-32,32 DOUBLE PRECISION LD23(K8), LD14(K8), LD7(K8), LDI5(K8) INTEGER K20, K21 K21 = KL + K4 K20 = MIN0 (KL + K4 + 32 - 1, KU) - (KL + K4) + 1 LD23(1:$K20) = A(J,K21:$K20) LD14(1:$K20) = B(J,K21:$K20) - LD23(1:$K20) * X(J-2,K21 X :$K20) LD7(1:$K20) = C(J,K21:$K20) - (LD23(1:$K20) * Y(J-2,K21 X :$K20) + LD14(1:$K20) * X(J-1,K21:$K20)) LDI5(1:$K20) = 1. / LD7(1:$K20) F(J,K21:$K20) = (F(J,K21:$K20) - LD23(1:$K20) * F(J-2, X K21:$K20) - LD14(1:$K20) * F(J-1,K21:$K20)) * LDI5(1:$ X K20) 13 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2196) K22 = (KU - KL + 32) / 32 d call trace_enter (2197) CDOALL 14 K5=0,32*K22-32,32 INTEGER K23, K24 K24 = KL + K5 K23 = MIN0 (KL + K5 + 32 - 1, KU) - (KL + K5) + 1 F(JU,K24:$K23) = F(JU,K24:$K23) F(JU-1,K24:$K23) = F(JU-1,K24:$K23) - X(JU-1,K24:$K23) * F(JU, X K24:$K23) 14 CONTINUE C d call trace_exit (2198) d call trace_enter (2199) DO 4 J = 2,JU-JL JX = JU-J K25 = (KU - KL + 32) / 32 CDOALL 15 K6=0,32*K25-32,32 INTEGER K26, K27 K27 = KL + K6 K26 = MIN0 (KL + K6 + 32 - 1, KU) - (KL + K6) + 1 F(JU-J,K27:$K26) = F(JU-J,K27:$K26) - X(JU-J,K27:$K26) * F(JU X -J+1,K27:$K26) - Y(JU-J,K27:$K26) * F(JU-J+2,K27:$K26) 15 CONTINUE 4 CONTINUE C d call trace_exit (2200) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C********************************************************************* C********************************************************************* C************ VECTORIZED IN K SCALAR PERIDC PENTA SOLVER ********* C********************************************************************* C********************************************************************* SUBROUTINE XPENTP(JDIM,KDIM,A,B,C,D,E,U1,U2,U3,U4,F, * JL,JU,KL,KU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXJ=289, MAXK=89, MAXJK=MAXJ*MAXK ) DOUBLE PRECISION LD,LD1,LD2,LDI,LDD,LU DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM) DIMENSION E(JDIM,KDIM),F(JDIM,KDIM),U1(JDIM,KDIM) DIMENSION U2(JDIM,KDIM),U3(JDIM,KDIM),U4(JDIM,KDIM) C COMMON/WORKSP/LDA(MAXJ),LDB(MAXJ),LUS(MAXJ), * LDS(MAXJ),WORK(MAXJ,88) DOUBLE PRECISION LDA,LDB,LUS,LDS INTEGER K1, K2, K3, K4, K5, K6, K7, K8, K9, K10, K11, K12, K13, X K14, K15, K16, K17, K20, K23, K26, K29, K32, K35, K38, K41, X K44, K47, K50, K53 C C ! START FORWARD GENERATION PROCESS AND SWEEP C I = JL K14 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2201) CDOALL 10 K1=0,32*K14-32,32 DOUBLE PRECISION LD3(K17), LDI1(K17) INTEGER K18, K19 K19 = KL + K1 K18 = MIN0 (KL + K1 + 32 - 1, KU) - (KL + K1) + 1 LD3(1:$K18) = C(I,K19:$K18) LDI1(1:$K18) = 1. / LD3(1:$K18) U1(I,K19:$K18) = D(I,K19:$K18) * LDI1(1:$K18) U2(I,K19:$K18) = E(I,K19:$K18) * LDI1(1:$K18) U3(I,K19:$K18) = A(I,K19:$K18) * LDI1(1:$K18) U4(I,K19:$K18) = B(I,K19:$K18) * LDI1(1:$K18) F(I,K19:$K18) = F(I,K19:$K18) * LDI1(1:$K18) 10 CONTINUE C d call trace_exit (2202) I = JL+1 K20 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2203) CDOALL 20 K2=0,32*K20-32,32 DOUBLE PRECISION LD11(K17), LD4(K17), LDI2(K17) INTEGER K21, K22 K22 = KL + K2 K21 = MIN0 (KL + K2 + 32 - 1, KU) - (KL + K2) + 1 LD11(1:$K21) = B(I,K22:$K21) LD4(1:$K21) = C(I,K22:$K21) - LD11(1:$K21) * U1(I-1,K22:$K21) LDI2(1:$K21) = 1. / LD4(1:$K21) U1(I,K22:$K21) = (D(I,K22:$K21) - LD11(1:$K21) * U2(I-1,K22:$K21 X )) * LDI2(1:$K21) U2(I,K22:$K21) = E(I,K22:$K21) * LDI2(1:$K21) U3(I,K22:$K21) = -LD11(1:$K21) * U3(I-1,K22:$K21) * LDI2(1:$K21) U4(I,K22:$K21) = (A(I,K22:$K21) - LD11(1:$K21) * U4(I-1,K22:$K21 X )) * LDI2(1:$K21) F(I,K22:$K21) = (F(I,K22:$K21) - LD11(1:$K21) * F(I-1,K22:$K21)) X * LDI2(1:$K21) 20 CONTINUE d call trace_exit (2204) K15 = KU - KL + 1 d call trace_enter (2205) CDOALL 11 K=KL,KL+K15-1 INTEGER I DOUBLE PRECISION LDI, LD2, LD1, LD C DO 1 I=JL+2,JU-4 LD2 = A(I,K) LD1 = B(I,K) - LD2 * U1(I-2,K) LD = C(I,K) - (LD2 * U2(I-2,K) + LD1 * U1(I-1,K)) LDI = 1. / LD U1(I,K) = (D(I,K) - LD1 * U2(I-1,K)) * LDI U2(I,K) = E(I,K) * LDI U3(I,K) = (-LD2 * U3(I-2,K) - LD1 * U3(I-1,K)) * LDI U4(I,K) = (-LD2 * U4(I-2,K) - LD1 * U4(I-1,K)) * LDI F(I,K) = (F(I,K) - LD2 * F(I-2,K) - LD1 * F(I-1,K)) * X LDI 1 CONTINUE 11 CONTINUE C d call trace_exit (2206) I = JU-3 K23 = (KU - KL + 32) / 32 K15 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2207) CDOALL 12 K3=0,32*K23-32,32 DOUBLE PRECISION LD22(K15), LD13(K15), LD6(K15), LDI4( X K15) INTEGER K24, K25 K25 = KL + K3 K24 = MIN0 (KL + K3 + 32 - 1, KU) - (KL + K3) + 1 LD22(1:$K24) = A(I,K25:$K24) LD13(1:$K24) = B(I,K25:$K24) - LD22(1:$K24) * U1(I-2,K25 X :$K24) LD6(1:$K24) = C(I,K25:$K24) - (LD22(1:$K24) * U2(I-2,K25 X :$K24) + LD13(1:$K24) * U1(I-1,K25:$K24)) LDI4(1:$K24) = 1. / LD6(1:$K24) U1(I,K25:$K24) = (D(I,K25:$K24) - LD13(1:$K24) * U2(I-1, X K25:$K24)) * LDI4(1:$K24) U3(I,K25:$K24) = (E(I,K25:$K24) - LD22(1:$K24) * U3(I-2, X K25:$K24) - LD13(1:$K24) * U3(I-1,K25:$K24)) * LDI4(1 X :$K24) U4(I,K25:$K24) = (-LD22(1:$K24) * U4(I-2,K25:$K24) - X LD13(1:$K24) * U4(I-1,K25:$K24)) * LDI4(1:$K24) F(I,K25:$K24) = (F(I,K25:$K24) - LD22(1:$K24) * F(I-2, X K25:$K24) - LD13(1:$K24) * F(I-1,K25:$K24)) * LDI4(1:$ X K24) 12 CONTINUE C d call trace_exit (2208) I = JU-2 K26 = (KU - KL + 32) / 32 K15 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2209) CDOALL 22 K4=0,32*K26-32,32 DOUBLE PRECISION LD23(K15), LD14(K15), LD7(K15), LDI5(K15) INTEGER K27, K28 K28 = KL + K4 K27 = MIN0 (KL + K4 + 32 - 1, KU) - (KL + K4) + 1 LD23(1:$K27) = A(I,K28:$K27) LD14(1:$K27) = B(I,K28:$K27) - LD23(1:$K27) * U1(I-2,K28:$K27) LD7(1:$K27) = C(I,K28:$K27) - (LD23(1:$K27) * U2(I-2,K28:$K27) + X LD14(1:$K27) * U1(I-1,K28:$K27)) LDI5(1:$K27) = 1. / LD7(1:$K27) U3(I,K28:$K27) = (D(I,K28:$K27) - LD23(1:$K27) * U3(I-2,K28:$K27 X ) - LD14(1:$K27) * U3(I-1,K28:$K27)) * LDI5(1:$K27) U4(I,K28:$K27) = (E(I,K28:$K27) - LD23(1:$K27) * U4(I-2,K28:$K27 X ) - LD14(1:$K27) * U4(I-1,K28:$K27)) * LDI5(1:$K27) F(I,K28:$K27) = (F(I,K28:$K27) - LD23(1:$K27) * F(I-2,K28:$K27) X - LD14(1:$K27) * F(I-1,K28:$K27)) * LDI5(1:$K27) 22 CONTINUE C d call trace_exit (2210) I = JU-1 K29 = (KU - KL + 32) / 32 d call trace_enter (2211) CDOALL 32 K5=0,32*K29-32,32 INTEGER K30, K31 K31 = KL + K5 K30 = MIN0 (KL + K5 + 32 - 1, KU) - (KL + K5) + 1 LDA(K31:$K30) = E(I,K31:$K30) LDB(K31:$K30) = -LDA(K31:$K30) * U1(JL,K31:$K30) F(I,K31:$K30) = F(I,K31:$K30) - LDA(K31:$K30) * F(JL,K31:$K30) - X LDB(K31:$K30) * F(JL+1,K31:$K30) LDS(K31:$K30) = LDA(K31:$K30) * U3(JL,K31:$K30) + LDB(K31:$K30) X * U3(JL+1,K31:$K30) U4(I,K31:$K30) = LDA(K31:$K30) * U4(JL,K31:$K30) + LDB(K31:$K30) X * U4(JL+1,K31:$K30) 32 CONTINUE d call trace_exit (2212) K16 = KU - KL + 1 d call trace_enter (2213) CDOALL 14 K=KL,KL+K16-1 INTEGER M DOUBLE PRECISION LDD DO 3 M=JL+2,JU-4 LDD = -LDA(K) * U2(M-2,K) - LDB(K) * U1(M-1,K) LDA(K) = LDB(K) LDB(K) = LDD LDS(K) = LDS(K) + LDD * U3(M,K) U4(I,K) = U4(I,K) + LDD * U4(M,K) F(I,K) = F(I,K) - LDD * F(M,K) 3 CONTINUE 14 CONTINUE C d call trace_exit (2214) M = JU-3 K32 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2215) CDOALL 15 K6=0,32*K32-32,32 DOUBLE PRECISION LDD2(K16) INTEGER K33, K34 K34 = KL + K6 K33 = MIN0 (KL + K6 + 32 - 1, KU) - (KL + K6) + 1 LDD2(1:$K33) = A(I,K34:$K33) - LDA(K34:$K33) * U2(M-2,K34:$K33) X - LDB(K34:$K33) * U1(M-1,K34:$K33) LDA(K34:$K33) = LDB(K34:$K33) LDB(K34:$K33) = LDD2(1:$K33) LDS(K34:$K33) = LDS(K34:$K33) + LDD2(1:$K33) * U3(M,K34:$K33) U4(I,K34:$K33) = U4(I,K34:$K33) + LDD2(1:$K33) * U4(M,K34:$K33) F(I,K34:$K33) = F(I,K34:$K33) - LDD2(1:$K33) * F(M,K34:$K33) 15 CONTINUE d call trace_exit (2216) M = JU-2 K35 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2217) CDOALL 25 K7=0,32*K35-32,32 DOUBLE PRECISION LDD3(K16) INTEGER K36, K37 K37 = KL + K7 K36 = MIN0 (KL + K7 + 32 - 1, KU) - (KL + K7) + 1 LDD3(1:$K36) = B(I,K37:$K36) - LDA(K37:$K36) * U2(M-2,K37:$K36) X - LDB(K37:$K36) * U1(M-1,K37:$K36) LDA(K37:$K36) = LDB(K37:$K36) LDB(K37:$K36) = LDD3(1:$K36) LDS(K37:$K36) = LDS(K37:$K36) + LDD3(1:$K36) * U3(M,K37:$K36) U4(I,K37:$K36) = U4(I,K37:$K36) + LDD3(1:$K36) * U4(M,K37:$K36) F(I,K37:$K36) = F(I,K37:$K36) - LDD3(1:$K36) * F(M,K37:$K36) 25 CONTINUE C d call trace_exit (2218) K38 = (KU - KL + 32) / 32 K16 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2219) CDOALL 35 K8=0,32*K38-32,32 DOUBLE PRECISION LDI6(K16) INTEGER K39, K40 K40 = KL + K8 K39 = MIN0 (KL + K8 + 32 - 1, KU) - (KL + K8) + 1 LDS(K40:$K39) = C(I,K40:$K39) - LDS(K40:$K39) LDI6(1:$K39) = 1. / LDS(K40:$K39) U4(I,K40:$K39) = (D(I,K40:$K39) - U4(I,K40:$K39)) * LDI6(1:$K39) F(I,K40:$K39) = F(I,K40:$K39) * LDI6(1:$K39) 35 CONTINUE C C THIS BREAKUP OF K LOOPS IS VERY IMPORTANT C d call trace_exit (2220) I = JU K41 = (KU - KL + 32) / 32 d call trace_enter (2221) CDOALL 16 K9=0,32*K41-32,32 INTEGER K42, K43 K43 = KL + K9 K42 = MIN0 (KL + K9 + 32 - 1, KU) - (KL + K9) + 1 LDA(K43:$K42) = D(I,K43:$K42) LDB(K43:$K42) = E(I,K43:$K42) - LDA(K43:$K42) * U1(JL,K43:$K42) LUS(K43:$K42) = LDA(K43:$K42) * U3(JL,K43:$K42) + LDB(K43:$K42) X * U3(JL+1,K43:$K42) LDS(K43:$K42) = LDA(K43:$K42) * U4(JL,K43:$K42) + LDB(K43:$K42) X * U4(JL+1,K43:$K42) F(I,K43:$K42) = F(I,K43:$K42) - LDA(K43:$K42) * F(JL,K43:$K42) - X LDB(K43:$K42) * F(JL+1,K43:$K42) C 16 CONTINUE d call trace_exit (2222) K17 = KU - KL + 1 d call trace_enter (2223) CDOALL 17 K=KL,KL+K17-1 INTEGER M DOUBLE PRECISION LDD C DO 4 M=JL+2,JU-3 LDD = -LDA(K) * U2(M-2,K) - LDB(K) * U1(M-1,K) LDA(K) = LDB(K) LDB(K) = LDD LUS(K) = LUS(K) + LDD * U3(M,K) LDS(K) = LDS(K) + LDD * U4(M,K) F(I,K) = F(I,K) - LDD * F(M,K) 4 CONTINUE 17 CONTINUE C d call trace_exit (2224) M = JU-2 K44 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2225) CDOALL 18 K10=0,32*K44-32,32 DOUBLE PRECISION LDD5(K17) INTEGER K45, K46 K46 = KL + K10 K45 = MIN0 (KL + K10 + 32 - 1, KU) - (KL + K10) + 1 LDD5(1:$K45) = A(I,K46:$K45) - LDA(K46:$K45) * U2(M-2,K46:$K45) X - LDB(K46:$K45) * U1(M-1,K46:$K45) LDA(K46:$K45) = LDB(K46:$K45) LDB(K46:$K45) = LDD5(1:$K45) LUS(K46:$K45) = LUS(K46:$K45) + LDD5(1:$K45) * U3(M,K46:$K45) LDS(K46:$K45) = LDS(K46:$K45) + LDD5(1:$K45) * U4(M,K46:$K45) F(I,K46:$K45) = F(I,K46:$K45) - LDD5(1:$K45) * F(M,K46:$K45) 18 CONTINUE d call trace_exit (2226) M = JU-1 K47 = (KU - KL + 32) / 32 K17 = MAX0 (MIN0 (KL + 31, KU) - KL + 1, 1) d call trace_enter (2227) CDOALL 28 K11=0,32*K47-32,32 DOUBLE PRECISION LDD6(K17), LDI7(K17) INTEGER K48, K49 K49 = KL + K11 K48 = MIN0 (KL + K11 + 32 - 1, KU) - (KL + K11) + 1 LDD6(1:$K48) = B(I,K49:$K48) - LUS(K49:$K48) F(I,K49:$K48) = F(I,K49:$K48) - LDD6(1:$K48) * F(M,K49:$K48) LDS(K49:$K48) = C(I,K49:$K48) - LDD6(1:$K48) * U4(M,K49:$K48) - X LDS(K49:$K48) LDI7(1:$K48) = 1. / LDS(K49:$K48) F(I,K49:$K48) = F(I,K49:$K48) * LDI7(1:$K48) 28 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2228) K50 = (KU - KL + 32) / 32 d call trace_enter (2229) CDOALL 19 K12=0,32*K50-32,32 INTEGER K51, K52 K52 = KL + K12 K51 = MIN0 (KL + K12 + 32 - 1, KU) - (KL + K12) + 1 F(JU,K52:$K51) = F(JU,K52:$K51) F(JU-1,K52:$K51) = F(JU-1,K52:$K51) - U4(JU-1,K52:$K51) * F(JU,K52 X :$K51) F(JU-2,K52:$K51) = F(JU-2,K52:$K51) - U4(JU-2,K52:$K51) * F(JU,K52 X :$K51) - U3(JU-2,K52:$K51) * F(JU-1,K52:$K51) F(JU-3,K52:$K51) = F(JU-3,K52:$K51) - U4(JU-3,K52:$K51) * F(JU,K52 X :$K51) - U3(JU-3,K52:$K51) * F(JU-1,K52:$K51) - U1(JU-3,K52:$K51 X ) * F(JU-2,K52:$K51) 19 CONTINUE C d call trace_exit (2230) d call trace_enter (2231) DO 2 I = 4,JU-JL IX = JU-I K53 = (KU - KL + 32) / 32 CDOALL 50 K13=0,32*K53-32,32 INTEGER K54, K55 K55 = KL + K13 K54 = MIN0 (KL + K13 + 32 - 1, KU) - (KL + K13) + 1 F(JU-I,K55:$K54) = F(JU-I,K55:$K54) - U4(JU-I,K55:$K54) * F(JU X ,K55:$K54) - U3(JU-I,K55:$K54) * F(JU-1,K55:$K54) - U1(JU-I,K55 X :$K54) * F(JU-I+1,K55:$K54) - U2(JU-I,K55:$K54) * F(JU-I+2,K55:$ X K54) 50 CONTINUE 2 CONTINUE C d call trace_exit (2232) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C*********************************************************************** C********************** METRIC CALCULATIONS **************************** C*********************************************************************** SUBROUTINE XYMETS(JDIM,KDIM,X,Y,XY,XIT,ETT,XYJ,WORK) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C COMMON/BASE/ 1 JMAX, KMAX, JM, KM, JBEGIN, JEND, 1 KBEGIN, KEND, JPLUS(999), JMINU (999), JLOW, JUP, 1 KLOW, KUP, PERIDC , NP, DT, CPUSEC, 1 FSMACH, ALPHA, GAMMA, GAMI, PI, 1 DIS2X, DIS2Y, DIS4X, DIS4Y, PHIDT, 1 THETAD , JACDT, RESID, IPRINT, NPCP, 1 JTAIL1, JTAIL2, NUMITE , ISTART, NSTEPS LOGICAL PERIDC COMMON/GRID/DYM,YMAX,XMIN,XMAX,THICK C DIMENSION XY(JDIM*KDIM*4),XYJ(JDIM*KDIM) DIMENSION XIT(JDIM*KDIM),ETT(JDIM*KDIM) DIMENSION X(JDIM*KDIM),Y(JDIM*KDIM) DIMENSION WORK(JDIM*KDIM*4) C C C COMPUTE XI DERIVATIVES OF X, Y C CALL XIDIF(JDIM,KDIM,X,Y,XY) C C COMPUTE ETA DERIVATIVES OF X, Y C CALL ETADIF(JDIM,KDIM,X,Y,XY) C C FORM METRICS AND JACOBIAN C CALL CALCME(JDIM,KDIM,XY,XIT,ETT,XYJ) C RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C***************************************************************** C***************************************************************** C ******** VECTORIZED IN J SCALAR PENTA SOLVER *************** C TWO F'S DONE AT A TIME *************** C***************************************************************** C***************************************************************** SUBROUTINE YPENT2(JDIM,KDIM,A,B,C,D,E,X,Y,F,KL,KU,JL,JU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LD,LD1,LD2,LDI DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM), * E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM,2) INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J13, J16, J19, X J22, J25 C C ! START FORWARD GENERATION PROCESS AND SWEEP C K = KL J7 = (JU - JL + 32) / 32 J10 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2233) CDOALL 10 J1=0,32*J7-32,32 DOUBLE PRECISION LD3(J10), LDI1(J10) INTEGER J11, J12 J12 = JL + J1 J11 = MIN0 (JL + J1 + 32 - 1, JU) - (JL + J1) + 1 LD3(1:$J11) = C(J12:$J11,K) LDI1(1:$J11) = 1. / LD3(1:$J11) F(J12:$J11,K,1) = F(J12:$J11,K,1) * LDI1(1:$J11) F(J12:$J11,K,2) = F(J12:$J11,K,2) * LDI1(1:$J11) X(J12:$J11,K) = D(J12:$J11,K) * LDI1(1:$J11) Y(J12:$J11,K) = E(J12:$J11,K) * LDI1(1:$J11) 10 CONTINUE C d call trace_exit (2234) K = KL+1 J13 = (JU - JL + 32) / 32 J10 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2235) CDOALL 20 J2=0,32*J13-32,32 DOUBLE PRECISION LD11(J10), LD4(J10), LDI2(J10) INTEGER J14, J15 J15 = JL + J2 J14 = MIN0 (JL + J2 + 32 - 1, JU) - (JL + J2) + 1 LD11(1:$J14) = B(J15:$J14,K) LD4(1:$J14) = C(J15:$J14,K) - LD11(1:$J14) * X(J15:$J14,K-1) LDI2(1:$J14) = 1. / LD4(1:$J14) F(J15:$J14,K,1) = (F(J15:$J14,K,1) - LD11(1:$J14) * F(J15:$J14,K X -1,1)) * LDI2(1:$J14) F(J15:$J14,K,2) = (F(J15:$J14,K,2) - LD11(1:$J14) * F(J15:$J14,K X -1,2)) * LDI2(1:$J14) X(J15:$J14,K) = (D(J15:$J14,K) - LD11(1:$J14) * Y(J15:$J14,K-1)) X * LDI2(1:$J14) Y(J15:$J14,K) = E(J15:$J14,K) * LDI2(1:$J14) 20 CONTINUE d call trace_exit (2236) J8 = JU - JL + 1 d call trace_enter (2237) CDOALL 11 J=JL,JL+J8-1 INTEGER K DOUBLE PRECISION LDI, LD2, LD1, LD C DO 1 K=KL+2,KU-2 LD2 = A(J,K) LD1 = B(J,K) - LD2 * X(J,K-2) LD = C(J,K) - (LD2 * Y(J,K-2) + LD1 * X(J,K-1)) LDI = 1. / LD F(J,K,1) = (F(J,K,1) - LD2 * F(J,K-2,1) - LD1 * F(J,K-1, X 1)) * LDI F(J,K,2) = (F(J,K,2) - LD2 * F(J,K-2,2) - LD1 * F(J,K-1, X 2)) * LDI X(J,K) = (D(J,K) - LD1 * Y(J,K-1)) * LDI Y(J,K) = E(J,K) * LDI 1 CONTINUE 11 CONTINUE C d call trace_exit (2238) K = KU-1 J16 = (JU - JL + 32) / 32 J8 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2239) CDOALL 12 J3=0,32*J16-32,32 DOUBLE PRECISION LD22(J8), LD13(J8), LD6(J8), LDI4(J8) INTEGER J17, J18 J18 = JL + J3 J17 = MIN0 (JL + J3 + 32 - 1, JU) - (JL + J3) + 1 LD22(1:$J17) = A(J18:$J17,K) LD13(1:$J17) = B(J18:$J17,K) - LD22(1:$J17) * X(J18:$J17,K-2) LD6(1:$J17) = C(J18:$J17,K) - (LD22(1:$J17) * Y(J18:$J17,K-2) + X LD13(1:$J17) * X(J18:$J17,K-1)) LDI4(1:$J17) = 1. / LD6(1:$J17) F(J18:$J17,K,1) = (F(J18:$J17,K,1) - LD22(1:$J17) * F(J18:$J17,K X -2,1) - LD13(1:$J17) * F(J18:$J17,K-1,1)) * LDI4(1:$J17) F(J18:$J17,K,2) = (F(J18:$J17,K,2) - LD22(1:$J17) * F(J18:$J17,K X -2,2) - LD13(1:$J17) * F(J18:$J17,K-1,2)) * LDI4(1:$J17) X(J18:$J17,K) = (D(J18:$J17,K) - LD13(1:$J17) * Y(J18:$J17,K-1)) X * LDI4(1:$J17) 12 CONTINUE C d call trace_exit (2240) K = KU J19 = (JU - JL + 32) / 32 J8 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2241) CDOALL 22 J4=0,32*J19-32,32 DOUBLE PRECISION LD23(J8), LD14(J8), LD7(J8), LDI5(J8) INTEGER J20, J21 J21 = JL + J4 J20 = MIN0 (JL + J4 + 32 - 1, JU) - (JL + J4) + 1 LD23(1:$J20) = A(J21:$J20,K) LD14(1:$J20) = B(J21:$J20,K) - LD23(1:$J20) * X(J21:$J20 X ,K-2) LD7(1:$J20) = C(J21:$J20,K) - (LD23(1:$J20) * Y(J21:$J20 X ,K-2) + LD14(1:$J20) * X(J21:$J20,K-1)) LDI5(1:$J20) = 1. / LD7(1:$J20) F(J21:$J20,K,1) = (F(J21:$J20,K,1) - LD23(1:$J20) * F( X J21:$J20,K-2,1) - LD14(1:$J20) * F(J21:$J20,K-1,1)) * LDI5(1:$ X J20) F(J21:$J20,K,2) = (F(J21:$J20,K,2) - LD23(1:$J20) * F( X J21:$J20,K-2,2) - LD14(1:$J20) * F(J21:$J20,K-1,2)) * LDI5(1:$ X J20) 22 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2242) J22 = (JU - JL + 32) / 32 d call trace_enter (2243) CDOALL 32 J5=0,32*J22-32,32 INTEGER J23, J24 J24 = JL + J5 J23 = MIN0 (JL + J5 + 32 - 1, JU) - (JL + J5) + 1 F(J24:$J23,KU,1) = F(J24:$J23,KU,1) F(J24:$J23,KU,2) = F(J24:$J23,KU,2) F(J24:$J23,KU-1,1) = F(J24:$J23,KU-1,1) - X(J24:$J23,KU-1) * F( X J24:$J23,KU,1) F(J24:$J23,KU-1,2) = F(J24:$J23,KU-1,2) - X(J24:$J23,KU-1) * F( X J24:$J23,KU,2) 32 CONTINUE C d call trace_exit (2244) d call trace_enter (2245) DO 2 K = 2,KU-KL KX = KU-K J25 = (JU - JL + 32) / 32 CDOALL 13 J6=0,32*J25-32,32 INTEGER J26, J27 J27 = JL + J6 J26 = MIN0 (JL + J6 + 32 - 1, JU) - (JL + J6) + 1 F(J27:$J26,KU-K,1) = F(J27:$J26,KU-K,1) - X(J27:$J26,KU-K) * X F(J27:$J26,KU-K+1,1) - Y(J27:$J26,KU-K) * F(J27:$J26,KU-K+2,1) F(J27:$J26,KU-K,2) = F(J27:$J26,KU-K,2) - X(J27:$J26,KU-K) * X F(J27:$J26,KU-K+1,2) - Y(J27:$J26,KU-K) * F(J27:$J26,KU-K+2,2) 13 CONTINUE 2 CONTINUE C d call trace_exit (2246) RETURN END C Cedar Fortran KAP (Fri, 14 Jun 1991) o3,r1,d3 16 Sep 1991 16:59:34 C***************************************************************** C***************************************************************** C ******** VECTORIZED IN J SCALAR PENTA SOLVER *************** C***************************************************************** C***************************************************************** SUBROUTINE YPENTA(JDIM,KDIM,A,B,C,D,E,X,Y,F,KL,KU,JL,JU) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LD,LD1,LD2,LDI DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM), * E(JDIM,KDIM),X(JDIM,KDIM),Y(JDIM,KDIM),F(JDIM,KDIM) INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J13, J16, J19, X J22, J25 C C ! START FORWARD GENERATION PROCESS AND SWEEP C K = KL J7 = (JU - JL + 32) / 32 J10 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2247) CDOALL 10 J1=0,32*J7-32,32 DOUBLE PRECISION LD3(J10), LDI1(J10) INTEGER J11, J12 J12 = JL + J1 J11 = MIN0 (JL + J1 + 32 - 1, JU) - (JL + J1) + 1 LD3(1:$J11) = C(J12:$J11,K) LDI1(1:$J11) = 1. / LD3(1:$J11) F(J12:$J11,K) = F(J12:$J11,K) * LDI1(1:$J11) X(J12:$J11,K) = D(J12:$J11,K) * LDI1(1:$J11) Y(J12:$J11,K) = E(J12:$J11,K) * LDI1(1:$J11) 10 CONTINUE C d call trace_exit (2248) K = KL+1 J13 = (JU - JL + 32) / 32 J10 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2249) CDOALL 20 J2=0,32*J13-32,32 DOUBLE PRECISION LD11(J10), LD4(J10), LDI2(J10) INTEGER J14, J15 J15 = JL + J2 J14 = MIN0 (JL + J2 + 32 - 1, JU) - (JL + J2) + 1 LD11(1:$J14) = B(J15:$J14,K) LD4(1:$J14) = C(J15:$J14,K) - LD11(1:$J14) * X(J15:$J14,K-1) LDI2(1:$J14) = 1. / LD4(1:$J14) F(J15:$J14,K) = (F(J15:$J14,K) - LD11(1:$J14) * F(J15:$J14,K-1)) X * LDI2(1:$J14) X(J15:$J14,K) = (D(J15:$J14,K) - LD11(1:$J14) * Y(J15:$J14,K-1)) X * LDI2(1:$J14) Y(J15:$J14,K) = E(J15:$J14,K) * LDI2(1:$J14) 20 CONTINUE d call trace_exit (2250) J8 = JU - JL + 1 d call trace_enter (2251) CDOALL 11 J=JL,JL+J8-1 INTEGER K DOUBLE PRECISION LDI, LD2, LD1, LD C DO 1 K=KL+2,KU-2 LD2 = A(J,K) LD1 = B(J,K) - LD2 * X(J,K-2) LD = C(J,K) - (LD2 * Y(J,K-2) + LD1 * X(J,K-1)) LDI = 1. / LD F(J,K) = (F(J,K) - LD2 * F(J,K-2) - LD1 * F(J,K-1)) * X LDI X(J,K) = (D(J,K) - LD1 * Y(J,K-1)) * LDI Y(J,K) = E(J,K) * LDI 1 CONTINUE 11 CONTINUE C d call trace_exit (2252) K = KU-1 J16 = (JU - JL + 32) / 32 J8 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2253) CDOALL 12 J3=0,32*J16-32,32 DOUBLE PRECISION LD22(J8), LD13(J8), LD6(J8), LDI4(J8) INTEGER J17, J18 J18 = JL + J3 J17 = MIN0 (JL + J3 + 32 - 1, JU) - (JL + J3) + 1 LD22(1:$J17) = A(J18:$J17,K) LD13(1:$J17) = B(J18:$J17,K) - LD22(1:$J17) * X(J18:$J17,K-2) LD6(1:$J17) = C(J18:$J17,K) - (LD22(1:$J17) * Y(J18:$J17,K-2) + X LD13(1:$J17) * X(J18:$J17,K-1)) LDI4(1:$J17) = 1. / LD6(1:$J17) F(J18:$J17,K) = (F(J18:$J17,K) - LD22(1:$J17) * F(J18:$J17,K-2) X - LD13(1:$J17) * F(J18:$J17,K-1)) * LDI4(1:$J17) X(J18:$J17,K) = (D(J18:$J17,K) - LD13(1:$J17) * Y(J18:$J17,K-1)) X * LDI4(1:$J17) 12 CONTINUE C d call trace_exit (2254) K = KU J19 = (JU - JL + 32) / 32 J8 = MAX0 (MIN0 (JL + 31, JU) - JL + 1, 1) d call trace_enter (2255) CDOALL 22 J4=0,32*J19-32,32 DOUBLE PRECISION LD23(J8), LD14(J8), LD7(J8), LDI5(J8) INTEGER J20, J21 J21 = JL + J4 J20 = MIN0 (JL + J4 + 32 - 1, JU) - (JL + J4) + 1 LD23(1:$J20) = A(J21:$J20,K) LD14(1:$J20) = B(J21:$J20,K) - LD23(1:$J20) * X(J21:$J20 X ,K-2) LD7(1:$J20) = C(J21:$J20,K) - (LD23(1:$J20) * Y(J21:$J20 X ,K-2) + LD14(1:$J20) * X(J21:$J20,K-1)) LDI5(1:$J20) = 1. / LD7(1:$J20) F(J21:$J20,K) = (F(J21:$J20,K) - LD23(1:$J20) * F(J21:$ X J20,K-2) - LD14(1:$J20) * F(J21:$J20,K-1)) * LDI5(1:$ X J20) 22 CONTINUE C C ! BACK SWEEP SOLUTION C d call trace_exit (2256) J22 = (JU - JL + 32) / 32 d call trace_enter (2257) CDOALL 32 J5=0,32*J22-32,32 INTEGER J23, J24 J24 = JL + J5 J23 = MIN0 (JL + J5 + 32 - 1, JU) - (JL + J5) + 1 F(J24:$J23,KU) = F(J24:$J23,KU) F(J24:$J23,KU-1) = F(J24:$J23,KU-1) - X(J24:$J23,KU-1) * F(J24:$ X J23,KU) 32 CONTINUE C d call trace_exit (2258) d call trace_enter (2259) DO 2 K = 2,KU-KL KX = KU-K J25 = (JU - JL + 32) / 32 CDOALL 13 J6=0,32*J25-32,32 INTEGER J26, J27 J27 = JL + J6 J26 = MIN0 (JL + J6 + 32 - 1, JU) - (JL + J6) + 1 F(J27:$J26,KU-K) = F(J27:$J26,KU-K) - X(J27:$J26,KU-K) * F( X J27:$J26,KU-K+1) - Y(J27:$J26,KU-K) * F(J27:$J26,KU-K+2) 13 CONTINUE 2 CONTINUE C d call trace_exit (2260) RETURN END