Commit 7837c6bb9f8a375624f50165e25268d578b41428
1 parent
a50b7398
Exists in
master
second modif
Showing
2 changed files
with
301 additions
and
2 deletions
Show diff stats
Sources/Prop_code/sub_slct_p.f
| ... | ... | @@ -161,9 +161,11 @@ c---angle selection |
| 161 | 161 | |
| 162 | 162 | c---set time |
| 163 | 163 | dtout=touts*dtr |
| 164 | - call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(:,1))) | |
| 164 | +!!!! elena : jsoutrev(1:nang,1) | |
| 165 | + call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(1:nang,1))) | |
| 165 | 166 | call ymds2js(iy0,im0,id0+int((idprop+1)*0.5),0.d0,jswout0) |
| 166 | - call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(:,nstep))) | |
| 167 | +!!!! elena : jsoutrev(1:nang,nstep) | |
| 168 | + call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(1:nang,nstep))) | |
| 167 | 169 | call ymds2js(iy0,im0,id0+int((idprop-1)*0.5),0.d0,jswout1) |
| 168 | 170 | |
| 169 | 171 | nwout=abs(jswout1-jswout0)/dtout | ... | ... |
| ... | ... | @@ -0,0 +1,297 @@ |
| 1 | +c======================================================================| | |
| 2 | + subroutine sub_slct(fnin,fnout,fntmp2,angref6,nang | |
| 3 | + 1 ,touts,dtr,idprop,ifnpo) | |
| 4 | +c======================================================================| | |
| 5 | + implicit none | |
| 6 | + integer,intent(in) :: nang,ifnpo | |
| 7 | + real*8,intent(in) :: angref6(nang),dtr,touts | |
| 8 | + character*100 ,intent(in):: fntmp2(nang),fnin,fnout | |
| 9 | + real*8,allocatable :: jsout(:,:),ddout(:,:,:) | |
| 10 | +!!!! elena : para(13) !!! not para(12) | |
| 11 | + real*8 :: angref0,para(13),js1,sec0 | |
| 12 | + integer :: i,j,ifn0,ifn1,ifn2,iang,nstep | |
| 13 | + integer :: iy0,im0,id0,jdd,jdd2,ihr0,imin0,isec0 | |
| 14 | + integer :: nwout,nnin | |
| 15 | + real*8 :: dtout,jswout0,jswout1,dtin | |
| 16 | + real*8 :: pdy0 | |
| 17 | + real*8,allocatable :: jswout(:),ddwout(:,:,:),ddwout1(:,:) | |
| 18 | + real*8,allocatable :: datain(:,:),jsin(:),jsinref(:,:) | |
| 19 | + real*8,allocatable :: dang2(:,:),dang1(:,:) | |
| 20 | + real*8,allocatable :: dangi(:),dango(:),ang0(:) | |
| 21 | + real*8,allocatable :: jsoutrev(:,:),ddoutrev(:,:,:) | |
| 22 | + integer :: idprop | |
| 23 | + real*8,allocatable :: angwin(:),angwout(:),dangw(:),angw(:) | |
| 24 | + real*8,allocatable :: jstmp(:) | |
| 25 | + real*8,allocatable :: dangi2(:) | |
| 26 | + real*8 :: jstmp0 ,dangtmp0,bd1,bd2 | |
| 27 | + integer :: idstep,idstep2 | |
| 28 | + integer,allocatable :: islct(:) | |
| 29 | + real*8,allocatable :: dangtmp(:,:,:),ddslct(:,:),jsslct(:) | |
| 30 | + real*8,allocatable :: jsinref2(:) | |
| 31 | + | |
| 32 | + ifn1=43 | |
| 33 | + ifn2=44 | |
| 34 | +!---read input data | |
| 35 | + call count_lines(fnin,nnin) | |
| 36 | + allocate(datain(13,nnin),jsin(nnin),dangi(nnin),dangi2(nnin)) | |
| 37 | + open(ifn2,file=fnin,status='old') | |
| 38 | + do i=1,nnin | |
| 39 | + read(ifn2,*) js1,(para(j),j=1,13) | |
| 40 | + datain(:,i)=para(1:13) | |
| 41 | + jsin(i)=js1 | |
| 42 | + enddo | |
| 43 | + dtin=jsin(2)-jsin(1) | |
| 44 | + close(ifn2,status='keep') | |
| 45 | + | |
| 46 | +!---count line & read simulation data | |
| 47 | + iang=1 | |
| 48 | + angref0=angref6(iang) | |
| 49 | + call count_lines(fntmp2(iang),nstep) | |
| 50 | +!!!! elena : jsout(nang,nstep), jsoutrev(nang,nstep) => jsout(6,nstep), jsoutrev(6,nstep) | |
| 51 | + allocate(jsout(6,nstep),ddout(nang,10,nstep) | |
| 52 | + 1 ,jsinref(nang,nstep),ang0(nstep),dango(nstep) | |
| 53 | + 1 ,dang1(nang,nstep),dang2(nang,nstep),jsoutrev(6,nstep)) | |
| 54 | + allocate(jstmp(nstep)) | |
| 55 | + allocate(ddslct(10,nstep),jsslct(nstep)) | |
| 56 | + allocate(islct(nstep),dangtmp(2,nang,nstep),jsinref2(nstep)) | |
| 57 | + | |
| 58 | + do iang=1,nang | |
| 59 | + angref0=angref6(iang) | |
| 60 | + open(ifn1,file=fntmp2(iang),status='old') | |
| 61 | + do i=1,nstep | |
| 62 | + read(ifn1,*) js1,(para(j),j=1,10) | |
| 63 | + jsout(iang,i)=js1 | |
| 64 | + ddout(iang,:,i)=para(1:10) | |
| 65 | + enddo | |
| 66 | + close(ifn1,status='keep') | |
| 67 | + enddo | |
| 68 | + jsinref(:,:)=ddout(:,10,:) | |
| 69 | + | |
| 70 | + write(ifnpo,*) nnin,nstep | |
| 71 | +c---time shift for radial difference | |
| 72 | + jsoutrev(:,:)=jsout(:,:) | |
| 73 | + | |
| 74 | +c---angle evaluation | |
| 75 | + do iang=1,nang !---6 ref.angles | |
| 76 | + angref0=angref6(iang) | |
| 77 | + | |
| 78 | + dangi(:)=datain(10,:)-angref0 !--input-ref | |
| 79 | + dangi2(1)=dangi(1) | |
| 80 | + do i=2,nnin | |
| 81 | + dangtmp0=dangi(i)-dangi(i-1) | |
| 82 | + if(abs(dangtmp0).ge.300.) dangtmp0=0. | |
| 83 | + dangi2(i)=dangi2(i-1)+dangtmp0 | |
| 84 | + enddo | |
| 85 | + bd1=dangi2(1); bd2=dangi2(nnin) | |
| 86 | + call convim(jsin,jsinref(iang,:),dangi2,dang1(iang,:) | |
| 87 | + 1 ,nnin,nstep,bd1,bd2) | |
| 88 | + | |
| 89 | + | |
| 90 | + dangi(:)=datain(12,:) !--input-ref | |
| 91 | + dangi2(1)=dangi(1) | |
| 92 | + do i=2,nnin | |
| 93 | + dangtmp0=dangi(i)-dangi(i-1) | |
| 94 | + if(abs(dangtmp0).ge.300.) dangtmp0=0. | |
| 95 | + dangi2(i)=dangi2(i-1)+dangtmp0 | |
| 96 | + enddo | |
| 97 | + bd1=dangi2(1); bd2=dangi2(nnin) | |
| 98 | + call convim(jsin,jsoutrev(iang,:),dangi2,ang0(:) | |
| 99 | + 1 ,nnin,nstep,bd1,bd2) | |
| 100 | + | |
| 101 | + dango(:)=ang0(:)-angref0 !--output-ref | |
| 102 | + | |
| 103 | +c---planeatry rotation effect (2nd step) | |
| 104 | + do idstep=1,3 | |
| 105 | + call angrev(dango,dango,nstep) | |
| 106 | + jstmp=jsoutrev(iang,:)+dango/360.*26.*86400 | |
| 107 | + do i=2,nstep | |
| 108 | + if(abs(jstmp(i)-jstmp(i-1)).gt.(maxval(jstmp)-minval(jstmp)) | |
| 109 | + 1 /float(nstep)*5.) jstmp(i)=jstmp(i-1) | |
| 110 | + enddo | |
| 111 | + | |
| 112 | + bd1=dangi2(1); bd2=dangi2(nnin) | |
| 113 | + call convim(jsin,jstmp,dangi2,ang0(:) | |
| 114 | + 1 ,nnin,nstep,bd1,bd2) | |
| 115 | + dango(:)=ang0(:)-angref0 !--output-ref | |
| 116 | + enddo | |
| 117 | + dang2(iang,:)=dango(:) | |
| 118 | + enddo | |
| 119 | +c---angle selection | |
| 120 | + do iang=1,nang | |
| 121 | + call angrev(dang1(iang,:),dang1(iang,:),nstep) | |
| 122 | + call angrev(dang2(iang,:),dang2(iang,:),nstep) | |
| 123 | + enddo | |
| 124 | + dangtmp(2,:,:)=abs(dang1(:,:))+abs(dang2(:,:)) | |
| 125 | + dangtmp(1,:,:)=dang1(:,:)+dang2(:,:) | |
| 126 | + | |
| 127 | + do iang=1,nang | |
| 128 | + idstep=0 | |
| 129 | + do i=2,nstep-1 | |
| 130 | + jsout(1,i)=jsinref(iang,i+1)-jsinref(iang,i-1) | |
| 131 | + if((idprop.eq.-1.and.jsout(1,i).ge.86400.*20.).or. | |
| 132 | + 1 (idprop.eq.1.and.jsout(1,i).le.-86400.*20.)) then | |
| 133 | + idstep=idstep+1 | |
| 134 | + jsout(2,idstep)=jsinref(iang,i) | |
| 135 | + jsout(3,idstep)=float(i) | |
| 136 | + endif | |
| 137 | + enddo | |
| 138 | + if(idstep.ge.1)then | |
| 139 | + do i=1,idstep/2. | |
| 140 | + jsout(4,i)=minval(minloc( | |
| 141 | + 1 abs(jsinref(iang,1:int(jsout(3,i)-1))-jsout(2,i*2)))) | |
| 142 | + dangtmp(2,iang,int(jsout(4,i)):int(jsout(3,i*2)))=901. | |
| 143 | + enddo | |
| 144 | + endif | |
| 145 | + do i=2,nstep-1 | |
| 146 | + jsout(5,i)=dang2(iang,i+1)-dang2(iang,i-1) | |
| 147 | + if(abs(jsout(5,i)).lt.1.e-10) dangtmp(2,iang,i)=902. | |
| 148 | + if((idprop.eq.-1.and.ddout(iang,10,i).ge.jsin(1)) | |
| 149 | + 1 .or.(idprop.eq.1.and.ddout(iang,10,i).le.jsin(1))) | |
| 150 | + 1 dangtmp(2,iang,i)=999. | |
| 151 | + enddo | |
| 152 | + enddo | |
| 153 | + | |
| 154 | + do i=1,nstep | |
| 155 | + islct(i)=minval(minloc(dangtmp(2,:,i))) | |
| 156 | + dango(i)=dangtmp(1,islct(i),i) | |
| 157 | + jsslct(i)=jsoutrev(islct(i),i)+ddout(islct(i),9,i) | |
| 158 | + 1 +dang2(islct(i),i)/360.*26.*86400. | |
| 159 | + ddslct(:,i)=ddout(islct(i),:,i) | |
| 160 | + enddo | |
| 161 | + | |
| 162 | +c---set time | |
| 163 | + dtout=touts*dtr | |
| 164 | +!!!! elena : jsoutrev(1:nang,1) | |
| 165 | + call js2ymds(iy0,im0,id0,sec0,maxval(jsoutrev(1:nang,1))) | |
| 166 | + call ymds2js(iy0,im0,id0+int((idprop+1)*0.5),0.d0,jswout0) | |
| 167 | +!!!! elena : jsoutrev(1:nang,nstep) | |
| 168 | + call js2ymds(iy0,im0,id0,sec0,minval(jsoutrev(1:nang,nstep))) | |
| 169 | + call ymds2js(iy0,im0,id0+int((idprop-1)*0.5),0.d0,jswout1) | |
| 170 | + | |
| 171 | + nwout=abs(jswout1-jswout0)/dtout | |
| 172 | + write(ifnpo,*) nwout | |
| 173 | + | |
| 174 | + allocate(jswout(nwout),ddwout(nang,10,nwout) | |
| 175 | + 1 ,ddwout1(11,nwout)) | |
| 176 | + do i=1,nwout | |
| 177 | + if(idprop.eq.1) jswout(i)=jswout0+dble(i)*dtout !--REV0318 | |
| 178 | + if(idprop.eq.-1) jswout(i)=jswout1+dble(i)*abs(dtout) !--REV0318 | |
| 179 | + enddo | |
| 180 | + | |
| 181 | +c---data and angle conversion | |
| 182 | + do i=1,10 | |
| 183 | + bd1=ddslct(i,1); bd2=ddslct(i,nstep) | |
| 184 | + if(idprop.eq.-1)then | |
| 185 | + bd1=ddslct(i,nstep); bd2=ddslct(i,1) | |
| 186 | + endif | |
| 187 | + call convim(jsslct(:),jswout(:),ddslct(i,:),ddwout1(i,:) | |
| 188 | + 1 ,nstep,nwout,bd1,bd2) | |
| 189 | + enddo | |
| 190 | +c---data lack | |
| 191 | + call convim(jsin,ddwout1(10,:),datain(13,:),ddwout1(11,:) | |
| 192 | + 1 ,nnin,nwout,datain(13,1),datain(13,nnin)) | |
| 193 | + do i=2,nwout | |
| 194 | + if(abs(ddwout1(3,i)-ddwout1(3,i-1)).le.abs(ddwout1(3,i))*1.e-10) | |
| 195 | + 1 ddwout1(11,i)=1. | |
| 196 | + enddo | |
| 197 | +c---angle evaluation: in-out | |
| 198 | + allocate(angwin(nwout),angwout(nwout),dangw(nwout),angw(nwout)) | |
| 199 | + | |
| 200 | + dangi(:)=datain(10,:) !--input-ref | |
| 201 | + dangi2(1)=dangi(1) | |
| 202 | + do i=2,nnin | |
| 203 | + dangtmp0=dangi(i)-dangi(i-1) | |
| 204 | + if(abs(dangtmp0).ge.300.) dangtmp0=0. | |
| 205 | + dangi2(i)=dangi2(i-1)+dangtmp0 | |
| 206 | + enddo | |
| 207 | + dangw=ddwout1(10,:) | |
| 208 | + do i=101,nwout-100 | |
| 209 | + dangw(i)=sum(ddwout1(10,i-100:i+100))/201. | |
| 210 | + enddo | |
| 211 | + | |
| 212 | + bd1=dangi2(1); bd2=dangi2(nnin) | |
| 213 | + if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif | |
| 214 | + call convim(jsin,dangw,dangi2,angwin(:) | |
| 215 | + 1 ,nnin,nwout,bd1,bd2) | |
| 216 | +! write(6,1115) real(dangi(1:60)) | |
| 217 | +! write(6,*) 'dangi2' | |
| 218 | +! write(6,1115) real(dangi2) | |
| 219 | +! write(6,*) 'angwin' | |
| 220 | +! write(6,1115) real(angwin) | |
| 221 | + | |
| 222 | + dangi(:)=datain(12,:) !--output-ref | |
| 223 | + dangi2(1)=dangi(1) | |
| 224 | + do i=2,nnin | |
| 225 | + dangtmp0=dangi(i)-dangi(i-1) | |
| 226 | + if(abs(dangtmp0).ge.300.) dangtmp0=0. | |
| 227 | + dangi2(i)=dangi2(i-1)+dangtmp0 | |
| 228 | + enddo | |
| 229 | + bd1=dangi2(1); bd2=dangi2(nnin) | |
| 230 | + if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif | |
| 231 | + call convim(jsin,jswout,dangi2,angwout(:) | |
| 232 | + 1 ,nnin,nwout,bd1,bd2) | |
| 233 | +! write(6, *) 'angwout', nnin, nwout | |
| 234 | +! write(6,1115) real(angwout(1:60)) | |
| 235 | +1115 format(15f7.1) | |
| 236 | + dangw=angwin-angwout !--difference | |
| 237 | + call angrev(dangw,dangw,nwout) | |
| 238 | +! write(6,1115) real(dangw) | |
| 239 | + | |
| 240 | + idstep=0; idstep2=0 | |
| 241 | + do i=2,nwout | |
| 242 | + angw(i)=dangw(i)-dangw(i-1) | |
| 243 | + if(angw(i).ge.300.) then | |
| 244 | + idstep=idstep+1 | |
| 245 | + jsout(2,idstep)=jswout(i) | |
| 246 | + jsout(3,idstep)=float(i) | |
| 247 | + endif | |
| 248 | + if(angw(i).le.-300.) then | |
| 249 | + idstep2=idstep2+1 | |
| 250 | + jsout(4,idstep2)=jswout(i) | |
| 251 | + jsout(5,idstep2)=float(i) | |
| 252 | + endif | |
| 253 | + enddo | |
| 254 | + | |
| 255 | + if(idprop.eq.1.and.idstep.ge.1.and.idstep*2.eq.idstep2)then | |
| 256 | + dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep))) | |
| 257 | + 1 =dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))+360. | |
| 258 | +!!!!!! elena : dimensions mismatch | |
| 259 | + dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2))) | |
| 260 | + 1 =dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2)))-360. | |
| 261 | + ddwout1(11,int(jsout(5,idstep*2-1)):int(jsout(5,idstep*2)))=1. | |
| 262 | + endif | |
| 263 | + | |
| 264 | + if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then | |
| 265 | + dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2))) | |
| 266 | + 1 =dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))-360. | |
| 267 | +!!!!!! elena : dimensions mismatch | |
| 268 | + dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2))) | |
| 269 | + 1 =dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))+360. | |
| 270 | + endif | |
| 271 | + | |
| 272 | + ddwout1(9,:)=dangw !separation angle | |
| 273 | +! write(6,1115) real(dangw) | |
| 274 | +c---output | |
| 275 | + open(ifn2,file=fnout,status='unknown',form='formatted') | |
| 276 | + do i=1,nwout,1 | |
| 277 | + call js2ymds(iy0,im0,id0,sec0,jswout(i)) | |
| 278 | + ihr0=int(sec0/3600.) | |
| 279 | + imin0=int((sec0-ihr0*3600.)/60.) | |
| 280 | + isec0=int(sec0-(ihr0*3600.+imin0*60.)) | |
| 281 | + pdy0=1.67e-27*ddwout1(1,i)*1.e6 | |
| 282 | + 1 *(ddwout1(3,i)**2+ddwout1(4,i)**2)*1.e6 | |
| 283 | + write(ifn2,'(i4.4,5(a1,i2.2),12(1pe13.5))') | |
| 284 | + 1 iy0,' ',im0,' ',id0,' ',ihr0,':',imin0,':',isec0 | |
| 285 | + 1 ,ddwout1(1:4,i),ddwout1(7,i) | |
| 286 | + 2 ,pdy0*1.e9,ddwout1(9,i),ddwout1(11,i) | |
| 287 | +c write(62,'(i12,12e18.9)') int(jswout(i)),(ddwout1(j,i),j=1,11) | |
| 288 | + enddo | |
| 289 | + | |
| 290 | +c---finish | |
| 291 | + deallocate(datain,jsin,dangi,dangi2 | |
| 292 | + 1 ,jsout,ddout,jsinref,ang0,dango,dang1,jsoutrev | |
| 293 | + 1 ,jswout,ddwout,ddwout1) | |
| 294 | + return | |
| 295 | + end | |
| 296 | +c================================================== | |
| 297 | + | ... | ... |