diff --git a/input/beam.data b/input/beam.data old mode 100755 new mode 100644 diff --git a/input/gray_params.data b/input/gray_params.data old mode 100755 new mode 100644 diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 old mode 100755 new mode 100644 diff --git a/src/graycore.f90 b/src/graycore.f90 index 7349ff3..505b6ae 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -8,7 +8,7 @@ contains eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & psipol0,chipol0, dpdv,jcd,pabs,icd, outp,rtrp,hcdp,ierr, rhout) - use const_and_precisions, only : zero, one + use const_and_precisions, only : zero, one, degree use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff use dispersion, only : expinit use gray_params, only : eqparam_type, prfparam_type, outparam_type, & @@ -59,7 +59,7 @@ contains real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx - real(wp_), dimension(3) :: xv,anv0,anv + real(wp_), dimension(3) :: xv,anv0,anv,bv real(wp_), dimension(:,:), allocatable :: yw,ypw,gri real(wp_), dimension(:,:,:), allocatable :: xc,du1,ggri integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1 @@ -115,6 +115,8 @@ contains ! print ne, Te, q, Jphi versus psi, rhop, rhot call print_bres(bres) call print_prof + call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), & + sin(beta0*degree)) ! ======= pre-proc prints END ====== ! ======= main loop BEGIN ====== @@ -152,7 +154,7 @@ contains xv = yw(1:3,jk) anv = yw(4:6,jk) call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & - psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn) + psinv,dens,btot,bv,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn) ! update global error code and print message if (ierrn/=0) then ierr = ior(ierr,ierrn) @@ -201,8 +203,8 @@ contains dids0(jk)=dids ccci0(jk)=ccci(jk,i) - call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr, & - anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) + call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,bv,ak0,anpl,anpr, & + anv,anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) end do @@ -596,8 +598,9 @@ contains ddr = anx**2 + any**2 + anz**2 - an20 ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) - call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,ak0,zero,zero,zero, & - zero,zero,zero,zero,zero,0,0,1,ddr,ddi) ! st=0, index_rt=1, Btot=0, psin=-1 + call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,(/zero,zero,zero/), & + ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, & + 0,0,1,ddr,ddi) ! st=0, index_rt=1, B=0, N=0, psin=-1 end do end subroutine ic_gb @@ -666,7 +669,7 @@ contains subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & - xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) + bv,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use errcodes, only : pnpl @@ -679,10 +682,11 @@ contains real(wp_), dimension(6), intent(out) :: dery real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm + real(wp_), dimension(3), intent(out) :: bv integer, intent(out) :: ierr ! local variables real(wp_) :: gr2,ajphi - real(wp_), dimension(3) :: dgr2,bv,derxg,deryg + real(wp_), dimension(3) :: dgr2,derxg,deryg real(wp_), dimension(3,3) :: derbv real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ @@ -1648,8 +1652,8 @@ bb: do write(uprj0+1,*) ' #sst j k xt yt zt rt' write(uwbm,*) ' #sst w1 w2' write(udisp,*) ' #sst Dr_Nr Di_Nr' - write(ucenr,*) ' #sst R z phi psin rhot ne Te Btot Nperp Npl ki '// & - 'alpha tau Pt dIds nhmax iohkw index_rt ddr' + write(ucenr,*) ' #sst R z phi psin rhot ne Te Btot Bx By Bx Nperp Npl '// & + 'Nx Ny Nz ki alpha tau Pt dIds nhmax iohkw index_rt ddr' write(uoutr,*) ' #i k sst x y R z psin tau Npl alpha index_rt' write(upec,*) ' #rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt' write(usumm,*) ' #Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // & @@ -1742,6 +1746,59 @@ bb: do + subroutine print_maps(bres,xgcn,r0,anpl0) + use const_and_precisions, only : wp_ + use gray_params, only : iequil + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, & + equinum_fpol, nq + use coreprofiles, only : density, temp + use units, only : umaps + implicit none +! arguments + real(wp_), intent(in) :: bres,xgcn,r0,anpl0 +! local variables + integer :: j,k + real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl + real(wp_), dimension(nq) :: r, z + + + dr = (rmxm-rmnm)/(nq-1) + dz = (zmxm-zmnm)/(nq-1) + do j=1,nq + r(j) = rmnm + dr*(j-1) + z(j) = zmnm + dz*(j-1) + end do + + write(umaps,*)'#R z psin Br Bphi Bz Btot ne Te X Y Npl' + do j=1,nq + rj=r(j) + anpl=anpl0*r0/rj + do k=1,nq + zk=z(k) + if (iequil < 2) then + call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br) + else + call equinum_psi(rj,zk,psinv=psin,dpsidr=bz,dpsidz=br) + call equinum_fpol(psin,fpolv=bphi) + end if + br = -br/rj + bphi = bphi/rj + bz = bz/rj + btot = sqrt(br**2+bphi**2+bz**2) + yg = btot/bres + te = temp(psin) + call density(psin,ne,dne) + xg = xgcn*ne + write(umaps,'(12(x,e12.5))') rj,zk,psin,br,bphi,bz,btot,ne,te,xg,yg,anpl + enddo + write(umaps,*) + enddo + + end subroutine print_maps + + + + subroutine print_surfq(qval) use const_and_precisions, only : wp_, one use equilibrium, only : psinr,nq,fq,frhotor,rmaxis,zmaxis, & @@ -1844,8 +1901,8 @@ bb: do - subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, & - dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) + subroutine print_output(i,jk,st,qj,xv,psinv,btot,bv,ak0,anpl,anpr,anv, & + anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) use const_and_precisions, only : degree,zero,one use equilibrium, only : frhotor use gray_params, only : istpl0 @@ -1854,7 +1911,7 @@ bb: do implicit none ! arguments integer, intent(in) :: i,jk,nhf,iokhawa,index_rt - real(wp_), dimension(3), intent(in) :: xv + real(wp_), dimension(3), intent(in) :: xv,bv,anv real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi ! local variables @@ -1879,8 +1936,8 @@ bb: do pt=exp(-tau) didsn=dids*1.0e2_wp_/qj - write(ucenr,'(16(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, & - psinv,rhot,dens,tekev,btot,anpr,anpl,akim,alpha,tau,pt,didsn, & + write(ucenr,'(22(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, & + psinv,rhot,dens,tekev,btot,bv,anpr,anpl,anv,akim,alpha,tau,pt,didsn, & nhf,iokhawa,index_rt,ddr end if diff --git a/src/units.f90 b/src/units.f90 index 90580bd..ce30d71 100644 --- a/src/units.f90 +++ b/src/units.f90 @@ -1,13 +1,15 @@ module units ! STANDARD -! integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99 -! integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71 -! integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33 -! integer, parameter :: udisp = 17, upec = 48, usumm = 7 + integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99 + integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71 + integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33 + integer, parameter :: udisp = 17, upec = 48, usumm = 7, umaps = 72 + integer, parameter :: uhead = 78 ! JINTRAC - integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644 - integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631 - integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633 - integer, parameter :: udisp =617, upec =648, usumm =607, uhead =638 +! integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644 +! integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631 +! integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633 +! integer, parameter :: udisp =617, upec =648, usumm =607, umaps =632 +! integer, parameter :: uhead =638 end module units \ No newline at end of file