diff --git a/Sources/Prop_code/sub_slct_p.f b/Sources/Prop_code/sub_slct_p.f index e7e91d0..183a6e7 100644 --- a/Sources/Prop_code/sub_slct_p.f +++ b/Sources/Prop_code/sub_slct_p.f @@ -161,9 +161,11 @@ c---angle selection c---set time dtout=touts*dtr - call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(:,1))) +!!!! elena : jsoutrev(1:nang,1) + call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(1:nang,1))) call ymds2js(iy0,im0,id0+int((idprop+1)*0.5),0.d0,jswout0) - call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(:,nstep))) +!!!! elena : jsoutrev(1:nang,nstep) + call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(1:nang,nstep))) call ymds2js(iy0,im0,id0+int((idprop-1)*0.5),0.d0,jswout1) nwout=abs(jswout1-jswout0)/dtout diff --git a/Sources/Prop_code/sub_slct_p.f.my.new b/Sources/Prop_code/sub_slct_p.f.my.new new file mode 100644 index 0000000..183a6e7 --- /dev/null +++ b/Sources/Prop_code/sub_slct_p.f.my.new @@ -0,0 +1,297 @@ +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) +!!!! 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(6,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 +!!!! elena : jsoutrev(1:nang,1) + call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(1:nang,1))) + call ymds2js(iy0,im0,id0+int((idprop+1)*0.5),0.d0,jswout0) +!!!! elena : jsoutrev(1:nang,nstep) + call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(1:nang,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================================================== + -- libgit2 0.21.2