Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Molecular CRM model implementation into develop #4

Merged
merged 2 commits into from
Jan 11, 2024
Merged

Molecular CRM model implementation into develop #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
103852c
Molecular CRM model implementation
holm10 May 23, 2023
4fea923
Merge branch 'upstream-develop-h2crm' into h2crm
holm10 Jan 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 48 additions & 0 deletions aph/aph.v
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@ sgvcxc real [m^2/s] /2.e-14/# const. value of sigv_cx for issgvcxc=1;
# const. sig for issgvcxc=2, sig-v=sig*sqrt(T/mi)
# with units now [m^2]
isaphdir integer /1/ #=1 uses aphdir; =0 uses explicit rate file names
crmnfname character*120 /"crumpet_nrates.dat"/ # Name of CRM rate data file
crmefname character*120 /"crumpet_Erates.dat"/ # Name of CRM E-sink data file
aphdir character*120 # name of directory containing data files
data_directory character*120 # another dirname containing data files. This is to be be passed in

Expand Down Expand Up @@ -38,6 +40,34 @@ htlcx(0:htnt,0:htnn,0:htns-1) _real
htlqa(0:htnt,0:htnn,0:htns-1) _real
# log(htqa) where htqa is ADPAK rate parameter data for energy loss


***** Rtcrumpet:
# data from CRUMPET rate tables
cmpe integer /60/ # number of energy points for crm params
cmpd integer /15/ # number of density points for crm params
crmdiss(cmpe,cmpd) _real [m**3/s] # Mol diss rate vs. temp and dens
crmarate(cmpe,cmpd) _real [m**3/s] # a-create rate from mol vs tem, dens
crmselm(cmpe,cmpd) _real [W m**3] # el e-change due to mol diss vs temp, dens
crmsiam(cmpe,cmpd) _real [W m**3] # ia e-change due to mol diss vs temp, dens
crmspotm(cmpe,cmpd) _real [W m**3] # Epot change due to mol diss vs temp, dens
crmsrada(cmpe,cmpd) _real [W m**3] # Radiation due to mol diss vs temp, dens
crmsradm(cmpe,cmpd) _real [W m**3] # Radiation due to mol diss vs temp, dens
cekpt(cmpe) _real # natural log of temperature (eV)
crlemin real # minimum of ekpt
crlemax real # maximum of ekpt
cerefmin real [eV] # minimum temperature in table data
cerefmax real [eV] # maximum temperature in table data
cdelekpt real # interval size for ekpt in table data
cdkpt(cmpd) _real # log10 of density (/m**3)
crldmin real # minimum of dkpt
crldmax real # maximum of dkpt
cdrefmin real [/m**3] # minimum density in table data
cdrefmax real [/m**3] # maximum density in table data
cdeldkpt real # dkpt interval in atomic data tables
ctaumin real # minimum tau in rtau opt-dep table data
ctaumax real # maximum tau in trau opt-dep table data
cdeltau real # log int. size of rtau in opt-dep table dat

***** Rtdegas:
# data from DEGAS rate tables eh.dat, ehr1.dat, nwfits and atmc.dat
mpe integer /48/ # number of energy points
Expand Down Expand Up @@ -142,6 +172,23 @@ readehr1(fname:string) subroutine
readehr2(fname:string) subroutine
# read 2004 data (n=2-9) files thin.dat, thickLyA.dat, thickAllLy.dat
# in fname filename to read
readcrumpete(fname:string) subroutine
# read CRUMPET energy sink/source file
# in fname filename to read
readcrumpetn(fname:string) subroutine
# read CRUMPET rate file
# in fname filename to read
sv_crumpet(te:real,ne:real,rate:integer) real function [X m**3/s]
# calculate CRUMPET atomic creation form molecules
# in te=electron temperature [J]
# in ne=electron density [/m**3]
# in rate=data table to be read
# 10: mol diss rate, [X]=1
# 11: mol->atom rate, [X]=1
# 20: el E change due to mol interaction, [X]=J
# 21: i/a E change due to mol interaction, [X]=J
# 22: Epot change due to mol interaction, [X]=J
# 23: radiation due to mol interaction, [X]=J

***** Aphwrk:
# working arrays for 2-d spline interpolation
Expand All @@ -166,5 +213,6 @@ rqacoef(1:nxdata_aph,1:nydata_aph) _real # spline coeff's for line emission
***** Subs:
# Subroutines that can be called from the parser
aphread subroutine
crumpetread subroutine


152 changes: 152 additions & 0 deletions aph/aphrates.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1336,3 +1336,155 @@ real function svdiss (te)
end

c----------------------------------------------------------------------c
real function sv_crumpet (te, ne, rate)
implicit none
real te, ne
integer rate

Use(Dim)
Use(UEpar)
Use(Share) # istabon
Use(Physical_constants) # ev
Use(Data_input)
Use(Rtdata)
Use(Rtcrumpet)

c local variables --
integer je,jd
real zloge,zlogd,rle,rld,fje,fjd,mins,maxs,epsilon,
. svd_crm11,svd_crm12,svd_crm21,svd_crm22,svd_crm1,svd_crm2


c Compute molecular dissociation rate - A. Holm's CRUMPET
c te [J] = electron temperature
c ne [/m**3] = electron density
c rate = rate to access
c 10: mol. diss, X=[]
c 11: atom source from mol dis, X=[]
c 20: el. energy change from mol interactions, X=[J]
c 21: i/a energy change from mol interactions, X=[J]
c 22: Epot (binding E) change from mol interactions, X=[J]
c 23: Atom radiation source from mol interactionsm X=[J]
c 24: Mol. radiation source from mol interactionsm X=[J]
c svma_crumpet [X/s] = radiation rate

c----------------------------------------------------------------------c
c Logarithmic interpolation as in Stotler-95



if (ishymol .eq. 0) then
sv_crumpet=0
elseif (ismolcrm .eq. 0) then
sv_crumpet=0
else
c compute abscissae --
zloge=log(te/ev)
rle=max(crlemin, min(zloge,crlemax))
zlogd=log10(ne)
rld=max(crldmin, min(zlogd,crldmax))
c table indicies for interpolation --
je=int((rle-crlemin)/cdelekpt) + 1
je=min(je,cmpe-1)
jd=int((rld-crldmin)/cdeldkpt) + 1
jd=min(jd,cmpd-1)
c fractional parts of intervals (je,je+1) and (jd,jd+1) --
fje=(rle-cekpt(je))/(cekpt(je+1)-cekpt(je))
fjd=(rld-cdkpt(jd))/(cdkpt(jd+1)-cdkpt(jd))

c pick the appropriate matrix to interpolate from
if (rate .eq. 10) then
svd_crm11=crmdiss(je,jd)
svd_crm12=crmdiss(je,jd+1)
svd_crm21=crmdiss(je+1,jd)
svd_crm22=crmdiss(je+1,jd+1)
elseif (rate .eq. 11) then
svd_crm11=crmarate(je,jd)
svd_crm12=crmarate(je,jd+1)
svd_crm21=crmarate(je+1,jd)
svd_crm22=crmarate(je+1,jd+1)
elseif (rate .eq. 20) then
svd_crm11=crmselm(je,jd)
svd_crm12=crmselm(je,jd+1)
svd_crm21=crmselm(je+1,jd)
svd_crm22=crmselm(je+1,jd+1)
elseif (rate .eq. 21) then
svd_crm11=crmsiam(je,jd)
svd_crm12=crmsiam(je,jd+1)
svd_crm21=crmsiam(je+1,jd)
svd_crm22=crmsiam(je+1,jd+1)
elseif (rate .eq. 22) then
svd_crm11=crmspotm(je,jd)
svd_crm12=crmspotm(je,jd+1)
svd_crm21=crmspotm(je+1,jd)
svd_crm22=crmspotm(je+1,jd+1)
elseif (rate .eq. 23) then
svd_crm11=crmsrada(je,jd)
svd_crm12=crmsrada(je,jd+1)
svd_crm21=crmsrada(je+1,jd)
svd_crm22=crmsrada(je+1,jd+1)
elseif (rate .eq. 24) then
svd_crm11=crmsradm(je,jd)
svd_crm12=crmsradm(je,jd+1)
svd_crm21=crmsradm(je+1,jd)
svd_crm22=crmsradm(je+1,jd+1)
else
call remark("Rate option not recognized!")
endif


c ... Check whether we have all positive, all negative, or zero crossing
mins=min(svd_crm11,svd_crm12,svd_crm21,svd_crm22)
maxs=max(svd_crm11,svd_crm12,svd_crm21,svd_crm22)
if (mins .gt. 0) then
c ... Minimum is positive, normal interpolation
svd_crm11=log(svd_crm11)
svd_crm12=log(svd_crm12)
svd_crm21=log(svd_crm21)
svd_crm22=log(svd_crm22)
svd_crm1=svd_crm11+fjd*(svd_crm12-svd_crm11)
svd_crm2=svd_crm21+fjd*(svd_crm22-svd_crm21)
sv_crumpet = exp( svd_crm1 + fje*(svd_crm2-svd_crm1) )

elseif (maxs .lt. 0) then
c ... Maximum is negative, interpolate of negative values
svd_crm11=log(-svd_crm11)
svd_crm12=log(-svd_crm12)
svd_crm21=log(-svd_crm21)
svd_crm22=log(-svd_crm22)
svd_crm1=svd_crm11+fjd*(svd_crm12-svd_crm11)
svd_crm2=svd_crm21+fjd*(svd_crm22-svd_crm21)
sv_crumpet = -exp( svd_crm1 + fje*(svd_crm2-svd_crm1) )

else
c ... Zero-crossing – linear interpolation
jd=jd-1
je=je-1
fjd=((ne/1e6)-(10**(10+0.5*jd)))/ (10**(10+0.5*jd)*(10**0.5 - 1) )
fje=(te/ev - (10**(-1.2+0.1*je)))/( (10**(-1.2+0.1*je))*(10**0.1-1))

svd_crm1=svd_crm11+fjd*(svd_crm12-svd_crm11)
svd_crm2=svd_crm21+fjd*(svd_crm22-svd_crm21)
sv_crumpet = svd_crm1 + fje*(svd_crm2-svd_crm1)


endif

endif


c----------------------------------------------------------------------c
return
end












Loading
Loading