return to old selection with modifs

Elena.Budnik
1 parent af16ae5a
Showing 3 changed files with 627 additions and 38 deletions Show diff stats
Sources/Prop_code/sub_slct_p.f
Sources/Prop_code/sub_slct_p.f.Tao.new
Sources/Prop_code/sub_slct_p.f.my
@@ -6,7 +6,8 @@ c======================================================================|
       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,allocatable :: jsout(:,:),ddout(:,:,:)
+!!!! elena : para(13) !!! not para(12)
       real*8 :: angref0,para(13),js1,sec0
       integer :: i,j,ifn0,ifn1,ifn2,iang,nstep
       integer :: iy0,im0,id0,jdd,jdd2,ihr0,imin0,isec0
@@ -30,7 +31,6 @@ c======================================================================|
  
       ifn1=43
       ifn2=44
-
 !---read input data
       call count_lines(fnin,nnin)
       allocate(datain(13,nnin),jsin(nnin),dangi(nnin),dangi2(nnin))
@@ -47,13 +47,13 @@ c======================================================================|
       iang=1
       angref0=angref6(iang)
       call count_lines(fntmp2(iang),nstep)
-      allocate(jsout(nang,nstep),ddout(nang,10,nstep)
+!!!! elena :   jsout(nang,nstep),  jsoutrev(nang,nstep)   => jsout(6,nstep), jsoutrev(6,nstep)
+      allocate(jsout(6,nstep),ddout(nang,10,nstep)
      1  ,jsinref(nang,nstep),ang0(nstep),dango(nstep)
-     1  ,dang1(nang,nstep),dang2(nang,nstep),jsoutrev(nang,nstep))
+     1  ,dang1(nang,nstep),dang2(nang,nstep),jsoutrev(6,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)
@@ -85,7 +85,8 @@ c---angle evaluation
         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
@@ -96,6 +97,7 @@ c---angle evaluation
         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)
@@ -125,24 +127,24 @@ c---angle selection
       do iang=1,nang
        idstep=0
        do i=2,nstep-1
-        jstmp2(1,i)=jsinref(iang,i+1)-jsinref(iang,i-1)
+        jsout(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)
+           jsout(2,idstep)=jsinref(iang,i)
+           jsout(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.
+         jsout(4,i)=minval(minloc(
+     1     abs(jsinref(iang,1:int(jsout(3,i)-1))-jsout(2,i*2))))
+         dangtmp(2,iang,int(jsout(4,i)):int(jsout(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.
+        jsout(5,i)=dang2(iang,i+1)-dang2(iang,i-1)
+        if(abs(jsout(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.
@@ -209,7 +211,12 @@ c---angle evaluation: in-out
        if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif
        call convim(jsin,dangw,dangi2,angwin(:)
      1       ,nnin,nwout,bd1,bd2)
-
+!        write(6,1115) real(dangi(1:60))
+!        write(6,*) 'dangi2'
+!        write(6,1115) real(dangi2)
+!        write(6,*) 'angwin'
+!        write(6,1115) real(angwin)
+       
        dangi(:)=datain(12,:) !--output-ref
        dangi2(1)=dangi(1)
        do i=2,nnin
@@ -221,41 +228,47 @@ c---angle evaluation: in-out
        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
+!       write(6, *) 'angwout', nnin, nwout
+!       write(6,1115) real(angwout(1:60))
+1115   format(15f7.1)
+       dangw=angwin-angwout !--difference    
        call angrev(dangw,dangw,nwout)
+!       write(6,1115) real(dangw)
+       
        idstep=0; idstep2=0
        do i=2,nwout
         angw(i)=dangw(i)-dangw(i-1)
-        if(angw(i).ge.300.) then
+        if(angw(i).ge.300.) then          
           idstep=idstep+1
-          jstmp2(2,idstep)=jswout(i)
-          jstmp2(3,idstep)=float(i)
+          jsout(2,idstep)=jswout(i)
+          jsout(3,idstep)=float(i)  
         endif
         if(angw(i).le.-300.) then
           idstep2=idstep2+1
-          jstmp2(4,idstep2)=jswout(i)
-          jstmp2(5,idstep2)=float(i)
+          jsout(4,idstep2)=jswout(i)
+          jsout(5,idstep2)=float(i)          
         endif
-       enddo
-
+       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.
+        dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))
+     1    =dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))+360.
+!!!!!! elena : dimensions mismatch  
+        dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2)))
+     1    =dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2)))-360.
+        ddwout1(11,int(jsout(5,idstep*2-1)):int(jsout(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.
+       
+       if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then           
+        dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))
+     1    =dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))-360.
+!!!!!! elena : dimensions mismatch
+        dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))
+     1    =dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))+360.
        endif
+       
        ddwout1(9,:)=dangw !separation angle
-
+!       write(6,1115) real(dangw) 
 c---output
       open(ifn2,file=fnout,status='unknown',form='formatted')
       do i=1,nwout,1
@@ -275,7 +288,7 @@ c       write(62,&#39;(i12,12e18.9)&#39;) int(jswout(i)),(ddwout1(j,i),j=1,11)
 c---finish
       deallocate(datain,jsin,dangi,dangi2
      1  ,jsout,ddout,jsinref,ang0,dango,dang1,jsoutrev
-     1  ,jswout,ddwout,ddwout1,jstmp2)
+     1  ,jswout,ddwout,ddwout1)
       return
       end
 c==================================================
@@ -0,0 +1,282 @@
+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==================================================
+
@@ -0,0 +1,294 @@
+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(:,:,:)
+!!!! elena : para(13) !!! not para(12)
+      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))
+
+      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
+        jsout(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
+           jsout(2,idstep)=jsinref(iang,i)
+           jsout(3,idstep)=float(i)
+        endif
+       enddo
+       if(idstep.ge.1)then
+        do i=1,idstep/2.
+         jsout(4,i)=minval(minloc(
+     1     abs(jsinref(iang,1:int(jsout(3,i)-1))-jsout(2,i*2))))
+         dangtmp(2,iang,int(jsout(4,i)):int(jsout(3,i*2)))=901.
+        enddo
+       endif
+       do i=2,nstep-1
+        jsout(5,i)=dang2(iang,i+1)-dang2(iang,i-1)
+        if(abs(jsout(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)
+!        write(6,1115) real(dangi(1:60))
+!        write(6,*) 'dangi2'
+!        write(6,1115) real(dangi2)
+!        write(6,*) 'angwin'
+!        write(6,1115) real(angwin)
+       
+       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)
+!       write(6, *) 'angwout', nnin, nwout
+!       write(6,1115) real(angwout(1:60))
+1115   format(15f7.1)
+       dangw=angwin-angwout !--difference    
+       call angrev(dangw,dangw,nwout)
+!       write(6,1115) real(dangw)
+       
+       idstep=0; idstep2=0
+       do i=2,nwout
+        angw(i)=dangw(i)-dangw(i-1)
+        if(angw(i).ge.300.) then          
+          idstep=idstep+1
+          jsout(2,idstep)=jswout(i)
+          jsout(3,idstep)=float(i)  
+        endif
+        if(angw(i).le.-300.) then
+          idstep2=idstep2+1
+          jsout(4,idstep2)=jswout(i)
+          jsout(5,idstep2)=float(i)          
+        endif
+       enddo             
+       
+       if(idprop.eq.1.and.idstep.ge.1.and.idstep*2.eq.idstep2)then
+        dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))
+     1    =dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))+360.
+!!!!!! elena : dimensions mismatch  
+        dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2)))
+     1    =dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2)))-360.
+        ddwout1(11,int(jsout(5,idstep*2-1)):int(jsout(5,idstep*2)))=1.
+       endif
+       
+       if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then           
+        dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))
+     1    =dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))-360.
+!!!!!! elena : dimensions mismatch
+        dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))
+     1    =dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))+360.
+       endif
+       
+       ddwout1(9,:)=dangw !separation angle
+!       write(6,1115) real(dangw) 
+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)
+      return
+      end
+c==================================================
+