diff --git a/basalt_buff_tunespin_bisec_v2.py b/basalt_buff_tunespin_bisec_v2.py index b018698..f8a7219 100644 --- a/basalt_buff_tunespin_bisec_v2.py +++ b/basalt_buff_tunespin_bisec_v2.py @@ -26,7 +26,7 @@ phnorm_pw = False use_local_storage = True -# use_local_storage = False +use_local_storage = False cec=float(sys.argv[3]) # targetpH = 6.0 @@ -38,12 +38,12 @@ dep_sample = 0.15 added_sp = 'gbas' -added_sp = 'cc' +# added_sp = 'cc' # mxing style imix = 2 # homogeneous -# imix = 3 # tilling +imix = 3 # tilling ztot_lab = 0.5 @@ -54,7 +54,7 @@ water_frac = 5. water_frac = 1. -water_frac = 2.5 +# water_frac = 2.5 catlist = ['ca','mg','k','na'] @@ -90,10 +90,9 @@ runname_lab = expid+'_'+added_sp+'_lab_tpH'+sys.argv[1].replace('.','p')+'_tau'+sys.argv[2].replace('.','p') # runname = 'chk_iter_incl2nd_basalt_tpH'+sys.argv[1].replace('.','p')+'_tau'+sys.argv[2].replace('.','p') -# outdir='/storage/scratch1/0/ykanzaki3/scepter_output/' -outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' +outdir = '../scepter_output/' if use_local_storage: - outdir_src = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir_src = '../scepter_output/' outdir = os.environ['TMPDIR'] + '/scepter_output/' for runname in [spinup_field,spinup_lab]: @@ -131,7 +130,8 @@ # fdust = 10000 exename = 'scepter' -exename_src = 'scepter_test' +# exename_src = 'scepter_test' +exename_src = 'scepter' to = ' ' where = '/' @@ -148,8 +148,8 @@ with open(dst, 'w') as file: file.writelines(data) -if added_sp == 'gbas': dustsrc = 'dust_gbasalt.in' -if added_sp == 'cc': dustsrc = 'dust_lime.in' +if added_sp == 'gbas': dustsrc = './data/dust_gbasalt.in' +if added_sp == 'cc': dustsrc = './data/dust_lime.in' dustdst = 'dust.in' os.system('cp ' + dustsrc + to + outdir + runname_field + where + dustdst) @@ -229,7 +229,7 @@ file.writelines(data) filename = '2ndslds.in' -srcfile = '/storage/coda1/p-creinhard3/0/ykanzaki3/PyWeath/data/2ndslds_def.in' +srcfile = './data/2ndslds_def.in' make_inputs.get_input_sld_properties( outdir=outdir ,runname=runname_field @@ -245,11 +245,11 @@ # compile # os.system('make') -os.system('make --file=makefile_test') -for runname in [runname_field,runname_lab]: - if not os.path.exists( outdir + runname) : os.system('mkdir -p ' + outdir + runname) - # os.system('cp ' + exename + to + outdir + runname) - os.system('cp ' + exename_src + to + outdir + runname + where + exename) +# os.system('make --file=makefile_test') +# for runname in [runname_field,runname_lab]: + # if not os.path.exists( outdir + runname) : + # os.system('mkdir -p ' + outdir + runname) + # os.system('cp ' + exename_src + to + outdir + runname + where + exename) maxiter = 50 diff --git a/get_inputs.py b/get_inputs.py new file mode 100644 index 0000000..8bbbc52 --- /dev/null +++ b/get_inputs.py @@ -0,0 +1,216 @@ +import os +import numpy as np +from shutil import copyfile +import time + +def get_input_frame(outdir,runname): + infile = outdir + runname + '/frame.in' + + with open(infile) as f: + lines = [line for line in f] + + ztot = float(lines[1].split()[0]) + nz = int(lines[2].split()[0]) + ttot = float(lines[3].split()[0]) + temp = float(lines[4].split()[0]) + fdust = float(lines[5].split()[0]) + fdust2 = float(lines[6].split()[0]) + taudust = float(lines[7].split()[0]) + omrain = float(lines[8].split()[0]) + zom = float(lines[9].split()[0]) + poro = float(lines[10].split()[0]) + moistsrf = float(lines[11].split()[0]) + zwater = float(lines[12].split()[0]) + zdust = float(lines[13].split()[0]) + w = float(lines[14].split()[0]) + q = float(lines[15].split()[0]) + p = float(lines[16].split()[0]) + nstep = int(lines[17].split()[0]) + rstrt = (lines[18].split()[0]) + runid = (lines[20].split()[0]) + + return ztot,nz,ttot,temp,fdust,fdust2,taudust,omrain,zom,poro,moistsrf,zwater,zdust,w,q,p,nstep,rstrt,runid + +def get_input_switches(outdir,runname): + infile = outdir + runname + '/switches.in' + + with open(infile) as f: + lines = [line for line in f] + + w_scheme = int(lines[1].split()[0]) + mix_scheme = int(lines[2].split()[0]) + poro_iter = (lines[3].split()[0]) + sldmin_lim = (lines[4].split()[0]) + display = (lines[5].split()[0]) + disp_lim = (lines[6].split()[0]) + restart = (lines[7].split()[0]) + rough = (lines[8].split()[0]) + act_ON = (lines[9].split()[0]) + dt_fix = (lines[10].split()[0]) + cec_on = (lines[11].split()[0]) + dz_fix = (lines[12].split()[0]) + close_aq = (lines[13].split()[0]) + poro_evol = (lines[14].split()[0]) + sa_evol_1 = (lines[15].split()[0]) + sa_evol_2 = (lines[16].split()[0]) + psd_bulk = (lines[17].split()[0]) + psd_full = (lines[18].split()[0]) + season = (lines[19].split()[0]) + + return w_scheme,mix_scheme,poro_iter,sldmin_lim,display,disp_lim,restart,rough,act_ON,dt_fix \ + ,cec_on,dz_fix,close_aq,poro_evol,sa_evol_1,sa_evol_2,psd_bulk,psd_full,season + +def get_input_tracers(outdir,runname): + + sld_list = np.genfromtxt(outdir + runname + '/slds.in',dtype='str',skip_header=1) + aq_list = np.genfromtxt(outdir + runname + '/solutes.in',dtype='str',skip_header=1) + gas_list = np.genfromtxt(outdir + runname + '/gases.in',dtype='str',skip_header=1) + exrxn_list = np.genfromtxt(outdir + runname + '/extrxns.in',dtype='str',skip_header=1) + + try: + sld_list = list(sld_list) + except: + sld_list = [str(sld_list)] + try: + aq_list = list(aq_list) + except: + aq_list = [aq_list[0]] + try: + gas_list = list(gas_list) + except: + print(gas_list.shape) + gas_list = [str(gas_list)] + try: + exrxn_list = list(exrxn_list) + except: + exrxn_list = [exrxn_list[0]] + + return sld_list,aq_list,gas_list,exrxn_list + +def get_input_tracer_bounds(**kwargs): + outdir = kwargs.get('outdir', '../scepter_output/') + runname = kwargs.get('runname', 'test_input') + pr_list = kwargs.get('pr_list', []) + rain_list = kwargs.get('rain_list', []) + atm_list = kwargs.get('atm_list', [('pco2',3.16e-4),('po2',0.21),('pnh3',1e-50),('pn2o',1e-50)]) + + notes = [ + '** parent rock wt fraction (e.g., "ab 0.2" in one line and "ka 0.001" in the next) (if not specified assumed 1e-20)' + ,'** solute concs. of rain in mol/L (if not specified assumed 1e-20)' + ,'** atmospheric composition in atm (if not specified assumed 1 PAL)' + ] + + pr_list.insert(0,'') + rain_list.insert(0,'') + atm_list.insert(0,'') + + values_lists = [ + pr_list + ,rain_list + ,atm_list + ] + + filenames = [ + 'parentrock.in' + ,'rain.in' + ,'atm.in' + ] + + nn = 3 + + if not os.path.exists(outdir + runname): os.makedirs(outdir + runname) + + for k in range(nn): + values = values_lists[k] + filename = filenames[k] + n = len(values) + input_text = '' + for i in range(n): + if i==0: + input_text += notes[k] + '\n' + else: + try: + input_text += values[i][0] + '\t' + str(values[i][1]) + '\n' + except: + print(n,i,values) + time.sleep(100) + + input_file = outdir + runname + '/' + filename + + with open(input_file, 'w') as file: + file.write(input_text) + + print(input_text) + + + del pr_list[0] + del rain_list[0] + del atm_list[0] + +def get_input_sld_properties(outdir,runname,filename): + + + infile = outdir + runname + '/' + filename + + sld_data_list = [] + + + with open(infile) as f: + cnt2 = 0 + for line in f: + if cnt2 ==0: + cnt2 += 1 + continue + else: + cnt2 += 1 + a = line.split() + tmplist = [] + cnt = 0 + for item in a: + if cnt ==0: tmplist.append(item) + else: tmplist.append(float(item)) + cnt += 1 + + sld_data_list.append(tmplist) + + return sld_data_list + +def main(): + + outdir = '../scepter_output/' + # runname = 'US_cropland_311_sph_N_spintuneup_field' + runname = 'test_Pot7' + + ztot,nz,ttot,temp,fdust,fdust2,taudust,omrain,zom,poro,moistsrf,zwater,zdust,w,q,p,nstep,rstrt,runid = get_input_frame(outdir,runname) + print( + ztot,nz,ttot,temp,fdust,fdust2,taudust,omrain,zom,poro,moistsrf,zwater,zdust,w,q,p,nstep,rstrt,runid + ) + + w_scheme,mix_scheme,poro_iter,sldmin_lim,display,disp_lim,restart,rough,act_ON,dt_fix \ + ,cec_on,dz_fix,close_aq,poro_evol,sa_evol_1,sa_evol_2,psd_bulk,psd_full,season = get_input_switches(outdir,runname) + print( + w_scheme,mix_scheme,poro_iter,sldmin_lim,display,disp_lim,restart,rough,act_ON,dt_fix + ,cec_on,dz_fix,close_aq,poro_evol,sa_evol_1,sa_evol_2,psd_bulk,psd_full,season + ) + + sld_list,aq_list,gas_list,exrxn_list = get_input_tracers(outdir,runname) + print( + sld_list,aq_list,gas_list,exrxn_list + ) + + filename = 'cec.in' + sld_data_list = get_input_sld_properties(outdir,runname,filename) + print( + sld_data_list, + sld_data_list[0][1] + ) + + sld_data_list = get_input_sld_properties('./','data','dust_gbasalt.in') + print( + sld_data_list + ) + +if __name__ == '__main__': + main() + + \ No newline at end of file diff --git a/get_int_prof.py b/get_int_prof.py index 42f7911..00bdb26 100644 --- a/get_int_prof.py +++ b/get_int_prof.py @@ -1042,7 +1042,7 @@ def get_dis_frac(outdir,runname,sp): def main(): - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' runname = 'US_cropland_251_pph_N_cacl2_2p5_homo2_mgo_field_tpH7p0_tau10' dep_sample = 0.15 get_intph_int_site(outdir,runname,dep_sample) diff --git a/get_int_prof_time.py b/get_int_prof_time.py index 5ae961e..2124edd 100644 --- a/get_int_prof_time.py +++ b/get_int_prof_time.py @@ -26,11 +26,6 @@ def linave(dep,phdep,ztot): return a def get_ph_int_site(outdir,runname,dep_sample,i): - # outdir = '../pyweath_output/' - # runname = '' - # runname = 'Sheldon_A_fick_noiter_test_psdfullpbe_w1_fit_r2000_sig0p2_loop_q0p2_chkall' - - # dep_sample = float(sys.argv[1]) infile = outdir+runname+'/prof/prof_aq-{:03d}.txt'.format(i) data = np.loadtxt(infile,skiprows=1) @@ -49,11 +44,6 @@ def get_ph_int_site(outdir,runname,dep_sample,i): return phintval def get_ac_int_site(outdir,runname,dep_sample): - # outdir = '../pyweath_output/' - # runname = '' - # runname = 'Sheldon_A_fick_noiter_test_psdfullpbe_w1_fit_r2000_sig0p2_loop_q0p2_chkall' - - # dep_sample = float(sys.argv[1]) infile = outdir+runname+'/prof/prof_aq(ads%cec)-020.txt' data = np.loadtxt(infile,skiprows=1) diff --git a/get_int_prof_time_dep.py b/get_int_prof_time_dep.py index 8d0208d..ef5cbc5 100644 --- a/get_int_prof_time_dep.py +++ b/get_int_prof_time_dep.py @@ -26,11 +26,6 @@ def linave(dep,phdep,ztot): return a def get_ph_int_site(outdir,runname,idep,i): - # outdir = '../pyweath_output/' - # runname = '' - # runname = 'Sheldon_A_fick_noiter_test_psdfullpbe_w1_fit_r2000_sig0p2_loop_q0p2_chkall' - - # dep_sample = float(sys.argv[1]) infile = outdir+runname+'/prof/prof_aq-{:03d}.txt'.format(i) data = np.loadtxt(infile,skiprows=1) @@ -42,11 +37,6 @@ def get_ph_int_site(outdir,runname,idep,i): return pH_dep def get_ac_int_site(outdir,runname,dep_sample): - # outdir = '../pyweath_output/' - # runname = '' - # runname = 'Sheldon_A_fick_noiter_test_psdfullpbe_w1_fit_r2000_sig0p2_loop_q0p2_chkall' - - # dep_sample = float(sys.argv[1]) infile = outdir+runname+'/prof/prof_aq(ads%cec)-020.txt' data = np.loadtxt(infile,skiprows=1) diff --git a/get_soilpH_time.py b/get_soilpH_time.py index de45c09..9ddce9c 100644 --- a/get_soilpH_time.py +++ b/get_soilpH_time.py @@ -62,8 +62,7 @@ def calc_soilpH(outdir, runname_field, dep_sample, itime, **kwargs): temp_lab=25 - # outdir = '/storage/scratch1/0/ykanzaki3/scepter_output/' - # outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + # outdir = '../scepter_output/' exename = 'scepter' to = ' ' @@ -485,7 +484,7 @@ def calc_soilpH(outdir, runname_field, dep_sample, itime, **kwargs): def timeseries_soilpH(): - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' runname_field = 'US_cropland_297_sph_N_DIC_v2_spintuneup_field' runname_field = 'US_cropland_311_sph_N_DIC_v2_spintuneup_field' runname_field = 'US_cropland_311_sph_N_DIC_cacl2_2p5_all_spintuneup_field' @@ -526,7 +525,7 @@ def CECseries_bufferpH(act_ON,database_name): cec_list = list(range(21)) cec_list = list(range(13)) - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' res_list = [] @@ -575,7 +574,7 @@ def S2006_bufferpH( use_local_storage, ): - # outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + # outdir = '../scepter_output/' # runname_field = 'US_cropland_297_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_311_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_297_sph_N_DIC_act_all_spintuneup_field' @@ -643,7 +642,7 @@ def MK2010_series_soilpH( use_local_storage, ): - # outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + # outdir = '../scepter_output/' # runname_field = 'US_cropland_297_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_311_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_297_sph_N_DIC_act_all_spintuneup_field' @@ -752,7 +751,7 @@ def MK2010_series_soilpH_2( use_local_storage, ): - # outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + # outdir = '../scepter_output/' # runname_field = 'US_cropland_297_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_311_sph_N_DIC_v2_spintuneup_field' # runname_field = 'US_cropland_297_sph_N_DIC_act_all_spintuneup_field' @@ -844,12 +843,11 @@ def main(): # timeseries_soilpH() + # exit("end run") - # act_ON = False # act_ON = True # database_name = 'Sikora' - # database_name = 'Goldberg' # CECseries_bufferpH(act_ON,database_name) # exit("end run") @@ -857,10 +855,10 @@ def main(): - i_parallel = int(sys.argv[1]) - n_parallel = int(sys.argv[2]) + # i_parallel = int(sys.argv[1]) + # n_parallel = int(sys.argv[2]) - name_base = sys.argv[3] + # name_base = sys.argv[3] # name_base = 'test_Pot7_25C_v2_alpha' name_base = 'test_Pot7_25C_v2_poro0p4_alpha' name_base = 'test_Pot7_25C_v2_densqrtz_alpha' @@ -885,8 +883,9 @@ def main(): name_base = 'test_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') name_base = 'chk_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') + name_base = 'chk2_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' alpha_list = list(np.linspace(0.1,5.5,19)) alpha_list = [3.4] @@ -897,7 +896,7 @@ def main(): alpha_list = [1.1] alpha_list = [1.2] alpha_list = [1.3] - alpha_list = [1.4] + # alpha_list = [1.4] # alpha_list = [1.7] # alpha_list = [5.5] # alpha_list = [(0, 6)] @@ -990,122 +989,6 @@ def main(): - - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' - runname_field_list = ['US_cropland_{:d}_sph_N_DIC_act_DA_1p1_spintuneup_field'.format(i) for i in [296,297,311,371]] - dep_sample = 0.15 - include_N = True - include_Cl = False - include_Al = False - use_local_storage = True - - cnt = 0 - for runname_field in runname_field_list: - for include_DIC in [False,True]: - if cnt%n_parallel!=i_parallel: - cnt += 1 - continue - else: - cnt += 1 - - MK2010_series_soilpH_2( - outdir, - runname_field, - dep_sample, - include_N, - include_Cl, - include_Al, - include_DIC, - use_local_storage, - ) - - - exit("end run") - - i_parallel = int(sys.argv[1]) - n_parallel = int(sys.argv[2]) - - name_base = sys.argv[3] - - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' - omrain_list = [100,300,600,900] - # ca_list = [1e-5,1e-4,1e-3,1e-2] - ca_list = [1e-6,1e-5,1e-4,1e-3] - cec_list = [5,10,20,30] - # cec_list = [0.01,0.1,1] - logkh_list = [5,6,7] - nacl_list = [1e-4,1e-3,1e-2,1e-1] - nacl_list = [1e-3,1e-2,1e-1,1e0] # 0.01 0.1 1 10 kg/ha --> 1e-3, 1e-2, 1e-1, 1e0 g/m2 - q_list = [1e-6,1e-4,1e-2,1e-1,1e0] - cnt = 0 - for omrain in omrain_list: - for ca in ca_list: - for cec in cec_list: - # for logkh in logkh_list: - # for nacl in nacl_list: - for q in q_list: - - if cnt%n_parallel!=i_parallel: - cnt += 1 - continue - else: - cnt += 1 - - # runname_field = 'test_inrt-g2-nacl-pco2_{:d}resp_{:d}lognacl_{:d}cec_{:d}logkh'.format( - # runname_field = 'test_inrt-g2-ca-pco2_{:d}resp_{:d}logca_{:d}cec_{:d}logkh'.format( - # int(omrain), - # int(-np.log10(ca)), - # int(cec), - # int(logkh), - # ) - - # if cec == int(cec): - # runname_field = 'test_inrt-g2-ca-pco2_{:d}resp_{:d}logca_{:d}cec_{:d}logkh'.format( - # int(omrain), - # int(-np.log10(ca)), - # int(cec), - # int(logkh), - # ) - # else : - # runname_field = 'test_inrt-g2-ca-pco2_{:d}resp_{:d}logca_{}cec_{:d}logkh'.format( - # int(omrain), - # int(-np.log10(ca)), - # str(cec).replace('.','p'), - # int(logkh), - # ) - - # runname_field = name_base + '_{:d}resp_{:d}logca_{:d}lognacl_{:d}cec'.format( - # int(omrain), - # int(-np.log10(ca)), - # int(-np.log10(nacl)), - # int(cec), - # ) - - runname_field = name_base + '_{:d}resp_{:d}logca_{:d}logq_{:d}cec'.format( - int(omrain), - int(-np.log10(ca)), - int(-np.log10(q)), - int(cec), - ) - - dep_sample = 0.15 - include_N = False - include_Cl = True - include_Al = False - use_local_storage = True - for include_DIC in [False,True]: - MK2010_series_soilpH( - outdir, - runname_field, - dep_sample, - include_N, - include_Cl, - include_Al, - include_DIC, - use_local_storage, - ) - - if __name__ == '__main__': diff --git a/get_soilpH_time_dep.py b/get_soilpH_time_dep.py index 474b933..6f1fd0e 100644 --- a/get_soilpH_time_dep.py +++ b/get_soilpH_time_dep.py @@ -38,8 +38,7 @@ def calc_soilpH(inputlist): temp_lab=25 - # outdir = '/storage/scratch1/0/ykanzaki3/scepter_output/' - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' exename = 'scepter' to = ' ' @@ -297,7 +296,7 @@ def main(): field_dir = sys.argv[1] - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' for itime in range(n_time): file = outdir + field_dir + '/prof/soil_ph_{:03d}.txt'.format(itime+1) diff --git a/make_inputs.py b/make_inputs.py index 08190fb..f96a886 100644 --- a/make_inputs.py +++ b/make_inputs.py @@ -4,7 +4,7 @@ import time def get_input_frame(**kwargs): - outdir = kwargs.get('outdir', '/storage/scratch1/0/ykanzaki3/scepter_output/') + outdir = kwargs.get('outdir', '../scepter_output/') runname = kwargs.get('runname', 'test_input') ztot = kwargs.get('ztot', 0.5) nz = kwargs.get('nz', 30) @@ -99,7 +99,7 @@ def get_input_frame(**kwargs): print(input_text) def get_input_switches(**kwargs): - outdir = kwargs.get('outdir', '/storage/scratch1/0/ykanzaki3/scepter_output/') + outdir = kwargs.get('outdir', '../scepter_output/') runname = kwargs.get('runname', 'test_input') w_scheme = kwargs.get('w_scheme', 1) mix_scheme = kwargs.get('mix_scheme', 1) @@ -189,7 +189,7 @@ def get_input_switches(**kwargs): print(input_text) def get_input_tracers(**kwargs): - outdir = kwargs.get('outdir', '/storage/scratch1/0/ykanzaki3/scepter_output/') + outdir = kwargs.get('outdir', '../scepter_output/') runname = kwargs.get('runname', 'test_input') sld_list = kwargs.get('sld_list', []) aq_list = kwargs.get('aq_list', []) @@ -250,7 +250,7 @@ def get_input_tracers(**kwargs): del exrxn_list[0] def get_input_tracer_bounds(**kwargs): - outdir = kwargs.get('outdir', '/storage/scratch1/0/ykanzaki3/scepter_output/') + outdir = kwargs.get('outdir', '../scepter_output/') runname = kwargs.get('runname', 'test_input') pr_list = kwargs.get('pr_list', []) rain_list = kwargs.get('rain_list', []) @@ -310,7 +310,7 @@ def get_input_tracer_bounds(**kwargs): del atm_list[0] def get_input_sld_properties(**kwargs): - outdir = kwargs.get('outdir', '/storage/scratch1/0/ykanzaki3/scepter_output/') + outdir = kwargs.get('outdir', '../scepter_output/') runname = kwargs.get('runname', 'test_input') sld_varlist = kwargs.get('sld_varlist', []) filename = kwargs.get('filename', 'kinspc.in') @@ -448,8 +448,7 @@ def get_input_climate_temp(**kwargs): def main(): - # outdir = '/storage/scratch1/0/ykanzaki3/scepter_output/' - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' runname = 'test_input' ztot=0.5 @@ -563,7 +562,7 @@ def main(): ,atm_list=atm_list ) filename = 'dust.in' - srcfile = '/storage/coda1/p-creinhard3/0/ykanzaki3/PyWeath/data/dust_gbasalt.in' + srcfile = './data/dust_gbasalt.in' sld_varlist =[] get_input_sld_properties( outdir=outdir @@ -589,8 +588,8 @@ def main(): ,sld_varlist=sld_varlist ) filename = '2ndslds.in' - srcfile = '/storage/coda1/p-creinhard3/0/ykanzaki3/PyWeath/data/2ndslds_def.in' - srcfile = '/storage/coda1/p-creinhard3/0/ykanzaki3/PyWeath/data/2ndslds_rm_al2o3.in' + srcfile = './data/2ndslds_def.in' + srcfile = './data/2ndslds_rm_al2o3.in' sld_varlist =[] get_input_sld_properties( outdir=outdir diff --git a/readme.txt b/readme.txt index e4b8350..503bed0 100644 --- a/readme.txt +++ b/readme.txt @@ -1,16 +1,68 @@ scepter -scepter.f90 - source code -makefile - compile (data dir needs to be specified) +scepter.f90 - source code +makefile - compile (data dir in L50 needs to be specified) scripts: -get_int_prof.py - contain functions to get output data -make_inputs.py - contain functions to make input data -tunespin_3_newton_inert_buff.py - run 3 variable iterations (output dir needs to be specified) -tunespin_3_newton_inert_buff_v2.py - same as above but with an option to use soil pH with CaCl2 -water_amb.py - sample porewater and leave it in the beaker in lab (output dir needs to be specified; no use for pipeline?) -basalt_buff_tunespin_bisec.py - run basalt application for a target pH etc. (output dir needs to be specified) -basalt_buff_tunespin_bisec_v2.py - same as above but with an option to use soil pH with CaCl2 -sub_jobs.py - submit multiple jobs (for spinup and basalt exp) -run_a_shell.sbatch - submit a single job (for GT cluster) \ No newline at end of file +get_int_prof.py - contain functions to get output data +get_int_prof_time.py - same above but with specifing the model time +get_int_prof_time_dep.py + - same above but with specifing the model depth +get_soilpH_time.py - contain functions to calculate soil pH +get_soilpH_time_dep.py - same above but with specifing the model depth +get_inputs.py - contain functions to retrieve input data +make_inputs.py - contain functions to make input data +tunespin_3_newton_inert_buff_v2_clean.py + - conduct 3 variable field-run iterations + (output dir needs to be specified) +basalt_buff_tunespin_bisec_v2.py + - conduct field-run iteration to get to a target soil/pw pH + (output dir needs to be specified) +spinup.py - run a spin-up run +spinup_inert.py - run series of run with bulk speces varying CECs +spinup_inrt2.py - run series of run with bulk and OM speces varying CECs etc. + + + +$$ running example simulations in GMD paper $$ + +1. test for Sikora buffer used in in Section 3 +1.1. in-silico field samples + a. modify outdir in L20 of spinup_inert.py + b. type: python3 spinup_inert.py +1.1. buffer pH for field samples in silico + a. modify outdir in L528 and + undo comments-out in L849-852 of get_soilpH_time.py + b. type: python3 get_soilpH_time.py + +2. mesocosm experiment in Section 3 +2.1. field simulation + a. modify outdir in L159 of spinup_inrt2.py + b. type: python3 spinup_inrt2.py + +2.2. laboratory simulation + a. modify outdir in L528, L888 of get_soilpH_time.py + b. type: python3 get_soilpH_time.py + +3. alkalinity requiremnt for ERW in Section 4 +3.1. field + laboratory simulation +3.1.1. spin/tune-up + a. modify outdir in L148 of tunespin_3_newton_inert_buff_v2_clean.py + b. type: python3 tunespin_3_newton_inert_buff_v2_clean.py spinup_run_name + 21.103289732688683 6.058006742238197 20.980309042371502 2.0516666666666667 + 8.222189843654622 0.282726679550165 0.35136107875550837 0.0010131311683626316 + 1.005952418781816 + [ note that the runtime inputs are to specify: run ID, CEC (cmol/kg), target soil pH, + target exchange acidity (%CEC), target soil OM (wt%), temperature (oC), + moisture, runoff (m/yr), erosion rate (m/yr), nitrification rate (gN/m2/yr) ] +3.1.2. basalt application + a. modify outdir in L93 of basalt_buff_tunespin_bisec_v2.py and + option of using soil/porewater pH (phnorm_pw=False/True, L25-26) + b. type: python3 basalt_buff_tunespin_bisec_v2.py 6.2 1 21.103289732688683 + basalt_run_name spinup_run_name + [ note that the runtime inputs are to specify: + target pH, duration of run, CEC (cmol/kg), basalt run name, spin/tune-up run name ] +3.2. calculating soil pH prodiles + a. modify outdir in L49, L299 of get_soilpH_time_dep.py + b. type: get_soilpH_time_dep.py basalt_run_name \ No newline at end of file diff --git a/spinup.py b/spinup.py index 3a026a7..54d7cf6 100644 --- a/spinup.py +++ b/spinup.py @@ -78,11 +78,12 @@ def run_a_scepter_run( # compile exename = 'scepter' - exename_src = 'scepter_test' + exename_src = 'scepter' + # exename_src = 'scepter_test' to = ' ' where = '/' - # os.system('make') - os.system('make --file=makefile_test') + os.system('make') + # os.system('make --file=makefile_test') if not os.path.exists( outdir + runname) : os.system('mkdir -p ' + outdir + runname) os.system('cp ' + exename_src + to + outdir + runname + where + exename) @@ -208,7 +209,7 @@ def run_a_scepter_run( def main(): - outdir_src = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir_src = '../scepter_output/' runname = 'test' # >>>> input variables of interests diff --git a/spinup_inert.py b/spinup_inert.py index b0032d6..8f29b26 100644 --- a/spinup_inert.py +++ b/spinup_inert.py @@ -5,7 +5,7 @@ def run_closed_inert_cec(cec,act_ON_IN,data): use_local_storage = True - # use_local_storage = False + use_local_storage = False # cec = 0 @@ -17,7 +17,7 @@ def run_closed_inert_cec(cec,act_ON_IN,data): alpha = 3.4 - outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir = '../scepter_output/' if use_local_storage: outdir = os.environ['TMPDIR'] + '/scepter_output/' # runname = 'test_inert_spinup_cec_{:d}'.format(cec) if act_ON_IN: runname = 'test_inert_spinup_act_cec_{:d}'.format(cec) + '_'+ data @@ -25,11 +25,12 @@ def run_closed_inert_cec(cec,act_ON_IN,data): # compile exename = 'scepter' - exename_src = 'scepter_test' + exename_src = 'scepter' + # exename_src = 'scepter_test' to = ' ' where = '/' - # os.system('make') - os.system('make --file=makefile_test') + os.system('make') + # os.system('make --file=makefile_test') if not os.path.exists( outdir + runname) : os.system('mkdir -p ' + outdir + runname) os.system('cp ' + exename_src + to + outdir + runname + where + exename) @@ -154,7 +155,7 @@ def run_closed_inert_cec(cec,act_ON_IN,data): ) filename = 'dust.in' - srcfile = '/storage/coda1/p-creinhard3/0/ykanzaki3/PyWeath/data/dust_gbasalt.in' + srcfile = './data/dust_gbasalt.in' sld_varlist =[] make_inputs.get_input_sld_properties( outdir=outdir @@ -188,7 +189,7 @@ def run_closed_inert_cec(cec,act_ON_IN,data): if use_local_storage: src = outdir + runname - dst = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + runname + dst = '../scepter_output/' + runname if not os.path.exists(dst): shutil.copytree(src, dst) @@ -198,13 +199,11 @@ def run_closed_inert_cec(cec,act_ON_IN,data): def main(): - # cec_list = [0,1,2,4,8,16,32,64] cec_list = list(range(21)) act_ON_IN = False act_ON_IN = True data = 'Sikora' - # data = 'Goldberg' for cec in cec_list: run_closed_inert_cec(cec,act_ON_IN,data) diff --git a/spinup_inrt2.py b/spinup_inrt2.py index 145dcad..a882143 100644 --- a/spinup_inrt2.py +++ b/spinup_inrt2.py @@ -15,7 +15,6 @@ def spin_inert(outdir_src,runname, act_ON, ): - # outdir_src = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' # runname = 'test' # cec = 10.0 @@ -157,7 +156,7 @@ def spin_inert(outdir_src,runname, def main(): - outdir_src = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir_src = '../scepter_output/' runname = 'test_Pot7_25C_alpha0p0' runname = 'test_Pot7_25C_alpha3p1' @@ -188,7 +187,7 @@ def main(): cl_tmp = 2.9e-4 # double CEC OM - 120, pwpH=6.1, alpha = 1.1 cl_tmp = 2.8e-4 # double CEC OM - 120, pwpH=6.1, alpha = 1.2 cl_tmp = 2.68e-4 # double CEC OM - 120, pwpH=6.1, alpha = 1.3 - cl_tmp = 2.55e-4 # double CEC OM - 120, pwpH=6.1, alpha = 1.4 + # cl_tmp = 2.55e-4 # double CEC OM - 120, pwpH=6.1, alpha = 1.4 ph_pw = 6.68 ph_pw = 6.3 @@ -269,16 +268,17 @@ def main(): if cec_1<0: exit("cec_1<0") - i_parallel = int(sys.argv[1]) - n_parallel = int(sys.argv[2]) + # i_parallel = int(sys.argv[1]) + # n_parallel = int(sys.argv[2]) - name_base = sys.argv[3] + # name_base = sys.argv[3] # name_base = 'test_Pot7_25C_v2_CEC10p0' name_base = 'test_Pot7_25C_v2_noact_alpha' name_base = 'test_Pot7_25C_v2_act_alpha' name_base = 'test_Pot7_25C_v2_act_2alpha' name_base = 'test_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') name_base = 'chk_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') + name_base = 'chk2_Pot7_25C_v2_act_2CEC{:.1f}-{:.1f}_pH{:.1f}_alpha'.format(cec_1,cec_2,ph_pw).replace('.','p') # name_base = 'test_Pot7_25C_v3_noact_alpha' act_ON = 'true' @@ -291,7 +291,7 @@ def main(): alpha_list = [(1.1, 1.1)] alpha_list = [(1.2, 1.2)] alpha_list = [(1.3, 1.3)] - alpha_list = [(1.4, 1.4)] + # alpha_list = [(1.4, 1.4)] for i in range(len(alpha_list)): diff --git a/sub_jobs.py b/sub_jobs.py deleted file mode 100644 index 999b008..0000000 --- a/sub_jobs.py +++ /dev/null @@ -1,177 +0,0 @@ -import os -import numpy as np -import time - -input_dir = './data/' -# input_dir = './' -input_file = 'friendship.txt' -input_file = 'GAdata.txt' -input_file = 'SL004810.txt' -skip_file = 'GA_fail.txt' -skip_file = 'SL004810_sph_fail.txt' -skip_file = 'SL004810_FAIL_sph_N_spintuneup_field.txt' - -data = np.loadtxt(input_dir+input_file,skiprows=1) -skipdata = np.loadtxt(input_dir+skip_file) -skipdata = [int(i) for i in skipdata] - -n_sample = data.shape[0] - -runs_chosen = [ - # 2,6,7 - # 7 - ] - -runtype = 'spinup' -runtype = 'basalt' - -for i in range(n_sample): - - - runtime = '12:20:00' - soft = 'python3' - - if 'friendship' in input_file: - - runid = data[i,0] - soilpH = data[i,1] - om = data[i,2] - cec = data[i,3] - acid = data[i,4] - buffpH = data[i,5] - - if acid ==0:acid = 0.1 - if soilpH ==5.6: soilpH = 5.7 - - elif 'GAdata' in input_file: - runid = data[i,0] - soilpH = data[i,8] - om = 5 # no data given - cec = data[i,1] - acid = data[i,6] - buffpH = 7 # no data given - - elif 'SL004810' in input_file: - - runid = data[i,0] - soilpH = data[i,1] - om = data[i,2] - cec = data[i,3] - acid = data[i,4] - buffpH = 7 # no data given - - - if len(runs_chosen)!=0 and runid not in runs_chosen: continue - - if runtype == 'spinup': - - code = 'tunespin_3_newton_inert_buff.py' - code = 'tunespin_3_newton_inert_buff_v2.py' - # code = 'tunespin_3b_newton_inert_buff.py' - - if code == 'tunespin_3b_newton_inert_buff.py' and np.isnan(buffpH): continue - - # runname = 'potato_'+str(int(runid)) - # runname = 'potato_pw_'+str(int(runid)) - # runname = 'potato_fert_'+str(int(runid)) - # runname = 'potato_fert_pw_'+str(int(runid)) - runname = 'potato_buff_'+str(int(runid)) - runname = 'GA_OM5_sph_'+str(int(runid)) - runname = input_file.replace('.txt','')+'_'+str(int(runid)) - - # runname += '_sph' - runname += '_sph_N' - - jobname = runname - - inputs = [ - runname - ,str(cec) - ,str(soilpH) - ,str(acid) - ,str(om) - ] - - if code == 'tunespin_3b_newton_inert_buff.py': - inputs = [ - runname - ,str(cec) - ,str(soilpH) - ,str(buffpH) - ,str(om) - ] - - input_runtime = ' '.join(inputs) - - command = 'sbatch --output='+jobname+'.out --job-name='+jobname+' --time='+runtime + ' run_a_shell.sbatch ' + soft + ' ' + code \ - + ' ' + input_runtime - - print(command) - os.system(command) - - time.sleep(5) - - elif runtype == 'basalt': - - if acid <= 0.1 : continue - - # if 'GAdata' in input_file and int(runid) in skipdata: continue - if 'SL004810' in input_file and int(runid) in skipdata: continue - - code = 'basalt_buff_tunespin_bisec.py' - code = 'basalt_buff_tunespin_bisec_v2.py' - - # runname = 'potato_bas_'+str(int(runid)) - # runname = 'potato_bas_pw_'+str(int(runid)) - # runname = 'potato_bas_fert_50u_'+str(int(runid)) - # runname = 'potato_bas_50u_'+str(int(runid)) - # runname = 'potato_bas_50u_'+str(int(runid)) - # runname = 'GA_OM5_sph_bas_10u_'+str(int(runid)) - runname = input_file.replace('.txt','')+'_'+str(int(runid)) - - # spinname = 'potato_'+str(int(runid)) - # spinname = 'potato_pw_'+str(int(runid)) - # spinname = 'potato_fert_'+str(int(runid)) - spinname = 'GA_OM5_sph_'+str(int(runid)) - spinname = input_file.replace('.txt','')+'_'+str(int(runid)) - - # spinname += '_N' - # runname += '_N_10u_bas' - # spinname += '_h' - # runname += '_10u_bas' - # spinname += '_sph' - # runname += '_sph_10u_bas' - spinname += '_sph_N' - runname += '_sph_N_10u_bas' - - targetpHs = [6.5, 6.8] - targetpHs = [5.4, 5.6, 5.8] - targettaus = [1] - - for targetpH in targetpHs: - - if soilpH >=targetpH: continue - - for targettau in targettaus: - - jobname = runname+'_'+'_'.join([str(targetpH),str(targettau)]).replace('.','p') - - - inputs = [ - str(targetpH) - ,str(targettau) - ,str(cec) - ,runname - ,spinname - ] - - input_runtime = ' '.join(inputs) - - command = 'sbatch --output='+jobname+'.out --job-name='+jobname+' --time='+runtime + ' run_a_shell.sbatch ' + soft + ' ' + code \ - + ' ' + input_runtime - - print(command) - os.system(command) - - time.sleep(5) - # break \ No newline at end of file diff --git a/tunespin_3_newton_inert_buff_v2_clean.py b/tunespin_3_newton_inert_buff_v2_clean.py index 6415523..01478a6 100644 --- a/tunespin_3_newton_inert_buff_v2_clean.py +++ b/tunespin_3_newton_inert_buff_v2_clean.py @@ -25,7 +25,7 @@ alphase = 'amal' use_CaCl2 = False -use_CaCl2 = True +# use_CaCl2 = True phnorm_pw = True phnorm_pw = False @@ -37,10 +37,10 @@ use_local_storage = False make_initial_geuss = True -# make_initial_geuss = False +make_initial_geuss = False liming = True -# liming = False +liming = False limesp = 'cao' limemwt = 56.079 @@ -50,7 +50,7 @@ water_frac = 5. water_frac = 1. -water_frac = 2.5 +# water_frac = 2.5 iter_max = 120 iter_max = 3000 @@ -83,7 +83,7 @@ make_initial_geuss = True if make_initial_geuss: - outdir_guess = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + outdir_guess = '../scepter_output/' try: # filename = '/iteration_tmp.res' @@ -145,8 +145,6 @@ # dep_sample = 0.25 -outdir = '/storage/scratch1/0/ykanzaki3/scepter_output/' -outdir = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' outdir = '../scepter_output/' if use_local_storage: outdir = os.environ['TMPDIR'] + '/scepter_output/' simid = sys.argv[1] @@ -155,11 +153,12 @@ # compile exename = 'scepter' -exename_src = 'scepter_test' +exename_src = 'scepter' +# exename_src = 'scepter_test' to = ' ' where = '/' -# os.system('make') -os.system('make --file=makefile_test') +os.system('make') +# os.system('make --file=makefile_test') prev_iter_exist = False if not os.path.exists( outdir + runname_field) : os.system('mkdir -p ' + outdir + runname_field) @@ -383,10 +382,13 @@ filename = 'psdrain.in' srcfile = './data/psdrain_320um.in' +sld_varlist = [ (5e-6,0.2,1), (5e-6,0.2,1), (5e-6,0.2,1), (5e-6,0.2,1), ] +# sld_varlist = [ (1e-6,0.2,1), (2e-6,0.2,1), (5e-6,0.2,1), (20e-6,0.2,1), ] make_inputs.get_input_sld_properties( outdir=outdir ,runname=runname_field ,filename = filename + # ,sld_varlist=sld_varlist ,srcfile = srcfile ) @@ -551,6 +553,7 @@ error = 1e4 tol = 1e-4 tol = 1e-3 +# tol = 1e-2 cnt = 1 res_list = [] @@ -1081,7 +1084,7 @@ def run_a_single_set(ca,logkh,omrain_field): if use_local_storage: for runname in [runname_field,runname_lab]: src = outdir + runname - dst = '/storage/coda1/p-creinhard3/0/ykanzaki3/scepter_output/' + runname + dst = '../scepter_output/' + runname if not os.path.exists(dst): shutil.copytree(src, dst)