Monte Carlo module addition + PHO network

.gitignore

@@ -0,0 +1,2 @@
+GJ1214_b/*
+output/*

README.md

@@ -22,6 +22,7 @@ VULCAN requires the following python packages:
 - matplotlib
 - PIL/Pillow (optional: for interactive plotting)
 and the embeded [FastChem](https://github.com/exoclime/FastChem) requires a standard C++ compiler, like g++ or Clang.
+- numba (optional: for Monte Carlo RT module)
 
 If any of the python packages are missing, you can install the full SciPy Stack via Pip, e.g.
 ```bash
@@ -87,6 +88,7 @@ for a weaker vertical mixing (K<sub>zz</sub>) and carbon rich (C/O=1) run. Set u
 │   ├── op.py
 │   ├── phy_const.py
 │   ├── store.py
+│   ├── gCMCRT.py
 │   ├── vulcan.py
 │   ├── vulcan_cfg.py
 ```
@@ -104,7 +106,8 @@ for a weaker vertical mixing (K<sub>zz</sub>) and carbon rich (C/O=1) run. Set u
 `NCHO_photo_netowrk.txt`: the default N-C-H-O photochemical kinetics network  
 `op.py`: all the modules for the numerical operations, e.g. computing reaction rates, ODE solvers etc.    
 `make_chem_funs.py`: the routine that runs first to produce the required `chem_funs.py` based on the assigned chemical network    
-`phy_const.py`: physical constants  
+`phy_const.py`: physical constants 
+`gCMCRT.py`: 1D Monte Carlo radiative-transfer module 
 `store.py`: modules to store all the variables  
 `vulcan.py`: the top-level main script of VULCAN  
 `vulcan_cfg.py`: the configuration file for VULCAN  
@@ -157,6 +160,11 @@ will read vulcan output (.vul files) can plot the species profiles. Species shou
 The script of ```plot_vulcan.py``` should also serve as a good example of how to access to output files. The first step is to use "pickle.load" to unpack the binary files. The main variables are stored in three basic classes: data['variable'], data['atm'], and data['parameter'].
 You can also find all the names of variables and the class structure in ```store.py```.
 
+### Using the Monte Carlo radiative-transfer (MCRT) module ###
+To enable the MCRT module edit the use_gCMCRT variable to True in vulcan_cfg.py. The number of photon packets is given by the variable Nph in the same file.
+The MCRT will run slow at first past as it compiles and caches, but then will run a lot faster. We suggest lowering the number of packets if runtimes are too slow, but no smaller than 1000 packets. Alternativly, lower the number of wavelength bins for the photochemistry. 
+We include timing output of the MCRT module step to aid this determination.
+
 ### Troubleshooting ###
 - `xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun`
 

atm/atm_GJ1214b_Kzz.txt

@@ -0,0 +1,118 @@
+#(dyne/cm2) (K)          (cm2/s)
+ Pressure        Temp     Kzz
+896151000.0 1258.19 270418.62518778135
+719686000.0 1214.18 301755.83042297367
+577969000.0 1179.72 336724.6613278543
+464159000.0 1154.42 375745.5184963396
+372759000.0 1135.97 419288.6712632056
+299358000.0 1122.3 467877.07161151076
+240410000.0 1112.59 522096.68192223064
+193070000.0 1105.98 582599.0861696698
+155052000.0 1101.44 650112.4919445151
+124520000.0 1098.49 725450.2611340531
+100000000.0 1096.84 809519.2813847678
+80308600.0 1095.92 903329.4514514874
+64494700.0 1095.26 1008010.7968839619
+51794700.0 1094.79 1124823.9021855702
+41595600.0 1094.46 1255172.7462258216
+33404800.0 1094.23 1400627.8559314867
+26827000.0 1094.06 1562935.8132243967
+21544300.0 1093.94 1744058.320403745
+17302000.0 1093.85 1946162.2468942143
+13895000.0 1093.79 2171690.4513717936
+11158800.0 1093.75 2423362.8756560283
+8961510.0 1093.72 2704186.251877813
+7196860.0 1093.69 3017558.3042297373
+5779690.0 1093.68 3367246.613278543
+4641590.0 1093.66 3757455.1849633963
+3727590.0 1093.63 4192886.7126320554
+2993580.0 1093.43 4678770.716115107
+2404100.0 1092.62 5220966.819222306
+1930700.0 1090.38 5825990.861696697
+1550520.0 1085.72 6501124.919445151
+1245200.0 1077.63 7254502.611340532
+1000000.0 1065.27 8095192.813847679
+803086.0 1047.75 9033294.514514875
+644947.0 1023.91 10080107.968839617
+517947.0 993.02 11248239.021855703
+415956.0 956.145 12551727.462258214
+334048.0 915.988 14006278.559314867
+268270.0 875.857 15629358.132243967
+215443.0 838.701 17440583.20403745
+173020.0 805.452 19461622.468942147
+138950.0 775.858 21716904.513717934
+111588.0 750.461 24233628.756560277
+89615.1 729.133 27041862.518778134
+71968.6 710.102 30175583.04229737
+57796.9 693.34 33672466.13278543
+46415.9 678.918 37574551.849633954
+37275.9 666.26 41928867.12632056
+29935.8 655.391 46787707.161151074
+24041.0 646.555 52209668.19222306
+19307.0 639.154 58259908.61696697
+15505.2 632.522 65011249.1944515
+12452.0 627.217 72545026.1134053
+10000.0 623.304 80951928.13847679
+8030.86 619.473 90332945.14514875
+6449.47 616.11 100801079.68839617
+5179.47 613.423 112482390.21855702
+4159.56 611.153 125517274.62258212
+3340.48 609.111 140062785.59314868
+2682.7 607.555 156293581.3224397
+2154.43 606.438 174405832.04037452
+1730.2 605.179 194616224.68942145
+1389.5 604.089 217169045.13717937
+1115.88 603.265 242336287.56560278
+896.151 602.356 270418625.1877813
+719.686 601.221 301755830.4229737
+577.969 600.14 336724661.32785434
+464.159 599.055 375745518.49633956
+372.759 597.823 419288671.26320565
+299.358 596.504 467877071.6115107
+240.41 595.079 522096681.9222306
+193.07 593.436 582599086.1696696
+155.052 591.557 650112491.9445151
+124.52 589.537 725450261.1340531
+100.0 587.409 809519281.3847679
+80.3086 585.174 903329451.4514874
+64.4947 582.978 1008010796.8839619
+51.7947 580.848 1124823902.1855702
+41.5956 578.767 1255172746.2258215
+33.4048 576.686 1400627855.9314866
+26.827 574.507 1562935813.2243967
+21.5443 572.13 1744058320.403745
+17.302 569.465 1946162246.8942142
+13.895 566.5 2171690451.3717937
+11.1588 563.245 2423362875.656028
+8.96151 559.674 2704186251.8778133
+7.19686 555.864 3017558304.2297363
+5.77969 551.937 3367246613.278543
+4.64159 547.976 3757455184.9633965
+3.72759 544.164 4192886712.632056
+2.99358 540.619 4678770716.115108
+2.4041 537.335 5220966819.222305
+1.9307 534.306 5825990861.696697
+1.55052 531.511 6501124919.445151
+1.2452 528.951 7254502611.34053
+1.0 526.637 8095192813.847678
+0.803086 524.583 9033294514.514875
+0.644947 522.677 10080107968.839619
+0.517947 520.932 11248239021.855701
+0.415956 519.343 12551727462.258215
+0.334048 517.911 14006278559.314867
+0.26827 516.629 15629358132.243967
+0.215443 515.487 17440583204.037453
+0.17302 514.475 19461622468.942142
+0.13895 513.582 21716904513.717937
+0.111588 512.801 24233628756.56028
+0.0896151 512.115 27041862518.77813
+0.0719686 511.519 30175583042.29737
+0.0577969 511.0 33672466132.78543
+0.0464159 510.543 37574551849.63395
+0.0372759 510.145 41928867126.320564
+0.0299358 509.797 46787707161.15108
+0.024041 509.495 52209668192.22305
+0.019307 509.235 58259908616.966965
+0.0155052 509.011 65011249194.451515
+0.012452 508.818 72545026113.4053
+0.01 508.66 80951928138.47678

atm/atm_GJ1214b_Kzz_aer.txt

@@ -0,0 +1,118 @@
+#(dyne/cm2) (K)          (cm2/s)    (cm-3)   (um)   (cm2)
+Pressure        Temp     Kzz   nd   rm   sig2
+896151000.0 1258.19 270418.62518778135 0.0 0.0 0.0
+719686000.0 1214.18 301755.83042297367 0.0 0.0 0.0
+577969000.0 1179.72 336724.6613278543 0.0 0.0 0.0
+464159000.0 1154.42 375745.5184963396 0.0 0.0 0.0
+372759000.0 1135.97 419288.6712632056 0.0 0.0 0.0
+299358000.0 1122.3 467877.07161151076 0.0 0.0 0.0
+240410000.0 1112.59 522096.68192223064 0.0 0.0 0.0
+193070000.0 1105.98 582599.0861696698 0.0 0.0 0.0
+155052000.0 1101.44 650112.4919445151 0.0 0.0 0.0
+124520000.0 1098.49 725450.2611340531 0.0 0.0 0.0
+100000000.0 1096.84 809519.2813847678 0.0 0.0 0.0
+80308600.0 1095.92 903329.4514514874 0.0 0.0 0.0
+64494700.0 1095.26 1008010.7968839619 0.0 0.0 0.0
+51794700.0 1094.79 1124823.9021855702 0.0 0.0 0.0
+41595600.0 1094.46 1255172.7462258216 0.0 0.0 0.0
+33404800.0 1094.23 1400627.8559314867 0.0 0.0 0.0
+26827000.0 1094.06 1562935.8132243967 0.0 0.0 0.0
+21544300.0 1093.94 1744058.320403745 0.0 0.0 0.0
+17302000.0 1093.85 1946162.2468942143 0.0 0.0 0.0
+13895000.0 1093.79 2171690.4513717936 0.0 0.0 0.0
+11158800.0 1093.75 2423362.8756560283 0.0 0.0 0.0
+8961510.0 1093.72 2704186.251877813 0.0 0.0 0.0
+7196860.0 1093.69 3017558.3042297373 0.0 0.0 0.0
+5779690.0 1093.68 3367246.613278543 0.0 0.0 0.0
+4641590.0 1093.66 3757455.1849633963 0.0 0.0 0.0
+3727590.0 1093.63 4192886.7126320554 0.0 0.0 0.0
+2993580.0 1093.43 4678770.716115107 0.0 0.0 0.0
+2404100.0 1092.62 5220966.819222306 0.0 0.0 0.0
+1930700.0 1090.38 5825990.861696697 0.0 0.0 0.0
+1550520.0 1085.72 6501124.919445151 0.0 0.0 0.0
+1245200.0 1077.63 7254502.611340532 0.0 0.0 0.0
+1000000.0 1065.27 8095192.813847679 0.0 0.0 0.0
+803086.0 1047.75 9033294.514514875 0.0 0.0 0.0
+644947.0 1023.91 10080107.968839617 0.0 0.0 0.0
+517947.0 993.02 11248239.021855703 0.0 0.0 0.0
+415956.0 956.145 12551727.462258214 0.0 0.0 0.0
+334048.0 915.988 14006278.559314867 0.0 0.0 0.0
+268270.0 875.857 15629358.132243967 0.0 0.0 0.0
+215443.0 838.701 17440583.20403745 0.0 0.0 0.0
+173020.0 805.452 19461622.468942147 0.0 0.0 0.0
+138950.0 775.858 21716904.513717934 0.0 0.0 0.0
+111588.0 750.461 24233628.756560277 0.0 0.0 0.0
+89615.1 729.133 27041862.518778134 1000000.0 1.0 1.0
+71968.6 710.102 30175583.04229737 1000000.0 1.0 1.0
+57796.9 693.34 33672466.13278543 1000000.0 1.0 1.0
+46415.9 678.918 37574551.849633954 1000000.0 1.0 1.0
+37275.9 666.26 41928867.12632056 1000000.0 1.0 1.0
+29935.8 655.391 46787707.161151074 1000000.0 1.0 1.0
+24041.0 646.555 52209668.19222306 1000000.0 1.0 1.0
+19307.0 639.154 58259908.61696697 1000000.0 1.0 1.0
+15505.2 632.522 65011249.1944515 1000000.0 1.0 1.0
+12452.0 627.217 72545026.1134053 1000000.0 1.0 1.0
+10000.0 623.304 80951928.13847679 1000000.0 1.0 1.0
+8030.86 619.473 90332945.14514875 1000000.0 1.0 1.0
+6449.47 616.11 100801079.68839617 1000000.0 1.0 1.0
+5179.47 613.423 112482390.21855702 1000000.0 1.0 1.0
+4159.56 611.153 125517274.62258212 1000000.0 1.0 1.0
+3340.48 609.111 140062785.59314868 1000000.0 1.0 1.0
+2682.7 607.555 156293581.3224397 1000000.0 1.0 1.0
+2154.43 606.438 174405832.04037452 1000000.0 1.0 1.0
+1730.2 605.179 194616224.68942145 1000000.0 1.0 1.0
+1389.5 604.089 217169045.13717937 1000000.0 1.0 1.0
+1115.88 603.265 242336287.56560278 1000000.0 1.0 1.0
+896.151 602.356 270418625.1877813 1000000.0 1.0 1.0
+719.686 601.221 301755830.4229737 1000000.0 1.0 1.0
+577.969 600.14 336724661.32785434 1000000.0 1.0 1.0
+464.159 599.055 375745518.49633956 1000000.0 1.0 1.0
+372.759 597.823 419288671.26320565 1000000.0 1.0 1.0
+299.358 596.504 467877071.6115107 1000000.0 1.0 1.0
+240.41 595.079 522096681.9222306 1000000.0 1.0 1.0
+193.07 593.436 582599086.1696696 1000000.0 1.0 1.0
+155.052 591.557 650112491.9445151 1000000.0 1.0 1.0
+124.52 589.537 725450261.1340531 1000000.0 1.0 1.0
+100.0 587.409 809519281.3847679 1000000.0 1.0 1.0
+80.3086 585.174 903329451.4514874 1000000.0 1.0 1.0
+64.4947 582.978 1008010796.8839619 1000000.0 1.0 1.0
+51.7947 580.848 1124823902.1855702 1000000.0 1.0 1.0
+41.5956 578.767 1255172746.2258215 1000000.0 1.0 1.0
+33.4048 576.686 1400627855.9314866 1000000.0 1.0 1.0
+26.827 574.507 1562935813.2243967 1000000.0 1.0 1.0
+21.5443 572.13 1744058320.403745 1000000.0 1.0 1.0
+17.302 569.465 1946162246.8942142 1000000.0 1.0 1.0
+13.895 566.5 2171690451.3717937 1000000.0 1.0 1.0
+11.1588 563.245 2423362875.656028 1000000.0 1.0 1.0
+8.96151 559.674 2704186251.8778133 1000000.0 1.0 1.0
+7.19686 555.864 3017558304.2297363 1000000.0 1.0 1.0
+5.77969 551.937 3367246613.278543 1000000.0 1.0 1.0
+4.64159 547.976 3757455184.9633965 1000000.0 1.0 1.0
+3.72759 544.164 4192886712.632056 1000000.0 1.0 1.0
+2.99358 540.619 4678770716.115108 1000000.0 1.0 1.0
+2.4041 537.335 5220966819.222305 1000000.0 1.0 1.0
+1.9307 534.306 5825990861.696697 1000000.0 1.0 1.0
+1.55052 531.511 6501124919.445151 1000000.0 1.0 1.0
+1.2452 528.951 7254502611.34053 1000000.0 1.0 1.0
+1.0 526.637 8095192813.847678 1000000.0 1.0 1.0
+0.803086 524.583 9033294514.514875 1000000.0 1.0 1.0
+0.644947 522.677 10080107968.839619 1000000.0 1.0 1.0
+0.517947 520.932 11248239021.855701 1000000.0 1.0 1.0
+0.415956 519.343 12551727462.258215 1000000.0 1.0 1.0
+0.334048 517.911 14006278559.314867 1000000.0 1.0 1.0
+0.26827 516.629 15629358132.243967 1000000.0 1.0 1.0
+0.215443 515.487 17440583204.037453 1000000.0 1.0 1.0
+0.17302 514.475 19461622468.942142 1000000.0 1.0 1.0
+0.13895 513.582 21716904513.717937 1000000.0 1.0 1.0
+0.111588 512.801 24233628756.56028 1000000.0 1.0 1.0
+0.0896151 512.115 27041862518.77813 1000000.0 1.0 1.0
+0.0719686 511.519 30175583042.29737 1000000.0 1.0 1.0
+0.0577969 511.0 33672466132.78543 1000000.0 1.0 1.0
+0.0464159 510.543 37574551849.63395 1000000.0 1.0 1.0
+0.0372759 510.145 41928867126.320564 1000000.0 1.0 1.0
+0.0299358 509.797 46787707161.15108 1000000.0 1.0 1.0
+0.024041 509.495 52209668192.22305 1000000.0 1.0 1.0
+0.019307 509.235 58259908616.966965 1000000.0 1.0 1.0
+0.0155052 509.011 65011249194.451515 1000000.0 1.0 1.0
+0.012452 508.818 72545026113.4053 1000000.0 1.0 1.0
+0.01 508.66 80951928138.47678 1000000.0 1.0 1.0

atm/make_atm_Kzz_cld.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+line_1 = '#(dyne/cm2) (K)          (cm2/s)    (cm-3)   (um)   (cm2)'
+line_2 = 'Pressure        Temp     Kzz   nd   rm   sig2'
+
+fname = 'atm_GJ1214b_Kzz.txt'
+
+
+data = np.loadtxt(fname,skiprows=2)
+
+p = data[:,0]
+T = data[:,1]
+Kzz = data[:,2]
+
+nz = len(p)
+
+# Set up aerosol component
+
+nd = np.zeros(nz)
+rm = np.zeros(nz)
+sig2 = np.zeros(nz)
+
+
+for i in range(nz):
+  if (p[i]/1e6 < 1e-1):
+    nd[i] = 1e6
+    rm[i] = 1e0
+    sig2[i] = 1.0
+
+fout = 'atm_GJ1214b_Kzz_aer.txt'
+f = open(fout, 'w')
+
+f.write(line_1 + '\n')
+f.write(line_2 + '\n')
+for i in range(nz):
+  f.write(str(p[i]) + ' ' + str(T[i]) + ' ' + str(Kzz[i]) + ' ' + str(nd[i]) + ' ' + str(rm[i]) + ' ' + str(sig2[i]) + '\n')
+
+
+

build_atm.py

@@ -337,6 +337,9 @@ def load_TPK(self, data_atm):
             if self.Kzz_prof == 'file':
                 atm_table = np.genfromtxt(vulcan_cfg.atm_file, names=True, dtype=None, skip_header=1)
                 p_file, T_file, Kzz_file = atm_table['Pressure'], atm_table['Temp'], atm_table['Kzz']
+
+                if (vulcan_cfg.use_aer == True):
+                  nd_file, rm_file, sig2_file = atm_table['nd'], atm_table['rm'], atm_table['sig2']
 
             else:     
                 atm_table = np.genfromtxt(vulcan_cfg.atm_file, names=True, dtype=None, skip_header=1)
@@ -370,6 +373,32 @@ def load_TPK(self, data_atm):
 
             elif self.Kzz_prof == 'const': data_atm.Kzz = np.repeat(self.const_Kzz,nz-1)
 
+            if (vulcan_cfg.use_aer == True):
+                PTK_fun['pnd'] = interpolate.interp1d(p_file, nd_file, assume_sorted = False, bounds_error=False, fill_value=(nd_file[np.argmin(p_file)], nd_file[np.argmax(p_file)]) )
+                # store aerosol data in data_atm
+                try:
+                    data_atm.nd = PTK_fun['pnd'](data_atm.pco)
+                # for SciPy earlier than v0.18.0
+                except ValueError:
+                    PTK_fun['pnd'] = interpolate.interp1d(p_file, nd_file, assume_sorted = False, bounds_error=False, fill_value=nd_file[np.argmin(p_file)] )
+                    data_atm.nd = PTK_fun['pnd'](data_atm.pco)
+                PTK_fun['prm'] = interpolate.interp1d(p_file, rm_file, assume_sorted = False, bounds_error=False, fill_value=(rm_file[np.argmin(p_file)], rm_file[np.argmax(p_file)]) )
+                # store aerosol data in data_atm
+                try:
+                    data_atm.rm = PTK_fun['prm'](data_atm.pco)
+                # for SciPy earlier than v0.18.0
+                except ValueError:
+                    PTK_fun['prm'] = interpolate.interp1d(p_file, rm_file, assume_sorted = False, bounds_error=False, fill_value=rm_file[np.argmin(p_file)] )
+                    data_atm.rm = PTK_fun['prm'](data_atm.pco)
+                PTK_fun['psig2'] = interpolate.interp1d(p_file, sig2_file, assume_sorted = False, bounds_error=False, fill_value=(sig2_file[np.argmin(p_file)], sig2_file[np.argmax(p_file)]) )
+                # store aerosol data in data_atm
+                try:
+                    data_atm.sig2 = PTK_fun['psig2'](data_atm.pco)
+                # for SciPy earlier than v0.18.0
+                except ValueError:
+                    PTK_fun['psig2'] = interpolate.interp1d(p_file, sig2_file, assume_sorted = False, bounds_error=False, fill_value=sig2_file[np.argmin(p_file)] )
+                    data_atm.sig2 = PTK_fun['psig2'](data_atm.pco)
+
         elif self.type == 'vulcan_ini':
             print ("Initializing PT from the prvious run " + vulcan_cfg.vul_ini)
             with open(vulcan_cfg.vul_ini, 'rb') as handle: