Commit 5c44cc6c923470fbb1f5ac9b128328133f25a774
1 parent
6ed9c6e8
Exists in
master
bug
Showing
3 changed files
with
70 additions
and
50 deletions
Show diff stats
Sources/Prop_code/main_sw1da_p.f
| ... | ... | @@ -8,32 +8,33 @@ c======================================================================| |
| 8 | 8 | real*8 :: angref0,touts,dtr |
| 9 | 9 | real*8 :: angref(30) |
| 10 | 10 | integer :: instop,instopfin,ifnpo,idprop |
| 11 | - character*100 :: fnin,fnout,fnin1,fnout1,fntmp(10),fntmp2(10) | |
| 11 | + character*100 :: fnin,fnout,fnin1,fnout1,fntmp(12),fntmp2(12) | |
| 12 | 12 | character*100 :: fdirtmp |
| 13 | 13 | real*8,allocatable :: angref6(:) |
| 14 | 14 | real*8 :: xmin1,xmax1,xear |
| 15 | 15 | integer :: ix |
| 16 | - | |
| 17 | -c---setting | |
| 16 | +!!!!! elena : namelist should be before executable statements | |
| 18 | 17 | namelist /INPARA1/instop,fnin,fnout,ifnpo,dtr,touts,idprop |
| 19 | 18 | 1 ,fdirtmp,angref,idp_in,idp_prop,idp_out,xmin1,xmax1 |
| 19 | + | |
| 20 | +c---setting | |
| 21 | + angref(:)=-1.e10 | |
| 22 | + | |
| 20 | 23 | open(unit=50,file='namelist',status='old',form='formatted') |
| 21 | 24 | read(50,nml=INPARA1) |
| 22 | 25 | close(50) |
| 23 | - | |
| 24 | - angref(:)=-1.e10 | |
| 26 | + | |
| 25 | 27 | c----set angle |
| 26 | 28 | do i=1,30 |
| 27 | 29 | if(angref(i).eq.-1.e10) goto 319 |
| 28 | 30 | enddo |
| 29 | 31 | 319 continue |
| 30 | 32 | nang=i-1 |
| 31 | - allocate(angref6(nang)) | |
| 32 | 33 | |
| 33 | - ! angref6(1:nang)=angref | |
| 34 | + allocate(angref6(nang)) | |
| 34 | 35 | !!!! elena : dimensions mismatch 6/30 |
| 35 | 36 | angref6(1:nang)=angref(1:nang) |
| 36 | - | |
| 37 | + | |
| 37 | 38 | c----set filename |
| 38 | 39 | do iang=1,nang |
| 39 | 40 | angref0=angref6(iang) |
| ... | ... | @@ -52,7 +53,7 @@ c---prepare dataset at angref0 |
| 52 | 53 | if(instop.eq.0) call count_lines(fntmp(1),instop) |
| 53 | 54 | |
| 54 | 55 | c---solar wind propagation along angref0 |
| 55 | - do iang=1,nang | |
| 56 | + do iang=1,nang | |
| 56 | 57 | angref0=angref6(iang) |
| 57 | 58 | fnin1=fntmp(iang) |
| 58 | 59 | fnout1=fntmp2(iang) |
| ... | ... | @@ -83,12 +84,12 @@ c======================================================================| |
| 83 | 84 | character*100 :: fnin,fnout |
| 84 | 85 | !!!!! elena |
| 85 | 86 | real*8 jjs |
| 86 | - | |
| 87 | - namelist /INPARA2/npout,bx0at1au,gm | |
| 87 | + | |
| 88 | + namelist /INPARA2/npout,bx0at1au,gm | |
| 88 | 89 | open(unit=50,file='namelist',status='old',form='formatted') |
| 89 | 90 | read(50,nml=INPARA2) |
| 90 | 91 | close(50) |
| 91 | - | |
| 92 | + | |
| 92 | 93 | touts=touts0*0.5 |
| 93 | 94 | ix=int((xmax1-xmin1)*400.+1.e-10) |
| 94 | 95 | xear=(1.-xmin1)*400. |
| ... | ... | @@ -177,16 +178,15 @@ c obtain time spacing |
| 177 | 178 | else |
| 178 | 179 | write(fnpo,*) 'Time Spacing is NOT Correct! : dt0 > dt' |
| 179 | 180 | write(fnpo,*) 'dt=',dt |
| 180 | - write(fnpo,*) 'dt0=',dt0 | |
| 181 | + write(fnpo,*) 'dt0=',dt0 | |
| 181 | 182 | endif |
| 182 | 183 | if (merr.ne.0) goto 9999 |
| 183 | 184 | |
| 184 | 185 | |
| 185 | 186 | c----data read & input---------- |
| 187 | + read(fnum,*) jjs,(parain(i),i=1,13) | |
| 186 | 188 | !!!!! elena : js is defined as REAL in the file |
| 187 | - read(fnum,*) jjs,(parain(i),i=1,13) | |
| 188 | 189 | js=int(jjs) |
| 189 | -! read(fnum,*) js,(parain(i),i=1,13) | |
| 190 | 190 | parain(3)=parain(3)*float(idprop) !--REV0318 |
| 191 | 191 | xin=minval(minloc(abs(parain(9)-x))) |
| 192 | 192 | if(parain(9).lt.x(xin).and.idprop.eq.1) xin=xin-1 !--REV0318 |
| ... | ... | @@ -196,7 +196,7 @@ c----data read & input---------- |
| 196 | 196 | parainref(1:2,ns)=parain(9:10) |
| 197 | 197 | parainref(3,ns)=x(xin) |
| 198 | 198 | parainref(4,ns)=parain(3) !km/s |
| 199 | - | |
| 199 | + | |
| 200 | 200 | call datin_varied3(ro,pr,vx,vy,vz,bx,by,bz |
| 201 | 201 | & ,vstar,ns,ix,parain,xin,paranm ,idprop,x) !--REV0318 |
| 202 | 202 | |
| ... | ... | @@ -218,20 +218,26 @@ c solve hydrodynamic equations |
| 218 | 218 | & ,gx,dsc,scm,dv,rr,rrm,drr,dx,ix) |
| 219 | 219 | call bnd3_2(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix,idprop,x) !--REV0318 |
| 220 | 220 | c----------------------------------------------------------------------| |
| 221 | -c data output | |
| 222 | - | |
| 221 | +c data output | |
| 222 | + | |
| 223 | 223 | c---monitoring |
| 224 | - if(minval(pr).lt.0.) write(6,*) ns,minval(pr) | |
| 225 | - if(minval(pr).lt.0.) stop 'Pr<0' | |
| 224 | + if(minval(pr).lt.0.) write(6,*) ns | |
| 225 | + if(minval(pr).lt.0.) write(6,*) minval(pr) | |
| 226 | + if(minval(pr).lt.0.) stop 'Pr<0' | |
| 226 | 227 | if(idprop.eq.1.and.minval(vx).lt.0.) stop 'Vx opposite' |
| 227 | - if(idprop.eq.-1.and.maxval(vx).gt.0.) stop 'Vx opposite' | |
| 228 | - | |
| 228 | + if(idprop.eq.-1.and.maxval(vx).gt.0.) then | |
| 229 | + write(6,1115) vx | |
| 230 | + stop 'Vx opposite' | |
| 231 | + endif | |
| 232 | +1115 format(15f10.7) | |
| 229 | 233 | c----first dt estimation |
| 230 | 234 | if (mod(int(ns),int(touts)).eq.0) then |
| 231 | 235 | js=js0+ns*dtr *float(idprop) !--REV0318 |
| 232 | 236 | xout=minval(minloc(abs(parain(11)-x))) |
| 233 | 237 | dtp1=(parain(11)-parain(9))*1.5d11/400./1.d3 |
| 234 | - | |
| 238 | + | |
| 239 | +! write(6,*) 'first dt estimation', dtp1 | |
| 240 | + | |
| 235 | 241 | if((js-dtp1.le.js0.and.idprop.eq.1) !--REV0318 |
| 236 | 242 | 1 .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then !--REV0318 |
| 237 | 243 | write(fnum2,704) js,1.,1.,100.,1.,1.,1.d-5,1.d-5,1.d-5 |
| ... | ... | @@ -239,44 +245,45 @@ c----first dt estimation |
| 239 | 245 | else |
| 240 | 246 | itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref))) !--REV0318 |
| 241 | 247 | |
| 242 | -!!!!! elena | |
| 248 | + !!!!! elena | |
| 243 | 249 | if(parainref(4,itin).eq.0.0)then |
| 244 | 250 | write(6,*) '!!!second dt',itin, ns, parain(12), parain(10) |
| 245 | 251 | itin = ns |
| 246 | 252 | endif |
| 247 | 253 | |
| 248 | 254 | c----second dt estimation |
| 249 | - dtp1=(parain(11)-parain(9))*1.5d11/parainref(4,itin)/1.d3 | |
| 250 | - | |
| 255 | + dtp1=(parain(11)-parain(9))*1.5d11/parainref(4,itin)/1.d3 | |
| 256 | + | |
| 251 | 257 | if((js-dtp1.le.js0.and.idprop.eq.1) |
| 252 | 258 | 1 .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then |
| 253 | 259 | write(fnum2,704) js,1.,1.,100.,1.,1.,1.d-5,1.d-5,1.d-5 |
| 254 | 260 | & ,0.d0,float(js0) |
| 255 | 261 | goto 1001 |
| 256 | 262 | endif |
| 263 | + | |
| 257 | 264 | itin_tmp=itin |
| 258 | 265 | itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref))) |
| 259 | 266 | |
| 260 | 267 | !!!!! elena |
| 261 | 268 | if (parainref(4,itin).eq.0.0) itin=itin_tmp |
| 262 | - | |
| 263 | 269 | c----output |
| 264 | 270 | dang0=-(parain(12)-parainref(2,itin)) |
| 265 | 271 | dang0=mod(dang0,360.) |
| 266 | 272 | if(dang0.ge.180.) dang0=dang0-360. |
| 273 | + | |
| 267 | 274 | if(angref0.ge.0.) |
| 268 | 275 | & dtp2=(parain(11)-x(xout))*1.5d11/vx(xout)/400.d3 |
| 269 | 276 | & -(parainref(1,itin)-parainref(3,itin))*1.5d11 !x_datain-x_in |
| 270 | - & /parainref(4,itin)/1.d3 | |
| 271 | - !!!! elena | |
| 277 | + & /parainref(4,itin)/1.d3 | |
| 278 | + !!!! elena | |
| 272 | 279 | if (parainref(4,itin).eq.0.0) then ! NaN |
| 273 | 280 | write(6,*) '!!!!!! NaN in dtp2 ', angref0 |
| 274 | 281 | write(6,*) parain(11), x(xout),vx(xout) |
| 275 | 282 | write(6,*) ns, itin, xout |
| 276 | 283 | write(6,*) parain(12), parain(10) |
| 277 | 284 | stop |
| 278 | - endif | |
| 279 | - | |
| 285 | + endif | |
| 286 | + | |
| 280 | 287 | write(fnum2,704) js,ro(xout)*paranm(1) |
| 281 | 288 | & ,pr(xout)*1937./ro(xout)*paranm(2) |
| 282 | 289 | & ,vx(xout)*paranm(3)*float(idprop) !--REV0318 |
| ... | ... | @@ -285,6 +292,7 @@ c----output |
| 285 | 292 | & ,bz(xout)*paranm(8)*16.34 |
| 286 | 293 | & ,dtp2*float(idprop),jsinref(itin)-parainref(6,itin) |
| 287 | 294 | endif |
| 295 | + | |
| 288 | 296 | 704 format(i15.1, 9e18.8,f18.1) |
| 289 | 297 | |
| 290 | 298 | 1001 continue |
| ... | ... | @@ -670,7 +678,7 @@ c----------------------------------------------------------------------| |
| 670 | 678 | pi = acos(-1.0d0) |
| 671 | 679 | pi4i=0.25/pi |
| 672 | 680 | dtmin=2.0e-10 |
| 673 | -! safety=0.8 !cf. 0.4 --original | |
| 681 | +! safety=0.8 !cf. 0.4 --original | |
| 674 | 682 | safety=0.4 |
| 675 | 683 | c----------------------------------------------------------------------| |
| 676 | 684 | dt=1.e20 |
| ... | ... | @@ -689,6 +697,7 @@ c----------------------------------------------------------------------| |
| 689 | 697 | dt=dtq(i) |
| 690 | 698 | endif |
| 691 | 699 | enddo |
| 700 | +! write(6,*) minval(pr) | |
| 692 | 701 | c----------------------------------------------------------------------| |
| 693 | 702 | c write the point where dt is smaller than critical value |
| 694 | 703 | c----------------------------------------------------------------------| | ... | ... |
Sources/Prop_code/sub_indat_p.f
| ... | ... | @@ -29,9 +29,8 @@ c======================================================================| |
| 29 | 29 | open(ifn1,file=fnin,status='old') |
| 30 | 30 | do i=1,nnin |
| 31 | 31 | read(ifn1,*) js1,(parain(j),j=1,12) |
| 32 | - !!!!!!! elena : Dimensions mismatch (13/12) | |
| 33 | - datain(1:12,i)=parain(:) | |
| 34 | -! datain(:,i)=parain(:) | |
| 32 | +!!!!!!! elena : Dimensions mismatch (13/12) | |
| 33 | + datain(1:12,i)=parain(:) | |
| 35 | 34 | jsin(i)=js1 |
| 36 | 35 | enddo |
| 37 | 36 | dtin=(jsin(2)-jsin(1))*float(idprop) !--REV0318 | ... | ... |
Sources/Prop_code/sub_slct_p.f
| ... | ... | @@ -7,7 +7,7 @@ c======================================================================| |
| 7 | 7 | real*8,intent(in) :: angref6(nang),dtr,touts |
| 8 | 8 | character*100 ,intent(in):: fntmp2(nang),fnin,fnout |
| 9 | 9 | real*8,allocatable :: jsout(:,:),ddout(:,:,:) |
| 10 | -!!!! elena : para(13) !!! not para(12) | |
| 10 | +!!!! elena : para(13) !!! not para(12) | |
| 11 | 11 | real*8 :: angref0,para(13),js1,sec0 |
| 12 | 12 | integer :: i,j,ifn0,ifn1,ifn2,iang,nstep |
| 13 | 13 | integer :: iy0,im0,id0,jdd,jdd2,ihr0,imin0,isec0 |
| ... | ... | @@ -84,7 +84,8 @@ c---angle evaluation |
| 84 | 84 | bd1=dangi2(1); bd2=dangi2(nnin) |
| 85 | 85 | call convim(jsin,jsinref(iang,:),dangi2,dang1(iang,:) |
| 86 | 86 | 1 ,nnin,nstep,bd1,bd2) |
| 87 | - | |
| 87 | + | |
| 88 | + | |
| 88 | 89 | dangi(:)=datain(12,:) !--input-ref |
| 89 | 90 | dangi2(1)=dangi(1) |
| 90 | 91 | do i=2,nnin |
| ... | ... | @@ -95,6 +96,7 @@ c---angle evaluation |
| 95 | 96 | bd1=dangi2(1); bd2=dangi2(nnin) |
| 96 | 97 | call convim(jsin,jsoutrev(iang,:),dangi2,ang0(:) |
| 97 | 98 | 1 ,nnin,nstep,bd1,bd2) |
| 99 | + | |
| 98 | 100 | dango(:)=ang0(:)-angref0 !--output-ref |
| 99 | 101 | |
| 100 | 102 | c---planeatry rotation effect (2nd step) |
| ... | ... | @@ -208,7 +210,12 @@ c---angle evaluation: in-out |
| 208 | 210 | if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif |
| 209 | 211 | call convim(jsin,dangw,dangi2,angwin(:) |
| 210 | 212 | 1 ,nnin,nwout,bd1,bd2) |
| 211 | - | |
| 213 | +! write(6,1115) real(dangi(1:60)) | |
| 214 | +! write(6,*) 'dangi2' | |
| 215 | +! write(6,1115) real(dangi2) | |
| 216 | +! write(6,*) 'angwin' | |
| 217 | +! write(6,1115) real(angwin) | |
| 218 | + | |
| 212 | 219 | dangi(:)=datain(12,:) !--output-ref |
| 213 | 220 | dangi2(1)=dangi(1) |
| 214 | 221 | do i=2,nnin |
| ... | ... | @@ -220,41 +227,46 @@ c---angle evaluation: in-out |
| 220 | 227 | if(idprop.eq.-1)then;bd1=dangi2(nnin);bd2=dangi2(1);endif |
| 221 | 228 | call convim(jsin,jswout,dangi2,angwout(:) |
| 222 | 229 | 1 ,nnin,nwout,bd1,bd2) |
| 223 | - | |
| 224 | - dangw=angwin-angwout !--difference | |
| 230 | +! write(6, *) 'angwout', nnin, nwout | |
| 231 | +! write(6,1115) real(angwout(1:60)) | |
| 232 | +1115 format(15f7.1) | |
| 233 | + dangw=angwin-angwout !--difference | |
| 225 | 234 | call angrev(dangw,dangw,nwout) |
| 235 | +! write(6,1115) real(dangw) | |
| 236 | + | |
| 226 | 237 | idstep=0; idstep2=0 |
| 227 | 238 | do i=2,nwout |
| 228 | 239 | angw(i)=dangw(i)-dangw(i-1) |
| 229 | - if(angw(i).ge.300.) then | |
| 240 | + if(angw(i).ge.300.) then | |
| 230 | 241 | idstep=idstep+1 |
| 231 | 242 | jsout(2,idstep)=jswout(i) |
| 232 | - jsout(3,idstep)=float(i) | |
| 243 | + jsout(3,idstep)=float(i) | |
| 233 | 244 | endif |
| 234 | 245 | if(angw(i).le.-300.) then |
| 235 | 246 | idstep2=idstep2+1 |
| 236 | 247 | jsout(4,idstep2)=jswout(i) |
| 237 | - jsout(5,idstep2)=float(i) | |
| 248 | + jsout(5,idstep2)=float(i) | |
| 238 | 249 | endif |
| 239 | - enddo | |
| 240 | - | |
| 250 | + enddo | |
| 251 | + | |
| 241 | 252 | if(idprop.eq.1.and.idstep.ge.1.and.idstep*2.eq.idstep2)then |
| 242 | 253 | dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep))) |
| 243 | 254 | 1 =dangw(int(jsout(5,idstep*2-1)):int(jsout(3,idstep)))+360. |
| 244 | - !!! elena dimensions mismatch | |
| 245 | - dangw(int(jsout(3,idstep))-1:int(jsout(5,idstep*2))) | |
| 255 | + dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2))) | |
| 246 | 256 | 1 =dangw(int(jsout(3,idstep))-1:int(jsout(5,idstep*2)))-360. |
| 247 | 257 | ddwout1(11,int(jsout(5,idstep*2-1)):int(jsout(5,idstep*2)))=1. |
| 248 | 258 | endif |
| 249 | - if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then | |
| 259 | + | |
| 260 | + if(idprop.eq.-1.and.idstep2.ge.1.and.idstep2*2.eq.idstep)then | |
| 250 | 261 | dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2))) |
| 251 | 262 | 1 =dangw(int(jsout(3,idstep2*2-1)):int(jsout(5,idstep2)))-360. |
| 252 | - !!! elena dimensions mismatch | |
| 263 | +!!!!!! elena : dimensions mismatch | |
| 253 | 264 | dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2))) |
| 254 | 265 | 1 =dangw(int(jsout(5,idstep2))-1:int(jsout(3,idstep2*2)))+360. |
| 255 | 266 | endif |
| 267 | + | |
| 256 | 268 | ddwout1(9,:)=dangw !separation angle |
| 257 | - | |
| 269 | +! write(6,1115) real(dangw) | |
| 258 | 270 | c---output |
| 259 | 271 | open(ifn2,file=fnout,status='unknown',form='formatted') |
| 260 | 272 | do i=1,nwout,1 | ... | ... |