From 900a51a08cd0d3a1af98f9ac9f58b2839968f9bf Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Fri, 10 Jul 2015 15:29:44 +0000 Subject: [PATCH] modified density smoothing at the boundary --- src/gray.f | 109 ++++++++++++++++++++++-------------------- src/interp_eqprof.f90 | 2 +- 2 files changed, 58 insertions(+), 53 deletions(-) diff --git a/src/gray.f b/src/gray.f index 9f11f9d..644f106 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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_ + psnpp=psrad(npp) + 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 + print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv + print*,psdbnd end if -c + end if -c + 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 diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90 index babfa50..2f7bbe7 100644 --- a/src/interp_eqprof.f90 +++ b/src/interp_eqprof.f90 @@ -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