istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

modified density smoothing at the boundary

Browse Source
This commit is contained in:
Daniela Farina 2015-07-10 15:29:44 +00:00
parent 87de4c9cc2
commit 900a51a08c
2 changed files with 58 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
src/gray.f
View File

@ -1972,10 +1972,10 @@ c
c
c
subroutine profdata
use const_and_precisions, only : wp_
use const_and_precisions, only : wp_,zero
use graydata_flags, only : iprof,iscal,nprof
use graydata_par, only : psdbnd,factb,factt,factn
use interp_eqprof, only : nsfd,npp,psnpp,aa,bb,cc,
use interp_eqprof, only : nsfd,npp,psnpp,denpp,ddenpp,d2denpp,
. psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec
use simplespline, only : difcsn
use dierckx, only : curfit,splev,splder
@ -1984,7 +1984,8 @@ c local variables
integer :: ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2
integer, dimension(:), allocatable :: iwrkf
real(wp_) :: aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi,
. teradi,deradi,zfci,xb,xe,sspl,dnpp,ddnpp,d2dnpp,dpsb,fp
. teradi,deradi,zfci,xb,xe,sspl,dpsb,fp
real(wp_) :: xnv,ynv
real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd,densi,
. ddensi,d2densi
c
@ -2033,54 +2034,52 @@ c
zfc(i)=zfci
wf(i)=1.0_wp_
end do
c
c spline approximation of temperature and Zeff
c
iopt=0
call difcsn(psrad,terad,npp,npp,iopt,ct,ier)
c
iopt=0
call difcsn(psrad,zfc,npp,npp,iopt,cz,ier)
c
c spline approximation of density
c
iopt=0
xb=0.0_wp_
xe=psrad(npp)
kspl=3
sspl=.001_wp_
c
call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd,
. tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
c
call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
nu=1
call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
dnpp=densi(npp)
ddnpp=ddensi(npp)
c
nu=2
call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
d2dnpp=d2densi(npp)
if(derad(npp).eq.0.0_wp_) then
psdbnd=psrad(npp)
else
psnpp=psrad(npp)
dpsb=-psnpp+psdbnd
nn=3
nn1=nn+1
nn2=nn+2
aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2)
aa=aa/(-dpsb)**nn/2.0_wp_
bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2)
bb=bb/(-dpsb)**nn1
cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2)
cc=cc/(-dpsb)**nn2/2.0_wp_
denpp=densi(npp)
ddenpp=ddensi(npp)
d2denpp=d2densi(npp)
psdbnd=psnpp
if(ddenpp.lt.0.0_wp_) then
xnv=psnpp-ddenpp/d2denpp
ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
if(d2denpp.gt.zero.and.ynv.ge.zero) then
psdbnd=xnv
else
psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
end if
c
print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
print*,psdbnd
end if
c
end if
deallocate(iwrkf,wrkf,wf,densi,ddensi,d2densi,wrkfd)
return
@ -4778,7 +4777,7 @@ c
use graydata_flags, only : iprof
use graydata_par, only : psdbnd
use graydata_anequil, only : dens0,aln1,aln2
use interp_eqprof, only : psnpp,aad=>aa,bbd=>bb,ccd=>cc,tfn,nsfd,
use interp_eqprof, only : psnpp,denpp,ddenpp,d2denpp,tfn,nsfd,
. cfn,npp
use dierckx, only : splev,splder
implicit none
@ -4786,7 +4785,7 @@ c arguments
real(wp_) :: arg
c local variables
integer :: ier,nn,nn1,nn2,nu
real(wp_) :: profd,dprofd,dpsib
real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
real(wp_), dimension(1) :: xxs,ffs
real(wp_), dimension(npp+4) :: wrkfd
c common/external functions/variables
@ -4809,16 +4808,20 @@ c
ddens=dens0*dprofd
else
if(arg.le.psdbnd.and.arg.gt.psnpp) then
c
c cubic interpolation for 1 < psi < psdbnd
c
nn=3
nn1=nn+1
nn2=nn+2
dpsib=arg-psdbnd
dens=aad*dpsib**nn+bbd*dpsib**nn1+ccd*dpsib**nn2
ddens=nn*aad*dpsib**(nn-1)+nn1*bbd*dpsib**nn
. +nn2*ccd*dpsib**nn1
c smooth interpolation for psnpp < psi < psdbnd
c dens = fp * fh
c fp: parabola matched at psi=psnpp with given profile density
c fh=(1-t)^3(1+3t+6t^2) is a smoothing function:
c fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1
fp=denpp+(arg-psnpp)*ddenpp+0.5_wp_*(arg-psnpp)**2*d2denpp
dfp=ddenpp+(arg-psnpp)*d2denpp
tt=(arg-psnpp)/(psdbnd-psnpp)
fh=(1.0_wp_-tt)**3*(1.0_wp_+3.0_wp_*tt+6.0_wp_*tt*tt)
dfh=-30.0_wp_*(1.0_wp_-tt)**2*tt*tt/(psdbnd-psnpp)
dens=fp*fh
ddens=dfp*fh+fp*dfh
else
xxs(1)=arg
ier=0
@ -5971,7 +5974,9 @@ c calculation of dP and dI over radial grid
if(j > 1) kkk=nrayth
do k=1,kkk
ii=iiv(j,k)
if (ii < nstep .and. psjki(j,k,ii+1) /= zero) ii=ii+1
if (ii < nstep ) then
if(psjki(j,k,ii+1) /= zero) ii=ii+1
end if
xxi=zero
ypt=zero
yamp=zero
@ -6195,7 +6200,7 @@ c arguments
real(wp_), dimension(nd) :: xx,yy
real(wp_), intent(out) :: xpk,ypk,dxxe
c local variables
integer :: imn,imx,ipk,ie1,ie2
integer :: imn,imx,ipk,ie
real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2
real(wp_) :: ypkp,ypkm
c
@ -6222,15 +6227,15 @@ c
xpk=xpkp
ypk=ypkp
yye=ypk*emn1
call locatex(yy,nd,1,ipk,yye,ie1)
if(ie1.gt.0.and.ie1.lt.nd) then
call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1)
call locatex(yy,nd,1,ipk,yye,ie)
if(ie.gt.0.and.ie.lt.nd) then
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1)
else
rte1=0.0_wp_
end if
call locatex(yy,nd,ipk,nd,yye,ie2)
if(ie2.gt.0.and.ie2.lt.nd) then
call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2)
call locatex(yy,nd,ipk,nd,yye,ie)
if(ie.gt.0.and.ie.lt.nd) then
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
else
rte2=0.0_wp_
end if
@ -6240,8 +6245,8 @@ c
ypk=yy(2)
rte1=0.0_wp_
yye=ypk*emn1
call locate(yy,nd,yye,ie2)
call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
call locate(yy,nd,yye,ie)
call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2)
end if
dxxe=rte2-rte1
if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe

2
src/interp_eqprof.f90
View File

@ -16,7 +16,7 @@ module interp_eqprof
! profdata
INTEGER, SAVE :: npp,nsfd
REAL(wp_), SAVE :: psnpp,aa,bb,cc
REAL(wp_), SAVE :: psnpp,denpp,ddenpp,d2denpp
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psrad,derad,terad,zfc,tfn,cfn
REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz