program peskin_tecplot ! to read in *.fl *.mk *.vp *.xf files from peskin's code output ! now >256 marker locations are ignored implicit none character(30) :: ZoneTitle(0:35) ! fl related integer, parameter :: valves=4,rad_arms=98,web_cir=24 integer :: klok,valve_no,rad_arms_no,web_cir_no,status,tmp_i1,tmp_i2,i,j,k character(10) :: tmp_c1 real :: x_fl(rad_arms+1,web_cir,valves),y_fl(rad_arms+1,web_cir,valves),z_fl(rad_arms+1,web_cir,valves) ! +1 to close circle real :: u_fl(rad_arms+1,web_cir,valves),v_fl(rad_arms+1,web_cir,valves),w_fl(rad_arms+1,web_cir,valves),p_fl(rad_arms+1,web_cir,valves) real :: x0(valves),y0(valves),z0(valves),u0(valves),v0(valves),w0(valves) ! mk related integer, parameter :: clouds=19 integer :: ng,mk_in_cloud(clouds),cloud_no integer :: max_mk_in_cloud,marker,size_x,size_y,size_z,mii,mjj,mkk integer, allocatable :: mi(:,:),mj(:,:),mk(:,:) ! indicate which cube marker belongs to character(25), allocatable :: tmp_c2(:,:,:) character(16), allocatable :: tmp_c3(:,:,:) character :: tmp_c4(29),tmp_c5(17),tmp_c6(11),tmp_c7(5) real :: cube_length,length,xx,yy,zz real(8) :: time real, allocatable :: marker_x(:,:),marker_y(:,:),marker_z(:,:),marker_u(:,:),marker_v(:,:),marker_w(:,:),marker_p(:,:) real, allocatable :: x_mk(:,:,:),y_mk(:,:,:),z_mk(:,:,:),u_mk(:,:,:),v_mk(:,:,:),w_mk(:,:,:),p_mk(:,:,:) ! xf related integer :: layers,layer_no,gp_no,fibre_no,fibres,pts,pt_no integer, allocatable :: gp_in_layer(:),grp(:,:),nfg(:,:),npf(:,:),fibre_in_layer(:),pt_in_layer(:),pt_in_layer2(:) real, allocatable :: x_xf(:,:),y_xf(:,:),z_xf(:,:),dummy_xf(:,:) character(20) :: hrt_fl,hrt_mk,hrt_vp,hrt_xf,hrt_out,hrt_no integer :: starting,stopping,interval,file_no ! Pls enter start,end timestep and interval starting=83456; stopping=88064; interval=512 !starting=55296; stopping=55296; interval=512 !write (tmp_c1,'(i10)') starting !tmp_c1=adjustl(tmp_c1) !arg_buffer(1)='hrt'//trim(tmp_c1)//'.txt' do file_no=starting,stopping,interval write (hrt_no,'(i10)') file_no hrt_no=adjustl(hrt_no) !hrt_out='hrt'//trim(hrt_no)//'.plt' !if filename is hrt512.plt if (file_no<1000) then hrt_out='hrt_000'//trim(hrt_no)//'.plt' !if filename is hrt000512.plt hrt_fl='hrt1113_000'//trim(hrt_no)//'.fl'; hrt_mk='hrt1113_000'//trim(hrt_no)//'.mk' hrt_vp='hrt1113_000'//trim(hrt_no)//'.vp'; hrt_xf='hrt1113_000'//trim(hrt_no)//'.xf' else if (file_no<10000) then hrt_out='hrt_00'//trim(hrt_no)//'.plt' hrt_fl='hrt1113_00'//trim(hrt_no)//'.fl'; hrt_mk='hrt1113_00'//trim(hrt_no)//'.mk' hrt_vp='hrt1113_00'//trim(hrt_no)//'.vp'; hrt_xf='hrt1113_00'//trim(hrt_no)//'.xf' else if (file_no<100000) then hrt_out='hrt_0'//trim(hrt_no)//'.plt' hrt_fl='hrt1113_0'//trim(hrt_no)//'.fl'; hrt_mk='hrt1113_0'//trim(hrt_no)//'.mk' hrt_vp='hrt1113_0'//trim(hrt_no)//'.vp'; hrt_xf='hrt1113_0'//trim(hrt_no)//'.xf' else hrt_out='hrt_'//trim(hrt_no)//'.plt' hrt_fl='hrt1113_'//trim(hrt_no)//'.fl'; hrt_mk='hrt1113_'//trim(hrt_no)//'.mk' hrt_vp='hrt1113_'//trim(hrt_no)//'.vp'; hrt_xf='hrt1113_'//trim(hrt_no)//'.xf' end if ! fl related ZoneTitle(32)="Aortic valve"; ZoneTitle(33)="Mitral valve" ZoneTitle(34)="Tricuspid valve"; ZoneTitle(35)="Pulmonic valve" open (unit = 1 , file = hrt_fl , form='formatted', status = "old", iostat = status) if (status/=0) print *, "file not found" klok=file_no read (1,*) tmp_c1 !read (1,*) klok read (1,*) tmp_c1 read (1,*) tmp_i1, tmp_c1 do i=1,valves read (1,*) tmp_i1, tmp_i2 ! values already known, so just scroll thru end do do valve_no=1,valves read (1,*) tmp_i1 read (1,*) read (1,*) u0(valve_no),v0(valve_no),w0(valve_no),x0(valve_no),y0(valve_no),z0(valve_no) ! meant for web circle 0 !read (1,*) u0(valve_no),v0(valve_no),tmp_c4,tmp_c5,tmp_c6,tmp_c7 ! meant for web circle 0 !if (len(tmp_c4)==29) then ! read ( tmp_c4(1:11), '(e11.11)' ) w0(valve_no) ! read ( tmp_c4(12:17), '(e6.6)' ) x0(valve_no) !read ( tmp_c4(18:23), '(e8.8)' ) y0(valve_no) !read ( tmp_c4(24:29), '(e8.8)' ) z0(valve_no) !end if do web_cir_no=1,web_cir do rad_arms_no=1,rad_arms read (1,*) u_fl(rad_arms_no,web_cir_no,valve_no),v_fl(rad_arms_no,web_cir_no,valve_no), & & w_fl(rad_arms_no,web_cir_no,valve_no),x_fl(rad_arms_no,web_cir_no,valve_no),y_fl(rad_arms_no,web_cir_no,valve_no), & & z_fl(rad_arms_no,web_cir_no,valve_no) !tmp_c4 !read ( tmp_c4(1:11), '(e11.11)' ) w_fl(rad_arms_no,web_cir_no,valve_no) !read ( tmp_c4(12:17), '(e6.6)' ) x_fl(rad_arms_no,web_cir_no,valve_no) !read ( tmp_c4(18:23), '(e8.8)' ) y_fl(rad_arms_no,web_cir_no,valve_no) !read ( tmp_c4(24:29), '(e8.8)' ) z_fl(rad_arms_no,web_cir_no,valve_no) end do end do end do close (1) u_fl(rad_arms+1,:,:)=u_fl(1,:,:); v_fl(rad_arms+1,:,:)=v_fl(1,:,:); w_fl(rad_arms+1,:,:)=w_fl(1,:,:) x_fl(rad_arms+1,:,:)=x_fl(1,:,:); y_fl(rad_arms+1,:,:)=y_fl(1,:,:); z_fl(rad_arms+1,:,:)=z_fl(1,:,:) ! mk related ZoneTitle(13)="Aortic ring"; ZoneTitle(14)="Mitral ring"; ZoneTitle(15)="Tricuspid ring"; ZoneTitle(16)="Pulmonic ring" ZoneTitle(17)="Left vent"; ZoneTitle(18)="Right vent"; ZoneTitle(19)="SVC source"; ZoneTitle(20)="IVC source" ZoneTitle(21)="PV source"; ZoneTitle(22)="PA source"; ZoneTitle(23)="Aorta source"; ZoneTitle(24)="Left atrium" ZoneTitle(25)="Right atrium"; ZoneTitle(26)="Aorta artery"; ZoneTitle(27)="Pulmonary artery"; ZoneTitle(28)="Pressure tap" ZoneTitle(29)="Marker 17"; ZoneTitle(30)="Marker 18"; ZoneTitle(31)="Marker 19" open (unit = 1 , file = hrt_mk , form='formatted', status = "old", iostat = status) if (status/=0) print *, "file not found" read (1,*) tmp_c1 read (1,*) ng read (1,*) tmp_i1 do cloud_no=1,clouds read (1,*) mk_in_cloud(cloud_no) end do !read (1,*) tmp_i1 read (1,*) tmp_c1 read (1,*) time if (file_no==starting) then max_mk_in_cloud=maxval(mk_in_cloud) allocate (marker_x(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_y(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_z(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (mi(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (mj(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (mk(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_u(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_v(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_w(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (marker_p(max_mk_in_cloud,clouds), STAT=status) if (status/=0) STOP "Cannot allocate memory" !size_x=ng/4;size_y=ng/4;size_z=ng/4 size_x=ng/8;size_y=ng/8;size_z=ng/8 allocate (x_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (y_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (z_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (u_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (v_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (w_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (p_mk(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (tmp_c2(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (tmp_c3(size_x,size_y,size_z), STAT=status) if (status/=0) STOP "Cannot allocate memory" end if marker_x=0.;marker_y=0.;marker_z=0.;marker_u=0.;marker_v=0.;marker_w=0.;marker_p=0. mi=0;mj=0;mk=0 x_mk=0.;y_mk=0.;z_mk=0.;u_mk=0.;v_mk=0.;w_mk=0.;p_mk=0. do cloud_no=1,clouds do marker=1,mk_in_cloud(cloud_no) read (1,*) marker_x(marker,cloud_no),marker_y(marker,cloud_no),marker_z(marker,cloud_no) end do end do close (1) ! vp section open (unit = 1 , file = hrt_vp , form='formatted', status = "old", iostat = status) if (status/=0) print *, "file not found" cube_length=64.*10./(2.*3.1415926536*5.85); length=cube_length/ng !read (1,*), tmp_i1 read (1,*), tmp_c1 read (1,*) do k=1,size_z do j=1,size_y do i=1,size_x read (1,*) u_mk(i,j,k),v_mk(i,j,k),w_mk(i,j,k),tmp_c2(i,j,k) tmp_c3(i,j,k)=tmp_c2(i,j,k) end do end do end do close(1) do k=1,size_z do j=1,size_y do i=1,size_x read ( tmp_c3(i,j,k), '(e16.16)' ) p_mk(i,j,k) end do end do end do do k=1,size_z do j=1,size_y do i=1,size_x !x(i,j,k)=length*i; y(i,j,k)=length*j; z(i,j,k)=length*(k-1) x_mk(i,j,k)=4.*i; y_mk(i,j,k)=4.*j; z_mk(i,j,k)=4.*(k-1) end do end do end do do cloud_no=1,clouds do marker=1,mk_in_cloud(cloud_no) xx=marker_x(marker,cloud_no) yy=marker_y(marker,cloud_no) zz=marker_z(marker,cloud_no) if (xx>256. .or. yy>256. .or. zz>256. .or. xx<4. .or. yy<4. .or. zz<4.) then print *, xx,yy,zz marker_u(marker,cloud_no)=0. marker_v(marker,cloud_no)=0. marker_w(marker,cloud_no)=0. marker_p(marker,cloud_no)=0. else mi(marker,cloud_no)=int(xx) ! guess mj(marker,cloud_no)=int(yy) ! guess mk(marker,cloud_no)=int(zz) ! guess call hunt(x_mk(:,1,1),size_x,xx,mi(marker,cloud_no)) call hunt(y_mk(1,:,1),size_y,yy,mj(marker,cloud_no)) call hunt(z_mk(1,1,:),size_z,zz,mk(marker,cloud_no)) mi(marker,cloud_no)=mi(marker,cloud_no)/2 !adjust for ng=256 mj(marker,cloud_no)=mj(marker,cloud_no)/2 mk(marker,cloud_no)=mk(marker,cloud_no)/2 mii=mi(marker,cloud_no);mjj=mj(marker,cloud_no);mkk=mk(marker,cloud_no) do k=0,1 do j=0,1 do i=0,1 marker_u(marker,cloud_no)=marker_u(marker,cloud_no)+0.125*u_mk(mii+i,mjj+j,mkk+k) marker_v(marker,cloud_no)=marker_v(marker,cloud_no)+0.125*v_mk(mii+i,mjj+j,mkk+k) marker_w(marker,cloud_no)=marker_w(marker,cloud_no)+0.125*w_mk(mii+i,mjj+j,mkk+k) marker_p(marker,cloud_no)=marker_p(marker,cloud_no)+0.125*p_mk(mii+i,mjj+j,mkk+k) end do end do end do end if end do end do ! xf related ZoneTitle(0)="Unknown layer" ZoneTitle(1)="Outer/Inner layer"; ZoneTitle(2)="RV IN/LV OUT layer"; ZoneTitle(3)="Lv cyl largest layer"; ZoneTitle(4)="Lv cyl layer" ZoneTitle(5)="Lv cyl layer"; ZoneTitle(6)="Lv cyl smallest layer"; ZoneTitle(7)="Aortic valve layer"; ZoneTitle(8)="Mitral valve layer" ZoneTitle(9)="Tricuspid valve layer"; ZoneTitle(10)="Pulmonic valve layer"; ZoneTitle(11)="Left atrium layer"; ZoneTitle(12)="Right atrium layer" open (unit = 1 , file = hrt_xf , form='formatted', status = "old", iostat = status) if (status/=0) print *, "file not found" read (1,*) tmp_i1 read (1,*) tmp_c1, layers if (file_no==starting) then allocate (gp_in_layer(0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (fibre_in_layer(0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (pt_in_layer(0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (pt_in_layer2(0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (grp(31,0:layers), STAT=status) ! 31 cos is max gp in layer, should be fixed in this case if (status/=0) STOP "Cannot allocate memory" allocate (nfg(31,0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (npf(31,0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" end if gp_in_layer=0; fibre_in_layer=0; pt_in_layer=0; pt_in_layer2=0 grp=0; nfg=0; npf=0 do layer_no=0,layers read (1,*) tmp_c1, gp_in_layer(layer_no), tmp_i1 do gp_no=1,gp_in_layer(layer_no) read (1,*) tmp_c1, grp(gp_no,layer_no), nfg(gp_no,layer_no), npf(gp_no,layer_no) fibre_in_layer(layer_no)=fibre_in_layer(layer_no)+nfg(gp_no,layer_no) pt_in_layer(layer_no)=pt_in_layer(layer_no)+nfg(gp_no,layer_no)*npf(gp_no,layer_no) end do end do if (file_no==starting) then allocate (x_xf(maxval(pt_in_layer),0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (y_xf(maxval(pt_in_layer),0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (z_xf(maxval(pt_in_layer),0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" allocate (dummy_xf(maxval(pt_in_layer),0:layers), STAT=status) if (status/=0) STOP "Cannot allocate memory" end if x_xf=0.; y_xf=0.; z_xf=0.; dummy_xf=0. read (1,*) tmp_c1 !read (1,*) tmp_i1 read (1,*) tmp_c1 read (1,*) tmp_i1 read (1,*) fibres if (fibres==sum(fibre_in_layer)) then print *, "total fibres tally" else print *, "total fibres do not tally" stop end if do layer_no=0,layers do fibre_no=1,fibre_in_layer(layer_no) read (1,*) tmp_i1, pts do pt_no=1,pts ! last no. repeat pt_in_layer2(layer_no)=pt_in_layer2(layer_no)+1 read (1,*) x_xf(pt_in_layer2(layer_no),layer_no),y_xf(pt_in_layer2(layer_no),layer_no),z_xf(pt_in_layer2(layer_no),layer_no) end do end do end do close (1) do layer_no=0,layers if (pt_in_layer2(layer_no)/=pt_in_layer(layer_no)) then print *, "total pts in layer", layer_no, "do not tally" stop end if end do print *, "total pts in layer tally" !call tecplot_fl("valves.plt") !call tecplot_mk_vp("markers_uvw.plt") !call tecplot_xf("fibres.plt") call tecplot_xf_mk_vp_fl(hrt_out) end do deallocate (mi, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (mj, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (tmp_c2, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (tmp_c3, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_x, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_y, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_z, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_u, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_v, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_w, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (marker_p, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (x_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (y_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (z_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (u_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (v_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (w_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (p_mk, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (gp_in_layer, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (grp, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (nfg, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (npf, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (fibre_in_layer, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (pt_in_layer, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (pt_in_layer2, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (x_xf, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (y_xf, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (z_xf, STAT=status) if (status/=0) STOP "Cannot deallocate memory" deallocate (dummy_xf, STAT=status) if (status/=0) STOP "Cannot deallocate memory" contains subroutine tecplot_xf_mk_vp_fl(filename) !to output node value binary tecplot file implicit none include 'tecio.f90' integer :: i,j,imax,jmax,kmax,debug,ier,itot,visdouble,disdouble,ICellMax,JCellMax,KCellMax,Filetype integer :: ZoneType,StrandID,ParentZn,IsBlock,NFConns,FNMode,ShrConn,varlist(7),varloc(7),sharevar(7) integer :: TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBoundaryConnections !integer :: clouds,valves character*1 :: nulchar character*20, intent(in) :: filename nulchar = char(0); debug = 0; visdouble = 0; disdouble = 0; varloc=1; sharevar=0 imax = 1; jmax = 1; kmax = 1; ICellMax = 0; JCellMax = 0; KCellMax = 0; Filetype = 0 ZoneType = 0; StrandID = 0; ParentZn = 0 ;IsBlock = 1; NFConns = 0; FNMode = 0; ShrConn = 0 TotalNumFaceNodes = 1; TotalNumBndryFaces = 1; TotalNumBoundaryConnections = 1 varlist(1:3)=0; varlist(4:7)=1 !varlist=1 !clouds=19; valves=4 ! xf related ! Open the file and write the tecplot datafile header information. ! !ier = tecini('flowfield'//nulchar,'x,y,u_node,v_node'//nulchar,'output.plt'//nulchar,'.'//nulchar,debug,visdouble) !ier = TecIni111('Flowfield'//nulchar,'x,y,z,u,v,w,p'//NULCHaR,filename//nulchar,'.'//nulchar,Filetype,debug,disdouble) ier = TecIni142('Flowfield'//nulchar,'x,y,z,u,v,w,p'//NULCHaR,filename//nulchar,'.'//nulchar,1,Filetype,debug,disdouble) do layer_no=0,layers ! Write the zone header information. StrandID = layer_no+1 !ParentZn = layer_no imax = pt_in_layer(layer_no) ier = TecZne142(ZoneTitle(layer_no)//nulchar, ZoneType,IMax,JMax,KMax,ICellMax,JCellMax,KCellMax,time,StrandID,ParentZn, & IsBlock,NFConns,FNMode,TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBoundaryConnections,varlist,varloc,sharevar,ShrConn) ! Write out the field data. itot = (imax)*(jmax)*kmax ier = tecdat142(itot,real(x_xf(:,layer_no)),disdouble) ier = tecdat142(itot,real(y_xf(:,layer_no)),disdouble) ier = tecdat142(itot,real(z_xf(:,layer_no)),disdouble) ! ier = tecdat142(itot,real(dummy_xf(:,layer_no)),disdouble) ! ier = tecdat142(itot,real(dummy_xf(:,layer_no)),disdouble) ! ier = tecdat142(itot,real(dummy_xf(:,layer_no)),disdouble) ! ier = tecdat142(itot,real(dummy_xf(:,layer_no)),disdouble) end do ! mk related varlist=0 do cloud_no=1,clouds ! Write the zone header information. StrandID = cloud_no+layers+1 !ParentZn = cloud_no+layers !StrandID = 1 imax = mk_in_cloud(cloud_no) ier = TecZne142(ZoneTitle(StrandID-1)//nulchar, ZoneType,IMax,JMax,KMax,ICellMax,JCellMax,KCellMax,time,StrandID,ParentZn, & IsBlock,NFConns,FNMode,TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBoundaryConnections,varlist,varloc,sharevar,ShrConn) ! Write out the field data. itot = mk_in_cloud(cloud_no) ier = tecdat142(itot,real(marker_x(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_y(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_z(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_u(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_v(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_w(1:itot,cloud_no)),disdouble) ier = tecdat142(itot,real(marker_p(1:itot,cloud_no)),disdouble) end do ! fl related imax = rad_arms+1; jmax = web_cir; kmax = 1 do valve_no=1,valves ! Write the zone header information. StrandID = valve_no+clouds+layers+1 !ParentZn = valve_no+clouds+layers ier = TecZne142(ZoneTitle(StrandID-1)//nulchar, ZoneType,IMax,JMax,KMax,ICellMax,JCellMax,KCellMax,time,StrandID,ParentZn, & IsBlock,NFConns,FNMode,TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBoundaryConnections,varlist,varloc,sharevar,ShrConn) ! Write out the field data. itot = (imax)*(jmax)*kmax ier = tecdat142(itot,real(x_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(y_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(z_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(u_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(v_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(w_fl(:,:,valve_no)),disdouble) ier = tecdat142(itot,real(p_fl(:,:,valve_no)),disdouble) end do ! vp related !imax = ng/4; jmax = ng/4; kmax = ng/4 imax = ng/8; jmax = ng/8; kmax = ng/8 ! Write the zone header information. StrandID = valves+clouds+layers+2 !ParentZn = valve_no+clouds+layers ier = TecZne142("UVWP field"//nulchar, ZoneType,IMax,JMax,KMax,ICellMax,JCellMax,KCellMax,time,StrandID,ParentZn, & IsBlock,NFConns,FNMode,TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBoundaryConnections,varlist,varloc,sharevar,ShrConn) ! Write out the field data. itot = (imax)*(jmax)*kmax ier = tecdat142(itot,real(2.*x_mk),disdouble) ier = tecdat142(itot,real(2.*y_mk),disdouble) ier = tecdat142(itot,real(2.*z_mk),disdouble) ier = tecdat142(itot,real(u_mk),disdouble) ier = tecdat142(itot,real(v_mk),disdouble) ier = tecdat142(itot,real(w_mk),disdouble) ier = tecdat142(itot,real(p_mk),disdouble) ! close file ier = tecend142() end subroutine tecplot_xf_mk_vp_fl SUBROUTINE hunt(xx,n,x,jlo) INTEGER jlo,n REAL x,xx(n) !Given an array xx(1:n), and given a value x, returns a value jlo such that x is between !xx(jlo) and xx(jlo+1). xx(1:n) must be monotonic, either increasing or decreasing. !jlo=0 or jlo=n is returned to indicate that x is out of range. jlo on input is taken as !the initial guess for jlo on output. INTEGER inc,jhi,jm LOGICAL ascnd ascnd=xx(n).ge.xx(1) !True if ascending order of table, false otherwise. if(jlo.le.0.or.jlo.gt.n)then !Input guess not useful. Go immediately to bisection. jlo=0 jhi=n+1 goto 3 endif inc=1 !Set the hunting increment. if(x.ge.xx(jlo).eqv.ascnd)then !Hunt up: 1 jhi=jlo+inc if(jhi.gt.n)then !Done hunting, since off end of table. jhi=n+1 else if(x.ge.xx(jhi).eqv.ascnd)then !Not done hunting, jlo=jhi inc=inc+inc !so double the increment goto 1 !and try again. endif !Done hunting, value bracketed. else !Hunt down: jhi=jlo 2 jlo=jhi-inc if(jlo.lt.1)then !Done hunting, since off end of table. jlo=0 else if(x.lt.xx(jlo).eqv.ascnd)then !Not done hunting, jhi=jlo inc=inc+inc !so double the increment goto 2 !and try again. endif !Done hunting, value bracketed. endif !Hunt is done, so begin the final bisection phase: 3 if(jhi-jlo.eq.1)then if(x.eq.xx(n))jlo=n-1 if(x.eq.xx(1))jlo=1 return endif jm=(jhi+jlo)/2 if(x.ge.xx(jm).eqv.ascnd)then jlo=jm else jhi=jm endif goto 3 END SUBROUTINE hunt end program peskin_tecplot