diff --git a/Makefile b/Makefile index 1e94bc7..c57d9a0 100644 --- a/Makefile +++ b/Makefile @@ -1,112 +1,110 @@ -# Makefile for JETTO fortran 90 library 'gray' -# Author: Lorenzo Figini (figini@ifp.cnr.it) +# Makefile for JETTO fortran 90 library 'gray' +# Author: Lorenzo Figini (figini@ifp.cnr.it) # # Derived from -# Makefile for JETTO fortran 77 library 'frantic' -# Author: Derek Harting (d.harting@fz-juelich.de) +# Makefile for JETTO fortran 77 library 'frantic' +# Author: Derek Harting (d.harting@fz-juelich.de) # and # Makefile for ITCequ Fortran 90 library - definitions # P. Strand, elfps@chalmers.se # G. Corrigan, gcor@jet.uk # D. Harting, d.harting@fz-juelich.de - -# Set the environment from the top-level Makefile of JETTO -# -------------------------------------------------------- -include ../include.mk +# Set the environment from the top-level Makefile of JETTO +# -------------------------------------------------------- +include ../include.mk # library name -# ------------ +# ------------ LIBNAME=$(JLIBDIR)/libgray.a -# List source and object files -# ---------------------------- -FSRC=$(wildcard *.f *.f90) -FOBJ0=$(FSRC:.f=.o) -FOBJ=$(FOBJ0:.f90=.o) - -######################################## -# Set up compiler specific environment # -######################################## +# List source and object files +# ---------------------------- +FSRC=$(wildcard *.f *.f90) +FOBJ0=$(FSRC:.f=.o) +FOBJ=$(FOBJ0:.f90=.o) +######################################## +# Set up compiler specific environment # +######################################## # set default compiler # -------------------- FC=$(F90) -# Set compiler debug and optimization flags -# ----------------------------------------- -ifeq ("$(DBG)","") - FFLAGS= $(AUTO) -O3 -else - FFLAGS= $(AUTO) -O0 -g -endif - -# Set include directories -# ----------------------- -#INC=-I$(JCOMMON_STD) +# Set compiler debug and optimization flags +# ----------------------------------------- +ifeq ("$(DBG)","") + FFLAGS= $(AUTO) -O3 +else + FFLAGS= $(AUTO) -O0 -g +endif + +# Set include directories +# ----------------------- +#INC=-I$(JCOMMON_STD) INC= - -##################################### -# Set up RULES specific environment # -##################################### + +##################################### +# Set up RULES specific environment # +##################################### # suffixes -# -------- +# -------- .SUFFIXES: .f90 .f .o .mod -# PHONY targets -# ------------- -.PHONY: all clean realclean $(MASTER_RULES) gray - -# default target -# -------------- -all: - @($(MAKE) -C ../ all) || exit $$? - -# Set rules passed to top-level Makefile -# -------------------------------------- -$(MASTER_RULES): - @($(MAKE) -C ../ $@) || exit $$? - -# rule for libgray.a -# -------------------- -$(LIBNAME): $(FOBJ) - ar vr $(LIBNAME) $? - -# synonym for libgray.a -# ----------------------- -gray: $(LIBNAME) - -# create rule for compiling f90 files -# ----------------------------------- -%.o: %.f90 - $(FC) -c $(FFLAGS) $(INC) $< - -# create rule for compiling f files -# --------------------------------- -%.o: %.f - $(FC) -c $(FFLAGS) $(INC) $< - -# make 'realclean' option -# ----------------------- -realclean: - rm -f *.o *.mod $(LIBNAME) - -# make 'clean' option -# ------------------- -clean: - rm -f *.o +# PHONY targets +# ------------- +.PHONY: all clean realclean $(MASTER_RULES) gray + +# default target +# -------------- +all: + @($(MAKE) -C ../ all) || exit $$? + +# Set rules passed to top-level Makefile +# -------------------------------------- +$(MASTER_RULES): + @($(MAKE) -C ../ $@) || exit $$? + +# rule for libgray.a +# -------------------- +$(LIBNAME): $(FOBJ) + ar vr $(LIBNAME) $? + +# synonym for libgray.a +# ----------------------- +gray: $(LIBNAME) + +# create rule for compiling f90 files +# ----------------------------------- +%.o: %.f90 + $(FC) -c $(FFLAGS) $(INC) $< + +# create rule for compiling f files +# --------------------------------- +%.o: %.f + $(FC) -c $(FFLAGS) $(INC) $< + +# make 'realclean' option +# ----------------------- +realclean: + rm -f *.o *.mod $(LIBNAME) + +# make 'clean' option +# ------------------- +clean: + rm -f *.o # Dependencies # ------------ gray_main.o: gray-externals.o itm_types.o -gray-externals.o: green_func_p.o reflections.o +gray-externals.o: green_func_p.o reflections.o scatterspl.o green_func_p.o: const_and_precisions.o const_and_precisions.o: itm_types.o itm_constants.o itm_constants.o: itm_types.o +scatterspl.o: itm_types.o # Special rule to preprocess source file and include svn revision # --------------------------------------------------------------- DIRECTIVES = -DREVISION="'$(shell svnversion)'" gray-externals.o:gray-externals.f - $(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC) $< - + $(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC) $< diff --git a/Makefile.standalone b/Makefile.standalone index 131d69d..775a962 100644 --- a/Makefile.standalone +++ b/Makefile.standalone @@ -3,8 +3,8 @@ EXE=gray LIBFILE=lib$(EXE).a # Objects list -OBJ=gray_main.o gray-externals.o grayl.o reflections.o green_func_p.o \ - const_and_precisions.o itm_constants.o itm_types.o +OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \ + green_func_p.o const_and_precisions.o itm_constants.o itm_types.o # Alternative search paths vpath %.f90 src @@ -30,10 +30,11 @@ $(LIBFILE): $(OBJ) # Dependencies on modules main.o: gray_main.o gray_main.o: gray-externals.o itm_types.o -gray-externals.o: green_func_p.o reflections.o +gray-externals.o: green_func_p.o reflections.o scatterspl.o green_func_p.o: const_and_precisions.o const_and_precisions.o: itm_types.o itm_constants.o itm_constants.o: itm_types.o +scatterspl.o: itm_types.o # General object compilation command %.o: %.f90 diff --git a/src/gray-externals.f b/src/gray-externals.f index 8f2347d..e4629bb 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -1023,6 +1023,7 @@ c subroutine equidata(ijetto,mr,mz,r,z,psin2d,psiaxis,psiedge, . rax,zax,nbnd,rbnd,zbnd,mrho,psijet,fpjet,qjet) use itm_constants, only : pi=>itm_pi + use reflections, only : inside implicit real*8 (a-h,o-z) integer ijetto,mr,mz,nbnd,mrho real*8 r(mr),z(mz),psin2d(mr,mz) @@ -1042,6 +1043,7 @@ c dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh) dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) + dimension rv1d(nnw*nnh),zv1d(nnw*nnh),wpsi(nnw*nnh) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh) @@ -1129,49 +1131,95 @@ c psi function c c do j=1,nz c do i=1,nr -c write(620,2021) rv(i),zv(j),psi(i,j) +c write(620,2021) rv(i),zv(j),psin(i,j) c enddo c write(620,*) ' ' c enddo - - do j=1,nz - do i=1,nr - psi(i,j)=psin(i,j)*psia - ij=nz*(i-1)+j - fvpsi(ij)=psin(i,j) - enddo - enddo c c spline fitting/interpolation of psin(i,j) and derivatives c - iopt=0 - call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, - . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, - . wrk,lwrk,iwrk,liwrk,ier) -c if ier=-1 data are fitted using sspl=0 - if(ier.eq.-1) then - sspl=0.0d0 - call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, - . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, - . wrk,lwrk,iwrk,liwrk,ier) - end if - nsrt=nsr - nszt=nsz - if (sspl.gt.0.0d0) then - call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, - . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) -c + if (ijetto.eq.1) then +c valid data on the whole grid do j=1,nz do i=1,nr +c psi(i,j)=psin(i,j)*psia ij=nz*(i-1)+j - psin(i,j)=ffvpsi(ij) - psi(i,j)=psin(i,j)*psia -c write(619,2021) rv(i),zv(j),psin(i,j) + fvpsi(ij)=psin(i,j) enddo -c write(619,*) ' ' enddo + iopt=0 + call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, + . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, + . wrk,lwrk,iwrk,liwrk,ier) +c if ier=-1 data are fitted using sspl=0 + if(ier.eq.-1) then + sspl=0.0d0 + call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, + . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, + . wrk,lwrk,iwrk,liwrk,ier) + end if + nsrt=nsr + nszt=nsz + else +c valid data only inside boundary + ij=0 + do j=1,nz + do i=1,nr + if (.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle + ij=ij+1 + rv1d(ij)=rv(i) + zv1d(ij)=zv(j) + fvpsi(ij)=psin(i,j) + wpsi(ij)=1.0d0 + end do + end do +c add boundary to the fit + if (ij+nbnd.le.nnw*nnh) then + do i=1,nbnd + ij=ij+1 + rv1d(ij)=rbnd(i) + zv1d(ij)=zbnd(i) + fvpsi(ij)=1.0d0 + wpsi(ij)=1.0d0 + end do + end if + nrz=ij +c +c fit as a scattered set of points +c + nsr=nr+4 + nsz=nz+4 + call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, + . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) +!!!!!!!!!!!! Sistemare anche dipendenze Makefile scatterspl:... +c if ier=-1 data are fitted using sspl=0 + if(ier.eq.-1) then + sspl=0.0d0 + nsr=nr+4 + nsz=nz+4 + call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, + . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) + end if + nsrt=nsr + nszt=nsz end if -2021 format(5(1x,e16.9)) +c +c re-evaluate psi on the grid using the spline (only for debug) +c +c call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, +c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) +c +c do j=1,nz +c do i=1,nr +c ij=nz*(i-1)+j +c psin(i,j)=ffvpsi(ij) +cc psi(i,j)=psin(i,j)*psia +c write(619,2021) rv(i),zv(j),psin(i,j) +c enddo +c write(619,*) ' ' +c enddo +c +c2021 format(5(1x,e16.9)) c nur=1 nuz=0 diff --git a/src/gray.f b/src/gray.f index bbb2a77..d742660 100644 --- a/src/gray.f +++ b/src/gray.f @@ -18,7 +18,7 @@ c input arguments c c ijetto Equilibrium source (1 EFIT, 2 ESCO) c If IJETTO=2, then PSIN values are valid only inside -c plasma boudary (PSIN=0 outside) +c plasma boudary (PSIN=1 outside) c mr Size of flux map grid in R direction c mz Size of flux map grid in Z direction c r R coordinates of flux map grid points [m] diff --git a/src/reflections.f90 b/src/reflections.f90 index 6e04214..1fd73d6 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -188,12 +188,17 @@ function inside(xc,yc,n,x,y) integer, dimension(n) :: jint real(r8), dimension(n) :: xint real(r8), dimension(n+1) :: xclosed,yclosed - integer :: i,nj + integer :: i,np,nj xclosed(1:n)=xc(1:n) yclosed(1:n)=yc(1:n) - xclosed(n+1)=xc(1) - yclosed(n+1)=yc(1) - call locate_unord(yclosed,n+1,y,jint,n,nj) + if (xc(1)==xc(n) .and. yc(1)==yc(n)) then + np=n + else + np=n+1 + xclosed(np)=xc(1) + yclosed(np)=yc(1) + end if + call locate_unord(yclosed(1:np),np,y,jint,n,nj) inside=.false. if (nj==0) return do i=1,nj diff --git a/src/scatterspl.f90 b/src/scatterspl.f90 new file mode 100644 index 0000000..025ab65 --- /dev/null +++ b/src/scatterspl.f90 @@ -0,0 +1,50 @@ +subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & + tx,nknt_x,ty,nknt_y,coeff,ierr) + use itm_types, only : r8 + implicit none + integer :: n + real(r8), dimension(n) :: x, y, z + real(r8), dimension(n) :: w + integer :: kspl + real(r8) :: sspl + real(r8) :: xmin, xmax, ymin, ymax + real(r8), dimension(nknt_x) :: tx + real(r8), dimension(nknt_y) :: ty + integer :: nknt_x, nknt_y + real(r8), dimension(nknt_x*nknt_y) :: coeff + integer :: ierr + + integer :: iopt + real(r8), parameter :: eps_r8=epsilon(1._r8) + real(r8) :: resid + integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest + real(r8), dimension(:), allocatable :: tx_tmp,ty_tmp + real(r8), dimension(:), allocatable :: wrk1, wrk2 + integer, dimension(:), allocatable :: iwrk + + nxest=nknt_x + nyest=nknt_y + ne = max(nxest,nyest) + allocate(tx_tmp(ne),ty_tmp(ne)) + + km = kspl+1 + u = nxest-km + v = nyest-km + b1 = kspl*min(u,v)+kspl+1 + b2 = (kspl+1)*min(u,v)+1 + lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1 + lwrk2 = u*v*(b2+1)+b2 + kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1) + allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)) + + iopt=0 + call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl,sspl, & + nxest,nyest,ne,eps_r8,nknt_x,tx_tmp,nknt_y,ty_tmp, & + coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr) + tx(1:nknt_x)=tx_tmp(1:nknt_x) + ty(1:nknt_y)=ty_tmp(1:nknt_y) + + deallocate(wrk1,wrk2,iwrk) + deallocate(tx_tmp,ty_tmp) + +end subroutine scatterspl