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 hrtxp_859.f second-order Immersed Boundary Method version 1 c hrtxp_863.f second-order Immersed Boundary Method version 2 c hrtxq_1113.f changed Cray parallel directives to OpenMP directives c hrtxq_1113a.f fixed bug related to restart 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=8,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=3*(2*NG+256)) 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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/HR/FU(0:NB,0:NB,0:NGM1) COMMON/HR/FV(0:NB,0:NB,0:NGM1) COMMON/HR/FW(0:NB,0:NB,0:NGM1) COMPLEX FU,FV,FW 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,2) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ tsrc common/source/ qsrc(nsrcs, 2) common/source/ psrc(nsrcspx,2) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:4095, nsrcspx) common/source/indxps(0:4095, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1,2) 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,2) 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,2) 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 REAL xmk_old_val(3),xmk_old_val2(3) !old values of xmk, used to store prev time steps 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(14),ttimer(14),rtimer(14) REAL*4 SECOND 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 = 1113 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 = 512 IF (DENOVA) THEN c NSTEP = 256*locref NSTEP = 2048 + 1 !2048+1 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,1),nsize*ng) call zero(v (0,1,0,1),nsize*ng) call zero(w (0,1,0,1),nsize*ng) call zero(u (0,1,0,2),nsize*ng) call zero(v (0,1,0,2),nsize*ng) call zero(w (0,1,0,2),nsize*ng) ELSE CALL RESTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,1,NREFOLD,NREF) NSTEP = 2048 + 1 !2048+1 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) CALL COPYX3D(XFN,NFSIZE) 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 COPYX3D(XMK,NMSIZE) 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),1),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 do 12 ist=1,2 do 12 isr=1,nsrcs psrc (isr,ist) = 0.0 qsrc (isr,ist) = 0.0 12 continue 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 source/sink inertial time constant tsrc = 400.0 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,2) xatapex(2) = xfn(2,katapex,2) xatapex(3) = xfn(3,katapex,2) c correct initial data is now stored in the checkpoint file. c provide correct initial data for COMPLEX QR(,,,2) and c provide correct initial data for COMPLEX UR(,,,2) following a restart. c if (KLOK0 .gt. 0) then c ISTEP0 = 0 c call sourceup(ISTEP0) c call initur c end if 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) c STOP END IF C C WRITE OUT SOME OUTPUT FOR TIME T-- C klokout = mod(klok,4) T = 0. + FLOAT(KLOK0)*FLOAT(NREF)*TD IF (KLOK0 .EQ. 0) THEN CALL WRVPXF(EXPNUM,KLOK,T,XFN,NFG,NPF,NGROUPS,KSTART,NFSTART, C NEST,NFIBER,XMK,NCLOUD,NMARKS,NFSKIP,NPSKIP, C FRC,STF,REST) END IF C C MAIN LOOP -- !TLAST = SECOND() call cpu_time(tlast) DO 5 KLOK=KLOK1,KLOKEND INTWR = 512 !512 INTOUT = 512 KLOKWR = MOD(KLOK, INTWR) KLOKOUT = MOD(KLOK, INTOUT) DO 4 NR=1,NREF T = T+TD c c set ISTEP as if this iteration of DO 5 follows from the end of the c previous iteration, even if this is the first iteration. ISTEP = 2 c c sort the linked lists using xfn(,,2) rtxf1 = rtc() !txfl1 = second() call cpu_time(txfl1) CALL XFLIST(XFN,NEXTN,ISTEP) !txfl2 = second() call cpu_time(txfl2) rtxf2 = rtc() timer(1) = txfl2 - txfl1 rtimer(1) = (rtxf2 - rtxf1)*8.0E-07 c c begin the first step of the second-order method ISTEP = 1 c c compute xfn(,,1) rtmov1 = rtc() !tmove1 = second() call cpu_time(tmove1) CALL MOVE(XFN,NEXTN,NFG,NPF,NGROUPS,NBUNCH, c katapex,xatapex,ISTEP) !tmove2 = second() call cpu_time(tmove2) rtmov2 = rtc() timer(2) = tmove2 - tmove1 rtimer(2) = (rtmov2 - rtmov1)*8.0E-07 c c sort the linked lists using xmk(,,2) c compute xmk(,,1) rtmvm1 = rtc() !tmovm1 = second() call cpu_time(tmovm1) CALL MOVEMK(XMK,NEXTM,NMARKT,ISTEP) !tmovm2 = second() call cpu_time(tmovm2) rtmvm2 = rtc() timer(3) = tmovm2 - tmovm1 rtimer(3) = (rtmvm2 - rtmvm1)*8.0E-07 c c sort the linked lists using xfn(,,1) rtxf1 = rtc() !txfl1 = second() call cpu_time(txfl1) CALL XFLIST(XFN,NEXTN,ISTEP) !txfl2 = second() call cpu_time(txfl2) rtxf2 = rtc() timer(4) = txfl2 - txfl1 rtimer(4) = (rtxf2 - rtxf1)*8.0E-07 c c compute the locations of the sources. call locsrc(nmsrc,nmtrl,nrvmp,npulm,xatapex,ISTEP) C C UPDATE THE MUSCLE ACTIVATION FUNCTIONS CALL ACTIVATE(KLOK,1) C C EVALUATE FORCE DENSITY rtfi1 = rtc() !tfib1 = second() call cpu_time(tfib1) C C EXPLICIT FORCE CALCULATION c since ISTEP is 1, use the first half of XFN in FIBERX and in FBEAMS 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() call cpu_time(tfib2) rtfi2 = rtc() timer(5) = tfib2 - tfib1 rtimer(5) = (rtfi2 - rtfi1)*8.0E-07 C C APPLY FORCE DENSITY TO FLUID rtpus1 = rtc() !tpush1 = second() call cpu_time(tpush1) c c compute the lattice forces from the fiber forces c since ISTEP is 1, use the first half of XFN in FORCES CALL FORCES(XFN,FRC,NEXTN,NFG,NPF,NGROUPS,NBUNCH) c c initialize (UR(,,,1),VR(,,,1),WR(,,,1)) with (U(,,,2),V(,,,2),W(,,,2)) CALL SKEWSM(ISTEP) c c update (UR(,,,1),VR(,,,1),WR(,,,1)) with (delta_t/rho)*(FU,FV,FW) CALL ADDFRC(ISTEP) !tpush2 = second() call cpu_time(tpush2) rtpus2 = rtc() timer(6) = tpush2 - tpush1 rtimer(6) = (rtpus2 - rtpus1)*8.0E-07 C C ADVANCE FLUID VELOCITY ONE TIME STEP -- rtfl1 = rtc() !tflu1 = second() call cpu_time(tflu1) CALL FLUID(nref,ISTEP) !tflu2 = second() call cpu_time(tflu2) rtfl2 = rtc() timer(7) = tflu2 - tflu1 rtimer(7) = (rtfl2 - rtfl1)*8.0E-07 c c begin the second step of the second-order method ISTEP = 2 c c compute xfn(,,2) [the linked lists are already sorted for this] rtmov1 = rtc() !tmove1 = second() call cpu_time(tmove1) CALL MOVE(XFN,NEXTN,NFG,NPF,NGROUPS,NBUNCH, c katapex,xatapex,ISTEP) !tmove2 = second() call cpu_time(tmove2) rtmov2 = rtc() timer(8) = tmove2 - tmove1 rtimer(8) = (rtmov2 - rtmov1)*8.0E-07 c c sort the linked lists using xmk(,,1) c compute xmk(,,2) rtmvm1 = rtc() !tmovm1 = second() call cpu_time(tmovm1) CALL MOVEMK(XMK,NEXTM,NMARKT,ISTEP) !tmovm2 = second() call cpu_time(tmovm2) rtmvm2 = rtc() timer(9) = tmovm2 - tmovm1 rtimer(9) = (rtmvm2 - rtmvm1)*8.0E-07 c c compute the locations of the sources. call locsrc(nmsrc,nmtrl,nrvmp,npulm,xatapex,ISTEP) C C APPLY FORCE DENSITY TO FLUID rtpus1 = rtc() !tpush1 = second() call cpu_time(tpush1) c c initialize (UR(,,,2),VR(,,,2),WR(,,,2)) with (U(,,,2),V(,,,2),W(,,,2)) CALL VISCON c c update (UR(,,,2),VR(,,,2),WR(,,,2)) with (UR(,,,1),VR(,,,1),WR(,,,1)) CALL SKEWSM(ISTEP) c c update (UR(,,,2),VR(,,,2),WR(,,,2)) with (2*delta_t/rho)*(FU,FV,FW) c use the same values of (FU,FV,FW) as computed by SUBROUTINE FORCES above. CALL ADDFRC(ISTEP) !tpush2 = second() call cpu_time(tpush2) rtpus2 = rtc() timer(10) = tpush2 - tpush1 rtimer(10) = (rtpus2 - rtpus1)*8.0E-07 c c solve for (UR(,,,2),VR(,,,2),WR(,,,2)) rtfl1 = rtc() !tflu1 = second() call cpu_time(tflu1) CALL FLUID(NREF,ISTEP) !tflu2 = second() call cpu_time(tflu2) rtfl2 = rtc() timer(11) = tflu2 - tflu1 rtimer(11) = (rtfl2 - rtfl1)*8.0E-07 c c check the maximum values of the new velocity field rtchk1 = rtc() !tchks1 = second() call cpu_time(tchks1) call chkspeed(speedm,ISTEP) !tchks2 = second() call cpu_time(tchks2) rtchk2 = rtc() timer(12) = tchks2 - tchks1 rtimer(12) = (rtchk2 - rtchk1)*8.0E-07 4 CONTINUE C C COMPUTE THE FLOWS THROUGH THE VALVE RINGS if (klokout .eq. 0) then !tflo1 = second() call cpu_time(tflo1) CALL FLOWS(xmk(1,1,ISTEP),nextw(1),nextw(nfwby2),nmarks,ncircs, c expnum,klok,klokout,nref,ISTEP) !tflo2 = second() call cpu_time(tflo2) timer(13) = tflo2 - tflo1 c c compute the flow out of each of the sources do 50 isrc=1,5 call srcflow(isrc,xsrc(1,isrc,2),nextw(1),nextw(nfwby2), c klok,klokout,nref,ISTEP) 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) ELSEIF (KLOKWR .EQ. 1) THEN C WRITE PSRC FOR A POSSIBLE FUTURE RESTART WITH TIMESTEP REFINEMENT MULTWR = KLOK/INTWR IUNIT = 2 - MOD(MULTWR,2) IUNITP = IUNIT + 30 CALL WRITPSRC(KLOK,IUNITP) END IF !twr1 = second() call cpu_time(twr1) IF (KLOKOUT .EQ. 0) THEN C WRITE OUT SOME OUTPUT CALL WRVPXF(EXPNUM,KLOK,T,XFN(1,1,ISTEP), C NFG,NPF,NGROUPS,KSTART,NFSTART, C NEST,NFIBER,XMK(1,1,ISTEP), C NCLOUD,NMARKS,NFSKIP,NPSKIP, C FRC,STF,REST) END IF !twr2 = second() call cpu_time(twr2) timer(14) = twr2 - twr1 !TTHIS = SECOND() call cpu_time(tthis) 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,'(a18,1x,a14)') ' second',' rtc' write(6,'(a10,f8.4,1x,f14.9)')'xflist ',timer( 1),rtimer( 1) write(6,'(a10,f8.4,1x,f14.9)')'move ',timer( 2),rtimer( 2) write(6,'(a10,f8.4,1x,f14.9)')'movemk ',timer( 3),rtimer( 3) write(6,'(a10,f8.4,1x,f14.9)')'xflist ',timer( 4),rtimer( 4) write(6,'(a10,f8.4,1x,f14.9)')'fiberx ',timer( 5),rtimer( 5) write(6,'(a10,f8.4,1x,f14.9)')'push up ',timer( 6),rtimer( 6) write(6,'(a10,f8.4,1x,f14.9)')'fluid ',timer( 7),rtimer( 7) write(6,'(a10,f8.4,1x,f14.9)')'move ',timer( 8),rtimer( 8) write(6,'(a10,f8.4,1x,f14.9)')'movemk ',timer( 9),rtimer( 9) write(6,'(a10,f8.4,1x,f14.9)')'push skew ',timer(10),rtimer(10) write(6,'(a10,f8.4,1x,f14.9)')'fluid ',timer(11),rtimer(11) write(6,'(a10,f8.4,1x,f14.9)')'chkspeed ',timer(12),rtimer(12) write(6,'(a10,f8.4)') 'flows ',timer(13) write(6,'(a10,f8.4)') 'wrvpxf ',timer(14) 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 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=8,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*17 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) 301 format(a3,i4,a6,i1,a3) 302 format(a3,i4,a5,i2,a3) 303 format(a3,i4,a4,i3,a3) 304 format(a3,i4,a3,i4,a3) 305 format(a3,i4,a2,i5,a3) 306 format(a3,i4,a1,i6,a3) 307 format(a3,i4, i7,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 elseif (expnum .lt. 1000) then 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 else if ( klok .lt. 10 ) then write(filemk,301) 'hrt' ,expnum,'_00000',klok,'.mk' write(filevp,301) 'hrt' ,expnum,'_00000',klok,'.vp' write(filexf,301) 'hrt' ,expnum,'_00000',klok,'.xf' elseif (klok .lt. 100 ) then write(filemk,302) 'hrt' ,expnum,'_0000' ,klok,'.mk' write(filevp,302) 'hrt' ,expnum,'_0000' ,klok,'.vp' write(filexf,302) 'hrt' ,expnum,'_0000' ,klok,'.xf' elseif (klok .lt. 1000 ) then write(filemk,303) 'hrt' ,expnum,'_000' ,klok,'.mk' write(filevp,303) 'hrt' ,expnum,'_000' ,klok,'.vp' write(filexf,303) 'hrt' ,expnum,'_000' ,klok,'.xf' elseif (klok .lt. 10000 ) then write(filemk,304) 'hrt' ,expnum,'_00' ,klok,'.mk' write(filevp,304) 'hrt' ,expnum,'_00' ,klok,'.vp' write(filexf,304) 'hrt' ,expnum,'_00' ,klok,'.xf' elseif (klok .lt. 100000 ) then write(filemk,305) 'hrt' ,expnum,'_0' ,klok,'.mk' write(filevp,305) 'hrt' ,expnum,'_0' ,klok,'.vp' write(filexf,305) 'hrt' ,expnum,'_0' ,klok,'.xf' elseif (klok .lt. 1000000) then write(filemk,306) 'hrt' ,expnum,'_' ,klok,'.mk' write(filevp,306) 'hrt' ,expnum,'_' ,klok,'.vp' write(filexf,306) 'hrt' ,expnum,'_' ,klok,'.xf' else write(filemk,307) 'hrt' ,expnum, klok,'.mk' write(filevp,307) 'hrt' ,expnum, klok,'.vp' write(filexf,307) '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,a17 ,/, 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=8,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 SUBROUTINE WRSTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,IUNIT) PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE,2) 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,2) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/CNST/VSC,ACNST,BCNST COMMON/NPM/NP1(NG),NM1(NG) common/source/xsrc(3,nsrcspx,2) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ tsrc common/source/ qsrc(nsrcs, 2) common/source/ psrc(nsrcspx,2) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:4095, nsrcspx) common/source/indxps(0:4095, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1,2) 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,tsrc,qsrc,psrc write(iunit) qr 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 WRITE(IUNIT) P WRITE(IUNIT) UR WRITE(IUNIT) VR WRITE(IUNIT) WR WRITE(IUNIT) PR 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 SUBROUTINE RESTART(NSTEP,LCUBE,NU,RHO,TD,ESTOP,ITMAX, C KLOK,index,NREFOLD,NREF) save iunit PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE,2) 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,2) COMMON/MAR/ NEXTM(NMSIZE) COMMON/MAR/ NMARKS(NCLMAX) COMMON/CNST/VSC,ACNST,BCNST COMMON/NPM/NP1(NG),NM1(NG) common/source/xsrc(3,nsrcspx,2) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ tsrc common/source/ qsrc(nsrcs, 2) common/source/ psrc(nsrcspx,2) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:4095, nsrcspx) common/source/indxps(0:4095, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1,2) 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,tsrc,qsrc,psrc read(iunit) qr 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 READ(IUNIT) P READ(IUNIT) UR READ(IUNIT) VR READ(IUNIT) WR READ(IUNIT) PR 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) call interpsrc(NREF,NREFOLD,psrc) IF (NREF .NE. NREFOLD) CALL REFINE(NREFOLD,NREF) end if RETURN END SUBROUTINE REFINE(NREFOLD,NREF) PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) COMMON/XAR/ XF (3,NFSIZE), XFN (3,NFSIZE,2) 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,2) common/source/resist(nsrcs) common/source/prsrvr(nsrcs) common/source/ tsrc common/source/ qsrc(nsrcs, 2) common/source/ psrc(nsrcspx,2) c common/source/ aps(nsrcs,0:nsrcs) c common/source/ pe(0:nb,1:ng,0:1,nsrcspx) common/source/ ps(0:4095, nsrcspx) common/source/indxps(0:4095, nsrcspx) common/source/ qr(0:nb,0:nb,0:ngm1,2) 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 ip=1,2 do 9 nsrc=1,nsrcspx psrc(nsrc,ip) = psrc(nsrc,ip)*FACTSQ 9 continue do 10 nsrc=1,nsrcs resist(nsrc) = resist(nsrc)*FACT prsrvr(nsrc) = prsrvr(nsrc)*FACTSQ 10 continue do 11 iq=1,2 do 11 nsrc=1,nsrcs qsrc(nsrc,iq) = qsrc(nsrc,iq)*FACT 11 continue tsrc = tsrc/FACT DO 20 L=1,2 DO 20 K=0,NGM1 DO 20 J=1,NG DO 20 I=1,NG U (I,J,K,L) = U (I,J,K,L)*FACT V (I,J,K,L) = V (I,J,K,L)*FACT W (I,J,K,L) = W (I,J,K,L)*FACT P (I,J,K,L) = P (I,J,K,L)*FACTSQ 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=8,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=8,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=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) character*17 filevp WRITE(IUNIT,9)KLOK,filevp 9 FORMAT( I5,10X,' = KLOK',1x,a17, 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,2),V(I,J,K,2),W(I,J,K,2), C P(I,J,K,2), I,J,K 100 CONTINUE RETURN END 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 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 6 FORMAT(1x,a28,4f10.6) RETURN END 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 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=8,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 Recompute REST0 C DO 1225 I=ISTART,ISTOP K = KSTART(I) CALL REREST(XF(1,K),REST0(K),NFG(I),NPF(I)) 1225 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 COPYX3D(XDATA,NSIZE) C C COPY VALUES FROM XDATA(,,1) TO XDATA(,,2) C DIMENSION XDATA(3,NSIZE,2) DO 10 N=1,NSIZE XDATA(1,N,2) = XDATA(1,N,1) XDATA(2,N,2) = XDATA(2,N,1) XDATA(3,N,2) = XDATA(3,N,1) 10 CONTINUE 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=8,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 c do 308 np=1,npf c if (laflag(nf,np) .eq. 1) then c write(6,"(3f7.2,2i4,a)")xf(1,nf,np),xf(2,nf,np),xf(3,nf,np), c c nf,np," RV FREE WALL" c else c write(6,"(3f7.2,2i4,a)")xf(1,nf,np),xf(2,nf,np),xf(3,nf,np), c c nf,np," SEPTUM" c end if c 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,rcchr)-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 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) logical first first = .true. xmax = xlvout + alvout xmin = xright - aright 1 x = (xmax+xmin)/2. if (.not.first) then if (x .eq. xold) then go to 2 end if end if first = .false. xold = x write(6,*) "in intrsct, x=",x 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 REREST(XF,REST,NFG,NPF) DIMENSION XF(3,NFG,NPF) DIMENSION REST(NFG,NPF) DO 200 NP=1,NPF 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 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 INITUR c c Provide correct initial data for COMPLEX UR(,,,2) following a restart c PARAMETER(L2NG=8,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP1=NG+1) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/AR/ U(0:NB,1:NG,0:NGM1,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) DO 100 K=0,NGM1 DO 100 J=1,NG DO 100 I=1,NG UR(I,J,K,2) = U(I,J,K,2) VR(I,J,K,2) = V(I,J,K,2) WR(I,J,K,2) = W(I,J,K,2) 100 CONTINUE RETURN END C**************************************************************************** SUBROUTINE SKEWSM(ISTEP) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP1=NG+1) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/AR/ U(0:NB,1:NG,0:NGM1,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) IF (ISTEP .EQ. 1) THEN C$OMP PARALLEL DO C$OMP1SHARED(U, V, W, UR, VR, WR) C$OMP2PRIVATE(K, KM1, KP1, J, JM1, JP1, I) C$OMP+SCHEDULE(RUNTIME) DO 101 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 101 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,2) = U (NG,J,K,2) U (NGP1,J,K,2) = U ( 1,J,K,2) V ( 0,J,K,2) = V (NG,J,K,2) V (NGP1,J,K,2) = V ( 1,J,K,2) W ( 0,J,K,2) = W (NG,J,K,2) W (NGP1,J,K,2) = W ( 1,J,K,2) DO 101 I=1,NG UR(I,J,K,1) = U (I ,J ,K ,2)*2.0 C -( ( U (I ,J ,K ,2)*( U (I+1,J ,K ,2) C -U (I-1,J ,K ,2)) C +V (I ,J ,K ,2)*( U (I ,JP1,K ,2) C -U (I ,JM1,K ,2)) C +W (I ,J ,K ,2)*( U (I ,J ,KP1,2) C -U (I ,J ,KM1,2)))/2.0 C +( U (I+1,J ,K ,2)* U (I+1,J ,K ,2) C -U (I-1,J ,K ,2)* U (I-1,J ,K ,2) C +V (I ,JP1,K ,2)* U (I ,JP1,K ,2) C -V (I ,JM1,K ,2)* U (I ,JM1,K ,2) C +W (I ,J ,KP1,2)* U (I ,J ,KP1,2) C -W (I ,J ,KM1,2)* U (I ,J ,KM1,2) )/2.0 C )/2.0 VR(I,J,K,1) = V (I ,J ,K ,2)*2.0 C -( ( U (I ,J ,K ,2)*( V (I+1,J ,K ,2) C -V (I-1,J ,K ,2)) C +V (I ,J ,K ,2)*( V (I ,JP1,K ,2) C -V (I ,JM1,K ,2)) C +W (I ,J ,K ,2)*( V (I ,J ,KP1,2) C -V (I ,J ,KM1,2)))/2.0 C +( U (I+1,J ,K ,2)* V (I+1,J ,K ,2) C -U (I-1,J ,K ,2)* V (I-1,J ,K ,2) C +V (I ,JP1,K ,2)* V (I ,JP1,K ,2) C -V (I ,JM1,K ,2)* V (I ,JM1,K ,2) C +W (I ,J ,KP1,2)* V (I ,J ,KP1,2) C -W (I ,J ,KM1,2)* V (I ,J ,KM1,2) )/2.0 C )/2.0 WR(I,J,K,1) = W (I ,J ,K ,2)*2.0 C -( ( U (I ,J ,K ,2)*( W (I+1,J ,K ,2) C -W (I-1,J ,K ,2)) C +V (I ,J ,K ,2)*( W (I ,JP1,K ,2) C -W (I ,JM1,K ,2)) C +W (I ,J ,K ,2)*( W (I ,J ,KP1,2) C -W (I ,J ,KM1,2)))/2.0 C +( U (I+1,J ,K ,2)* W (I+1,J ,K ,2) C -U (I-1,J ,K ,2)* W (I-1,J ,K ,2) C +V (I ,JP1,K ,2)* W (I ,JP1,K ,2) C -V (I ,JM1,K ,2)* W (I ,JM1,K ,2) C +W (I ,J ,KP1,2)* W (I ,J ,KP1,2) C -W (I ,J ,KM1,2)* W (I ,J ,KM1,2) )/2.0 C )/2.0 101 CONTINUE ELSE C$OMP PARALLEL DO C$OMP1SHARED(UR, VR, WR) C$OMP2PRIVATE(K, KM1, KP1, J, JM1, JP1, I) C$OMP+SCHEDULE(RUNTIME) 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 UR( 0,J,K,1) = UR(NG,J,K,1) UR(NGP1,J,K,1) = UR( 1,J,K,1) VR( 0,J,K,1) = VR(NG,J,K,1) VR(NGP1,J,K,1) = VR( 1,J,K,1) WR( 0,J,K,1) = WR(NG,J,K,1) WR(NGP1,J,K,1) = WR( 1,J,K,1) DO 105 I=1,NG UR(I,J,K,2) = UR(I ,J ,K ,2) C -( ( UR(I ,J ,K ,1)*( UR(I+1,J ,K ,1) C -UR(I-1,J ,K ,1)) C +VR(I ,J ,K ,1)*( UR(I ,JP1,K ,1) C -UR(I ,JM1,K ,1)) C +WR(I ,J ,K ,1)*( UR(I ,J ,KP1,1) C -UR(I ,J ,KM1,1)))/2.0 C +( UR(I+1,J ,K ,1)* UR(I+1,J ,K ,1) C -UR(I-1,J ,K ,1)* UR(I-1,J ,K ,1) C +VR(I ,JP1,K ,1)* UR(I ,JP1,K ,1) C -VR(I ,JM1,K ,1)* UR(I ,JM1,K ,1) C +WR(I ,J ,KP1,1)* UR(I ,J ,KP1,1) C -WR(I ,J ,KM1,1)* UR(I ,J ,KM1,1) )/2.0 C ) VR(I,J,K,2) = VR(I ,J ,K ,2) C -( ( UR(I ,J ,K ,1)*( VR(I+1,J ,K ,1) C -VR(I-1,J ,K ,1)) C +VR(I ,J ,K ,1)*( VR(I ,JP1,K ,1) C -VR(I ,JM1,K ,1)) C +WR(I ,J ,K ,1)*( VR(I ,J ,KP1,1) C -VR(I ,J ,KM1,1)))/2.0 C +( UR(I+1,J ,K ,1)* VR(I+1,J ,K ,1) C -UR(I-1,J ,K ,1)* VR(I-1,J ,K ,1) C +VR(I ,JP1,K ,1)* VR(I ,JP1,K ,1) C -VR(I ,JM1,K ,1)* VR(I ,JM1,K ,1) C +WR(I ,J ,KP1,1)* VR(I ,J ,KP1,1) C -WR(I ,J ,KM1,1)* VR(I ,J ,KM1,1) )/2.0 C ) WR(I,J,K,2) = WR(I ,J ,K ,2) C -( ( UR(I ,J ,K ,1)*( WR(I+1,J ,K ,1) C -WR(I-1,J ,K ,1)) C +VR(I ,J ,K ,1)*( WR(I ,JP1,K ,1) C -WR(I ,JM1,K ,1)) C +WR(I ,J ,K ,1)*( WR(I ,J ,KP1,1) C -WR(I ,J ,KM1,1)))/2.0 C +( UR(I+1,J ,K ,1)* WR(I+1,J ,K ,1) C -UR(I-1,J ,K ,1)* WR(I-1,J ,K ,1) C +VR(I ,JP1,K ,1)* WR(I ,JP1,K ,1) C -VR(I ,JM1,K ,1)* WR(I ,JM1,K ,1) C +WR(I ,J ,KP1,1)* WR(I ,J ,KP1,1) C -WR(I ,J ,KM1,1)* WR(I ,J ,KM1,1) )/2.0 C ) 105 CONTINUE END IF RETURN END C**************************************************************************** SUBROUTINE UPWIND c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP1=NG+1) COMMON/AR/ U(0:NB,1:NG,0:NGM1,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR C$OMP PARALLEL DO C$OMP1SHARED(U, V, W, UR, VR, WR) C$OMP2PRIVATE(K, KM1, KP1, J, JM1, JP1, I) C$OMP+SCHEDULE(RUNTIME) 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,2) = U(NG,J,K,2) U(NGP1,J,K,2) = U( 1,J,K,2) V( 0,J,K,2) = V(NG,J,K,2) V(NGP1,J,K,2) = V( 1,J,K,2) W( 0,J,K,2) = W(NG,J,K,2) W(NGP1,J,K,2) = W( 1,J,K,2) DO 105 I=1,NG UR(I,J,K,1) = U(I,J,K,2)*2.0 C -U(I,J,K,2)*CVMGP(U(I ,J ,K ,2)-U(I-1,J ,K ,2), C U(I+1,J ,K ,2)-U(I ,J ,K ,2), C U(I ,J ,K ,2)) C -V(I,J,K,2)*CVMGP(U(I ,J ,K ,2)-U(I ,JM1,K ,2), C U(I ,JP1,K ,2)-U(I ,J ,K ,2), C V(I ,J ,K ,2)) C -W(I,J,K,2)*CVMGP(U(I ,J ,K ,2)-U(I ,J ,KM1,2), C U(I ,J ,KP1,2)-U(I ,J ,K ,2), C W(I ,J ,K ,2)) VR(I,J,K,1) = V(I,J,K,2)*2.0 C -U(I,J,K,2)*CVMGP(V(I ,J ,K ,2)-V(I-1,J ,K ,2), C V(I+1,J ,K ,2)-V(I ,J ,K ,2), C U(I ,J ,K ,2)) C -V(I,J,K,2)*CVMGP(V(I ,J ,K ,2)-V(I ,JM1,K ,2), C V(I ,JP1,K ,2)-V(I ,J ,K ,2), C V(I ,J ,K ,2)) C -W(I,J,K,2)*CVMGP(V(I ,J ,K ,2)-V(I ,J ,KM1,2), C V(I ,J ,KP1,2)-V(I ,J ,K ,2), C W(I ,J ,K ,2)) WR(I,J,K,1) = W(I,J,K,2)*2.0 C -U(I,J,K,2)*CVMGP(W(I ,J ,K ,2)-W(I-1,J ,K ,2), C W(I+1,J ,K ,2)-W(I ,J ,K ,2), C U(I ,J ,K ,2)) C -V(I,J,K,2)*CVMGP(W(I ,J ,K ,2)-W(I ,JM1,K ,2), C W(I ,JP1,K ,2)-W(I ,J ,K ,2), C V(I ,J ,K ,2)) C -W(I,J,K,2)*CVMGP(W(I ,J ,K ,2)-W(I ,J ,KM1,2), C W(I ,J ,KP1,2)-W(I ,J ,K ,2), C W(I ,J ,K ,2)) 105 CONTINUE RETURN END C**************************************************************************** SUBROUTINE VISCON c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) PARAMETER(NGP1=NG+1) COMMON/AR/ U(0:NB,1:NG,0:NGM1,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/AR/ P(0:NB,1:NG,0:NGM1,2) COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR COMMON/CNST/VSC,ACNST,BCNST C THIS ROUTINE MAKES A VISCOUS-TERM CONTRIBUTION TO U, V, W C AND ALSO MAKES A PRESSURE-TERM CONTRIBUTION TO U, V, W C FOR ALL I, J, K: C UR(I,J,K,2)= 2.0*U(I,J,K) C + VSC*( U(I+1,J,K)+U(I-1,J,K) -2.0*U(I,J,K) C +U(I,J+1,K)+U(I,J-1,K) -2.0*U(I,J,K) C +U(I,J,K+1)+U(I,J,K-1) -2.0*U(I,J,K) ) C C = (2.0-6.0*VSC)*U(I,J,K) C + VSC*( U(I+1,J,K)+U(I-1,J,K) C +U(I,J+1,K)+U(I,J-1,K) C +U(I,J,K+1)+U(I,J,K-1) ) C C AND SIMILARLY FOR VR(I,J,K,2) AND WR(I,J,K,2). C VCNST = 2.0 - 6.0*VSC C$OMP PARALLEL DO C$OMP1SHARED(U, V, W, P, UR, VR, WR, VCNST, VSC) C$OMP2PRIVATE(K, KM1, KP1, J, JM1, JP1, I) C$OMP+SCHEDULE(RUNTIME) 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,2) = U(NG,J,K,2) U(NGP1,J,K,2) = U( 1,J,K,2) V( 0,J,K,2) = V(NG,J,K,2) V(NGP1,J,K,2) = V( 1,J,K,2) W( 0,J,K,2) = W(NG,J,K,2) W(NGP1,J,K,2) = W( 1,J,K,2) DO 105 I=1,NG UR(I,J,K,2) = VCNST*U(I ,J ,K ,2) C +VSC*( U(I+1,J ,K ,2) + U(I-1,J ,K ,2) C +U(I ,JP1,K ,2) + U(I ,JM1,K ,2) C +U(I ,J ,KP1,2) + U(I ,J ,KM1,2)) VR(I,J,K,2) = VCNST*V(I ,J ,K ,2) C +VSC*( V(I+1,J ,K ,2) + V(I-1,J ,K ,2) C +V(I ,JP1,K ,2) + V(I ,JM1,K ,2) C +V(I ,J ,KP1,2) + V(I ,J ,KM1,2)) WR(I,J,K,2) = VCNST*W(I ,J ,K ,2) C +VSC*( W(I+1,J ,K ,2) + W(I-1,J ,K ,2) C +W(I ,JP1,K ,2) + W(I ,JM1,K ,2) C +W(I ,J ,KP1,2) + W(I ,J ,KM1,2)) 105 CONTINUE RETURN END C**************************************************************************** SUBROUTINE FORCES(XFN,FRC,NEXTN,NFG,NPF,NGROUPS,NBUNCH) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,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/HR/FU(0:NB,0:NB,0:NGM1) COMMON/HR/FV(0:NB,0:NB,0:NGM1) COMMON/HR/FW(0:NB,0:NB,0:NGM1) COMPLEX FU,FV,FW 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 FULIN(-15:16*NGP2),FVLIN(-15:16*NGP2),FWLIN(-15:16*NGP2) DIMENSION INDUVW( 16*NG ) DIMENSION FLKZP1(NPTMAX) dimension lindx(NPTMAX) dimension timer(10) REAL*4 SECOND C C THIS ROUTINE SPREADS THE FIBER FORCES ONTO FLUID-LATTICE FORCE ARRAYS. 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 ZERO OUT THE FLUID-LATTICE FORCE ARRAYS FU, FV, FW. 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 FU(I,J,K) = FU(I,J,K) + D * FRC(1,L) C FV(I,J,K) = FV(I,J,K) + D * FRC(2,L) C FW(I,J,K) = FW(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.) C$OMP PARALLEL DO C$OMP1SHARED(FU, FV, FW) C$OMP2PRIVATE(K, J, I) C$OMP+SCHEDULE(RUNTIME) DO 105 K=0,NGM1 DO 105 J=0,NB DO 105 I=0,NB FU(I,J,K) = (0.0,0.0) FV(I,J,K) = (0.0,0.0) FW(I,J,K) = (0.0,0.0) 105 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 C$OMP PARALLEL DO C$OMP1SHARED(FIRSTN, FRC, MSHIFT, NEXTN, NUMBER, XFN, C$OMP2 IZERO, JZERO, FU, FV, FW) C$OMP3PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, C$OMP4 FLIZP1, FLJZP1, FLKZP1, FORCE1, FORCE2, FORCE3, C$OMP5 IJ, I, I0, I3D, II, INDUVW, IZ, J, J3D, JJ, JZ, K, KZ, L, C$OMP6 LINDX, M, MZERO, NPOINTS, NPT, NUMREM, RAD1, RAD2, RAD3, C$OMP7 XFN1OLD, XFN2OLD, XFN3OLD, C$OMP8 FULIN, FVLIN, FWLIN) C$OMP+SCHEDULE(RUNTIME) 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 FULIN(M) = 0.0 FVLIN(M) = 0.0 FWLIN(M) = 0.0 41 CONTINUE CALL GATHER(MSHIFT,FULIN(1),FU(0,0,0),INDUVW(1)) CALL GATHER(MSHIFT,FVLIN(1),FV(0,0,0),INDUVW(1)) CALL GATHER(MSHIFT,FWLIN(1),FW(0,0,0),INDUVW(1)) DO 42 M=16*NG+1,16*NGP2 FULIN(M) = 0.0 FVLIN(M) = 0.0 FWLIN(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 if ((1.+4.*ARG2*(1.-ARG2))<0.) then print *, "arg2 1sqrt error", NPT,XFN2OLD(NPT),FLJZP1 end if 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 FULIN(M+MZERO) = FULIN(M+MZERO) + DELTA(M,NPT)*FORCE1(NPT) FVLIN(M+MZERO) = FVLIN(M+MZERO) + DELTA(M,NPT)*FORCE2(NPT) FWLIN(M+MZERO) = FWLIN(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 FULIN(M) = FULIN(M) + FULIN(M+MSHIFT) FVLIN(M) = FVLIN(M) + FVLIN(M+MSHIFT) FWLIN(M) = FWLIN(M) + FWLIN(M+MSHIFT) 14 CONTINUE DO 15 M=16*(NG-1)+1,16*NG FULIN(M) = FULIN(M) + FULIN(M-MSHIFT) FVLIN(M) = FVLIN(M) + FVLIN(M-MSHIFT) FWLIN(M) = FWLIN(M) + FWLIN(M-MSHIFT) 15 CONTINUE ctim t2 = second() ctim timer(9) = timer(9) + t2 - t1 ctim t1 = t2 CALL SCATTER(MSHIFT,FU(0,0,0),INDUVW(1),FULIN(1)) CALL SCATTER(MSHIFT,FV(0,0,0),INDUVW(1),FVLIN(1)) CALL SCATTER(MSHIFT,FW(0,0,0),INDUVW(1),FWLIN(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 SUBROUTINE ADDFRC(ISTEP) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,NG=2**L2NG,NB=NG+2) PARAMETER(NSIZE=(NB+1)*NG) PARAMETER(NGM1=NG-1) COMMON/HR/FU(0:NB,0:NB,0:NGM1) COMMON/HR/FV(0:NB,0:NB,0:NGM1) COMMON/HR/FW(0:NB,0:NB,0:NGM1) COMPLEX FU,FV,FW COMMON/BR/UR(0:NB,0:NB,0:NGM1,2) COMMON/BR/VR(0:NB,0:NB,0:NGM1,2) COMMON/BR/WR(0:NB,0:NB,0:NGM1,2) COMMON/BR/PR(0:NB,0:NB,0:NGM1) COMPLEX UR,VR,WR,PR STEPSIZ = FLOAT(ISTEP) C$OMP PARALLEL DO C$OMP1SHARED(ISTEP, STEPSIZ, FU, FV, FW, UR, VR, WR) C$OMP2PRIVATE(K, J, I) C$OMP+SCHEDULE(RUNTIME) DO 100 K=0,NGM1 DO 100 J=1,NG DO 100 I=1,NG UR(I,J,K,ISTEP) = UR(I,J,K,ISTEP) + STEPSIZ*FU(I,J,K) VR(I,J,K,ISTEP) = VR(I,J,K,ISTEP) + STEPSIZ*FV(I,J,K) WR(I,J,K,ISTEP) = WR(I,J,K,ISTEP) + STEPSIZ*FW(I,J,K) 100 CONTINUE RETURN END C**************************************************************************** SUBROUTINE XFLIST(XFN,NEXTN,INDEX) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,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/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) INTEGER FIRSTN DIMENSION XFN(3,NFSIZE,2),NEXTN(*) dimension timer(10) REAL*4 SECOND MODNG(K)=MOD(K+NG,NG) 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 C$OMP PARALLEL DO C$OMP1SHARED(FIRSTN, NEXTN, NUMBER, IZERO, JZERO, XFN, INDEX) C$OMP2PRIVATE(II, IJ, IN, JJ, JN, N, NEXTNOLD, NPREV) C$OMP+SCHEDULE(RUNTIME) 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,INDEX)+FLNG)-1) + 1 IN = MODNG(INT(XFN(1,N,INDEX)+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 RETURN END C********************************************************************** SUBROUTINE MOVE(XFN,NEXTN,NFG,NPF,NGROUPS,NBUNCH, c katapex,xatapex,ISTEP) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) COMMON/HIST/ FIRSTN(1:NG,1:NG) COMMON/HIST/ NUMBER(1:NG,1:NG) INTEGER FIRSTN DIMENSION XFN(3,NFSIZE,2),NEXTN(*) DIMENSION NFG(*),NPF(*),NGROUPS(*) DIMENSION XFN1OLD(NPTMAX),XFN2OLD(NPTMAX),XFN3OLD(NPTMAX) DIMENSION XFN1DEL(NPTMAX),XFN2DEL(NPTMAX),XFN3DEL(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) REAL*4 SECOND 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 THE INTEGERS IXFNNEW,IXFNDEL ARE USED TO SELECT THE VALUE C OF THE LAST ARRAY INDEX OF XFN, U, V, AND W, C AND STEPSIZ IS USED TO SET THE SIZE OF THE TIME STEP THUSLY: C C XFN(,,IXFNNEW) = XFN(,,0) + STEPSIZ*D(XFN(,,IXFNDEL))*U(,,,IXFNDEL) C C ISTEP 1 2 C IXFNNEW 1 2 C IXFNDEL 2 1 C STEPSIZ 0.5 1.0 C IXFNNEW = ISTEP IXFNDEL = 3 - ISTEP STEPSIZ = FLOAT(ISTEP)/2.0 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 C$OMP PARALLEL DO C$OMP1SHARED(FIRSTN, NEXTN, NUMBER, U, V, W, XFN, C$OMP1 ISTEP, IXFNNEW, IXFNDEL, STEPSIZ) C$OMP2PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, C$OMP3 FLIZP1, FLJZP1, FLKZP1, C$OMP4 IJ, I, I0, I3D, II, IZ, J, J3D, JJ, JZ, K, K3D, KZ, L, C$OMP5 LINDX, M, MZERO, NPOINTS, NUMREM, NPT, RAD1, RAD2, RAD3, C$OMP6 UINT, ULIN, VINT, C$OMP7 VLIN, WINT, WLIN, XFN1OLD, XFN2OLD, XFN3OLD, C$OMP8 XFN1DEL, XFN2DEL, XFN3DEL) C$OMP+SCHEDULE(RUNTIME) 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 THE VALUE OF THE LAST INDEX OF THE 3D VELOCITY ARRAY DEPENDS ON THE STEP. 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,IXFNDEL) VLIN(I0+I) = V(I3D,J3D,K3D,IXFNDEL) WLIN(I0+I) = W(I3D,J3D,K3D,IXFNDEL) 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 C C ALWAYS GATHER XFNOLD FROM XFN(,,2) C IN STEP 1, GATHER XFNDEL FROM XFN(,,2) (ACTUALLY, COPY FROM XFNOLD) C IN STEP 2, GATHER XFNDEL FROM XFN(,,1) C CALL GATHER(NPOINTS,XFN1OLD,XFN(1,1,2),LINDX) CALL GATHER(NPOINTS,XFN2OLD,XFN(2,1,2),LINDX) CALL GATHER(NPOINTS,XFN3OLD,XFN(3,1,2),LINDX) IF (ISTEP .EQ. 1) THEN DO 602 NPT=1,NPOINTS XFN1DEL(NPT) = XFN1OLD(NPT) XFN2DEL(NPT) = XFN2OLD(NPT) XFN3DEL(NPT) = XFN3OLD(NPT) 602 CONTINUE ELSE CALL GATHER(NPOINTS,XFN1DEL,XFN(1,1,1),LINDX) CALL GATHER(NPOINTS,XFN2DEL,XFN(2,1,1),LINDX) CALL GATHER(NPOINTS,XFN3DEL,XFN(3,1,1),LINDX) ENDIF IZ = INT(XFN1DEL( 1) - 1. + FLNG) - NG JZ = INT(XFN2DEL( 1) - 1. + FLNG) - NG FLIZP1 = IZ+1 FLJZP1 = JZ+1 DO 6 NPT=1,NPOINTS KZ = INT(XFN3DEL(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 = XFN3DEL(NPT) - FLKZP1(NPT) ARG2 = XFN2DEL(NPT) - FLJZP1 ARG1 = XFN1DEL(NPT) - FLIZP1 if ((1.+4.*ARG2*(1.-ARG2))<0.) then print *, "arg2 2sqrt error", NPT,XFN2DEL(NPT),FLJZP1 end if 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(XFN3DEL(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 C C ALWAYS UPDATE XFNOLD. C THE VELOCITY DATA WITH WHICH XFNOLD IS UPDATED DEPENDS ON THE STEP. C ctim t1 = t2 DO 124 M=1,64 DO 123 NPT=1,NPOINTS XFN1OLD(NPT) = XFN1OLD(NPT) + STEPSIZ*UINT(M,NPT) XFN2OLD(NPT) = XFN2OLD(NPT) + STEPSIZ*VINT(M,NPT) XFN3OLD(NPT) = XFN3OLD(NPT) + STEPSIZ*WINT(M,NPT) 123 CONTINUE 124 CONTINUE ctim t2 = second() ctim timer(9) = timer(9) + t2 - t1 C C THE VALUE OF THE LAST INDEX OF XFN DEPENDS ON THE STEP, C BUT THE COORDINATE VALUES BEING WRITTEN ARE ALWAYS IN XFNOLD. C ctim t1 = t2 CALL SCATTER(NPOINTS,XFN(1,1,IXFNNEW),LINDX,XFN1OLD) CALL SCATTER(NPOINTS,XFN(2,1,IXFNNEW),LINDX,XFN2OLD) CALL SCATTER(NPOINTS,XFN(3,1,IXFNNEW),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,ixfnnew) xatapex(2) = xfn(2,katapex,ixfnnew) xatapex(3) = xfn(3,katapex,ixfnnew) RETURN END SUBROUTINE CHKSPEED(SPEEDM,ISTEP) c c This subroutine contains microtasking directives. c PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) 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 C$OMP PARALLEL DO C$OMP1SHARED(U, V, W, UMAX, VMAX, WMAX, SMAX, ISTEP) C$OMP2PRIVATE( ULIN , VLIN , WLIN , SLIN , C$OMP3 UINDEX, VINDEX, WINDEX, SINDEX, C$OMP4 UI0 , VI0 , WI0 , SI0 , C$OMP5 I , J , K , S ) C$OMP+SCHEDULE(RUNTIME) DO 100 K=0,NGM1 C C FOR EACH PLANE K, CREATE THE ARRAY S WHOSE MAXIMUM IS TO BE FOUND C DO 5 J=1,NG DO 5 I=1,NG S(I,J) = ABS(U(I,J,K,ISTEP)) 2 + ABS(V(I,J,K,ISTEP)) 3 + ABS(W(I,J,K,ISTEP)) 5 CONTINUE C C FIND THE ROW INDEX OF THE MAXIMUM VALUE IN EACH COLUMN C DO 10 J=1,NG UINDEX(J) = IDAMAX(NG,U(1,J,K,ISTEP),1) VINDEX(J) = IDAMAX(NG,V(1,J,K,ISTEP),1) WINDEX(J) = IDAMAX(NG,W(1,J,K,ISTEP),1) SINDEX(J) = IDAMAX(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,ISTEP),UINDEX(1)) CALL GATHER(NG,VLIN(1),V(0,1,K,ISTEP),VINDEX(1)) CALL GATHER(NG,WLIN(1),W(0,1,K,ISTEP),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 = IDAMAX(NG,ULIN(1),1) VI0 = IDAMAX(NG,VLIN(1),1) WI0 = IDAMAX(NG,WLIN(1),1) SI0 = IDAMAX(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 = IDAMAX(NG,UMAX(0),1) - 1 VIM = IDAMAX(NG,VMAX(0),1) - 1 WIM = IDAMAX(NG,WMAX(0),1) - 1 SIM = IDAMAX(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,ISTEP) c c This subroutine contains microtasking directives. c USE, INTRINSIC :: IEEE_ARITHMETIC PARAMETER(L2NG=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) 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,NMSIZE,2),NEXTN(*) DIMENSION XMK1OLD(NPTMAX),XMK2OLD(NPTMAX),XMK3OLD(NPTMAX) DIMENSION XMK1DEL(NPTMAX),XMK2DEL(NPTMAX),XMK3DEL(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.) C THE INTEGERS IXMKNEW,IXMKDEL ARE USED TO SELECT THE VALUE C OF THE LAST ARRAY INDEX OF XMK, U, V, AND W, C AND STEPSIZ IS USED TO SET THE SIZE OF THE TIME STEP THUSLY: C C XMK(,,IXMKNEW) = XMK(,,0) + STEPSIZ*D(XMK(,,IXMKDEL))*U(,,,IXMKDEL) C C ISTEP 1 2 C IXMKNEW 1 2 C IXMKDEL 2 1 C STEPSIZ 0.5 1.0 C IXMKNEW = ISTEP IXMKDEL = 3 - ISTEP STEPSIZ = FLOAT(ISTEP)/2.0 C C UPDATE THE LINKED LISTS WHICH SORT THE XMK DATA BY (X,Y) COORDINATE C SORT BY THE DATA IN XMK(,,IXMKDEL) C DO 89 IZERO=1,4 DO 89 JZERO=1,4 C$OMP PARALLEL DO C$OMP1SHARED(FIRSTN, NEXTN, NUMBER, IZERO, JZERO, XMK, IXMKDEL) C$OMP2PRIVATE(II, IJ, IN, JJ, JN, N, NEXTNOLD, NPREV) C$OMP+SCHEDULE(RUNTIME) 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,IXMKDEL)+FLNG)-1) + 1 IN = MODNG(INT(XMK(1,N,IXMKDEL)+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 C$OMP PARALLEL DO C$OMP1SHARED(FIRSTN, NEXTN, NUMBER, U, V, W, XMK, C$OMP1 ISTEP, IXMKNEW, IXMKDEL, STEPSIZ) C$OMP2PRIVATE(ARG1, ARG2, ARG3, D1, D12, D2, D3, DELTA, C$OMP3 FLIZP1, FLJZP1, FLKZP1, C$OMP4 IJ, I, I0, I3D, II, IZ, J, J3D, JJ, JZ, K, K3D, KZ, L, C$OMP5 LINDX, M, MZERO, NPOINTS, NUMREM, NPT, RAD1, RAD2, RAD3, C$OMP6 UINT, ULIN, VINT, C$OMP7 VLIN, WINT, WLIN, XMK1OLD, XMK2OLD, XMK3OLD, C$OMP8 XMK1DEL, XMK2DEL, XMK3DEL) C$OMP+SCHEDULE(RUNTIME) 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 THE VALUE OF THE LAST INDEX OF THE 3D VELOCITY ARRAY DEPENDS ON THE STEP. 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,IXMKDEL) VLIN(I0+I) = V(I3D,J3D,K3D,IXMKDEL) WLIN(I0+I) = W(I3D,J3D,K3D,IXMKDEL) 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 C C ALWAYS GATHER XMKOLD FROM XMK(,,2) C IN STEP 1, GATHER XMKDEL FROM XMK(,,2) (ACTUALLY, COPY FROM XMKOLD) C IN STEP 2, GATHER XMKDEL FROM XMK(,,1) C CALL GATHER(NPOINTS,XMK1OLD,XMK(1,1,2),LINDX) CALL GATHER(NPOINTS,XMK2OLD,XMK(2,1,2),LINDX) CALL GATHER(NPOINTS,XMK3OLD,XMK(3,1,2),LINDX) do ii=1,3 if (ieee_is_nan(xmk(ii,1,2))==.true.) then print *, "gather error",ii end if end do IF (ISTEP .EQ. 1) THEN DO 602 NPT=1,NPOINTS XMK1DEL(NPT) = XMK1OLD(NPT) XMK2DEL(NPT) = XMK2OLD(NPT) XMK3DEL(NPT) = XMK3OLD(NPT) if (ieee_is_nan(XMK1DEL(NPT))==.true.) then print *, "XMK1DEL error" end if if (ieee_is_nan(XMK2DEL(NPT))==.true.) then print *, "XMK2DEL error" end if if (ieee_is_nan(XMK3DEL(NPT))==.true.) then print *, "XMK3DEL error" end if 602 CONTINUE ELSE CALL GATHER(NPOINTS,XMK1DEL,XMK(1,1,1),LINDX) CALL GATHER(NPOINTS,XMK2DEL,XMK(2,1,1),LINDX) CALL GATHER(NPOINTS,XMK3DEL,XMK(3,1,1),LINDX) do ii=1,3 if (ieee_is_nan(xmk(ii,1,1))==.true.) then print *, "gather xmk error" end if end do ENDIF IZ = INT(XMK1DEL( 1) - 1. + FLNG) - NG JZ = INT(XMK2DEL( 1) - 1. + FLNG) - NG FLIZP1 = IZ+1 FLJZP1 = JZ+1 DO 6 NPT=1,NPOINTS KZ = INT(XMK3DEL(NPT) - 1. + FLNG) - NG FLKZP1(NPT) = KZ+1 if (ieee_is_nan(FLIZP1)==.true.) then print *, "FLIZP1 error" print *, FLIZP1,iz,XMK1DEL,flng,ng end if if (ieee_is_nan(FLJZP1)==.true.) then print *, "FLJZP1 error" print *, FLJZP1,jz,XMK2DEL,flng,ng end if if (ieee_is_nan(FLKZP1(NPT))==.true.) then print *, "FLKZP1(NPT) error" print *, FLKZP1(NPT),kz,flng,ng end if 6 CONTINUE DO 7 NPT=1,NPOINTS ARG3 = XMK3DEL(NPT) - FLKZP1(NPT) ARG2 = XMK2DEL(NPT) - FLJZP1 ARG1 = XMK1DEL(NPT) - FLIZP1 !print *, "XMK1DEL(2)",XMK1DEL(2),FLIZP1 !print *, "XMK2DEL(2)",XMK2DEL(2),FLJZP1 !print *, "XMK3DEL(2)",XMK3DEL(2),FLKZP1(NPT) if ((1.+4.*ARG2*(1.-ARG2))<0.) then print *, "arg2 3sqrt error", NPT,XMK2DEL(NPT) print *, FLJZP1,JZ,XMK2DEL( 1),XMK1DEL(NPT),FLIZP1 print *, IZ,XMK1DEL( 1),XMK3DEL(NPT) print *, FLKZP1(NPT),KZ,XMK3DEL(NPT),FLNG !ARG2 = XMK2DEL(NPT) - int(XMK2DEL(NPT)) end if RAD3 = SQRT(1.+4.*ARG3*(1.-ARG3)) RAD2 = SQRT(1.+4.*ARG2*(1.-ARG2)) RAD1 = SQRT(1.+4.*ARG1*(1.-ARG1)) if ((1.+4.*ARG2*(1.-ARG2))<0.) then print *, "arg2 sqrt error", NPT,XMK2DEL(NPT),FLJZP1 end if 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. do ii=0,3 if (ieee_is_nan(D1(NPT,ii))==.true.) then print *, "d1 error" print *, ii,npt,D1(NPT,ii),rad1 end if if (ieee_is_nan(D2(NPT,ii))==.true.) then print *, "d2 error" print *, ii,npt,D2(NPT,ii),rad2 end if if (ieee_is_nan(D3(NPT,ii))==.true.) then print *, "d3 error" print *, ii,npt,D3(NPT,ii),rad3 end if end do if (ieee_is_nan(arg1)==.true.) then print *, "arg1 error" print *, arg1,XMK1DEL(NPT),FLIZP1 end if if (ieee_is_nan(arg2)==.true.) then print *, "arg2 error" print *, arg2,XMK2DEL(NPT),FLJZP1 end if if (ieee_is_nan(arg3)==.true.) then print *, "arg3 error" print *, arg1,XMK3DEL(NPT),FLKZP1(NPT) end if 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) if (ieee_is_nan(D12(NPT,I,J))==.true.) then print *, "d12 error b4" print *, i,j,npt,D12(NPT,I,J) end if 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 if (ieee_is_nan(DELTA(M,NPT))==.true.) then print *, "delta error b4" print *, M,NPT,DELTA(M,NPT) end if DELTA(M,NPT) = D12(NPT,I,J) * D3(NPT,K) if (ieee_is_nan(DELTA(M,NPT))==.true.) then print *, "delta error after" print *, M,NPT,delta(m,npt),D12(NPT,I,J),D3(NPT,K) end if 10 CONTINUE DO 122 NPT=1,NPOINTS MZERO = 16*(INT(XMK3DEL(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) if (ieee_is_nan(UINT(M,NPT))==.true.) then print *, "UINT error" print *, M,NPT end if if (ieee_is_nan(VINT(M,NPT))==.true.) then print *, "VINT error" print *, M,NPT end if if (ieee_is_nan(WINT(M,NPT))==.true.) then print *, "WINT error" print *, M,NPT end if 121 CONTINUE 122 CONTINUE C C ALWAYS UPDATE XMKOLD. C THE VELOCITY DATA WITH WHICH XMKOLD IS UPDATED DEPENDS ON THE STEP. C DO 124 M=1,64 DO 123 NPT=1,NPOINTS XMK1OLD(NPT) = XMK1OLD(NPT) + STEPSIZ*UINT(M,NPT) XMK2OLD(NPT) = XMK2OLD(NPT) + STEPSIZ*VINT(M,NPT) XMK3OLD(NPT) = XMK3OLD(NPT) + STEPSIZ*WINT(M,NPT) if (ieee_is_nan(XMK1OLD(NPT))==.true.) then print *, "XMK1OLD error" print *, NPT,M end if if (ieee_is_nan(XMK2OLD(NPT))==.true.) then print *, "XMK2OLD error" print *, NPT,M end if if (ieee_is_nan(XMK3OLD(NPT))==.true.) then print *, "XMK3OLD error" print *, NPT,M end if 123 CONTINUE 124 CONTINUE C C THE VALUE OF THE LAST INDEX OF XMK DEPENDS ON THE STEP, C BUT THE COORDINATE VALUES BEING WRITTEN ARE ALWAYS IN XMKOLD. C CALL SCATTER(NPOINTS,XMK(1,1,IXMKNEW),LINDX,XMK1OLD) CALL SCATTER(NPOINTS,XMK(2,1,IXMKNEW),LINDX,XMK2OLD) CALL SCATTER(NPOINTS,XMK(3,1,IXMKNEW),LINDX,XMK3OLD) do ii=1,NPOINTS if (ieee_is_nan(XMK(1,1,IXMKNEW))==.true.) then print *, "XMK_1 error",ii end if if (ieee_is_nan(XMK(2,1,IXMKNEW))==.true.) then print *, "XMK_2 error",ii end if if (ieee_is_nan(XMK(3,1,IXMKNEW))==.true.) then print *, "XMK_3 error",ii end if end do IF (NUMREM .NE. 0) GO TO 5 20 CONTINUE RETURN END 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=8,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 SUBROUTINE INHIST PARAMETER(L2NG=8,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,2) 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,2) 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) ISTEP = 2 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,ISTEP)+FLNG)-1) + 1 II=MODNG(INT(XFN(1,N,ISTEP)+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,ISTEP)+FLNG)-1) + 1 II=MODNG(INT(XMK(1,N,ISTEP)+FLNG)-1) + 1 NEXTM(N)=FIRSTM(II,JJ) FIRSTM(II,JJ)=N NUMBEM(II,JJ)=NUMBEM(II,JJ)+1 4 CONTINUE RETURN END subroutine flows (xmk,nextnm,nextnc,nmarks,ncircs, c expnum,klok,klokout,nref,istep) 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=8,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*17 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) 301 format(a3,i4,a6,i1,a3) 302 format(a3,i4,a5,i2,a3) 303 format(a3,i4,a4,i3,a3) 304 format(a3,i4,a3,i4,a3) 305 format(a3,i4,a2,i5,a3) 306 format(a3,i4,a1,i6,a3) 307 format(a3,i4, i7,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 elseif (expnum .lt. 1000) then 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 else if ( klok .lt. 10 ) then write(filefl,301) 'hrt' ,expnum,'_00000',klok,'.fl' elseif (klok .lt. 100 ) then write(filefl,302) 'hrt' ,expnum,'_0000' ,klok,'.fl' elseif (klok .lt. 1000 ) then write(filefl,303) 'hrt' ,expnum,'_000' ,klok,'.fl' elseif (klok .lt. 10000 ) then write(filefl,304) 'hrt' ,expnum,'_00' ,klok,'.fl' elseif (klok .lt. 100000 ) then write(filefl,305) 'hrt' ,expnum,'_0' ,klok,'.fl' elseif (klok .lt. 1000000) then write(filefl,306) 'hrt' ,expnum,'_' ,klok,'.fl' else write(filefl,307) '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,istep) 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,a17, 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 ,3f9.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 subroutine srcflow(isrc,xsrc,nextnm,nextnc,klok,klokout,nref, c istep) 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=8,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,istep) 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 subroutine emflow(xmk,nextnm,nextnc,nmarks,ncircs,xcen, c rad,xnorm,axnorm,xweb,uweb,aweb,acmp, c flow,iflag,istep) 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=8,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,2) COMMON/AR/ V(0:NB,1:NG,0:NGM1,2) COMMON/AR/ W(0:NB,1:NG,0:NGM1,2) 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) if ((1.+4.*ARG2*(1.-ARG2))<0.) then print *, "arg2 4sqrt error", NPT,XM2OLD(NPT),FLJZP1 end if 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,ISTEP),INDEX(0,NPT)) CALL GATHER(64,VLIN(0),V(0,1,0,ISTEP),INDEX(0,NPT)) CALL GATHER(64,WLIN(0),W(0,1,0,ISTEP),INDEX(0,NPT)) UWEB(1,NM,NC) = DDOT(64,DELTA(0,NPT),1,ULIN(0),1) UWEB(2,NM,NC) = DDOT(64,DELTA(0,NPT),1,VLIN(0),1) UWEB(3,NM,NC) = DDOT(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 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=8,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 C$OMP PARALLEL DO C$OMP1SHARED(XFN, STF0, REST0, FRC, LAFLAG, ARFLAG, NFG, NPF, C$OMP2 STF , REST , UNSTBL, C$OMP3 ISTART, ISTOP, KGROUP) C$OMP4PRIVATE(I, K, KOUNT) C$OMP+SCHEDULE(RUNTIME) 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 C$OMP PARALLEL DO C$OMP1SHARED(XFN, STF0, REST0, FRC, LAFLAG, ARFLAG, NFG, NPF, C$OMP2 XF0, STF , REST , C$OMP3 ISTART, ISTOP, KGROUP, KGRPS) C$OMP4PRIVATE(I, K, K0, KOUNT) C$OMP+SCHEDULE(RUNTIME) 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 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= 0. b is returned if c < 0 real, intent(in) :: a,b,c cvmgp = b if (c>=0.) cvmgp=a end function cvmgp real(8) function ddot(n,x,incx,y,incy) !Computes a dot product (inner product) of two real vectors integer, intent(in) :: n,incx,incy real(8), intent(in) :: x((n-1)*abs(incx)+1),y((n-1)*abs(incy)+1) integer :: i ddot=0. do i=1,n,incx ddot=x(i)*y(i)+ddot end do end function ddot integer function idamax(n,x,incx) !IDAMAX Searches a vector for the first occurrence of the maximum absolute value integer, intent(in) :: n,incx real(8), intent(in) :: x((n-1)*abs(incx)+1) integer :: i real :: value value=0. do i=1,n,incx if (abs(x(i))>value) then value=abs(x(i)) idamax=i end if end do end function idamax