#!/usr/bin/env python from __future__ import print_function, division, absolute_import from builtins import range import sys import os import string import numpy as np import h5py from ncutils import read_nc from scipy.interpolate import interp1d from itertools import product dir_pyartdeco = "/home/mathieu/work/RTM/artdeco/pyartdeco/" dir_artdeco = "/home/mathieu/work/RTM/artdeco/fortran/" os.environ["DIR_PYARTDECO"] = dir_pyartdeco os.environ["DIR_ARTDECO"] = dir_artdeco sys.path.append(dir_pyartdeco+'tools/') import pyartdeco_utils as ad dir_opt = '/rfs/proj/artdeco_lib/opt/' from scipy.integrate import simps dir_runlib_log = "/tmp/runlib_log/" from runlib.condor import CondorPool as Pool #runlib_method = "none" runlib_method = "condor" delta_sumf11 = 1e-3 trunc_model ='dm' baum_dir = '/rfs/data/baum_opt/' cdffile = 'GeneralHabitMixture_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' ang_s = read_nc(baum_dir+cdffile,'phase_angles') ang_s = ang_s.astype(np.float64) ang_store = np.append(ang_s[ang_s <= 80.0], np.linspace(80.5, 175.5, num=191, endpoint=True)) ang_store = np.append(ang_store, ang_s[ang_s >= 176.0]) ang_store = np.sort(ang_store)[::-1] mu_store = np.cos(ang_store*np.pi/180.0) nang_store = len(ang_store) wtol = 0.001/100. rho_liq_wat = 1.0e-12 # g micron-3 ######################################################################################################### def gamma_size_dist(r, reff, veff): # as defined in mie3 # reff is the effective radius, reff must be positive # veff is the effective variance, veff must be positive and less than 0.5 # if veff is larger than 1/3, n(r) is singular at r = 0 #from scipy.special import gammaln from scipy.special import loggamma alpha = 1.0/veff - 3.0 b = 1.0/(veff*reff) alpha1= alpha+1.0 #logC = alpha1 * np.log(b) - gammaln(alpha1) logC = alpha1 * np.log(b) - loggamma(alpha1) n = np.exp( logC + alpha * np.log(r) - b * r) return n def logn_size_dist(r, reff, veff): # as defined in mie3 root2p = np.sqrt(np.pi*2.0) rg = reff/(1.0+veff)**2.5 flogrg = np.log(rg) flogsi = np.sqrt(np.log(1.0+veff)) C = 1.0/(root2p*flogsi) fac = -0.50/(flogsi*flogsi) n = C * np.exp( fac*(np.log(r)-flogrg)**2.0 ) / r return n ######################################################################################### def compute_opt(wvl, reff, veff, gamma): if gamma: ptcle_name = 'liq_cl_gamma_veff_%.2f'%veff+'_reff_%.2f'%reff else: ptcle_name = 'liq_cl_veff_%.2f'%veff+'_reff_%.2f'%reff saveroot = ptcle_name keywords = ' verbose opt_only nowarn' artdeco_in=[] artdeco_in.append('# \n') artdeco_in.append('# MAIN INPUT FILE FOR ARTDECO PROGRAM \n') artdeco_in.append('# \n') artdeco_in.append('######################## \n') artdeco_in.append('# keywords (Ex: none, verbose,...) \n') artdeco_in.append(keywords+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# outfiles root name \n') artdeco_in.append(saveroot+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# mode \n') artdeco_in.append(' mono \n') artdeco_in.append('#######################\n') artdeco_in.append('# filters\n') artdeco_in.append(' none \n') artdeco_in.append('######################## \n') artdeco_in.append('# Wavelengths (microns) \n') s= '' for iwl in range(len(wvl)): s = s+' %.9f'%wvl[iwl] artdeco_in.append(s+' \n') artdeco_in.append('###################### \n') artdeco_in.append('# Particles \n') artdeco_in.append('# type interp. H-G Tau_550 vertical distribution type vdist parameters (km) \n') artdeco_in.append(ptcle_name+' no \n') artdeco_in.append('# \n') # write artdeco_in.dat f = open(dir_artdeco+'input/artdeco_in_'+saveroot, 'w') for j in range(len(artdeco_in)): f.write(artdeco_in[j]) f.close() prop_in=[] prop_in.append('# \n') prop_in.append('# ARTDECO : File containig properties for new particles (aerosols or hygrosols) \n') prop_in.append('# needed to obtain optical properties (file opt_*.dat) \n') prop_in.append('# \n') prop_in.append('# Material \n') prop_in.append('lwat \n') prop_in.append('# optical model \n') prop_in.append('mie3 \n') prop_in.append('#################################################### \n') prop_in.append('# INPUT FOR THE MEERHOFF MIE PROGRAM VERSION 3.0 size distribution \n') prop_in.append('# truncation of Mie sum (delta) \n') prop_in.append('1.00E-12 \n') prop_in.append('# cutoff value of the size distribution (cutoff) \n') prop_in.append('1.00E-06 \n') prop_in.append('# index size distribution (idis) \n') if gamma: prop_in.append('2 \n') else: prop_in.append('5 \n') prop_in.append('# number of subintervals for r (nsubr) \n') prop_in.append('20 \n') prop_in.append('# number of Gauss points in sub (ngaur) \n') prop_in.append('500 \n') prop_in.append('# PAR1 \n') prop_in.append('%.2f'%reff+' \n') prop_in.append('# PAR2 \n') prop_in.append('%.2f'%veff+' \n') prop_in.append('# PAR3 \n') prop_in.append('0.0 \n') prop_in.append('# \n') # write f = open(dir_artdeco+'lib/prop/prop_'+ptcle_name+'.dat', 'w') for j in range(len(prop_in)): f.write(prop_in[j]) f.close() # run ARTDECO one time os.system(dir_artdeco.strip()+"src/artdeco artdeco_in_"+saveroot) os.system("rm -rf "+dir_artdeco+'input/artdeco_in_'+saveroot) os.system("rm -rf "+dir_artdeco+'out/'+saveroot+'_unknown') return ################################################### def computeBetal(wvl, reff, veff, nbetal, gamma): keywords = ' betal_only print_betal print_recomp verbose' if gamma: ptcle = 'liq_cl_gamma_veff_%.2f'%veff+'_reff_%.2f'%reff else: ptcle = 'liq_cl_veff_%.2f'%veff+'_reff_%.2f'%reff saveroot = "get_betal_"+ptcle+"_%i"%nbetal if not os.path.isfile(dir_opt+'/opt_'+ptcle+'.dat'): print('No file'+dir_opt+'/opt_'+ptcle+'.dat !!!') return artdeco_in=[] artdeco_in.append('# \n') artdeco_in.append('# MAIN INPUT FILE FOR ARTDECO PROGRAM \n') artdeco_in.append('# \n') artdeco_in.append('######################## \n') artdeco_in.append('# keywords (Ex: none, verbose,...) \n') artdeco_in.append(keywords+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# outfiles root name \n') artdeco_in.append(saveroot+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# mode \n') artdeco_in.append(' mono \n') artdeco_in.append('#######################\n') artdeco_in.append('# filters\n') artdeco_in.append(' none \n') artdeco_in.append('######################## \n') artdeco_in.append('# Wavelengths (microns) \n') s= '' for iwl in range(len(wvl)): s = s+' %.9f'%wvl[iwl] artdeco_in.append(s+' \n') artdeco_in.append('###################### \n') artdeco_in.append('# Particles \n') artdeco_in.append('# type interp. H-G Tau_550 vertical distribution type vdist parameters (km) \n') artdeco_in.append(ptcle+' no \n') artdeco_in.append('# \n') artdeco_in.append('########################## \n') artdeco_in.append('# Truncation method (none, dfit, DM, Potter) \n') artdeco_in.append(trunc_model+' \n') artdeco_in.append('########################## \n') artdeco_in.append('# Number of Betal \n') artdeco_in.append('%i'%nbetal+' \n') artdeco_in.append('# \n') # write artdeco_in.dat f = open(dir_artdeco+'input/artdeco_in_'+saveroot+'.dat', 'w') for j in range(len(artdeco_in)): f.write(artdeco_in[j]) f.close() # run ARTDECO one time os.system(dir_artdeco.strip()+"src/artdeco artdeco_in_"+saveroot+".dat") os.system("rm -rf "+dir_artdeco+'input/artdeco_in_'+saveroot+'.dat') return ################################################### def doit(wvl, reff, veff, nbetal, res_file, gamma=False, write_library=True): nwvl = len(wvl) n_nbetal = len(nbetal) nbetalmax = np.max(nbetal) nreff = len(reff) nveff = len(veff) r4vol = np.linspace(1e-3, 2000., num=2000) # create empty hdf5 file to be filled # as result come if (not os.path.isfile(dir_opt+res_file)) and (write_library): f = h5py.File(dir_opt+res_file, 'a') grp = f.require_group("liquid_cloud") subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (nwvl,), dtype='float64') dset = subgrp.require_dataset("reff", (nreff,), dtype='float64') dset = subgrp.require_dataset("veff", (nveff,), dtype='float64') dset = subgrp.require_dataset("phase_angles", (nang_store,), dtype='float64') dset = subgrp.require_dataset("mu", (nang_store,), dtype='float64') dset = subgrp.require_dataset("nbetal", (n_nbetal,), dtype='i') dset = subgrp.require_dataset("ibetal", (nbetalmax+1 ,), dtype='i') dset[...] = np.arange(nbetalmax+1) subgrp = grp.require_group("data") dset = subgrp.require_dataset("single_scattering_albedo", (nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("reff,veff,wavelengths")) dset = subgrp.require_dataset("Cext", (nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("reff,veff,wavelengths")) dset.attrs.create("unit", np.string_("micron^2")) dset = subgrp.require_dataset("extinction_coefficient_over_wc", (nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("reff,veff,wavelengths")) dset.attrs.create("unit", np.string_("m^2g-1")) dset = subgrp.require_dataset("integrate_p11", (nreff,nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("reff,veff,wavelengths")) dset = subgrp.require_dataset("integrate_p11_recomp", (n_nbetal,nreff,nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("nbetal,reff,veff,wavelengths")) dset = subgrp.require_dataset("p11_phase_function", (nang_store,nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,reff,veff,wavelengths")) dset = subgrp.require_dataset("p21_phase_function", (nang_store,nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,reff,veff,wavelengths")) dset = subgrp.require_dataset("p34_phase_function", (nang_store,nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,reff,veff,wavelengths")) dset = subgrp.require_dataset("p44_phase_function", (nang_store,nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,reff,veff,wavelengths")) dset = subgrp.require_dataset("truncation_coefficient", (n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha1_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha2_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha3_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha4_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("beta1_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("beta2_betal", (nbetalmax+1, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p11_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p22_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p33_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p21_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p34_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p44_phase_function", (nang_store, n_nbetal, nreff, nveff,nwvl), dtype='float64') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,veff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) f.close() print('Compute optical properties ...') params = [] for (reff_tmp, veff_tmp) in product(reff, veff): params.append((wvl, reff_tmp, veff_tmp, gamma)) # we revert the param list to compute the bigger particles first params = params[::-1] print(" Number of compute_opt() call is", len(params)) if runlib_method == 'condor': pool = Pool(log=dir_runlib_log+'compute_opt_liqcl', progressbar=True, prio = -5) it = pool.map(compute_opt, *zip(*params)) else: it = map(compute_opt, *zip(*params)) for res in it: toto=0 if write_library: print('Compute Betal ...') params = [] for (reff_tmp, veff_tmp, nbetal_tmp) in product(reff, veff, nbetal): params.append((wvl, reff_tmp, veff_tmp, nbetal_tmp, gamma)) print(" Number of compute_betal() call is", len(params)) if runlib_method == 'condor': pool = Pool(log=dir_runlib_log+'compute_betal_liqcl', progressbar=True, prio = -5) it = pool.map(computeBetal, *zip(*params)) else: it = map(computeBetal, *zip(*params)) for res in it: toto=0 for iveff in range(nveff): for ireff in range(nreff): if gamma: ptcle_name = 'liq_cl_gamma_veff_%.2f'%veff[iveff]+'_reff_%.2f'%reff[ireff] else: ptcle_name = 'liq_cl_veff_%.2f'%veff[iveff]+'_reff_%.2f'%reff[ireff] print('Read opt, Betal and write HDF for ', ptcle_name) # open file during computation f = h5py.File(dir_opt+res_file, 'r') test = f['liquid_cloud']['data']['integrate_p11'][ireff,iveff,:] if not 0.0 in test[:]: print('particle already computed ', ptcle_name) # value already computed and stored continue f.close() trunc_coeff = None alpha1 = None alpha2 = None alpha3 = None alpha4 = None beta1 = None beta2 = None p11_recomp = None p22_recomp = None p33_recomp = None p44_recomp = None p21_recomp = None p34_recomp = None p11_recomp_integrate = None p11_integrate = None g = None ssa = None Cext = None p11 = None p44 = None p21 = None p34 = None Sext_o_WC = None wvl_tmp, ang_tmp, p11_tmp, p44_tmp, p21_tmp, p34_tmp, \ Cext_tmp, ssa_tmp, g_tmp = ad.get_opt_all(dir_opt, ptcle_name) # read optical properties for iwvl in range(nwvl): iwvl_tmp = np.argwhere( (wvl[iwvl]*(1+wtol) > wvl_tmp) & (wvl[iwvl]*(1-wtol) < wvl_tmp) )[0][0] if iwvl==0 : nang = len(ang_tmp) ang = np.zeros(nang) ang[:] = ang_tmp[:] trunc_coeff = np.zeros((n_nbetal, nwvl)) alpha1 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) alpha2 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) alpha3 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) alpha4 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) beta1 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) beta2 = np.zeros((nbetalmax+1, n_nbetal, nwvl)) p11_recomp = np.zeros((nang, n_nbetal, nwvl)) p22_recomp = np.zeros((nang, n_nbetal, nwvl)) p33_recomp = np.zeros((nang, n_nbetal, nwvl)) p44_recomp = np.zeros((nang, n_nbetal, nwvl)) p21_recomp = np.zeros((nang, n_nbetal, nwvl)) p34_recomp = np.zeros((nang, n_nbetal, nwvl)) p11_recomp_integrate = np.zeros((n_nbetal, nwvl)) p11_integrate = np.zeros(nwvl) g = np.zeros(nwvl) ssa = np.zeros(nwvl) Cext = np.zeros(nwvl) p11 = np.zeros((nang, nwvl)) p44 = np.zeros((nang, nwvl)) p21 = np.zeros((nang, nwvl)) p34 = np.zeros((nang, nwvl)) Sext_o_WC = np.zeros(nwvl) Cext[iwvl] = Cext_tmp[iwvl_tmp] ssa[iwvl] = ssa_tmp[iwvl_tmp] g[iwvl] = g_tmp[iwvl_tmp] if gamma: n_r = gamma_size_dist(r4vol,reff[ireff],veff[iveff]) else: n_r = logn_size_dist(r4vol, reff[ireff], veff[iveff]) Vol = np.trapz( n_r * 4./3.*np.pi*pow(r4vol,3.0), x=r4vol) / np.trapz( n_r, x=r4vol) Sext_o_WC[iwvl] = Cext[iwvl] / (Vol*rho_liq_wat) *1e-12 # microns2/g --> m2/g p11[:,iwvl] = p11_tmp[iwvl_tmp,:] p44[:,iwvl] = p44_tmp[iwvl_tmp,:] p21[:,iwvl] = p21_tmp[iwvl_tmp,:] p34[:,iwvl] = p34_tmp[iwvl_tmp,:] # check if angle def is the same for all wavelength and all particle size dist. if (nang!=len(ang_tmp)): print('') print('There is a problem with nang') print('') exit(0) for iang in range(nang): if ang_tmp[iang]!=ang[iang]: print('') print('There is a problem with ang') print('') exit(0) p11_integrate[iwvl] = np.trapz(p11[:,iwvl], x=np.cos(ang*np.pi/180.0)) for i_nbetal in range(n_nbetal): saveroot = "get_betal_"+ptcle_name+"_%i"%nbetal[i_nbetal] for iwvl in range(nwvl): alpha1_tmp, alpha2_tmp, alpha3_tmp, alpha4_tmp, beta1_tmp, beta2_tmp, nbetal_tmp, trunc_coef_tmp =\ ad.get_betal(ptcle_name, wvl[iwvl], dir_artdeco+'out/'+saveroot+'/Betal_'+ptcle_name+'_'+trunc_model+'.dat') if nbetal_tmp!=nbetal[i_nbetal]: print('') print('There is a problem with nbetal') print('') exit(0) alpha1[0:nbetal_tmp+1, i_nbetal, iwvl] = alpha1_tmp alpha2[0:nbetal_tmp+1, i_nbetal, iwvl] = alpha2_tmp alpha3[0:nbetal_tmp+1, i_nbetal, iwvl] = alpha3_tmp alpha4[0:nbetal_tmp+1, i_nbetal, iwvl] = alpha4_tmp beta1[0:nbetal_tmp+1, i_nbetal, iwvl] = beta1_tmp beta2[0:nbetal_tmp+1, i_nbetal, iwvl] = beta2_tmp trunc_coeff[i_nbetal,iwvl] = trunc_coef_tmp ang_tmp, p11_tmp, p22_tmp, p33_tmp, p44_tmp, p21_tmp, p34_tmp, \ Cext_tmp, ssa_tmp, g_tmp, \ nbetal_tmp, nteta_tmp, trunc_coef_tmp, sumf11_tmp = \ ad.get_opt(dir_opt, ptcle_name, wvl[iwvl], recomp = dir_artdeco+'out/'+saveroot+'/opt_recomp_'+ptcle_name+'_'+trunc_model+'.dat') # check if angle def is the same for all wavelength and all particle size dist. if (nang!=len(ang_tmp)): print('') print('There is a problem with nang') print('') exit(0) for iang in range(nang): if ang_tmp[iang]!=ang[iang]: print('') print('There is a problem with ang') print('') exit(0) if abs(sumf11_tmp - 2.00) > delta_sumf11 : print('sumf11_tmp != 2.0, ', ptcle_name, nbetal[i_nbetal], sumf11_tmp) p11_recomp[:, i_nbetal,iwvl] = p11_tmp p22_recomp[:, i_nbetal,iwvl] = p22_tmp p33_recomp[:, i_nbetal,iwvl] = p33_tmp p44_recomp[:, i_nbetal,iwvl] = p44_tmp p21_recomp[:, i_nbetal,iwvl] = p21_tmp p34_recomp[:, i_nbetal,iwvl] = p34_tmp p11_recomp_integrate[i_nbetal, iwvl] = np.trapz(p11_tmp, x= np.cos(ang*np.pi/180.0)) #print 'p11_recomp_integrate[i_nbetal,iwvl] = ', p11_recomp_integrate[i_nbetal, iwvl] ### Normalisation de Cext Cext_norm = np.zeros((nwvl)) indrefwvl = int(np.argwhere(wvl == 0.550)) Cext_norm[:] = Cext_norm[:] / Cext[indrefwvl] # re-interpolation des matrice de phase print(" interpolate phase matrices on lighter angular grid") p11_store = np.zeros((nang_store,nwvl)) p44_store = np.zeros((nang_store,nwvl)) p21_store = np.zeros((nang_store,nwvl)) p34_store = np.zeros((nang_store,nwvl)) p11_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) p22_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) p33_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) p44_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) p21_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) p34_recomp_store = np.zeros((nang_store,n_nbetal,nwvl)) for iwl in range(nwvl): fp11 = interp1d(ang,p11[:,iwl]) p11_store[:,iwl] = fp11(ang_store) fp44 = interp1d(ang,p44[:,iwl]) p44_store[:,iwl] = fp44(ang_store) fp21 = interp1d(ang,p21[:,iwl]) p21_store[:,iwl] = fp21(ang_store) fp34 = interp1d(ang,p34[:,iwl]) p34_store[:,iwl] = fp34(ang_store) for i_nbetal in range(n_nbetal): fp11 = interp1d(ang,p11_recomp[:,i_nbetal,iwl]) p11_recomp_store[:,i_nbetal,iwl] = fp11(ang_store) fp22 = interp1d(ang,p22_recomp[:,i_nbetal,iwl]) p22_recomp_store[:,i_nbetal,iwl] = fp22(ang_store) fp33 = interp1d(ang,p33_recomp[:,i_nbetal,iwl]) p33_recomp_store[:,i_nbetal,iwl] = fp33(ang_store) fp44 = interp1d(ang,p44_recomp[:,i_nbetal,iwl]) p44_recomp_store[:,i_nbetal,iwl] = fp44(ang_store) fp21 = interp1d(ang,p21_recomp[:,i_nbetal,iwl]) p21_recomp_store[:,i_nbetal,iwl] = fp21(ang_store) fp34 = interp1d(ang,p34_recomp[:,i_nbetal,iwl]) p34_recomp_store[:,i_nbetal,iwl] = fp34(ang_store) # suppress ARTDECO outputs for i_nbetal in range(n_nbetal): saveroot = "get_betal_"+ptcle_name+"_%i"%nbetal[i_nbetal] print(dir_artdeco+'out/'+saveroot) os.system("rm -rf "+dir_artdeco+'out/'+saveroot) # write result f = h5py.File(dir_opt+res_file, 'a') print("write result in hdf5 file") if 0.0 in f["liquid_cloud"]["axis"]["wavelengths"][:]: dset = f["liquid_cloud"]["axis"]["wavelengths"] dset[...] = wvl[:] dset = f["liquid_cloud"]["axis"]["reff"] dset[...] = reff[:] dset = f["liquid_cloud"]["axis"]["veff"] dset[...] = veff[:] dset = f["liquid_cloud"]["axis"]["phase_angles"] dset[...] = ang_store[:] dset = f["liquid_cloud"]["axis"]["mu"] dset[...] = mu_store[:] dset = f["liquid_cloud"]["axis"]["nbetal"] dset[...] = nbetal[:] dset = f["liquid_cloud"]["data"]["single_scattering_albedo"] dset[ireff,iveff,:] = ssa[:] dset = f["liquid_cloud"]["data"]["Cext"] dset[ireff,iveff,:] = Cext[:] dset = f["liquid_cloud"]["data"]["extinction_coefficient_over_wc"] dset[ireff,iveff,:] = Sext_o_WC[:] dset = f["liquid_cloud"]["data"]["integrate_p11"] dset[ireff,iveff,:] = p11_integrate[:] dset = f["liquid_cloud"]["data"]["integrate_p11_recomp"] dset[:,ireff,iveff,:] = p11_recomp_integrate[:,:] dset = f["liquid_cloud"]["data"]["p11_phase_function"] dset[:,ireff,iveff,:] = p11_store[:,:] dset = f["liquid_cloud"]["data"]["p21_phase_function"] dset[:,ireff,iveff,:] = p21_store[:,:] dset = f["liquid_cloud"]["data"]["p34_phase_function"] dset[:,ireff,iveff,:] = p34_store[:,:] dset = f["liquid_cloud"]["data"]["p44_phase_function"] dset[:,ireff,iveff,:] = p44_store[:,:] dset = f["liquid_cloud"]["data"]["truncation_coefficient"] dset[:,ireff,iveff,:] = trunc_coeff[:,:] dset = f["liquid_cloud"]["data"]["alpha1_betal"] dset[:,:,ireff,iveff,:] = alpha1[:,:,:] dset = f["liquid_cloud"]["data"]["alpha2_betal"] dset[:,:,ireff,iveff,:] = alpha2[:,:,:] dset = f["liquid_cloud"]["data"]["alpha3_betal"] dset[:,:,ireff,iveff,:] = alpha3[:,:,:] dset = f["liquid_cloud"]["data"]["alpha4_betal"] dset[:,:,ireff,iveff,:] = alpha4[:,:,:] dset = f["liquid_cloud"]["data"]["beta1_betal"] dset[:,:,ireff,iveff,:] = beta1[:,:,:] dset = f["liquid_cloud"]["data"]["beta2_betal"] dset[:,:,ireff,iveff,:] = beta2[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p11_phase_function"] dset[:,:,ireff,iveff,:] = p11_recomp_store[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p22_phase_function"] dset[:,:,ireff,iveff,:] = p22_recomp_store[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p33_phase_function"] dset[:,:,ireff,iveff,:] = p33_recomp_store[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p21_phase_function"] dset[:,:,ireff,iveff,:] = p21_recomp_store[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p34_phase_function"] dset[:,:,ireff,iveff,:] = p34_recomp_store[:,:,:] dset = f["liquid_cloud"]["data"]["recomposed_p44_phase_function"] dset[:,:,ireff,iveff,:] = p44_recomp_store[:,:,:] f.close() if __name__=='__main__': nbetal = np.array([8,16,32], dtype = int) baum_dir = '/rfs/data/baum_opt/' cdffile = 'GeneralHabitMixture_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' wvl = read_nc(baum_dir+cdffile,'wavelengths') wvl=wvl[1:] wvl = wvl[(wvl>0.5)&(wvl<20.0)] # nreff = 23 # min_reff = 5.0 # max_reff = 60.0 # dreff = (max_reff - min_reff) / (nreff-1) # reff = np.zeros(nreff) # for j in range(nreff): # reff[j] = min_reff + (j * dreff) nreff = 16 min_reff = 1.0 max_reff = 31.0 dreff = (max_reff - min_reff) / (nreff-1) reff = np.zeros(nreff) for j in range(nreff): reff[j] = min_reff + (j * dreff) nveff = 1 veff = np.array([0.11]) file_res="opt_gamma_liquid.h5" doit(wvl, reff, veff, nbetal, file_res, gamma=True, write_library=True) # ####### # # FORUM # nbetal = np.array([4, 8], dtype = int) # nwvl = 450 # wvl = np.logspace(np.log10(0.2), np.log10(100.0), num=nwvl, endpoint=True) # good = np.where( (wvl>=5.0) ) # wvl = wvl[good] # wvl = np.append(wvl,0.55) # wvl = np.sort(wvl) # nreff = 26 # min_reff = 5.0 # max_reff = 30.0 # dreff = (max_reff - min_reff) / (nreff-1) # reff = np.zeros(nreff) # for j in range(nreff): # reff[j] = min_reff + (j * dreff) # nveff = 1 # veff = np.array([0.11]) # file_res="opt_liquid_forum.h5" # doit(wvl, reff, veff, nbetal, file_res, gamma=True, write_library=True) # import matplotlib.pyplot as plt # reff = 30.0 # veff = 0.11 # r = np.linspace(1e-3, 2000., num=2000) # n_r = gamma_size_dist(r,reff,veff) # #n_r = logn_size_dist(r, reff, veff) # print(np.trapz( n_r, x=r) ) # print( np.trapz( n_r * 4./3.*np.pi*pow(r,3.0), x=r) / # np.trapz( n_r, x=r) ) # print( np.trapz( n_r * np.pi*pow(r,2.0), x=r) / # np.trapz( n_r, x=r) ) # print( np.trapz( n_r * pow(r,3.0), x=r) / # np.trapz( n_r * pow(r,2.0), x=r) ) # plt.plot(r, n_r) # plt.show()