PROGRAM PULSE3D C********************************************************************** C (WRITES DATA FILE FOR RESTART VERSION SUITABLE FOR THE CRAY-2) C (WHOLE HEART, WITH VALVES) C (FIBER ACTIVATION: AN INTEGER ARRAY, LAFLAG, FLAGS ACTIVE LINKS, [A AND V]) c fixed bug in subroutine fluidup: c changed: c WR(I,J,K) = (WR(I,J,K) + PR(I,J,K)*PRFACT(K))/QRFACT(I,J,K) c to: c WR(I,J,K) = (WR(I,J,K) + PR(I,J,K)*PRFACT(K+1))/QRFACT(I,J,K) c (Combined upwind difference [hrtxp_61u6.f] and microtasking [hrtxp_61m3e.f]) C C OUTPUT FILES: C FORT.7 VELOCITY AND PRESSURE ON A COARSE GRID C FORT.8 VELOCITY AND PRESSURE ON A FINE GRID C FORT.9 VELOCITY AND LOCATION ON THE VALVE FLOWMETER WEBS C FORT.10 FIBER DATA C FORT.11 MARKER DATA C C INPUT PARAMETERS: C LCUBE=LENGTH OF SIDE OF CUBE (CM) C NU=KINEMATIC VISCOSITY (CM**2/SEC) C RHO=FLUID DENSITY (GM/CM**3) C (MU=NU*RHO=DYNAMIC VISCOSITY) C NG=NUMBER OF GRID POINTS IN EACH DIRECTION C TD=TIME STEP (SEC) C C PARAMETERS: C PLANES ARE NG X NG C BUT ARE STORED IN ARRAYS WITH DIMENSIONS C (0:NB,1:NG) WHERE NB=NG+2 C C A FIBER IS COMPOSED OF POINTS C A GROUP IS COMPOSED OF FIBERS HAVING THE SAME NUMBER OF POINTS C A BUNCH IS COMPOSED OF GROUPS C C NBUNCH = NUMBER OF BUNCHES IN THE ENTIRE STRUCTURE C NGROUPS(J)= NUMBER OF GROUPS IN BUNCH J, J=1,...,NBUNCH (NGROUPS(0)=0) C NFG(I) = NUMBER OF FIBERS IN GROUP I, I=NGROUPS(J-1)+1,...,NGROUPS(J) C NPF(I) = NUMBER OF POINTS IN A FIBER IN GROUP I C IMAX = SUM(J=1,NBUNCH):NGROUPS(J) C IMAX IS THE TOTAL NUMBER OF GROUPS IN THE ENTIRE STRUCTURE. C C NFSIZE = MAX(J=1,NBUNCH):SUM(I=ISTART(J),ISTOP(J)):NFG(I)*NPF(I) C ISTART(J) = NGROUPS(J-1)+1 C ISTOP (J) = NGROUPS(J) C NGROUPS(0)=0 C C ASSUMED VALUE MAX(I=1,IMAX):NFG(I) = 64 C ASSUMED VALUE MAX(I=1,IMAX):NPF(I) = 530 C ASSUMED VALUE MAX(I=1,IMAX):NFG(I)*NPF(I) = 33920 C C WARNING: THE ASSUMED VALUES GIVEN ABOVE C AND THE CONSTANT NFSIZE=606638 GIVEN C BELOW WERE DETERMINED EMPIRICALLY FOR A C PARTICULAR HEART ANATOMY. ANY CHANGE TO C THIS ANATOMY MAY REQUIRE ALTERATION OF C OF THESE VALUES. C A SIMILAR WARNING APPLIES TO NMSIZE. C c nsrcs = the number of sources in the heart, reasonably 5: c (1) superior vena cava c (2) inferior vena cava c (3) pulmonary veins (left atrium) c (4) pulmonary artery (normally thought of as a sink) c (5) aorta (normally thought of as a sink) c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(NFWBY2=NFSIZE/2) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) PARAMETER(NCONES=12) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) C SIZES OF CFFT3D WORK ARRAYS PARAMETER(NTABLE=100+2*3*NG) PARAMETER(NWORK=4*NG*NG+1 ) C C outflow valve leaflet beam parameters C PARAMETER(KGRPS=24) PARAMETER(NPF1=62,NPF2=80,NPF3=1,NPF4=20) PARAMETER(NPFINV=(NPF2-NPF1+1)+(NPF4-NPF3+1)) PARAMETER(NFGINV=256) PARAMETER(NLEAFS=6) PARAMETER(NPNFNL = NPFINV*NFGINV*NLEAFS) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) COMMON/BR/UR(0:NB,0:NB,0:NGM1) COMMON/BR/VR(0:NB,0:NB,0:NGM1) COMMON/BR/WR(0:NB,0:NB,0:NGM1) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/FFTFACT/ PRDENO( 0:NB , 0:NB ,0:NGM1) COMMON/FFTFACT/ QRFACT( 0:NB , 0:NB ,0:NGM1) COMMON/FFTFACT/ VRFACT( 0:NB ) COMMON/FFTFACT/ PRFACT( 0:NB ) REAL PRDENO,QRFACT COMPLEX VRFACT,PRFACT C WORK SPACE REQUIRED BY CFFT3D COMMON/FFTWORK/ FSCALE,INC1X,INC2X,INC3X COMMON/FFTWORK/ TABLE(NTABLE),WORK(NWORK,16) COMMON/KKGWORK/ TRIGX(2*NG),TRIGY(2*NG),TRIGZ(2*NG) COMMON/KKGWORK/ IFAX(19) , IFAY(19) , IFAZ(19) REAL FSCALE INTEGER INC1X,INC2X,INC3X REAL TABLE ,WORK REAL TRIGX ,TRIGY ,TRIGZ INTEGER IFAX , IFAY , IFAZ common/source/xsrc(3,nsrcspx) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ qsrc(nsrcs) common/source/ psrc(nsrcspx) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:511, nsrcspx) common/source/indxps(0:511, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1) complex qr COMMON/BSTAR/TD,TIME,USTAR,PSTAR,FSTAR,RSTAR COMMON/BSTAR/LCUBE,NU,MU,MASS,LENGTH REAL LCUBE,NU,MU,MASS,LENGTH COMMON/CNST/VSC,ACNST,BCNST COMMON/NPM/NP1(NG),NM1(NG) COMMON/WKSP/WKSP1(0:NB,1:NG),WKSP2(0:NB,1:NG) COMMON/BULGE/FLNG2,HZ,HSCALE,XSHIFT,YSHIFT,ZSHIFT COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) COMMON/XAR/ XF (3*NFSIZE), XFN (3,NFSIZE) COMMON/XAR/ STF0 ( NFSIZE) COMMON/XAR/ REST0 ( NFSIZE) COMMON/XAR/ FRC (3,NFSIZE) COMMON/XAR/ NEXTN ( NFSIZE),LAFLAG ( NFSIZE),ARFLAG ( NFSIZE) COMMON/XAR/ LAYER(IMAX*NFGMAX),NFIBER(0:NFGMAX,IMAX,0:NCONES) COMMON/XAR/ NGROUPS(NBUNCH) COMMON/XAR/ NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) COMMON/XAR/ MRAMP(IMAX),MFLAT LOGICAL ARFLAG COMMON/MAR/ XMK(3,NMSIZE) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) COMMON/HIST/ FIRSTM(1:NG,1:NG) COMMON/HIST/ NUMBEM(1:NG,1:NG) INTEGER FIRSTN INTEGER FIRSTM dimension nfskip (0:ncones) dimension npskip (0:ncones) dimension xatapex(3) C C IN THE FOLLOWING (WORKSPACE) ARRAYS, THE LAST DIMENSION, C SHOULD BE MAX(I=1,IMAX):NFG(I)*NPF(I). THIS IS THE LARGEST C NUMBER OF POINTS IN ANY ONE GROUP. C DIMENSION STF ( NFSIZE),REST( NFSIZE) DIMENSION BFIB(3,3,NPFGMX),CFIB(3,3,NPFGMX) DIMENSION XFIB(3, NPFGMX),XM (3,3,NPFGMX) LOGICAL UNSTBL( NFSIZE) C C outflow valve beam structures C DIMENSION XFINV(3,NPFINV,NFGINV,NLEAFS) DIMENSION FRCINV(3,NPFINV,NFGINV,NLEAFS) DIMENSION DX02( NPFINV,NFGINV,NLEAFS) DIMENSION DX0 ( NPFINV,NFGINV,NLEAFS) DIMENSION DDX0( NPFINV,NFGINV,NLEAFS) DIMENSION CSQ0( NPFINV,NFGINV,NLEAFS) DIMENSION ARANTI( NPFINV,NFGINV ) DIMENSION XF0(3 ) EQUIVALENCE( XFINV(1,1,1,1),XF( 1)) EQUIVALENCE(FRCINV(1,1,1,1),XF( 3*NPNFNL+1)) EQUIVALENCE( DX02( 1,1,1),XF( 6*NPNFNL+1)) EQUIVALENCE( DX0 ( 1,1,1),XF( 7*NPNFNL+1)) EQUIVALENCE( DDX0( 1,1,1),XF( 8*NPNFNL+1)) EQUIVALENCE( CSQ0( 1,1,1),XF( 9*NPNFNL+1)) EQUIVALENCE(ARANTI( 1,1 ),XF(10*NPNFNL+1)) EQUIVALENCE( XF0(1 ),XF(10*NPNFNL+1+NPFINV*NFGINV)) DIMENSION SBEND(NLEAFS) C C IN THE FOLLOWING (WORKSPACE) ARRAYS, THE LAST DIMENSION C SHOULD BE MAX(I=1,IMAX):NFG(I). THIS IS THE LARGEST C NUMBER OF FIBERS IN ANY ONE GROUP C DIMENSION CSAVE(3,3,NFGMAX),XCON(3,NFGMAX), D(3,NFGMAX) DIMENSION T1( NFGMAX),T2( NFGMAX) dimension nextw(NFSIZE) dimension timer(9),ttimer(9),rtimer(9) LOGICAL DENOVA INTEGER EXPNUM character*18 srcnam(nsrcs) data srcnam/ 1 'superior vena cava', 2 'inferior vena cava', 3 'pulmonary vein ', 4 'pulmonary artery ', 5 'aorta '/ C C DENOVA = .TRUE. EXPNUM = 800 NREFOLD = 1 NREF = 1 BETA = NREF**2 C C INPUT PARAMETERS -- C PI=4.*ATAN(1.) c c LCUBE is based on a domain of 64 meshwidths and a mitral ring radius c radius of 5.85 meshwidths which should represent the 10 cm c circumference of the human mitral ring c TD is based on dividing the period of the human heartbeat (0.8 sec) c into 256 steps c MFLAT = 16 DO 10 I=1,IMAX MRAMP(I) = 16 10 CONTINUE locref = 128 IF (DENOVA) THEN c NSTEP = 256*locref NSTEP = 64 LCUBE = float(64)*10.0/(2.0*pi*5.85) NU = 0.04 RHO = 1. TD = 0.8/float(256*locref) ESTOP = 2.**(-24) ITMAX = 128 KLOK = 0 call zero(u (0,1,0),nsize*ng) call zero(v (0,1,0),nsize*ng) call zero(w (0,1,0),nsize*ng) ELSE CALL RESTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,1,NREFOLD,NREF) NSTEP = 2048 END IF KLOK0 = KLOK KLOK1 = KLOK + 1 KLOKEND = KLOK + NSTEP MU = NU*RHO H = LCUBE/NG C C ESTABLISH CONVERSION FACTORS FROM CGS UNITS TO PROGRAM UNITS C MASS = RHO*H**3 LENGTH = H TIME = TD USTAR = LENGTH/TIME PSTAR = (MASS*LENGTH/TIME**2)/LENGTH**2 FSTAR = (MASS*LENGTH/TIME**2)/LENGTH**3 rstar = pstar*time/length**3 VSC = NU/((LENGTH**2)/TIME) H64 = LCUBE/64. TD64 = 0.000244140625 STFSCAL = ((TD/TD64)**2)/(RHO*(H/H64)**4) CALL INFIBER(XF , XFN , C STF0 , C REST0 ,laflag,arflag, C NGROUPS,NFG,NPF,KSTART,NFSTART, C NEST,LAYER,NFIBER,STFSCAL,katapex,xatapex) C NOTE THAT THE CONTENTS OF SEVERAL OF THE ARRAY ARGUMENTS ABOVE WILL BE C OVER-WRITTEN BY RESTART IF RESTART IS CALLED CALL INMARK(XMK,ncloud,nmarks,ncircs,nmarkt) CALL INHIST CALL INBEAM(XFN,FRC,XFINV,FRCINV,DX02,DX0,DDX0,CSQ0,ARANTI, C NFG,NPF,NPFINV,NFGINV,NLEAFS, C NPF1,NPF2,NPF3,NPF4,KGRPS) CALL INANCH(XF0,XFN(1,KSTART(25)),NFG(25),NPF(25),2) c c prsrvr(i) is the pressure in the reservoir connected to source i. c prsrvr in mm.Hg prsrvr(1) = 100.0 prsrvr(2) = 100.0 prsrvr(3) = 15.0 prsrvr(4) = 5.0 prsrvr(5) = 80.0 psrc (1) = 0.0 psrc (2) = 0.0 psrc (3) = 0.0 psrc (4) = 0.0 psrc (5) = 0.0 qsrc (1) = 0.0 qsrc (2) = 0.0 qsrc (3) = 0.0 qsrc (4) = 0.0 qsrc (5) = 0.0 c c prsrvr in dynes/cm**2 prconv = 1333.0 do 101 isrc=1,nsrcs prsrvr(isrc) = prsrvr(isrc)*prconv 101 continue resfact = 1.00 c resist(i) is the resistance between source i and its reservoir c resist in mm.Hg/(liters/min.) resist(1) = 40.0/resfact resist(2) = 40.0/resfact resist(3) = 2.0/resfact resist(4) = 2.0/resfact resist(5) = 1.5/resfact c c resist in (dynes/cm**2)/(cm**3/sec)=dynes*sec/cm**5 reconv = 1333.0*60.0/1000.0 do 102 isrc=1,nsrcs resist(isrc) = resist(isrc)*reconv 102 continue C C INITIALIZE THE FLUID SOLVER IF (DENOVA) THEN CALL INFLUIDU ELSE CALL RESTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,2,NREFOLD,NREF) CALL INFLUIDU xatapex(1) = xfn(1,katapex) xatapex(2) = xfn(2,katapex) xatapex(3) = xfn(3,katapex) END IF C INITIALIZE THE ACTIVATION VARIABLES -- IF (DENOVA) CALL INMUSCLE(TD,KLOK0,NREF) WRITE(6,*) NSTEP ,' = NSTEP' WRITE(6,*) LCUBE ,' CM = LCUBE' WRITE(6,*) NU ,' CM**2/SEC = NU' WRITE(6,*) RHO ,' GM/CM**3 = RHO' WRITE(6,*) NG,L2NG,' = NG L2NG' WRITE(6,*) TD ,' SEC = TD' WRITE(6,*) PI,' = PI' WRITE(6,*) H ,' CM = H ' WRITE(6,*) MU,' (GM/CM**3)*(CM**2/SEC) = MU' if (denova) then c print prsrvr in cgs units and then convert prsrvr to program units do 111 isrc=1,nsrcs write(6,131) prsrvr(isrc),srcnam(isrc) prsrvr(isrc) = prsrvr(isrc)/pstar 111 continue c print resist in cgs units and then convert resist to program units do 112 isrc=1,nsrcs write(6,132) resist(isrc),srcnam(isrc) resist(isrc) = resist(isrc)/rstar 112 continue else c print prsrvr in cgs units do 121 isrc=1,nsrcs write(6,131) prsrvr(isrc)*pstar,srcnam(isrc) 121 continue c print resist in cgs units do 122 isrc=1,nsrcs write(6,132) resist(isrc)*rstar,srcnam(isrc) 122 continue end if 131 format(1x,f9.1,' (gm*cm/sec**2)/cm**2 = ', c a18,' reservoir pressure') 132 format(1x,f7.1,' ((gm*cm/sec**2)/cm**2)/(cm**3/sec) = ', c a18,' resistance') WRITE(6,*)'CONVERSION FACTORS --' WRITE(6,*) MASS ,' = MASS' WRITE(6,*) LENGTH,' = LENGTH' WRITE(6,*) TIME ,' = TIME' WRITE(6,*) USTAR ,' = USTAR' WRITE(6,*) PSTAR ,' = PSTAR' WRITE(6,*) FSTAR ,' = FSTAR' WRITE(6,*) rstar ,' = rstar' WRITE(6,*) VSC ,' = VSC = KINEMATIC VISCOSITY IN PROGRAM UNITS' WRITE(6,*) STFSCAL,' = STFSCAL' WRITE(6,*) WRITE(6,*) ESTOP,' = ESTOP' WRITE(6,*) ITMAX,' = ITMAX' c c compute the index of the first source marker nmsrc = 1 do 1 icloud=1,6 nmsrc = nmsrc + nmarks(icloud) 1 continue nmtrl = nmarks(1) + 1 write(6,*)'the first source marker is number ',nmsrc write(6,*)'the first mitral marker is number ',nmtrl nrvmp = 1 do 141 icloud=1,15 nrvmp = nrvmp + nmarks(icloud) 141 continue npulm = 1 do 142 icloud=1,3 npulm = npulm + nmarks(icloud) 142 continue write(6,*)"the first rv mid marker is number ",nrvmp write(6,*)"the first pulmon marker is number ",npulm WRITE(6,*)'N,NM1,NP1=' DO 2 N=1,NG WRITE(6,*)N,NM1(N),NP1(N) 2 CONTINUE IF (DENOVA) THEN c write a restart file for klok=0 into file fort.2 CALL WRSTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C 0,2) STOP END IF C C WRITE OUT SOME OUTPUT FOR TIME T-- C klokout = mod(klok,4) T = 0. + FLOAT(KLOK0)*FLOAT(NREF)*TD CALL WRVPXF(EXPNUM,KLOK,T,XFN,NFG,NPF,NGROUPS,KSTART,NFSTART, C NEST,NFIBER,XMK,NCLOUD,NMARKS,NFSKIP,NPSKIP, C FRC) C C MAIN LOOP -- TLAST = SECOND() DO 5 KLOK=KLOK1,KLOKEND INTWR = 2048 INTOUT = 64 KLOKWR = MOD(KLOK, INTWR) KLOKOUT = MOD(KLOK, INTOUT) DO 4 NR=1,NREF T = T+TD c c compute the locations of the sources. c nmsrc is the index of the first source marker, c which is located in cloud 7. c the last argument is the number of source. call movesrc(xmk(1,nmsrc),nmarks(7),xsrc(1,1),nsrcs) c c compute the location of the left ventricular pressure tap. c nmtrl is the index of the first mitral ring marker, c which is located in cloud 2. c the last argument is the number of lv pressure taps c positioned by subroutine movesrc. call movesrc(xmk(1,nmtrl),nmarks(2),xsrc(1,nsrcs+1),1) write(6,45) nsrcs+1,(xsrc(ii,nsrcs+1),ii=1,3) do 44 is=2,5 fr = float(is-1)/5. xsrc(1,nsrcs+is) = xsrc(1,nsrcs+1)+fr*(xatapex(1)-xsrc(1,nsrcs+1)) xsrc(2,nsrcs+is) = xsrc(2,nsrcs+1)+fr*(xatapex(2)-xsrc(2,nsrcs+1)) xsrc(3,nsrcs+is) = xsrc(3,nsrcs+1)+fr*(xatapex(3)-xsrc(3,nsrcs+1)) write(6,45) nsrcs+is,(xsrc(ii,nsrcs+is),ii=1,3) 44 continue 45 format(' lv tap ',i2,3f6.2) 59 format(a8,i2,3f6.2) c compute the locations of the left ventricular outflow tract pressure taps c ktaps counts the number of additional pressure taps (max[ktaps]=ntapx) c nartc is the index of the first aortic ring marker, c which is in fact the very first marker (and in cloud 1). c The last argument is the number of ventricular outflow taps c positioned by subroutine movesrc. ktapsold = 5 ktaps = ktapsold + 3 nartc = 1 call movesrc(xmk(1,nartc),nmarks(1),xsrc(1,nsrcs+ktaps),1) do 46 is=ktapsold+1,ktaps-1 fr = float(is-ktapsold)/float(ktaps-ktapsold) xsrc(1,nsrcs+is) = xsrc(1,nsrcs+3) & +fr*(xsrc(1,nsrcs+ktaps)-xsrc(1,nsrcs+3)) xsrc(2,nsrcs+is) = xsrc(2,nsrcs+3) & +fr*(xsrc(2,nsrcs+ktaps)-xsrc(2,nsrcs+3)) xsrc(3,nsrcs+is) = xsrc(3,nsrcs+3) & +fr*(xsrc(3,nsrcs+ktaps)-xsrc(3,nsrcs+3)) c write(6,59) "lvo tap ",nsrcs+is,(xsrc(ii,nsrcs+is),ii=1,3) 46 continue c compute the locations of the aorta (artery) pressure taps c There is already a tap at the aortic sink which has index nsrcs. ktapsold = ktaps ktaps = ktapsold + 9 do 47 is=ktapsold+1,ktaps fr = float(is-ktapsold)/float(ktaps+1-ktapsold) xsrc(1,nsrcs+is) = xsrc(1,nsrcs+ktapsold) & +fr*(xsrc(1,nsrcs)-xsrc(1,nsrcs+ktapsold)) xsrc(2,nsrcs+is) = xsrc(2,nsrcs+ktapsold) & +fr*(xsrc(2,nsrcs)-xsrc(2,nsrcs+ktapsold)) xsrc(3,nsrcs+is) = xsrc(3,nsrcs+ktapsold) & +fr*(xsrc(3,nsrcs)-xsrc(3,nsrcs+ktapsold)) c write(6,59) "aor tap ",nsrcs+is,(xsrc(ii,nsrcs+is),ii=1,3) 47 continue c compute the locations of the right ventricular midplane pressure taps c nrvmp is the index of the first rv midplane marker, c which is located in cloud 16. c The last argument is the number of rv midplane pressure taps ktapsold = ktaps ktaps = ktapsold + 4 call movesrc(xmk(1,nrvmp),nmarks(16),xsrc(1,nsrcs+ktapsold+1),4) c compute the locations of the right ventricular outflow tract pressure taps c npulm is the index of the first pulmonic valve ring marker c which is in cloud 4. ktapsold = ktaps ktapsom1 = ktapsold - 1 ktaps = ktapsold + 3 call movesrc(xmk(1,npulm),nmarks(4),xsrc(1,nsrcs+ktaps),1) c the rest of the rv outflow taps are located along a line c which starts at the NEXT TO LAST rv midplane tap (i.e., ktapsold-1) c and ends at the tap in the center of the pulmonic ring (i.e., ktaps) do 48 is=ktapsold+1,ktaps-1 fr = float(is-ktapsold)/float(ktaps-ktapsold) xsrc(1,nsrcs+is) = xsrc(1,nsrcs+ktapsom1) & +fr*(xsrc(1,nsrcs+ktaps)-xsrc(1,nsrcs+ktapsom1)) xsrc(2,nsrcs+is) = xsrc(2,nsrcs+ktapsom1) & +fr*(xsrc(2,nsrcs+ktaps)-xsrc(2,nsrcs+ktapsom1)) xsrc(3,nsrcs+is) = xsrc(3,nsrcs+ktapsom1) & +fr*(xsrc(3,nsrcs+ktaps)-xsrc(3,nsrcs+ktapsom1)) c write(6,59) "rvo tap ",nsrcs+is,(xsrc(ii,nsrcs+is),ii=1,3) 48 continue c compute the locations of the pulmonary artery pressure taps c There is already a tap at the pulmonic sink which has index nsrcs-1. ktapsold = ktaps ktaps = ktapsold + 9 do 49 is=ktapsold+1,ktaps fr = float(is-ktapsold)/float(ktaps+1-ktapsold) xsrc(1,nsrcs+is) = xsrc(1,nsrcs+ktapsold) & +fr*(xsrc(1,nsrcs-1)-xsrc(1,nsrcs+ktapsold)) xsrc(2,nsrcs+is) = xsrc(2,nsrcs+ktapsold) & +fr*(xsrc(2,nsrcs-1)-xsrc(2,nsrcs+ktapsold)) xsrc(3,nsrcs+is) = xsrc(3,nsrcs+ktapsold) & +fr*(xsrc(3,nsrcs-1)-xsrc(3,nsrcs+ktapsold)) c write(6,59) "pul tap ",nsrcs+is,(xsrc(ii,nsrcs+is),ii=1,3) 49 continue C timer(1) = 0.0 C C UPDATE THE MUSCLE ACTIVATION FUNCTIONS CALL ACTIVATE(KLOK,1) C C EVALUATE FORCE DENSITY rtfi1 = rtc() tfib1 = second() C C EXPLICIT FORCE CALCULATION CALL FIBERX(XF0,STF0,REST0,FRC,LAFLAG,ARFLAG, C XFN,STF ,REST ,UNSTBL, C NFG,NPF,NGROUPS) C outflow valve leaflet beam forces SBEND(1) = 1.0 SBEND(2) = 1.0 SBEND(3) = 1.0 SBEND(4) = 1.0/6.0 SBEND(5) = 1.0/6.0 SBEND(6) = 1.0/6.0 DELTAU = 1.00 S1 = 0.0480000*STF(1) S2 = 0.0120000*STF(1) CALL FBEAMS(XFN,FRC,XFINV,FRCINV,DX02,DX0,DDX0,CSQ0,ARANTI, C NFG,NPF,NPFINV,NFGINV,NLEAFS, C NPF1,NPF2,NPF3,NPF4,KGRPS,DELTAU,S1,S2,SBEND) tfib2 = second() rtfi2 = rtc() timer(2) = tfib2 - tfib1 rtimer(2) = (rtfi2 - rtfi1)*4.167E-09 C C APPLY FORCE DENSITY TO FLUID c ttpus1 = TSECND() rtpus1 = rtc() tpush1 = second() CALL PUSHUP(XFN,FRC,NEXTN,NFG,NPF,NGROUPS,NBUNCH) tpush2 = second() rtpus2 = rtc() c ttpus2 = TSECND() timer(3) = tpush2 - tpush1 rtimer(3) = (rtpus2 - rtpus1)*4.167E-09 c ttimer(3) = ttpus2 - ttpus1 C C ADVANCE FLUID VELOCITY ONE TIME STEP -- rtfl1 = rtc() tflu1 = second() CALL FLUIDUP(nref) tflu2 = second() rtfl2 = rtc() timer(4) = tflu2 - tflu1 rtimer(4) = (rtfl2 - rtfl1)*4.167E-09 c c check the maximum values of the new velocity field rtchk1 = rtc() tchks1 = second() call chkspeed(speedm) tchks2 = second() rtchk2 = rtc() timer(5) = tchks2 - tchks1 rtimer(5) = (rtchk2 - rtchk1)*4.167E-09 C C MOVE THE BOUNDARY POINTS IN THE NEW VELOCITY FIELD c ttmov1 = TSECND() rtmov1 = rtc() tmove1 = second() CALL MOVE(XFN, NEXTN,NFG,NPF,NGROUPS,NBUNCH,katapex,xatapex) tmove2 = second() rtmov2 = rtc() c ttmov2 = TSECND() timer(6) = tmove2 - tmove1 rtimer(6) = (rtmov2 - rtmov1)*4.167E-09 c ttimer(6) = ttmov2 - ttmov1 C C MOVE THE FLUID MARKERS IN THE NEW VELOCITY FIELD rtmvm1 = rtc() tmovm1 = second() CALL MOVEMK(XMK,NEXTM,NMARKT) tmovm2 = second() rtmvm2 = rtc() timer(7) = tmovm2 - tmovm1 rtimer(7) = (rtmvm2 - rtmvm1)*4.167E-09 4 CONTINUE C C COMPUTE THE FLOWS THROUGH THE VALVE RINGS timer(8) = 0.0 if (klokout .eq. 0) then tflo1 = second() CALL FLOWS(xmk,nextw(1),nextw(nfwby2),nmarks,ncircs, c expnum,klok,klokout,nref) tflo2 = second() timer(8) = tflo2 - tflo1 c c compute the flow out of each of the sources do 50 isrc=1,5 call srcflow(isrc,xsrc(1,isrc),nextw(1),nextw(nfwby2), c klok,klokout,nref) 50 continue end if IF (KLOKWR .EQ. 0) THEN C WRITE RESTART DATA MULTWR = KLOK/INTWR IUNIT = 2 - MOD(MULTWR,2) CALL WRSTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,IUNIT) END IF twr1 = second() IF (KLOKOUT .EQ. 0) THEN C WRITE OUT SOME OUTPUT CALL WRVPXF(EXPNUM,KLOK,T,XFN,NFG,NPF,NGROUPS,KSTART,NFSTART, C NEST,NFIBER,XMK,NCLOUD,NMARKS,NFSKIP,NPSKIP, C FRC) END IF twr2 = second() timer(9) = twr2 - twr1 TTHIS = SECOND() TUSED = TTHIS - TLAST WRITE(6,999)TLAST,TTHIS,TUSED TLAST = TTHIS 999 FORMAT(' TIME IN = ',F8.2,' OUT = ',F8.2,' USED = ',F8.2) write(6,'(a17,1x,a14)') ' second',' rtc' write(6,'(a9,f8.4,1x,f14.9)')'fiberx ',timer(2),rtimer(2) write(6,'(a9,f8.4,1x,f14.9)')'pushup ',timer(3),rtimer(3) write(6,'(a9,f8.4,1x,f14.9)')'fluidup ',timer(4),rtimer(4) write(6,'(a9,f8.4,1x,f14.9)')'chkspeed ',timer(5),rtimer(5) write(6,'(a9,f8.4,1x,f14.9)')'move ',timer(6),rtimer(6) write(6,'(a9,f8.4,1x,f14.9)')'movemk ',timer(7),rtimer(7) write(6,'(a9,f8.4)') 'flows ',timer(8) write(6,'(a9,f8.4)') 'wrvpxf ',timer(9) WRITE(6,*)' KLOK: ',KLOK, ' ; TIME: ',T if (speedm .gt. 1.0) then write(6,*) write(6,*) "speedm=",speedm call exit(1) endif WRITE(6,*) 5 CONTINUE END SUBROUTINE CHEKFL(FL,NFG,NPF,NGROUPS,NBUNCH) PARAMETER(IMAX=63) DIMENSION FL(*),NFG(*),NPF(*),NGROUPS(*) DIMENSION FLMIN(IMAX),FLMAX(IMAX) ISTOP = 0 ISTART = 1 C DO OVER ALL BUNCHES DO 2000 NJ=1,NBUNCH ISTOP = ISTOP + NGROUPS(NJ) C DO OVER ALL GROUPS IN THE BUNCH K = 1 DO 1000 I=ISTART,ISTOP CALL FLIMIT(FL(K),NFG(I),NPF(I),FLMIN(I),FLMAX(I)) CDEBUGWRITE(0,100)I,FLMIN(I),FLMAX(I) 100 FORMAT(' I=',I2,' FLMIN=',E12.6,' FLMAX=',E12.6) K = K + NFG(I)*NPF(I) 1000 CONTINUE ISTART = ISTOP + 1 2000 CONTINUE FLMINMN = FLMIN(1) FLMAXMX = FLMAX(1) DO 3000 I=2,ISTOP IF (FLMIN(I) .LT. FLMINMN) FLMINMN = FLMIN(I) IF (FLMAX(I) .GT. FLMAXMX) FLMAXMX = FLMAX(I) 3000 CONTINUE WRITE(6,3001) FLMINMN/(27./512.),FLMAXMX/(27./512.) 3001 FORMAT(' MINIMUM RATIO FL/(27./512.) = ',E10.4, C /, ' MAXIMUM RATIO FL/(27./512.) = ',E10.4) RETURN END SUBROUTINE FLIMIT(FL,NFG,NPF,FLMIN,FLMAX) DIMENSION FL(NFG,NPF) FLMIN = FL(1,1) FLMAX = FL(1,1) DO 100 NP=1,NPF DO 100 NF=1,NFG IF (FL(NF,NP) .LT. FLMIN) THEN FLMIN = FL(NF,NP) ELSEIF (FL(NF,NP) .GT. FLMAX) THEN FLMAX = FL(NF,NP) END IF 100 CONTINUE RETURN END SUBROUTINE WRVPXF(EXPNUM, C KLOK,T,XFN,NFG,NPF,NGROUPS,KSTART,NFSTART, C NEST,NFIBER,XMK,NCLOUD,NMARKS,NFSKIP,NPSKIP, C FRC,STF,REST) save writ12,writ13 PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) DIMENSION NGROUPS(NBUNCH) DIMENSION NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) DIMENSION XFN(3,NFSIZE) DIMENSION CROSS(3,NPFGMX),VOLUME(IMAX) DIMENSION XMK(3,*) DIMENSION NMARKS(NCLOUD) DIMENSION NFIBER(0:NFGMAX,IMAX,0:NEST) DIMENSION NFSKIP(0:NEST) DIMENSION NPSKIP(0:NEST) DIMENSION FRC(3,NFSIZE) DIMENSION STF( NFSIZE) DIMENSION REST( NFSIZE) logical writ12,writ13 INTEGER EXPNUM character*14 filemk,filevp,filexf character*64 command do 1 icom=1,64 command(icom:icom) = ' ' 1 continue writ12 = .false. writ13 = .false. DO 10 NL=0,NEST NFSKIP(NL)= 2**(L2NG-5) NPSKIP(NL)= 2**(L2NG-5) 10 CONTINUE DO 20 NL=7,8 NFSKIP(NL)= 1 20 CONTINUE NPSKIP(7) = 1 NPSKIP( 8) = 1 NPSKIP( 9) = 1 NFSKIP( 9) = 1 NPSKIP(10) = 1 NFSKIP(10) = 1 NSKIP = 2**(L2NG-5) 101 format(a4,i2,a4,i1,a3) 102 format(a4,i2,a3,i2,a3) 103 format(a4,i2,a2,i3,a3) 104 format(a4,i2,a1,i4,a3) 105 format(a4,i2, i5,a3) 201 format(a3,i3,a4,i1,a3) 202 format(a3,i3,a3,i2,a3) 203 format(a3,i3,a2,i3,a3) 204 format(a3,i3,a1,i4,a3) 205 format(a3,i3, i5,a3) if (expnum .lt. 100) then if ( klok .lt. 10 ) then write(filemk,101) 'hrtx',expnum,'k000',klok,'.mk' write(filevp,101) 'hrtx',expnum,'k000',klok,'.vp' write(filexf,101) 'hrtx',expnum,'k000',klok,'.xf' elseif (klok .lt. 100 ) then write(filemk,102) 'hrtx',expnum,'k00' ,klok,'.mk' write(filevp,102) 'hrtx',expnum,'k00' ,klok,'.vp' write(filexf,102) 'hrtx',expnum,'k00' ,klok,'.xf' elseif (klok .lt. 1000 ) then write(filemk,103) 'hrtx',expnum,'k0' ,klok,'.mk' write(filevp,103) 'hrtx',expnum,'k0' ,klok,'.vp' write(filexf,103) 'hrtx',expnum,'k0' ,klok,'.xf' elseif (klok .lt. 10000) then write(filemk,104) 'hrtx',expnum,'k' ,klok,'.mk' write(filevp,104) 'hrtx',expnum,'k' ,klok,'.vp' write(filexf,104) 'hrtx',expnum,'k' ,klok,'.xf' else write(filemk,105) 'hrtx',expnum, klok,'.mk' write(filevp,105) 'hrtx',expnum, klok,'.vp' write(filexf,105) 'hrtx',expnum, klok,'.xf' end if else if ( klok .lt. 10 ) then write(filemk,201) 'hrt' ,expnum,'k000',klok,'.mk' write(filevp,201) 'hrt' ,expnum,'k000',klok,'.vp' write(filexf,201) 'hrt' ,expnum,'k000',klok,'.xf' elseif (klok .lt. 100 ) then write(filemk,202) 'hrt' ,expnum,'k00' ,klok,'.mk' write(filevp,202) 'hrt' ,expnum,'k00' ,klok,'.vp' write(filexf,202) 'hrt' ,expnum,'k00' ,klok,'.xf' elseif (klok .lt. 1000 ) then write(filemk,203) 'hrt' ,expnum,'k0' ,klok,'.mk' write(filevp,203) 'hrt' ,expnum,'k0' ,klok,'.vp' write(filexf,203) 'hrt' ,expnum,'k0' ,klok,'.xf' elseif (klok .lt. 10000) then write(filemk,204) 'hrt' ,expnum,'k' ,klok,'.mk' write(filevp,204) 'hrt' ,expnum,'k' ,klok,'.vp' write(filexf,204) 'hrt' ,expnum,'k' ,klok,'.xf' else write(filemk,205) 'hrt' ,expnum, klok,'.mk' write(filevp,205) 'hrt' ,expnum, klok,'.vp' write(filexf,205) 'hrt' ,expnum, klok,'.xf' end if end if IF (KLOK .GT. 0) THEN open(7,file=filevp,form='formatted') CALL VPOUT(filevp,KLOK,NSKIP,7) close(7) c write(command,'(a20,i2,a1,a14)') c c 'cfs -r5 store hrtin_',expnum,'/',filevp c write(6,'(a37)') command(1:37) c istat = iexec(command) c write(6,*) 'istat after ',command(1:37),' = ',istat END IF open(10,file=filexf,form='formatted') ISTOP = 0 ISTART = 1 WRITE(10,9910) NG,NEST if (writ12) CWRITE(12,9910) NG,NEST if (writ13) CWRITE(13,9910) NG,NEST DO 2001 NJ=1,NBUNCH ISTOP = ISTOP + NGROUPS(NJ) NFTOT = 0 DO 502 NL=0,NEST NGRLAY = 0 DO 401 I=ISTART,ISTOP IF (NFIBER(0,I,NL) .GT. 0) THEN NGRLAY = NGRLAY + 1 NFTOT = NFTOT + (1+(NFIBER(0,I,NL)-1)/NFSKIP(NL)) END IF 401 CONTINUE WRITE(10,9911) NGRLAY,NL if (writ12) CWRITE(12,9911) NGRLAY,NL if (writ13) CWRITE(13,9911) NGRLAY,NL DO 501 I=ISTART,ISTOP IF (NFIBER(0,I,NL) .GT. 0) THEN WRITE(10,9912) I,1+(NFIBER(0,I,NL)-1)/NFSKIP(NL), C 2+(NPF(I) -1)/NPSKIP(NL) if (writ12) C WRITE(12,9912) I,1+(NFIBER(0,I,NL)-1)/NFSKIP(NL), C 2+(NPF(I) -1)/NPSKIP(NL) if (writ13) C WRITE(13,9912) I,1+(NFIBER(0,I,NL)-1)/NFSKIP(NL), C 2+(NPF(I) -1)/NPSKIP(NL) END IF 501 CONTINUE 502 CONTINUE WRITE(10,9913) if (writ12) CWRITE(12,9913) if (writ13) CWRITE(13,9913) WRITE(10,9914) KLOK,filexf,T if (writ12) CWRITE(12,9914) KLOK,filexf,T if (writ13) CWRITE(13,9914) KLOK,filexf,T WRITE(10,9915) NFTOT if (writ12) CWRITE(12,9915) NFTOT if (writ13) CWRITE(13,9915) NFTOT 9910 FORMAT( 1X, I4,10X,' = NG' ,/, C 'C', I4,10X,' = MAXIMUM-LAYER-NUMBER, STARTING WITH 0') 9911 FORMAT('C',2I4, 6X,' = NUMBER-OF-GROUPS LAYER-NUMBER') 9912 FORMAT('C' 3I4, 2X,' = GRP NFG NPF' ) 9913 FORMAT('C*') 9914 FORMAT( I5,10X,' = KLOK',1x,a14 ,/, C 1X, F14.9 ,' = TIME' ) 9915 FORMAT( 1X, I4,10X,' = NUMBER OF FIBERS IN THIS DATA FILE') DO 1001 NL=0,NEST DO 1001 I=ISTART,ISTOP K = KSTART(I) CVOL CALL COMPVOL(XFN(1,K),CROSS,NFG(I),NPF(I),VOLUME(I),I) CALL XFOUT(KLOK,T,XFN(1,K),NFG(I),NPF(I),I,NL,NFIBER(0,I,NL), C NFSTART(I),NFSKIP(NL),NPSKIP(NL), C FRC(1,K),STF(K),REST(K),writ12,writ13) 1001 CONTINUE CVOL WRITE(6,*) CVOL WRITE(6,9920) VOLUME(ISTART),VOLUME(ISTOP), CVOL C VOLUME(ISTART)-VOLUME(ISTOP) ISTART = ISTOP + 1 2001 CONTINUE 9920 FORMAT(' OUTER VOL = ',F10.7, C ' INNER VOL = ',F10.7, C ' DIFFERENCE = ',F10.7) close(10) c write(command,'(a20,i2,a1,a14)') c c 'cfs -r5 store hrtin_',expnum,'/',filexf c write(6,'(a37)') command(1:37) c istat = iexec(command) c write(6,*) 'istat after ',command(1:37),' = ',istat open(11,file=filemk,form='formatted') WRITE(11,9930) NG,NCLOUD,(NMARKS(NC),NC,NC=1,NCLOUD) WRITE(11,9914) KLOK,filemk,T 9930 FORMAT('V3 ',' = FORMAT VERSION' , C /, I4 ,' = NG' , C /, I4 ,' = NUMBER OF CLOUDS ' , C 19(/, I4 ,' = NMARKS IN CLOUD ',I2) ) CALL XMKOUT(XMK,NCLOUD,NMARKS) close(11) write(command,'(a20,i2,a1,a14)') c 'cfs -r5 store hrtin_',expnum,'/',filemk c write(6,'(a37)') command(1:37) c istat = iexec(command) c write(6,*) 'istat after ',command(1:37),' = ',istat RETURN END SUBROUTINE COMPVOL(XF,CROSS,NFG,NPF,VOLUME,I) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) C C THE PURPOSE OF THE ROUTINE IS TO COMPUTE THE VOLUME IN LAYER I. C THE VOLUME IN THE LAYER IS NORMALIZED BY THE VOLUME OF THE DOMAIN. C DIMENSION XF(3,NFG,NPF),CROSS(3,NFG,NPF) DO 100 NP=1,NPF NPP1 = MOD(NP ,NPF)+1 NPM1 = MOD(NP +NPF-2,NPF)+1 NPM2 = MOD(NPM1+NPF-2,NPF)+1 DO 100 NF=1,NFG NFP1 = MOD(NF ,NFG)+1 NFM1 = MOD(NF +NFG-2,NFG)+1 IF (I .EQ. 1) THEN CROSS(1,NF,NP) = 1 + (XF(2,NFM1,NPP1)-XF(2,NF,NP)) * (XF(3,NFM1,NP )-XF(3,NF,NP)) 2 - (XF(3,NFM1,NPP1)-XF(3,NF,NP)) * (XF(2,NFM1,NP )-XF(2,NF,NP)) 3 + (XF(2,NFM1,NP )-XF(2,NF,NP)) * (XF(3,NF ,NPM1)-XF(3,NF,NP)) 4 - (XF(3,NFM1,NP )-XF(3,NF,NP)) * (XF(2,NF ,NPM1)-XF(2,NF,NP)) CROSS(2,NF,NP) = 1 + (XF(3,NFM1,NPP1)-XF(3,NF,NP)) * (XF(1,NFM1,NP )-XF(1,NF,NP)) 2 - (XF(1,NFM1,NPP1)-XF(1,NF,NP)) * (XF(3,NFM1,NP )-XF(3,NF,NP)) 3 + (XF(3,NFM1,NP )-XF(3,NF,NP)) * (XF(1,NF ,NPM1)-XF(1,NF,NP)) 4 - (XF(1,NFM1,NP )-XF(1,NF,NP)) * (XF(3,NF ,NPM1)-XF(3,NF,NP)) CROSS(3,NF,NP) = 1 + (XF(1,NFM1,NPP1)-XF(1,NF,NP)) * (XF(2,NFM1,NP )-XF(2,NF,NP)) 2 - (XF(2,NFM1,NPP1)-XF(2,NF,NP)) * (XF(1,NFM1,NP )-XF(1,NF,NP)) 3 + (XF(1,NFM1,NP )-XF(1,NF,NP)) * (XF(2,NF ,NPM1)-XF(2,NF,NP)) 4 - (XF(2,NFM1,NP )-XF(2,NF,NP)) * (XF(1,NF ,NPM1)-XF(1,NF,NP)) ELSE CROSS(1,NF,NP) = 1 + (XF(2,NFM1,NPM1)-XF(2,NF,NP)) * (XF(3,NFM1,NPM2)-XF(3,NF,NP)) 2 - (XF(3,NFM1,NPM1)-XF(3,NF,NP)) * (XF(2,NFM1,NPM2)-XF(2,NF,NP)) 3 + (XF(2,NFM1,NPM2)-XF(2,NF,NP)) * (XF(3,NF ,NPM1)-XF(3,NF,NP)) 4 - (XF(3,NFM1,NPM2)-XF(3,NF,NP)) * (XF(2,NF ,NPM1)-XF(2,NF,NP)) CROSS(2,NF,NP) = 1 + (XF(3,NFM1,NPM1)-XF(3,NF,NP)) * (XF(1,NFM1,NPM2)-XF(1,NF,NP)) 2 - (XF(1,NFM1,NPM1)-XF(1,NF,NP)) * (XF(3,NFM1,NPM2)-XF(3,NF,NP)) 3 + (XF(3,NFM1,NPM2)-XF(3,NF,NP)) * (XF(1,NF ,NPM1)-XF(1,NF,NP)) 4 - (XF(1,NFM1,NPM2)-XF(1,NF,NP)) * (XF(3,NF ,NPM1)-XF(3,NF,NP)) CROSS(3,NF,NP) = 1 + (XF(1,NFM1,NPM1)-XF(1,NF,NP)) * (XF(2,NFM1,NPM2)-XF(2,NF,NP)) 2 - (XF(2,NFM1,NPM1)-XF(2,NF,NP)) * (XF(1,NFM1,NPM2)-XF(1,NF,NP)) 3 + (XF(1,NFM1,NPM2)-XF(1,NF,NP)) * (XF(2,NF ,NPM1)-XF(2,NF,NP)) 4 - (XF(2,NFM1,NPM2)-XF(2,NF,NP)) * (XF(1,NF ,NPM1)-XF(1,NF,NP)) END IF 100 CONTINUE VOLUME = 0.0 DO 200 NP=1,NPF DO 200 NF=1,NFG VOLUME = VOLUME - CROSS(1,NF,NP)*XF(1,NF,NP) 2 - CROSS(2,NF,NP)*XF(2,NF,NP) 3 - CROSS(3,NF,NP)*XF(3,NF,NP) 200 CONTINUE VOLUME = VOLUME/(6.0*(NG**3)) RETURN END C*********************************************************************** SUBROUTINE WRSTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,IUNIT) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) PARAMETER(NCONES=12) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE) COMMON/XAR/ STF0 ( NFSIZE) COMMON/XAR/ REST0 ( NFSIZE) COMMON/XAR/ FRC (3,NFSIZE) COMMON/XAR/ NEXTN ( NFSIZE),LAFLAG ( NFSIZE),ARFLAG ( NFSIZE) COMMON/XAR/ LAYER(IMAX*NFGMAX),NFIBER(0:NFGMAX,IMAX,0:NCONES) COMMON/XAR/ NGROUPS(NBUNCH) COMMON/XAR/ NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) COMMON/XAR/ MRAMP(IMAX),MFLAT LOGICAL ARFLAG COMMON/MAR/ XMK(3,NMSIZE) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/CNST/VSC,ACNST,BCNST COMMON/NPM/NP1(NG),NM1(NG) common/source/xsrc(3,nsrcspx) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ qsrc(nsrcs) common/source/ psrc(nsrcspx) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:511, nsrcspx) common/source/indxps(0:511, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1) complex qr COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) COMMON/HIST/ FIRSTM(1:NG,1:NG) COMMON/HIST/ NUMBEM(1:NG,1:NG) INTEGER FIRSTN INTEGER FIRSTM REAL LCUBE,NU,MU,MASS,LENGTH C C INPUT PARAMETERS -- C REWIND(IUNIT) WRITE(IUNIT) NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX,MRAMP,MFLAT,KLOK WRITE(IUNIT) VSC,ACNST,BCNST WRITE(IUNIT) NP1,NM1 WRITE(IUNIT) NGROUPS,NFG,NPF write(iunit) xsrc,resist,prsrvr,qsrc,psrc WRITE(IUNIT) ACTPA ,ACTPV ,ACTPC ,ACTPR, C ALFAA ,ALFAV ,ALFAC ,ALFAR, C AA , AV , AC , AR, C TRAALF,TRVALF,TRCALF,TRRALF, C TRABET,TRVBET,TRCBET,TRRBET, C TFAALF,TFVALF,TFCALF,TFRALF, C TFABET,TFVBET,TFCBET,TFRBET, C NTA , NTV, C NDA , NDV, C NTAZ, NTVZ WRITE(IUNIT) FIRSTN,NUMBER,FIRSTM,NUMBEM WRITE(IUNIT) U WRITE(IUNIT) V WRITE(IUNIT) W C WRITE(IUNIT) XF WRITE(IUNIT) XFN WRITE(IUNIT) STF0 WRITE(IUNIT) REST0 WRITE(IUNIT) NEXTN C WRITE(IUNIT) FRC WRITE(IUNIT) XMK WRITE(IUNIT) NEXTM REWIND(IUNIT) RETURN END C*********************************************************************** SUBROUTINE RESTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,index,NREFOLD,NREF) save iunit PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) PARAMETER(NCONES=12) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE) COMMON/XAR/ STF0 ( NFSIZE) COMMON/XAR/ REST0 ( NFSIZE) COMMON/XAR/ FRC (3,NFSIZE) COMMON/XAR/ NEXTN ( NFSIZE),LAFLAG ( NFSIZE),ARFLAG ( NFSIZE) COMMON/XAR/ LAYER(IMAX*NFGMAX),NFIBER(0:NFGMAX,IMAX,0:NCONES) COMMON/XAR/ NGROUPS(NBUNCH) COMMON/XAR/ NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) COMMON/XAR/ MRAMP(IMAX),MFLAT LOGICAL ARFLAG COMMON/MAR/ XMK(3,NMSIZE) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/CNST/VSC,ACNST,BCNST COMMON/NPM/NP1(NG),NM1(NG) common/source/xsrc(3,nsrcspx) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ qsrc(nsrcs) common/source/ psrc(nsrcspx) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:511, nsrcspx) common/source/indxps(0:511, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1) complex qr COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) COMMON/HIST/ FIRSTM(1:NG,1:NG) COMMON/HIST/ NUMBEM(1:NG,1:NG) INTEGER FIRSTN INTEGER FIRSTM REAL LCUBE,NU,MU,MASS,LENGTH logical unit1,unit2 if ((index .ne. 1) .and. (index .ne. 2)) then write(6,*) 'BADINDEX' call exit(1) elseif (index .eq. 1) then inquire(file='fort.1',exist=unit1) inquire(file='fort.2',exist=unit2) if (unit1 .and. unit2) then write(6,*)' both fort.1 and fort.2 exist' rewind 1 read(1) NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX,MRAMP,MFLAT,KLOK1 rewind 2 read(2) NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX,MRAMP,MFLAT,KLOK2 if (klok1 .eq. klok2) then write(6,*)' fort.1 and fort.2 represent the same time-step' call exit(1) else if (klok1 .gt. klok2) then write(6,*)' fort.1 is the later of the two' iunit = 1 else write(6,*)' fort.2 is the later of the two' iunit = 2 end if elseif((.not.unit1) .and. (.not.unit2)) then write(6,*)' neither fort.1 nor fort.2 exists' call exit(1) elseif (unit1) then iunit = 1 elseif (unit2) then iunit = 2 end if C C INPUT PARAMETERS -- C REWIND(IUNIT) READ(IUNIT) NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX,MRAMP,MFLAT,KLOK READ(IUNIT) VSC,ACNST,BCNST READ(IUNIT) NP1,NM1 READ(IUNIT) NGROUPS,NFG,NPF TD = TD*FLOAT(NREFOLD)/FLOAT(NREF) return else read(iunit) xsrc,resist,prsrvr,qsrc,psrc READ(IUNIT) ACTPA ,ACTPV ,ACTPC ,ACTPR, C ALFAA ,ALFAV ,ALFAC ,ALFAR, C AA , AV , AC , AR, C TRAALF,TRVALF,TRCALF,TRRALF, C TRABET,TRVBET,TRCBET,TRRBET, C TFAALF,TFVALF,TFCALF,TFRALF, C TFABET,TFVBET,TFCBET,TFRBET, C NTA , NTV, C NDA , NDV, C NTAZ, NTVZ READ(IUNIT) FIRSTN, NUMBER, FIRSTM, NUMBEM READ(IUNIT) U READ(IUNIT) V READ(IUNIT) W C READ(IUNIT) XF READ(IUNIT) XFN READ(IUNIT) STF0 READ(IUNIT) REST0 READ(IUNIT) NEXTN C READ(IUNIT) FRC READ(IUNIT) XMK READ(IUNIT) NEXTM REWIND(IUNIT) IF (NREF .NE. NREFOLD) CALL REFINE(NREFOLD,NREF) end if RETURN END C*********************************************************************** SUBROUTINE REFINE(NREFOLD,NREF) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NCONES=12) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) COMMON/BSTAR/TD,TIME,USTAR,PSTAR,FSTAR,RSTAR COMMON/BSTAR/LCUBE,NU,MU,MASS,LENGTH REAL LCUBE,NU,MU,MASS,LENGTH COMMON/CNST/VSC,ACNST,BCNST COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE) COMMON/XAR/ STF0 ( NFSIZE) COMMON/XAR/ REST0 ( NFSIZE) COMMON/XAR/ FRC (3,NFSIZE) COMMON/XAR/ NEXTN ( NFSIZE),LAFLAG ( NFSIZE),ARFLAG ( NFSIZE) COMMON/XAR/ LAYER(IMAX*NFGMAX),NFIBER(0:NFGMAX,IMAX,0:NCONES) COMMON/XAR/ NGROUPS(NBUNCH) COMMON/XAR/ NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) COMMON/XAR/ MRAMP(IMAX),MFLAT LOGICAL ARFLAG common/source/xsrc(3,nsrcspx) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ qsrc(nsrcs) common/source/ psrc(nsrcspx) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:511, nsrcspx) common/source/indxps(0:511, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1) complex qr COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) FACT = FLOAT(NREFOLD)/FLOAT(NREF) FACTSQ = FACT**2 C TD = TD*FACT TIME = TD USTAR = LENGTH/TIME PSTAR = (MASS*LENGTH/TIME**2)/LENGTH**2 FSTAR = (MASS*LENGTH/TIME**2)/LENGTH**3 rstar = pstar*time/length**3 VSC = NU/((LENGTH**2)/TIME) ACNST = VSC BCNST = 1.+2.*VSC TRAALF = TRAALF*FACT TRVALF = TRVALF*FACT TRCALF = TRCALF*FACT TRRALF = TRRALF*FACT TRABET = TRABET*FACT TRVBET = TRVBET*FACT TRCBET = TRCBET*FACT TRRBET = TRRBET*FACT TFAALF = TFAALF*FACT TFVALF = TFVALF*FACT TFCALF = TFCALF*FACT TFRALF = TFRALF*FACT TFABET = TFABET*FACT TFVBET = TFVBET*FACT TFCBET = TFCBET*FACT TFRBET = TFRBET*FACT do 9 nsrc=1,nsrcspx c qsrc(nsrc) = qsrc(nsrc)*FACT psrc(nsrc) = psrc(nsrc)*FACTSQ 9 continue do 10 nsrc=1,nsrcs resist(nsrc) = resist(nsrc)*FACT prsrvr(nsrc) = prsrvr(nsrc)*FACTSQ qsrc (nsrc) = qsrc (nsrc)*FACT 10 continue DO 20 K=0,NGM1 DO 20 J=1,NG DO 20 I=1,NG U (I,J,K) = U (I,J,K)*FACT V (I,J,K) = V (I,J,K)*FACT W (I,J,K) = W (I,J,K)*FACT 20 CONTINUE DO 30 N=1,NFSIZE STF0 (N) = STF0 (N)*FACTSQ 30 CONTINUE NREFOLD = NREF RETURN END SUBROUTINE XMKOUT(XMK,NCLOUD,NMARKS) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) DIMENSION XMK(3,*) DIMENSION NMARKS(NCLOUD) KSTOP = 0 KSTRT = 1 DO 200 NC=1,NCLOUD KSTOP = KSTOP + NMARKS(NC) DO 100 K=KSTRT,KSTOP WRITE(11,11) XMK(1,K),XMK(2,K),XMK(3,K),K,NC 100 CONTINUE KSTRT = KSTRT + NMARKS(NC) 200 CONTINUE 11 FORMAT(3F7.2,1X,I4,1X,I2) RETURN END SUBROUTINE XFOUT(KLOK,T,XF,NFG,NPF,NGR,LAYER,NFIBER,NFSTART, C NFSKIP,NPSKIP,FRC,STF,REST,writ12,writ13) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NFGMAX=64) DIMENSION XF(3,NFG,NPF) DIMENSION NFIBER(0:NFGMAX) DIMENSION FRC(3,NFG,NPF) DIMENSION STF( NFG,NPF) DIMENSION REST( NFG,NPF) logical writ12,writ13 C C THIS INCARNATION OF XFOUT WRITES OUT ALL THE FIBERS IN ONE GROUP C WHICH ARE ALSO IN ONE LAYER. THE ARRAY NFIBER CONTAINS THE INDICES OF C THESE FIBERS NUMBERED FROM 1 TO THE MAXIMUM NUMBER IN THE HEART. C NFSTART CONTAINS THE INDEX OF THE FIRST FIBER IN THE ARRAY XF; IT IS C NECESSARY TO SUBTRACT NFSTART-1 FROM THE ENTRIES IN NFIBER TO OBTAIN C THE CORRECT REGISTRATION WITH XF. C NFSHIFT = NFSTART - 1 DO 200 NFIB=1,NFIBER(0),NFSKIP WRITE(10,10) NFIBER(NFIB),2+(NPF-1)/NPSKIP if (writ12) CWRITE(12,10) NFIBER(NFIB),2+(NPF-1)/NPSKIP if (writ13) CWRITE(13,10) NFIBER(NFIB),2+(NPF-1)/NPSKIP 10 FORMAT(1X,2I4,6X,' = FIBER POINTS') NF = NFIBER(NFIB) - NFSHIFT DO 100 NP=1,NPF,NPSKIP WRITE(10,210) XF(1,NF,NP),XF(2,NF,NP),XF(3,NF,NP),NP,NF,NGR,LAYER if (writ12) CWRITE(12,212) FRC(1,NF,NP),FRC(2,NF,NP),FRC(3,NF,NP), C NP,NF,NGR,LAYER if (writ13) CWRITE(13,213) STF( NF,NP),REST( NF,NP),NP,NF,NGR,LAYER 100 CONTINUE WRITE(10,210) XF(1,NF, 1),XF(2,NF, 1),XF(3,NF, 1), 1,NF,NGR,LAYER if (writ12) CWRITE(12,212) FRC(1,NF, 1),FRC(2,NF, 1),FRC(3,NF, 1), C 1,NF,NGR,LAYER if (writ13) CWRITE(13,213) STF( NF, 1),REST( NF, 1), 1,NF,NGR,LAYER 200 CONTINUE 210 FORMAT(3F7.2,1X,I4,I3,I3,I3) 212 FORMAT(3E16.8,1X,I4,I3,I3,I3) 213 FORMAT(2E16.8,1X,I4,I3,I3,I3) RETURN END SUBROUTINE VPOUT(filevp,KLOK,NSKIP,IUNIT) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) character*14 filevp WRITE(IUNIT,9)KLOK,filevp 9 FORMAT( I5,10X,' = KLOK',1x,a14, C /,8X,'U',16X,'V',16X,'W',16X,'P',8X,' I J K') 10 FORMAT(4E17.9,3I3) DO 100 K=0,NGM1,NSKIP DO 100 J=NSKIP,NG,NSKIP DO 100 I=NSKIP,NG,NSKIP WRITE(IUNIT,10) U(I,J,K),V(I,J,K),W(I,J,K), C P(I,J,K),I,J,K 100 CONTINUE RETURN END C********************************************************************** SUBROUTINE INMUSCLE(TD,KLOK0,NREF) c c initialize the activation variables c COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) C C TD = LENGTH OF ONE TIME-STEP IN SECONDS) C KLOK0 = LAST TIME-STEP BEFORE RESTART BEGINS (= 0 FOR INITIAL START) C ALFAA = FIRST STAGE OF TWO-STEP ACTIVATION (ATRIAL) C ALFAV = FIRST STAGE OF TWO-STEP ACTIVATION (VENTRICULAR) C ALFAC = FIRST STAGE OF TWO-STEP ACTIVATION (PAPILLARY) C AA = ATRIAL ACTIVATION C AV = VENTRICULAR ACTIVATION C AC = CHORD (PAPILLARY MUSCLE) ACTIVATION C ACTPA = MAXIMUM POSSIBLE VALUE OF A C ACTPV = MAXIMUM POSSIBLE VALUE OF V C ACTPC = MAXIMUM POSSIBLE VALUE OF C C TDAM ,TDVM = DURATIONS OF SYSTOLE (ATRIAL, VENTRICULAR) C TRAM ,TRVM = EXCITATION TIME-CONSTANT (ATRIAL, VENTRICULAR) C TFAM ,TFVM = RELAXATION TIME-CONSTANT (ATRIAL, VENTRICULAR) C TRAALF,TRVALF,TRCALF,TRRALF = FIRST-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TRABET,TRVBET,TRCBET,TRRBET = SECOND-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TFABET,TFVBET,TFCBET,TFRBET = ONE-STEP TIME FACTORS (RELAXATION [A,V,C,R]) C TDELAY = ATRIO-VENTRICULAR EXCITATION DELAY (SECONDS) C NDELAY = ATRIO-VENTRICULAR EXCITATION DELAY (TIME-STEPS) C NTA,NTAZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (ATRIAL) C NTV,NTVZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (VENTRICULAR AND PAPILLARY) C ACTPA = 0.25 ACTPV = ACTPA ACTPC = ACTPV ACTPR = ACTPV C TIMES IN SECONDS TDAM = 0.150 c TRAM = 0.025 asyst_exp0, non-convergence in klok=218, group=12 c TRAM = 0.032 asyst_exp1, non-convergence in klok=218, group=36 c TRAM = 0.047 asyst_exp2, non-convergence in klok=218, group=36 c introduce two-step relaxation c TRAM = 0.047 asyst_exp3, ok thru klok=224 c (non-convergence in klok=226, group=33) c TRAM = 0.025 asyst_exp4, non-convergence in klok=219, group=30 c TRAM = 0.032 asyst_exp5, non-convergence in klok=221, group=35 c TRAM = 0.062 asyst_exp6, non-convergence in klok=231, group=35) TRAM = 0.094 c TFAM = 0.060 TFAM = 0.047 TDVM = 0.300 c TRVM = 0.015 TRVM = 0.141 c TFVM = 0.060 TFVM = 0.047 c TDELAY = 0.100 TDELAY = 0.16875 C DIMENSIONLESS RATIO OF BET VALUES TO ALF VALUES FTIME = 1.0 C TIME FACTORS "DELTA T OVER TAU" TRAALF = TD/TRAM TRVALF = TD/TRVM TRCALF = TRVALF TRRALF = TRVALF TRABET = TRAALF*FTIME TRVBET = TRVALF*FTIME TRCBET = TRCALF*FTIME TRRBET = TRRALF*FTIME WRITE(6,9)'TR{A,V,C,R}ALF=',TRAALF,TRVALF,TRCALF,TRRALF WRITE(6,9)'TR{A,V,C,R}BET=',TRABET,TRVBET,TRCBET,TRRBET TFAALF = TD/TFAM TFVALF = TD/TFVM TFCALF = TFVALF TFRALF = TFVALF TFABET = TFAALF*FTIME TFVBET = TFVALF*FTIME TFCBET = TFCALF*FTIME TFRBET = TFRALF*FTIME WRITE(6,9)'TF{A,V,C,R}ALF=',TFAALF,TFVALF,TFCALF,TFRALF WRITE(6,9)'TF{A,V,C,R}BET=',TFABET,TFVBET,TFCBET,TFRBET 9 format(1x,a28,4f10.4) C TIMES IN TIME-STEPS NTA = 24577 NDA = TDAM/(TD*FLOAT(NREF)) NTAZ = NTA + NDA NDELAY = TDELAY/(TD*FLOAT(NREF)) NTV = NTA + NDELAY NDV = TDVM/(TD*FLOAT(NREF)) NTVZ = NTV + NDV WRITE(6,10)'NTA NTAZ=',NTA,NTAZ WRITE(6,10)'NTV NTVZ=',NTV,NTVZ 10 format(1x, a9,2i6 ) TRAALFR = TRAALF TRVALFR = TRVALF TRCALFR = TRCALF TRRALFR = TRRALF TRABETR = TRABET TRVBETR = TRVBET TRCBETR = TRCBET TRRBETR = TRRBET TFAALFR = TFAALF TFVALFR = TFVALF TFCALFR = TFCALF TFRALFR = TFRALF TFABETR = TFABET TFVBETR = TFVBET TFCBETR = TFCBET TFRBETR = TFRBET TRAALF = TRAALF*FLOAT(NREF) TRVALF = TRVALF*FLOAT(NREF) TRCALF = TRCALF*FLOAT(NREF) TRRALF = TRRALF*FLOAT(NREF) TRABET = TRABET*FLOAT(NREF) TRVBET = TRVBET*FLOAT(NREF) TRCBET = TRCBET*FLOAT(NREF) TRRBET = TRRBET*FLOAT(NREF) TFAALF = TFAALF*FLOAT(NREF) TFVALF = TFVALF*FLOAT(NREF) TFCALF = TFCALF*FLOAT(NREF) TFRALF = TFRALF*FLOAT(NREF) TFABET = TFABET*FLOAT(NREF) TFVBET = TFVBET*FLOAT(NREF) TFCBET = TFCBET*FLOAT(NREF) TFRBET = TFRBET*FLOAT(NREF) NREFR = NREF NREF = 1 ALFAA = 0.0 ALFAV = 0.0 ALFAC = 0.0 ALFAR = 0.0 AA = 0.0 AV = 0.0 AC = 0.0 AR = 0.0 DO 100 KLOK=1,KLOK0 DO 100 NR=1,NREF C IF (KLOK .EQ. NTAZ) THEN C ALFAA = AA C END IF KASYST = (NTA-KLOK)*(KLOK-NTAZ) IF (KASYST .GE. 0) THEN ALFAA = (ACTPA*TRAALF+ALFAA)/(TRAALF+1.) ELSE ALFAA = ALFAA /(TFAALF+1.) END IF IF (ALFAA .GE. AA) THEN AA = (ALFAA*TRABET+ AA)/(TRABET+1.) ELSE AA = (ALFAA*TFABET+ AA)/(TFABET+1.) END IF C IF (KASYST .GE. 0) THEN C ALFAA = (ACTPA*TRAALF+ALFAA)/(TRAALF+1.) C AA = (ALFAA*TRABET+ AA)/(TRABET+1.) C ELSE C ALFAA = ALFAA /(TFAALF+1.) C AA = (ALFAA*TFABET+ AA)/(TFABET+1.) C END IF C IF (KLOK .EQ. NTVZ) THEN C ALFAV = AV C ALFAC = AC C END IF KVSYST = (NTV-KLOK)*(KLOK-NTVZ) IF (KVSYST .GE. 0) THEN ALFAV = (ACTPV*TRVALF+ALFAV)/(TRVALF+1.) ALFAC = (ACTPC*TRCALF+ALFAC)/(TRCALF+1.) ALFAR = (ACTPR*TRRALF+ALFAR)/(TRRALF+1.) ELSE ALFAV = ALFAV /(TFVALF+1.) ALFAC = ALFAC /(TFCALF+1.) ALFAR = ALFAR /(TFRALF+1.) END IF IF (ALFAV .GE. AV) THEN AV = (ALFAV*TRVBET+ AV)/(TRVBET+1.) AC = (ALFAC*TRCBET+ AC)/(TRCBET+1.) AR = (ALFAR*TRRBET+ AR)/(TRRBET+1.) ELSE AV = (ALFAV*TFVBET+ AV)/(TFVBET+1.) AC = (ALFAC*TFCBET+ AC)/(TFCBET+1.) AR = (ALFAR*TFRBET+ AR)/(TFRBET+1.) END IF C IF (KVSYST .GE. 0) THEN C ALFAV = (ACTPV*TRVALF+ALFAV)/(TRVALF+1.) C AV = (ALFAV*TRVBET+ AV)/(TRVBET+1.) C ALFAC = (ACTPC*TRCALF+ALFAC)/(TRCALF+1.) C AC = (ALFAC*TRCBET+ AC)/(TRCBET+1.) C ALFAR = (ACTPR*TRRALF+ALFAR)/(TRRALF+1.) C AR = (ALFAR*TRRBET+ AR)/(TRRBET+1.) C ELSE C ALFAV = ALFAV /(TFVALF+1.) C AV = (ALFAV*TFVBET+ AV)/(TFVBET+1.) C ALFAC = ALFAC /(TFCALF+1.) C AC = (ALFAC*TFCBET+ AC)/(TFCBET+1.) C ALFAR = ALFAR /(TFRALF+1.) C AR = (ALFAR*TFRBET+ AR)/(TFRBET+1.) C END IF 100 CONTINUE write(6,200) AA,AV,AC 200 format(1x,'IN> AA, AV, AC=',3f10.6) TRAALF = TRAALFR TRVALF = TRVALFR TRCALF = TRCALFR TRRALF = TRRALFR TRABET = TRABETR TRVBET = TRVBETR TRCBET = TRCBETR TRRBET = TRRBETR TFAALF = TFAALFR TFVALF = TFVALFR TFCALF = TFCALFR TFRALF = TFRALFR TFABET = TFABETR TFVBET = TFVBETR TFCBET = TFCBETR TFRBET = TFRBETR NREF = NREFR RETURN END C********************************************************************** SUBROUTINE ACTIVATE(KLOK,index) c c update SFACTA,SFACTV,SFACTC,SFACTR c update RFACTA,RFACTV,RFACTC,RFACTR c COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) C C KLOK = WHAT TIME(-STEP) IS IT C ALFAA = FIRST STAGE OF TWO-STEP ACTIVATION (ATRIAL) C ALFAV = FIRST STAGE OF TWO-STEP ACTIVATION (VENTRICULAR) C ALFAC = FIRST STAGE OF TWO-STEP ACTIVATION (PAPILLARY) C AA = ATRIAL ACTIVATION C AV = VENTRICULAR ACTIVATION C AC = CHORD (PAPILLARY MUSCLE) ACTIVATION C ACTPA = MAXIMUM POSSIBLE VALUE OF A C ACTPV = MAXIMUM POSSIBLE VALUE OF V C ACTPC = MAXIMUM POSSIBLE VALUE OF C C TRAALF,TRVALF,TRCALF,TRRALF = FIRST-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TRABET,TRVBET,TRCBET,TRRBET = SECOND-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TFABET,TFVBET,TFCBET,TFRBET = ONE-STEP TIME FACTORS (RELAXATION [A,V,C,R]) C NTA,NTAZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (ATRIAL) C NTV,NTVZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (VENTRICULAR AND PAPILLARY) C C UPDATE THE ACTIVATION FUNCTIONS C C SET THE RESTING LENGTH ACTIVATION MULTIPLIERS ARFMA = 1.0 ARFMV = 1.5 ARFMC = 3.0 ARFMR = 4.0 C EVALUATE RFACTP WHICH IS USED TO PRE-STRESS C THE AORTIC VALVE LAYER FIBERS KLOKPR = 128 KPREST = ( 0-KLOK)*(KLOK-KLOKPR) IF (KPREST .GE. 0) THEN RFACTP(1) = 1. - 0.00*(FLOAT(KLOK)/FLOAT(KLOKPR)) RFACTP(2) = 1. - 0.00*(FLOAT(KLOK)/FLOAT(KLOKPR)) ELSE RFACTP(1) = 1. - 0.00 RFACTP(2) = 1. - 0.00 END IF C IF (KLOK .EQ. NTAZ) THEN C ALFAA = AA C END IF KASYST = (NTA-KLOK)*(KLOK-NTAZ) IF (KASYST .GE. 0) THEN ALFAA = (ACTPA*TRAALF+ALFAA)/(TRAALF+1.) ELSE ALFAA = ALFAA /(TFAALF+1.) END IF IF (ALFAA .GE. AA) THEN AA = (ALFAA*TRABET+ AA)/(TRABET+1.) ELSE AA = (ALFAA*TFABET+ AA)/(TFABET+1.) END IF C IF (KASYST .GE. 0) THEN C ALFAA = (ACTPA*TRAALF+ALFAA)/(TRAALF+1.) C AA = (ALFAA*TRABET+ AA)/(TRABET+1.) C ELSE C ALFAA = ALFAA /(TFAALF+1.) C AA = (ALFAA*TFABET+ AA)/(TFABET+1.) C END IF SFACTA = 1./(1.-AA) RFACTA = 1.*(1.-AA*ARFMA) C IF (KLOK .EQ. NTVZ) THEN C ALFAV = AV C ALFAC = AC C END IF KVSYST = (NTV-KLOK)*(KLOK-NTVZ) IF (KVSYST .GE. 0) THEN ALFAV = (ACTPV*TRVALF+ALFAV)/(TRVALF+1.) ALFAC = (ACTPC*TRCALF+ALFAC)/(TRCALF+1.) ALFAR = (ACTPR*TRRALF+ALFAR)/(TRRALF+1.) ELSE ALFAV = ALFAV /(TFVALF+1.) ALFAC = ALFAC /(TFCALF+1.) ALFAR = ALFAR /(TFRALF+1.) END IF IF (ALFAV .GE. AV) THEN AV = (ALFAV*TRVBET+ AV)/(TRVBET+1.) AC = (ALFAC*TRCBET+ AC)/(TRCBET+1.) AR = (ALFAR*TRRBET+ AR)/(TRRBET+1.) ELSE AV = (ALFAV*TFVBET+ AV)/(TFVBET+1.) AC = (ALFAC*TFCBET+ AC)/(TFCBET+1.) AR = (ALFAR*TFRBET+ AR)/(TFRBET+1.) END IF C IF (KVSYST .GE. 0) THEN C ALFAV = (ACTPV*TRVALF+ALFAV)/(TRVALF+1.) C AV = (ALFAV*TRVBET+ AV)/(TRVBET+1.) C ALFAC = (ACTPC*TRCALF+ALFAC)/(TRCALF+1.) C AC = (ALFAC*TRCBET+ AC)/(TRCBET+1.) C ALFAR = (ACTPR*TRRALF+ALFAR)/(TRRALF+1.) C AR = (ALFAR*TRRBET+ AR)/(TRRBET+1.) C ELSE C ALFAV = ALFAV /(TFVALF+1.) C AV = (ALFAV*TFVBET+ AV)/(TFVBET+1.) C ALFAC = ALFAC /(TFCALF+1.) C AC = (ALFAC*TFCBET+ AC)/(TFCBET+1.) C ALFAR = ALFAR /(TFRALF+1.) C AR = (ALFAR*TFRBET+ AR)/(TFRBET+1.) C END IF SFACTV = 1./(1.-AV)**13 RFACTV = 1.*(1.-AV*ARFMV) SFACTC = 1./(1.-AC)**13 RFACTC = 1.*(1.-AC*ARFMC) SFACTR = 1./(1.-AR)**13 RFACTR = 1.*(1.-AR*ARFMR) WRITE(6,6) ' A{A,V,C,R}=', AA, AV, AC, AR if (index .eq. 0) return WRITE(6,6) 'SFACT{A,V,C,R}=',SFACTA,SFACTV,SFACTC,SFACTR WRITE(6,6) 'RFACT{A,V,C,R}=',RFACTA,RFACTV,RFACTC,RFACTR C WRITE(6,6) 'RFACTP =',RFACTP C WRITE(6,6) 'RFACTP =',RFACTP 6 FORMAT(1x,a28,4f10.6) RETURN END C********************************************************************** SUBROUTINE MUSCLE( STF0, STF, C REST0,REST,laflag,arflag,NFG,NPF) c c update STF and REST for the flagged layers c COMMON/BACT/ACTPA ,ACTPV ,ACTPC ,ACTPR COMMON/BACT/ALFAA ,ALFAV ,ALFAC ,ALFAR COMMON/BACT/ AA , AV , AC , AR COMMON/BACT/TRAALF,TRVALF,TRCALF,TRRALF COMMON/BACT/TRABET,TRVBET,TRCBET,TRRBET COMMON/BACT/TFAALF,TFVALF,TFCALF,TFRALF COMMON/BACT/TFABET,TFVBET,TFCBET,TFRBET COMMON/BACT/ NTA , NTV COMMON/BACT/ NDA , NDV COMMON/BACT/ NTAZ, NTVZ COMMON/BACT/SFACTA,SFACTV,SFACTC,SFACTR COMMON/BACT/RFACTA,RFACTV,RFACTC,RFACTR,RFACTP(2) DIMENSION STF0(NFG,NPF), STF (NFG,NPF) DIMENSION REST0(NFG,NPF),REST (NFG,NPF) DIMENSION laflag(NFG,NPF) logical arflag(NFG,NPF) LOGICAL FLAGA,FLAGV,FLAGC,FLAGR,FLAGP1,FLAGP2 C C KLOK = WHAT TIME(-STEP) IS IT C ALFAA = FIRST STAGE OF TWO-STEP ACTIVATION (ATRIAL) C ALFAV = FIRST STAGE OF TWO-STEP ACTIVATION (VENTRICULAR) C ALFAC = FIRST STAGE OF TWO-STEP ACTIVATION (PAPILLARY) C AA = ATRIAL ACTIVATION C AV = VENTRICULAR ACTIVATION C AC = CHORD (PAPILLARY MUSCLE) ACTIVATION C ACTPA = MAXIMUM POSSIBLE VALUE OF A C ACTPV = MAXIMUM POSSIBLE VALUE OF V C ACTPC = MAXIMUM POSSIBLE VALUE OF C C TRAALF,TRVALF,TRCALF,TRRALF = FIRST-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TRABET,TRVBET,TRCBET,TRRBET = SECOND-STEP TIME FACTORS (EXCITATION [A,V,C,R]) C TFABET,TFVBET,TFCBET,TFRBET = ONE-STEP TIME FACTORS (RELAXATION [A,V,C,R]) C NTA,NTAZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (ATRIAL) C NTV,NTVZ = FIRST,LAST KLOK VALUES FOR SYSTOLE (VENTRICULAR AND PAPILLARY) C DO 10 NP=1,NPF DO 10 NF=1,NFG FLAGA = (laflag(NF,NP).eq.11) .or. (laflag(NF,NP).eq.12) STF( NF,NP) = CVMGT(SFACTA* STF0(NF,NP), STF0(NF,NP),FLAGA) REST(NF,NP) = CVMGT(RFACTA*REST0(NF,NP),REST0(NF,NP),FLAGA) 10 CONTINUE DO 13 NP=1,NPF DO 13 NF=1,NFG C FLAGV = (laflag(NF,NP).ge.0) .and. (laflag(NF,NP).le.6) FLAGV = (laflag(NF,NP).eq.0) .or. C ((laflag(NF,NP).ge.2) .and. (laflag(NF,NP).le.6)) STF( NF,NP) = CVMGT(SFACTV* STF0(NF,NP), STF(NF,NP),FLAGV) REST(NF,NP) = CVMGT(RFACTV*REST0(NF,NP),REST(NF,NP),FLAGV) 13 CONTINUE DO 14 NP=1,NPF DO 14 NF=1,NFG FLAGR = laflag(NF,NP).eq.1 STF( NF,NP) = CVMGT(SFACTR* STF0(NF,NP), STF(NF,NP),FLAGR) REST(NF,NP) = CVMGT(RFACTR*REST0(NF,NP),REST(NF,NP),FLAGR) 14 CONTINUE DO 16 NP=1,NPF DO 16 NF=1,NFG FLAGC = (laflag(NF,NP).eq.8) .or. (laflag(NF,NP).eq.9) STF( NF,NP) = CVMGT(SFACTC* STF0(NF,NP), STF(NF,NP),FLAGC) REST(NF,NP) = CVMGT(RFACTC*REST0(NF,NP),REST(NF,NP),FLAGC) 16 CONTINUE C PRE-STRESS THE AORTIC VALVE LAYER FIBERS C BY DECREASING THEIR RESTING LENGTHS DO 18 NP=1,NPF DO 18 NF=1,NFG FLAGP1 = ((laflag(NF,NP).eq.7).and. .not.arflag(NF,NP)) REST(NF,NP) = CVMGT(RFACTP(1)*REST0(NF,NP),REST(NF,NP),FLAGP1) 18 CONTINUE DO 19 NP=1,NPF DO 19 NF=1,NFG FLAGP2 = ((laflag(NF,NP).eq.7).and. arflag(NF,NP)) REST(NF,NP) = CVMGT(RFACTP(2)*REST0(NF,NP),REST(NF,NP),FLAGP2) 19 CONTINUE RETURN END C********************************************************************** SUBROUTINE INFIBER(XF , XFN , C STF0 , C REST0 ,laflag,arflag, C NGROUPS,NFG,NPF,KSTART,NFSTART, C NEST,LAYER,NFIBER,STFSCAL, c katapex,xatapex) C C THIS SUBROUTINE DRIVES SUBROUTINE FIBGEN WHICH READS DATA FOR A C FIBER-WRAPPED HEART FROM A PREVIOUSLY CONSTRUCTED DATA FILE. C C TERMINOLOGY: C C A FIBER IS COMPOSED OF POINTS C A GROUP IS COMPOSED OF FIBERS HAVING THE SAME NUMBER OF POINTS C A BUNCH IS COMPOSED OF GROUPS C C NBUNCH = NUMBER OF BUNCHES IN THE ENTIRE STRUCTURE C NGROUPS(J)= NUMBER OF GROUPS IN BUNCH J, J=1,...,NBUNCH (NGROUPS(0)=0) C NFG(I) = NUMBER OF FIBERS IN GROUP I, I=NGROUPS(J-1)+1,...,NGROUPS(J) C NPF(I) = NUMBER OF POINTS IN A FIBER IN GROUP I C IMAX = SUM(J=1,NBUNCH):NGROUPS(J) C IMAX IS THE TOTAL NUMBER OF GROUPS IN THE ENTIRE STRUCTURE. C C NFW = 512*(1+NFSIZE/512) C NFSIZE = MAX(J=1,NBUNCH):SUM(I=ISTART(J),ISTOP(J)):NFG(I)*NPF(I) C ISTART(J) = NGROUPS(J-1)+1 C ISTOP (J) = NGROUPS(J) C NGROUPS(0)=0 C C ASSUMED VALUE MAX(I=1,IMAX):NFG(I) = 64 C ASSUMED VALUE MAX(I=1,IMAX):NPF(I) = 530 C ASSUMED VALUE MAX(I=1,IMAX):NFG(I)*NPF(I) = 33920 C PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NCONES=12) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) C COMMON/BULGE/FLNG2,HZ,HSCALE,XSHIFT,YSHIFT,ZSHIFT C DIMENSION NGROUPS(NBUNCH) DIMENSION NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) DIMENSION XMIN(IMAX),XMAX(IMAX) DIMENSION YMIN(IMAX),YMAX(IMAX) DIMENSION ZMIN(IMAX),ZMAX(IMAX) dimension nfatzmin(imax),npatzmin(imax) dimension xatapex(3) DIMENSION XF (3,NFSIZE), XFN (3,NFSIZE) DIMENSION STF0 ( NFSIZE) DIMENSION REST0 ( NFSIZE) DIMENSION laflag(nfsize) LOGICAL arflag(nfsize) DIMENSION LAYER(*),NFIBER(0:NFGMAX,IMAX,0:NCONES) DIMENSION AXVALV(4),BYVALV(4),XCVALV(4),YCVALV(4) DIMENSION AXVEIN(4),BYVEIN(4),XCVEIN(4),YCVEIN(4) DIMENSION AXVENA(2),BYVENA(2),XCVENA(2),YCVENA(2) DIMENSION ZSRC(NSRCS) DIMENSION VLVNAM(4) DIMENSION SRCNAM(NSRCS) CHARACTER*9 VLVNAM CHARACTER*18 SRCNAM DATA VLVNAM /'AORTIC ','MITRAL ','TRICUSPID','PULMONIC'/ DATA SRCNAM /'SUPERIOR VENA CAVA', 2 'INFERIOR VENA CAVA', 3 'PULMONARY VEIN ', 4 'PULMONARY ARTERY ', 5 'AORTA '/ C C THERE IS ONE BUNCH CONTAINING MANY GROUPS. FIBERS IN A GROUP HAVE C NOTHING NECESSARILY IN COMMON OTHER THAN THE SAME NUMBER OF POINTS. C THE LAYER STRUCTURE OF THE HEART IS CONTAINED WITHIN AUXILIARY C DATA (FIBER NUMBER AND LAYER NUMBER) IN THE INPUT FILE; THE LAYER C STRUCTURE IS IMPORTANT ONLY FOR OUTPUT PURPOSES. C c open(4,file='/m/h2/bio0001/mcqueen/heart.in', c c form='formatted',access='sequential',status='old') c open(4,file='/m/s1/bio0002/mcqueen/heart.in', c c form='formatted',access='sequential',status='old') open(4,file='heart.in', c form='formatted',access='sequential',status='old') READ(4,*) NGINPUT READ(4,*) HEIGHT READ(4,*) THETAT READ(4,*) THETAL READ(4,*) THETAW READ(4,*) THETAR READ(4,*) NEST READ(4,*) HRINGS DO 10 IV=1,4 READ(4,*) AXVALV(IV),BYVALV(IV),XCVALV(IV),YCVALV(IV) 10 CONTINUE read(4,*) hveins do 11 iv=1,4 read(4,*) axvein(iv),byvein(iv),xcvein(iv),ycvein(iv) 11 continue do 12 iv=1,2 read(4,*) axvena(iv),byvena(iv),xcvena(iv),ycvena(iv) 12 continue do 13 is=1,nsrcs read(4,*) zsrc(is) 13 continue WRITE(6,21) NGINPUT,'NG ' WRITE(6,22) HEIGHT ,'HEIGHT ' WRITE(6,22) THETAT ,'THETAT (DEGREES) ' WRITE(6,22) THETAL ,'THETAL (DEGREES) ' WRITE(6,22) THETAW ,'THETAW (DEGREES) ' WRITE(6,22) THETAR ,'THETAR (DEGREES) ' WRITE(6,21) NEST ,'MAXIMUM LAYER NUMBER' 21 FORMAT(I4,3X,' = ',A20) 22 FORMAT(F7.2 ,' = ',A20) WRITE(6,23) HRINGS,' = HEIGHT OF VALVE RINGS ABOVE EQUATOR' 23 FORMAT(21X,F7.2,A38) DO 24 IV=1,4 WRITE(6,25) AXVALV(IV),BYVALV(IV),XCVALV(IV),YCVALV(IV),VLVNAM(IV) 24 CONTINUE 25 FORMAT(4F7.2,' = AX BY XC YC ',A9,' VALVE RING') WRITE(6,23) HVEINS,' = HEIGHT OF VEINS ABOVE VALVE RINGS ' DO 26 IV=1,4 WRITE(6,27) AXVEIN(IV),BYVEIN(IV),XCVEIN(IV),YCVEIN(IV), c 'PULMONARY VEIN NO.',IV 26 CONTINUE 27 FORMAT(4F7.2,' = AX BY XC YC ',A18,1X,I1) WRITE(6,27) AXVENA( 1),BYVENA( 1),XCVENA( 1),YCVENA( 1), c 'SUPERIOR VENA CAVA' WRITE(6,27) AXVENA( 2),BYVENA( 2),XCVENA( 2),YCVENA( 2), c 'INFERIOR VENA CAVA' DO 28 IS=1,5 WRITE(6,29) ZSRC(IS),SRCNAM(IS) 28 CONTINUE 29 FORMAT(21X,F7.2,' = HEIGHT OF ',A18,' SOURCE/SINK') call geometry(nginput,height, 2 thetat,thetal,thetaw,thetar, 3 alvinn,blvinn,xlvinn,ylvinn, 4 alvout,blvout,xlvout,ylvout, 5 aright,bright,xright,yright, 6 thetai ,thetaid) xlvinn = xlvinn + 32.0 ylvinn = ylvinn + 32.0 zlvinn = height xlvout = xlvout + 32.0 ylvout = ylvout + 32.0 zlvout = height xright = xright + 32.0 yright = yright + 32.0 zright = height write(6,*) " a b x y c z" write(6,*) alvinn,blvinn,xlvinn,ylvinn,zlvinn write(6,*) alvout,blvout,xlvout,ylvout,zlvout write(6,*) aright,bright,xright,yright,zright c scale up the lengths to accomodate the new computational grid c xfscal is passed as an argument to fibgen xfscal = float(ng)/float(nginput) height = height*xfscal hrings = hrings*xfscal hveins = hveins*xfscal do 32 iv=1,4 axvalv(iv) = axvalv(iv)*xfscal byvalv(iv) = byvalv(iv)*xfscal xcvalv(iv) = xcvalv(iv)*xfscal ycvalv(iv) = ycvalv(iv)*xfscal 32 continue do 33 iv=1,4 axvein(iv) = axvein(iv)*xfscal byvein(iv) = byvein(iv)*xfscal xcvein(iv) = xcvein(iv)*xfscal ycvein(iv) = ycvein(iv)*xfscal 33 continue do 34 iv=1,4 axvena(iv) = axvena(iv)*xfscal byvena(iv) = byvena(iv)*xfscal xcvena(iv) = xcvena(iv)*xfscal ycvena(iv) = ycvena(iv)*xfscal 34 continue do 35 is=1,nsrcs zsrc(is) = zsrc(is)*xfscal 35 continue alvinn = alvinn*xfscal blvinn = blvinn*xfscal xlvinn = xlvinn*xfscal ylvinn = ylvinn*xfscal zlvinn = height alvout = alvout*xfscal blvout = blvout*xfscal xlvout = xlvout*xfscal ylvout = ylvout*xfscal zlvout = height aright = aright*xfscal blvinn = blvinn*xfscal xright = xright*xfscal yright = yright*xfscal zright = height write(6,*) " a b x y c z" write(6,*) alvinn,blvinn,xlvinn,ylvinn,zlvinn write(6,*) alvout,blvout,xlvout,ylvout,zlvout write(6,*) aright,bright,xright,yright,zright write(6,*) thetai,thetaid ctemp = cos(thetai) stemp = sin(thetai) rtemp = 1./sqrt((ctemp/alvout)**2 + (stemp/blvout)**2) xtemp = xlvout + rtemp*cos(thetai) ytemp = ylvout + rtemp*sin(thetai) ztemp = zlvout write(6,"(3f7.2,2i4,a)") xtemp,ytemp,ztemp,0,0," SUBTENDED ANGLE" xtemp = xlvout ytemp = ylvout ztemp = zlvout write(6,"(3f7.2,2i4,a)") xtemp,ytemp,ztemp,0,0," SUBTENDED ANGLE" xtemp = xlvout + rtemp*cos(thetai) ytemp = ylvout - rtemp*sin(thetai) ztemp = zlvout write(6,"(3f7.2,2i4,a)") xtemp,ytemp,ztemp,0,0," SUBTENDED ANGLE" THETAT = THETAT*ATAN(1.0)/45.0 FLNG2 = FLOAT(NG/2) C HZ = HEIGHT + HRINGS HZ = HEIGHT HSCALE = 1.0/(HZ/2.0)**2 C C READ THE NUMBER OF GROUPS IN THE BUNCH C READ(4,*) NGROUPS(1) WRITE(6,31) NGROUPS(1) 30 CONTINUE 31 FORMAT(I4,8X,' = NUMBER OF GROUPS OF FIBERS') C C READ THE NUMBER OF FIBERS PER GROUP C AND THE NUMBER OF POINTS PER FIBER C NFIBS = 0 ISTART = 1 ISTOP = NGROUPS(1) NPTOTAL = 0 C DO OVER ALL GROUPS DO 100 I=ISTART,ISTOP READ(4,*) ICHECK,NFG(I),NPF(I) WRITE(6,40) I,NFG(I),NPF(I) 40 FORMAT(3I4,' = GROUP FIBERS/GROUP POINTS/FIBER') NFIBS = NFIBS + NFG(I) NPTOTAL = NPTOTAL + NFG(I)*NPF(I) 100 CONTINUE WRITE(6,201) NPTOTAL,NFSIZE 201 FORMAT(I12,' = TOTAL NUMBER OF POINTS (NFSIZE = ',I6,')') C C INITIALIZE THE ARRAYS LAYER AND NFIBER TO SUITABLE OUT-OF-RANGE VALUES C -1 IS SUITABLE FOR LAYER BECAUSE LAYERS ARE NUMBERED FROM 0 ON UP C 0 IS SUITABLE FOR NFIBER BECAUSE FIBERS ARE NUMBERED FROM 1 ON UP C DO 300 N=1,NFIBS LAYER(N) = -1 300 CONTINUE DO 400 NL=0,NCONES DO 400 NGR=1,IMAX NFIBER(0,NGR,NL) = 0 400 CONTINUE C C READ THE COORDINATES OF POINTS ON THE FIBERS C K = 1 NF = 1 DO 1000 I=ISTART,ISTOP CDEBUGWRITE(0,*)'CALLING FIBGEN WITH I=',I KSTART(I) = K NFSTART(I) = NF CALL FIBGEN(I,XF(1,K),STF0(K),REST0(K),laflag(k),arflag(k), C NFG(I),NPF(I), C LAYER,NFIBER,HEIGHT,HRINGS,hveins,STFSCAL, C XMIN(I),XMAX(I),YMIN(I),YMAX(I),ZMIN(I),ZMAX(I), C THETAT, c nfatzmin(i),npatzmin(i),xfscal, c alvout,blvout,xlvout,ylvout,zlvout, c aright,bright,xright,yright,zright, c thetai) K = K + NFG(I)*NPF(I) NF = NF + NFG(I) IF (K .GT. (NFSIZE+1)) THEN WRITE(6,*) 'KINFIBER' call exit(1) END IF KSTOP = K - 1 1000 CONTINUE C C Grade the stiffness of the outflow valve leaflets (groups 1 thru 24) C I = 1 K = KSTART(I) CALL GRADOUT(STF0(K),NFG(I),NPF(I)) KGRPS = 24 DO 1010 I=1,KGRPS K = KSTART(I) CALL GRADER(STF0(K),NFG(I),NPF(I)) 1010 CONTINUE I = 1 K = KSTART(I) CALL GRADOUT(STF0(K),NFG(I),NPF(I)) C C COMPUTE THE OVERALL MINIMUM AND MAXIMUM X-DIRECTION VALUES OF XF c and then compute the index of the point at the apex. c This point is used to position pressure taps in the lv along a c line from the apex to the center of the mitral ring. The apical c point is in the first bunch. C XLO = XMIN(ISTART) XHI = XMAX(ISTART) YLO = YMIN(ISTART) YHI = YMAX(ISTART) ZLO = ZMIN(ISTART) ZHI = ZMAX(ISTART) iatzmin = istart DO 1100 I=ISTART+1,ISTOP if (zmin(i) .lt. zlo) iatzmin = i IF (XMIN(I) .LT. XLO) XLO = XMIN(I) IF (XMAX(I) .GT. XHI) XHI = XMAX(I) IF (YMIN(I) .LT. YLO) YLO = YMIN(I) IF (YMAX(I) .GT. YHI) YHI = YMAX(I) IF (ZMIN(I) .LT. ZLO) ZLO = ZMIN(I) IF (ZMAX(I) .GT. ZHI) ZHI = ZMAX(I) 1100 CONTINUE c the index of the apex point is the sum of: c the number of points on all the fibers in all the groups up to, c but not including, the group which contains the apex point c plus c (in the group which contains the apex point) c the number of points on all the fibers up to, but not including, c the fiber which contains the apex point c plus c (in the group which contains the apex point, c on the fiber containing the apex point) c the number of points up to and including the apex point katapex = 0 do 1150 i=istart,iatzmin-1 katapex = katapex + nfg(i)*npf(i) 1150 continue katapex = katapex + (nfatzmin(iatzmin)-1)*npf(iatzmin) katapex = katapex + npatzmin(iatzmin) C C TRANSLATE THE HEART TO THE ORIGIN OF THE DOMAIN C XSHIFT = -(XLO + XHI)/2.0 YSHIFT = -(YLO + YHI)/2.0 ZSHIFT = -(ZLO + ZHI)/2.0 DO 1200 K=KSTART(ISTART),KSTOP XF(1,K) = XF(1,K) + XSHIFT XF(2,K) = XF(2,K) + YSHIFT XF(3,K) = XF(3,K) + ZSHIFT 1200 CONTINUE C C ROTATE THE HEART ABOUT THE ORIGIN SO THAT UNIT VECTORS FIXED TO C THE HEART AND ORIGINALLY POINTING IN THE X, Y AND Z DIRECTIONS C END UP POINTING IN THE Z, X AND Y DIRECTIONS, RESPECTIVELY. C SINCE THE (ORIGINAL) X EXTENT OF THE HEART IS LESS THAN THE C Y EXTENT OR THE Z EXTENT, THIS SHOULD TEND TO KEEP THE HEART C AWAY FROM THE COMPRESSIBLE FLOW FIELD OF THE BOUNDARY SINK. C DO 1210 K=KSTART(ISTART),KSTOP XTMP = XF(1,K) YTMP = XF(2,K) ZTMP = XF(3,K) XF(1,K) = YTMP XF(2,K) = ZTMP XF(3,K) = XTMP 1210 CONTINUE C C TRANSLATE THE HEART TO THE CENTER OF THE DOMAIN C THE ADDITIONAL 0.5 IN THE Z DIRECTION IS TO POSITION THE HEART C MIDWAY BETWEEN THE SINK ON PLANE 1 AND THE SINK ON PLANE NG C (I.E., THE PERIODIC IMAGE OF THE SINK ON PLANE 0) C XADD = FLOAT(NG)/2.0 YADD = FLOAT(NG)/2.0 ZADD = FLOAT(NG)/2.0 + 0.5 DO 1220 K=KSTART(ISTART),KSTOP XF(1,K) = XF(1,K) + XADD XF(2,K) = XF(2,K) + YADD XF(3,K) = XF(3,K) + ZADD 1220 CONTINUE C C Change REST0 on redundant fibers (groups 1,5,9,13,17,21) C DO 1226 I=1,21,4 K = KSTART(I) CALL NEWREST(REST0(K),NFG(I),NPF(I)) 1226 CONTINUE c c Make note of the position of the point at the apex. c xatapex(1) = xf(1,katapex) xatapex(2) = xf(2,katapex) xatapex(3) = xf(3,katapex) c c Initialize the array xfn c do 1300 k=1,nfsize xfn( 1,k) = xf( 1,k) xfn( 2,k) = xf( 2,k) xfn( 3,k) = xf( 3,k) 1300 continue close(4) RETURN END SUBROUTINE FIBGEN(I,XF,STF,REST,laflag,arflag,NFG,NPF, C LAYER,NFIBER,HEIGHT,HRINGS,hveins,STFSCAL, C XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX, C THETAT, c nfatzmin,npatzmin,xfscal, c alvout,blvout,xlvout,ylvout,zlvout, c aright,bright,xright,yright,zright, c thetai) PARAMETER(L2NG=7,NG=2**L2NG) PARAMETER(NFGMAX=64,NPFMAX=530) PARAMETER(IMAX=63) PARAMETER(NCONES=12) COMMON/BULGE/FLNG2,HZ,HSCALE,XSHIFT,YSHIFT,ZSHIFT DIMENSION XF (3,NFG,NPF) DIMENSION STF ( NFG,NPF) DIMENSION REST( NFG,NPF) DIMENSION LAFLAG(NFG,NPF) LOGICAL ARFLAG(NFG,NPF) DIMENSION COSPHI( NFGMAX),SINPHI( NFGMAX) DIMENSION LAYER(*),NFIBER(0:NFGMAX,IMAX,0:NCONES) DIMENSION LAYTMP(NFGMAX) dimension xfrot(3,npfmax),xcomp(3) dimension stfness(0:ncones) C C The line in the input file corresponding to the first point on each C fiber contains some auxiliary information in addition to the point C coordinates. Some of this information is redundant and/or not used, C placed in the file primarily to provide visual reference points C in the event that the data is perused by a human reader. Some of the C auxiliary information is used, however. C C NP1 = point number on the fiber (not used, first point always 1) C NF1 = fiber number (1..total number of fibers in the heart) C NGR = group number C NL0 = layer number in original layered structure C NF0 = fiber number in layer NL0 (not used) C NF1, NGR and NL0 are used for output purposes as follows: C LAYER(NF1) = layer number of fiber NF1 C NFIBER(K,NGR,NL0) = Kth fiber in group NGR in layer NL0 C NFIBER(0,NGR,NL0) = the number of fibers in group NGR in layer NL0; thus C fiber numbers (i.e., the values of NF1) are stored in C NFIBER(1,NGR,NL0)..NFIBER(NFNL0,NGR,NL0), C where NFNL0 = NFIBER(0,NGR,NL0) C DO 1 NF=1,NFG CDEBUGWRITE(0,*)'IN FIBGEN, DO 1 NF=',NF CDEBUGIF (I .EQ. 16) WRITE(0,*)'IN FIBGEN, NF=',NF READ(4,*) XF(1,NF, 1),XF(2,NF, 1),XF(3,NF, 1),NP1,NF1,NGR,NF0,NL0 xf(1,nf, 1) = xf(1,nf, 1)*xfscal xf(2,nf, 1) = xf(2,nf, 1)*xfscal xf(3,nf, 1) = xf(3,nf, 1)*xfscal LAYER(NF1) = NL0 LAYTMP(NF) = NL0 NFNL0 = NFIBER(0,NGR,NL0) + 1 NFIBER(NFNL0,NGR,NL0) = NF1 NFIBER( 0,NGR,NL0) = NFNL0 arflag(NF,1) = ((NL0 .eq. 7) .or. (NL0 .eq. 10)) .and. c (NF0 .gt. 768) DO 1 NP=2,NPF CDEBUGIF ((I .EQ. 16) .AND.(NF .EQ. 40)) WRITE(0,*)'IN FIBGEN, NP=',NP READ(4,*) XF(1,NF,NP),XF(2,NF,NP),XF(3,NF,NP) xf(1,nf,np) = xf(1,nf,np)*xfscal xf(2,nf,np) = xf(2,nf,np)*xfscal xf(3,nf,np) = xf(3,nf,np)*xfscal arflag(NF,NP) = arflag(NF,1) 1 CONTINUE c c stfness is a local mutiplier to increase (or decrease) stf c c stfness = 2.0 c stfness = 3.0 c fibers in the ventricles do 2 lay=0,6 stfness(lay) = stfscal*3.0/16.0 2 continue stfness(1) = stfness(1)*25.0 c fibers in the valve layers do 3 lay=7,10 stfness(lay) = stfscal*24.0 3 continue stfness(7) = stfness(7)*10.0*80.0/22.5 stfness(7) = stfness(7)*0.5 stfness(8) = stfness(8)*8.0 stfness(10) = stfness(10)/6.0 c fibers in the atria do 4 lay=11,12 stfness(lay) = stfscal*6.0*128.0 4 continue C C COMPUTE THE RESTING LENGTHS OF THE LINKS C DO 200 NP=1,NPF CDEBUGWRITE(0,*)'IN FIBGEN, DO 200 NP=',NP NP1 = MOD(NP,NPF) + 1 DO 200 NF=1,NFG REST(NF,NP) = SQRT( (XF(1,NF,NP) - XF(1,NF,NP1))**2 2 + (XF(2,NF,NP) - XF(2,NF,NP1))**2 3 + (XF(3,NF,NP) - XF(3,NF,NP1))**2) 200 CONTINUE C C COMPUTE THE STIFFNESS OF THE LINKS C DO 300 NF=1,NFG CDEBUGWRITE(0,*)'IN FIBGEN, DO 300 NF=',NF CDEBUGWRITE(0,*)'LAYTMP(NF) =',LAYTMP(NF) stiff = stfness(LAYTMP(NF)) DO 299 NP=1,NPF CDEBUGWRITE(0,*)'NP REST(NF,NP) =',NP,(XF(IC,NF,NP),IC=1,3) 299 CONTINUE DO 300 NP=1,NPF CDEBUGWRITE(0,*)'NP REST(NF,NP) =',NP,REST(NF,NP) STF (NF,NP) = stiff/REST(NF,NP) 300 CONTINUE c c Flag the links which are to be made contractile. In the left atrium c (layer 11) these are all links. (The left atrial appendage and the c pulmonary veins are to be made contractile.) In the c right atrium (layer 12) these are links which are simultaneously below c the height where the superior vena cava begins and above the height c where the inferior vena cava begins. c The flag value for contractile fibers is the layer number. The flag c value for non-contractile fibers is -1. Layers begin with layer 0. c This mechanism will make it possible to have contractile fibers in c other layers at a later stage of the project. c 'height' is the height of the equator above the apex (0.0) c 'hrings' is the height of the valve rings above the equator c 'hveins' is the height of the veins above the valve rings c hatrium = height + hrings + hveins do 305 np=1,npf do 305 nf=1,nfg laflag(nf,np) = -1 305 continue do 350 nf=1,nfg if (laytmp(nf) .eq. 1) then c c thetlv is the angle relative to the LV OUTER ellipse c thetrv is the angle relative to the RV ellipse c radilv is the actual distance from LV OUTER ellipse "center" c radirv is the actual distance from RV ellipse "center" c distlv is the separation from the LV OUTER ellipse c distrv is the separation from the RV ellipse c "separation" is measured along the line joining the point to the c ellipse in question. It would perhaps be better to measure separation c along the normal to the ellipse, but that is hard to compute and it c is expected that, in practice, the points will be VERY CLOSE to one c of the two ellipses except where the RV and the SEPTUM come together c and the small error should be tolerable. c do 306 np=1,npf c c provisionally, assign to each laflag(nf,np) the value of 2 c laflag = 1 is being redefined here to mean the RV FREE WALL only. laflag(nf,np) = 2 zscal = xf(3,nf,np)/height zscal = amin1(zscal,1.0) alv = alvout * zscal blv = blvout * zscal xlv = 64.0 + (xlvout-64.0)*zscal ylv = 64.0 + (ylvout-64.0)*zscal arv = aright * zscal brv = bright * zscal xrv = 64.0 + (xright-64.0)*zscal yrv = 64.0 + (yright-64.0)*zscal thetlv = atan2(xf(2,nf,np)-ylv, xf(1,nf,np)-xlv ) thetrv = atan2(xf(2,nf,np)-yrv, xf(1,nf,np)-xrv ) radilv = sqrt((xf(2,nf,np)-ylv)**2 + (xf(1,nf,np)-xlv)**2) radirv = sqrt((xf(2,nf,np)-yrv)**2 + (xf(1,nf,np)-xrv)**2) c = cos(thetai) s = sin(thetai) radius = 1./sqrt((c/alv)**2 + (s/blv)**2) xinter = xlv + c*radius yinthi = ylv + s*radius yintlo = ylv - s*radius c c check whether the point is within the angle subtended by the septum. if ((thetlv .ge. thetai) .or. (thetlv .le. -thetai)) then c c compute the location of the point at thetrv on the RV ellipse c = cos(thetrv) s = sin(thetrv) radius = 1./sqrt((c/arv)**2 + (s/brv)**2) distrv = abs(radius-radirv) c c compute the location of the point at thetlv on the LV ellipse c = cos(thetlv) s = sin(thetlv) radius = 1./sqrt((c/alv)**2 + (s/blv)**2) distlv = abs(radius-radilv) c c compare the distances to the two ellipses c note: radius is computed for the LV ellipse epsilon = 0.10 if (radilv .gt. radius+epsilon) then laflag(nf,np) = 1 elseif ((xf(2,nf,np) .gt. yinthi) .or. c (xf(2,nf,np) .lt. yintlo)) then laflag(nf,np) = 1 c elseif (distrv .lt. distlv ) then c laflag(nf,np) = 1 c else c write(6,"(3f7.2,2i4,a,f7.2)") c c xf(1,nf,np),xf(2,nf,np),xf(3,nf,np), c c nf,np," SEPTUM",yinter end if end if 306 continue c eliminate isolated points with laflag(nf,np)=1 do 307 np=1,npf npp1 = np + 1 if (np .eq. npf) npp1 = 1 npm1 = np - 1 if (np .eq. 1) npm1 = npf if (laflag(nf,np) .eq. 1) then if ((laflag(nf,npp1).ne.1).and.(laflag(nf,npm1).ne.1)) then c write(6,*) "isolated nf np = ",nf,np laflag(nf,np) = 2 end if end if 307 continue do 308 np=1,npf if (laflag(nf,np) .eq. 1) then write(6,"(3f7.2,2i4,a)")xf(1,nf,np),xf(2,nf,np),xf(3,nf,np), c nf,np," RV FREE WALL" else write(6,"(3f7.2,2i4,a)")xf(1,nf,np),xf(2,nf,np),xf(3,nf,np), c nf,np," SEPTUM" end if 308 continue elseif ((laytmp(nf) .ge. 0) .and. (laytmp(nf) .le. 6)) then do 310 np=1,npf laflag(nf,np) = laytmp(nf) 310 continue elseif (laytmp(nf) .eq. 7) then do 315 np=1,npf laflag(nf,np) = laytmp(nf) 315 continue elseif (laytmp(nf) .eq. 11) then do 320 np=1,npf laflag(nf,np) = laytmp(nf) 320 continue elseif (laytmp(nf) .eq. 12) then do 330 np=1,npf if ((hatrium-xf(3,nf,np))*(xf(3,nf,np)-height) .ge. 0.0) then laflag(nf,np) = laytmp(nf) end if 330 continue end if 350 continue C C THE AXIS OF SYMMETRY OF THE LV MAKES AN ANGLE THETAT WITH RESPECT TO C THE Z-AXIS. C APPLY TO EACH DATA POINT A RADIAL EXPANSION ABOUT THE INTERSECTION OF C THE AXIS OF SYMMETRY AND THE PLANE Z=Z-COORDINATE OF THE DATA POINT C WITH MAGNITUDE DETERMINED BY THE Z-COORDINATE OF THE DATA. C THE DATA IS CONTAINED BETWEEN Z=0.0 AND Z=HZ. C HZ, HSCALE AND FLNG2 HAVE BEEN SET IN SUBROUTINE INFIBER. C THE LAYTMP TEST INSURES THAT THE EXPANSION APPLIES ONLY TO CERTAIN C LAYERS. LAYERS > 6 ARE COMPOSED OF VALVES AND GREAT VESSELS. c the expansion mentioned above has been commented out with the symbol CN c which stands for 'Comment No expansion'. To re-enable this expansion, c simply replace CN by two blank spaces. c C CN TTH = TAN(THETAT) CN DO 400 NP=1,NPF CN DO 400 NF=1,NFG CN IF ((LAYTMP(NF) .LE. 6) .OR. (LAYTMP(NF) .EQ. 9)) THEN CN IF (XF(3,NF,NP) .LE. HZ) THEN CN XCEN = FLNG2 + TTH*XF(3,NF,NP) CN YCEN = FLNG2 CN SCALE = 1.0 + HSCALE* XF(3,NF,NP)*(HZ-XF(3,NF,NP)) CN XF(1,NF,NP) = XCEN + SCALE*(XF(1,NF,NP)-XCEN) CN XF(2,NF,NP) = YCEN + SCALE*(XF(2,NF,NP)-YCEN) CN END IF CN END IF CN400 CONTINUE c c Apply the expansion artfully to the mitral valve (layer 8) c c To understand the following, bear in mind that the fibers of the inflow c valves consist of three parts: a lower portion emanating from the apex c which is actually embedded in the lv inner layer, a middle portion which c connects the lv inner wall to an 'anchor point' at the bottom of the c valve leaflet, and an upper portion which goes from the 'anchor point' c up to the top of the leaflet. (Since the fiber returns to the apex one c might say there are two lower, two middle, and two upper portions.) c The middle portion constitutes the papillary muscle and the upper portion c constitutes the chord and the leaflet. Separating 'chord' from 'leaflet' c at this place in the overall design-to-test process would be difficult. c 'Leaflet' is the region where various fibers criss-cross to give the c appearance of a leaflet; 'chord' is the region where fibers run c side-by-side like the ropes descending from the parachute silk. c c Points on the inflow valve fibers are ordered: from lowest xf(3,,) c up to the valve leaflet and then back to lowest xf(3,,). To isolate c the chordae of any fiber: c (i) begin at point 1 and search with increasing index for a point c at or just above the equator; these are the two end points of c one chord of that fiber. c (ii) begin at point npf and search with decreasing index for a point c at or just above the equator; these are the two end points of c the other chord of that fiber. c For this reason the following loop, which is intended to isolate the c chordal portion of the mitral valve fibers has the fiber index as the c outer loop counter. c pi = 4.0*atan(1.0) conv = pi/180. thetat = 15.00*conv cthetat = cos(thetat) sthetat = sin(thetat) fpapil = 0.375 papz1 = 13.50*xfscal hzsrch = fpapil*hz hanchr = hz + hrings - papz1 CN restfac = 0.6 CN stffac = 4.0 c c Rotate each point of the fiber about the apex as if the lv axis of c symmetry were to be made vertical. This is an angle thetat about c the -y axis. The fiber kink is the point that is at or just above c hzsrch=fpapil*hz in the rotated coordinates. The end of the chord c is the point that is at or just above hz in the unrotated coordinates. c For this reason, it is not necessary to rotate the entire fiber before c searching, however since this rotation would vectorize in the absence c of if tests, rotate the entire fiber before searching. c Notice that some of the search loops employ the rotated fiber and other c search loops employ the unrotated fiber. Make modifications with care. c do 590 nf=1,nfg CDEBUGWRITE(0,*)'IN FIBGEN, DO 590 NF=',NF if (laytmp(nf) .eq. 8) then c rotate the entire fiber do 510 np=1,npf xfrot(1,np) = flng2 + cthetat*(xf(1,nf,np) - flng2) c - sthetat*(xf(3,nf,np) ) xfrot(2,np) = xf(2,nf,np) xfrot(3,np) = + sthetat*(xf(1,nf,np) - flng2) c + cthetat*(xf(3,nf,np) ) 510 continue c c search xfrot with increasing np for hzsrch do 520 np=1,npf if (xfrot(3,np) .gt. hzsrch) go to 521 nkink = np 520 continue c search xf with increasing np for hz 521 do 522 np=nkink+1,npf nanchr = np if (xf(3,nf,np) .ge. hanchr) go to 523 522 continue 523 continue c mark the links which are embedded in the wall or below the anchor point do 524 np=1,nanchr-1 laflag(nf,np) = laytmp(nf) laflag(nf,np) = 0 524 continue c apply the expansion to the lower part of the chord CN do 525 np=1,nkink CN XCEN = FLNG2 + TTH*XF(3,NF,NP) CN YCEN = FLNG2 CN SCALE = 1.0 + HSCALE* XF(3,NF,NP)*(HZ-XF(3,NF,NP)) CN XF(1,NF,NP) = XCEN + SCALE*(XF(1,NF,NP)-XCEN) CN XF(2,NF,NP) = YCEN + SCALE*(XF(2,NF,NP)-YCEN) CN525 continue c redistribute the points on the upper part of the chord CN xcomp(1) = xf(1,nf,nanchr)-xf(1,nf,nkink) CN xcomp(2) = xf(2,nf,nanchr)-xf(2,nf,nkink) CN xcomp(3) = xf(3,nf,nanchr)-xf(3,nf,nkink) CN do 526 np=nkink+1,nanchr CN frac = float(np-nkink)/float(nanchr-nkink) CN xf(1,nf,np) = xf(1,nf,nkink) + frac*xcomp(1) CN xf(2,nf,np) = xf(2,nf,nkink) + frac*xcomp(2) CN xf(3,nf,np) = xf(3,nf,nkink) + frac*xcomp(3) CN526 continue c recompute the resting lengths and stiffnesses of the chord links CN do 527 np=1,nkink-1 CN NP1 = NP + 1 CN REST(NF,NP) = restfac*SQRT( (XF(1,NF,NP) - XF(1,NF,NP1))**2 CN 2 + (XF(2,NF,NP) - XF(2,NF,NP1))**2 CN 3 + (XF(3,NF,NP) - XF(3,NF,NP1))**2) CN STF (NF,NP) = stffac*STFSCAL/REST(NF,NP) CN527 continue CN do 528 np=nkink,nanchr-1 CN NP1 = NP + 1 CN REST(NF,NP) = SQRT( (XF(1,NF,NP) - XF(1,NF,NP1))**2 CN 2 + (XF(2,NF,NP) - XF(2,NF,NP1))**2 CN 3 + (XF(3,NF,NP) - XF(3,NF,NP1))**2) CN STF (NF,NP) = stffac*STFSCAL/REST(NF,NP) CN REST(NF,NP) = 0.0 CN528 continue c c search xfrot with decreasing np for hzsrch do 530 np=npf,1,-1 if (xfrot(3,np) .gt. hzsrch) go to 531 nkink = np 530 continue c search xf with decreasing np for hz 531 do 532 np=nkink-1,1,-1 nanchr = np if (xf(3,nf,np) .ge. hanchr) go to 533 532 continue 533 continue c mark the links which are embedded in the wall or below the anchor point do 534 np=nanchr,npf laflag(nf,np) = laytmp(nf) laflag(nf,np) = 0 534 continue c The mitral valve is constructed by fans of chordae whose regions of c overlap are taken to be the leaflets. The chordae fan out from two c points onto an interpolating surface. On each mitral valve fiber, c define the papillary muscle to be the portions of the fiber which c lie below the fork onto the interpolating surface. The links on c these portions can be identified by the fact that their endpoints c are superposed on one another. Simultaneously search forward from c point 2 and backward from point npf, flagging links whose endpoints c are coincident. do 551 np=2,(npf+1)/2 distsq = (xf(1,nf,np)-xf(1,nf,npf-(np-2)))**2 c +(xf(2,nf,np)-xf(2,nf,npf-(np-2)))**2 c +(xf(3,nf,np)-xf(3,nf,npf-(np-2)))**2 if ((distsq .eq. 0.0 ) .and. c ((hrings+hz-xf(3,nf,np)) .ge. 16.0)) then if (laytmp(nf) .ne. 0) then laflag(nf,np-1) = laytmp(nf) laflag(nf,npf-(np-2)) = laytmp(nf) end if end if 551 continue do 552 np=2,npf-1 if ( laflag(nf,np ) .ne. laytmp(nf) ) then if ((laflag(nf,np-1) .eq. laytmp(nf)) .and. c (laflag(nf,np+1) .eq. laytmp(nf))) then laflag(nf,np ) = laytmp(nf) end if end if 552 continue c apply the expansion to the lower part of the chord CN do 535 np=npf,nkink,-1 CN XCEN = FLNG2 + TTH*XF(3,NF,NP) CN YCEN = FLNG2 CN SCALE = 1.0 + HSCALE* XF(3,NF,NP)*(HZ-XF(3,NF,NP)) CN XF(1,NF,NP) = XCEN + SCALE*(XF(1,NF,NP)-XCEN) CN XF(2,NF,NP) = YCEN + SCALE*(XF(2,NF,NP)-YCEN) CN535 continue c redistribute the points on the upper part of the chord CN xcomp(1) = xf(1,nf,nanchr)-xf(1,nf,nkink) CN xcomp(2) = xf(2,nf,nanchr)-xf(2,nf,nkink) CN xcomp(3) = xf(3,nf,nanchr)-xf(3,nf,nkink) CN do 536 np=nkink-1,nanchr,-1 CN frac = float(np-nkink)/float(nanchr-nkink) CN xf(1,nf,np) = xf(1,nf,nkink) + frac*xcomp(1) CN xf(2,nf,np) = xf(2,nf,nkink) + frac*xcomp(2) CN xf(3,nf,np) = xf(3,nf,nkink) + frac*xcomp(3) CN536 continue c recompute the resting lengths and stiffnesses of the chord links CN do 537 np=nanchr,nkink-1 CN NP1 = NP + 1 CN REST(NF,NP) = SQRT( (XF(1,NF,NP) - XF(1,NF,NP1))**2 CN 2 + (XF(2,NF,NP) - XF(2,NF,NP1))**2 CN 3 + (XF(3,NF,NP) - XF(3,NF,NP1))**2) CN STF (NF,NP) = stffac*STFSCAL/REST(NF,NP) CN REST(NF,NP) = 0.0 CN537 continue CN do 538 np=nkink,npf-1 CN NP1 = MOD(NP,NPF) + 1 CN REST(NF,NP) = restfac*SQRT( (XF(1,NF,NP) - XF(1,NF,NP1))**2 CN 2 + (XF(2,NF,NP) - XF(2,NF,NP1))**2 CN 3 + (XF(3,NF,NP) - XF(3,NF,NP1))**2) CN STF (NF,NP) = stffac*STFSCAL/REST(NF,NP) CN538 continue endif 590 continue C C COMPUTE THE MINIMUM AND MAXIMUM X-DIRECTION VALUES, XMIN AND XMAX. C COMPUTE THE MINIMUM AND MAXIMUM Y-DIRECTION VALUES, YMIN AND YMAX. C COMPUTE THE MINIMUM AND MAXIMUM Z-DIRECTION VALUES, ZMIN AND ZMAX. C THESE WILL BE USED SUBSEQUENTLY TO CENTER THE HEART IN THE DOMAIN. C XMIN = XF(1,1,1) XMAX = XF(1,1,1) YMIN = XF(2,1,1) YMAX = XF(2,1,1) ZMIN = XF(3,1,1) ZMAX = XF(3,1,1) nfatzmin = 1 npatzmin = 1 DO 600 NP=1,NPF CDEBUGWRITE(0,*)'IN FIBGEN, DO 600 NP=',NP DO 600 NF=1,NFG IF (XF(1,NF,NP) .LT. XMIN) THEN XMIN = XF(1,NF,NP) ELSEIF (XF(1,NF,NP) .GT. XMAX) THEN XMAX = XF(1,NF,NP) END IF IF (XF(2,NF,NP) .LT. YMIN) THEN YMIN = XF(2,NF,NP) ELSEIF (XF(2,NF,NP) .GT. YMAX) THEN YMAX = XF(2,NF,NP) END IF IF (XF(3,NF,NP) .LT. ZMIN) THEN ZMIN = XF(3,NF,NP) nfatzmin = nf npatzmin = np ELSEIF (XF(3,NF,NP) .GT. ZMAX) THEN ZMAX = XF(3,NF,NP) END IF 600 CONTINUE RETURN END C**************************************************************************** subroutine geometry(ng,height, 2 thetatd,thetald,thetawd,thetard, 3 alvinn ,blvinn ,xlvinn ,ylvinn , 4 alvout ,blvout ,xlvout ,ylvout , 5 aright ,bright ,xright ,yright , 6 thetai ,thetaid) c c constants c pi = 4.*atan(1.0) twopi = 2.*pi piby2 = pi/2. dphi = twopi/float(nsmax) radians = pi/180. FLNG = NG FLNG2 = NG*2 FLNG3 = NG*3 FLNGBY2 = NG/2 c c convert theta in degress to theta in radians c thetat = thetatd*radians thetal = thetald*radians thetaw = thetawd*radians thetar = thetard*radians alpha = tan(thetat) beta = height/cos(thetat) c innermost lv cone r = 1./tan(thetal/2.) alvinn = sqrt(1.+alpha**2)*beta*r/(r**2-alpha**2) blvinn = beta/sqrt(r**2-alpha**2) xlvinn = alpha*height*(r**2+1.)/(r**2-alpha**2) ylvinn = 0.0 arealv = pi*alvinn*blvinn c outermost lv cone r = 1./tan(thetaw+thetal/2.) alvout = sqrt(1.+alpha**2)*beta*r/(r**2-alpha**2) blvout = beta/sqrt(r**2-alpha**2) xlvout = alpha*height*(r**2+1.)/(r**2-alpha**2) ylvout = 0.0 ratio = blvout/alvout theta0 = thetat-thetal/2.-thetaw-thetar x0 = height*tan(theta0) call rvcoord(ratio,aright,xright,alvout,xlvout,omegar,omegal, c x0 ,arealv) bright = aright*ratio yright = 0.0 write(6,*) "lv inner a b x = ",alvinn,blvinn,xlvinn write(6,*) "lv outer a b x = ",alvout,blvout,xlvout write(6,*) "lv right a b x = ",aright,bright,xright call intrsct(alvout,blvout,xlvout,aright,bright,xright, c xinter,yinter) c c measure the location of (xinter,yinter) by the angle relative c to the "center" of the LV OUTER ellipse (xlvout,0.0). c thetai = atan2(yinter-0.0,xinter-xlvout) thetaid = thetai/radians write(6,*) "thetai = ",thetaid," degrees" return end subroutine rvcoord(ratio,aright,xright,a1,x1,omegar,omegal, c x0 ,arealv) parameter(pi=3.141592653589793, piby2=pi/2.) c c ratio = ratio of minor-axis to major-axis for the r.v. ellipse c and for the outermost ellipse of the l.v. cylinder c aright = major axis of r.v. ellipse c xright = center of r.v. ellipse c a1 = major axis of the outermost ellipse of the l.v. cylinder c x1 = center of the outermost ellipse of the l.v. cylinder c x0 = desired x-coordinate of the right tip of the r.v. ellipse, c where "the right tip" is furthest from the l.v. c arealv = the area enclosed by the innermost ellipse of the l.v. c c a2max ,a2min = max, min possible values of aright c xa2max,xa2min = values of xright corresponding to a2max,a2min c a2min = 0.5*(x1 - a1 - x0) xa2min = x0 + a2min a2max = 0.5*(x1 + a1 - x0) xa2max = x0 + a2max do 100 k=1,100 aright = (a2max+a2min)/2. xright = x0 + aright area = astar(ratio,aright,xright,a1,x1,xstar) write(6,*) "in rvcoord, k area arealv =",k,area,arealv if ((aright .eq. a2min) .or. (aright .eq. a2max)) then omegar = acos((xstar-xright)/aright) omegal = acos((xstar-x1 )/a1 ) return end if if (area .gt. arealv) then a2max = aright elseif(area .lt. arealv)then a2min = aright else omegar = acos((xstar-xright)/aright) omegal = acos((xstar-x1 )/a1 ) return end if 100 continue return end function astar(ratio,a2,x2,a1,x1,xstar) parameter(pi=3.141592653589793) xint(a1,a2,x1,x2) = 0.5*((x2+x1) - (a2**2 - a1**2)/(x2-x1)) yint(ratio,a1,xstar,x1) = ratio*sqrt(a1**2 - (xstar-x1)**2) piby2 = pi/2.0 xstar = xint(a1,a2,x1,x2) ystar = yint(ratio,a1,xstar,x1) area1 = piby2*ratio*a1**2 c + (xstar-x1)*ystar c + ratio*(a1**2)*asin((xstar-x1)/a1) area2 = piby2*ratio*a2**2 c - (xstar-x2)*ystar c - ratio*(a2**2)*asin((xstar-x2)/a2) astar = pi*ratio*a2**2 - (area1+area2) return end subroutine intrsct(alvout,blvout,xlvout,aright,bright,xright, c xinter,yinter) xmax = xlvout + alvout xmin = xright - aright 1 x = (xmax+xmin)/2. ylvout = blvout*sqrt(1.-((x-xlvout)/alvout)**2) yright = bright*sqrt(1.-((x-xright)/aright)**2) write(6,*) xmin,x,xmax,ylvout,yright if (ylvout .gt. yright) then xmax = x go to 1 else if (ylvout .lt. yright) then xmin = x go to 1 else go to 2 end if 2 xinter = x yinter = ylvout return end C**************************************************************************** SUBROUTINE GRADER(STF0,NFG,NPF) DIMENSION STF0(NFG,NPF) NPF1 = 61 NPF2 = 80 NPF3 = 1 NPF4 = 20 NPNUM = NPF1 NPDEN = NPF2-NPF1 DO 10 NP=NPF1,NPF2 WEIGHT = 1.0 - 0.9*(1.0 - (FLOAT(NP-NPNUM)/FLOAT(NPDEN))**4) DO 10 NF=1,NFG STF0(NF,NP) = STF0(NF,NP)*WEIGHT 10 CONTINUE NPNUM = NPF4 NPDEN = NPF4-NPF3 DO 20 NP=NPF3,NPF4 WEIGHT = 1.0 - 0.9*(1.0 - (FLOAT(NPNUM-NP)/FLOAT(NPDEN))**4) DO 20 NF=1,NFG STF0(NF,NP) = STF0(NF,NP)*WEIGHT 20 CONTINUE RETURN END SUBROUTINE NEWREST(REST0,NFG,NPF) DIMENSION REST0(NFG,NPF) DO 10 NP=1,NPF DO 10 NF=1,NFG/2 REST0(NF,NP) = REST0(NF,NP)/0.95 10 CONTINUE RETURN END SUBROUTINE GRADOUT(STF0,NFG,NPF) DIMENSION STF0(NFG,NPF) DO 10 NP=1,NPF WRITE(6,*) "STF0(1,",NP,")=",STF0(1,NP) 10 CONTINUE RETURN END C**************************************************************************** SUBROUTINE INANCH(XF0,XFN,NFG,NPF,NGRPS) DIMENSION XF0(3,*),XFN(3,*) DIMENSION NFG(NGRPS),NPF(NGRPS) K = 1 DO 10 NGRP=1,NGRPS CALL ANCHOR0(XF0(1,K),XFN(1,K),NFG(NGRP),NPF(NGRP)) K = K + NFG(NGRP)*NPF(NGRP) 10 CONTINUE RETURN END C**************************************************************************** SUBROUTINE ANCHOR0(XF0,XFN,NFG,NPF) DIMENSION XF0(3,NFG,NPF),XFN(3,NFG,NPF) DO 10 NP=1,NPF DO 10 NF=1,NFG XF0(1,NF,NP) = XFN(1,NF,NP) XF0(2,NF,NP) = XFN(2,NF,NP) XF0(3,NF,NP) = XFN(3,NF,NP) 10 CONTINUE RETURN END C**************************************************************************** SUBROUTINE ANCHORZ(XF0,XF,STF,REST,NFG,NPF,FRC) DIMENSION XF0(3,NFG,NPF),XF(3,NFG,NPF),FRC(3,NFG,NPF) DIMENSION STF(NFG,NPF),REST(NFG,NPF) DIMENSION F(3),D(3) LOGICAL RFLAG C CALL ZERO(FRC,3*NFG*NPF) SFACTR = 0.00000002 DO 1 K=1,NPF DO 1 N=1,NFG R= SQRT( (XF(1,N,K)-XF0(1,N,K))**2 C +(XF(2,N,K)-XF0(2,N,K))**2 C +(XF(3,N,K)-XF0(3,N,K))**2) IF (R .EQ. 0.0) GO TO 1 T = SFACTR*STF(N,K)*REST(N,K)*R D(1) = (XF0(1,N,K)-XF(1,N,K))/R D(2) = (XF0(2,N,K)-XF(2,N,K))/R D(3) = (XF0(3,N,K)-XF(3,N,K))/R F(1) = T*D(1) F(2) = T*D(2) F(3) = T*D(3) FRC(1,N,K )= FRC(1,N,K ) +F(1) FRC(2,N,K )= FRC(2,N,K ) +F(2) FRC(3,N,K )= FRC(3,N,K ) +F(3) 1 CONTINUE RETURN END C**************************************************************************** SUBROUTINE PUSHUP(XFN,FRC,NEXTN,NFG,NPF,NGROUPS,NBUNCH) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFSIZE=606638) PARAMETER(NPTMAX=256) PARAMETER(FLNG=NG) PARAMETER(NGP1=NG+1) PARAMETER(NGP2=NG+2) PARAMETER(NGBY4=NG/4) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/BR/UR(0:NB,0:NB,0:NGM1) COMMON/BR/VR(0:NB,0:NB,0:NGM1) COMMON/BR/WR(0:NB,0:NB,0:NGM1) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) INTEGER FIRSTN DIMENSION XFN(3,*),FRC(3,*),NEXTN(*) DIMENSION NFG(*),NPF(*),NGROUPS(*) DIMENSION XFN1OLD(NPTMAX),XFN2OLD(NPTMAX),XFN3OLD(NPTMAX) DIMENSION FORCE1(NPTMAX), FORCE2(NPTMAX), FORCE3(NPTMAX) DIMENSION IMOD(0:4), JMOD(0:4), KMOD(0:4) DIMENSION D1(NPTMAX,0:3), D2(NPTMAX,0:3), D3(NPTMAX,0:3) DIMENSION D12(NPTMAX,0:3,0:3) DIMENSION DELTA(0:64,NPTMAX),INDEX(0:64) DIMENSION ULIN(-15:16*NGP2),VLIN(-15:16*NGP2),WLIN(-15:16*NGP2) DIMENSION INDUVW( 16*NG ) DIMENSION FLKZP1(NPTMAX) dimension lindx(NPTMAX) dimension timer(10) C C THIS ROUTINE APPLIES THE FIBER FORCES TO THE FLUID. C C FIBER POINT L INFLUENCES FLUID POINT (I,J,K) IFF C ALL OF THE FOLLOWING CONDITIONS ARE SATISFIED: C C ABS(I-XFN(1,L)) < 2. C ABS(J-XFN(2,L)) < 2. C ABS(K-XFN(3,L)) < 2. C C USING FORTRAN THESE VALUES OF (I,J,K) ARE FOUND AS FOLLOWS: C C IZ= XFN(1,L) -1. C JZ= XFN(2,L) -1. C KZ= XFN(3,L) -1. C C THEN C C I = IZ,IZ+3 C J = JZ,JZ+3 C K = KZ,KZ+3 C C IF THESE CONDITIONS ARE SATISFIED, THE COEFFICIENT THAT LINKS C FIBER POINT L TO THE FLUID POINT (I,J,K) IS: C C D = DEL(I-XFN(1,L)) * DEL(J-XFN(2,L)) * DEL(K-XFN(3,L)) C C WHERE C DEL(R) = (1. + COS((PI/2.)*R)) / 4. C C ALGORITHM IS AS FOLLOWS: C C LET L=1,2, ... ,NP C FOR EACH L, LET (I,J,K) COVER RANGE OF VALUES DEFINED ABOVE. C C FOR EACH VALUE OF (L,I,J,K) C C U(I,J,K) = U(I,J,K) + D * FRC(1,L) C V(I,J,K) = V(I,J,K) + D * FRC(2,L) C W(I,J,K) = W(I,J,K) + D * FRC(3,L) C C WHERE D WAS DEFINED ABOVE. C MODNG(K)=MOD(K+NG,NG) DEL(R) = (1. + COS((PI/2.)*R))/4. PI = 4. * ATAN(1.) CMIC$ DO ALL CMIC$1SHARED(U, V, W, UR, VR, WR) CMIC$2PRIVATE(K, KM1, KP1, J, JM1, JP1, I) DO 105 K=0,NGM1 KM1 = K - 1 KP1 = K + 1 IF (K .EQ. 0) THEN KM1 = NGM1 ELSEIF (K .EQ. NGM1) THEN KP1 = 0 END IF DO 105 J=1,NG JM1 = J - 1 JP1 = J + 1 IF (J .EQ. 1) THEN JM1 = NG ELSEIF (J .EQ. NG) THEN JP1 = 1 END IF U( 0,J,K) = U(NG,J,K) U(NGP1,J,K) = U( 1,J,K) V( 0,J,K) = V(NG,J,K) V(NGP1,J,K) = V( 1,J,K) W( 0,J,K) = W(NG,J,K) W(NGP1,J,K) = W( 1,J,K) DO 105 I=1,NG UR(I,J,K) = U(I,J,K) C -U(I,J,K)*CVMGP(U(I ,J ,K )-U(I-1,J ,K ), C U(I+1,J ,K )-U(I ,J ,K ), C U(I ,J ,K )) C -V(I,J,K)*CVMGP(U(I ,J ,K )-U(I ,JM1,K ), C U(I ,JP1,K )-U(I ,J ,K ), C V(I ,J ,K )) C -W(I,J,K)*CVMGP(U(I ,J ,K )-U(I ,J ,KM1), C U(I ,J ,KP1)-U(I ,J ,K ), C W(I ,J ,K )) VR(I,J,K) = V(I,J,K) C -U(I,J,K)*CVMGP(V(I ,J ,K )-V(I-1,J ,K ), C V(I+1,J ,K )-V(I ,J ,K ), C U(I ,J ,K )) C -V(I,J,K)*CVMGP(V(I ,J ,K )-V(I ,JM1,K ), C V(I ,JP1,K )-V(I ,J ,K ), C V(I ,J ,K )) C -W(I,J,K)*CVMGP(V(I ,J ,K )-V(I ,J ,KM1), C V(I ,J ,KP1)-V(I ,J ,K ), C W(I ,J ,K )) WR(I,J,K) = W(I,J,K) C -U(I,J,K)*CVMGP(W(I ,J ,K )-W(I-1,J ,K ), C W(I+1,J ,K )-W(I ,J ,K ), C U(I ,J ,K )) C -V(I,J,K)*CVMGP(W(I ,J ,K )-W(I ,JM1,K ), C W(I ,JP1,K )-W(I ,J ,K ), C V(I ,J ,K )) C -W(I,J,K)*CVMGP(W(I ,J ,K )-W(I ,J ,KM1), C W(I ,J ,KP1)-W(I ,J ,K ), C W(I ,J ,K )) 105 CONTINUE ctim do 1 nt=1,10 ctim timer(nt) = 0.0 ctim1 continue ctim t1 = second() cser DO 2 JJ=1,NG cser DO 2 II=1,NG cser FIRSTN(II,JJ)=0 cser NUMBER(II,JJ)=0 cser2 CONTINUE ctim t2 = second() ctim timer(1) = t2 - t1 C C SORT THE XFN DATA BY X-COORDINATE AND Y-COORDINATE USING LINKED LISTS C ctim t1 = t2 cser DO 3 N=1,NFSIZE cser JJ=MODNG(INT(XFN(2,N)+FLNG)-1) + 1 cser II=MODNG(INT(XFN(1,N)+FLNG)-1) + 1 cser NEXTN(N)=FIRSTN(II,JJ) cser FIRSTN(II,JJ)=N cser NUMBER(II,JJ)=NUMBER(II,JJ)+1 cser3 CONTINUE ctim t2 = second() ctim timer(2) = t2 - t1 C C UPDATE THE LINKED LISTS WHICH SORT THE XFN DATA BY (X,Y) COORDINATE C DO 89 IZERO=1,4 DO 89 JZERO=1,4 CMIC$ DOALL CMIC$1SHARED(FIRSTN, NEXTN, NUMBER, IZERO, JZERO, XFN) CMIC$2PRIVATE(II, IJ, IN, JJ, JN, N, NEXTNOLD, NPREV) DO 85 IJ=0,NGBY4**2-1 JJ = JZERO + (IJ/NGBY4)*4 II = IZERO + (MOD(IJ,NGBY4))*4 IF (NUMBER(II,JJ) .EQ. 0) GO TO 85 NPREV = 0 N = FIRSTN(II,JJ) 82 JN = MODNG(INT(XFN(2,N)+FLNG)-1) + 1 IN = MODNG(INT(XFN(1,N)+FLNG)-1) + 1 IF ((IN.NE.II) .OR. (JN.NE.JJ)) THEN c c point N is in the wrong linked list c remember the pointer to the next one in the present linked list c NEXTNOLD = NEXTN(N) c c add point N to the correct linked list c NEXTN(N) = FIRSTN(IN,JN) FIRSTN(IN,JN) = N NUMBER(IN,JN) = NUMBER(IN,JN) + 1 c c remove point N from the present linked list c IF (NPREV .EQ. 0) THEN FIRSTN(II,JJ) = NEXTNOLD ELSE NEXTN(NPREV) = NEXTNOLD END IF NUMBER(II,JJ) = NUMBER(II,JJ) - 1 c c advance to the remembered next point in the present linked list c the previous-point pointer is unchanged c N = NEXTNOLD c ELSE c c point N is in the correct linked list c advance the previous-point pointer and then the current-point pointer c NPREV = N N = NEXTN(N) c END IF IF (N .NE. 0) GO TO 82 85 CONTINUE 89 CONTINUE C C FIRSTN(I,J) CONTAINS THE INDEX OF THE FIRST POINT WHICH SIMULTANEOUSLY IS C BETWEEN PLANES I AND I+1 AND IS BETWEEN PLANES J AND J+1. C THIS POINT HAS COORDINATES: C (XFN(1,FIRSTN(I,J)),XFN(2,FIRSTN(I,J)),XFN(3,FIRSTN(I,J))). C NEXTN(FIRSTN(I,J)) CONTAINS THE INDEX OF THE SECOND SUCH POINT C NEXTN(NEXTN(FIRSTN(I,J))) CONTAINS THE INDEX OF THE THIRD SUCH POINT C ETC. C IF FIRSTN(I,J) CONTAINS THE VALUE 0, THERE ARE NO SUCH POINTS. C cser DO 20 JJ=1,NG cser DO 20 II=1,NG MSHIFT = 16*NG DO 25 JZERO=1,4 DO 25 IZERO=1,4 CMIC$ DO ALL CMIC$1SHARED(FIRSTN, FRC, MSHIFT, NEXTN, NUMBER, UR, VR, WR, XFN, CMIC$2 IZERO, JZERO) CMIC$3PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, CMIC$4 FLIZP1, FLJZP1, FLKZP1, FORCE1, FORCE2, FORCE3, CMIC$5 IJ, I, I0, I3D, II, INDUVW, IZ, J, J3D, JJ, JZ, K, KZ, L, CMIC$6 LINDX, M, MZERO, NPOINTS, NPT, NUMREM, RAD1, RAD2, RAD3, CMIC$7 ULIN, VLIN, WLIN, XFN1OLD, XFN2OLD, XFN3OLD) DO 20 IJ=0,NGBY4**2-1 JJ = JZERO + (IJ/NGBY4)*4 II = IZERO + (MOD(IJ,NGBY4))*4 L = FIRSTN(II,JJ) IF (L.EQ.0) GO TO 20 C C CONSTRUCT THE INDEX ARRAY FOR SUBSEQUENT GATHERING AND SCATTERING OF THE C VELOCITY DATA FOR THE APPROPRIATE 4 X 4 X NG COLUMN OF THE 3D GRID INTO C AND OUT OF THE APPROPRIATE LINEAR ARRAYS, ONE FOR EACH DIRECTION. C ZERO-FILL THE LINEAR ARRAYS. NOTICE THAT THE LINEAR DATA ARRAYS ARE LARGER C THAN THE INDEX ARRAY; THIS FACILITATES SATISFYING THE PERIODIC B.C. C ctim t1 = t2 DO 4 K=0,NGM1 DO 4 J=-1,2 J3D = MODNG(JJ+J-1)+1 I0 = K*16+(J+1)*4 + 2 DO 4 I=-1,2 I3D = MODNG(II+I-1)+1 INDUVW(I0+I) = (K*(NB+1)**2 + J3D*(NB+1) + I3D)*2 + 1 4 CONTINUE DO 41 M=-15,0 ULIN(M) = 0.0 VLIN(M) = 0.0 WLIN(M) = 0.0 41 CONTINUE CALL GATHER(MSHIFT,ULIN(1),UR(0,0,0),INDUVW(1)) CALL GATHER(MSHIFT,VLIN(1),VR(0,0,0),INDUVW(1)) CALL GATHER(MSHIFT,WLIN(1),WR(0,0,0),INDUVW(1)) DO 42 M=16*NG+1,16*NGP2 ULIN(M) = 0.0 VLIN(M) = 0.0 WLIN(M) = 0.0 42 CONTINUE ctim t2 = second() ctim timer(3) = timer(3) + t2 - t1 NUMREM = NUMBER(II,JJ) 5 IF (NUMREM .GE. NPTMAX) THEN NPOINTS = NPTMAX NUMREM = NUMREM - NPTMAX ELSE NPOINTS = NUMREM NUMREM = 0 END IF ctim t1 = t2 DO 601 NPT=1,NPOINTS LINDX(NPT) = (L-1)*3 + 1 L = NEXTN(L) 601 CONTINUE CALL GATHER(NPOINTS,XFN1OLD,XFN(1,1),LINDX) CALL GATHER(NPOINTS,XFN2OLD,XFN(2,1),LINDX) CALL GATHER(NPOINTS,XFN3OLD,XFN(3,1),LINDX) CALL GATHER(NPOINTS, FORCE1,FRC(1,1),LINDX) CALL GATHER(NPOINTS, FORCE2,FRC(2,1),LINDX) CALL GATHER(NPOINTS, FORCE3,FRC(3,1),LINDX) IZ = INT(XFN1OLD( 1) - 1. + FLNG) - NG JZ = INT(XFN2OLD( 1) - 1. + FLNG) - NG FLIZP1 = IZ+1 FLJZP1 = JZ+1 DO 6 NPT=1,NPOINTS KZ = INT(XFN3OLD(NPT) - 1. + FLNG) - NG FLKZP1(NPT) = KZ+1 6 CONTINUE ctim t2 = second() ctim timer(4) = timer(4) + t2 - t1 ctim t1 = t2 DO 7 NPT=1,NPOINTS ARG3 = XFN3OLD(NPT) - FLKZP1(NPT) ARG2 = XFN2OLD(NPT) - FLJZP1 ARG1 = XFN1OLD(NPT) - FLIZP1 RAD3 = SQRT(1.+4.*ARG3*(1.-ARG3)) RAD2 = SQRT(1.+4.*ARG2*(1.-ARG2)) RAD1 = SQRT(1.+4.*ARG1*(1.-ARG1)) D3(NPT,3) = (1.+2.*ARG3-RAD3)/8. D2(NPT,3) = (1.+2.*ARG2-RAD2)/8. D1(NPT,3) = (1.+2.*ARG1-RAD1)/8. D3(NPT,2) = (1.+2.*ARG3+RAD3)/8. D2(NPT,2) = (1.+2.*ARG2+RAD2)/8. D1(NPT,2) = (1.+2.*ARG1+RAD1)/8. D3(NPT,1) = (3.-2.*ARG3+RAD3)/8. D2(NPT,1) = (3.-2.*ARG2+RAD2)/8. D1(NPT,1) = (3.-2.*ARG1+RAD1)/8. D3(NPT,0) = (3.-2.*ARG3-RAD3)/8. D2(NPT,0) = (3.-2.*ARG2-RAD2)/8. D1(NPT,0) = (3.-2.*ARG1-RAD1)/8. 7 CONTINUE ctim t2 = second() ctim timer(5) = timer(5) + t2 - t1 ctim t1 = t2 DO 9 J=0,3 DO 9 I=0,3 DO 9 NPT=1,NPOINTS D12(NPT,I,J) = D1(NPT,I)*D2(NPT,J) 9 CONTINUE ctim t2 = second() ctim timer(6) = timer(6) + t2 - t1 ctim t1 = t2 DO 10 K=0,3 DO 10 J=0,3 DO 10 I=0,3 M = K*16 + J*4 + I + 1 DO 10 NPT=1,NPOINTS DELTA(M,NPT) = D12(NPT,I,J) * D3(NPT,K) 10 CONTINUE ctim t2 = second() ctim timer(7) = timer(7) + t2 - t1 ctim t1 = t2 DO 12 NPT=1,NPOINTS MZERO = 16*(INT(XFN3OLD(NPT) - 1. + FLNG) - NG) DO 11 M=1,64 ULIN(M+MZERO) = ULIN(M+MZERO) + DELTA(M,NPT)*FORCE1(NPT) VLIN(M+MZERO) = VLIN(M+MZERO) + DELTA(M,NPT)*FORCE2(NPT) WLIN(M+MZERO) = WLIN(M+MZERO) + DELTA(M,NPT)*FORCE3(NPT) 11 CONTINUE 12 CONTINUE ctim t2 = second() ctim timer(8) = timer(8) + t2 - t1 IF (NUMREM .NE. 0) GO TO 5 ctim t1 = t2 DO 14 M=1,32 ULIN(M) = ULIN(M) + ULIN(M+MSHIFT) VLIN(M) = VLIN(M) + VLIN(M+MSHIFT) WLIN(M) = WLIN(M) + WLIN(M+MSHIFT) 14 CONTINUE DO 15 M=16*(NG-1)+1,16*NG ULIN(M) = ULIN(M) + ULIN(M-MSHIFT) VLIN(M) = VLIN(M) + VLIN(M-MSHIFT) WLIN(M) = WLIN(M) + WLIN(M-MSHIFT) 15 CONTINUE ctim t2 = second() ctim timer(9) = timer(9) + t2 - t1 ctim t1 = t2 CALL SCATTER(MSHIFT,UR(0,0,0),INDUVW(1),ULIN(1)) CALL SCATTER(MSHIFT,VR(0,0,0),INDUVW(1),VLIN(1)) CALL SCATTER(MSHIFT,WR(0,0,0),INDUVW(1),WLIN(1)) ctim t2 = second() ctim timer(10) = timer(10) + t2 - t1 20 CONTINUE 25 CONTINUE ctim write(6,'(a8,i14 )') 'PUSH ' ctim write(6,'(a8,f14.8)') ' DO 2 ',timer(1) ctim write(6,'(a8,f14.8)') ' DO 3 ',timer(2) ctim write(6,'(a8,f14.8)') ' DO 4 ',timer(3) ctim write(6,'(a8,f14.8)') ' DO 6 ',timer(4) ctim write(6,'(a8,f14.8)') ' DO 7 ',timer(5) ctim write(6,'(a8,f14.8)') ' DO 9 ',timer(6) ctim write(6,'(a8,f14.8)') ' DO 10 ',timer(7) ctim write(6,'(a8,f14.8)') ' UPDATE ',timer(8) ctim write(6,'(a8,f14.8)') ' FOLD ',timer(9) ctim write(6,'(a8,f14.8)') ' SCATTER',timer(10) RETURN END C********************************************************************** SUBROUTINE MOVE(XFN,NEXTN,NFG,NPF,NGROUPS,NBUNCH, c katapex,xatapex) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFSIZE=606638) PARAMETER(NPTMAX=256) PARAMETER(FLNG=NG) PARAMETER(NGP2=NG+2) PARAMETER(NGBY4=NG/4) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) INTEGER FIRSTN DIMENSION XFN(3,*),NEXTN(*) DIMENSION NFG(*),NPF(*),NGROUPS(*) DIMENSION XFN1OLD(NPTMAX),XFN2OLD(NPTMAX),XFN3OLD(NPTMAX) DIMENSION IMOD(0:4) , JMOD(0:4) , KMOD(0:4) DIMENSION D1( NPTMAX,0:3),D2(NPTMAX,0:3),D3(NPTMAX,0:3) DIMENSION D12(NPTMAX,0:3,0:3) DIMENSION DELTA(0:64,NPTMAX),INDEX(0:64) DIMENSION ULIN(-15:16*NGP2),VLIN(-15:16*NGP2),WLIN(-15:16*NGP2) DIMENSION UINT(0:64,NPTMAX),VINT(0:64,NPTMAX),WINT(0:64,NPTMAX) dimension xatapex(3) DIMENSION FLKZP1(NPTMAX) dimension lindx(NPTMAX) dimension timer(10) C C THIS ROUTINE MOVES THE FIBER POINTS AT THE LOCAL FLUID VELOCITY. C C FIBER POINT L IS INFLUENCED BY FLUID POINT (I,J,K) IFF C ALL THE FOLLOWING CONDITIONS ARE SATISFIED: C C ABS(I - XFN(1,L)) < 2. C ABS(J - XFN(2,L)) < 2. C ABS(K - XFN(3,L)) < 2. C C USING FORTRAN, THESE VALUES OF (I,J,K) ARE FOUND AS FOLLOWS: C C IZ = XFN(1,L) -1. C JZ = XFN(2,L) -1. C KZ = XFN(3,L) -1. C C THEN C I = IZ,IZ+3 C J = JZ,JZ+3 C K = KZ,KZ+3 C C IF THESE CONDITIONS ARE SATISFIED, THE COEFFICIENT THAT LINKS C FIBER POINT L TO FLUID POINT (I,J,K) IS: C C D = DEL(I-XFN(1,L)) * DEL(J-XFN(2,L)) * DEL(K-XFN(3,L)) C C WHERE C C DEL(R) = (1. + COS((PI/2.)*R))/4. C C ALGORITHM IS AS FOLLOWS: C C LET L=1,2, ... ,NP C FOR EACH L, LET (I,J,K) COVER RANGE OF VALUES DEFINED ABOVE. C C FOR EACH VALUE OF (L,I,J,K) C C XFN(1,L) = XFN(1,L) + D * U(I,J,K) C XFN(2,L) = XFN(2,L) + D * V(I,J,K) C XFN(3,L) = XFN(3,L) + D * W(I,J,K) C C WHERE D WAS DEFINED ABOVE. C MODNG(K)=MOD(K+NG,NG) DEL(R) = (1. + COS((PI/2.)*R))/4. PI = 4. * ATAN(1.) c rewind(41) c rewind(42) ctim do 1 nt=1,10 ctim timer(nt) = 0.0 ctim1 continue ctim t1 = second() cser DO 2 JJ=1,NG cser DO 2 II=1,NG cser FIRSTN(II,JJ)=0 cser NUMBER(II,JJ)=0 cser2 CONTINUE ctim t2 = second() ctim timer(1) = t2 - t1 C C SORT THE XFN DATA BY X-COORDINATE AND Y-COORDINATE USING LINKED LISTS C ctim t1 = t2 cser DO 3 N=1,NFSIZE cser JJ=MODNG(INT(XFN(2,N)+FLNG)-1) + 1 cser II=MODNG(INT(XFN(1,N)+FLNG)-1) + 1 cser NEXTN(N)=FIRSTN(II,JJ) cser FIRSTN(II,JJ)=N cser NUMBER(II,JJ)=NUMBER(II,JJ)+1 cser3 CONTINUE c call statist ctim t2 = second() ctim timer(2) = t2 - t1 C C FIRSTN(I,J) CONTAINS THE INDEX OF THE FIRST POINT WHICH SIMULTANEOUSLY IS C BETWEEN PLANES I AND I+1 AND IS BETWEEN PLANES J AND J+1. C THIS POINT HAS COORDINATES: C (XFN(1,FIRSTN(I,J)),XFN(2,FIRSTN(I,J)),XFN(3,FIRSTN(I,J))). C NEXTN(FIRSTN(I,J)) CONTAINS THE INDEX OF THE SECOND SUCH POINT C NEXTN(NEXTN(FIRSTN(I,J))) CONTAINS THE INDEX OF THE THIRD SUCH POINT C ETC. C IF FIRSTN(I,J) CONTAINS THE VALUE 0, THERE ARE NO SUCH POINTS. C cser DO 20 JJ=1,NG cser DO 20 II=1,NG CMIC$ DO ALL CMIC$1SHARED(FIRSTN, NEXTN, NUMBER, U, V, W, XFN) CMIC$2PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, CMIC$3 FLIZP1, FLJZP1, FLKZP1, CMIC$4 IJ, I, I0, I3D, II, IZ, J, J3D, JJ, JZ, K, K3D, KZ, L, CMIC$5 LINDX, M, MZERO, NPOINTS, NUMREM, NPT, RAD1, RAD2, RAD3, CMIC$6 UINT, ULIN, VINT, CMIC$7 VLIN, WINT, WLIN, XFN1OLD, XFN2OLD, XFN3OLD) DO 20 IJ=0,NG**2-1 JJ = 1 + IJ/NG II = 1 + MOD(IJ,NG) L = FIRSTN(II,JJ) IF (L.EQ.0) GO TO 20 C C COPY THE VELOCITY DATA FOR THE APPROPRIATE 4 X 4 X NG COLUMN OF THE 3D GRID C INTO THE APPROPRIATE LINEAR ARRAYS, ONE FOR EACH DIRECTION. NOTICE THAT K C DELIBERATELY COVERS A LARGER RANGE THAN K3D TO ENSURE PERIODICITY OF THE C LINEAR ARRAY VELOCITY DATA. C ctim t1 = t2 DO 4 K=-1,NG+1 K3D = MODNG(K) DO 4 J=-1,2 J3D = MODNG(JJ+J-1)+1 I0 = K*16+(J+1)*4 + 2 DO 4 I=-1,2 I3D = MODNG(II+I-1)+1 ULIN(I0+I) = U(I3D,J3D,K3D) VLIN(I0+I) = V(I3D,J3D,K3D) WLIN(I0+I) = W(I3D,J3D,K3D) 4 CONTINUE ctim t2 = second() ctim timer(3) = timer(3) + t2 - t1 NUMREM = NUMBER(II,JJ) 5 IF (NUMREM .GE. NPTMAX) THEN NPOINTS = NPTMAX NUMREM = NUMREM - NPTMAX ELSE NPOINTS = NUMREM NUMREM = 0 END IF ctim t1 = t2 DO 601 NPT=1,NPOINTS LINDX(NPT) = (L-1)*3 + 1 L = NEXTN(L) 601 CONTINUE CALL GATHER(NPOINTS,XFN1OLD,XFN(1,1),LINDX) CALL GATHER(NPOINTS,XFN2OLD,XFN(2,1),LINDX) CALL GATHER(NPOINTS,XFN3OLD,XFN(3,1),LINDX) IZ = INT(XFN1OLD( 1) - 1. + FLNG) - NG JZ = INT(XFN2OLD( 1) - 1. + FLNG) - NG FLIZP1 = IZ+1 FLJZP1 = JZ+1 DO 6 NPT=1,NPOINTS KZ = INT(XFN3OLD(NPT) - 1. + FLNG) - NG FLKZP1(NPT) = KZ+1 6 CONTINUE ctim t2 = second() ctim timer(4) = timer(4) + t2 - t1 ctim t1 = t2 DO 7 NPT=1,NPOINTS ARG3 = XFN3OLD(NPT) - FLKZP1(NPT) ARG2 = XFN2OLD(NPT) - FLJZP1 ARG1 = XFN1OLD(NPT) - FLIZP1 RAD3 = SQRT(1.+4.*ARG3*(1.-ARG3)) RAD2 = SQRT(1.+4.*ARG2*(1.-ARG2)) RAD1 = SQRT(1.+4.*ARG1*(1.-ARG1)) D3(NPT,3) = (1.+2.*ARG3-RAD3)/8. D2(NPT,3) = (1.+2.*ARG2-RAD2)/8. D1(NPT,3) = (1.+2.*ARG1-RAD1)/8. D3(NPT,2) = (1.+2.*ARG3+RAD3)/8. D2(NPT,2) = (1.+2.*ARG2+RAD2)/8. D1(NPT,2) = (1.+2.*ARG1+RAD1)/8. D3(NPT,1) = (3.-2.*ARG3+RAD3)/8. D2(NPT,1) = (3.-2.*ARG2+RAD2)/8. D1(NPT,1) = (3.-2.*ARG1+RAD1)/8. D3(NPT,0) = (3.-2.*ARG3-RAD3)/8. D2(NPT,0) = (3.-2.*ARG2-RAD2)/8. D1(NPT,0) = (3.-2.*ARG1-RAD1)/8. 7 CONTINUE ctim t2 = second() ctim timer(5) = timer(5) + t2 - t1 ctim t1 = t2 DO 9 J=0,3 DO 9 I=0,3 DO 9 NPT=1,NPOINTS D12(NPT,I,J) = D1(NPT,I)*D2(NPT,J) 9 CONTINUE ctim t2 = second() ctim timer(6) = timer(6) + t2 - t1 ctim t1 = t2 DO 10 K=0,3 DO 10 J=0,3 DO 10 I=0,3 M = K*16 + J*4 + I + 1 DO 10 NPT=1,NPOINTS DELTA(M,NPT) = D12(NPT,I,J) * D3(NPT,K) 10 CONTINUE ctim t2 = second() ctim timer(7) = timer(7) + t2 - t1 ctim t1 = t2 DO 122 NPT=1,NPOINTS MZERO = 16*(INT(XFN3OLD(NPT) - 1. + FLNG) - NG) DO 121 M=1,64 UINT(M,NPT) = ULIN(MZERO+M)*DELTA(M,NPT) VINT(M,NPT) = VLIN(MZERO+M)*DELTA(M,NPT) WINT(M,NPT) = WLIN(MZERO+M)*DELTA(M,NPT) 121 CONTINUE 122 CONTINUE ctim t2 = second() ctim timer(8) = timer(8) + t2 - t1 ctim t1 = t2 DO 124 M=1,64 DO 123 NPT=1,NPOINTS XFN1OLD(NPT) = XFN1OLD(NPT) + UINT(M,NPT) XFN2OLD(NPT) = XFN2OLD(NPT) + VINT(M,NPT) XFN3OLD(NPT) = XFN3OLD(NPT) + WINT(M,NPT) 123 CONTINUE 124 CONTINUE ctim t2 = second() ctim timer(9) = timer(9) + t2 - t1 ctim t1 = t2 CALL SCATTER(NPOINTS,XFN(1,1),LINDX,XFN1OLD) CALL SCATTER(NPOINTS,XFN(2,1),LINDX,XFN2OLD) CALL SCATTER(NPOINTS,XFN(3,1),LINDX,XFN3OLD) ctim t2 = second() ctim timer(10) = timer(10) + t2 - t1 IF (NUMREM .NE. 0) GO TO 5 20 CONTINUE ctim write(6,'(a8 )') 'MOVE ' ctim write(6,'(a8,f14.8)') ' DO 2 ',timer(1) ctim write(6,'(a8,f14.8)') ' DO 3 ',timer(2) ctim write(6,'(a8,f14.8)') ' DO 4 ',timer(3) ctim write(6,'(a8,f14.8)') ' DO 6 ',timer(4) ctim write(6,'(a8,f14.8)') ' DO 7 ',timer(5) ctim write(6,'(a8,f14.8)') ' DO 9 ',timer(6) ctim write(6,'(a8,f14.8)') ' DO 10 ',timer(7) ctim write(6,'(a8,f14.8)') ' INTERP ',timer(8) ctim write(6,'(a8,f14.8)') ' MOVE ',timer(9) ctim write(6,'(a8,f14.8)') ' SCATTER',timer(10) C C XFN IS NOW THE FIBER CONFIGURATION AT THE END OF THE CURRENT TIME STEP C THE FIBER FORCES DURING THE NEXT TIME STEP C xatapex(1) = xfn(1,katapex) xatapex(2) = xfn(2,katapex) xatapex(3) = xfn(3,katapex) RETURN END subroutine chkuvw PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) umax = 0.0 vmax = 0.0 wmax = 0.0 do 20 kk=0,ngm1 call uvwmax(u(0,1,kk),v(0,1,kk),w(0,1,kk),umax,vmax,wmax) 20 continue call speedmax(smax) write(6,*)'umax,vmax,wmax = ',umax,vmax,wmax write(6,*)'smax = ',smax write(6,*) return end SUBROUTINE UVWMAX(U,V,W,UMAX,VMAX,WMAX) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2,NBP1=NB+1) PARAMETER(NSIZE=(NB+1)*NG) DIMENSION U(0:NB,1:NG) DIMENSION V(0:NB,1:NG) DIMENSION W(0:NB,1:NG) DIMENSION ULIN (NG),VLIN (NG),WLIN (NG) INTEGER UINDEX(NG),VINDEX(NG),WINDEX(NG) INTEGER UI0 ,VI0 ,WI0 C C FOR EACH VELOCITY COMPONENT, U, V, W: C C FIND THE ROW INDEX OF THE MAXIMUM VALUE IN EACH COLUMN C DO 10 J=1,NG UINDEX(J) = ISAMAX(NG,U(1,J),1) VINDEX(J) = ISAMAX(NG,V(1,J),1) WINDEX(J) = ISAMAX(NG,W(1,J),1) 10 CONTINUE C C CONVERT EACH ROW INDEX INTO A LOCATION IN THE (0:NB,1:NG) ARRAY C DO 20 J=1,NG UINDEX(J) = UINDEX(J) + 1 + NBP1*(J-1) VINDEX(J) = VINDEX(J) + 1 + NBP1*(J-1) WINDEX(J) = WINDEX(J) + 1 + NBP1*(J-1) 20 CONTINUE C C GATHER THE MAXIMUM VALUE FROM EACH COLUMN INTO A LINEAR ARRAY C CALL GATHER(NG,ULIN(1),U(0,1),UINDEX(1)) CALL GATHER(NG,VLIN(1),V(0,1),VINDEX(1)) CALL GATHER(NG,WLIN(1),W(0,1),WINDEX(1)) C C FIND THE INDEX OF THE MAXIMUM VALUE IN THE LINEAR ARRAY C UI0 = ISAMAX(NG,ULIN(1),1) VI0 = ISAMAX(NG,VLIN(1),1) WI0 = ISAMAX(NG,WLIN(1),1) C C UPDATE THE OVERALL MAXIMUM VALUE C IF (ABS(ULIN(UI0)) .GT. UMAX) UMAX = ABS(ULIN(UI0)) IF (ABS(VLIN(VI0)) .GT. VMAX) VMAX = ABS(VLIN(VI0)) IF (ABS(WLIN(WI0)) .GT. WMAX) WMAX = ABS(WLIN(WI0)) RETURN END SUBROUTINE CHKSPEED(SPEEDM) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2,NBP1=NB+1) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) DIMENSION S(0:NB,1:NG) DIMENSION UMAX(0:NGM1),VMAX(0:NGM1),WMAX(0:NGM1),SMAX(0:NGM1) DIMENSION ULIN (NG) ,VLIN (NG) ,WLIN (NG) ,SLIN (NG) INTEGER UINDEX(NG) ,VINDEX(NG) ,WINDEX(NG) ,SINDEX(NG) INTEGER UI0 ,VI0 ,WI0 ,SI0 INTEGER UIM ,VIM ,WIM ,SIM CMIC$ DO ALL CMIC$1SHARED(U, V, W, UMAX, VMAX, WMAX, SMAX) CMIC$2PRIVATE( ULIN , VLIN , WLIN , SLIN , CMIC$3 UINDEX, VINDEX, WINDEX, SINDEX, CMIC$4 UI0 , VI0 , WI0 , SI0 , CMIC$5 I , J , K , S ) DO 100 K=0,NGM1 C C FOR EACH PLANE K, CREATE THE ARRAY S WHOSE MAXIMUMM IS TO BE FOUND C DO 5 J=1,NG DO 5 I=1,NG S(I,J) = ABS(U(I,J,K)) + ABS(V(I,J,K)) + ABS(W(I,J,K)) 5 CONTINUE C C FIND THE ROW INDEX OF THE MAXIMUM VALUE IN EACH COLUMN C DO 10 J=1,NG UINDEX(J) = ISAMAX(NG,U(1,J,K),1) VINDEX(J) = ISAMAX(NG,V(1,J,K),1) WINDEX(J) = ISAMAX(NG,W(1,J,K),1) SINDEX(J) = ISAMAX(NG,S(1,J) ,1) 10 CONTINUE C C CONVERT EACH ROW INDEX INTO A LOCATION IN THE (0:NB,1:NG) ARRAY C DO 20 J=1,NG UINDEX(J) = UINDEX(J) + 1 + NBP1*(J-1) VINDEX(J) = VINDEX(J) + 1 + NBP1*(J-1) WINDEX(J) = WINDEX(J) + 1 + NBP1*(J-1) SINDEX(J) = SINDEX(J) + 1 + NBP1*(J-1) 20 CONTINUE C C GATHER THE MAXIMUM VALUE FROM EACH COLUMN INTO A LINEAR ARRAY C CALL GATHER(NG,ULIN(1),U(0,1,K),UINDEX(1)) CALL GATHER(NG,VLIN(1),V(0,1,K),VINDEX(1)) CALL GATHER(NG,WLIN(1),W(0,1,K),WINDEX(1)) CALL GATHER(NG,SLIN(1),S(0,1 ),SINDEX(1)) C C FIND THE INDEX OF THE MAXIMUM VALUE IN THE LINEAR ARRAY C UI0 = ISAMAX(NG,ULIN(1),1) VI0 = ISAMAX(NG,VLIN(1),1) WI0 = ISAMAX(NG,WLIN(1),1) SI0 = ISAMAX(NG,SLIN(1),1) C C STORE THE MAXIMUM VALUE IN THE LOCATION CORRESPONDING TO PLANE K C UMAX(K) = ABS(ULIN(UI0)) VMAX(K) = ABS(VLIN(VI0)) WMAX(K) = ABS(WLIN(WI0)) SMAX(K) = ABS(SLIN(SI0)) C 100 CONTINUE C C FIND THE INDEX OF THE OVERALL MAXIMUM VALUE, NOTING DIMENSIONS (0:NGM1) C UIM = ISAMAX(NG,UMAX(0),1) - 1 VIM = ISAMAX(NG,VMAX(0),1) - 1 WIM = ISAMAX(NG,WMAX(0),1) - 1 SIM = ISAMAX(NG,SMAX(0),1) - 1 C C PRINT OUT THE MAXIMUM VALUES C write(6,*)'umax,vmax,wmax = ',umax(uim),vmax(vim),wmax(wim) write(6,*)'smax = ',smax(sim) speedm = smax(sim) write(6,*) RETURN END SUBROUTINE MOVEMK(XMK,NEXTN,NMARKT) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP2=NG+2) PARAMETER(NMSIZE=7593) PARAMETER(NPTMAX=256) PARAMETER(NCLMAX=19) PARAMETER(FLNG=NG) PARAMETER(NGBY4=NG/4) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/HIST/ FDUMMY(1:NG,1:NG) COMMON/HIST/ NDUMMY(1:NG,1:NG) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) INTEGER FDUMMY INTEGER FIRSTN DIMENSION XMK(3,*),NEXTN(*) DIMENSION XMK1OLD(NPTMAX),XMK2OLD(NPTMAX),XMK3OLD(NPTMAX) DIMENSION D1( NPTMAX,0:3),D2(NPTMAX,0:3),D3(NPTMAX,0:3) DIMENSION D12(NPTMAX,0:3,0:3) DIMENSION DELTA(0:64,NPTMAX),INDEX(0:64) DIMENSION ULIN(-15:16*NGP2),VLIN(-15:16*NGP2),WLIN(-15:16*NGP2) DIMENSION UINT(0:64,NPTMAX),VINT(0:64,NPTMAX),WINT(0:64,NPTMAX) DIMENSION FLKZP1(NPTMAX) dimension lindx(NPTMAX) C C THIS ROUTINE MOVES THE MARKER POINTS AT THE LOCAL FLUID VELOCITY. C C MARKER POINT L IS INFLUENCED BY FLUID POINT (I,J,K) IFF C ALL THE FOLLOWING CONDITIONS ARE SATISFIED: C C ABS(I - XMK(1,L)) < 2. C ABS(J - XMK(2,L)) < 2. C ABS(K - XMK(3,L)) < 2. C C USING FORTRAN, THESE VALUES OF (I,J,K) ARE FOUND AS FOLLOWS: C C IZ = XMK(1,L) -1. C JZ = XMK(2,L) -1. C KZ = XMK(3,L) -1. C C THEN C I = IZ,IZ+3 C J = JZ,JZ+3 C K = KZ,KZ+3 C C IF THESE CONDITIONS ARE SATISFIED, THE COEFFICIENT THAT LINKS C MARKER POINT L TO FLUID POINT (I,J,K) IS: C C D = DEL(I-XMK(1,L)) * DEL(J-XMK(2,L)) * DEL(K-XMK(3,L)) C C WHERE C C DEL(R) = (1. + COS((PI/2.)*R))/4. C C ALGORITHM IS AS FOLLOWS: C C LET L=1,2, ... ,NP C FOR EACH L, LET (I,J,K) COVER RANGE OF VALUES DEFINED ABOVE. C C FOR EACH VALUE OF (L,I,J,K) C C XMK(1,L) = XMK(1,L) + D * U(I,J,K) C XMK(2,L) = XMK(2,L) + D * V(I,J,K) C XMK(3,L) = XMK(3,L) + D * W(I,J,K) C C WHERE D WAS DEFINED ABOVE. C MODNG(K)=MOD(K+NG,NG) DEL(R) = (1. + COS((PI/2.)*R))/4. PI = 4. * ATAN(1.) cser DO 2 KK=0,NGM1 cser FIRSTN(KK)=0 cser NUMBER(KK)=0 cser2 CONTINUE C C SORT THE XMK DATA BY Z-COORDINATE USING LINKED LISTS C cser DO 3 N=1,NMARKT C KK=XMK(3,N) cser KK=MODNG(INT(XMK(3,N)+FLNG)) cser NEXTN(N)=FIRSTN(KK) cser FIRSTN(KK)=N cser NUMBER(KK)=NUMBER(KK)+1 cser3 CONTINUE C C UPDATE THE LINKED LISTS WHICH SORT THE XMK DATA BY (X,Y) COORDINATE C DO 89 IZERO=1,4 DO 89 JZERO=1,4 CMIC$ DOALL CMIC$1SHARED(FIRSTN, NEXTN, NUMBER, IZERO, JZERO, XMK) CMIC$2PRIVATE(II, IJ, IN, JJ, JN, N, NEXTNOLD, NPREV) DO 85 IJ=0,NGBY4**2-1 JJ = JZERO + (IJ/NGBY4)*4 II = IZERO + (MOD(IJ,NGBY4))*4 IF (NUMBER(II,JJ) .EQ. 0) GO TO 85 NPREV = 0 N = FIRSTN(II,JJ) 82 JN = MODNG(INT(XMK(2,N)+FLNG)-1) + 1 IN = MODNG(INT(XMK(1,N)+FLNG)-1) + 1 IF ((IN.NE.II) .OR. (JN.NE.JJ)) THEN c c point N is in the wrong linked list c remember the pointer to the next one in the present linked list c NEXTNOLD = NEXTN(N) c c add point N to the correct linked list c NEXTN(N) = FIRSTN(IN,JN) FIRSTN(IN,JN) = N NUMBER(IN,JN) = NUMBER(IN,JN) + 1 c c remove point N from the present linked list c IF (NPREV .EQ. 0) THEN FIRSTN(II,JJ) = NEXTNOLD ELSE NEXTN(NPREV) = NEXTNOLD END IF NUMBER(II,JJ) = NUMBER(II,JJ) - 1 c c advance to the remembered next point in the present linked list c the previous-point pointer is unchanged c N = NEXTNOLD c ELSE c c point N is in the correct linked list c advance the previous-point pointer and then the current-point pointer c NPREV = N N = NEXTN(N) c END IF IF (N .NE. 0) GO TO 82 85 CONTINUE 89 CONTINUE C C FIRSTN(I,J) CONTAINS THE INDEX OF THE FIRST POINT WHICH SIMULTANEOUSLY IS C BETWEEN PLANES I AND I+1 AND IS BETWEEN PLANES J AND J+1. C THIS POINT HAS COORDINATES: C (XMK(1,FIRSTN(I,J)),XMK(2,FIRSTN(I,J)),XMK(3,FIRSTN(I,J))). C NEXTN(FIRSTN(I,J)) CONTAINS THE INDEX OF THE SECOND SUCH POINT C NEXTN(NEXTN(FIRSTN(I,J))) CONTAINS THE INDEX OF THE THIRD SUCH POINT C ETC. C IF FIRSTN(I,J) CONTAINS THE VALUE 0, THERE ARE NO SUCH POINTS. C CMIC$ DO ALL CMIC$1SHARED(FIRSTN, NEXTN, NUMBER, U, V, W, XMK) CMIC$2PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, CMIC$3 FLIZP1, FLJZP1, FLKZP1, CMIC$4 IJ, I, I0, I3D, II, IZ, J, J3D, JJ, JZ, K, K3D, KZ, L, CMIC$5 LINDX, M, MZERO, NPOINTS, NUMREM, NPT, RAD1, RAD2, RAD3, CMIC$6 UINT, ULIN, VINT, CMIC$7 VLIN, WINT, WLIN, XMK1OLD, XMK2OLD, XMK3OLD) DO 20 IJ=0,NG**2-1 JJ = 1 + IJ/NG II = 1 + MOD(IJ,NG) L = FIRSTN(II,JJ) IF (L.EQ.0) GO TO 20 C C COPY THE VELOCITY DATA FOR THE APPROPRIATE 4 X 4 X NG COLUMN OF THE 3D GRID C INTO THE APPROPRIATE LINEAR ARRAYS, ONE FOR EACH DIRECTION. NOTICE THAT K C DELIBERATELY COVERS A LARGER RANGE THAN K3D TO ENSURE PERIODICITY OF THE C LINEAR ARRAY VELOCITY DATA. C DO 4 K=-1,NG+1 K3D = MODNG(K) DO 4 J=-1,2 J3D = MODNG(JJ+J-1)+1 I0 = K*16+(J+1)*4 + 2 DO 4 I=-1,2 I3D = MODNG(II+I-1)+1 ULIN(I0+I) = U(I3D,J3D,K3D) VLIN(I0+I) = V(I3D,J3D,K3D) WLIN(I0+I) = W(I3D,J3D,K3D) 4 CONTINUE NUMREM = NUMBER(II,JJ) 5 IF (NUMREM .GE. NPTMAX) THEN NPOINTS = NPTMAX NUMREM = NUMREM - NPTMAX ELSE NPOINTS = NUMREM NUMREM = 0 END IF DO 601 NPT=1,NPOINTS LINDX(NPT) = (L-1)*3 + 1 L = NEXTN(L) 601 CONTINUE CALL GATHER(NPOINTS,XMK1OLD,XMK(1,1),LINDX) CALL GATHER(NPOINTS,XMK2OLD,XMK(2,1),LINDX) CALL GATHER(NPOINTS,XMK3OLD,XMK(3,1),LINDX) IZ = INT(XMK1OLD( 1) - 1. + FLNG) - NG JZ = INT(XMK2OLD( 1) - 1. + FLNG) - NG FLIZP1 = IZ+1 FLJZP1 = JZ+1 DO 6 NPT=1,NPOINTS KZ = INT(XMK3OLD(NPT) - 1. + FLNG) - NG FLKZP1(NPT) = KZ+1 6 CONTINUE DO 7 NPT=1,NPOINTS ARG3 = XMK3OLD(NPT) - FLKZP1(NPT) ARG2 = XMK2OLD(NPT) - FLJZP1 ARG1 = XMK1OLD(NPT) - FLIZP1 RAD3 = SQRT(1.+4.*ARG3*(1.-ARG3)) RAD2 = SQRT(1.+4.*ARG2*(1.-ARG2)) RAD1 = SQRT(1.+4.*ARG1*(1.-ARG1)) D3(NPT,3) = (1.+2.*ARG3-RAD3)/8. D2(NPT,3) = (1.+2.*ARG2-RAD2)/8. D1(NPT,3) = (1.+2.*ARG1-RAD1)/8. D3(NPT,2) = (1.+2.*ARG3+RAD3)/8. D2(NPT,2) = (1.+2.*ARG2+RAD2)/8. D1(NPT,2) = (1.+2.*ARG1+RAD1)/8. D3(NPT,1) = (3.-2.*ARG3+RAD3)/8. D2(NPT,1) = (3.-2.*ARG2+RAD2)/8. D1(NPT,1) = (3.-2.*ARG1+RAD1)/8. D3(NPT,0) = (3.-2.*ARG3-RAD3)/8. D2(NPT,0) = (3.-2.*ARG2-RAD2)/8. D1(NPT,0) = (3.-2.*ARG1-RAD1)/8. 7 CONTINUE DO 9 J=0,3 DO 9 I=0,3 DO 9 NPT=1,NPOINTS D12(NPT,I,J) = D1(NPT,I)*D2(NPT,J) 9 CONTINUE DO 10 K=0,3 DO 10 J=0,3 DO 10 I=0,3 M = K*16 + J*4 + I + 1 DO 10 NPT=1,NPOINTS DELTA(M,NPT) = D12(NPT,I,J) * D3(NPT,K) 10 CONTINUE DO 122 NPT=1,NPOINTS MZERO = 16*(INT(XMK3OLD(NPT) - 1. + FLNG) - NG) DO 121 M=1,64 UINT(M,NPT) = ULIN(MZERO+M)*DELTA(M,NPT) VINT(M,NPT) = VLIN(MZERO+M)*DELTA(M,NPT) WINT(M,NPT) = WLIN(MZERO+M)*DELTA(M,NPT) 121 CONTINUE 122 CONTINUE DO 124 M=1,64 DO 123 NPT=1,NPOINTS XMK1OLD(NPT) = XMK1OLD(NPT) + UINT(M,NPT) XMK2OLD(NPT) = XMK2OLD(NPT) + VINT(M,NPT) XMK3OLD(NPT) = XMK3OLD(NPT) + WINT(M,NPT) 123 CONTINUE 124 CONTINUE CALL SCATTER(NPOINTS,XMK(1,1),LINDX,XMK1OLD) CALL SCATTER(NPOINTS,XMK(2,1),LINDX,XMK2OLD) CALL SCATTER(NPOINTS,XMK(3,1),LINDX,XMK3OLD) IF (NUMREM .NE. 0) GO TO 5 20 CONTINUE RETURN END C********************************************************************** SUBROUTINE INMARK(XMK,ncloud,nmarks,ncircs,nmarkt) c c This subroutine reads in the fluid markers from a named external file c and initializes some variables used by the flowmeter routines. c PARAMETER(L2NG=7,NG=2**L2NG) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) parameter(nmarmx=98,ncirmx=25) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) COMMON/BULGE/FLNG2,HZ,HSCALE,XSHIFT,YSHIFT,ZSHIFT dimension xmk(3,*) dimension nmarks(nclmax) dimension axvalv(4),byvalv(4),xcvalv(4),ycvalv(4) dimension axvein(4),byvein(4),xcvein(4),ycvein(4) dimension axvena(2),byvena(2),xcvena(2),ycvena(2) dimension zsrc(nsrcs) c open(4,file='/m/h2/bio0001/mcqueen/markers.in', c c form='formatted',access='sequential',status='old') c open(4,file='/m/s1/bio0002/mcqueen/markers.in', c c form='formatted',access='sequential',status='old') open(4,file='markers.in', c form='formatted',access='sequential',status='old') read(4,*) ngdata read(4,*) height read(4,*) thetat read(4,*) thetal read(4,*) thetaw read(4,*) thetar read(4,*) nest read(4,*) hrings do 10 iv=1,4 read(4,*) axvalv(iv),byvalv(iv),xcvalv(iv),ycvalv(iv) 10 continue read(4,*) hveins do 11 iv=1,4 read(4,*) axvein(iv),byvein(iv),xcvein(iv),ycvein(iv) 11 continue do 12 iv=1,2 read(4,*) axvena(iv),byvena(iv),xcvena(iv),ycvena(iv) 12 continue do 13 is=1,nsrcs read(4,*) zsrc(is) 13 continue read(4,*) ncloud c check memory allocation for array nmarks if (ncloud .gt. nclmax) then write(6,*) ' number of clouds = ',ncloud,' too large in INMARK' call exit(1) end if read(4,*) scale radmax = axvalv(1) if (byvalv( 1) .gt. radmax) radmax = byvalv( 1) do 15 iv=2,4 if (axvalv(iv) .gt. radmax) radmax = axvalv(iv) if (byvalv(iv) .gt. radmax) radmax = byvalv(iv) 15 continue ncircs = radmax*scale + 1. c check memory allocation in subroutine flows for "flowmeter arrays" if (ncircs .gt. ncirmx) then write(6,*) ' ncircs = ',ncircs,' too large in INMARK' call exit(1) end if c scale up the lengths to accomodate the new computational grid xmscal = float(ng)/float(ngdata) height = height*xmscal hrings = hrings*xmscal hveins = hveins*xmscal do 32 iv=1,4 axvalv(iv) = axvalv(iv)*xmscal byvalv(iv) = byvalv(iv)*xmscal xcvalv(iv) = xcvalv(iv)*xmscal ycvalv(iv) = ycvalv(iv)*xmscal 32 continue do 33 iv=1,4 axvein(iv) = axvein(iv)*xmscal byvein(iv) = byvein(iv)*xmscal xcvein(iv) = xcvein(iv)*xmscal ycvein(iv) = ycvein(iv)*xmscal 33 continue do 34 iv=1,4 axvena(iv) = axvena(iv)*xmscal byvena(iv) = byvena(iv)*xmscal xcvena(iv) = xcvena(iv)*xmscal ycvena(iv) = ycvena(iv)*xmscal 34 continue do 35 is=1,nsrcs zsrc(is) = zsrc(is)*xmscal 35 continue thetat = thetat*atan(1.0)/45.0 kstop = 0 kstrt = 1 do 30 nc=1,ncloud read(4,*) nmarks(nc) if (nc .le. 4) then c check memory allocation in subroutine flows for "flowmeter arrays" if (nmarks(nc) .gt. nmarmx) then write(6,*)' number of marks = ',nmarks(nc), c ' on valve ring = ',nc,' too large in INMARK' call exit(1) end if end if read(4,*) kstop = kstop + nmarks(nc) if (kstop .gt. nmsize) then write(6,*)' overall number of markers = ',kstop, c ' exceeds the maximum nmsize = ',nmsize, c ' beginning with cloud ',nc call exit(1) end if do 20 k=kstrt,kstop read(4,*) xmk(1,k),xmk(2,k),xmk(3,k) xmk(1,k) = xmk(1,k)*xmscal xmk(2,k) = xmk(2,k)*xmscal xmk(3,k) = xmk(3,k)*xmscal 20 continue kstrt = kstrt + nmarks(nc) 30 continue c nmarkt = actual total number of markers (used in subroutine movemk) nmarkt = kstop C C APPLY TO EACH DATA POINT A RADIAL EXPANSION ABOUT THE NOW-VERTICAL AXIS C OF SYMMETRY WITH MAGNITUDE DETERMINED BY THE THIRD (Z) COMPONENT C OF THE DATA. THE DATA IS CONTAINED BETWEEN Z=0.0 AND Z=HZ. C HZ, HSCALE AND FLNG2 HAVE BEEN SET IN SUBROUTINE INFIBER. C c TTH = TAN(THETAT) c DO 400 K=1,KSTOP c IF (XMK(3,K) .LE. HZ) THEN c XCEN = FLNG2 + TTH*XMK(3,K) c YCEN = FLNG2 c SCALE = 1.0 + HSCALE* XMK(3,K)*(HZ-XMK(3,K)) c XMK(1,K) = XCEN + SCALE*(XMK(1,K)-XCEN) c XMK(2,K) = YCEN + SCALE*(XMK(2,K)-YCEN) c END IF c 400 CONTINUE C C XSHIFT, YSHIFT AND ZSHIFT HAVE BEEN SET IN SUBROUTINE INFIBER. C TRANSLATE THE MARKERS TO THE CENTER OF THE DOMAIN. C DO 600 K=1,KSTOP XMK(1,K) = XMK(1,K) + XSHIFT XMK(2,K) = XMK(2,K) + YSHIFT XMK(3,K) = XMK(3,K) + ZSHIFT 600 CONTINUE C C ROTATE THE MARKERS ABOUT THE ORIGIN SO THAT UNIT VECTORS FIXED TO C THE MARKERS AND ORIGINALLY POINTING IN THE X, Y AND Z DIRECTIONS C END UP POINTING IN THE Z, X AND Y DIRECTIONS, RESPECTIVELY. C DO 610 K=1,KSTOP XTMP = XMK(1,K) YTMP = XMK(2,K) ZTMP = XMK(3,K) XMK(1,K) = YTMP XMK(2,K) = ZTMP XMK(3,K) = XTMP 610 CONTINUE C C TRANSLATE THE MARKERS TO THE CENTER OF THE DOMAIN C THE ADDITIONAL 0.5 IN THE Z DIRECTION IS TO POSITION THE HEART C MIDWAY BETWEEN THE SINK ON PLANE 1 AND THE SINK ON PLANE NG C (I.E., THE PERIODIC IMAGE OF THE SINK ON PLANE 0) C XADD = FLOAT(NG)/2.0 YADD = FLOAT(NG)/2.0 ZADD = FLOAT(NG)/2.0 + 0.5 DO 620 K=1,KSTOP XMK(1,K) = XMK(1,K) + XADD XMK(2,K) = XMK(2,K) + YADD XMK(3,K) = XMK(3,K) + ZADD 620 CONTINUE close(4) RETURN END C********************************************************************** SUBROUTINE INHIST PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) PARAMETER(NMSIZE=7593) PARAMETER(NCLMAX=19) PARAMETER(NCONES=12) PARAMETER(FLNG=NG) COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE) COMMON/XAR/ STF0 ( NFSIZE) COMMON/XAR/ REST0 ( NFSIZE) COMMON/XAR/ FRC (3,NFSIZE) COMMON/XAR/ NEXTN ( NFSIZE),LAFLAG ( NFSIZE),ARFLAG ( NFSIZE) COMMON/XAR/ LAYER(IMAX*NFGMAX),NFIBER(0:NFGMAX,IMAX,0:NCONES) COMMON/XAR/ NGROUPS(NBUNCH) COMMON/XAR/ NFG(IMAX),NPF(IMAX),KSTART(IMAX),NFSTART(IMAX) COMMON/XAR/ MRAMP(IMAX),MFLAT LOGICAL ARFLAG COMMON/MAR/ XMK(3,NMSIZE) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) COMMON/HIST/ FIRSTM(1:NG,1:NG) COMMON/HIST/ NUMBEM(1:NG,1:NG) INTEGER FIRSTN INTEGER FIRSTM MODNG(K)=MOD(K+NG,NG) DO 2 JJ=1,NG DO 2 II=1,NG FIRSTN(II,JJ)=0 NUMBER(II,JJ)=0 FIRSTM(II,JJ)=0 NUMBEM(II,JJ)=0 2 CONTINUE DO 3 N=1,NFSIZE JJ=MODNG(INT(XFN(2,N)+FLNG)-1) + 1 II=MODNG(INT(XFN(1,N)+FLNG)-1) + 1 NEXTN(N)=FIRSTN(II,JJ) FIRSTN(II,JJ)=N NUMBER(II,JJ)=NUMBER(II,JJ)+1 3 CONTINUE DO 4 N=1,NMSIZE JJ=MODNG(INT(XMK(2,N)+FLNG)-1) + 1 II=MODNG(INT(XMK(1,N)+FLNG)-1) + 1 NEXTM(N)=FIRSTM(II,JJ) FIRSTM(II,JJ)=N NUMBEM(II,JJ)=NUMBEM(II,JJ)+1 4 CONTINUE RETURN END C********************************************************************** subroutine flows (xmk,nextnm,nextnc,nmarks,ncircs, c expnum,klok,klokout,nref) c c This subroutine acts as a driver for the flowmeter, c calling subroutine emflow for each of the 4 valve rings. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NGM1=NG-1) PARAMETER(NMSIZE=7593) parameter(nmarmx=98,ncirmx=25) dimension xmk(3,*) dimension nextnm(*) dimension nextnc(*) dimension nmarks(*) c "flowmeter arrays" dimension rad( 3,nmarmx) dimension xnorm( 3,nmarmx) dimension axnorm( nmarmx) dimension xweb( 3,nmarmx,0:ncirmx) dimension uweb( 3,nmarmx,0:ncirmx) dimension aweb( nmarmx,0:ncirmx) dimension acmp( 3,nmarmx,0:ncirmx) dimension xcen( 3,4) dimension flow( 4) integer expnum character*14 filefl character*64 command do 1 icom=1,64 command(icom:icom) = ' ' 1 continue flnref = nref 101 format(a4,i2,a4,i1,a3) 102 format(a4,i2,a3,i2,a3) 103 format(a4,i2,a2,i3,a3) 104 format(a4,i2,a1,i4,a3) 105 format(a4,i2, i5,a3) 201 format(a3,i3,a4,i1,a3) 202 format(a3,i3,a3,i2,a3) 203 format(a3,i3,a2,i3,a3) 204 format(a3,i3,a1,i4,a3) 205 format(a3,i3, i5,a3) if (expnum .lt. 100) then if ( klok .lt. 10 ) then write(filefl,101) 'hrtx',expnum,'k000',klok,'.fl' elseif (klok .lt. 100 ) then write(filefl,102) 'hrtx',expnum,'k00' ,klok,'.fl' elseif (klok .lt. 1000 ) then write(filefl,103) 'hrtx',expnum,'k0' ,klok,'.fl' elseif (klok .lt. 10000) then write(filefl,104) 'hrtx',expnum,'k' ,klok,'.fl' else write(filefl,105) 'hrtx',expnum, klok,'.fl' end if else if ( klok .lt. 10 ) then write(filefl,201) 'hrt' ,expnum,'k000',klok,'.fl' elseif (klok .lt. 100 ) then write(filefl,202) 'hrt' ,expnum,'k00' ,klok,'.fl' elseif (klok .lt. 1000 ) then write(filefl,203) 'hrt' ,expnum,'k0' ,klok,'.fl' elseif (klok .lt. 10000) then write(filefl,204) 'hrt' ,expnum,'k' ,klok,'.fl' else write(filefl,205) 'hrt' ,expnum, klok,'.fl' end if end if if (klokout .eq. 0) then open(9,file=filefl,form='formatted') end if k = 1 do 100 iv=1,4 call emflow(xmk(1,k),nextnm,nextnc,nmarks(iv),ncircs,xcen(1,iv), c rad,xnorm,axnorm,xweb,uweb,aweb,acmp, c flow(iv),0) k = k + nmarks(iv) if (klokout .ne. 0) go to 100 c c write out the velocities and spatial locations on the flowmeter web if (iv.eq.1) write(9,60) klok,filefl,4,(ncircs,nmarks(i),i,i=1,4) 60 format(' VELOCITIES AND LOCATIONS ON FLOWMETER WEB', c /,i5,5x,' = KLOK',1x,a14, c /,i5,5x,' = NUMBER OF VALVES', c4(/,i5,i5,' = NUMBER OF WEB CIRCLES AND RADIAL ARMS IN VALVE',i2)) write(9,71) iv nc = 0 nm = 1 write(9,72) (flnref*uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3),nc do 70 nc=1,ncircs nm = 1 write(9,72) (flnref*uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3),nc do 70 nm=2,nmarks(iv) write(9,72) (flnref*uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3) 70 continue 71 format(i5,' = VALVE NUMBER', c /,6x,'U',5x,6x,'V',5x,6x'W',5x,3x,'X',2x,3x,'Y',2x,3x,'Z') 72 format(3e12.4 ,3f6.2,i3) 100 continue write(6,200) (flnref*flow(iv),iv=1,4) 200 format(' VALVE FLOWS = ',4E15.6) if (klokout .eq. 0) then close(9) write(command,'(a20,i2,a1,a14)') c 'cfs -r5 store hrtin_',expnum,'/',filefl c write(6,'(a37)') command(1:37) c istat = iexec(command) c write(6,*) 'istat after ',command(1:37),' = ',istat end if return end C********************************************************************** subroutine srcflow(isrc,xsrc,nextnm,nextnc,klok,klokout,nref) save icrnr,srclabel c c This subroutine acts as a driver for the flowmeter, c calling subroutine emflow for each of the 6 faces of the cube c surrounding xsrc. c c xcrnr contains the coordinates of the 8 corners of a cube, 12 meshwidths c on a side, surrounding the source located at xsrc. The cube is c located midway between the mesh in all 3 directions. xcrnr must c be recomputed each time-step. c icrnr contains the indices of the corners in the order necessary to c compute the outward flow through each of the 6 faces of the cube. c The first dimension of icrnr is 5 to allow the first index to be c the last index so that the face is 'closed'. icrnr need only be c computed once. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NGM1=NG-1) PARAMETER(NMSIZE=7593) parameter(nmarmx=128,ncirmx=23) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) parameter(nunits=12+nsrcs) dimension xsrc(3) dimension nextnm(*) dimension nextnc(*) c "flowmeter arrays" dimension xsquar(3,nmarmx) dimension rad( 3,nmarmx) dimension xnorm( 3,nmarmx) dimension axnorm( nmarmx) dimension xweb( 3,nmarmx,0:ncirmx) dimension uweb( 3,nmarmx,0:ncirmx) dimension aweb( nmarmx,0:ncirmx) dimension acmp( 3,nmarmx,0:ncirmx) dimension nmarks( 6) dimension xcen( 3,6) dimension flow( 6) dimension xcrnr( 3,8) dimension icrnr( 5,6) character*6 srclabel(5) data icrnr/1,2,3,4,1, 2 5,6,7,8,5, 3 1,5,8,2,1, 4 3,7,6,4,3, 5 1,4,6,5,1, 6 2,8,7,3,2/ data srclabel/'SVCSRC','IVCSRC',' PVSRC',' PASRC',' AoSRC'/ do 8 i1=1,3 i3 = (i1-1)*2 + 1 do 4 i4=1,4 i2 = icrnr(i4,i3) xcrnr(i1,i2) = xsrc(i1) - 6.0 4 continue i3 = (i1-1)*2 + 2 do 6 i4=1,4 i2 = icrnr(i4,i3) xcrnr(i1,i2) = xsrc(i1) + 6.0 6 continue 8 continue ncircs = ncirmx flowsum = 0.0 imkmx = nmarmx/4 do 100 iv=1,6 k = 0 do 20 i4=1,4 i1 = icrnr(i4 ,iv) i2 = icrnr(i4+1,iv) do 10 imk=1,imkmx frac = float(imk)/float(imkmx) xsquar(1,k+imk) = xcrnr(1,i1) + frac*(xcrnr(1,i2) - xcrnr(1,i1)) xsquar(2,k+imk) = xcrnr(2,i1) + frac*(xcrnr(2,i2) - xcrnr(2,i1)) xsquar(3,k+imk) = xcrnr(3,i1) + frac*(xcrnr(3,i2) - xcrnr(3,i1)) 10 continue k = k + imkmx 20 continue nmarks(iv) = k call emflow(xsquar(1,1),nextnm,nextnc,nmarks(iv),ncircs, c xcen(1,iv),rad,xnorm,axnorm,xweb,uweb,aweb,acmp, c flow(iv),1) flowsum = flowsum + flow(iv) c if (klokout .ne. 0) go to 100 c c write out the velocities and spatial locations on the flowmeter web c if (iv .eq. 1) write(9,60) klok,4,(ncircs,nmarks(i),i,i=1,4) c 60 format(' VELOCITIES AND LOCATIONS ON FLOWMETER WEB', c c /,i5,5x,' = KLOK', c c /,i5,5x,' = NUMBER OF VALVES', c c4(/,i5,i5,' = NUMBER OF WEB CIRCLES AND RADIAL ARMS IN VALVE',i2)) c write(9,71) iv c nc = 0 c nm = 1 c write(9,72) (uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3),nc c do 70 nc=1,ncircs c nm = 1 c write(9,72) (uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3),nc c do 70 nm=2,nmarks(iv) c write(9,72) (uweb(i,nm,nc),i=1,3),(xweb(i,nm,nc),i=1,3) c 70 continue c 71 format(i5,' = VALVE NUMBER', c c /,6x,'U',5x,6x,'V',5x,6x'W',5x,3x,'X',2x,3x,'Y',2x,3x,'Z') c 72 format(3e12.4 ,3f6.2,i3) 100 continue flnref = nref write(6,200) srclabel(isrc), c (flnref*flow(iv),iv=1,6,2), c srclabel(isrc), c (flnref*flow(iv),iv=2,6,2), c flnref*flowsum 200 format(A6,' FLOWS = ',3E15.6, c /, A6,' FLOWS = ',3E15.6,E14.6,' =SUM') return end C********************************************************************** subroutine emflow(xmk,nextnm,nextnc,nmarks,ncircs,xcen, c rad,xnorm,axnorm,xweb,uweb,aweb,acmp, c flow,iflag) c c This routine computes the flow through one of the 4 valve rings. c The identity of the ring is not relevant and is not known to the routine. c Each ring is marked by nmarks 'fluid markers' whose coordinates are c passed to the routine in the array xmk. These markers are used to c establish a spider-web-like mesh in the valve ring, and the lattice c velocities are interpolated to the intersections of the webs. The c interpolated velocities and the areas of the web subdivisions are used to c compute the flows through the rings. c c nmarks = number of markers on the valve ring. This is the c number of radial arms in the spider-web. c ncircs = number of circular tracks in each spider-web, c not including the center of the web which is track 0. c xweb = 3-space coordinates of the intersections on the spider-web. c uweb = 3-space velocities at the intersections on the spider-web. c xcen = 3-space coordinates of the "center" of the spider-web. c uweb = 3-space velocity at the intersections on the spider-web. c rad = 3-space vectors along the radial arms of the spider-web. c xnorm = 3-space unit normals to the subdivisions of the spider-web. c aweb = areas of the subdivisions of the spider-web. c axnorm = auxiliary variables used to compute xnorm. c acmp = auxiliary variables used to compute aweb. c iflag = indicator of whether the motion of the meter is to be c subtracted from the motion of the fluid (iflag=0) c or not (iflag=1) c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NMSIZE=7593) PARAMETER(FLNG=NG) COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) dimension xmk(3,*),rad(3,*),xnorm(3,*),axnorm(*) dimension xweb(3,nmarks,0:ncircs),xcen(3) dimension uweb(3,nmarks,0:ncircs) dimension aweb( nmarks,0:ncircs) dimension acmp(3,nmarks,0:ncircs) dimension nextnm(nmarks,0:ncircs) dimension nextnc(nmarks,0:ncircs) DIMENSION FRSTNM(0:NGM1),FRSTNC(0:NGM1),NUMBER(0:NGM1) INTEGER FRSTNM ,FRSTNC DIMENSION XM1OLD(64), XM2OLD(64), XM3OLD(64) DIMENSION IZ(64), JZ(64), KZ(64) DIMENSION IMOD(0:4,64), JMOD(0:4,64), KMOD(0:4,64) DIMENSION D1(0:4,64), D2(0:4,64), D3(0:4,64) DIMENSION DELTA(0:64,64),INDEX(0:64,64) DIMENSION ULIN(0:63), VLIN(0:63), WLIN(0:63) DIMENSION FLIZP1(64), FLJZP1(64), FLKZP1(64) C C MARKER POINT L IS INFLUENCED BY FLUID POINT (I,J,K) IFF C ALL THE FOLLOWING CONDITIONS ARE SATISFIED: C C ABS(I - XMK(1,L)) < 2. C ABS(J - XMK(2,L)) < 2. C ABS(K - XMK(3,L)) < 2. C C USING FORTRAN, THESE VALUES OF (I,J,K) ARE FOUND AS FOLLOWS: C C IZ = XMK(1,L) -1. C JZ = XMK(2,L) -1. C KZ = XMK(3,L) -1. C C THEN C I = IZ,IZ+3 C J = JZ,JZ+3 C K = KZ,KZ+3 C C IF THESE CONDITIONS ARE SATISFIED, THE COEFFICIENT THAT LINKS C FIBER POINT L TO FLUID POINT (I,J,K) IS: C C D = DEL(I-XMK(1,L)) * DEL(J-XMK(2,L)) * DEL(K-XMK(3,L)) C C WHERE C C DEL(R) = (1. + COS((PI/2.)*R))/4. C C ALGORITHM IS AS FOLLOWS: C C LET L=1,2, ... ,NP C FOR EACH L, LET (I,J,K) COVER RANGE OF VALUES DEFINED ABOVE. C C FOR EACH VALUE OF (L,I,J,K) C C XMK(1,L) = XMK(1,L) + D * U(I,J,K) C XMK(2,L) = XMK(2,L) + D * V(I,J,K) C XMK(3,L) = XMK(3,L) + D * W(I,J,K) C C WHERE D WAS DEFINED ABOVE. C MODNG(K)=MOD(K+NG,NG) DEL(R) = (1. + COS((PI/2.)*R))/4. PI = 4. * ATAN(1.) c c compute the center of the spider-web c xc = 0.0 yc = 0.0 zc = 0.0 do 1 n=1,nmarks xc = xc + xmk(1,n) yc = yc + xmk(2,n) zc = zc + xmk(3,n) 1 continue xcen(1) = xc/float(nmarks) xcen(2) = yc/float(nmarks) xcen(3) = zc/float(nmarks) c c compute the coordinates of the spider-web mesh points c c points on the center of the web (track 0) nc = 0 do 3 nm=1,nmarks xweb(1,nm,nc) = xcen(1) xweb(2,nm,nc) = xcen(2) xweb(3,nm,nc) = xcen(3) 3 continue c points on the perimeter of the web (track ncircs) nc = ncircs do 4 nm=1,nmarks xweb(1,nm,nc) = xmk(1,nm) xweb(2,nm,nc) = xmk(2,nm) xweb(3,nm,nc) = xmk(3,nm) 4 continue c points on all the intermediate circular tracks do 5 nc=1,ncircs-1 wt1 = float( nc)/float(ncircs) wt2 = float(ncircs-nc)/float(ncircs) do 5 nm=1,nmarks xweb(1,nm,nc) = wt1*xweb(1,nm,ncircs) + wt2*xweb(1,nm,0) xweb(2,nm,nc) = wt1*xweb(2,nm,ncircs) + wt2*xweb(2,nm,0) xweb(3,nm,nc) = wt1*xweb(3,nm,ncircs) + wt2*xweb(3,nm,0) 5 continue C C INITIALIZE THE LINKED LIST VARIABLES C DO 12 KK=0,NGM1 FRSTNM(KK) = 0 FRSTNC(KK) = -1 NUMBER(KK) = 0 12 CONTINUE C C SORT THE xweb DATA BY Z-COORDINATE USING LINKED LISTS C do 13 nc=0,ncircs do 13 nm=1,nmarks C KK = xweb(3,nm,nc) KK = MODNG(INT(xweb(3,nm,nc)+NG)) NEXTNM(NM,NC) = FRSTNM(KK) NEXTNC(NM,NC) = FRSTNC(KK) FRSTNM(KK) = nm FRSTNC(KK) = nc NUMBER(KK) = NUMBER(KK)+1 13 CONTINUE c c Establish the minimum and maximum values of kk. This will minimize the c number of planes of velocity data which are read later. c do 14 kk=0,ngm1 kkmin = kk if (number(kk) .gt. 0) go to 16 14 continue c c Start loop at kkmin to cover the likely possibility that kkmax=kkmin. c Because the previous loop established that number(kkmin) > 0, c the following loop is guaranteed to assign a value to kkmax c 16 do 17 kk=kkmin,ngm1 if (number(kk) .eq. 0) go to 18 kkmax = kk 17 continue 18 continue C C FRSTNM(K) CONTAINS THE NM INDEX AND FRSTNC(K) CONTAINS THE NC INDEX C OF THE FIRST POINT BETWEEN PLANES K AND K+1, WHICH HAS COORDINATES C C (XWEB(1,FRSTNM(K),FRSTNC(K)), C XWEB(2,FRSTNM(K),FRSTNC(K)), C XWEB(3,FRSTNM(K),FRSTNC(K))). C C NEXTNM(FRSTNM(K),FRSTNC(K)) CONTAINS THE NM INDEX OF C THE SECOND SUCH POINT, C NEXTNC(FRSTNM(K),FRSTNC(K)) CONTAINS THE NC INDEX OF C THE SECOND SUCH POINT, C C NEXTNM(NEXTNM(FRSTNM(K),FRSTNC(K)), CONTAINS THE NM INDEX OF C NEXTNC(FRSTNM(K),FRSTNC(K))) THE THIRD SUCH POINT, C NEXTNC(NEXTNM(FRSTNM(K),FRSTNC(K)), CONTAINS THE NC INDEX OF C NEXTNC(FRSTNM(K),FRSTNC(K))) THE THIRD SUCH POINT, C C ETC. C DO 120 KK=kkmin,kkmax KXM1=MODNG(KK-1) KX =MODNG(KK ) KXP1=MODNG(KK+1) KXP2=MODNG(KK+2) KXP3=MODNG(KK+3) NM = FRSTNM(KK) NC = FRSTNC(KK) IF (NM .EQ. 0) GO TO 120 105 IF (NUMBER(KK) .GE. 64) THEN NPOINTS = 64 NUMBER(KK) = NUMBER(KK) - 64 ELSE NPOINTS = NUMBER(KK) NUMBER(KK) = 0 END IF NM0 = NM NC0 = NC DO 106 NPT=1,NPOINTS XM1OLD(NPT) = XWEB(1,NM,NC) XM2OLD(NPT) = XWEB(2,NM,NC) XM3OLD(NPT) = XWEB(3,NM,NC) IZ(NPT) = INT(XM1OLD(NPT) -1. + FLNG) - NG JZ(NPT) = INT(XM2OLD(NPT) -1. + FLNG) - NG KZ(NPT) = INT(XM3OLD(NPT) -1. + FLNG) - NG FLIZP1(NPT) = IZ(NPT)+1 FLJZP1(NPT) = JZ(NPT)+1 FLKZP1(NPT) = KZ(NPT)+1 NMTEMP = NM NM = NEXTNM(NMTEMP,NC) NC = NEXTNC(NMTEMP,NC) 106 CONTINUE DO 107 NPT=1,NPOINTS ARG3 = XM3OLD(NPT) - FLKZP1(NPT) ARG2 = XM2OLD(NPT) - FLJZP1(NPT) ARG1 = XM1OLD(NPT) - FLIZP1(NPT) RAD3 = SQRT(1.+4.*ARG3*(1.-ARG3)) RAD2 = SQRT(1.+4.*ARG2*(1.-ARG2)) RAD1 = SQRT(1.+4.*ARG1*(1.-ARG1)) D3(3,NPT) = (1.+2.*ARG3-RAD3)/8. D2(3,NPT) = (1.+2.*ARG2-RAD2)/8. D1(3,NPT) = (1.+2.*ARG1-RAD1)/8. D3(2,NPT) = (1.+2.*ARG3+RAD3)/8. D2(2,NPT) = (1.+2.*ARG2+RAD2)/8. D1(2,NPT) = (1.+2.*ARG1+RAD1)/8. D3(1,NPT) = (3.-2.*ARG3+RAD3)/8. D2(1,NPT) = (3.-2.*ARG2+RAD2)/8. D1(1,NPT) = (3.-2.*ARG1+RAD1)/8. D3(0,NPT) = (3.-2.*ARG3-RAD3)/8. D2(0,NPT) = (3.-2.*ARG2-RAD2)/8. D1(0,NPT) = (3.-2.*ARG1-RAD1)/8. 107 CONTINUE DO 108 M=0,3 DO 108 NPT=1,NPOINTS KMOD(M,NPT) = MODNG(KZ(NPT)+M) JMOD(M,NPT) = MODNG(JZ(NPT)+M-1)+1 IMOD(M,NPT) = MODNG(IZ(NPT)+M-1)+1 108 CONTINUE DO 110 K=0,3 DO 110 J=0,3 DO 110 I=0,3 M = K*16 + J*4 + I DO 110 NPT=1,NPOINTS DELTA(M,NPT) = D1(I,NPT) * D2(J,NPT) * D3(K,NPT) INDEX(M,NPT) = KMOD(K,NPT) *(NB+1)*NG C + (JMOD(J,NPT)-1)*(NB+1) C + IMOD(I,NPT) + 1 110 CONTINUE NM = NM0 NC = NC0 DO 112 NPT=1,NPOINTS CALL GATHER(64,ULIN(0),U(0,1,0),INDEX(0,NPT)) CALL GATHER(64,VLIN(0),V(0,1,0),INDEX(0,NPT)) CALL GATHER(64,WLIN(0),W(0,1,0),INDEX(0,NPT)) UWEB(1,NM,NC) = SDOT(64,DELTA(0,NPT),1,ULIN(0),1) UWEB(2,NM,NC) = SDOT(64,DELTA(0,NPT),1,VLIN(0),1) UWEB(3,NM,NC) = SDOT(64,DELTA(0,NPT),1,WLIN(0),1) NMTEMP = NM NM = NEXTNM(NMTEMP,NC) NC = NEXTNC(NMTEMP,NC) 112 CONTINUE IF (NUMBER(KK) .NE. 0) GO TO 105 120 CONTINUE c c UWEB now contains the velocity interpolated to the points of the spider-web. c UWEB( ,ncircs,nm),nm=1,nmarks contains the velocity at the locations of the c fluid markers, i.e., the velocity of the "flowmeter". c if (iflag .ne. 0) go to 399 c c compute the average 3-space velocity of the markers c umkavg = 0.0 vmkavg = 0.0 wmkavg = 0.0 do 200 nm=1,nmarks umkavg = umkavg + uweb(1,nm,ncircs) vmkavg = vmkavg + uweb(2,nm,ncircs) wmkavg = wmkavg + uweb(3,nm,ncircs) 200 continue umkavg = umkavg/float(nmarks) vmkavg = vmkavg/float(nmarks) wmkavg = wmkavg/float(nmarks) c c subtract the linearly-interpolated velocity of the flowmeter from the c previously computed velocities on the web. This gives the velocities c relative to the web. c do 300 nc=0,ncircs-1 wt1 = float( nc)/float(ncircs) wt2 = float(ncircs-nc)/float(ncircs) do 300 nm=1,nmarks uweb(1,nm,nc) = uweb(1,nm,nc) - wt1*uweb(1,nm,ncircs) - wt2*umkavg uweb(2,nm,nc) = uweb(2,nm,nc) - wt1*uweb(2,nm,ncircs) - wt2*vmkavg uweb(3,nm,nc) = uweb(3,nm,nc) - wt1*uweb(3,nm,ncircs) - wt2*wmkavg 300 continue do 310 nm=1,nmarks uweb(1,nm,ncircs) = 0.0 uweb(2,nm,ncircs) = 0.0 uweb(3,nm,ncircs) = 0.0 310 continue c c UWEB now contains the velocities at the web intersections relative to the web. c 399 continue c c Compute the unit normals (xnorm()) to the triangular areas defined by the c average marker location (xcen()) and two neighboring markers. Store the c normals in an array indexed by the smaller of the marker indices (except for c the last which uses markers index=nmarks and index=1 and is itself indexed c nmarks). These are also the unit normals to any sub-area of the triangles such c as the trapezoidal areas contained within neighboring radial arms and c neighboring circular tracks of the web. c c radial vectors do 400 nm=1,nmarks rad(1,nm) = xmk(1,nm) - xcen(1) rad(2,nm) = xmk(2,nm) - xcen(2) rad(3,nm) = xmk(3,nm) - xcen(3) 400 continue c c normal vectors via cross product do 410 nm=1,nmarks-1 xnorm(3,nm) = rad(1,nm)*rad(2,nm+1) - rad(2,nm)*rad(1,nm+1) xnorm(1,nm) = rad(2,nm)*rad(3,nm+1) - rad(3,nm)*rad(2,nm+1) xnorm(2,nm) = rad(3,nm)*rad(1,nm+1) - rad(1,nm)*rad(3,nm+1) 410 continue nm = nmarks nmp1 = 1 xnorm(3,nm) = rad(1,nm)*rad(2,nmp1) - rad(2,nm)*rad(1,nmp1) xnorm(1,nm) = rad(2,nm)*rad(3,nmp1) - rad(3,nm)*rad(2,nmp1) xnorm(2,nm) = rad(3,nm)*rad(1,nmp1) - rad(1,nm)*rad(3,nmp1) c c absolute values (i.e., lengths) of normal vectors do 420 nm=1,nmarks axnorm(nm) = sqrt(xnorm(1,nm)**2 +xnorm(2,nm)**2 +xnorm(3,nm)**2) 420 continue c c unit normal vectors do 430 nm=1,nmarks xnorm(1,nm) = xnorm(1,nm)/axnorm(nm) xnorm(2,nm) = xnorm(2,nm)/axnorm(nm) xnorm(3,nm) = xnorm(3,nm)/axnorm(nm) 430 continue c c Compute the subareas of the web. c Each subarea is a trapezoid (the subareas at the center of the web are c triangles which are simply trapezoids with tops of length 0). c Let A be a vector from left to right along the bottom. c Let C be a vector from left to right along the top. c Let B be a vector from the tail end of A to the tail end of C. c Let cross denote the vector cross product operator. c The area of the trapezoid is given by one-half the magnitude of the vector: c c (A + C) cross B c c (Note that A and C are parallel, thus | A | + | C | = | A + C |) c The array acmp( ,nm,nc) contains the 3 components of this vector for c each of the nm*nc trapezoids. ("acmp" stands for "area components".) c do 500 nc=0,ncircs-1 do 500 nm=1,nmarks-1 acmp(3,nm,nc) = ( xweb(1,nm+1,nc ) - xweb(1,nm,nc ) c + xweb(1,nm+1,nc+1) - xweb(1,nm,nc+1) ) c *( xweb(2,nm ,nc ) - xweb(2,nm,nc+1) ) c - ( xweb(2,nm+1,nc ) - xweb(2,nm,nc ) c + xweb(2,nm+1,nc+1) - xweb(2,nm,nc+1) ) c *( xweb(1,nm ,nc ) - xweb(1,nm,nc+1) ) acmp(1,nm,nc) = ( xweb(2,nm+1,nc ) - xweb(2,nm,nc ) c + xweb(2,nm+1,nc+1) - xweb(2,nm,nc+1) ) c *( xweb(3,nm ,nc ) - xweb(3,nm,nc+1) ) c - ( xweb(3,nm+1,nc ) - xweb(3,nm,nc ) c + xweb(3,nm+1,nc+1) - xweb(3,nm,nc+1) ) c *( xweb(2,nm ,nc ) - xweb(2,nm,nc+1) ) acmp(2,nm,nc) = ( xweb(3,nm+1,nc ) - xweb(3,nm,nc ) c + xweb(3,nm+1,nc+1) - xweb(3,nm,nc+1) ) c *( xweb(1,nm ,nc ) - xweb(1,nm,nc+1) ) c - ( xweb(1,nm+1,nc ) - xweb(1,nm,nc ) c + xweb(1,nm+1,nc+1) - xweb(1,nm,nc+1) ) c *( xweb(3,nm ,nc ) - xweb(3,nm,nc+1) ) 500 continue nm = nmarks nmp1 = 1 do 510 nc=0,ncircs-1 acmp(3,nm,nc) = ( xweb(1,nmp1,nc ) - xweb(1,nm,nc ) c + xweb(1,nmp1,nc+1) - xweb(1,nm,nc+1) ) c *( xweb(2,nm ,nc ) - xweb(2,nm,nc+1) ) c - ( xweb(2,nmp1,nc ) - xweb(2,nm,nc ) c + xweb(2,nmp1,nc+1) - xweb(2,nm,nc+1) ) c *( xweb(1,nm ,nc ) - xweb(1,nm,nc+1) ) acmp(1,nm,nc) = ( xweb(2,nmp1,nc ) - xweb(2,nm,nc ) c + xweb(2,nmp1,nc+1) - xweb(2,nm,nc+1) ) c *( xweb(3,nm ,nc ) - xweb(3,nm,nc+1) ) c - ( xweb(3,nmp1,nc ) - xweb(3,nm,nc ) c + xweb(3,nmp1,nc+1) - xweb(3,nm,nc+1) ) c *( xweb(2,nm ,nc ) - xweb(2,nm,nc+1) ) acmp(2,nm,nc) = ( xweb(3,nmp1,nc ) - xweb(3,nm,nc ) c + xweb(3,nmp1,nc+1) - xweb(3,nm,nc+1) ) c *( xweb(1,nm ,nc ) - xweb(1,nm,nc+1) ) c - ( xweb(1,nmp1,nc ) - xweb(1,nm,nc ) c + xweb(1,nmp1,nc+1) - xweb(1,nm,nc+1) ) c *( xweb(3,nm ,nc ) - xweb(3,nm,nc+1) ) 510 continue do 520 nc=0,ncircs-1 do 520 nm=1,nmarks aweb(nm,nc) = 0.5*sqrt( acmp(1,nm,nc)**2 c + acmp(2,nm,nc)**2 c + acmp(3,nm,nc)**2) 520 continue c c Compute the flow through the web. c Let area denote the area of trapezoid (nm,nc) c Let vel denote the avg. of the velocity vectors at that trapezoid's corners. c Let norm denote the unit normal vector to that trapezoid. c Let dot denote the vector dot product operator. c Flow is the sum over all the trapezoids of: c c area*(vel dot norm) c flow = 0.0 do 600 nc=0,ncircs-1 do 600 nm=1,nmarks-1 flow = flow + aweb(nm,nc)*( ( uweb(1,nm,nc )+uweb(1,nm+1,nc ) c +uweb(1,nm,nc+1)+uweb(1,nm+1,nc+1)) c * xnorm(1,nm) c + ( uweb(2,nm,nc )+uweb(2,nm+1,nc ) c +uweb(2,nm,nc+1)+uweb(2,nm+1,nc+1)) c * xnorm(2,nm) c + ( uweb(3,nm,nc )+uweb(3,nm+1,nc ) c +uweb(3,nm,nc+1)+uweb(3,nm+1,nc+1)) c * xnorm(3,nm) ) 600 continue nm = nmarks nmp1 = 1 do 610 nc=0,ncircs-1 flow = flow + aweb(nm,nc)*( ( uweb(1,nm,nc )+uweb(1,nmp1,nc ) c +uweb(1,nm,nc+1)+uweb(1,nmp1,nc+1)) c * xnorm(1,nm) c + ( uweb(2,nm,nc )+uweb(2,nmp1,nc ) c +uweb(2,nm,nc+1)+uweb(2,nmp1,nc+1)) c * xnorm(2,nm) c + ( uweb(3,nm,nc )+uweb(3,nmp1,nc ) c +uweb(3,nm,nc+1)+uweb(3,nmp1,nc+1)) c * xnorm(3,nm) ) 610 continue flow = 0.25*flow RETURN END C********************************************************************** SUBROUTINE FIBERX(XF0,STF0,REST0,FRC,LAFLAG,ARFLAG, C XFN,STF ,REST ,UNSTBL, C NFG,NPF,NGROUPS) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) PARAMETER(IMAX=63,NBUNCH=1) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) parameter(nunits=12+nsrcs) C C A FIBER IS COMPOSED OF POINTS C A GROUP IS COMPOSED OF FIBERS HAVING THE SAME NUMBER OF POINTS C A BUNCH IS COMPOSED OF GROUPS C C NBUNCH =NUMBER OF BUNCHES IN THE ENTIRE STRUCTURE C NGROUPS(J)=NUMBER OF GROUPS IN BUNCH J, J=1,...,NBUNCH (NGROUPS(0)=0) C NFG(I) =NUMBER OF FIBERS IN GROUP I, I=NGROUPS(J-1)+1,...,NGROUPS(J) C NPF(I) =NUMBER OF POINTS IN A FIBER IN GROUP I C DIMENSION NGROUPS(NBUNCH) C C IN THE FOLLOWING COMMENTS, LET IMAX= SUM(J=1,NBUNCH):NGROUPS(J). C IMAX IS THE TOTAL NUMBER OF GROUPS IN THE ENTIRE STRUCTURE. C C IN THE FOLLOWING ARRAYS, THE LAST DIMENSION, WHICH IS UNSPECIFIED C HERE, SHOULD BE IMAX C DIMENSION NFG(*),NPF(*) C C IN THE FOLLOWING ARRAYS WHICH HOLD INPUT OR OUPUT TO C THIS SUBROUTINE, THE LAST DIMENSION, NOT SPECIFIED HERE, C SHOULD BE NFSIZE C WHERE NFSIZE = MAX(J=1,NBUNCH):SUM(I=ISTART(J),ISTOP(J)):NFG(I)*NPF(I) C ISTART(J)=NGROUPS(J-1)+1 C ISTOP (J)=NGROUPS(J) C NGROUPS(0)=0 C DIMENSION XF0 (3,*) DIMENSION FRC (3,*) DIMENSION STF0 (*) DIMENSION REST0 (*) DIMENSION LAFLAG(*) LOGICAL ARFLAG(*) DIMENSION XFN (3,*) DIMENSION STF ( *) DIMENSION REST( *) LOGICAL UNSTBL(*) DIMENSION KGROUP(IMAX) C C THE CALLED ROUTINES WORK WITH ALL OF THE FIBERS (A GROUP AT A TIME). C FIBFRCZ COMPUTES FRC(XF). C C THE POINTS OF FIBER I ARE DESIGNATED BY THE INDICES C FROM NPFSUM(I-1)+1 THROUGH NPFSUM(I) C WHERE NPFSUM(I)=SUM(I=1,N):NFG(I)*NPF(I) C C IN THE FOLLOWING COMMENTS K RUNS THROUGH ALL OF THE INDICES C OF THE FIBER POINTS. THAT IS K=1,NPFSUM(NF): C INPUT VARIABLES: C NF=NUMBER OF FIBERS C NPF(I)=NUMBER OF POINTS IN FIBER I, I=1,NF C XF( ,K)= INITIAL GUESS FOR POSITION OF FIBER POINT K C STF(K)= STIFFNESS OF FIBER LINK K C REST(K)= RESTING LENGTH OF FIBER LINK K C FL= GIVEN PARAMETER C ITMAX= MAXIMUM NUMBER OF ITERATIONS IN FIBNWT C ESTOP= STOPPING CRITERION FOR FIBNWT C OUTPUT VARIABLES: C XF( ,K)=POSITION OF FIBER POINT K AS DETERMINED BY FIBNWT C FRC( ,K)=FORCE APPLIED TO THE FLUID BY FIBER POINT K C ISTART = 1 ISTOP = NGROUPS(1) K = 1 DO 10 I=ISTART,ISTOP KGROUP(I) = K K = K + NFG(I)*NPF(I) 10 CONTINUE NMLNKS = K - 1 CMIC$ DO ALL CMIC$1SHARED(XFN, STF0, REST0, FRC, LAFLAG, ARFLAG, NFG, NPF, CMIC$2 STF , REST , UNSTBL, CMIC$3 ISTART, ISTOP, KGROUP) CMIC$4PRIVATE(I, K, KOUNT) DO 1 I=ISTART,ISTOP K = KGROUP(I) CALL MUSCLE( STF0(K), STF(K), C REST0(K),REST(K),laflag(K),arflag(k),NFG(I),NPF(I)) CALL FIBFRCZ(XFN(1,K),STF(K),REST(K),NFG(I),NPF(I),FRC(1,K), C UNSTBL(K) ) 1 CONTINUE CMIC$ DO ALL CMIC$1SHARED(XFN, STF0, REST0, FRC, LAFLAG, ARFLAG, NFG, NPF, CMIC$2 XF0, STF , REST , CMIC$3 ISTART, ISTOP, KGROUP, KGRPS) CMIC$4PRIVATE(I, K, K0, KOUNT) DO 3 I=25,26 K = KGROUP(I) K0 = K - KGROUP(25) + 1 CALL ANCHORZ(XF0(1,K0), C XFN(1,K),STF(K),REST(K),NFG(I),NPF(I),FRC(1,K)) 3 CONTINUE KOUNT = ILSUM(NMLNKS,UNSTBL,1) IF (KOUNT .GT. 0) THEN DO 2 I=ISTART,ISTOP K = KGROUP(I) KOUNT = ILSUM(NFG(I)*NPF(I),UNSTBL(K),1) WRITE(6,*) 'THE NUMBER OF LONG LINKS IN GROUP ',I,' = ',KOUNT 2 CONTINUE CALL EXIT(1) END IF RETURN END C********************************************************************** SUBROUTINE FIBERS(FL , XF , XFSAVE, 1 XFZ, XFZOLD, XFZNEW, 2 STF, STFOLD, STFNEW, STF0, 3 REST,RESTOLD,RESTNEW,REST0,LAFLAG,NFG,NPF, 4 B,C,X,XM, 5 CSAVE,XCON,D,T1,T2, 6 ITMAX,ESTOP,MRAMP,MFLAT,FRC,NGROUPS,NBUNCH) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NFGMAX=64,NPFMAX=530,NPFGMX=NFGMAX*NPFMAX) PARAMETER(NFSIZE=606638) parameter(nsrcs=5,ntapx=33,nsrcspx=nsrcs+ntapx) parameter(nunits=12+nsrcs) C C A FIBER IS COMPOSED OF POINTS C A GROUP IS COMPOSED OF FIBERS HAVING THE SAME NUMBER OF POINTS C A BUNCH IS COMPOSED OF GROUPS C C NBUNCH =NUMBER OF BUNCHES IN THE ENTIRE STRUCTURE C NGROUPS(J)=NUMBER OF GROUPS IN BUNCH J, J=1,...,NBUNCH (NGROUPS(0)=0) C NFG(I) =NUMBER OF FIBERS IN GROUP I, I=NGROUPS(J-1)+1,...,NGROUPS(J) C NPF(I) =NUMBER OF POINTS IN A FIBER IN GROUP I C DIMENSION NGROUPS(NBUNCH) C C IN THE FOLLOWING COMMENTS, LET IMAX= SUM(J=1,NBUNCH):NGROUPS(J). C IMAX IS THE TOTAL NUMBER OF GROUPS IN THE ENTIRE STRUCTURE. C C IN THE FOLLOWING ARRAYS, THE LAST DIMENSION, WHICH IS UNSPECIFIED C HERE, SHOULD BE IMAX C DIMENSION NFG(*),NPF(*) DIMENSION MRAMP(*) C C IN THE FOLLOWING (WORKSPACE) ARRAYS, THE LAST DIMENSION, WHICH IS C UNSPECIFIED HERE, SHOULD BE MAX(I=1,IMAX):NFG(I)*NPF(I) C DIMENSION B(3,3,*),C(3,3,*),X(3,*),XM(3,3,*) C C IN THE FOLLOWING (WORKSPACE) ARRAYS, THE LAST DIMENSION, WHICH IS C UNSPECIFIED HERE, SHOULD BE MAX(I=1,IMAX):NFG(I) C DIMENSION CSAVE(3,3,*),XCON(3,*),D(3,*),T1(*),T2(*) C C IN THE FOLLOWING ARRAYS WHICH HOLD INPUT OR OUPUT TO C THIS SUBROUTINE, THE LAST DIMENSION, NOT SPECIFIED HERE, C SHOULD BE NFSIZE C WHERE NFSIZE = MAX(J=1,NBUNCH):SUM(I=ISTART(J),ISTOP(J)):NFG(I)*NPF(I) C ISTART(J)=NGROUPS(J-1)+1 C ISTOP (J)=NGROUPS(J) C NGROUPS(0)=0 C DIMENSION XF(3,*),XFZ(3,*),STF(*),REST(*),FRC(3,*) DIMENSION FL(*) DIMENSION XFZNEW(3,*) DIMENSION XFZOLD(3,*) DIMENSION STFNEW(*) DIMENSION STFOLD(*) DIMENSION STF0 (*) DIMENSION RESTNEW(*) DIMENSION RESTOLD(*) DIMENSION REST0 (*) DIMENSION LAFLAG(*) C C THE CALLED ROUTINES WORK WITH ALL OF THE FIBERS (A GROUP AT A TIME). C FOR EACH FIBER, FIBNWT SOLVES FOR A FIBER CONFIGURATION C XF THAT SATISFIES C XF = XFZ + FL*FRC(XF). C THEN FIBFRCZ COMPUTES FRC(XF). C C THE POINTS OF FIBER I ARE DESIGNATED BY THE INDICES C FROM NPFSUM(I-1)+1 THROUGH NPFSUM(I) C WHERE NPFSUM(I)=SUM(I=1,N):NFG(I)*NPF(I) C C IN THE FOLLOWING COMMENTS K RUNS THROUGH ALL OF THE INDICES C OF THE FIBER POINTS. THAT IS K=1,NPFSUM(NF): C INPUT VARIABLES: C NF=NUMBER OF FIBERS C NPF(I)=NUMBER OF POINTS IN FIBER I, I=1,NF C XF( ,K)= INITIAL GUESS FOR POSITION OF FIBER POINT K C XFZ( ,K)= REFERENCE POSITION OF FIBER POINT K C STF(K)= STIFFNESS OF FIBER LINK K C REST(K)= RESTING LENGTH OF FIBER LINK K C FL= GIVEN PARAMETER C ITMAX= MAXIMUM NUMBER OF ITERATIONS IN FIBNWT C ESTOP= STOPPING CRITERION FOR FIBNWT C OUTPUT VARIABLES: C XF( ,K)=POSITION OF FIBER POINT K AS DETERMINED BY FIBNWT C FRC( ,K)=FORCE APPLIED TO THE FLUID BY FIBER POINT K C ISTART = 1 ISTOP = NGROUPS(1) c rewind(43) c rewind(44) K = 1 DO 1 I=ISTART,ISTOP write(6,*)' GROUP=',I CALL MUSCLE( STF0(K), STFOLD, STFNEW(K), C REST0(K),RESTOLD,RESTNEW(K),laflag(K),NFG(I),NPF(I)) CALL FIBNWT(FL(K), XF( 1,K), XFSAVE, 1 XFZ , XFZOLD(1,K), XFZNEW(1,K), 2 STF , STFOLD , STFNEW( K), 3 REST ,RESTOLD ,RESTNEW( K),NFG(I),NPF(I), 4 B,C,X,XM,ITMAX,ESTOP,MRAMP(I),MFLAT, 5 CSAVE,XCON,D,T1,T2) C CALL FIBFRCZ(XF(1,K),STFNEW(K),RESTNEW(K),NFG(I),NPF(I),FRC(1,K)) CALL FIBFRCZ(XF(1,K),STF ,REST ,NFG(I),NPF(I),FRC(1,K)) K=K+NFG(I)*NPF(I) 1 CONTINUE RETURN END C********************************************************************** SUBROUTINE FIBFRC(FL,XF,XFZ,STF,REST,NFG,NPF,B,C,X, C XCON,D,T1,T2) C C THIS ROUTINE COMPUTES FIBER FORCES AND DERIVATIVES C OF FIBER FORCES. C IT INITIALIZES FIBTRI TO SOLVE THE FOLLOWING C EQUATIONS FOR DELTAXF: C C (I-FL*(DFRC/DXF))*DELTAXF= XFZ+FL*FRC-XF C C WHERE I=IDENTITY. C C INPUT VARIABLES: C C FL= GIVEN PARAMETER C XF(K)= POSITION OF THE K-TH BOUNDARY POINT(CURRENT GUESS) C XFZ(K)= REFERENCE POSITION OF THE K-TH BOUNDARY POINT C STF(K)= STIFFNESS OF LINK K C REST(K)= RESTING LENGTH OF LINK K C NPF= NUMBER OF POINTS IN FIBER(K=1,2,...,NPF) C NFG= NUMBER OF FIBERS IN A GROUP. C EACH FIBER IN A GROUP HAS THE SAME NUMBER OF POINTS. C FIBERS IN A DIFFERENT GROUP MAY HAVE A DIFFERENT C NUMBER OF POINTS. C C NOTE THAT LINK K JOINS POINT K AND POINT K+1. C FIBER IS PERIODIC: LINK NPF JOINS POINT NPF AND POINT 1. C C OUTPUT VARIABLES: (LET J,J1,J2=1,2,3) C C B(J1,J2,K)= DELTA(J1,J2) -FL*DFRC(J1,K)/DXF(J2,K) C C(J1,J2,K)= FL*DFRC(J1,K)/DXF(J2,K+1) C X(J,K)= XFZ(J,K) -XF(J,K) +FL*FRC(J,K) C C C IN THE FOREGOING, THE QUANTITIES FRC AND DFRC/DXF C NEED NOT BE STORED.THE ARRAYS C,X ARE COMPUTED BY C SUMMING OVER LINKS AS FOLLOWS: C C INITIALIZATION : X= XFZ-XF C C FOR K=1,2,...NPF, EVALUATE THE CONTRIBUTION OF LINK K C AS FOLLOWS: C C LET R= DIST(XF( ,K),XF( ,K+1)) C IF R< REST(K), IGNORE LINK K.OTHERWISE C LET T= STF(K)*((R-REST(K))**2)/R C = STF(K)*R*(1.-REST(K)/R)**2 C dT/dR= STF(K)*(1.-(REST(K)/R)**2) C XCON(J)= FL*T*(XF(J,K+1)-XF(J,K))/R C X(J,K)= X(J,K) +XCON(J) C X(J,K+1)= X(J,K+1) -XCON(J) C T1= T/R C T2= dT/dR-T1 C D(J)= (XF(J,K+1)-XF(J,K))/R C C C(J1,J2,K)= FL*(T1*DELTA(J1,J2) +T2*D(J1)*D(J2)) C C THIS COMPLETES THE LOOP OVER K. C C ONCE C IS DEFINED, B IS COMPUTED AS FOLLOWS: C B(J1,J2,K)= DELTA(J1,J2) +C(J1,J2,K) +C(J1,J2,K-1) C C C C NOTE: FUNCTION CVMGT(X1,X2,FLAG) C LOGICAL FLAG C IF (FLAG) THEN C CVMGT=X1 C ELSE C CVMGT=X2 C END IF C RETURN C END C DIMENSION XF(3,NFG,NPF),XFZ(3,NFG,NPF) DIMENSION STF(NFG,NPF),REST(NFG,NPF) DIMENSION FL(NFG,NPF) DIMENSION B(3,3,NFG,NPF),C(3,3,NFG,NPF),X(3,NFG,NPF) DIMENSION D(3,NFG),XCON(3,NFG) DIMENSION T1(NFG),T2(NFG) LOGICAL RFLAG DO 1 K=1,NPF DO 1 N=1,NFG X(1,N,K)= (XFZ(1,N,K)-XF(1,N,K))/FL(N,K) X(2,N,K)= (XFZ(2,N,K)-XF(2,N,K))/FL(N,K) X(3,N,K)= (XFZ(3,N,K)-XF(3,N,K))/FL(N,K) 1 CONTINUE DO 10 K=1,NPF KP1= K+1 IF(K.EQ.NPF) KP1=1 DO 4 N=1,NFG R= SQRT((XF(1,N,K)-XF(1,N,KP1))**2 C +(XF(2,N,K)-XF(2,N,KP1))**2 C +(XF(3,N,K)-XF(3,N,KP1))**2) clin RATIO = REST(N,K)/R RATIO = R/REST(N,K) - 1.0 cexp ERATIO= EXP(RATIO) RFLAG= R.LT.REST(N,K) clin T =CVMGT(0., STF(N,K)*R*(1.-RATIO)**2, RFLAG) cexp T = CVMGT(0., 2.*STF(N,K)*REST(N,K)*(ERATIO-(1.+RATIO)), RFLAG) T = CVMGT(0., STF(N,K)*REST(N,K)*RATIO**2, RFLAG) T1(N)= T/R clin T2(N)=CVMGT(0., STF(N,K)*(1.-RATIO**2)-T/R, RFLAG) cexp T2(N)= CVMGT(0., 2.*STF(N,K)*(ERATIO-1.)-T/R, RFLAG) T2(N)= CVMGT(0., 2.*STF(N,K)*RATIO-T/R, RFLAG) XCON(1,N)= T*(XF(1,N,KP1)-XF(1,N,K))/R XCON(2,N)= T*(XF(2,N,KP1)-XF(2,N,K))/R XCON(3,N)= T*(XF(3,N,KP1)-XF(3,N,K))/R X(1,N,K )= X(1,N,K ) +XCON(1,N) X(2,N,K )= X(2,N,K ) +XCON(2,N) X(3,N,K )= X(3,N,K ) +XCON(3,N) X(1,N,KP1)= X(1,N,KP1) -XCON(1,N) X(2,N,KP1)= X(2,N,KP1) -XCON(2,N) X(3,N,KP1)= X(3,N,KP1) -XCON(3,N) D(1,N)= (XF(1,N,KP1) -XF(1,N,K))/R D(2,N)= (XF(2,N,KP1) -XF(2,N,K))/R D(3,N)= (XF(3,N,KP1) -XF(3,N,K))/R 4 CONTINUE DO 2 J1=1,3 DO 2 J2=1,3 DO 2 N =1,NFG C(J1,J2,N,K)= T2(N)*D(J1,N)*D(J2,N) 2 CONTINUE DO 3 J=1,3 DO 3 N=1,NFG C(J,J,N,K)= C(J,J,N,K)+T1(N) 3 CONTINUE 10 CONTINUE DO 20 K=1,NPF KM1= K-1 IF(K.EQ.1) KM1=NPF DO 17 J1=1,3 DO 17 J2=1,3 DO 17 N =1,NFG B(J1,J2,N,K)= C(J1,J2,N,K) +C(J1,J2,N,KM1) 17 CONTINUE DO 18 J=1,3 DO 18 N=1,NFG B(J,J,N,K)= B(J,J,N,K) +1./FL(N,K) 18 CONTINUE 20 CONTINUE RETURN END C********************************************************************** SUBROUTINE FIBFRCZ(XF,STF,REST,NFG,NPF,FRC,UNSTBL) DIMENSION XF(3,NFG,NPF),FRC(3,NFG,NPF) DIMENSION STF(NFG,NPF),REST(NFG,NPF) DIMENSION F(3),D(3) LOGICAL RFLAG LOGICAL UNSTBL(NFG,NPF) C C THIS ROUTINE COMPUTES FIBER FORCES AS IN FIBFRC. C HERE,HOWEVER,THE RESULTS ARE STORED IN THE ARRAY FRC. C C INPUT VARIABLES: C C XF( ,K)= POSITION OF POINT K C STF(K)= STIFFNESS OF LINK K C REST(K)= RESTING LENGTH OF LINK K C NPF= NUMBER OF POINTS IN A FIBER C NFG= NUMBER OF FIBERS IN A GROUP C C OUTPUT VARIABLE: C C FRC( ,K)= FORCE APPLIED BY POINT K TO FLUID C C NOTE THAT LINK K JOINS POINT K AND POINT K+1. C FIBER IS PERIODIC: LINK NPF JOINS POINT NPF AND POINT 1. C C FORCES ARE COMPUTED BY SUMMING OVER LINKS,AS FOLLOWS: C INITIALIZE FRC( , )=0. THEN,FOR K=1,2,...,NPF C LET R= DIST(XF( ,K),XF( ,K+1)) C IF R j do 30 i=j+1,n sum = a(i,j) do 20 k=1,j-1 sum = sum - a(i,k)*a(j,k) 20 continue a(i,j) = sum/a(j,j) 30 continue c 100 continue c c forward solution for y: lower*y=-b (Note: y is stored in x) c do 200 i=1,n sum = -b(i) do 110 j=1,i-1 sum = sum - a(i,j)*x(j) 110 continue x(i) = sum/a(i,i) 200 continue c c backward solution for x: transpose(lower)*x=x c do 300 i=n,1,-1 sum = x(i) do 210 j=i+1,n sum = sum - a(j,i)*x(j) 210 continue x(i) = sum/a(i,i) 300 continue return end C*********************************************************************** subroutine movesrc(xmk,nmarks,xsrc,nsrcs) c c this subroutine computes the locations of each source by averaging c the locations of the markers in the cloud associated with the source. c dimension xmk(3,*) dimension nmarks(nsrcs) dimension xsrc(3,nsrcs) dimension xsum(3) nstop = 0 do 200 isrc=1,nsrcs nstrt = nstop + 1 nstop = nstop + nmarks(isrc) number = nmarks(isrc) do 100 i=1,3 xsrc(i,isrc) = (ssum(number,xmk(i,nstrt),3))/float(number) 100 continue 200 continue c do 300 isrc=1,nsrcs c write(6,250) isrc,(xsrc(i,isrc),i=1,3) c 250 format(' xsrc([1..3],',i1,') = ',3f6.2) c 300 continue return end C**************************************************************************** SUBROUTINE FFT2D(A,B,ISIGN) C*****TWO-DIMENSIONAL FOURIER TRANSFORM OF A COMPLEX ARRAY A+IB. C (THAT IS, THE REAL PART IS STORED IN A AND THE IMMAGINARY C PART IS STORED IN B.) C THE TRANSFORM IS COMPUTED IN PLACE, OVERWRITING A AND B. C VECTORIZATION IS ACHIEVED AS FOLLOWS: C WHEN COMPUTING TRANSFORM IN X-DIRECTION, VECTORIZE OVER Y; C WHEN COMPUTING TRANSFORM IN Y-DIRECTION, VECTORIZE OVER X. C C THE DATA A AND B ARE N X N, C WHERE N IS SPECIFIED IN THE PARAMETER STATEMENT BELOW. C THESE N X N DATA ARE STORED IN ARRAYS C WITH DIMENSIONS 0...NB, 1...NG C C THE ARGUMENT ISIGN MAY BE ANY POSITIVE OR NEGATIVE INTEGER (NOT 0). C ITS SIGN DETERMINES THE SIGN OF THE EXPONENT IN THE FOURIER TRANSFORM. C C NORMALIZATION IS PERFORMED FOR EITHER VALUE OF ISIGN. C THE NORMALIZATION IS DIVISION BY N. (NOTE THAT N=SQRT(N**2).) C C ALTHOUGH THE PROGRAM COMPUTES THE COMPLEX FOURIER TRANSFORM, THE C ARITHMETIC IS WRITTEN OUT IN TERMS OF THE REAL AND IMAGINARY PARTS. C C THE CODE IS A VECTOR VERSION OF THE ONE-DIMENSIONAL FFT C PROGRAM OF COOLEY, LEWIS, AND WELCH, IEEE TRANSACTIONS E-12, C NO. 1 (MARCH 1965), CITED BY DAHLQUIST & BJORK, NUMERICAL METHODS, C PRENTICE HALL, ENGLEWOOD CLIFFS, NEW JERSEY, P416 C PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(M=L2NG,N=2**M) DIMENSION A(0:NB,1:NG),B(0:NB,1:NG) DIMENSION T1(N),T2(N),T3(N),T4(N) SIGN=ISIGN/IABS(ISIGN) NV2=N/2 NM1=N-1 C C*****FOURIER TRANSFORM IN X-DIRECTION: C*****BIT-REVERSAL PERMUTATION: C J=1 DO7I=1,NM1 IF(I.GE.J)GOTO5 DO91 IV=1,N T1(IV)=A(J,IV) T3(IV)=B(J,IV) T2(IV)=A(I,IV) 91 T4(IV)=B(I,IV) DO92IV=1,N A(I,IV)=T1(IV) B(I,IV)=T3(IV) A(J,IV)=T2(IV) 92 B(J,IV)=T4(IV) 5 K=NV2 6 IF(K.GE.J)GO TO 7 J=J-K K=K/2 GO TO 6 7 J=J+K C C*****MAIN LOOP: C DO20L=1,M LE=2**L LE1=LE/2 U1=1. U2=0. ANG=SIGN*3.14159265358979/LE1 W1=COS(ANG) W2=SIN(ANG) DO20J=1,LE1 DO10I=J,N,LE IP=I+LE1 DO93IV=1,N T1(IV)=A(I,IV) T3(IV)=B(I,IV) T2(IV)=A(IP,IV)*U1-B(IP,IV)*U2 93 T4(IV)=A(IP,IV)*U2+B(IP,IV)*U1 DO94IV=1,N A(IP,IV)=T1(IV)-T2(IV) B(IP,IV)=T3(IV)-T4(IV) A(I ,IV)=T1(IV)+T2(IV) 94 B(I ,IV)=T3(IV)+T4(IV) 10 CONTINUE TU1=U1*W1-U2*W2 U2 =U1*W2+U2*W1 20 U1=TU1 C C*****FOURIER TRANSFORM IN Y-DIRECTION: C*****BIT-REVERSAL PERMUTATION: C J=1 DO107I=1,NM1 IF(I.GE.J)GOTO105 DO191 IV=1,N T1(IV)=A(IV,J) T3(IV)=B(IV,J) T2(IV)=A(IV,I) 191 T4(IV)=B(IV,I) DO192IV=1,N A(IV,I)=T1(IV) B(IV,I)=T3(IV) A(IV,J)=T2(IV) 192 B(IV,J)=T4(IV) 105 K=NV2 106 IF(K.GE.J)GO TO 107 J=J-K K=K/2 GO TO 106 107 J=J+K C C*****MAIN LOOP: C DO120L=1,M LE=2**L LE1=LE/2 U1=1. U2=0. ANG=SIGN*3.14159265358979/LE1 W1=COS(ANG) W2=SIN(ANG) DO120J=1,LE1 DO110I=J,N,LE IP=I+LE1 DO193IV=1,N T1(IV)=A(IV,I) T3(IV)=B(IV,I) T2(IV)=A(IV,IP)*U1-B(IV,IP)*U2 193 T4(IV)=A(IV,IP)*U2+B(IV,IP)*U1 DO194IV=1,N A(IV,IP)=T1(IV)-T2(IV) B(IV,IP)=T3(IV)-T4(IV) A(IV,I )=T1(IV)+T2(IV) 194 B(IV,I )=T3(IV)+T4(IV) 110 CONTINUE TU1=U1*W1-U2*W2 U2 =U1*W2+U2*W1 120 U1=TU1 C C NORMALIZATION LOOP: C DO200I=1,N DO200IV=1,N A(I,IV)=A(I,IV)/N B(I,IV)=B(I,IV)/N 200 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE INFLUIDU PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP1=NG+1) PARAMETER(FLNG=NG) C SIZES OF CFFT3D WORK ARRAYS PARAMETER(NTABLE=100+2*3*NG) PARAMETER(NWORK=4*NG*NG+1 ) C COMMON/AR/ U(0:NB,1:NG,0:NGM1) COMMON/AR/ V(0:NB,1:NG,0:NGM1) COMMON/AR/ W(0:NB,1:NG,0:NGM1) COMMON/AR/ P(0:NB,1:NG,0:NGM1) C COMMON/BR/UR(0:NB,0:NB,0:NGM1) COMMON/BR/VR(0:NB,0:NB,0:NGM1) COMMON/BR/WR(0:NB,0:NB,0:NGM1) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR C COMMON/NPM/NP1(NG),NM1(NG) C COMMON/CNST/VSC,ACNST,BCNST C COMMON/BSTAR/TD,TIME,USTAR,PSTAR,FSTAR,RSTAR COMMON/BSTAR/LCUBE,NU,MU,MASS,LENGTH REAL LCUBE,NU,MU,MASS,LENGTH C C NOTE THAT QRFACT IS DIMENSIONED (-1:NGP1,-1:NGP1,0:NGM1) C WHICH IS THE EQUIVALENT OF ( 0:NB , 0:NB ,0:NGM1) C COMMON/FFTFACT/ PRDENO(-1:NGP1,-1:NGP1,0:NGM1) COMMON/FFTFACT/ QRFACT(-1:NGP1,-1:NGP1,0:NGM1) COMMON/FFTFACT/ VRFACT(-1:NGP1) COMMON/FFTFACT/ PRFACT(-1:NGP1) REAL PRDENO,QRFACT COMPLEX VRFACT,PRFACT DIMENSION SINSQ (-1:NGP1) C WORK SPACE REQUIRED BY CFFT3D COMMON/FFTWORK/ FSCALE,INC1X,INC2X,INC3X COMMON/FFTWORK/ TABLE(NTABLE),WORK(NWORK,16) COMMON/KKGWORK/ TRIGX(2*NG),TRIGY(2*NG),TRIGZ(2*NG) COMMON/KKGWORK/ IFAX(19) , IFAY(19) , IFAZ(19) REAL FSCALE INTEGER INC1X,INC2X,INC3X REAL TABLE ,WORK REAL TRIGX ,TRIGY ,TRIGZ INTEGER IFAX , IFAY , IFAZ C MODNG(K) = MOD(K+NG,NG) C C INITIALIZE CONSTANTS FOR TRIDIAGONAL SYSTEMS ACNST = VSC BCNST = 1.+2.*VSC C C INITIALIZE THE ARRAYS THAT ARE USED FOR C PERIODIC SUBSCRIPT ARITHMETIC: DO 5 I=1,NG NP1(I) = I+1 NM1(I) = I-1 5 CONTINUE NP1(NG) = 1 NM1( 1) = NG C C INITIALIZE THE CONSTANT WAVE-NUMBER ARRAYS USED IN THE FLUID SOLVER. PI = 4.*ATAN(1.) TWOPI = 2.*PI FOURNU = 4.*VSC MIDPT = NG/2 DO 10 KS=0,NGM1 VRFACT(KS) = CMPLX(0.,-SIN(TWOPI*KS/FLNG)) PRFACT(KS) = VRFACT(KS) 10 CONTINUE VRFACT(MIDPT) = (0.0,0.0) PRFACT(MIDPT) = VRFACT(MIDPT) DO 20 KS=0,NGM1 SINSQ (KS) = (SIN(PI*KS/FLNG))**2 20 CONTINUE DO 25 K=0,NGM1 DO 25 J=0,NGM1 DO 25 I=0,NGM1 QRFACT(I,J,K) = 1. + FOURNU*(SINSQ(I) + SINSQ(J) + SINSQ(K)) 25 CONTINUE DO 30 KS=0,NGM1 SINSQ(KS) = (SIN(TWOPI*KS/FLNG))**2 30 CONTINUE SINSQ(MIDPT) = 0.0 DO 35 K=0,NGM1 DO 35 J=0,NGM1 DO 35 I=0,NGM1 PRDENO(I,J,K) = (SINSQ(I) + SINSQ(J) + SINSQ(K)) 35 CONTINUE C C PRDENO(-1:NGP1,-1:NGP1,0:NGM1) PRDENO(0:NB,0:NB,0:NGM1) C C I J K C 0, 0, 0 1, 1, 0 C NG/2, 0, 0 NG/2+1, 1, 0 C 0, NG/2, 0 1, NG/2+1, 0 C NG/2, NG/2, 0 NG/2+1, NG/2+1, 0 C 0, 0, NG/2 1, 1, NG/2 C NG/2, 0, NG/2 NG/2+1, 1, NG/2 C 0, NG/2, NG/2 1, NG/2+1, NG/2 C NG/2, NG/2, NG/2 NG/2+1, NG/2+1, NG/2 C PRDENO( 0, 0, 0) = 1.0 PRDENO(MIDPT, 0, 0) = 1.0 PRDENO( 0,MIDPT, 0) = 1.0 PRDENO(MIDPT,MIDPT, 0) = 1.0 PRDENO( 0, 0,MIDPT) = 1.0 PRDENO(MIDPT, 0,MIDPT) = 1.0 PRDENO( 0,MIDPT,MIDPT) = 1.0 PRDENO(MIDPT,MIDPT,MIDPT) = 1.0 FSCALE = (1./FLNG**3) INC1X = 1 INC2X = NB+1 INC3X = (NB+1)**2 C INITIALIZE THE ARRAYS TABLE AND WORK CALL CFFT3D( 0,NG,NG,NG,FSCALE, C UR(1,1,0),INC1X,INC2X,INC3X, C UR(1,1,0),INC1X,INC2X,INC3X, C TABLE,NTABLE,WORK,NWORK) RETURN END SUBROUTINE MYCFFT3D(ISIGN, NSIZ1, NSIZ2, NSIZ3, FSCALE, C A, INC1A, INC2A, INC3A, C B, INC1B, INC2B, INC3B, C U, NTABLE, WORK, NWORK ) PARAMETER(L2NG=7,NG=2**L2NG,NB=NG+2,NGM1=NG-1) PARAMETER(NBP1=NB+1) PARAMETER(M=L2NG,N=2**M) COMPLEX A(1:NBP1,1:NBP1,1:NG) COMPLEX B(1:NBP1,1:NBP1,1:NG) COMPLEX U(2*N,-1:1) COMPLEX T(1:N),W IF (ISIGN .EQ. 0) THEN PI = 3.141592653589793238462643 CMIC$ DO ALL CMIC$1 SHARED(PI, U) CMIC$2 PRIVATE(ANG, IS, J, L, LE, LE1, SIGN, W) DO 3 IS=-1,1,2 SIGN = IS DO 2 L=1,M LE = 2**L LE1 = LE/2 U(LE1,IS) = (1.,0.) ANG = SIGN*PI/LE1 W = CMPLX(COS(ANG),SIN(ANG)) IF (LE1 .EQ. 1) THEN W = (-1., 0.) ELSEIF (LE1 .EQ. 2) THEN W = CMPLX(0., SIGN) END IF DO 2 J=1,LE1-1 U(LE1+J,IS) = U(LE1+J-1,IS)*W 2 CONTINUE 3 CONTINUE RETURN END IF SIGN = ISIGN/IABS(ISIGN) IS = SIGN NV2 = N/2 NM1 = N-1 CMIC$ DO ALL CMIC$1SHARED(A, B, FSCALE) CMIC$2PRIVATE(IW, IV, I) DO 1 IW=1,N DO 1 IV=1,N DO 1 I=1,N B(I,IV,IW) = FSCALE*A(I,IV,IW) 1 CONTINUE CDIR$ IVDEP CMIC$ DO ALL CMIC$1SHARED(A, B, IS, NM1, NV2, U) CMIC$2PRIVATE(IW, IV, IP, I, J, K, L, LE, LE1, T) DO 30 IW=1,N J = 1 DO 7 I=1,NM1 IF (I.LT.J) THEN DO 91 IV=1,N T(IV) = B(J,IV,IW) 91 CONTINUE DO 92 IV=1,N B(J,IV,IW) = B(I,IV,IW) B(I,IV,IW) = T( IV) 92 CONTINUE END IF 5 K = NV2 6 IF (K.LT.J) THEN J = J-K K = K/2 GO TO 6 END IF J = J+K 7 CONTINUE DO 20 L=1,M LE = 2**L LE1 = LE/2 C U = (1.,0.) C ANG = SIGN*3.14159265358979/LE1 C W = CMPLX(COS(ANG),SIN(ANG)) DO 20 J=1,LE1 DO 10 I=J,N,LE IP = I+LE1 DO 95 IV=1,N T(IV)= B(IP,IV,IW)*U(LE1+J-1,IS) 95 CONTINUE DO 96 IV=1,N B(IP,IV,IW) = B(I,IV,IW)-T(IV) B( I,IV,IW) = B(I,IV,IW)+T(IV) 96 CONTINUE 10 CONTINUE C U = U*W 20 CONTINUE 30 CONTINUE CDIR$ IVDEP CMIC$ DO ALL CMIC$1SHARED(A, B, IS, NM1, NV2, U) CMIC$2PRIVATE(IW, IV, IP, I, J, K, L, LE, LE1, T) DO 130 IW=1,N J = 1 DO 107 I=1,NM1 IF (I.LT.J) THEN DO 191 IV=1,N T(IV) = B(IW,J,IV) 191 CONTINUE DO 192 IV=1,N B(IW,J,IV) = B(IW,I,IV) B(IW,I,IV) = T( IV) 192 CONTINUE END IF 105 K = NV2 106 IF (K.LT.J) THEN J = J-K K = K/2 GO TO 106 END IF J = J+K 107 CONTINUE DO 120 L=1,M LE = 2**L LE1 = LE/2 C U = (1.,0.) C ANG = SIGN*3.14159265358979/LE1 C W = CMPLX(COS(ANG),SIN(ANG)) DO 120 J=1,LE1 DO 110 I=J,N,LE IP = I+LE1 DO 195 IV=1,N T(IV)= B(IW,IP,IV)*U(LE1+J-1,IS) 195 CONTINUE DO 196 IV=1,N B(IW,IP,IV) = B(IW,I ,IV)-T(IV) B(IW,I ,IV) = B(IW,I ,IV)+T(IV) 196 CONTINUE 110 CONTINUE C U = U*W 120 CONTINUE 130 CONTINUE CDIR$ IVDEP CMIC$ DO ALL CMIC$1SHARED(A, B, IS, NM1, NV2, U) CMIC$2PRIVATE(IW, IV, IP, I, J, K, L, LE, LE1, T) DO 230 IW=1,N J = 1 DO 207 I=1,NM1 IF (I.LT.J) THEN DO 291 IV=1,N T(IV) = B(IV,IW,J) 291 CONTINUE DO 292 IV=1,N B(IV,IW,J) = B(IV,IW,I) B(IV,IW,I) = T(IV) 292 CONTINUE END IF 205 K = NV2 206 IF (K.LT.J) THEN J = J-K K = K/2 GO TO 206 END IF J = J+K 207 CONTINUE DO 220 L=1,M LE = 2**L LE1 = LE/2 C U = (1.,0.) C ANG = SIGN*3.14159265358979/LE1 C W = CMPLX(COS(ANG),SIN(ANG)) DO 220 J=1,LE1 DO 210 I=J,N,LE IP = I+LE1 DO 295 IV=1,N T(IV)= B(IV,IW,IP)*U(LE1+J-1,IS) 295 CONTINUE DO 296 IV=1,N B(IV,IW,IP) = B(IV,IW,I )-T(IV) B(IV,IW,I ) = B(IV,IW,I )+T(IV) 296 CONTINUE 210 CONTINUE C U = U*W 220 CONTINUE 230 CONTINUE RETURN END subroutine cfft3d (isign,nx,ny,nz,fscale, c z ,idum1x,idum1y,idum1z, c zdum,idum2x,idum2y,idum2z, c tabdum,ntabdum,work,nwork) c subroutine f3dinc (nx,ny,nz,z,work,trigx,trigy,trigz, c c ifax,ifay,ifaz,isign,is,ns,nwork,fscale) c c...... Description: c c F3DINC is a procedure to compute an in-core 3-D complex FFT. c c c...... Arguments: c c nx,ny,nz : size of the FFT (integers) c c z : Input FFT declared as complex z(nx+is,ny+is,nz) c c work : Real array of size work(nwork,ns) c c trigx: Real array of size (2 * Nx), initialized cftfax(nx,ifax,trigx) c c trigy: Real array of size (2 * Ny), initialized cftfax(ny,ifay,trigy) c c trigz: Real array of size (2 * Nz), initialized cftfax(nz,ifaz,trigz) c c ifax: integer array of size 19, initialized cftfax(nx,ifax,trigx) c c ifay: integer array of size 19, initialized cftfax(ny,ifay,trigy) c c ifaz: integer array of size 19, initialized cftfax(nz,ifaz,trigz) c c isign: integer +1 for forward transform or -1 for inverse c determines the sign of the exponent in c exp (isign*2*pi*j*k/N) * f (j,k) c c is: a pad to make leading dimensions odd c c ns: number of parallel chunks, typically number of available CPUs c c nwork: workspace required by each CPU, at least 4*n*n, n=max(nx,ny,nz) c c..... Method c c Do Nz 2-D FFTs of the X-Y planes c i.e. c (1) For all planes 1..Nz c Perform NY simultaneous 1-D FFTs of length NX c c and c c (2) For all planes 1..Nz c Perform NX simultaneous 1-D FFTs of length NY c c Finally, c (3) For all planes 1..NY c Perform NX simultaneous 1-D FFTs of NZ c c------------------------------ implicit none c------------------------------ c integer idum1x,idum1y,idum1z,idum2x,idum2y,idum2z real zdum real tabdum integer ntabdum integer is ! padding to give odd first dimensions integer ns ! # planes in a slab integer nx,ny,nz,n,nxb,nyb,m c parameter(ns=16) parameter(ns=8) parameter(is=3) INTEGER L2NG, NG PARAMETER(L2NG=7,NG=2**L2NG) COMMON/KKGWORK/ TRIGX(2*NG),TRIGY(2*NG),TRIGZ(2*NG) COMMON/KKGWORK/ IFAX(19) , IFAY(19) , IFAZ(19) REAL TRIGX ,TRIGY ,TRIGZ INTEGER IFAX , IFAY , IFAZ external second,cftfax,timef,cfftmlt real timef c integer nx,ny,nz,n,nxb,nyb,m c c integer is c integer ns ! # planes in a slab integer ioff integer isign c integer second real ss integer ix,iy,iz,ip integer incx,incy,incz,jumpx,jumpy,jumpz,lot c integer ifax (19), ifay (19), ifaz (19) integer nwork real work (nwork,ns) c real trigx (2*nx), trigy (2*ny), trigz (2*nz) real z (2,nx+is,ny+is,nz) real fscale integer isc,jsc,ksc c c---------------------------------------------------------------- c Executable code c---------------------------------------------------------------- c if (isign .eq. 0) then call cftfax (nx,ifax,trigx) call cftfax (ny,ifay,trigy) call cftfax (nz,ifaz,trigz) return end if c jumpx = 2*(nx + is) incx = 2 c incy = 2*(nx + is) jumpy = 2 c c ------------------------------------- c For all planes 1..Nz do a 2-D FFT c Perform NY simultaneous NX point FFT c ------------------------------------- c do iz = 1,nz,ns cmic$ doall autoscope cmic$. shared (z,nx,ny,nz,jumpx,jumpy,jumpz,incx,incy,incz,work, cmic$. isign,ns,trigx,trigy,trigz,ifax,ifay,ifaz) do ip = 1,min (ns, nz-iz+1) call cfftmlt (z(1,1,1,iz+ip-1),z(2,1,1,iz+ip-1), . work(1,ip),trigx,ifax,incx,jumpx,nx,ny,isign) enddo enddo c c ------------------------------------- c For all planes 1..Nz do a 2-D FFT c Perform nx simultaneous ny point FFT c ------------------------------------- c do iz = 1,nz,ns cmic$ doall autoscope cmic$. shared (z,nx,ny,nz,jumpx,jumpy,jumpz,incx,incy,incz,work, cmic$. isign,ns,trigx,trigy,trigz,ifax,ifay,ifaz) do ip = 1,min (ns, nz-iz+1) call cfftmlt (z(1,1,1,iz+ip-1),z(2,1,1,iz+ip-1), . work(1,ip),trigy,ifay,incy,jumpy,ny,nx,isign) enddo enddo c c---------------------------------------------------------------- c---------------------------------------------------------------- c c c c -------------------------------- c NX 1-D transforms of length NZ c Repeat for all planes 1..NY c -------------------------------- c incz = 2*(nx+is)*(ny+is) jumpz = 2 do iy = 1,ny,ns cmic$ doall autoscope cmic$. shared (z,nx,ny,nz,jumpx,jumpy,jumpz,incx,incy,incz,work, cmic$. isign,ns,trigx,trigy,trigz,ifax,ifay,ifaz) do ip = 1,min (ns, ny-iy+1) call cfftmlt (z(1,1,iy+ip-1,1),z(2,1,iy+ip-1,1), . work(1,ip),trigz,ifaz,incz,jumpz,nz,nx,isign) enddo enddo c---------------------------------------------------------------- c---------------------------------------------------------------- c C SCALE THE TRANSFORMED DATA CMIC$ DO ALL CMIC$1SHARED(FSCALE, NX, NY, NZ, Z) CMIC$2PRIVATE(ISC, JSC, KSC) DO 520 KSC=1,NZ DO 520 JSC=1,NY DO 520 ISC=1,NX Z(1,ISC,JSC,KSC) = FSCALE*Z(1,ISC,JSC,KSC) Z(2,ISC,JSC,KSC) = FSCALE*Z(2,ISC,JSC,KSC) 520 CONTINUE return end