#!/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/' dir_data = dir_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. ######################################################################################################### # aerosols only opac_aer_components = ['opac_inso', 'opac_miam', 'opac_micm', 'opac_minm', 'opac_mitr', 'opac_soot', 'opac_ssam00', 'opac_ssam50', 'opac_ssam70', 'opac_ssam80', 'opac_ssam90', 'opac_ssam95', 'opac_ssam98', 'opac_ssam99', 'opac_sscm00', 'opac_sscm50', 'opac_sscm70', 'opac_sscm80', 'opac_sscm90', 'opac_sscm95', 'opac_sscm98', 'opac_sscm99', 'opac_suso00', 'opac_suso50', 'opac_suso70', 'opac_suso80', 'opac_suso90', 'opac_suso95', 'opac_suso98', 'opac_suso99', 'opac_waso00', 'opac_waso50', 'opac_waso70', 'opac_waso80', 'opac_waso90', 'opac_waso95', 'opac_waso98', 'opac_waso99'] ################################################################################### def compute_aer_components(wvl, ptcle): nwvl = len(wvl) keywords = 'verbose opt_only' saveroot = 'opt_'+ptcle 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(nwvl): 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') # write artdeco_in.dat f = open(dir_artdeco+'input/artdeco_in_opt_'+ptcle+'.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_opt_"+ptcle+".dat") os.system("rm -rf "+dir_artdeco+'input/artdeco_in_opt_'+ptcle+'.dat') os.system("rm -rf "+dir_artdeco+'out/'+saveroot+'_unknown') return ptcle ################################################################################### def compute_mix_opt(wvl, fraction_in, ptcles, ptcle_name): print(" Compute properties for ", ptcle_name) nptcle = len(ptcles) nwvl = len(wvl) # normalization de fraction fraction = fraction_in / np.sum(fraction_in) # first compute all components optical properties Qext = np.zeros((nptcle, nwvl)) ssa = np.zeros((nptcle, nwvl)) g = np.zeros((nptcle, nwvl)) Qext_eff = np.zeros(nwvl) Qabs_eff = np.zeros(nwvl) Qsca_eff = np.zeros(nwvl) ssa_eff = np.zeros(nwvl) g_eff = np.zeros(nwvl) for iptcle in range(nptcle): print(" Read properties for ", ptcles[iptcle]) if not os.path.isfile(dir_opt+'/opt_'+ptcles[iptcle]+'.dat'): print('No file'+dir_opt+'/opt_'+ptcles[iptcle]+'.dat !!!') return wvl_tmp, ang_tmp, p11_tmp, p44_tmp, p21_tmp, p34_tmp, \ Qext_tmp, ssa_tmp, g_tmp = ad.get_opt_all(dir_opt, ptcles[iptcle]) # retrieve 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 (iptcle == 0) and (iwvl == 0) : nang = len(ang_tmp) ang = np.zeros(nang) p11 = np.zeros((nang, nptcle, nwvl)) p44 = np.zeros((nang, nptcle, nwvl)) p21 = np.zeros((nang, nptcle, nwvl)) p34 = np.zeros((nang, nptcle, nwvl)) p11_eff = np.zeros((nang, nwvl)) p44_eff = np.zeros((nang, nwvl)) p21_eff = np.zeros((nang, nwvl)) p34_eff = np.zeros((nang, nwvl)) ang[:] = ang_tmp Qext[iptcle,iwvl] = Qext_tmp[iwvl_tmp] ssa[iptcle,iwvl] = ssa_tmp[iwvl_tmp] g[iptcle,iwvl] = g_tmp[iwvl_tmp] p11[:,iptcle,iwvl] = p11_tmp[iwvl_tmp,:] p44[:,iptcle,iwvl] = p44_tmp[iwvl_tmp,:] p21[:,iptcle,iwvl] = p21_tmp[iwvl_tmp,:] p34[:,iptcle,iwvl] = p34_tmp[iwvl_tmp,:] Qsca = ssa * Qext Qabs = (1.0-ssa) * Qext print(" Compute it... ") for iwvl in range(nwvl): Qsca_eff[iwvl] = np.sum(Qsca[:,iwvl] * fraction) Qabs_eff[iwvl] = np.sum(Qabs[:,iwvl] * fraction) Qext_eff[iwvl] = np.sum(Qext[:,iwvl] * fraction) ssa_eff[iwvl] = Qsca_eff[iwvl] / Qext_eff[iwvl] frac_sca = fraction * Qsca[:,iwvl] frac_sca = frac_sca / np.sum(frac_sca) for iang in range(nang): p11_eff[iang,iwvl] = np.sum(frac_sca * p11[iang,:,iwvl]) p44_eff[iang,iwvl] = np.sum(frac_sca * p44[iang,:,iwvl]) p21_eff[iang,iwvl] = np.sum(frac_sca * p21[iang,:,iwvl]) p34_eff[iang,iwvl] = np.sum(frac_sca * p34[iang,:,iwvl]) #sumf11 = simps(p11_eff[:,iwvl], np.cos(ang*np.pi/180.0)) sumf11 = np.trapz(p11_eff[:,iwvl], x=np.cos(ang*np.pi/180.0)) if abs(sumf11 - 2.00) > delta_sumf11 : print('sumf 11 = ', sumf11, ' for wvl=', wvl[iwvl]) # Write opt_ file to be used for th betal expension fopt = open(dir_opt+'opt_'+ptcle_name+'.dat', 'w') fopt.write('# Optical properties to be used by ARTDECO \n') fopt.write('# \n') fopt.write('# Used model to obtain that properties is: \n') fopt.write(' unknown \n') fopt.write('# Used material is: \n') fopt.write(' unknown \n') fopt.write('# Number of phase matrix elements : \n') fopt.write(' 4 \n') fopt.write('# number of wavelengths \n') fopt.write('%i'%nwvl+' \n') fopt.write('# lambda(microns) nteta Cext (microns^2) SSA g \n') for iwvl in range(nwvl): s = ' %.7e'%wvl[iwvl]+' %i'%nang+' %.7e'%Qext_eff[iwvl]+' %.7e'%ssa_eff[iwvl]+' -32768 \n' fopt.write(s) fopt.write('# Phase matrix \n') for iwvl in range(nwvl): fopt.write('# lambda = %.4f'%wvl[iwvl]+' \n') fopt.write('# u F11 F44 F33 F21 F34 \n') for iang in range(nang): s = ' %16.8e'%(np.cos(ang[iang]*np.pi/180.))+' %15.7e'%p11_eff[iang, iwvl]+' %15.7e'%p44_eff[iang, iwvl]+' %15.7e'%p21_eff[iang, iwvl]+' %15.7e'%p34_eff[iang, iwvl]+' \n' fopt.write(s) fopt.close() return ################################################################################### def computeBetal(ptcle, wl, nbetal): keywords = ' betal_only print_betal print_recomp verbose' n_nbetal = len(nbetal) if not os.path.isfile(dir_opt+'/opt_'+ptcle+'.dat'): print('No file'+dir_opt+'/opt_'+ptcle+'.dat !!!') return for i_nbetal in range(n_nbetal): saveroot = "get_betal_"+ptcle+"_%.5f"%wl+"_%i"%nbetal[i_nbetal] 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('# artdeco_path \n') artdeco_in.append('/home/mathieu/work/RTM/artdeco/fortran/ \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') artdeco_in.append(' %.9f \n'%wl) 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[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 "+dir_artdeco+'input/artdeco_in_'+saveroot+'.dat') os.system("rm -rf "+dir_artdeco+'input/artdeco_in_'+saveroot+'.dat') return ################################################################################### def compute_components(wvl): wvl_list = [wvl]*len(opac_aer_components) if runlib_method == 'condor': pool = Pool(log=dir_runlib_log+'opt_prop', progressbar=True) pool.map(compute_aer_components, wvl_list, opac_aer_components) else: it = map(compute_aer_components, wvl_list, opac_aer_components) ################################################################################### def doit(aer, wvl, nbetal, res_file): ################################ # See OPAC (Hess, 1998 paper) if aer == 'cont_clean': ptcle_name = 'opac_cont_clean' fraction = np.array([ 2600, 0.15]) # this must be a fraction IN NUMBER #this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00', 'opac_inso'], ['opac_waso50', 'opac_inso'], ['opac_waso70', 'opac_inso'], ['opac_waso80', 'opac_inso'], ['opac_waso90', 'opac_inso'], ['opac_waso95', 'opac_inso'], ['opac_waso98', 'opac_inso'], ['opac_waso99', 'opac_inso']] elif aer == 'cont_avg': ptcle_name = 'opac_cont_avg' fraction = np.array([ 7000, 0.4, 8300]) # this must be a fraction IN NUMBER #this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00', 'opac_inso', 'opac_soot'], ['opac_waso50', 'opac_inso', 'opac_soot'], ['opac_waso70', 'opac_inso', 'opac_soot'], ['opac_waso80', 'opac_inso', 'opac_soot'], ['opac_waso90', 'opac_inso', 'opac_soot'], ['opac_waso95', 'opac_inso', 'opac_soot'], ['opac_waso98', 'opac_inso', 'opac_soot'], ['opac_waso99', 'opac_inso', 'opac_soot']] elif aer == 'cont_polluted': ptcle_name = 'opac_cont_polluted' fraction = np.array([ 15700, 0.6, 34300]) # this must be a fraction IN NUMBER #this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00', 'opac_inso', 'opac_soot'], ['opac_waso50', 'opac_inso', 'opac_soot'], ['opac_waso70', 'opac_inso', 'opac_soot'], ['opac_waso80', 'opac_inso', 'opac_soot'], ['opac_waso90', 'opac_inso', 'opac_soot'], ['opac_waso95', 'opac_inso', 'opac_soot'], ['opac_waso98', 'opac_inso', 'opac_soot'], ['opac_waso99', 'opac_inso', 'opac_soot']] elif aer == 'urban': ptcle_name = 'opac_urban' fraction = np.array([ 28000, 1.5, 130000]) # this must be a fraction IN NUMBER #this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00', 'opac_inso', 'opac_soot'], ['opac_waso50', 'opac_inso', 'opac_soot'], ['opac_waso70', 'opac_inso', 'opac_soot'], ['opac_waso80', 'opac_inso', 'opac_soot'], ['opac_waso90', 'opac_inso', 'opac_soot'], ['opac_waso95', 'opac_inso', 'opac_soot'], ['opac_waso98', 'opac_inso', 'opac_soot'], ['opac_waso99', 'opac_inso', 'opac_soot']] elif aer =='desert': ptcle_name = 'opac_desert' fraction = np.array([2000, 269.5, 30.5, 0.142]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00','opac_minm', 'opac_miam','opac_micm'], ['opac_waso50','opac_minm', 'opac_miam','opac_micm'], ['opac_waso70','opac_minm', 'opac_miam','opac_micm'], ['opac_waso80','opac_minm', 'opac_miam','opac_micm'], ['opac_waso90','opac_minm', 'opac_miam','opac_micm'], ['opac_waso95','opac_minm', 'opac_miam','opac_micm'], ['opac_waso98','opac_minm', 'opac_miam','opac_micm'], ['opac_waso99','opac_minm', 'opac_miam','opac_micm']] elif aer == 'maritime_clean': ptcle_name = 'opac_maritime_clean' fraction = np.array([1500, 20, 3.2e-3]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00','opac_ssam00','opac_sscm00'], ['opac_waso50','opac_ssam50','opac_sscm50'], ['opac_waso70','opac_ssam70','opac_sscm70'], ['opac_waso80','opac_ssam80','opac_sscm80'], ['opac_waso90','opac_ssam90','opac_sscm90'], ['opac_waso95','opac_ssam95','opac_sscm95'], ['opac_waso98','opac_ssam98','opac_sscm98'], ['opac_waso99','opac_ssam99','opac_sscm99']] elif aer == 'maritime_polluted': ptcle_name = 'opac_maritime_polluted' fraction = np.array([3800, 20, 3.2e-3, 5180]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00','opac_ssam00','opac_sscm00','opac_soot'], ['opac_waso50','opac_ssam50','opac_sscm50','opac_soot'], ['opac_waso70','opac_ssam70','opac_sscm70','opac_soot'], ['opac_waso80','opac_ssam80','opac_sscm80','opac_soot'], ['opac_waso90','opac_ssam90','opac_sscm90','opac_soot'], ['opac_waso95','opac_ssam95','opac_sscm95','opac_soot'], ['opac_waso98','opac_ssam98','opac_sscm98','opac_soot'], ['opac_waso99','opac_ssam99','opac_sscm99','opac_soot']] elif aer == 'maritime_tropical': ptcle_name = 'opac_maritime_tropical' fraction = np.array([590, 10, 1.3e-3]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00','opac_ssam00','opac_sscm00'], ['opac_waso50','opac_ssam50','opac_sscm50'], ['opac_waso70','opac_ssam70','opac_sscm70'], ['opac_waso80','opac_ssam80','opac_sscm80'], ['opac_waso90','opac_ssam90','opac_sscm90'], ['opac_waso95','opac_ssam95','opac_sscm95'], ['opac_waso98','opac_ssam98','opac_sscm98'], ['opac_waso99','opac_ssam99','opac_sscm99']] elif aer == 'arctic': ptcle_name = 'opac_arctic' fraction = np.array([1300, 0.01, 1.9, 5300]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00','opac_inso','opac_ssam00','opac_soot'], ['opac_waso50','opac_inso','opac_ssam50','opac_soot'], ['opac_waso70','opac_inso','opac_ssam70','opac_soot'], ['opac_waso80','opac_inso','opac_ssam80','opac_soot'], ['opac_waso90','opac_inso','opac_ssam90','opac_soot'], ['opac_waso95','opac_inso','opac_ssam95','opac_soot'], ['opac_waso98','opac_inso','opac_ssam98','opac_soot'], ['opac_waso99','opac_inso','opac_ssam99','opac_soot']] elif aer == 'antarctic': ptcle_name = 'opac_antarctic' fraction = np.array([42.9, 0.47e-1, 0.53e-2]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_suso00','opac_ssam00','opac_mitr'], ['opac_suso50','opac_ssam50','opac_mitr'], ['opac_suso70','opac_ssam70','opac_mitr'], ['opac_suso80','opac_ssam80','opac_mitr'], ['opac_suso90','opac_ssam90','opac_mitr'], ['opac_suso95','opac_ssam95','opac_mitr'], ['opac_suso98','opac_ssam98','opac_mitr'], ['opac_suso99','opac_ssam99','opac_mitr']] elif aer == "soot": ptcle_name = 'opac_component_soot' fraction = np.array([100]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_soot'], ['opac_soot'], ['opac_soot'], ['opac_soot'], ['opac_soot'], ['opac_soot'], ['opac_soot'], ['opac_soot']] elif aer == "waso": ptcle_name = 'opac_component_waso' fraction = np.array([100]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the code hygro = np.array([0,50,70,80,90,95,98,99]) ptcles = [['opac_waso00'], ['opac_waso50'], ['opac_waso70'], ['opac_waso80'], ['opac_waso90'], ['opac_waso95'], ['opac_waso98'], ['opac_waso99']] # check whether the opac particle is already present # in library file filename = dir_data+res_file if os.path.isfile(filename): f = h5py.File(filename, 'a') grpkeys = f.keys() f.close if ptcle_name in grpkeys: print(ptcle_name+' already present in base') return print("Compute opt...") params = [] for ihy in range(len(hygro)): params.append((wvl, fraction, ptcles[:][ihy], ptcle_name+'%02i'%hygro[ihy])) print(" Number of compute_mix_opt() call is", len(params)) if runlib_method == 'condor': pool = Pool(log=dir_runlib_log+'compute_mix_opt', progressbar=True, prio = -5) it = pool.map(compute_mix_opt, *zip(*params)) else: it = map(compute_mix_opt, *zip(*params)) for ii in it: print(ii) # compute BETAL print('Compute Legendre poly. exp. ...') params = [] for (hy, wl) in product(hygro,wvl): params.append((ptcle_name+'%02i'%hy, wl, nbetal)) print(" Number of computeBetal() call is", len(params)) if runlib_method == 'condor': pool = Pool(log=dir_runlib_log+'computeBetal', progressbar=True, prio = -5) it = pool.map(computeBetal, *zip(*params)) else: it = map(computeBetal, *zip(*params)) for ii in it: print(ii) # Retrieve OPT print("Retrieve opt...") nwvl = len(wvl) nhygro = len(hygro) g = np.zeros((nhygro,nwvl)) ssa = np.zeros((nhygro,nwvl)) Cext = np.zeros((nhygro,nwvl)) p11_integrate = np.zeros((nhygro, nwvl)) for ihy in range(nhygro): ptcle_tmp = ptcle_name+'%02i'%hygro[ihy] 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_tmp) # retrieve 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 (ihy==0) & (iwvl==0): nang = len(ang_tmp) ang = np.zeros(nang) p11 = np.zeros((nang, nhygro, nwvl)) p44 = np.zeros((nang, nhygro, nwvl)) p21 = np.zeros((nang, nhygro, nwvl)) p34 = np.zeros((nang, nhygro, nwvl)) ang[:] = ang_tmp Cext[ihy,iwvl] = Cext_tmp[iwvl_tmp] ssa[ihy,iwvl] = ssa_tmp[iwvl_tmp] g[ihy,iwvl] = g_tmp[iwvl_tmp] p11[:,ihy,iwvl] = p11_tmp[iwvl_tmp,:] p44[:,ihy,iwvl] = p44_tmp[iwvl_tmp,:] p21[:,ihy,iwvl] = p21_tmp[iwvl_tmp,:] p34[:,ihy,iwvl] = p34_tmp[iwvl_tmp,:] #p11_integrate[ihy,iwvl] = simps(p11[:,ihy,iwvl], np.cos(ang*np.pi/180.0)) p11_integrate[ihy,iwvl] = np.trapz(p11[:,ihy,iwvl], x=np.cos(ang*np.pi/180.0)) # Get BETAL print("Retrieve Betal...") n_nbetal = len(nbetal) nbetalmax = np.amax(nbetal) trunc_coeff = np.zeros((n_nbetal, nhygro, nwvl)) alpha1 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) alpha2 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) alpha3 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) alpha4 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) beta1 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) beta2 = np.zeros((nbetalmax+1, n_nbetal, nhygro, nwvl)) p11_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p22_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p33_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p44_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p21_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p34_recomp = np.zeros((nang, n_nbetal, nhygro, nwvl)) p11_recomp_integrate = np.zeros((n_nbetal, nhygro, nwvl)) for ihy in range(nhygro): ptcle_tmp = ptcle_name+'%02i'%hygro[ihy] for iwvl in range(nwvl): for i_nbetal in range(n_nbetal): saveroot = "get_betal_"+ptcle_tmp+"_%.5f"%wvl[iwvl]+"_%i"%nbetal[i_nbetal] alpha1_tmp, alpha2_tmp, alpha3_tmp, alpha4_tmp, beta1_tmp, beta2_tmp, nbetal_tmp, trunc_coef_tmp =\ ad.get_betal(ptcle_tmp, wvl[iwvl], dir_artdeco+'out/'+saveroot+'/Betal_'+ptcle_tmp+'_'+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, ihy, iwvl] = alpha1_tmp alpha2[0:nbetal_tmp+1, i_nbetal, ihy, iwvl] = alpha2_tmp alpha3[0:nbetal_tmp+1, i_nbetal, ihy, iwvl] = alpha3_tmp alpha4[0:nbetal_tmp+1, i_nbetal, ihy, iwvl] = alpha4_tmp beta1[0:nbetal_tmp+1, i_nbetal, ihy, iwvl] = beta1_tmp beta2[0:nbetal_tmp+1, i_nbetal, ihy, iwvl] = beta2_tmp trunc_coeff[i_nbetal, ihy, 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_tmp, wvl[iwvl], recomp = dir_artdeco+'out/'+saveroot+'/opt_recomp_'+ptcle_tmp+'_'+trunc_model+'.dat') os.system("rm -rf "+dir_artdeco+'out/'+saveroot) # 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_tmp, nbetal[i_nbetal], sumf11_tmp) p11_recomp[:, i_nbetal, ihy, iwvl] = p11_tmp p22_recomp[:, i_nbetal, ihy, iwvl] = p22_tmp p33_recomp[:, i_nbetal, ihy, iwvl] = p33_tmp p44_recomp[:, i_nbetal, ihy, iwvl] = p44_tmp p21_recomp[:, i_nbetal, ihy, iwvl] = p21_tmp p34_recomp[:, i_nbetal, ihy, iwvl] = p34_tmp #p11_recomp_integrate[i_nbetal,ihy, iwvl] = simps(p11_tmp, np.cos(ang*np.pi/180.0)) p11_recomp_integrate[i_nbetal,ihy, iwvl] = np.trapz(p11_tmp, x=np.cos(ang*np.pi/180.0)) #print 'p11_recomp_integrate[i_nbetal,ihy,iwvl] = ', p11_recomp_integrate[i_nbetal,ihy, iwvl] # re-interpolation des matrice de phase print("Interpolate phase matrices on lighter angular grid") p11_store = np.zeros((nang_store,nhygro,nwvl)) p44_store = np.zeros((nang_store,nhygro,nwvl)) p21_store = np.zeros((nang_store,nhygro,nwvl)) p34_store = np.zeros((nang_store,nhygro,nwvl)) p11_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) p22_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) p33_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) p44_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) p21_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) p34_recomp_store = np.zeros((nang_store,n_nbetal,nhygro,nwvl)) for iii in range(nhygro): for iwl in range(nwvl): fp11 = interp1d(ang,p11[:,iii,iwl]) p11_store[:,iii,iwl] = fp11(ang_store) fp44 = interp1d(ang,p44[:,iii,iwl]) p44_store[:,iii,iwl] = fp44(ang_store) fp21 = interp1d(ang,p21[:,iii,iwl]) p21_store[:,iii,iwl] = fp21(ang_store) fp34 = interp1d(ang,p34[:,iii,iwl]) p34_store[:,iii,iwl] = fp34(ang_store) for i_nbetal in range(n_nbetal): fp11 = interp1d(ang,p11_recomp[:,i_nbetal,iii,iwl]) p11_recomp_store[:,i_nbetal,iii,iwl] = fp11(ang_store) fp22 = interp1d(ang,p22_recomp[:,i_nbetal,iii,iwl]) p22_recomp_store[:,i_nbetal,iii,iwl] = fp22(ang_store) fp33 = interp1d(ang,p33_recomp[:,i_nbetal,iii,iwl]) p33_recomp_store[:,i_nbetal,iii,iwl] = fp33(ang_store) fp44 = interp1d(ang,p44_recomp[:,i_nbetal,iii,iwl]) p44_recomp_store[:,i_nbetal,iii,iwl] = fp44(ang_store) fp21 = interp1d(ang,p21_recomp[:,i_nbetal,iii,iwl]) p21_recomp_store[:,i_nbetal,iii,iwl] = fp21(ang_store) fp34 = interp1d(ang,p34_recomp[:,i_nbetal,iii,iwl]) p34_recomp_store[:,i_nbetal,iii,iwl] = fp34(ang_store) # Normalisation de Cext Cext_norm = np.zeros_like(Cext) indrefwvl = int(np.argwhere(wvl == 0.550)) for iii in range(nhygro): Cext_norm[iii,:] = Cext[iii,:] / Cext[iii,indrefwvl] ############################ # WRITE A RESULT FILE (HDF5) f = h5py.File(filename, 'a') grp = f.require_group(ptcle_name) subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (nwvl,), dtype='float64') dset[...] = wvl dset = subgrp.require_dataset("humidity", (nhygro,), dtype='float64') dset[...] = hygro dset = subgrp.require_dataset("phase_angles", (nang_store,), dtype='float64') dset[...] = ang_store dset = subgrp.require_dataset("mu", (nang_store,), dtype='float64') dset[...] = np.cos(ang_store*np.pi/180.0) dset = subgrp.require_dataset("nbetal", (n_nbetal,), dtype='i') dset[...] = nbetal 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", (nhygro,nwvl), dtype='float64') dset[...] = ssa dset.attrs.create("dimensions", np.string_("humidity,wavelengths")) dset = subgrp.require_dataset("Cext", (nhygro,nwvl), dtype='float64') dset[...] = Cext dset.attrs.create("dimensions", np.string_("humidity,wavelengths")) #dset.attrs.create("unit", np.string_("micron^2")) # dset = subgrp.require_dataset("normed_ext_coeff", (nhygro,nwvl), dtype='float64') # dset[...] = Cext_norm # dset.attrs.create("dimensions", np.string_("humidity,wavelengths") dset = subgrp.require_dataset("integrate_p11", (nhygro,nwvl), dtype='float64') dset[...] = p11_integrate dset.attrs.create("dimensions", np.string_("humidity,wavelengths")) dset = subgrp.require_dataset("p11_phase_function", (nang_store,nhygro,nwvl), dtype='float64') dset[...] = p11_store dset.attrs.create("dimensions", np.string_("mu,humidity,wavelengths")) dset = subgrp.require_dataset("p21_phase_function", (nang_store,nhygro,nwvl), dtype='float64') dset[...] = p21_store dset.attrs.create("dimensions", np.string_("mu,humidity,wavelengths")) dset = subgrp.require_dataset("p34_phase_function", (nang_store,nhygro,nwvl), dtype='float64') dset[...] = p34_store dset.attrs.create("dimensions", np.string_("mu,humidity,wavelengths")) dset = subgrp.require_dataset("p44_phase_function", (nang_store,nhygro,nwvl), dtype='float64') dset[...] = p44_store dset.attrs.create("dimensions", np.string_("mu,humidity,wavelengths")) dset = subgrp.require_dataset("truncation_coefficient", (n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = trunc_coeff dset.attrs.create("dimensions", np.string_("nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("integrate_p11_recomp", (n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p11_recomp_integrate dset.attrs.create("dimensions", np.string_("nbetal,humidity,wavelengths")) dset = subgrp.require_dataset("alpha1_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = alpha1 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha2_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = alpha2 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha3_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = alpha3 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("alpha4_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = alpha4 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("beta1_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = beta1 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("beta2_betal", (nbetalmax+1, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = beta2 dset.attrs.create("dimensions", np.string_("ibetal,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p11_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p11_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p22_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p22_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p33_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p33_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p44_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p44_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p21_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p21_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset = subgrp.require_dataset("recomposed_p34_phase_function", (nang_store, n_nbetal,nhygro,nwvl), dtype='float64') dset[...] = p34_recomp_store dset.attrs.create("dimensions", np.string_("mu,nbetal,humidity,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) f.close return if __name__=='__main__': nwvl = 450 wvl = np.logspace(np.log10(0.2), np.log10(100.0), num=nwvl, endpoint=True) good = np.where( (wvl>0.25) & (wvl<40.0) ) wvl = wvl[good] # good = np.where( (wvl>=1.0) ) # wvl = wvl[good] wvl = np.append(wvl,0.55) wvl = np.sort(wvl) print(wvl) print(len(wvl)) #compute_components(wvl) nbetal = np.array([8, 16, 32, 64], dtype = int) # res_file = "opt_opac.h5" # doit('cont_clean', wvl, nbetal, res_file) # doit('cont_avg', wvl, nbetal, res_file) # doit('cont_polluted', wvl, nbetal, res_file) # doit('urban', wvl, nbetal, res_file) # doit('desert', wvl, nbetal, res_file) # doit('maritime_clean', wvl, nbetal, res_file) # doit('maritime_polluted', wvl, nbetal, res_file) # doit('maritime_tropical', wvl, nbetal, res_file) # doit('arctic', wvl, nbetal, res_file) # doit('antarctic', wvl, nbetal, res_file) res_file = "opt_opac_component.h5" doit('soot', wvl, nbetal, res_file) doit('waso', wvl, nbetal, res_file)