added missing file
This commit is contained in:
parent
c6c1395cff
commit
dcc199b336
162
src/coreprofiles.f90
Normal file
162
src/coreprofiles.f90
Normal file
@ -0,0 +1,162 @@
|
||||
module coreprofiles
|
||||
use const_and_precisions, only : wp_
|
||||
implicit none
|
||||
|
||||
INTEGER, SAVE :: npp,nsfd
|
||||
REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
|
||||
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psrad,derad,terad,zfc
|
||||
REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn
|
||||
REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
|
||||
|
||||
contains
|
||||
|
||||
subroutine alloc_profvec(ier)
|
||||
implicit none
|
||||
integer, intent(out) :: ier
|
||||
integer :: npest
|
||||
|
||||
ier=0
|
||||
if(npp.le.0) then
|
||||
ier = -1
|
||||
return
|
||||
end if
|
||||
|
||||
npest=npp+4
|
||||
|
||||
call dealloc_profvec
|
||||
allocate(psrad(npp),terad(npp),derad(npp),zfc(npp),ct(npp,4), &
|
||||
cz(npp,4),tfn(npest),cfn(npest), &
|
||||
stat=ier)
|
||||
if (ier/=0) call dealloc_profvec
|
||||
end subroutine alloc_profvec
|
||||
|
||||
subroutine dealloc_profvec
|
||||
implicit none
|
||||
if(allocated(psrad)) deallocate(psrad)
|
||||
if(allocated(terad)) deallocate(terad)
|
||||
if(allocated(derad)) deallocate(derad)
|
||||
if(allocated(zfc)) deallocate(zfc)
|
||||
if(allocated(ct)) deallocate(ct)
|
||||
if(allocated(cz)) deallocate(cz)
|
||||
if(allocated(tfn)) deallocate(tfn)
|
||||
if(allocated(cfn)) deallocate(cfn)
|
||||
|
||||
end subroutine dealloc_profvec
|
||||
|
||||
subroutine density(psin,dens,ddens)
|
||||
use const_and_precisions, only : wp_
|
||||
use graydata_flags, only : iprof
|
||||
use graydata_anequil, only : dens0,aln1,aln2
|
||||
use dierckx, only : splev,splder
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_), intent(out) :: dens,ddens
|
||||
! local variables
|
||||
integer, parameter :: nn=3, nn1=nn+1, nn2=nn+2
|
||||
integer :: ier,nu
|
||||
real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
|
||||
real(wp_), dimension(1) :: xxs,ffs
|
||||
real(wp_), dimension(npp+4) :: wrkfd
|
||||
|
||||
!
|
||||
! computation of density [10^19 m^-3] and derivative wrt psi
|
||||
!
|
||||
dens=0.0_wp_
|
||||
ddens=0.0_wp_
|
||||
if(psin.ge.psdbnd.or.psin.lt.0.0_wp_) return
|
||||
!
|
||||
if(iprof.eq.0) then
|
||||
if(psin.gt.1.0_wp_) return
|
||||
profd=(1.0_wp_-psin**aln1)**aln2
|
||||
dens=dens0*profd
|
||||
dprofd=-aln1*aln2*psin**(aln1-1.0_wp_) &
|
||||
*(1.0_wp_-psin**aln1)**(aln2-1.0_wp_)
|
||||
ddens=dens0*dprofd
|
||||
else
|
||||
if(psin.gt.psnpp) then
|
||||
!
|
||||
! smooth interpolation for psnpp < psi < psdbnd
|
||||
! dens = fp * fh
|
||||
! fp: parabola matched at psi=psnpp with given profile density
|
||||
! fh=(1-t)^3(1+3t+6t^2) is a smoothing function:
|
||||
! fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1
|
||||
!
|
||||
dpsib=psin-psnpp
|
||||
fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
|
||||
dfp=ddenpp+dpsib*d2denpp
|
||||
tt=dpsib/(psdbnd-psnpp)
|
||||
fh=(1.0_wp_-tt)**3*(1.0_wp_+3.0_wp_*tt+6.0_wp_*tt**2)
|
||||
dfh=-30.0_wp_*(1.0_wp_-tt)**2*tt**2/(psdbnd-psnpp)
|
||||
dens=fp*fh
|
||||
ddens=dfp*fh+fp*dfh
|
||||
else
|
||||
xxs(1)=psin
|
||||
ier=0
|
||||
call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
|
||||
dens=ffs(1)
|
||||
nu=1
|
||||
ier=0
|
||||
call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
|
||||
ddens=ffs(1)
|
||||
if(ier.gt.0) print*,ier
|
||||
if(abs(dens).lt.1.0e-10_wp_) dens=0.0_wp_
|
||||
end if
|
||||
if(dens.lt.0.0_wp_) print*,' DENSITY NEGATIVE',dens
|
||||
end if
|
||||
end subroutine density
|
||||
|
||||
function temp(psin)
|
||||
use const_and_precisions, only : wp_,zero,one
|
||||
use graydata_flags, only : iprof
|
||||
use graydata_anequil, only : te0,dte0,alt1,alt2
|
||||
use utils, only : locate
|
||||
use simplespline, only :spli
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: temp
|
||||
! local variables
|
||||
integer :: k
|
||||
real(wp_) :: proft,dps
|
||||
|
||||
temp=zero
|
||||
if(psin.ge.one.or.psin.lt.zero) return
|
||||
if(iprof.eq.0) then
|
||||
proft=(1.0_wp_-psin**alt1)**alt2
|
||||
temp=(te0-dte0)*proft+dte0
|
||||
else
|
||||
call locate(psrad,npp,psin,k)
|
||||
k=max(1,min(k,npp-1))
|
||||
dps=psin-psrad(k)
|
||||
temp=spli(ct,npp,k,dps)
|
||||
endif
|
||||
end function temp
|
||||
|
||||
function fzeff(psin)
|
||||
use const_and_precisions, only : wp_,zero,one
|
||||
use graydata_flags, only : iprof
|
||||
use graydata_anequil, only : zeffan
|
||||
use utils, only : locate
|
||||
use simplespline, only :spli
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_) :: fzeff
|
||||
! local variables
|
||||
integer :: k
|
||||
real(wp_) :: dps
|
||||
|
||||
fzeff=one
|
||||
if(psin.gt.one.or.psin.lt.zero) return
|
||||
if(iprof.eq.0) then
|
||||
fzeff=zeffan
|
||||
else
|
||||
call locate(psrad,npp,psin,k)
|
||||
k=max(1,min(k,npp-1))
|
||||
dps=psin-psrad(k)
|
||||
fzeff=spli(cz,npp,k,dps)
|
||||
endif
|
||||
end function fzeff
|
||||
|
||||
end module coreprofiles
|
||||
Loading…
Reference in New Issue
Block a user