Refactor emissivity IMC-DDMC interface/boundary condition.

+ Add pure function with emissivity formula to transportmod.
+ Use in 1D/3D spherical, 3D cylindrical, and 3D planar geometry.

GRID/leakage_opacity1.f90

@@ -141,8 +141,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(1,l)=grd_opaclump(1,l)+(specval*speclump)*&
                 1.5d0*pp*grd_xarr(i)**2/(dx3(i)*thelp)
         else
@@ -165,8 +165,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(2,l)=grd_opaclump(2,l)+(specval*speclump)*&
                 1.5d0*pp*grd_xarr(i+1)**2/(dx3(i)*thelp)
         else
@@ -189,8 +189,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i)*dyac(j)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(xm(i)*dyac(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(3,l)=grd_opaclump(3,l)+(specval*speclump)*&
                 0.75d0*pp*dx2(i)*sqrt(1d0-grd_yarr(j)**2)/(dy(j)*dx3(i)*thelp)
         else
@@ -214,8 +214,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i)*dyac(j)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(xm(i)*dyac(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(4,l)=grd_opaclump(4,l)+(specval*speclump)*&
                 0.75d0*pp*dx2(i)*sqrt(1d0-grd_yarr(j+1)**2)/(dy(j)*dx3(i)*thelp)
         else
@@ -242,9 +242,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i)*ym(j) * &
-              dz(k)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(xm(i)*ym(j)*dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(5,l)=grd_opaclump(5,l)+(specval*speclump)*&
                 0.75d0*pp*dx2(i)*dyac(j)/(dy(j)*dx3(i)*dz(k)*thelp)
         else
@@ -269,9 +268,8 @@ subroutine leakage_opacity1
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i)*ym(j) * &
-              dz(k)*thelp
-           pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(xm(i)*ym(j)*dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(6,l)=grd_opaclump(6,l)+(specval*speclump)*&
                 0.75d0*pp*dx2(i)*dyac(j)/(dy(j)*dx3(i)*dz(k)*thelp)
         else

GRID/leakage_opacity11.f90

@@ -120,8 +120,8 @@ subroutine leakage_opacity11
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           ppl = 4d0/(3d0*help+6d0*pc_dext)
+           ppl = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(1,l)=grd_opaclump(1,l)+(specval*speclump)*&
                 1.5d0*ppl*(thelp*grd_xarr(i))**2/ &
                 (3d0*grd_vol(l)/pc_pi4)
@@ -145,8 +145,8 @@ subroutine leakage_opacity11
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           ppr = 4d0/(3d0*help+6d0*pc_dext)
+           ppr = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(2,l)=grd_opaclump(2,l)+(specval*speclump)*&
                 1.5d0*ppr*(thelp*grd_xarr(i+1))**2/ &
                 (3d0*grd_vol(l)/pc_pi4)

GRID/leakage_opacity2.f90

@@ -138,8 +138,8 @@ subroutine leakage_opacity2
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           ppl = 4d0/(3d0*help+6d0*pc_dext)
+           ppl = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(1,l)=grd_opaclump(1,l)+(specval*speclump)*&
                 ppl*(thelp*grd_xarr(i))/(dx2(i)*thelp**2)
         else
@@ -162,8 +162,8 @@ subroutine leakage_opacity2
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           ppr = 4d0/(3d0*help+6d0*pc_dext)
+           ppr = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(2,l)=grd_opaclump(2,l)+(specval*speclump)*&
                 ppr*(thelp*grd_xarr(i+1))/(dx2(i)*thelp**2)
         else
@@ -186,8 +186,8 @@ subroutine leakage_opacity2
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dy(j)*thelp
-           ppl = 4d0/(3d0*help+6d0*pc_dext)
+           ppl = ddmc_emiss_bc(dy(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(3,l)=grd_opaclump(3,l)+(specval*speclump)*&
                 0.5d0*ppl/(thelp*dy(j))
         else
@@ -210,8 +210,8 @@ subroutine leakage_opacity2
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dy(j)*thelp
-           ppr = 4d0/(3d0*help+6d0*pc_dext)
+           ppr = ddmc_emiss_bc(dy(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(4,l)=grd_opaclump(4,l)+(specval*speclump)*&
                 0.5d0*ppr/(thelp*dy(j))
         else
@@ -237,9 +237,8 @@ subroutine leakage_opacity2
 
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i) * &
-              dz(k)*thelp
-           ppl = 4d0/(3d0*help+6d0*pc_dext)
+           ppl = ddmc_emiss_bc(xm(i)*dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(5,l)=grd_opaclump(5,l)+(specval*speclump)*&
               0.5d0*ppl/(xm(i)*dz(k)*thelp)
         else
@@ -263,9 +262,8 @@ subroutine leakage_opacity2
 
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*xm(i) * &
-              dz(k)*thelp
-           ppr = 4d0/(3d0*help+6d0*pc_dext)
+           ppr = ddmc_emiss_bc(xm(i)*dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(6,l)=grd_opaclump(6,l)+(specval*speclump)*&
               0.5d0*ppr/(xm(i)*dz(k)*thelp)
         else

GRID/leakage_opacity3.f90

@@ -127,14 +127,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(1,l) = grd_opaclump(1,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dx(i))
         else
@@ -157,14 +151,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dx(i)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dx(i)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(2,l) = grd_opaclump(2,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dx(i))
         else
@@ -187,14 +175,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dy(j)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dy(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(3,l) = grd_opaclump(3,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dy(j))
         else
@@ -217,14 +199,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dy(j)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dy(j)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(4,l) = grd_opaclump(4,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dy(j))
         else
@@ -247,14 +223,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dz(k)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(5,l) = grd_opaclump(5,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dz(k))
         else
@@ -277,14 +247,8 @@ subroutine leakage_opacity3
 !
         if(lhelp) then
 !-- DDMC interface
-           help = (grd_cap(ig,l)+grd_sig(l))*dz(k)*thelp
-           alb = grd_fcoef(l)*grd_cap(ig,l)/ &
-                (grd_cap(ig,l)+grd_sig(l))
-           eps = (4d0/3d0)*sqrt(3d0*alb)/(1d0+pc_dext*sqrt(3d0*alb))
-           beta = 1.5d0*alb*help**2+sqrt(3d0*alb*help**2 + &
-                2.25d0*alb**2*help**4)
-           pp = 0.5d0*eps*beta/(beta-0.75*eps*help)
-!              pp = 4d0/(3d0*help+6d0*pc_dext)
+           pp = ddmc_emiss_bc(dz(k)*thelp, grd_fcoef(l), &
+                grd_cap(ig,l), grd_sig(l), pc_dext)
            grd_opaclump(6,l) = grd_opaclump(6,l)+(specval*speclump)*&
                 0.5d0*pp/(thelp*dz(k))
         else

SOURCE/interior_source.f90

@@ -107,11 +107,11 @@ subroutine interior_source
 !-- calculating direction cosine (comoving)
      call rnd_r(r1,rnd_state)
      mu0 = 1d0-2d0*r1
-     ptcl%mu = mu0 !overwrite when isvelocity
+     ptcl%mu = mu0
 !-- sampling azimuthal angle of direction
      call rnd_r(r1,rnd_state)
      om0 = pc_pi2*r1
-     ptcl%om = om0 !overwrite when isvelocity
+     ptcl%om = om0
 
 !
 !-- x position
@@ -146,7 +146,7 @@ subroutine interior_source
   enddo !j
   enddo !k
   if(ipart/=src_nsurf+src_nnonth) stop 'interior_source: n/=nexecsrc'
-  
+
 
 !-- Thermal volume particle instantiation: loop
   iimpi = -1