V3FIT
•All Classes Namespaces Files Functions Variables Enumerations Macros Pages
rdeqdsk.f
1  subroutine rdeqdsk(g1)
2  USE precision
3  USE gfile
4  USE mapg_mod
5 c-----------------------------------------------------------------------
6 c changed dpsi to dpsiv: rlm 7/3/96
7 c-----------------------------------------------------------------------
8  IMPLICIT NONE
9  TYPE(geqdsk) :: g1
10  CHARACTER*8 ntitle(5)
11  CHARACTER*6 dat
12  INTEGER :: ipestg, limitr, mbdry
13  REAL(rprec) :: xdim,zdim,rc,redge,zmid
14  REAL(rprec) :: btor
15  REAL(rprec) :: totcur,psimx(2),xax(2)
16  REAL(rprec) :: zax(2),psisep,xsep,zsep
17  REAL(rprec) :: dpsiv
18  REAL(rprec) :: zmin
19  INTEGER :: i,j,isave
20  neqdsk=118
21  ntitle=" "
22 c
23  open (neqdsk,file=filename,status='old', iostat=isave)
24  if(isave /= 0) then
25  print*,"iostat=",isave," opening file=",trim(filename),
26  . " on unit=",neqdsk
27  stop 'EQDSK'
28  endif
29 c read eqdsk file
30  write(6,'("begin read eqdsk")')
31  READ(neqdsk,200,err=211)(ntitle(i),i=1,5),dat,ipestg,nx,nz
32  GOTO 212
33  211 print*,' READ ERROR :ipestg,nx,nz=',ipestg,nx,nz
34  stop 'eqdsk 1st line'
35  212 continue
36  200 format(6a8,3i4)
37  call eqdsk_allocate(nx,nz)
38  read(neqdsk,300)xdim,zdim,rc,redge,zmid
39  read(neqdsk,300)xaxis,zaxis,psiaxis,psilim,btor
40  read(neqdsk,300)totcur,psimx(1),psimx(2),xax(1),xax(2)
41  read(neqdsk,300)zax(1),zax(2),psisep,xsep,zsep
42  read(neqdsk,300)(sf(i),i=1,nx)
43  read(neqdsk,300)(sp(i),i=1,nx)
44  read(neqdsk,300)(sffp(i),i=1,nx)
45  read(neqdsk,300)(spp(i),i=1,nx)
46  read(neqdsk,300)((psixz(i,j),i=1,nx),j=1,nz)
47  300 format(5e16.9)
48  kvtor=0
49  nmass=0
50  read(neqdsk,300,end=500) (qpsi(i),i=1,nx)
51  read(neqdsk,'(2i5)') nbndry,nlim
52  limitr=nlim; mbdry=nbndry
53  allocate(xbndry(nbndry),zbndry(nbndry))
54  read(neqdsk,300) (xbndry(i),zbndry(i),i=1,nbndry)
55  allocate(xlim(nlim),zlim(nlim))
56  read(neqdsk,300) (xlim(i),zlim(i),i=1,nlim)
57  write(6,'("finished eqdsk")')
58  if(rotate.ne.0) then
59  read(neqdsk,'(i5,e16.9,i5)',end=500,err=500) kvtor,rvtor,nmass
60  if(kvtor.gt.0) then
61  read(neqdsk,300) (pressw(i),i=1,nx)
62  read(neqdsk,300) (pwprim(i),i=1,nx)
63 c guard against negative or zero pressw
64  isave=1
65  do i=1,nx
66  if(pressw(i).le.0.) then
67  isave=i
68  go to 75
69  endif
70  end do
71  75 continue
72  if(isave.ne.1) then
73  do i=isave,nx
74  pressw(i)=pressw(isave-1)+(i-isave)/(nx-isave)*
75  $ (1.-pressw(isave-1))
76  pwprim(i)=(pressw(isave-1)-1.)/
77  $ (psilim-psiaxis)*(nx-isave)/(nx-1)
78  end do
79  endif
80  endif
81  if(nmass.gt.0) then
82  read(neqdsk,300) (rho0(i),i=1,nx)
83  dpsiv=(psilim-psiaxis)/(nx-1.)
84  do i=2,nx-1
85  rho0p(i)=(rho0(i+1)-rho0(i-1))/(2.*dpsiv)
86  end do
87  rho0p(1)=(-3.*rho0(1)+4.*rho0(2)-rho0(3))/(2.*dpsiv)
88  rho0p(nx)=(3.*rho0(nx)-4.*rho0(nx-1)+rho0(nx-2))/(2.*dpsiv)
89  endif
90  endif
91 c generate x,z gridd
92  500 continue
93  dx=xdim/(nx-1.)
94  dz=zdim/(nz-1.)
95  do i=0,nx - 1
96  xgrid(i + 1)=redge + dx*i
97 ! xgrid(i)=redge+(i-1.)*dx
98  end do
99 ! nz is an integer. Need the intger devide here so the floor comes for free.
100  zmin = zmid - zdim/2.0
101  do i=0,nz - 1
102  zgrid(i + 1) = zmin + dz*i
103 ! zgrid(i)=-0.5*zdim+(i-1.)*dz
104  end do
105 c if sf is negative change it's sign
106 c I don't know what the purpose of a negative B field is.
107  if(sf(nx).lt.0.) then
108  do i=1,nx
109  sf(i)=-sf(i)
110  end do
111  endif
112  close(neqdsk)
113  call getbpsq(psixz,nxd,nzd,xgrid,dx,dz,nx,nz,bpsq)
114 ! load g
115  write(6,'("begin loading geqdsk")')
116  call gfile_allocate(nx,nz,mbdry,limitr,g1)
117  read(ntitle(4)(3:8),fmt='(i10)')i
118  g1%shot=i
119  read(ntitle(5)(1:6),fmt='(i10)')j
120  g1%time=j
121  g1%source=trim(filename)
122  g1%xdim=xdim
123  g1%zdim=zdim
124  g1%mw=nx
125  g1%mh=nz
126  g1%rzero=rc
127  g1%rgrid1=redge
128  g1%zmid=zmid
129  g1%rmaxis=xaxis
130  g1%zmaxis=zaxis
131  g1%ssimag=psiaxis
132  g1%ssibry=psilim
133  g1%bcentr=btor
134  g1%cpasma=totcur
135  g1%fpol=sf
136  g1%pres=sp
137  g1%ffprim=sffp
138  g1%pprime=spp
139  g1%psirz=psixz(1:nx,1:nz)
140  g1%qpsi=qpsi
141  g1%nbdry=mbdry
142  g1%limitr=limitr
143  g1%rbdry(1:mbdry)=xbndry(1:mbdry)
144  g1%zbdry(1:mbdry)=zbndry(1:mbdry)
145  g1%xlim(1:limitr)=xlim(1:limitr)
146  g1%ylim(1:limitr)=zlim(1:limitr)
147  write(6,'("end loading geqdsk")')
148  end subroutine rdeqdsk
149 
150  subroutine eqdsk_allocate(mw,mh)
151  USE mapg_mod
152  INTEGER, INTENT(IN) :: mw,mh
153  IF(ALLOCATED(psixz))stop 'allocation error in reqdsk'
154  ALLOCATE(psixz(mw,mh),bpsq(mw,mh),stat=istat)
155  if(istat.ne.0)print*,"STAT=",istat
156  ALLOCATE(sp(mw),spp(mw),sf(mw),sffp(mw),qpsi(mw)
157  . ,stat=istat)
158  if(istat.ne.0)print*,"STAT=",istat
159  ALLOCATE(xgrid(mw),zgrid(mw))
160  ALLOCATE(pressw(mw),pwprim(mw),rho0(mw),rho0p(mw),stat=istat)
161  if(istat.ne.0)print*,"STAT=",istat
162 ! parameter(nxd=129,nzd=129,nxzd=6*(nxd+nzd),nh2=2*nzd,
163 ! $ nwk=2*nxd*nzd+nh2)
164 ! nxd=mw; nzd=mh;nxzd=6*(nxd+nzd);nh2=2*nzd;nwk=2*nxd*nzd+nh2
165  nxd=mw; nzd=mh;nxzd=2000;nh2=2*nzd;nwk=2*nxd*nzd+nh2
166 ! allocate splining arrays for eqdsk
167  ALLOCATE(csplpsi(4,nxd,nzd),stat=istat)
168  if(istat.ne.0)print*,"STAT=",istat
169 ! initialize
170  pressw=0;pwprim=0;rho0=0;rho0p=0;csplpsi=0;xgrid=0;zgrid=0
171  psixz=0;bpsq=0;sp=0;spp=0;sf=0;sffp=0;qpsi=0
172  end subroutine eqdsk_allocate
173 
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17