diff --git a/Sources/Prop_code/main_sw1da_p.f b/Sources/Prop_code/main_sw1da_p.f index 4a1cb3a..1b4700a 100644 --- a/Sources/Prop_code/main_sw1da_p.f +++ b/Sources/Prop_code/main_sw1da_p.f @@ -8,32 +8,33 @@ c======================================================================| real*8 :: angref0,touts,dtr real*8 :: angref(30) integer :: instop,instopfin,ifnpo,idprop - character*100 :: fnin,fnout,fnin1,fnout1,fntmp(10),fntmp2(10) + character*100 :: fnin,fnout,fnin1,fnout1,fntmp(12),fntmp2(12) character*100 :: fdirtmp real*8,allocatable :: angref6(:) real*8 :: xmin1,xmax1,xear integer :: ix - -c---setting +!!!!! elena : namelist should be before executable statements namelist /INPARA1/instop,fnin,fnout,ifnpo,dtr,touts,idprop 1 ,fdirtmp,angref,idp_in,idp_prop,idp_out,xmin1,xmax1 + +c---setting + angref(:)=-1.e10 + open(unit=50,file='namelist',status='old',form='formatted') read(50,nml=INPARA1) close(50) - - angref(:)=-1.e10 + c----set angle do i=1,30 if(angref(i).eq.-1.e10) goto 319 enddo 319 continue nang=i-1 - allocate(angref6(nang)) - ! angref6(1:nang)=angref + allocate(angref6(nang)) !!!! elena : dimensions mismatch 6/30 angref6(1:nang)=angref(1:nang) - + c----set filename do iang=1,nang angref0=angref6(iang) @@ -52,7 +53,7 @@ c---prepare dataset at angref0 if(instop.eq.0) call count_lines(fntmp(1),instop) c---solar wind propagation along angref0 - do iang=1,nang + do iang=1,nang angref0=angref6(iang) fnin1=fntmp(iang) fnout1=fntmp2(iang) @@ -83,12 +84,12 @@ c======================================================================| character*100 :: fnin,fnout !!!!! elena real*8 jjs - - namelist /INPARA2/npout,bx0at1au,gm + + namelist /INPARA2/npout,bx0at1au,gm open(unit=50,file='namelist',status='old',form='formatted') read(50,nml=INPARA2) close(50) - + touts=touts0*0.5 ix=int((xmax1-xmin1)*400.+1.e-10) xear=(1.-xmin1)*400. @@ -177,16 +178,15 @@ c obtain time spacing else write(fnpo,*) 'Time Spacing is NOT Correct! : dt0 > dt' write(fnpo,*) 'dt=',dt - write(fnpo,*) 'dt0=',dt0 + write(fnpo,*) 'dt0=',dt0 endif if (merr.ne.0) goto 9999 c----data read & input---------- + read(fnum,*) jjs,(parain(i),i=1,13) !!!!! elena : js is defined as REAL in the file - read(fnum,*) jjs,(parain(i),i=1,13) js=int(jjs) -! read(fnum,*) js,(parain(i),i=1,13) parain(3)=parain(3)*float(idprop) !--REV0318 xin=minval(minloc(abs(parain(9)-x))) if(parain(9).lt.x(xin).and.idprop.eq.1) xin=xin-1 !--REV0318 @@ -196,7 +196,7 @@ c----data read & input---------- parainref(1:2,ns)=parain(9:10) parainref(3,ns)=x(xin) parainref(4,ns)=parain(3) !km/s - + call datin_varied3(ro,pr,vx,vy,vz,bx,by,bz & ,vstar,ns,ix,parain,xin,paranm ,idprop,x) !--REV0318 @@ -218,20 +218,26 @@ c solve hydrodynamic equations & ,gx,dsc,scm,dv,rr,rrm,drr,dx,ix) call bnd3_2(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix,idprop,x) !--REV0318 c----------------------------------------------------------------------| -c data output - +c data output + c---monitoring - if(minval(pr).lt.0.) write(6,*) ns,minval(pr) - if(minval(pr).lt.0.) stop 'Pr<0' + if(minval(pr).lt.0.) write(6,*) ns + if(minval(pr).lt.0.) write(6,*) minval(pr) + if(minval(pr).lt.0.) stop 'Pr<0' if(idprop.eq.1.and.minval(vx).lt.0.) stop 'Vx opposite' - if(idprop.eq.-1.and.maxval(vx).gt.0.) stop 'Vx opposite' - + if(idprop.eq.-1.and.maxval(vx).gt.0.) then + write(6,1115) vx + stop 'Vx opposite' + endif +1115 format(15f10.7) c----first dt estimation if (mod(int(ns),int(touts)).eq.0) then js=js0+ns*dtr *float(idprop) !--REV0318 xout=minval(minloc(abs(parain(11)-x))) dtp1=(parain(11)-parain(9))*1.5d11/400./1.d3 - + +! write(6,*) 'first dt estimation', dtp1 + if((js-dtp1.le.js0.and.idprop.eq.1) !--REV0318 1 .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then !--REV0318 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 else itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref))) !--REV0318 -!!!!! elena + !!!!! elena if(parainref(4,itin).eq.0.0)then write(6,*) '!!!second dt',itin, ns, parain(12), parain(10) itin = ns endif c----second dt estimation - dtp1=(parain(11)-parain(9))*1.5d11/parainref(4,itin)/1.d3 - + dtp1=(parain(11)-parain(9))*1.5d11/parainref(4,itin)/1.d3 + if((js-dtp1.le.js0.and.idprop.eq.1) 1 .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then write(fnum2,704) js,1.,1.,100.,1.,1.,1.d-5,1.d-5,1.d-5 & ,0.d0,float(js0) goto 1001 endif + itin_tmp=itin itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref))) !!!!! elena if (parainref(4,itin).eq.0.0) itin=itin_tmp - c----output dang0=-(parain(12)-parainref(2,itin)) dang0=mod(dang0,360.) if(dang0.ge.180.) dang0=dang0-360. + if(angref0.ge.0.) & dtp2=(parain(11)-x(xout))*1.5d11/vx(xout)/400.d3 & -(parainref(1,itin)-parainref(3,itin))*1.5d11 !x_datain-x_in - & /parainref(4,itin)/1.d3 - !!!! elena + & /parainref(4,itin)/1.d3 + !!!! elena if (parainref(4,itin).eq.0.0) then ! NaN write(6,*) '!!!!!! NaN in dtp2 ', angref0 write(6,*) parain(11), x(xout),vx(xout) write(6,*) ns, itin, xout write(6,*) parain(12), parain(10) stop - endif - + endif + write(fnum2,704) js,ro(xout)*paranm(1) & ,pr(xout)*1937./ro(xout)*paranm(2) & ,vx(xout)*paranm(3)*float(idprop) !--REV0318 @@ -285,6 +292,7 @@ c----output & ,bz(xout)*paranm(8)*16.34 & ,dtp2*float(idprop),jsinref(itin)-parainref(6,itin) endif + 704 format(i15.1, 9e18.8,f18.1) 1001 continue @@ -670,7 +678,7 @@ c----------------------------------------------------------------------| pi = acos(-1.0d0) pi4i=0.25/pi dtmin=2.0e-10 -! safety=0.8 !cf. 0.4 --original +! safety=0.8 !cf. 0.4 --original safety=0.4 c----------------------------------------------------------------------| dt=1.e20 @@ -689,6 +697,7 @@ c----------------------------------------------------------------------| dt=dtq(i) endif enddo +! write(6,*) minval(pr) c----------------------------------------------------------------------| c write the point where dt is smaller than critical value c----------------------------------------------------------------------| diff --git a/Sources/Prop_code/sub_indat_p.f b/Sources/Prop_code/sub_indat_p.f index 2ce340f..0495ae1 100644 --- a/Sources/Prop_code/sub_indat_p.f +++ b/Sources/Prop_code/sub_indat_p.f @@ -29,9 +29,8 @@ c======================================================================| open(ifn1,file=fnin,status='old') do i=1,nnin read(ifn1,*) js1,(parain(j),j=1,12) - !!!!!!! elena : Dimensions mismatch (13/12) - datain(1:12,i)=parain(:) -! datain(:,i)=parain(:) +!!!!!!! elena : Dimensions mismatch (13/12) + datain(1:12,i)=parain(:) jsin(i)=js1 enddo dtin=(jsin(2)-jsin(1))*float(idprop) !--REV0318 diff --git a/Sources/Prop_code/sub_slct_p.f b/Sources/Prop_code/sub_slct_p.f index f9b0ace..21e0cf0 100644 --- a/Sources/Prop_code/sub_slct_p.f +++ b/Sources/Prop_code/sub_slct_p.f @@ -7,7 +7,7 @@ c======================================================================| 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) +!!!! 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 @@ -84,7 +84,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 @@ -95,6 +96,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) @@ -208,7 +210,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 @@ -220,41 +227,46 @@ 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 jsout(2,idstep)=jswout(i) - jsout(3,idstep)=float(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) + jsout(5,idstep2)=float(i) endif - enddo - + 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))-1:int(jsout(5,idstep*2))) + dangw(int(jsout(3,idstep)):int(jsout(5,idstep*2))) 1 =dangw(int(jsout(3,idstep))-1: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 + + 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 +!!!!!! 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 -- libgit2 0.21.2