#!/usr/bin/env python # This was developped by Mathieu Compiegne at HYGEOS # All routines in this file act # on the artdeco f2py object from __future__ import print_function, division, absolute_import from builtins import range import sys import os import string import numpy as np import h5py import collections import time compare = lambda x, y: collections.Counter(x) == collections.Counter(y) ################################################### def run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=None, verbose=False, save_layers=None, get_trans_atm=False, rho_surf_output=False): ''' save_layers is either None or a directory PATH and ROOT name save_layers = {"PATH":string, "ROOT":string} ''' import artdeco if verbose: print("") print("====================================") print("") print("") print("") print(" Enter run_artdeco ") print("") print("") print("") # set artdeco variables t0 = time.time() set_artdeco_var(artdeco, artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=verbose, save_layers=save_layers) t_set_var = time.time()-t0 if (artdeco_in.mode == 'kdis'): if len(atmos.gasabs) < 1: print("\n (run_artdeco) No absorbing gas although you are in Kdis mode !") exit() # adjust the KDIS gases to those of the atmosphere select_kdis_gas(atmos.gasabs, artdeco) if (artdeco.mcommon.nptcle > 0): t0 = time.time() # if not already available, we need to compute Betal artdeco.mget_betal.get_betal() if verbose: print("\n (run_artdeco) Get_betal time =", time.time()-t0) t0 = time.time() # set the atmosphere optical properties artdeco.mlayers.get_layers() t_get_layer = time.time()-t0 if save_layers: save_smartg_layers(artdeco, save_layers) # this unallocate all ARTDECO arrays deallocate(artdeco) return if get_trans_atm: n_call_rtes = np.zeros(artdeco.mcommon.nlambda, dtype=int) trans_mol_atm = np.zeros(artdeco.mcommon.nlambda) if verbose: print("\n (run_artdeco) Compute atm gas transmission...") t0 = time.time() for j in range(artdeco.mcommon.nlambda): artdeco.mlayers.get_mono_layers(j+1) # +1 cause this index is in fortran n_call_rtes[j] = artdeco.mcommon.n_aitaui_mono_layers trans_mol_atm[j] = np.sum( artdeco.mcommon.ai_mono_layers[0:artdeco.mcommon.n_aitaui_mono_layers] * np.exp( - np.sum( artdeco.mcommon.gas_dtau_mono_layers[:,0:artdeco.mcommon.n_aitaui_mono_layers] * (1.0 - artdeco.mcommon.gas_ssa_mono_layers[:,0:artdeco.mcommon.n_aitaui_mono_layers]), axis=0 ))) lamb=np.copy(artdeco.mcommon.wlambda) deallocate(artdeco) if verbose: print("\n (run_artdeco) Computed atm gas transmission in %.2e seconds"%(time.time()-t0)) return lamb, trans_mol_atm, n_call_rtes t0 = time.time() if artdeco_in.rt_model == 'doad': # call DOAD artdeco.mcall_doad.call_doad() elif artdeco_in.rt_model == 'disort': # call DISORT artdeco.mcall_od.call_od() elif artdeco_in.rt_model == 'sinsca': artdeco.msinsca.call_sinsca() if verbose: print("\n (run_artdeco) SetVar time = ", t_set_var) print("\n (run_artdeco) GetLayers time = ", t_get_layer) print("\n (run_artdeco) RT time = ", time.time()-t0) # returned values flux = np.zeros((artdeco.mcommon.nsza, 3, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) alt_flux = np.zeros(artdeco.mcommon.nalt_atm) lamb = np.zeros(artdeco.mcommon.nlambda) rho_surf = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nlambda)) radiance = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) radiance_levels = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) flux[:,:,:,:] = artdeco.mcommon.flux_out[:,:,:,:] alt_flux[:] = artdeco.mcommon.alt_atm[:] lamb[:] = artdeco.mcommon.wlambda[:] radiance[:,:,:,:,:] = artdeco.mcommon.svr[:,:,:,:,:] radiance_levels[:,:,:,:,:,:] = artdeco.mcommon.svr_levels[:,:,:,:,:,:] rho_surf[:,:,:,:] = artdeco.mcommon.rho_surf[:,:,:,:] # this unallocate all ARTDECO arrays deallocate(artdeco) if verbose: print("") print("") print("") print(" Done run_artdeco ") print("") print("") print("") print("====================================") print("") if rho_surf_output: return lamb, radiance, radiance_levels, flux, alt_flux, rho_surf else: return lamb, radiance, radiance_levels, flux, alt_flux ################################################################################################# def save_smartg_layers(artdeco, save_layers): print("(save_smartg_layers) Save layers for SMART-G") if get_char(artdeco.mcommon.mode) == "kdis": print("") print("(save_smartg_layers) ERROR ") print(" No possible save if ARTDECO used in kdis mode") print("") return # Save the gas absorption and diffusion opacity as a function of wavelength # and altitude cumul_tau_ray = np.zeros(( artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda )) cumul_tau_mol = np.zeros(( artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda )) cumul_tau_ray[:,:] = np.cumsum( np.append( np.zeros((1, artdeco.mcommon.nlambda)), artdeco.mcommon.gas_dtausca_layers[:,:], axis=0), axis = 0) cumul_tau_mol[:,:] = np.cumsum( np.append( np.zeros((1, artdeco.mcommon.nlambda)), artdeco.mcommon.gas_dtauabs_layers[:,0,0,:], axis=0), axis = 0) f = open(save_layers["PATH"]+'cTauRay_'+save_layers["ROOT"]+'.dat', 'w') f.write('# Number of wavelength \n') f.write('%i \n'%artdeco.mcommon.nlambda) f.write('# Number of altitude \n') f.write('%i \n'% artdeco.mcommon.nalt_atm) s ='' for i in range(artdeco.mcommon.nlambda): s = s + ' %.6f'%(artdeco.mcommon.wlambda[i]*1.e3) f.write('# Wavelengths\n ') f.write(s + '\n') f.write('# OD \n') for ialt in range(artdeco.mcommon.nalt_atm): s = ' %.2f'%artdeco.mcommon.alt_atm[ialt] for i in range(artdeco.mcommon.nlambda): s = s + ' %.6e'%cumul_tau_ray[ialt, i] f.write(s + '\n') f.close() f = open(save_layers["PATH"]+'cTauGas_'+save_layers["ROOT"]+'.dat', 'w') f.write('# Number of wavelength \n') f.write('%i \n'%artdeco.mcommon.nlambda) f.write('# Number of altitude \n') f.write('%i \n'% artdeco.mcommon.nalt_atm) s ='' for i in range(artdeco.mcommon.nlambda): s = s + ' %.6f'%(artdeco.mcommon.wlambda[i]*1.e3) f.write('# Wavelengths\n ') f.write(s + '\n') f.write('# OD \n') for ialt in range(artdeco.mcommon.nalt_atm): s = ' %.2f'%artdeco.mcommon.alt_atm[ialt] for i in range(artdeco.mcommon.nlambda): s = s + ' %.6e'%cumul_tau_mol[ialt, i] f.write(s + '\n') f.close() if artdeco.mcommon.nptcle > 0: # If particles in the atmosphere # Save the scattering and absorption opacity as a function # of wavelength and altitude dtau_sca_ptcle = np.zeros(( artdeco.mcommon.nlayers, artdeco.mcommon.nlambda )) dtau_abs_ptcle = np.zeros(( artdeco.mcommon.nlayers, artdeco.mcommon.nlambda )) dtau_abs_ptcle[:,:] = artdeco.mcommon.ptcle_dtau_layers * ( 1.0 - artdeco.mcommon.ptcle_ssa_layers) dtau_sca_ptcle[:,:] = artdeco.mcommon.ptcle_dtau_layers * ( artdeco.mcommon.ptcle_ssa_layers) cumul_tau_sca_ptcle = np.cumsum( np.append( np.zeros((1, artdeco.mcommon.nlambda)), dtau_sca_ptcle[:,:], axis=0), axis = 0) cumul_tau_abs_ptcle = np.cumsum( np.append( np.zeros((1, artdeco.mcommon.nlambda)), dtau_abs_ptcle[:,:], axis=0), axis = 0) f = open(save_layers["PATH"]+'cTauAbs_ptcle_'+save_layers["ROOT"]+'.dat', 'w') f.write('# Number of wavelength \n') f.write('%i \n'%artdeco.mcommon.nlambda) f.write('# Number of altitude \n') f.write('%i \n'% artdeco.mcommon.nalt_atm) s ='' for i in range(artdeco.mcommon.nlambda): s = s + ' %.6f'%(artdeco.mcommon.wlambda[i]*1.e3) f.write('# Wavelengths\n ') f.write(s + '\n') f.write('# OD \n') for ialt in range(artdeco.mcommon.nalt_atm): s = ' %.2f'%artdeco.mcommon.alt_atm[ialt] for i in range(artdeco.mcommon.nlambda): s = s + ' %.6e'%cumul_tau_abs_ptcle[ialt, i] f.write(s + '\n') f.close() f = open(save_layers["PATH"]+'cTauSca_ptcle_'+save_layers["ROOT"]+'.dat', 'w') f.write('# Number of wavelength \n') f.write('%i \n'%artdeco.mcommon.nlambda) f.write('# Number of altitude \n') f.write('%i \n'% artdeco.mcommon.nalt_atm) s ='' for i in range(artdeco.mcommon.nlambda): s = s + ' %.6f'%(artdeco.mcommon.wlambda[i]*1.e3) f.write('# Wavelengths\n ') f.write(s + '\n') f.write('# OD \n') for ialt in range(artdeco.mcommon.nalt_atm): s = ' %.2f'%artdeco.mcommon.alt_atm[ialt] for i in range(artdeco.mcommon.nlambda): s = s + ' %.6e'%cumul_tau_sca_ptcle[ialt, i] f.write(s + '\n') f.close() nang = artdeco.mcommon.ptcle_nang_phasemat_layers[0] # We also save a phase matrix # ptcle_phasemat_layers(nlayers, 6, nang_max,nlambda) # ptcle_u_phasemat_layers(nang_max, nlambda) # ptcle_nang_phasemat_layers(nlambda) if len(np.unique(artdeco.mcommon.ptcle_nang_phasemat_layers)) != 1: print("(save_smartg_layers) ERROR") print(" Phase matrix angular grid must be common over wavelengths") print(" ") exit() for iang in range(nang): if len(np.unique(artdeco.mcommon.ptcle_u_phasemat_layers[iang, :])) != 1: print("(save_smartg_layers) ERROR") print(" Phase matrix angular grid must be common over wavelengths") print(" ") exit() # for iang in range(nang): # for iwvl in range(artdeco.mcommon.nlambda): # if len(np.unique(artdeco.mcommon.ptcle_phasemat_layers[:,0,iang, iwvl])) != 1: # print "(save_smartg_layers) ERROR" # print " Phase matrix should be constant over layers " # print " " # exit() for iang in range(nang): iwvl = 0 if len(np.unique(artdeco.mcommon.ptcle_phasemat_layers[:,0,iang, iwvl])) != 1: print("(save_smartg_layers) ERROR") print(" Phase matrix should be constant over layers ") print(" ") exit() iwvl = int(artdeco.mcommon.nlambda / 2) if len(np.unique(artdeco.mcommon.ptcle_phasemat_layers[:,0,iang, iwvl])) != 1: print("(save_smartg_layers) ERROR") print(" Phase matrix should be constant over layers ") print(" ") exit() iwvl = artdeco.mcommon.nlambda - 1 if len(np.unique(artdeco.mcommon.ptcle_phasemat_layers[:,0,iang, iwvl])) != 1: print("(save_smartg_layers) ERROR") print(" Phase matrix should be constant over layers ") print(" ") exit() for iptcle in range(artdeco.mcommon.nptcle): if artdeco.mcommon.ptcle_nelem[iptcle] != 4: print("(save_smartg_layers) ERROR") print(" Used particules should all have nelem = 4") print(" ") exit() # F11 # F22 # F33 # F44 # F21 # F34 filename = save_layers["PATH"]+'phasemat_'+save_layers["ROOT"]+'.dat' f = open(filename, 'w') # SSA in SMARTG f.write(' %i'%nang+' \n') # select a layer with ptcles ilay = np.argwhere(artdeco.mcommon.ptcle_dtau_layers[:,0] > 0.0 )[0][0] iwvl = 0 f.write(' %f \n'%(artdeco.mcommon.wlambda[iwvl]*1e3)) f.write(' %f \n'%artdeco.mcommon.ptcle_ssa_layers[ilay,iwvl] ) M2 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] - artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M1 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] + artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M3 = artdeco.mcommon.ptcle_phasemat_layers[ilay,2,0:nang, iwvl] ang = np.arccos(artdeco.mcommon.ptcle_u_phasemat_layers[0:nang, iwvl]) * 180.0 / np.pi for i in range(nang): s = ' %e'%ang[i]+' %e'%M1[i]+' %e'%M2[i]+' %e'%M3[i]+' %e'%0.0 f.write(s+' \n') iwvl = int(artdeco.mcommon.nlambda / 2) f.write(' %f \n'%(artdeco.mcommon.wlambda[iwvl]*1e3)) f.write(' %f \n'%artdeco.mcommon.ptcle_ssa_layers[ilay,iwvl] ) M2 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] - artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M1 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] + artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M3 = artdeco.mcommon.ptcle_phasemat_layers[ilay,2,0:nang, iwvl] ang = np.arccos(artdeco.mcommon.ptcle_u_phasemat_layers[0:nang, iwvl]) * 180.0 / np.pi for i in range(nang): s = ' %e'%ang[i]+' %e'%M1[i]+' %e'%M2[i]+' %e'%M3[i]+' %e'%0.0 f.write(s+' \n') iwvl = artdeco.mcommon.nlambda - 1 f.write(' %f \n'%(artdeco.mcommon.wlambda[iwvl]*1e3)) f.write(' %f \n'%artdeco.mcommon.ptcle_ssa_layers[ilay,iwvl] ) M2 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] - artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M1 = artdeco.mcommon.ptcle_phasemat_layers[ilay,0,0:nang, iwvl] + artdeco.mcommon.ptcle_phasemat_layers[ilay,4,0:nang, iwvl] M3 = artdeco.mcommon.ptcle_phasemat_layers[ilay,2,0:nang, iwvl] ang = np.arccos(artdeco.mcommon.ptcle_u_phasemat_layers[0:nang, iwvl]) * 180.0 / np.pi for i in range(nang): s = ' %e'%ang[i]+' %e'%M1[i]+' %e'%M2[i]+' %e'%M3[i]+' %e'%0.0 f.write(s+' \n') f.write("# \n") f.close() ################################################################################################# def set_artdeco_var(artdeco, artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=None, verbose=False, save_layers=None): ''' ''' if verbose: print("(set_artdeco_var) start...") artdeco.mcommon.f2py_mode = True # init some values artdeco.mcommon.verbose = False artdeco.mcommon.opt_only = False artdeco.mcommon.betal_only = False artdeco.mcommon.thermal = False artdeco.mcommon.do_rt = True artdeco.mcommon.sinsca_corint = False artdeco.mcommon.sinsca_corint_brdf = False artdeco.mcommon.warning = True artdeco.mcommon.no_rayleigh = False artdeco.mcommon.thermal_only = False artdeco.mcommon.flux_only = False artdeco.mcommon.print_betal = False artdeco.mcommon.print_recomp = False artdeco.mcommon.print_aitaui = False artdeco.mcommon.od_no_check = False artdeco.mcommon.is_secsca = False artdeco.mcommon.bottom_f0 = False artdeco.mcommon.od_deltam = False artdeco.mcommon.od_corint = False artdeco.mcommon.od_do_secsca = True artdeco.mcommon.od_use_fisot = False artdeco.mcommon.flag_user_vdist = False set_char(artdeco.mcommon.mode, 'undef') set_char(artdeco.mcommon.trunc_method, 'undef') set_char(artdeco.mcommon.rt_model, 'undef') artdeco.mcommon.nmat = artdeco.mconstants.undef_i artdeco.mcommon.hg_flag = False artdeco.mcommon.flag_mie = False artdeco.mcommon.flag_hm = False artdeco.mcommon.flag_baum = False artdeco.mcommon.nptcle = 0 artdeco.mcommon.nbetal_in = -1 artdeco.mcommon.nmaxbetal = -1 artdeco.mcommon.lat = artdeco.mconstants.lat_ray_default artdeco.mcommon.lon = -1 if 'verbose' in artdeco_in.keywords: artdeco.mcommon.verbose = True if verbose: artdeco.mcommon.verbose = True if 'opt_only' in artdeco_in.keywords: print(" (set_artdeco_var) ERROR") print(" opt_only not allowed") exit() if 'betal_only' in artdeco_in.keywords: print(" (set_artdeco_var) ERROR") print(" betal_only not allowed") exit() if 'nowarn' in artdeco_in.keywords: artdeco.mcommon.warning = False if 'no_rayleigh' in artdeco_in.keywords: artdeco.mcommon.no_rayleigh = True if 'thermal_only' in artdeco_in.keywords: artdeco.mcommon.thermal_only = True if 'flux_only' in artdeco_in.keywords: artdeco.mcommon.flux_only = True if 'print_betal' in artdeco_in.keywords: artdeco.mcommon.print_betal = True if 'print_recomp' in artdeco_in.keywords: artdeco.mcommon.print_recomp = True if 'print_aitaui' in artdeco_in.keywords: artdeco.mcommon.print_aitaui = True if 'od_no_check' in artdeco_in.keywords: artdeco.mcommon.od_no_check = True if 'bottom_f0' in artdeco_in.keywords: artdeco.mcommon.bottom_f0 = True if artdeco_in.rt_model != 'disort': print(" (set_artdeco_var) ERROR") print(" with -bottom_f0- ") print(" rt_model must be -disort-") exit() if artdeco.mcommon.betal_only or artdeco.mcommon.opt_only: artdeco.mcommon.do_rt = False # mode if artdeco_in.mode not in ['mono', 'kdis', 'lbl']: print(" (set_artdeco_var) ERROR") print(" mode must be -mono-, -lbl- or -kdis-") exit() set_char(artdeco.mcommon.mode, artdeco_in.mode) # kdis model if artdeco_in.mode == 'kdis': set_char(artdeco.mcommon.kdis_model, artdeco_in.kdis_model) if artdeco_in.filters == 'none': artdeco.mcommon.filter = False artdeco.mcommon.nfilter = 0 else: print("(set_artdeco_var) ERROR") print(" filter use not implemented") exit() ################################ # set wlambda and kdis stuff if artdeco_in.mode == 'mono': artdeco.mcommon.nlambda = len(artdeco_in.wavel) artdeco.mcommon.wlambda = np.array(artdeco_in.wavel) elif artdeco_in.mode == 'lbl': artdeco.mcommon.nlambda = len(atmos.wvl_tauabsgas) artdeco.mcommon.wlambda = np.array(atmos.wvl_tauabsgas) elif artdeco_in.mode == 'kdis': if len(artdeco_in.wavel) == 1: artdeco.mcommon.lambda_bound[0] = artdeco_in.wavel[0] artdeco.mcommon.lambda_bound[1] = artdeco_in.wavel[0] elif len(artdeco_in.wavel) == 2: artdeco.mcommon.lambda_bound[:] = artdeco_in.wavel[:] # set wavelength variables if artdeco_in.wavel[0] == -1.0: artdeco.mcommon.lambda_bound[0] = np.amin(kdis.wvlband[1,:]) if artdeco_in.wavel[1] == -1.0: artdeco.mcommon.lambda_bound[1] = np.amax(kdis.wvlband[2,:]) kdis_flag = np.zeros(kdis.nwvl, dtype='int') nband = 0 if not artdeco.mcommon.filter: if len(artdeco_in.wavel) == 2: # find kdis band that are needed for ik in range(kdis.nwvl): if (kdis.wvlband[2,ik] > artdeco.mcommon.lambda_bound[0]) and (kdis.wvlband[1,ik] < artdeco.mcommon.lambda_bound[1]): kdis_flag[ik] = 1 nband = nband + 1 elif len(artdeco_in.wavel) == 1: if not artdeco.mcommon.lambda_bound[0] in kdis.wvlband[0,:]: print("(set_artdeco_var) ERROR") print(" The wvl you request is not present in your kdis (no kdis band central wvl corresponding)") print(" wvl requested = ", artdeco.mcommon.lambda_bound[0]) print(" kdis central wvl = ", kdis.wvlband[0,:]) exit() for ik in range(kdis.nwvl): if (kdis.wvlband[0,ik] == artdeco_in.wavel[0]): kdis_flag[ik] = 1 nband = nband + 1 artdeco.mcommon.kdis_nsp = kdis.nsp artdeco.mcommon.kdis_nsp_c = kdis.nsp_c artdeco.mcommon.kdis_c_desc = -1 if kdis.c_desc == "density": artdeco.mcommon.kdis_c_desc = 1 elif kdis.c_desc == "molar_fraction": artdeco.mcommon.kdis_c_desc = 2 if (kdis.nsp + kdis.nsp_c) > 0: artdeco.mcommon.kdis_nwvl = nband artdeco.mcommon.kdis_nt = kdis.nt artdeco.mcommon.kdis_np = kdis.np artdeco.mcommon.kdis_nmaxai = kdis.nmaxai artdeco.mcommon.kdis_t = kdis.t artdeco.mcommon.kdis_p = kdis.p artdeco.mcommon.kdis_wvlband = kdis.wvlband[:,kdis_flag==1] if solrad.solar_spec!='none': artdeco.mcommon.kdis_solrad_nwvl = nband artdeco.mcommon.kdis_solrad_solcste = solrad.kdis_solrad_solcste artdeco.mcommon.kdis_solrad_wvl = solrad.kdis_solrad_wvl[kdis_flag==1] artdeco.mcommon.kdis_solrad_F0 = solrad.kdis_solrad_F0[kdis_flag==1] else: artdeco.mcommon.kdis_solrad_nwvl = nband artdeco.mcommon.kdis_solrad_solcste = -32768 artdeco.mcommon.kdis_solrad_wvl = kdis.wvlband[0,kdis_flag==1] artdeco.mcommon.kdis_solrad_F0 = np.ones(nband) #print artdeco.mcommon.kdis_solrad_F0 if (kdis.nsp>0): artdeco.mcommon.kdis_nai = kdis.nai[:,kdis_flag==1] artdeco.mcommon.kdis_ki = kdis.ki[:,kdis_flag==1,:,:,:] artdeco.mcommon.kdis_xsect = kdis.xsect[:,kdis_flag==1] artdeco.mcommon.kdis_ai = kdis.ai[:,kdis_flag==1,:] artdeco.mcommon.kdis_fcont = kdis.fcont[:] artdeco.mf2py_utility.allocate_kdis_species(artdeco.mcommon.kdis_nsp) for isp in range(artdeco.mcommon.kdis_nsp): artdeco.mf2py_utility.set_kdis_species(isp+1, kdis.species[isp].strip()) if (kdis.nsp_c>0): artdeco.mcommon.kdis_nc = kdis.nc artdeco.mcommon.kdis_c = kdis.c artdeco.mcommon.kdis_nai_c = kdis.nai_c[:,kdis_flag==1] artdeco.mcommon.kdis_ki_c = kdis.ki_c[:,kdis_flag==1,:,:,:,:] artdeco.mcommon.kdis_xsect_c = kdis.xsect_c[:,kdis_flag==1] artdeco.mcommon.kdis_ai_c = kdis.ai_c[:,kdis_flag==1,:] artdeco.mcommon.kdis_fcont_c = kdis.fcont_c[:] artdeco.mf2py_utility.allocate_kdis_species_c(artdeco.mcommon.kdis_nsp_c) for isp in range(artdeco.mcommon.kdis_nsp_c): artdeco.mf2py_utility.set_kdis_species_c(isp+1, kdis.species_c[isp].strip()) artdeco.mcommon.wlambda = artdeco.mcommon.kdis_wvlband[0,:] artdeco.mcommon.nlambda = len(artdeco.mcommon.wlambda) # wlambda just match kdis_wvlband[0,:] so artdeco.mcommon.lambda_ikdis = np.arange(artdeco.mcommon.nlambda, dtype='int') # +1 cause index starts at 0 in fortran artdeco.mcommon.lambda_ikdis[:] = artdeco.mcommon.lambda_ikdis[:]+1 #print artdeco.mcommon.lambda_ikdis #print artdeco.mcommon.wlambda #print artdeco.mcommon.kdis_wvlband[0,:] else: print("(set_artdeco_var) ERROR") print(" filter use not implemented") exit() else: print("Error") exit() ############## # Geometry if ( len(geom.sza) > 1 ) and ( artdeco_in.rt_model == 'disort' ): print("(set_artdeco_var) ERROR") print(" Only one SZA can be provided at a time") exit() artdeco.mcommon.get_possol = False artdeco.mcommon.sza = geom.sza artdeco.mcommon.nsza = len(artdeco.mcommon.sza) artdeco.mcommon.vza = geom.vza artdeco.mcommon.nvza = len(artdeco.mcommon.vza) artdeco.mcommon.vaa = geom.vaa artdeco.mcommon.nvaa = len(artdeco.mcommon.vaa) #print artdeco.mcommon.vza ################## # set solar flux if artdeco.mcommon.thermal_only: artdeco.mcommon.f0 = np.zeros(artdeco.mcommon.nlambda) else: if artdeco_in.mode == 'mono': artdeco.mcommon.f0 = np.ones(artdeco.mcommon.nlambda) elif artdeco_in.mode == 'kdis': artdeco.mcommon.f0 = artdeco.mcommon.kdis_solrad_F0 elif artdeco_in.mode == 'lbl': if solrad.solar_spec != 'none': print("(set_artdeco_var) ERROR") print(" specific solrad ,ot implemented for lbl mode") exit() else: artdeco.mcommon.f0 = np.ones(artdeco.mcommon.nlambda) artdeco.mcommon.varsol_fact = np.ones(artdeco.mcommon.nsza) artdeco.mcommon.svin = np.array([1.0,0.0,0.0,0.0]) #print artdeco.mcommon.f0 ################# # Polarisation artdeco.mcommon.nmat = artdeco_in.nmat ############# # Thermal artdeco.mcommon.thermal = artdeco_in.thermal ############# # RT stuff if save_layers: set_char(artdeco.mcommon.rt_model, "sinsca") else: set_char(artdeco.mcommon.rt_model, artdeco_in.rt_model) if get_char(artdeco.mcommon.rt_model) not in ['disort','doad','sinsca']: print("(set_artdeco_var) ERROR") print(" rt_model must be -disort-, -sinsca- or -doad-") exit() # RT model if get_char(artdeco.mcommon.rt_model) == 'doad': nstreams = int(artdeco_in.nstreams/2) artdeco.mcommon.doad_eps = artdeco_in.doad_eps artdeco.mcommon.doad_nmug = nstreams artdeco.mcommon.doad_nfoumx = artdeco_in.doad_nfoumx artdeco.mcommon.doad_ncoef_min_brdf = artdeco_in.doad_ncoef_min_brdf #print artdeco.mcommon.doad_eps #print artdeco.mcommon.doad_nmug #print artdeco.mcommon.doad_nfoumx #print "artdeco.mcommon.doad_ncoef_min_brdf =", artdeco.mcommon.doad_ncoef_min_brdf if artdeco.mcommon.thermal: print(" ") print(" (set_artdeco_var) ERROR") print(" Thermal emission can not be accounted for with DOAD") print(" ") exit() elif get_char(artdeco.mcommon.rt_model) == 'disort': artdeco.mcommon.od_nstr = check_disort_nstr(artdeco_in.od_nstr, artdeco.mcommon.sza[0]) nstreams = artdeco.mcommon.od_nstr artdeco.mcommon.od_deltam = artdeco_in.od_deltam artdeco.mcommon.od_corint = artdeco_in.od_corint artdeco.mcommon.od_do_secsca = artdeco_in.od_do_secsca artdeco.mcommon.od_accur = artdeco_in.od_accur artdeco.mcommon.od_use_fisot = artdeco_in.od_use_fisot if artdeco.mcommon.nmat != 1 : print(" ") print(" (set_artdeco_var) ERROR") print(" (set_artdeco_var) nmat must be = 1 with DISORT") print(" ") exit() elif get_char(artdeco.mcommon.rt_model) == 'sinsca': nstreams = artdeco_in.nstreams if artdeco.mcommon.thermal: print(" ") print(" (set_artdeco_var) ERROR") print(" Thermal emission can not be accounted for with -sinsca-") print(" ") exit() # truncation method if save_layers: set_char(artdeco.mcommon.trunc_method, "none") else: set_char(artdeco.mcommon.trunc_method, artdeco_in.trunc_method) if get_char(artdeco.mcommon.trunc_method) not in ['none','dm','dfit','potter']: print("(set_artdeco_var) ERROR") print(" trunc_method must be -none-, -potter-, -dfit- or -dm-") exit() if (get_char(artdeco.mcommon.trunc_method) == 'potter') and (get_char(artdeco.mcommon.rt_model) == 'disort'): print("(set_artdeco_var) ERROR") print(" Potter truncation can not be used with DISORT") exit() if get_char(artdeco.mcommon.trunc_method) == 'dfit': artdeco.mcommon.dfit_thetac = artdeco_in.dfit_thetac artdeco.mcommon.dfit_fitall = artdeco_in.dfit_fitall if get_char(artdeco.mcommon.trunc_method) == 'potter': artdeco.mcommon.potter_theta_min = artdeco_in.potter_theta_min artdeco.mcommon.potter_theta_max = artdeco_in.potter_theta_max # set nbetal if (get_char(artdeco.mcommon.rt_model) == 'sinsca') and (get_char(artdeco.mcommon.trunc_method) == 'none'): artdeco.mcommon.nbetal_in = -1 elif get_char(artdeco.mcommon.rt_model) == 'doad': artdeco.mcommon.nbetal_in = nstreams * 2 else: artdeco.mcommon.nbetal_in = nstreams if artdeco.mcommon.nbetal_in > 0: artdeco.mcommon.nmaxbetal = artdeco.mcommon.nbetal_in else: artdeco.mcommon.nmaxbetal = 0 ################## # particles #print ptcle.__dict__.keys() artdeco.mcommon.nptcle = ptcle.nptcle if artdeco.mcommon.nptcle>0: artdeco.mf2py_utility.allocate_ptcle_char_arr(artdeco.mcommon.nptcle) for iptcle in range(artdeco.mcommon.nptcle): artdeco.mf2py_utility.set_ptcle_vdist(iptcle+1, ptcle.vdist_type[iptcle].strip()) artdeco.mf2py_utility.set_ptcle_type(iptcle+1, ptcle.type[iptcle].strip()) artdeco.mcommon.ptcle_h_vdist = ptcle.h_vdist artdeco.mcommon.ptcle_sd_gauss_vdist = ptcle.sd_gauss_vdist artdeco.mcommon.ptcle_hg = ptcle.hg # this is not used in pyartdeco artdeco.mcommon.ptcle_opt_interp = np.full(artdeco.mcommon.nptcle, False, dtype='bool') artdeco.mcommon.ptcle_tau_reflamb = ptcle.tau_reflamb artdeco.mcommon.reflamb = ptcle.wlref artdeco.mcommon.ptcle_cext_reflamb = np.zeros(artdeco.mcommon.nptcle) for iptcle in range(artdeco.mcommon.nptcle): artdeco.mcommon.ptcle_cext_reflamb[iptcle] = ptcle.optical_properties[iptcle]["cext_reflamb"] if 'user' in ptcle.vdist_type : artdeco.mcommon.ptcle_user_vdist_nalt = ptcle.user_vdist_nalt artdeco.mcommon.ptcle_user_vdist_alt = ptcle.user_vdist_alt artdeco.mcommon.ptcle_user_vdist = ptcle.user_vdist artdeco.mcommon.ptcle_nelem = np.zeros(artdeco.mcommon.nptcle, dtype="int") for iptcle in range(artdeco.mcommon.nptcle): artdeco.mcommon.ptcle_nelem[iptcle] = ptcle.optical_properties[iptcle]["nelem"] artdeco.mcommon.nang_max = ptcle.nang_phasemat artdeco.mcommon.ptcle_nang_phasemat = np.zeros((artdeco.mcommon.nptcle, artdeco.mcommon.nlambda), dtype='int') artdeco.mcommon.ptcle_nang_phasemat[:,:] = ptcle.nang_phasemat artdeco.mcommon.ptcle_opt = np.zeros((artdeco.mcommon.nptcle, 3, artdeco.mcommon.nlambda)) artdeco.mcommon.ptcle_u_phasemat = np.zeros((artdeco.mcommon.nptcle, artdeco.mcommon.nang_max, artdeco.mcommon.nlambda)) artdeco.mcommon.ptcle_phasemat = np.zeros((artdeco.mcommon.nptcle, 6, artdeco.mcommon.nang_max, artdeco.mcommon.nlambda)) artdeco.mcommon.ptcle_betal_flag = np.full((artdeco.mcommon.nptcle, artdeco.mcommon.nlambda), False, dtype=bool) artdeco.mcommon.ptcle_nbetal = np.zeros((artdeco.mcommon.nptcle, artdeco.mcommon.nlambda), dtype='int') artdeco.mf2py_utility.alloc_ptcle_betal() artdeco.mcommon.ptcle_trunccoeff = np.zeros((artdeco.mcommon.nptcle, artdeco.mcommon.nlambda)) artdeco.mcommon.ptcle_trunc_phasemat = np.zeros((artdeco.mcommon.nptcle, 6, artdeco.mcommon.nang_max, artdeco.mcommon.nlambda)) artdeco.mcommon.ptcle_trunc_normphasemat = np.zeros((artdeco.mcommon.nptcle, artdeco.mcommon.nlambda)) opt, u_phasemat, phasemat, betal, betal_flag, trunc_coeff, trunc_phasemat, trunc_normphasemat = \ ptcle.get_opt(artdeco.mcommon.wlambda, artdeco.mcommon.nbetal_in, get_char(artdeco.mcommon.trunc_method)) artdeco.mcommon.ptcle_betal_flag[:,:] = betal_flag[:,:] if artdeco.mcommon.nbetal_in>0 : artdeco.mcommon.ptcle_nbetal[:,:] = artdeco.mcommon.nbetal_in artdeco.mcommon.ptcle_betal[:,:,:,:] = betal[:,:,:,:] artdeco.mcommon.ptcle_trunccoeff[:,:] = trunc_coeff[:,:] artdeco.mcommon.ptcle_trunc_normphasemat[:,:] = trunc_normphasemat[:,:] artdeco.mcommon.ptcle_trunc_phasemat[:,:,:,:] = trunc_phasemat[:,:,:,:] artdeco.mcommon.ptcle_opt[:,:,:] = opt[:,:,:] artdeco.mcommon.ptcle_phasemat[:,:,:,:] = phasemat[:,:,:,:] for iwl in range(artdeco.mcommon.nlambda): for iptcle in range(artdeco.mcommon.nptcle): artdeco.mcommon.ptcle_u_phasemat[iptcle,:,iwl] = ptcle.u_phasemat[:] # for iwl in range(artdeco.mcommon.nlambda): # for iptcle in range(artdeco.mcommon.nptcle): # for iang in range(artdeco.mcommon.ptcle_nang_phasemat[iptcle,iwl]): # print(artdeco.mcommon.ptcle_u_phasemat[iptcle,iang,iwl], ",") # for iwl in range(artdeco.mcommon.nlambda): # for iptcle in range(artdeco.mcommon.nptcle): # for iang in range(artdeco.mcommon.ptcle_nang_phasemat[iptcle,iwl]): # print(artdeco.mcommon.ptcle_phasemat[iptcle,0,iang,iwl], ",") # exit() ####################################################################################### else: if verbose: print("(set_artdeco_var) no particle") ######### # TMS if save_layers : artdeco.mcommon.sinsca_corint = False artdeco.mcommon.sinsca_corint_brdf = False else: if get_char(artdeco.mcommon.rt_model) != 'sinsca' : artdeco.mcommon.sinsca_corint = artdeco_in.corint artdeco.mcommon.sinsca_corint_brdf = artdeco_in.corint_brdf else: artdeco.mcommon.sinsca_corint = False artdeco.mcommon.sinsca_corint_brdf = False ################## # atmosphere #print atmos.__dict__.keys() alt_atm, t_atm, p_atm, u_air_atm, u_atm = atmos.regrid() artdeco.mcommon.nalt_atm = len(alt_atm) artdeco.mcommon.nlayers = artdeco.mcommon.nalt_atm - 1 artdeco.mcommon.alt_atm = alt_atm artdeco.mcommon.p_atm = p_atm artdeco.mcommon.t_atm = t_atm artdeco.mcommon.u_air_atm = u_air_atm # for ialt in range(len(alt_atm)): # print(" %.3f"%artdeco.mcommon.alt_atm[ialt]+" %.3f"%artdeco.mcommon.p_atm[ialt]+" %.3f"%artdeco.mcommon.t_atm[ialt]) if artdeco_in.mode == 'lbl': artdeco.mcommon.lblod_atm = atmos.tauabsgas if atmos.ngas > 0: artdeco.mcommon.u_atm = u_atm artdeco.mcommon.ngas_u_atm = len(atmos.gas) artdeco.mf2py_utility.allocate_gas_u_atm(artdeco.mcommon.ngas_u_atm) for isp in range(artdeco.mcommon.ngas_u_atm): artdeco.mf2py_utility.set_gas_u_atm(isp+1, atmos.gas[isp].strip()) artdeco.mcommon.depol = np.zeros((artdeco.mcommon.nlambda)) depol = atmos.get_depol(artdeco.mcommon.wlambda) #print depol artdeco.mcommon.depol[:] = depol[:] artdeco.mcommon.depol[depol==-1] = artdeco.mconstants.undef_dp if len(atmos.tau_ray.keys())>0: artdeco.mcommon.ray_od = atmos.get_ray_od(artdeco.mcommon.p_atm, artdeco.mcommon.wlambda) ############# # Surface artdeco.mcommon.test_glitter = surface.test_glitter artdeco.mcommon.ocean_whitecaps = surface.ocean_whitecaps if surface.family not in ['brdf', 'lambert']: print("(set_artdeco_var) ERROR") print(" Surface family must be -brdf- or -lambert-") exit() #print surface.__dict__.keys() set_char(artdeco.mcommon.surface_family, surface.family) set_char(artdeco.mcommon.surface_type, surface.kind) if surface.family == 'lambert': artdeco.mcommon.surface_albedo = np.zeros(artdeco.mcommon.nlambda) albedo = surface.get_it(artdeco.mcommon.wlambda, False) artdeco.mcommon.surface_albedo[:] = albedo if artdeco.mcommon.bottom_f0 and np.count_nonzero(artdeco.mcommon.surface_albedo)!=0: print("(set_artdeco_var) ERROR") print(" No surface allowed if -bottom_f0-") exit() elif surface.family == 'brdf': if artdeco.mcommon.bottom_f0 : print("(set_artdeco_var) ERROR") print(" No surface allowed if -bottom_f0-") exit() if surface.kind == 'ocean': surface_brdf_model = 'ocean' set_char(artdeco.mcommon.surface_brdf_model, surface_brdf_model) artdeco.mcommon.wind_spd = surface.wdspd artdeco.mcommon.ocean_xsal = surface.xsal artdeco.mcommon.ocean_pcl = surface.pcl artdeco.mcommon.ocean_shadow = surface.shadow artdeco.mcommon.ocean_distr = surface.sdistr artdeco.mcommon.ocean_azw = surface.azw artdeco.mcommon.ocean_rsw = np.zeros(artdeco.mcommon.nlambda) artdeco.mcommon.ocean_rsw[:] = surface.get_Rsw(artdeco.mcommon.wlambda) #print("artdeco.mcommon.ocean_rsw[:] = ",artdeco.mcommon.ocean_rsw[:]) artdeco.mcommon.ocean_6s_glitter = surface._6s_glitter elif surface.kind == 'land': surface_brdf_model = 'lisp_ross' set_char(artdeco.mcommon.surface_brdf_model, surface_brdf_model) artdeco.mcommon.surface_n_par_brdf = 3 artdeco.mcommon.surface_par_brdf = np.zeros((artdeco.mcommon.surface_n_par_brdf, artdeco.mcommon.nlambda)) artdeco.mcommon.surface_n_par_bpdf = 3 artdeco.mcommon.surface_par_bpdf = np.zeros((artdeco.mcommon.surface_n_par_bpdf, artdeco.mcommon.nlambda)) if artdeco.mcommon.nmat==1: iso, vol, geo = surface.get_it(artdeco.mcommon.wlambda, False) #print iso #print vol #print geo artdeco.mcommon.surface_par_brdf[0, :] = iso artdeco.mcommon.surface_par_brdf[1, :] = vol / iso artdeco.mcommon.surface_par_brdf[2, :] = geo / iso else: iso, vol, geo, cv, nr, ni = surface.get_it(artdeco.mcommon.wlambda, True) # print iso # print vol # print geo # print cv # print nr # print ni artdeco.mcommon.surface_par_brdf[0, :] = iso artdeco.mcommon.surface_par_brdf[1, :] = vol / iso artdeco.mcommon.surface_par_brdf[2, :] = geo / iso artdeco.mcommon.surface_par_bpdf[0, :] = cv artdeco.mcommon.surface_par_bpdf[1, :] = nr artdeco.mcommon.surface_par_bpdf[2, :] = ni else: print("(set_artdeco_var) ERROR") print(" Surface kind must be -land- or -ocean-") print(" if family is -brdf-") exit() artdeco.mcommon.surface_temp = surface.temp if artdeco.mcommon.bottom_f0 and artdeco.mcommon.surface_temp>0: print("(set_artdeco_var) ERROR") print(" No surface allowed if -bottom_f0-") print(" (set surf temp = 0)") exit() artdeco.mcommon.rt_cputime = np.zeros( artdeco.mcommon.nlambda ) artdeco.mcommon.rho_surf = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nlambda)) artdeco.mcommon.svr = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) artdeco.mcommon.svr_levels = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) artdeco.mcommon.svr_sig = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) artdeco.mcommon.flux_out = np.zeros((artdeco.mcommon.nsza, 3, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) if artdeco.mcommon.sinsca_corint : artdeco.mcommon.svr_notms = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) artdeco.mcommon.svr_levels_notms = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) artdeco.mcommon.svr_sig_notms = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) if artdeco_in.mode =='kdis': artdeco.mcommon.fluxint_out = np.zeros((artdeco.mcommon.nsza, 3, artdeco.mcommon.nalt_atm)) if verbose: print("(set_artdeco_var) exit") ####################################################################################################################### def select_kdis_gas(species_in, artdeco): # This routines modify the ARTDECO CKD arrays for the KDIS # to select only gases that are listed in "species" # # This must take place righ after the call for artdeco.minout.read_data() # # artdeco.kdis varaibale must pre-exist for that # routine to be called # print 'old' # print artdeco.mcommon.kdis_nsp # for isp_kdis in range( artdeco.mcommon.kdis_nsp ): # print artdeco.mf2py_utility.print_kdis_species(isp_kdis+1).strip(), artdeco.mcommon.kdis_nai[isp_kdis, 0] # print "" # print "" if not 'all' in species_in: species = [] species_c = [] for isp in range( len(species_in) ): for isp_kdis in range(artdeco.mcommon.kdis_nsp): kdis_sp = artdeco.mf2py_utility.print_kdis_species(isp_kdis+1) if isinstance(kdis_sp, str): kdis_sp = kdis_sp.strip() elif isinstance(kdis_sp, bytes): kdis_sp = str(kdis_sp, 'utf-8').strip() if kdis_sp == species_in[isp].strip() : species.append(kdis_sp) break for isp_kdis in range(artdeco.mcommon.kdis_nsp_c): kdis_sp = artdeco.mf2py_utility.print_kdis_species_c(isp_kdis+1) if isinstance(kdis_sp, str): kdis_sp = kdis_sp.strip() elif isinstance(kdis_sp, bytes): kdis_sp = str(kdis_sp, 'utf-8').strip() if kdis_sp == species_in[isp].strip() : species_c.append(kdis_sp) break if (artdeco.mcommon.kdis_nsp > 0) : nspecies = len(species) old_index = np.zeros(nspecies, dtype=int) kdis_nai_old = np.zeros((artdeco.mcommon.kdis_nsp,artdeco.mcommon.kdis_nwvl), dtype='int') kdis_ai_old = np.zeros((artdeco.mcommon.kdis_nsp,artdeco.mcommon.kdis_nwvl,artdeco.mcommon.kdis_nmaxai)) kdis_ki_old = np.zeros((artdeco.mcommon.kdis_nsp,artdeco.mcommon.kdis_nwvl,artdeco.mcommon.kdis_nmaxai,artdeco.mcommon.kdis_np,artdeco.mcommon.kdis_nt)) kdis_xsect_old = np.zeros((artdeco.mcommon.kdis_nsp,artdeco.mcommon.kdis_nwvl)) kdis_fcont_old = np.zeros( artdeco.mcommon.kdis_nsp) kdis_nsp_old = artdeco.mcommon.kdis_nsp kdis_nai_old[:,:] = artdeco.mcommon.kdis_nai[:,:] # kdis_nai(kdis_nsp,kdis_nwvl) kdis_ai_old[:,:,:] = artdeco.mcommon.kdis_ai[:,:,:] # kdis_ai(kdis_nsp,kdis_nwvl,kdis_nmaxai) kdis_ki_old[:,:,:,:,:] = artdeco.mcommon.kdis_ki[:,:,:,:,:] # kdis_ki(kdis_nsp,kdis_nwvl,kdis_nmaxai,kdis_np,kdis_nt) kdis_xsect_old[:,:] = artdeco.mcommon.kdis_xsect[:,:] # kdis_xsect(kdis_nsp,kdis_nwvl) kdis_fcont_old[:] = artdeco.mcommon.kdis_fcont[:] # a weight to the continuum to be added # check if the gas species are defined for isp in range(nspecies): flag = False for isp_kdis in range(kdis_nsp_old): kdis_sp = artdeco.mf2py_utility.print_kdis_species(isp_kdis+1) if isinstance(kdis_sp, str): kdis_sp = kdis_sp.strip() elif isinstance(kdis_sp, bytes): kdis_sp = str(kdis_sp, 'utf-8').strip() if kdis_sp == species[isp].strip() : flag = True old_index[isp] = isp_kdis #print species[isp].strip(), kdis_sp.strip(), old_index[isp] if flag == False: print('') print('(py select_kdis_gas) Error') print('(py select_kdis_gas) Specie ', species[isp].strip(), ' is not defined in KDIS') print('') exit() # Reset ARTDECO arrays artdeco.mcommon.kdis_nsp = nspecies for isp_kdis in range(kdis_nsp_old): artdeco.mf2py_utility.set_kdis_species(isp_kdis+1, 'unkown') artdeco.mcommon.kdis_nai[:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_ai[:,:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_ki[:,:,:,:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_xsect[:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_fcont[:] = artdeco.mconstants.undef_dp for isp in range(nspecies): artdeco.mf2py_utility.set_kdis_species(isp+1, species[isp].strip()) artdeco.mcommon.kdis_nai[isp,:] = kdis_nai_old[old_index[isp], :] artdeco.mcommon.kdis_ai[isp,:,:] = kdis_ai_old[old_index[isp],:,:] artdeco.mcommon.kdis_ki[isp,:,:,:,:] = kdis_ki_old[old_index[isp],:,:,:,:] artdeco.mcommon.kdis_xsect[isp,:] = kdis_xsect_old[old_index[isp],:] artdeco.mcommon.kdis_fcont[isp] = kdis_fcont_old[old_index[isp]] if (artdeco.mcommon.kdis_nsp_c > 0) : nspecies_c = len(species_c) old_index_c = np.zeros(nspecies_c, dtype=int) kdis_nai_c_old = np.zeros((artdeco.mcommon.kdis_nsp_c,artdeco.mcommon.kdis_nwvl), dtype='int') kdis_ai_c_old = np.zeros((artdeco.mcommon.kdis_nsp_c,artdeco.mcommon.kdis_nwvl,artdeco.mcommon.kdis_nmaxai)) kdis_ki_c_old = np.zeros((artdeco.mcommon.kdis_nsp_c,artdeco.mcommon.kdis_nwvl,artdeco.mcommon.kdis_nmaxai,artdeco.mcommon.kdis_np,artdeco.mcommon.kdis_nt, artdeco.mcommon.kdis_nc)) kdis_xsect_c_old = np.zeros((artdeco.mcommon.kdis_nsp_c,artdeco.mcommon.kdis_nwvl)) kdis_fcont_c_old = np.zeros( artdeco.mcommon.kdis_nsp_c) kdis_nsp_c_old = artdeco.mcommon.kdis_nsp_c kdis_nai_c_old[:,:] = artdeco.mcommon.kdis_nai_c[:,:] # kdis_nai(kdis_nsp,kdis_nwvl) kdis_ai_c_old[:,:,:] = artdeco.mcommon.kdis_ai_c[:,:,:] # kdis_ai(kdis_nsp,kdis_nwvl,kdis_nmaxai) kdis_ki_c_old[:,:,:,:,:,:] = artdeco.mcommon.kdis_ki_c[:,:,:,:,:,:] # kdis_ki(kdis_nsp,kdis_nwvl,kdis_nmaxai,kdis_np,kdis_nt,kdis_nc) kdis_xsect_c_old[:,:] = artdeco.mcommon.kdis_xsect_c[:,:] # kdis_xsect(kdis_nsp,kdis_nwvl) kdis_fcont_c_old[:] = artdeco.mcommon.kdis_fcont_c[:] # a weight to the continuum to be added # check if the gas species are defined for isp in range(nspecies_c): flag = False for isp_kdis in range(kdis_nsp_c_old): kdis_sp = artdeco.mf2py_utility.print_kdis_species_c(isp_kdis+1) if isinstance(kdis_sp, str): kdis_sp = kdis_sp.strip() elif isinstance(kdis_sp, bytes): kdis_sp = str(kdis_sp, 'utf-8').strip() if kdis_sp == species_c[isp].strip() : flag = True old_index_c[isp] = isp_kdis if flag == False: print('') print('(py select_kdis_gas) Error') print('(py select_kdis_gas) Specie ', species_c[isp].strip(), ' is not defined in KDIS') print('') exit() # Reset ARTDECO arrays artdeco.mcommon.kdis_nsp_c = nspecies_c for isp_kdis in range(kdis_nsp_c_old): artdeco.mf2py_utility.set_kdis_species_c(isp_kdis+1, 'unkown') artdeco.mcommon.kdis_nai_c[:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_ai_c[:,:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_ki_c[:,:,:,:,:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_xsect_c[:,:] = artdeco.mconstants.undef_dp artdeco.mcommon.kdis_fcont_c[:] = artdeco.mconstants.undef_dp for isp in range(nspecies_c): artdeco.mf2py_utility.set_kdis_species_c(isp+1, species_c[isp].strip()) artdeco.mcommon.kdis_nai_c[isp,:] = kdis_nai_c_old[old_index_c[isp], :] artdeco.mcommon.kdis_ai_c[isp,:,:] = kdis_ai_c_old[old_index_c[isp],:,:] artdeco.mcommon.kdis_ki_c[isp,:,:,:,:,:] = kdis_ki_c_old[old_index_c[isp],:,:,:,:,:] artdeco.mcommon.kdis_xsect_c[isp,:] = kdis_xsect_c_old[old_index_c[isp],:] artdeco.mcommon.kdis_fcont_c[isp] = kdis_fcont_c_old[old_index_c[isp]] ################################################################################### def get_char(chartdeco): # print(chartdeco, chartdeco.__class__, chartdeco.dtype) # print(chartdeco) # exit() if len(chartdeco.shape) == 0: # after a change in f2py format for character string return str(chartdeco, 'utf-8').strip().strip("\x00") else: charout = '' for cc in chartdeco: if cc != ' ': charout = charout + cc else: return charout.strip() def cmp_char(chartdeco, charin): if len(chartdeco.shape) == 0: # after a change in f2py format for character string if str(chartdeco).strip() == charin: return True else: return False else: if type(charin) is list: tmp = charin else: tmp = list(charin) n = len(tmp) n2 = len(chartdeco) if ( n > n2): return False charf2py = np.chararray(n2) for i in range(n): charf2py[i] = tmp[i] for i in range(n, n2): charf2py[i] = ' ' if (charf2py == chartdeco).all(): return True else: return False def set_char(charf2py, charin): # Routine to set the charcater # string variable (numpy array of string) # of f2py module # print("(set_char)", charf2py) # print("(set_char)", hex(id(charf2py))) # print("(set_char)", hex(id(charin))) # exit() if len(charf2py.shape) == 0: # after a change in f2py format for character string if len(charin)> charf2py.dtype.itemsize: print(" (py) Pb in set_char()") exit() charf2py[...] = charin+( charf2py.dtype.itemsize - len(charin))*" " else: if type(charin) is list: tmp = charin else: tmp = list(charin) n = len(tmp) n2 = len(charf2py) if ( n > len(charf2py)): print(" (py) Pb in set_char()") exit() for i in range(n): charf2py[i] = tmp[i] for i in range(n, n2): charf2py[i] = ' ' # print("(set_char)", charf2py) # print("(set_char)", hex(id(charf2py))) return ################################################################################### def deallocate(artdeco): # in ordrer of occurance in artdeco_common.f90 artdeco.mcommon.rt_cputime = None artdeco.mcommon.rt_cputime_filter = None artdeco.mcommon.sza = None artdeco.mcommon.vza = None artdeco.mcommon.vaa = None artdeco.mcommon.day = None artdeco.mcommon.h_tu = None artdeco.mcommon.varsol_fact = None artdeco.mcommon.filter_name = None artdeco.mcommon.nlamb_filter = None artdeco.mcommon.bound_lamb_filter = None artdeco.mcommon.lamb_filter = None artdeco.mcommon.trans_filter = None artdeco.mcommon.solrad_wvl = None artdeco.mcommon.solrad_F0 = None artdeco.mcommon.f0 = None artdeco.mcommon.alt_atm = None artdeco.mcommon.p_atm = None artdeco.mcommon.t_atm = None artdeco.mcommon.u_air_atm = None artdeco.mcommon.gas_u_atm = None artdeco.mcommon.u_atm = None artdeco.mcommon.ray_od = None artdeco.mcommon.lblod_atm = None artdeco.mcommon.kdis_nt = 0 artdeco.mcommon.kdis_np = 0 artdeco.mcommon.kdis_nc = 0 artdeco.mcommon.kdis_nwvl = 0 artdeco.mcommon.kdis_t = None artdeco.mcommon.kdis_p = None artdeco.mcommon.kdis_c = None artdeco.mcommon.kdis_wvlband = None artdeco.mcommon.kdis_nsp = 0 artdeco.mcommon.kdis_species = None artdeco.mcommon.kdis_nai = None artdeco.mcommon.kdis_ai = None artdeco.mcommon.kdis_ki = None artdeco.mcommon.kdis_xsect = None artdeco.mcommon.kdis_fcont = None artdeco.mcommon.kdis_nsp_c = 0 artdeco.mcommon.kdis_species_c = None artdeco.mcommon.kdis_nai_c = None artdeco.mcommon.kdis_ai_c = None artdeco.mcommon.kdis_ki_c = None artdeco.mcommon.kdis_xsect_c = None artdeco.mcommon.kdis_fcont_c = None artdeco.mcommon.kdis_solrad_f0 = None artdeco.mcommon.kdis_solrad_wvl = None artdeco.mcommon.wlambda = None artdeco.mcommon.lambda_ikdis = None artdeco.mcommon.depol = None artdeco.mcommon.ptcle_type = None artdeco.mcommon.ptcle_hg = None artdeco.mcommon.ptcle_opt_interp = None artdeco.mcommon.ptcle_tau_reflamb = None artdeco.mcommon.ptcle_vdist_type = None artdeco.mcommon.ptcle_h_vdist = None artdeco.mcommon.ptcle_sd_gauss_vdist = None artdeco.mcommon.ptcle_user_vdist_nalt = None artdeco.mcommon.ptcle_user_vdist_alt = None artdeco.mcommon.ptcle_user_vdist = None artdeco.mcommon.ptcle_opt_model = None artdeco.mcommon.ptcle_material = None artdeco.mcommon.ptcle_param_mie3 = None artdeco.mcommon.ptcle_nphot_hm = None artdeco.mcommon.ptcle_param_hm = None artdeco.mcommon.ptcle_param_baum = None artdeco.mcommon.ptcle_opt_exist_reflamb = None artdeco.mcommon.ptcle_opt_exist = None artdeco.mcommon.ptcle_nelem = None artdeco.mcommon.ptcle_flag_store_opt = None artdeco.mcommon.ptcle_opt_store = None artdeco.mcommon.ptcle_nlamb_phasemat_store = None artdeco.mcommon.ptcle_nang_phasemat_store = None artdeco.mcommon.ptcle_lamb_phasemat_store = None artdeco.mcommon.ptcle_u_phasemat_store = None artdeco.mcommon.ptcle_phasemat_store = None artdeco.mcommon.ptcle_flag_new_opt = None artdeco.mcommon.ptcle_opt_new = None artdeco.mcommon.ptcle_nlamb_phasemat_new = None artdeco.mcommon.ptcle_nang_phasemat_new = None artdeco.mcommon.ptcle_lamb_phasemat_new = None artdeco.mcommon.ptcle_u_phasemat_new = None artdeco.mcommon.ptcle_phasemat_new = None artdeco.mcommon.ptcle_ilambda = None artdeco.mcommon.ptcle_ireflambda = None artdeco.mcommon.ptcle_cext_reflamb = None artdeco.mcommon.ptcle_opt = None artdeco.mcommon.ptcle_nang_phasemat = None artdeco.mcommon.ptcle_u_phasemat = None artdeco.mcommon.ptcle_phasemat = None artdeco.mcommon.ptcle_nbetal = None artdeco.mcommon.ptcle_betal = None artdeco.mcommon.ptcle_betal_flag = None artdeco.mcommon.ptcle_trunc_normphasemat = None artdeco.mcommon.ptcle_trunccoeff = None artdeco.mcommon.ptcle_trunc_phasemat = None artdeco.mcommon.material_name = None artdeco.mcommon.material_nlamb_refind = None artdeco.mcommon.material_refind = None artdeco.mcommon.nai_layers = None artdeco.mcommon.ai_layers = None artdeco.mcommon.gas_dtauabs_layers= None artdeco.mcommon.ptcle_dtauabs_layers= None artdeco.mcommon.tot_dtausca_layers= None artdeco.mcommon.nbetal_layers= None artdeco.mcommon.betal_layers= None artdeco.mcommon.ptcle_dtau_layers= None artdeco.mcommon.ptcle_ssa_layers= None artdeco.mcommon.ptcle_nang_phasemat_layers= None artdeco.mcommon.ptcle_u_phasemat_layers= None artdeco.mcommon.ptcle_phasemat_layers= None artdeco.mcommon.gas_dtausca_layers= None artdeco.mcommon.ptcle_dtau_layers_unscaled= None artdeco.mcommon.ptcle_ssa_layers_unscaled= None artdeco.mcommon.ptcle_ssa_layers_tms= None artdeco.mcommon.ptcle_phasemat_tms_layers= None artdeco.mcommon.ptcle_phasemat_actual_layers= None artdeco.mcommon.ptcle_phasemat_delta_layers= None artdeco.mcommon.ptcle_u_phasemat_delta_layers = None artdeco.mcommon.ptcle_secsca_mu_eq = None artdeco.mcommon.ptcle_secsca_neg_phas = None artdeco.mcommon.ptcle_secsca_norm_phas = None artdeco.mcommon.ai_mono_layers= None artdeco.mcommon.tot_dtau_mono_layers= None artdeco.mcommon.ssa_mono_layers= None artdeco.mcommon.gas_dtau_mono_layers= None artdeco.mcommon.gas_ssa_mono_layers= None artdeco.mcommon.surface_albedo = None artdeco.mcommon.surface_par_brdf = None artdeco.mcommon.surface_par_bpdf = None artdeco.mcommon.rho_surf = None artdeco.mcommon.svr = None artdeco.mcommon.svr_levels = None artdeco.mcommon.svr_sig = None artdeco.mcommon.svr_notms = None artdeco.mcommon.svr_levels_notms = None artdeco.mcommon.svr_sig_notms = None artdeco.mcommon.flux_out = None artdeco.mcommon.fluxint_out = None artdeco.mcommon.svr_filter = None artdeco.mcommon.svr_sig_filter = None artdeco.mcommon.flux_filter = None artdeco.mcommon.lamb_baum = None artdeco.mcommon.de_baum = None artdeco.mcommon.mu_baum = None artdeco.mcommon.cext_baum = None artdeco.mcommon.ssa_baum = None artdeco.mcommon.g_baum = None artdeco.mcommon.ptcle_phasemat_baum = None def check_disort_nstr(nstreams_init, sza): tol = 1e-4 valid = False nstreams = nstreams_init mu0 = np.cos(sza * np.pi /180.0) if mu0<0: print("") print("") print("(check_disort_nstr) ERROR") print("(check_disort_nstr) mu0 must be positive") print("") print("") exit() #print 'check_disort_nstr',mu0 while not (valid): GMT, wght = np.polynomial.legendre.leggauss(int(nstreams/2)) # convert from (-1, 1) to (0, 1) as in DISORT GMT = 0.5 * GMT + 0.5 #print 'check_disort_nstr',GMT if True in ((np.abs(mu0-GMT)/mu0)