c======================================================================| subroutine sub_slct(fnin,fnout,fntmp2,angref6,nang 1 ,touts,dtr,idprop,ifnpo) c======================================================================| implicit none integer,intent(in) :: nang,ifnpo real*8,intent(in) :: angref6(nang),dtr,touts character*100 ,intent(in):: fntmp2(nang),fnin,fnout real*8,allocatable :: jsout(:,:),ddout(:,:,:),jstmp2(:,:) real*8 :: angref0,para(13),js1,sec0 integer :: i,j,ifn0,ifn1,ifn2,iang,nstep integer :: iy0,im0,id0,jdd,jdd2,ihr0,imin0,isec0 integer :: nwout,nnin real*8 :: dtout,jswout0,jswout1,dtin real*8 :: pdy0 real*8,allocatable :: jswout(:),ddwout(:,:,:),ddwout1(:,:) real*8,allocatable :: datain(:,:),jsin(:),jsinref(:,:) real*8,allocatable :: dang2(:,:),dang1(:,:) real*8,allocatable :: dangi(:),dango(:),ang0(:) real*8,allocatable :: jsoutrev(:,:),ddoutrev(:,:,:) integer :: idprop real*8,allocatable :: angwin(:),angwout(:),dangw(:),angw(:) real*8,allocatable :: jstmp(:) real*8,allocatable :: dangi2(:) real*8 :: jstmp0 ,dangtmp0,bd1,bd2 integer :: idstep,idstep2 integer,allocatable :: islct(:) real*8,allocatable :: dangtmp(:,:,:),ddslct(:,:),jsslct(:) real*8,allocatable :: jsinref2(:) ifn1=43 ifn2=44 !---read input data call count_lines(fnin,nnin) allocate(datain(13,nnin),jsin(nnin),dangi(nnin),dangi2(nnin)) open(ifn2,file=fnin,status='old') do i=1,nnin read(ifn2,*) js1,(para(j),j=1,13) datain(:,i)=para(1:13) jsin(i)=js1 enddo dtin=jsin(2)-jsin(1) close(ifn2,status='keep') !---count line & read simulation data iang=1 angref0=angref6(iang) call count_lines(fntmp2(iang),nstep) allocate(jsout(nang,nstep),ddout(nang,10,nstep) 1 ,jsinref(nang,nstep),ang0(nstep),dango(nstep) 1 ,dang1(nang,nstep),dang2(nang,nstep),jsoutrev(nang,nstep)) allocate(jstmp(nstep)) allocate(ddslct(10,nstep),jsslct(nstep)) allocate(islct(nstep),dangtmp(2,nang,nstep),jsinref2(nstep)) allocate(jstmp2(6,nstep)) do iang=1,nang angref0=angref6(iang) open(ifn1,file=fntmp2(iang),status='old') do i=1,nstep read(ifn1,*) js1,(para(j),j=1,10) jsout(iang,i)=js1 ddout(iang,:,i)=para(1:10) enddo close(ifn1,status='keep') enddo jsinref(:,:)=ddout(:,10,:) write(ifnpo,*) nnin,nstep c---time shift for radial difference jsoutrev(:,:)=jsout(:,:) c---angle evaluation do iang=1,nang !---6 ref.angles angref0=angref6(iang) dangi(:)=datain(10,:)-angref0 !--input-ref dangi2(1)=dangi(1) do i=2,nnin dangtmp0=dangi(i)-dangi(i-1) if(abs(dangtmp0).ge.300.) dangtmp0=0. dangi2(i)=dangi2(i-1)+dangtmp0 enddo bd1=dangi2(1); bd2=dangi2(nnin) call convim(jsin,jsinref(iang,:),dangi2,dang1(iang,:) 1 ,nnin,nstep,bd1,bd2) dangi(:)=datain(12,:) !--input-ref dangi2(1)=dangi(1) do i=2,nnin dangtmp0=dangi(i)-dangi(i-1) if(abs(dangtmp0).ge.300.) dangtmp0=0. dangi2(i)=dangi2(i-1)+dangtmp0 enddo bd1=dangi2(1); bd2=dangi2(nnin) call convim(jsin,jsoutrev(iang,:),dangi2,ang0(:) 1 ,nnin,nstep,bd1,bd2) dango(:)=ang0(:)-angref0 !--output-ref c---planeatry rotation effect (2nd step) do idstep=1,3 call angrev(dango,dango,nstep) jstmp=jsoutrev(iang,:)+dango/360.*26.*86400 do i=2,nstep if(abs(jstmp(i)-jstmp(i-1)).gt.(maxval(jstmp)-minval(jstmp)) 1 /float(nstep)*5.) jstmp(i)=jstmp(i-1) enddo bd1=dangi2(1); bd2=dangi2(nnin) call convim(jsin,jstmp,dangi2,ang0(:) 1 ,nnin,nstep,bd1,bd2) dango(:)=ang0(:)-angref0 !--output-ref enddo dang2(iang,:)=dango(:) enddo c---angle selection do iang=1,nang call angrev(dang1(iang,:),dang1(iang,:),nstep) call angrev(dang2(iang,:),dang2(iang,:),nstep) enddo dangtmp(2,:,:)=abs(dang1(:,:))+abs(dang2(:,:)) dangtmp(1,:,:)=dang1(:,:)+dang2(:,:) do iang=1,nang idstep=0 do i=2,nstep-1 jstmp2(1,i)=jsinref(iang,i+1)-jsinref(iang,i-1) if((idprop.eq.-1.and.jsout(1,i).ge.86400.*20.).or. 1 (idprop.eq.1.and.jsout(1,i).le.-86400.*20.)) then idstep=idstep+1 jstmp2(2,idstep)=jsinref(iang,i) jstmp2(3,idstep)=float(i) endif enddo if(idstep.ge.1)then do i=1,idstep/2. jstmp2(4,i)=minval(minloc( 1 abs(jsinref(iang,1:int(jstmp2(3,i)-1))-jstmp2(2,i*2)))) dangtmp(2,iang,int(jstmp2(4,i)):int(jstmp2(3,i*2)))=901. enddo endif do i=2,nstep-1 jstmp2(5,i)=dang2(iang,i+1)-dang2(iang,i-1) if(abs(jstmp2(5,i)).lt.1.e-10) dangtmp(2,iang,i)=902. if((idprop.eq.-1.and.ddout(iang,10,i).ge.jsin(1)) 1 .or.(idprop.eq.1.and.ddout(iang,10,i).le.jsin(1))) 1 dangtmp(2,iang,i)=999. enddo enddo do i=1,nstep islct(i)=minval(minloc(dangtmp(2,:,i))) dango(i)=dangtmp(1,islct(i),i) jsslct(i)=jsoutrev(islct(i),i)+ddout(islct(i),9,i) 1 +dang2(islct(i),i)/360.*26.*86400. ddslct(:,i)=ddout(islct(i),:,i) enddo c---set time dtout=touts*dtr call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(:,1))) call ymds2js(iy0,im0,id0+int((idprop+1)*0.5),0.d0,jswout0) call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(:,nstep))) call ymds2js(iy0,im0,id0+int((idprop-1)*0.5),0.d0,jswout1) nwout=abs(jswout1-jswout0)/dtout write(ifnpo,*) nwout allocate(jswout(nwout),ddwout(nang,10,nwout) 1 ,ddwout1(11,nwout)) do i=1,nwout if(idprop.eq.1) jswout(i)=jswout0+dble(i)*dtout !--REV0318 if(idprop.eq.-1) jswout(i)=jswout1+dble(i)*abs(dtout) !--REV0318 enddo c---data and angle conversion do i=1,10 bd1=ddslct(i,1); bd2=ddslct(i,nstep) if(idprop.eq.-1)then bd1=ddslct(i,nstep); bd2=ddslct(i,1) endif call convim(jsslct(:),jswout(:),ddslct(i,:),ddwout1(i,:) 1 ,nstep,nwout,bd1,bd2) enddo c---data lack call convim(jsin,ddwout1(10,:),datain(13,:),ddwout1(11,:) 1 ,nnin,nwout,datain(13,1),datain(13,nnin)) do i=2,nwout if(abs(ddwout1(3,i)-ddwout1(3,i-1)).le.abs(ddwout1(3,i))*1.e-10) 1 ddwout1(11,i)=1. enddo c---angle evaluation: in-out allocate(angwin(nwout),angwout(nwout),dangw(nwout),angw(nwout)) dangi(:)=datain(10,:) !--input-ref dangi2(1)=dangi(1) do i=2,nnin dangtmp0=dangi(i)-dangi(i-1) if(abs(dangtmp0).ge.300.) dangtmp0=0. dangi2(i)=dangi2(i-1)+dangtmp0 enddo dangw=ddwout1(10,:) do i=101,nwout-100 dangw(i)=sum(ddwout1(10,i-100:i+100))/201. enddo bd1=dangi2(1); bd2=dangi2(nnin) if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif call convim(jsin,dangw,dangi2,angwin(:) 1 ,nnin,nwout,bd1,bd2) dangi(:)=datain(12,:) !--output-ref dangi2(1)=dangi(1) do i=2,nnin dangtmp0=dangi(i)-dangi(i-1) if(abs(dangtmp0).ge.300.) dangtmp0=0. dangi2(i)=dangi2(i-1)+dangtmp0 enddo bd1=dangi2(1); bd2=dangi2(nnin) if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif call convim(jsin,jswout,dangi2,angwout(:) 1 ,nnin,nwout,bd1,bd2) dangw=angwin-angwout !--difference call angrev(dangw,dangw,nwout) idstep=0; idstep2=0 do i=2,nwout angw(i)=dangw(i)-dangw(i-1) if(angw(i).ge.300.) then idstep=idstep+1 jstmp2(2,idstep)=jswout(i) jstmp2(3,idstep)=float(i) endif if(angw(i).le.-300.) then idstep2=idstep2+1 jstmp2(4,idstep2)=jswout(i) jstmp2(5,idstep2)=float(i) endif enddo if(idprop.eq.1.and.idstep.ge.1.and.idstep*2.eq.idstep2)then dangw(int(jstmp2(5,idstep*2-1)):int(jstmp2(3,idstep))) 1 =dangw(int(jstmp2(5,idstep*2-1)):int(jstmp2(3,idstep)))+360. dangw(int(jstmp2(3,idstep)):int(jstmp2(5,idstep*2))) 1 =dangw(int(jstmp2(3,idstep)):int(jstmp2(5,idstep*2)))-360. c 1 =dangw(int(jstmp2(3,idstep))-1:int(jstmp2(5,idstep*2)))-360. ddwout1(11,int(jstmp2(5,idstep*2-1)):int(jstmp2(5,idstep*2)))=1. endif if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then dangw(int(jstmp2(3,idstep2*2-1)):int(jstmp2(5,idstep2))) 1 =dangw(int(jstmp2(3,idstep2*2-1)):int(jstmp2(5,idstep2)))-360. c dangw(int(jstmp2(5,idstep2)):int(jstmp2(3,idstep2*2))) dangw(int(jstmp2(5,idstep2))-1:int(jstmp2(3,idstep2*2))) 1 =dangw(int(jstmp2(5,idstep2))-1:int(jstmp2(3,idstep2*2)))+360. endif ddwout1(9,:)=dangw !separation angle c---output open(ifn2,file=fnout,status='unknown',form='formatted') do i=1,nwout,1 call js2ymds(iy0,im0,id0,sec0,jswout(i)) ihr0=int(sec0/3600.) imin0=int((sec0-ihr0*3600.)/60.) isec0=int(sec0-(ihr0*3600.+imin0*60.)) pdy0=1.67e-27*ddwout1(1,i)*1.e6 1 *(ddwout1(3,i)**2+ddwout1(4,i)**2)*1.e6 write(ifn2,'(i4.4,5(a1,i2.2),12(1pe13.5))') 1 iy0,' ',im0,' ',id0,' ',ihr0,':',imin0,':',isec0 1 ,ddwout1(1:4,i),ddwout1(7,i) 2 ,pdy0*1.e9,ddwout1(9,i),ddwout1(11,i) c write(62,'(i12,12e18.9)') int(jswout(i)),(ddwout1(j,i),j=1,11) enddo c---finish deallocate(datain,jsin,dangi,dangi2 1 ,jsout,ddout,jsinref,ang0,dango,dang1,jsoutrev 1 ,jswout,ddwout,ddwout1,jstmp2) return end c==================================================