sub_slct_p.f 9.75 KB
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==================================================