Add optional emission opacity array for eventual coronal emissivity

GAS/gas_update.f

@@ -156,8 +156,17 @@ subroutine gas_update(it)
 c-- initialize opacity
         gas_sig = 0d0
         gas_cap = 0.0
+        if(in_doemiss) then
+          if(.not.allocated(gas_em_cap)) then
+            stop 'in_doemiss and gas_em_cap not allocated'
+          endif
+          if(.not.allocated(gas_em_capgrey)) then
+            stop 'in_doemiss and gas_em_capgrey not allocated'
+          endif
+          gas_em_cap = 0.0
+        endif
         if(.not.in_notbopac) then
-          call tabular_opacity
+          call tabular_opacity(in_doemiss)
         else
           call physical_opacity(dtempfrac*temp)
         endif
@@ -233,9 +242,18 @@ subroutine gas_update(it)
 c-- initialize opacity
         gas_sig = 0d0
         gas_cap = 0.0
+        if(in_doemiss) then
+          if(.not.allocated(gas_em_cap)) then
+            stop 'in_doemiss and gas_em_cap not allocated'
+          endif
+          if(.not.allocated(gas_em_capgrey)) then
+            stop 'in_doemiss and gas_em_capgrey not allocated'
+          endif
+          gas_em_cap = 0.0
+        endif
 c-- calculate opacities
         if(.not.in_notbopac) then
-          call tabular_opacity
+          call tabular_opacity(in_doemiss)
         else
           call physical_opacity(temp)
         endif

GAS/opacity_planckmean.f

@@ -30,6 +30,13 @@ subroutine opacity_planckmean
        call specintv(1/gas_temp(i),grp_ng,specarr)
        gas_capgrey(i) = sum(gas_cap(:,i)*specarr)
       enddo !i
+c-- compute Planck weighted integral of emission opacity
+      if (in_doemiss) then
+        do i=1,gas_ncell
+           call specintv(1/gas_temp(i),grp_ng,specarr)
+           gas_em_capgrey(i) = sum(gas_em_cap(:,i)*specarr)
+        enddo !i
+      endif
 c
 c
 cc-- old

GAS/tabular_opacity.f

@@ -6,7 +6,7 @@
 * The Government is granted for itself and others acting on its behalf a nonexclusive, paid-up,
 * irrevocable worldwide license in this material to reproduce, prepare. derivative works, distribute
 * copies to the public, perform publicly and display publicly, and to permit others to do so.
-      subroutine tabular_opacity
+      subroutine tabular_opacity(lemiss)
 c     --------------------------
       use tbxsmod
       use miscmod, only:binsrch
@@ -15,10 +15,11 @@ subroutine tabular_opacity
       use groupmod
       use physconstmod
       implicit none
+      logical, intent(in) :: lemiss
 ************************************************************************
 * Interpolate Fontes opacity table
 ************************************************************************
-      integer :: itemp,irho,i,iz,ig,j,l
+      integer :: itemp,irho,i,iz,ig,j,l,ll
       real*8 :: massfr, rhopart, temph, rhoh
       real*8 :: phi1t,phi2t,phi1r,phi2r
       real*8 :: sig1,sig2,sig3,sig4
@@ -149,6 +150,103 @@ subroutine tabular_opacity
        if(iand(j,8)==0) write(0,*) 'opacity_calc: all cap==inf'
       endif
 c
+c
+c
+c--interpolate emission opacity
+      if (lemiss) then
+c
+c--   calculate emission opacities
+        do i=1,gas_ncell
+          if(gas_mass(i)<=0d0) cycle
+c--   table abundance for cell
+          ll = 1
+          do l=1,tb_nelem
+            iz=tb_ielem(l)
+            if(gas_natom1fr(iz,i)<=0d0) cycle
+c--   mass fraction helper
+            massfr=gas_natom1fr(iz,i)*gas_natom(i) *
+     &           elem_data(iz)%m*pc_amu/gas_mass(i)
+c--   partial density for table interpolation
+            rhopart=massfr*gas_rho(i)
+c--   search 1d tb arrays for temp-rho point
+            irho=binsrch(rhopart,tb_rho,tb_nrho,.false.)
+            itemp=binsrch(gas_temp(i),tb_temp,tb_ntemp,.false.)
+c--   do not allow for extrapolation
+            rhoh=max(rhopart,tb_rho(irho))
+            rhoh=min(rhoh,tb_rho(irho+1))
+            temph=max(gas_temp(i),tb_temp(itemp))
+            temph=min(temph,tb_temp(itemp+1))
+c--   2d bilinear interpolation in temp-rho
+c--   temperature basis functions
+            phi1t=(tb_temp(itemp+1)-temph) /
+     &           (tb_temp(itemp+1)-tb_temp(itemp))
+            phi2t=(temph-tb_temp(itemp)) /
+     &           (tb_temp(itemp+1)-tb_temp(itemp))
+c--   density basis functions
+            phi1r=log(tb_rho(irho+1)/rhoh) /
+     &           log(tb_rho(irho+1)/tb_rho(irho))
+            phi2r=log(rhoh/tb_rho(irho)) /
+     &           log(tb_rho(irho+1)/tb_rho(irho))
+c--   emission
+            if (tb_ielem_em_mask(l,irho)) then
+              cap1=tb_em_cap(:,itemp,irho,ll)*sngl(phi1t*phi1r)
+              cap4=tb_em_cap(:,itemp+1,irho,ll)*sngl(phi2t*phi1r)
+            else
+              cap1=tb_cap(:,itemp,irho,l)*sngl(phi1t*phi1r)
+              cap4=tb_cap(:,itemp+1,irho,l)*sngl(phi2t*phi1r)
+            endif
+            if (tb_ielem_em_mask(l,irho+1)) then
+              cap2=tb_em_cap(:,itemp,irho+1,ll)*sngl(phi1t*phi2r)
+              cap3=tb_em_cap(:,itemp+1,irho+1,ll)*sngl(phi2t*phi2r)
+            else
+              cap2=tb_cap(:,itemp,irho+1,l)*sngl(phi1t*phi2r)
+              cap3=tb_cap(:,itemp+1,irho+1,l)*sngl(phi2t*phi2r)
+            endif
+            cap(:,i)=cap1+cap2+cap3+cap4
+c--   add macroscopic emission opacity to total
+            gas_em_cap(:,i)=gas_em_cap(:,i)+sngl(rhopart)*cap(:,i)
+c--   increment emission table sub-counter
+            if (tb_ielem_em_mask(l,irho).or.tb_ielem_em_mask(l,irho+1))
+     &           ll = ll + 1
+c
+          enddo
+        enddo
+c
+c--   sanity check (duplicated from physical_opacity)
+c--   emission
+        j = 0
+        do i=1,gas_ncell
+          if(gas_mass(i)<=0d0) cycle
+          do ig=1,grp_ng
+            if(gas_em_cap(ig,i)==0.) j = ior(j,1)
+            if(gas_em_cap(ig,i)<0.) j = ior(j,2)
+            if(gas_em_cap(ig,i)/=gas_em_cap(ig,i)) j = ior(j,4)
+            if(gas_em_cap(ig,i)>huge(gas_em_cap)) j = ior(j,8)
+          enddo                 !ig
+        enddo                   !i
+        if(iand(j,1)/=0) write(0,*) 'opacity_calc: some em_cap==0'
+        if(iand(j,2)/=0) write(0,*) 'opacity_calc: some em_cap<0'
+        if(iand(j,4)/=0) write(0,*) 'opacity_calc: some em_cap==NaN'
+        if(iand(j,8)/=0) write(0,*) 'opacity_calc: some em_cap==inf'
+c--   refine sanity check
+        if(j>0) then
+          j = 0
+          do i=1,gas_ncell
+            if(gas_mass(i)<=0d0) cycle
+            do ig=1,grp_ng
+              if(gas_em_cap(ig,i)/=0.) j = ior(j,1)
+              if(gas_em_cap(ig,i)>=0.) j = ior(j,2)
+              if(gas_em_cap(ig,i)==gas_em_cap(ig,i)) j = ior(j,4)
+              if(gas_em_cap(ig,i)<=huge(gas_em_cap)) j = ior(j,8)
+            enddo               !ig
+          enddo                 !i
+          if(iand(j,1)==0) write(0,*) 'opacity_calc: all em_cap==0'
+          if(iand(j,2)==0) write(0,*) 'opacity_calc: all em_cap<0'
+          if(iand(j,4)==0) write(0,*) 'opacity_calc: all em_cap==NaN'
+          if(iand(j,8)==0) write(0,*) 'opacity_calc: all em_cap==inf'
+        endif
+      endif
+c
 c-- remove opacity helpers from heap
       deallocate(sig)
       deallocate(cap,cap1,cap2,cap3,cap4)

Input/input.emiss.par

@@ -0,0 +1,75 @@
+&inputpars
+
+!-- misc
+ in_nomp = 16
+ in_io_grabstdout = f
+ in_puretran = f
+ in_noreadstruct = f
+ in_noeos = t
+ in_flx_noobservertime = f
+ in_trn_nolumpshortcut = t
+ in_doemiss = t
+
+!-- time step
+ in_tsp_tfirst = 4841952 ! 56 days
+ in_tsp_tlast = 4841952.25
+ in_tsp_gridtype = 'expo'
+ in_tsp_nt = 1
+ in_tsp_itrestart = 1
+
+!-- grid settings
+ in_grd_igeom = 11
+ in_ndim = 1, 1, 1
+
+!-- flat structures
+ in_gas_gastempinit = 0d0
+ in_gas_cvcoef = -1d0
+
+!-- flux settings
+ in_flx_ndim = 1024, 1, 1
+ in_flx_wlmin = 1d-5
+ in_flx_wlmax = 1.28d-3
+
+!-- group settings
+ in_grp_ng = 1024
+ in_grp_wlmin = 1d-5
+ in_grp_wlmax = 1.28d-3
+
+!-- source
+ in_novolsrc=t
+! in_srctype="tabl"
+! in_opcapgam=0.1d0
+! in_gas_srccoef=1.3d16
+! in_gas_srcrpwr=1d0
+! in_gas_srctimepwr=-1.3d0
+
+!-- opacity
+ in_opacanaltype="none"
+ in_gas_capcoef = 100
+ in_gas_caprpwr = 1d0
+
+!-- particle numbers
+ in_src_n2s  = 10
+ in_prt_n2max = 16
+ in_srcepwr = 0.5d0
+
+!-- imc tally method
+ in_trn_isimcanlog = f
+ in_trn_isddmcanlog = f
+!-- ddmc mfp per cell
+ in_trn_tauddmc = 10
+ in_taulump = 10
+
+!-- opacity debug switches
+ in_nobbopac = t
+ in_nobfopac = t
+ in_noffopac = t
+ in_nothmson = t
+ in_noplanckweighting = t
+ in_notbopac = f
+ in_notbbbopac = f
+ in_notbbfopac = f
+ in_notbffopac = f
+ in_notbthmson = f
+
+/

Input/input.str_emiss

@@ -0,0 +1,4 @@
+# spherical
+# 1 1 1 5 1
+#  rightvel        mass        temp          ye          nd
+ 7.4925e+09  2.0000e+30  3.0000e+04  3.9000e-01  1.0000e+00

gasmod.f

@@ -57,6 +57,7 @@ module gasmod
 c-- Line+Cont extinction coeff
       real*8,allocatable :: gas_capcoef(:) !(ncell)
       real*4,allocatable :: gas_cap(:,:) !(ng,ncell)
+      real*4,allocatable :: gas_em_cap(:,:) !(ng,ncell) -- emission opacity
 c-- leakage opacities
 c     real*8,allocatable :: dd_opacleak(:,:) !(6,ncell)
 c-- scattering coefficient
@@ -65,6 +66,7 @@ module gasmod
       real*8,allocatable :: gas_capgam(:) !(ncell)
 c-- Planck opacity (gray)
       real*8,allocatable :: gas_capgrey(:)!(ncell)
+      real*8,allocatable :: gas_em_capgrey(:)!(ncell) -- Planck-integrated emission opacity
 c-- Fleck factor
       real*8,allocatable :: gas_fcoef(:)  !(ncell)
 c
@@ -81,11 +83,12 @@ module gasmod
       contains
 c
 c
-      subroutine gasmod_init(ltalk,icell1,ncell,ngin)
+      subroutine gasmod_init(ltalk,icell1,ncell,ngin,lemiss)
 c----------------------------------------!{{{
       implicit none
       logical,intent(in) :: ltalk
       integer,intent(in) :: icell1,ncell,ngin
+      logical,intent(in) :: lemiss
 ************************************************************************
 * Allocate gas variables.
 *
@@ -149,6 +152,11 @@ subroutine gasmod_init(ltalk,icell1,ncell,ngin)
 c
 c-- ndim=2 alloc big
       allocate(gas_cap(ng,gas_ncell))
+c-- optional emissivity arrays
+      if (lemiss) then
+        allocate(gas_em_cap(ng,gas_ncell))
+        allocate(gas_em_capgrey(gas_ncell))
+      endif
 !}}}
       end subroutine gasmod_init
 c
@@ -180,7 +188,10 @@ subroutine gas_dealloc
       deallocate(gas_evolinit)
       deallocate(gas_natom1fr)
       deallocate(gas_natom0fr)
-      deallocate(gas_cap)!}}}
+      deallocate(gas_cap)
+      if (allocated(gas_em_cap)) deallocate(gas_em_cap)
+      if (allocated(gas_em_cap)) deallocate(gas_em_capgrey)
+!}}}
       end subroutine gas_dealloc
 c
       end module gasmod

gridmod.f

@@ -52,6 +52,7 @@ module gridmod
 
 c-- Line+Cont extinction coeff
       real*4,allocatable :: grd_cap(:,:) !(ng,ncell)
+      real*4,allocatable :: grd_em_cap(:,:) !(ng,ncell)
 
 c-- leakage opacities
       real*8,allocatable :: grd_opaclump(:,:) !(10,ncell) leak(6),speclump,caplump,igemitmax,doplump
@@ -60,6 +61,7 @@ module gridmod
       real*8,allocatable :: grd_sig(:) !(ncell) !grey scattering opacity
 c-- Planck opacity (gray)
       real*8,allocatable :: grd_capgrey(:) !(ncell)
+      real*8,allocatable :: grd_em_capgrey(:) !(ncell)
 c-- Fleck factor
       real*8,allocatable :: grd_fcoef(:)  !(ncell)
 
@@ -105,10 +107,10 @@ end function emitgroup
 c
       contains
 c
-      subroutine gridmod_init(ltalk,ngin,ncell,lvoid,idd1,ndd)
+      subroutine gridmod_init(ltalk,ngin,ncell,lvoid,idd1,ndd,lemiss)
 c     --------------------------------------------------!{{{
       implicit none
-      logical,intent(in) :: ltalk,lvoid
+      logical,intent(in) :: ltalk,lvoid,lemiss
       integer,intent(in) :: ngin
       integer,intent(in) :: ncell,idd1,ndd
 ************************************************************************
@@ -147,6 +149,9 @@ subroutine gridmod_init(ltalk,ngin,ncell,lvoid,idd1,ndd)
        n = int((int(grd_ncell,8)*(8*(12+6) + 5*4))/1024) !kB
        write(6,*) 'ALLOC grd      :',n,"kB",n/1024,"MB",n/1024**2,"GB"
        n = int((int(grd_ncell,8)*4*ng)/1024) !kB
+       if (lemiss) then
+         n = 2 * n
+       endif
        write(6,*) 'ALLOC grd_cap  :',n,"kB",n/1024,"MB",n/1024**2,"GB"
       endif
 c
@@ -180,6 +185,12 @@ subroutine gridmod_init(ltalk,ngin,ncell,lvoid,idd1,ndd)
       allocate(grd_emitprob(grd_nep,grd_ncell))
 c-- ndim=4 alloc
       allocate(grd_cap(ng,grd_ncell))
+c
+c-- optional emissivity arrays
+      if (lemiss) then
+        allocate(grd_em_cap(ng,grd_ncell))
+        allocate(grd_em_capgrey(grd_ncell))
+      endif
 c!}}}
       end subroutine gridmod_init
 c
@@ -215,7 +226,11 @@ subroutine grid_dealloc
       deallocate(grd_opaclump)
       deallocate(grd_emitprob)
 c-- ndim=4 alloc
-      deallocate(grd_cap)!}}}
+      deallocate(grd_cap)
+c-- optional arrays
+      if (allocated(grd_em_cap)) deallocate(grd_em_cap)
+      if (allocated(grd_em_cap)) deallocate(grd_em_capgrey)
+!}}}
       end subroutine grid_dealloc
 c
       end module gridmod