Convert UCLALES-SALSA output to DEPHY format#
Code to read the LES output files and write to DEPHY format (NetCDF)
Modified from “convert_comble_dephy_forcing_to_DHARMA_LES_and_ModelE3_SCM_forcing.ipynb”
Import libraries#
import xarray as xr
import pandas as pd
import numpy as np
import os
import datetime as dt
from netCDF4 import Dataset
Specify directory locations#
Model name: UCLALES-SALSA
Microphysics:
a) SB (i.e., without SALSA): diagnostic cloud (saturation adjustment + CDNC as an input) and 2M (number and mass) rain-ice/snow/graupel microphysics
b) SALSA: prognostic sectional aerosol-cloud-rain-ice microphysics
# Specify sources and targets
# Name of the LES simulation
#exp_name = 'COMBLE_SB_FixN'
#exp_name = 'COMBLE_SB_FixN_noice'
#exp_name = 'COMBLE_SB_FixN_def_z0'
#exp_name = 'COMBLE_SALSA_ProgNc_noice'
exp_name = 'COMBLE_SALSA_ProgNc'
#exp_name = 'COMBLE_SALSA_ProgNcNi'
# Github scratch directory where processed model output will be committed
my_gitdir = '../../output_les/uclales-salsa/sandbox/'
# Source directory (local disk)
my_rundir = '../../../outputs/'
# Output file name prefix
prefix = 'UCLALES-SALSA_Lx25_dx100' # Small domain
# Output file name
sb='SB' in exp_name # SB or SALSA microphysics
if sb:
# Seifert and Beheng (SB) microphysics
my_outfile = prefix+exp_name[9:] # Remove 'COMBLE_SB'
else:
# SALSA microphysics
my_outfile = prefix+exp_name[12:] # Remove 'COMBLE_SALSA'
# LES inputs
input_filename_ts=my_rundir+exp_name+'.ts.nc'
input_filename_ps=my_rundir+exp_name+'.ps.nc'
input_filename_cs=my_rundir+exp_name+'.cs.nc' # Optional
input_filename_an=my_rundir+exp_name+'.nc' # Optional
if not os.path.exists(input_filename_ts): print('Missing ts data!')
if not os.path.exists(input_filename_ps): print('Missing ps data!')
cs=os.path.exists(input_filename_cs)
an=os.path.exists(input_filename_an)
# Additional note about microphysics
if sb:
# SB
micro = 'Without SALSA, i.e., using the 2M bulk microphysics by Seifert and Beheng (SB). '+\
'This version has diagnostic clouds (saturation adjustment + fixed CDNC) and prognostic 2M rain, '+\
'ice, snow and graupel categories.'
else:
# SALSA
micro = 'UCLALES-SALSA with prognostic sectional aerosol, cloud and ice (called snow) categories.'
# Constants
p00=1.0e+05 # Reference pressure (Pa)
Rd=287.04 # Specific gas constant for dry air (R_specific=R/M), J/kg/K
cp=1005.0 # Specific heat for a constant pressure
alvl=2.5e+6 # Latent heat for water
alvi=2.834e+6 # Latent heat for ice
Prepare output file in DEPHY format#
Read requested variables list#
Variable description, naming, units, and dimensions.
# read list of requested variables
vars_mean_list = pd.read_excel('https://docs.google.com/spreadsheets/d/1Vl8jYGviet7EtXZuQiitrx4NSkV1x27aJAhxxjBb9zI/export?gid=0&format=xlsx',
sheet_name='Mean')
pd.set_option('display.max_rows', None)
vars_mean_list
| standard_name | variable_id | units | dimensions | comment (10-min average reported at endpoints, green=minimum) | |
|---|---|---|---|---|---|
| 0 | time | time | s | – | dimension, seconds since 2020-03-12 18:00:00 |
| 1 | height | zf | m | – | dimension, altitude of mid-level points above ... |
| 2 | layer_top_height | ze | m | – | dimension, altitude of layer top points above ... |
| 3 | surface_pressure | ps | Pa | time | – |
| 4 | surface_temperature | ts | K | time | – |
| 5 | surface_friction_velocity | ustar | m s-1 | time | – |
| 6 | surface_roughness_length_for_momentum_in_air | z0 | m | time | – |
| 7 | surface_roughness_length_for_heat_in_air | z0h | m | time | – |
| 8 | surface_roughness_length_for_humidity_in_air | z0q | m | time | – |
| 9 | surface_upward_sensible_heat_flux | hfss | W m-2 | time | – |
| 10 | surface_upward_latent_heat_flux | hfls | W m-2 | time | – |
| 11 | obukhov_length | ol | m | time | – |
| 12 | atmosphere_mass_content_of_liquid_cloud_water | lwpc | kg m-2 | time | default breakdown to cloud water and rain in a... |
| 13 | atmosphere_mass_content_of_rain_water | lwpr | kg m-2 | time | – |
| 14 | atmosphere_mass_content_of_ice_water | iwp | kg m-2 | time | all ice-phase hydrometeors |
| 15 | cloud_area_fraction | clt | 1 | time | fraction of atmospheric columns with total hyd... |
| 16 | optical_depth | od | 1 | time | scene (all sky); mid-visible, all hydrometeors... |
| 17 | optical_depth_of_liquid_cloud | odlc | 1 | time | as for optical_depth but cloud droplets only |
| 18 | precipitation_flux_at_surface | pr | kg m-2 s-1 | time | liquid and ice phase, all hydrometeors |
| 19 | precipitation_flux_at_surface_in_ice_phase | pri | kg m-2 s-1 | time | ice phase hydrometeors only |
| 20 | toa_outgoing_longwave_flux | rlut | W m-2 | time | at top of atmosphere |
| 21 | surface_downwelling_longwave_flux | rlds | W m-2 | time | – |
| 22 | surface_upwelling_longwave_flux | rlus | W m-2 | time | – |
| 23 | surface_sea_spray_number_flux | ssaf | m-2 s-1 | time | when using prognostic aerosol; emission only (... |
| 24 | air_pressure | pa | Pa | time, height | – |
| 25 | eastward_wind | ua | m s-1 | time, height | – |
| 26 | northward_wind | va | m s-1 | time, height | – |
| 27 | air_dry_density | rhoa | kg m-3 | time, height | – |
| 28 | air_temperature | ta | K | time, height | – |
| 29 | water_vapor_mixing_ratio | qv | kg kg-1 | time, height | per kg dry air |
| 30 | relative_humidity | hur | 1 | time, height | relative to liquid |
| 31 | relative_humidity_over_ice | huri | 1 | time, height | relative to ice |
| 32 | air_potential_temperature | theta | K | time, height | – |
| 33 | specific_turbulent_kinetic_energy_resolved | tke_res | m2 s-2 | time, height | resolved only; compute mean values horizontall... |
| 34 | specific_turbulent_kinetic_energy_sgs | tke_sgs | m2 s-2 | time, height | SGS only |
| 35 | mass_mixing_ratio_of_cloud_liquid_water_in_air | qlc | kg kg-1 | time, height | scene (all sky) per kg dry air; default breakd... |
| 36 | mass_mixing_ratio_of_rain_water_in_air | qlr | kg kg-1 | time, height | – |
| 37 | mass_mixing_ratio_of_cloud_ice_in_air | qic | kg kg-1 | time, height | default breakdown to cloud_ice, snow and graup... |
| 38 | mass_mixing_ratio_of_snow_in_air | qis | kg kg-1 | time, height | – |
| 39 | mass_mixing_ratio_of_graupel_in_air | qig | kg kg-1 | time, height | – |
| 40 | number_of_liquid_cloud_droplets_in_air | nlc | kg-1 | time, height | scene (all sky) per kg dry air; default breakd... |
| 41 | number_of_rain_drops_in_air | nlr | kg-1 | time, height | – |
| 42 | number_of_cloud_ice_crystals_in_air | nic | kg-1 | time, height | default breakdown to cloud_ice, snow and graup... |
| 43 | number_of_snow_crystals_in_air | nis | kg-1 | time, height | – |
| 44 | number_of_graupel_crystals_in_air | nig | kg-1 | time, height | – |
| 45 | number_of_total_aerosol_mode1 | na1 | kg-1 | time, height | when using prognostic aerosol; scene (all sky)... |
| 46 | number_of_total_aerosol_mode2 | na2 | kg-1 | time, height | accumulation mode |
| 47 | number_of_total_aerosol_mode3 | na3 | kg-1 | time, height | sea spray mode |
| 48 | number_of_liquid_cloud_droplets_in_cloud | nlcic | kg-1 | time, height | per kg dry air where cloud water > 0.01 g/kg; ... |
| 49 | number_of_ice_crystals_in_cloud | niic | kg-1 | time, height | total ice hydrometeors per kg dry air where cl... |
| 50 | dissipation_rate_of_turbulent_kinetic_energy | eps | m2 s-3 | time, height | total (including numerical diffusion contribut... |
| 51 | zonal_momentum_flux | uw | m2 s-2 | time, height | total (resolved plus SGS); for resolved compon... |
| 52 | meridional_momentum_flux | vw | m2 s-2 | time, height | total (resolved plus SGS); for resolved compon... |
| 53 | variance_of_upward_air_velocity | w2 | m2 s-2 | time, height | total (resolved plus SGS); for resolved compon... |
| 54 | vertical_flux_potential_temperature | wth | K m s-1 | time, height | total (resolved plus SGS); for resolved compon... |
| 55 | vertical_flux_liquid_ice_water_potential_tempe... | vf_thli | K m s-1 | time, height | total (resolved plus SGS); include sedimentati... |
| 56 | vertical_flux_water_vapor | wqv | kg kg-1 m s-1 | time, height | total (resolved plus SGS); for resolved compon... |
| 57 | vertical_flux_total_water | vf_qt | kg kg-1 m s-1 | time, height | total (resolved plus SGS); vapor+all liquid+al... |
| 58 | area_fraction_of_liquid_cloud | flc | 1 | time, height | fraction of cells with cloud water threshold o... |
| 59 | precipitation_flux_in_air | prf | kg m-2 s-1 | time, height | liquid and ice phase, all hydrometeors |
| 60 | precipitation_flux_in_air_in_ice_phase | prfi | kg m-2 s-1 | time, height | – |
| 61 | downwelling_longwave_flux_in_air | rld | W m-2 | time, height | – |
| 62 | upwelling_longwave_flux_in_air | rlu | W m-2 | time, height | – |
| 63 | tendency_of_air_potential_temperature_due_to_r... | dth_rad | K s-1 | time, height | – |
| 64 | tendency_of_air_potential_temperature_due_to_m... | dth_micro | K s-1 | time, height | including net condensation and precipitation |
| 65 | tendency_of_air_potential_temperature_due_to_m... | dth_turb | K s-1 | time, height | including surface fluxes |
| 66 | tendency_of_water_vapor_mixing_ratio_due_to_mi... | dq_micro | kg kg-1 s-1 | time, height | including net condensation and precipitation |
| 67 | tendency_of_water_vapor_mixing_ratio_due_to_mi... | dq_turb | kg kg-1 s-1 | time, height | including surface fluxes |
| 68 | tendency_of_aerosol_number_due_to_warm_microph... | dna_micro_warm | kg-1 s-1 | time, height | activated and unactivated aerosol (sum over al... |
| 69 | tendency_of_aerosol_number_due_to_cold_microph... | dna_micro_cold | kg-1 s-1 | time, height | activated and unactivated aerosol (sum over al... |
| 70 | tendency_of_aerosol_number_due_to_mixing | dna_turb | kg-1 s-1 | time, height | activated and unactivated aerosol (sum over al... |
| 71 | tendency_of_ice_number_due_to_heterogeneous_fr... | dni_het | kg-1 s-1 | time, height | sum of primary ice nucleation due to activatio... |
| 72 | tendency_of_ice_number_due_to_secondary_ice_pr... | dni_sip | kg-1 s-1 | time, height | sum of secondary ice formation processes (e.g.... |
| 73 | tendency_of_ice_number_due_to_homogeneous_free... | dni_hom | kg-1 s-1 | time, height | ice nucleation source due to homogoeneous free... |
Create DEPHY output file#
Write a single file to contain all domain-mean scalar and profile outputs. This code expects the write directory to be pre-existing (already created by the user). In the case that this output will be committed to the comble-mip GitHub repository, see above “Specify directory locations”.
Definitions for the requested outputs#
# Time variable conversions
def calc_ts(standard_name):
# Return -999 if no data
tmp=xr.DataArray(np.ones(les_ts['psrf'].shape))*-999.
#
# match to LES variable names and determine output
if standard_name=='surface_pressure':
return les_ts['psrf'],''
elif standard_name=='surface_temperature':
return les_ts['tsrf'],''
elif standard_name=='surface_friction_velocity':
return les_ts['ustar'],''
elif standard_name=='surface_roughness_length_for_momentum_in_air':
return les_ts['z0m'],''
elif standard_name=='surface_roughness_length_for_heat_in_air':
if np.all(les_ts['z0t']==-999.):
# The same value (z0m) for momentum, humidity and heat
return les_ts['z0m'],'The same as that for momentum'
else:
return les_ts['z0t'],''
elif standard_name=='surface_roughness_length_for_humidity_in_air':
if np.all(les_ts['z0t']==-999.):
# The same value (z0m) for momentum, humidity and heat
return les_ts['z0m'],'The same as that for momentum'
else:
# The same as that for heat
return les_ts['z0t'],'The same as that for heat'
elif standard_name=='surface_upward_sensible_heat_flux':
return les_ts['shf_bar'],''
elif standard_name=='surface_upward_latent_heat_flux':
return les_ts['lhf_bar'],''
elif standard_name=='obukhov_length':
return les_ts['obl'],''
elif standard_name=='atmosphere_mass_content_of_liquid_cloud_water':
# default breakdown to cloud water and rain in a bulk scheme or radius separated at
# 40 micrometers in a bin scheme; if additional categories are used in a bulk scheme
# (e.g., a drizzle class), provide additional standard_name variable_id (e.g., lwpd);
# all lwp* variables will be summed to obtain total liquid water path
if sb:
return les_ts['lwp_bar'],''
else:
return les_ts['lwp_bar'],'Includes aerosol water' # For SALSA this includes aerosol water
elif standard_name=='atmosphere_mass_content_of_rain_water':
return les_ts['rwp_bar'],''
elif standard_name=='atmosphere_mass_content_of_ice_water':
# all ice-phase hydrometeors
if sb:
return les_ts['iwp_bar']+les_ts['swp_bar']+les_ts['gwp_bar'],''
else:
return les_ts['swp_bar'],''
elif standard_name=='cloud_area_fraction':
# fraction of atmospheric columns with total hydrometeor mid-visible optical thickness > 2;
# may assume geometric scatterers; to be compared with satellite measurements
if 'frac_od' in les_ts.keys():
return les_ts['frac_od'],''
else:
return les_ts['cfrac'],'Based on cloud water' # Original cloud fraction
elif standard_name=='optical_depth':
# scene (all sky); mid-visible, all hydrometeors; may assume geometric scatterers
return les_ts['tau_liq']+les_ts['tau_ice'],'No gases or vapors'
elif standard_name=='optical_depth_of_liquid_cloud':
# as for optical_depth but cloud droplets only
return les_ts['tau_liq'],'No gases or vapors'
elif standard_name=='precipitation_flux_at_surface':
# liquid and ice phase, all hydrometeors
# Unit conversion for prcp: W/m2/(J/kg)=kg/m2/s
# Note: cloud water precipitation ignored for SB
if sb:
return (les_ts['prcp'])/alvl+(les_ts['iprcp']+les_ts['sprcp']+les_ts['gprcp'])/alvi,''
else:
return les_ts['rmH2Ocl']+les_ts['rmH2Opr']+les_ts['rmH2Osn'],''
elif standard_name=='precipitation_flux_at_surface_in_ice_phase':
# ice phase hydrometeors only
if sb:
return (les_ts['iprcp']+les_ts['sprcp']+les_ts['gprcp'])/alvi,''
else:
return les_ts['rmH2Osn'],''
elif standard_name=='toa_outgoing_longwave_flux':
return les_ts['toa_lwu'],''
elif standard_name=='surface_downwelling_longwave_flux':
return les_ts['srf_lwd'],''
elif standard_name=='surface_upwelling_longwave_flux':
return les_ts['srf_lwu'],''
elif standard_name=='surface_sea_spray_number_flux':
# when using prognostic aerosol; emission only (not including dry or wet deposition); total over all aerosol modes
if sb:
return tmp,'NA'
else:
return les_ts['flx_aer'],''
#
print('Undefined ts:'+standard_name)
return tmp,'NA'
# Profile variable conversions
def calc_ps(standard_name):
# Return -999 if no data
tmp=xr.DataArray(np.ones(les_ps['p'].shape))*-999.
#
# SALSA bin limits (lower limits as diameters in nm)
if not sb: rad=np.round(les_ps['B_Rd12a'].values*2e9)
#
# match to LES variable names and determine output
if standard_name=='air_pressure':
return les_ps['p'],''
elif standard_name=='eastward_wind':
return les_ps['u'],''
elif standard_name=='northward_wind':
return les_ps['v'],''
elif standard_name=='air_dry_density':
# Note: UCLALES-SALSA assumes the same densities for dry and moist air
return les_ps['rho_air'],'For dry air'
elif standard_name=='air_temperature':
# Calculate fron theta_L
return (les_ps['thl']+(les_ps['l']+les_ps['rr'])*(alvl/cp))*(les_ps['p']/p00)**(Rd/cp),'Calculated from theta_L'
#return les_ps['T_avg'],''
elif standard_name=='water_vapor_mixing_ratio':
# per kg dry air
return les_ps['rv'],''
elif standard_name=='relative_humidity':
# relative to liquid
return les_ps['SS_avg']/100.+1.,'' # Convert SS to RH/100%
elif standard_name=='relative_humidity_over_ice':
# relative to ice
return les_ps['SSi_avg']/100.+1.,'' # Convert SS to RH/100%
elif standard_name=='air_potential_temperature':
# Calculate fron theta_L
return les_ps['thl']+(les_ps['l']+les_ps['rr'])*(alvl/cp),'Calculated from theta_L'
elif standard_name=='specific_turbulent_kinetic_energy_resolved':
# resolved only
return les_ps['res_tke'],''
elif standard_name=='specific_turbulent_kinetic_energy_sgs':
# SGS only
return les_ps['sfs_tke'],''
elif standard_name=='mass_mixing_ratio_of_cloud_liquid_water_in_air':
# scene (all sky) per kg dry air; default breakdown to cloud water and rain in a bulk scheme
# or radius separated at 40 micrometers in a bin scheme; if additional categories are used
# in a bulk scheme (e.g., a drizzle class), provide additional standard_name variable_id
# (e.g., qld); all ql* variables will be summed to obtain total liquid mixing ratio
if sb:
return les_ps['l'],''
else:
return les_ps['l'],'Includes aerosol water'
elif standard_name=='mass_mixing_ratio_of_rain_water_in_air':
return les_ps['rr'],''
elif standard_name=='mass_mixing_ratio_of_cloud_ice_in_air':
# default breakdown to cloud_ice, snow and graupel; if other categories are used,
# provide additional standard_name and qiX variable_id; all qi* variables will be
# summed to obtain total ice mixing ratio
if sb:
return les_ps['ri'],''
else:
# SALSA has just ice (called snow) so set snow and graupel to NA
return les_ps['rs'],'Called snow in SALSA'
elif standard_name=='mass_mixing_ratio_of_snow_in_air':
if sb:
return les_ps['rs'],''
else:
return tmp,'NA'
elif standard_name=='mass_mixing_ratio_of_graupel_in_air':
if sb:
return les_ps['rg'],''
else:
return tmp,'NA'
elif standard_name=='number_of_liquid_cloud_droplets_in_air':
# scene (all sky) per kg dry air; default breakdown to cloud water and rain;
# if other categories are used, provide additional standard_name and nl* variable_id
return les_ps['tcNct'],''
elif standard_name=='number_of_rain_drops_in_air':
return les_ps['tcNrt'],''
elif standard_name=='number_of_cloud_ice_crystals_in_air':
# default breakdown to cloud_ice, snow and graupel; if other categories are used,
# provide additional standard_name and niX variable_id; all ni* variables will be summed
# to obtain total ice number mixing ratio
if sb:
return les_ps['tcNit'],''
else:
return les_ps['tcNst'],'Called snow in SALSA'
elif standard_name=='number_of_snow_crystals_in_air':
if sb:
return les_ps['tcNst'],''
else:
return tmp,'NA'
elif standard_name=='number_of_graupel_crystals_in_air':
if sb:
return les_ps['tcNgt'],''
else:
return tmp,'NA'
elif standard_name=='number_of_total_aerosol_mode1':
# when using prognostic aerosol; scene (all sky); activated + unactivated: Aitken mode
if sb:
return tmp,'NA'
else:
# SALSA Aitken mode: bins 0-2 (1a aerosol)
aa=les_ps['B_Naa'][:,0:2,:].sum(dim='B_Rd12a') #+les_ps['B_Nca'][:,0:-1,:].sum(dim='B_Rd2ab')
return aa,'SALSA 1a aerosol (d<'+str(rad[3])+' nm)'
elif standard_name=='number_of_total_aerosol_mode2':
# accumulation mode
if sb:
return tmp,'NA'
else:
# SALSA accumulation mode: bins 3-7
# The larger aerosol and cloud a-bins
aa=les_ps['B_Naa'][:,3:7,:].sum(dim='B_Rd12a')+les_ps['B_Nca'][:,0:4,:].sum(dim='B_Rd2ab')
return aa,'SALSA 2a aerosol+cloud ('+str(rad[3])+' nm < d <'+str(rad[8])+'nm)'
elif standard_name=='number_of_total_aerosol_mode3':
# sea spray mode
if sb:
return tmp,'NA'
else:
# SALSA sea spray mode: bins 8-
aa=les_ps['B_Naa'][:,8:,:].sum(dim='B_Rd12a')+les_ps['B_Nca'][:,5:,:].sum(dim='B_Rd2ab')
return aa,'SALSA 2a aerosol+cloud (d>'+str(rad[8])+'nm)'
elif standard_name=='number_of_liquid_cloud_droplets_in_cloud':
# per kg dry air where cloud water > 0.01 g/kg; breakdown to cloud water and rain in
# a bulk scheme or radius separated at 40 micrometers in a bin scheme
return les_ps['icNct'],''
elif standard_name=='number_of_ice_crystals_in_cloud':
# total ice hydrometeors per kg dry air where cloud water > 0.01 g/kg; within liquid
# cloud boundaries in order to check immersion freezing
if sb:
return les_ps['icNit']+les_ps['icNst']+les_ps['icNgt'],''
else:
return les_ps['icNst'],''
elif standard_name=='dissipation_rate_of_turbulent_kinetic_energy':
# total (including numerical diffusion contributions if relevant); report as negative
return les_ps['diss'],''
elif standard_name=='zonal_momentum_flux':
# total (resolved plus SGS)
return les_ps['tot_uw'],''
elif standard_name=='meridional_momentum_flux':
return les_ps['tot_vw'],''
elif standard_name=='variance_of_upward_air_velocity':
return les_ps['w_2'],''
elif standard_name=='vertical_flux_potential_temperature':
# total (resolved plus SGS)
# Not available: les_ps['tot_tw'] is for ice-liquid potential temperature
return tmp,'NA'
elif standard_name=='vertical_flux_liquid_ice_water_potential_temperature':
# total (resolved plus SGS); include sedimentation fluxes
# les_ps['tot_tw'] does not include sedimentation fluxes
return tmp,'NA'
elif standard_name=='vertical_flux_water_vapor':
# total (resolved plus SGS)
# SALSA has water vapor, SB total water
if sb:
return tmp,'NA'
else:
return les_ps['tot_qw'],''
elif standard_name=='vertical_flux_total_water':
# total (resolved plus SGS); vapor+all liquid+all ice; include sedimentation fluxes
if sb and 'hrate' in les_ps.keys():
return les_ps['tot_qw']+(les_ps['crate']+les_ps['rrate'])/alvl+(les_ps['irate']+les_ps['srate']+les_ps['grate']+les_ps['hrate'])/alvi,''
elif sb:
return les_ps['tot_qw']+(les_ps['crate']+les_ps['rrate'])/alvl+(les_ps['irate']+les_ps['srate']+les_ps['grate'])/alvi,''
else:
return tmp,'NA'
elif standard_name=='area_fraction_of_liquid_cloud':
# fraction of cells with cloud water threshold of 0.01 g/kg, not including rain class
return les_ps['frac_ic'],''
elif standard_name=='precipitation_flux_in_air':
# liquid and ice phase, all hydrometeors
# W/m2/(J/kg)=kg/m2/s
if sb:
return (les_ps['crate']+les_ps['rrate'])/alvl+(les_ps['irate']+les_ps['srate']+les_ps['grate'])/alvi,''
else:
return (les_ps['crate']+les_ps['rrate'])/alvl+(les_ps['srate'])/alvi,''
elif standard_name=='precipitation_flux_in_air_in_ice_phase':
if sb:
return (les_ps['irate']+les_ps['srate']+les_ps['grate'])/alvi,''
else:
return les_ps['srate']/alvi,''
elif standard_name=='downwelling_longwave_flux_in_air':
return les_ps['lw_down'],''
elif standard_name=='upwelling_longwave_flux_in_air':
return les_ps['lw_up'],''
elif standard_name=='tendency_of_air_potential_temperature_due_to_radiative_heating':
# There are no theta outputs
return tmp,'NA'
elif standard_name=='tendency_of_air_potential_temperature_due_to_microphysics':
return tmp,'NA'
elif standard_name=='tendency_of_air_potential_temperature_due_to_mixing':
return tmp,'NA'
elif standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_microphysics':
# including net condensation and precipitation
# This is available for SALSA
if sb:
return tmp,'NA'
else:
return les_ps['mcrp_rv'],'Without precipitation'
elif standard_name=='tendency_of_water_vapor_mixing_ratio_due_to_mixing':
# including surface fluxes
if sb:
return tmp,'NA'
else:
return les_ps['advf_rv']+les_ps['diff_rv'],''
elif standard_name=='tendency_of_aerosol_number_due_to_warm_microphysics':
# activated and unactivated aerosol (sum over all modes); losses to liquid-liquid hydrometeor collisions
# + any other liquid-phase sources/sinks that may be accounted for (e.g., from drop breakup);
# this quantity accounts for all microphysical aerosol source terms in a liquid-only simulation
# Separate warm and cold outputs are not available!
return tmp,'NA'
elif standard_name=='tendency_of_aerosol_number_due_to_cold_microphysics':
return tmp,'NA'
elif standard_name=='tendency_of_aerosol_number_due_to_mixing':
# activated and unactivated aerosol (sum over all modes); including surface fluxes (sea spray,
# dry deposition if included, etc.)
if sb:
return tmp,'NA'
else:
return les_ps['advf_na']+les_ps['advf_nc']+les_ps['srfc_na'],'No dry deposition'
elif standard_name=='tendency_of_ice_number_due_to_heterogeneous_freezing':
# sum of primary ice nucleation due to activation of heterogeneous INP; also, in a diagnostic run,
# the ice crystals introduced to meet the diagnostic target
if sb:
# Fixed primary ice
return les_ps['nucl_ni'],''
else:
# Fixed primary ice
return les_ps['nucf_ns'],''
elif standard_name=='tendency_of_ice_number_due_to_secondary_ice_production':
# sum of secondary ice formation processes (e.g., Hallet-Mossop plus any others)
return tmp,'NA' # Switched off
elif standard_name=='tendency_of_ice_number_due_to_homogeneous_freezing':
# ice nucleation source due to homogoeneous freezing
if sb:
return tmp,'NA' # Switched off
else:
return les_ps['nucm_ns'],''
#
print('Undefined ps:'+standard_name+'!')
return tmp,'NA'
Generate output#
# Open source files
les_ts = xr.open_dataset(input_filename_ts)
les_ps = xr.open_dataset(input_filename_ps)
# create DEPHY output file
dephy_filename = my_gitdir + my_outfile + '.nc'
if os.path.exists(dephy_filename):
try:
os.remove(dephy_filename)
print('The file ' + dephy_filename + ' has been deleted successfully')
except Exception as e:
print('Failed to delete '+dephy_filename+': ',e)
dephy_file = Dataset(dephy_filename,mode='w',format='NETCDF3_CLASSIC')
# create global attributes
dephy_file.title='UCLALES-SALSA results for COMBLE-MIP case'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Tomi Raatikainen (tomi.raatikainen@fmi.fi)'
dephy_file.source=exp_name
dephy_file.version=dt.datetime.now().strftime('%Y-%m-%d %H:%M:%S')
dephy_file.format_version='DEPHY SCM format version 1.6'
dephy_file.script='convert_UCLALES-SALSA_output_to_dephy_format.ipynb'
dephy_file.startDate='2020-03-12T22:00:00Z'
# Additional note about microphysics
dephy_file.note=micro
# create dimensions
# zf = dimension, altitude of mid-level points above sea surface = zf
nz = les_ps.sizes['zt']-1
zf = dephy_file.createDimension('zf', nz)
zf = dephy_file.createVariable('zf', np.float64, ('zf',))
zf.units = 'm'
zf.long_name = 'altitude'
zf[:] = les_ps['zt'][1:].data
# ze = dimension, altitude of layer top points above sea surface (used for vertical integrations) = zm
ze = dephy_file.createDimension('ze', nz)
ze = dephy_file.createVariable('ze', np.float64, ('ze',))
ze.units = 'm'
ze.long_name = 'altitude'
ze[:] = les_ps['zm'][1:].data
# time = seconds since 2020-03-12 18:00:00 - Note: here the start time is 2020-03-12 22:00:00
nt = les_ps.sizes['time']
time = dephy_file.createDimension('time', nt)
time = dephy_file.createVariable('time', np.float64, ('time',))
time.units = 'seconds since ' + dephy_file.startDate
time.long_name = 'time'
time[:] = 5.*np.floor(les_ps['time'].data/5.) # Model time may exceed target time by 1 s
# create variables
for index in vars_mean_list.index[2:]:
# Output
std_name = vars_mean_list.standard_name.iat[index]
var_name = vars_mean_list.variable_id.iat[index]
if vars_mean_list.dimensions.iat[index]=='time':
# 1D time series
new_sca = dephy_file.createVariable(var_name, np.float64, ('time'))
new_sca.units = vars_mean_list.units.iat[index]
new_sca.long_name = std_name
try:
out,note=calc_ts(std_name)
if note!='NA' and not np.all(out==-999.): # No data: do not fill the variable field
# Current ts-data is with 2 min output frequency => needs to be averaged
if np.any(out==-999.): print('Warning: averaging missing values for ',std_name)
time0=time[0]
for i in range(nt):
new_sca[i]=out.sel(time=slice(time0,time[i]+1.5)).mean('time')
time0=time[i]+1.5
if len(note)>0: new_sca.note=note
except Exception as e:
print('Failed: '+std_name,e)
elif vars_mean_list.dimensions.iat[index]=='time, height':
# 2D (z,time)
new_snd = dephy_file.createVariable(var_name, np.float64, ('time','zf'))
new_snd.units = vars_mean_list.units.iat[index]
new_snd.long_name = std_name
try:
arr,note = calc_ps(std_name)
if note!='NA' and not np.all(arr==-999.): # No data: do not fill the variable field
if np.any(arr==-999.): print('Warning: including missing values for ',std_name)
new_snd[:]=arr[:,1:] # Without z<0
if len(note)>0: new_snd.note=note
except Exception as e:
print('Failed: '+std_name,e)
print(dephy_file)
# Close all files
dephy_file.close()
les_ts.close()
les_ps.close()
The file ../../output_les/uclales-salsa/sandbox/UCLALES-SALSA_Lx25_dx100_ProgNc.nc has been deleted successfully
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
title: UCLALES-SALSA results for COMBLE-MIP case
reference: https://github.com/ARM-Development/comble-mip
authors: Tomi Raatikainen (tomi.raatikainen@fmi.fi)
source: COMBLE_SALSA_ProgNc
version: 2024-08-07 09:02:16
format_version: DEPHY SCM format version 1.6
script: convert_UCLALES-SALSA_output_to_dephy_format.ipynb
startDate: 2020-03-12T22:00:00Z
note: UCLALES-SALSA with prognostic sectional aerosol, cloud and ice (called snow) categories.
dimensions(sizes): zf(159), ze(159), time(121)
variables(dimensions): float64 zf(zf), float64 ze(ze), float64 time(time), float64 ps(time), float64 ts(time), float64 ustar(time), float64 z0(time), float64 z0h(time), float64 z0q(time), float64 hfss(time), float64 hfls(time), float64 ol(time), float64 lwpc(time), float64 lwpr(time), float64 iwp(time), float64 clt(time), float64 od(time), float64 odlc(time), float64 pr(time), float64 pri(time), float64 rlut(time), float64 rlds(time), float64 rlus(time), float64 ssaf(time), float64 pa(time, zf), float64 ua(time, zf), float64 va(time, zf), float64 rhoa(time, zf), float64 ta(time, zf), float64 qv(time, zf), float64 hur(time, zf), float64 huri(time, zf), float64 theta(time, zf), float64 tke_res(time, zf), float64 tke_sgs(time, zf), float64 qlc(time, zf), float64 qlr(time, zf), float64 qic(time, zf), float64 qis(time, zf), float64 qig(time, zf), float64 nlc(time, zf), float64 nlr(time, zf), float64 nic(time, zf), float64 nis(time, zf), float64 nig(time, zf), float64 na1(time, zf), float64 na2(time, zf), float64 na3(time, zf), float64 nlcic(time, zf), float64 niic(time, zf), float64 eps(time, zf), float64 uw(time, zf), float64 vw(time, zf), float64 w2(time, zf), float64 wth(time, zf), float64 vf_thli(time, zf), float64 wqv(time, zf), float64 vf_qt(time, zf), float64 flc(time, zf), float64 prf(time, zf), float64 prfi(time, zf), float64 rld(time, zf), float64 rlu(time, zf), float64 dth_rad(time, zf), float64 dth_micro(time, zf), float64 dth_turb(time, zf), float64 dq_micro(time, zf), float64 dq_turb(time, zf), float64 dna_micro_warm(time, zf), float64 dna_micro_cold(time, zf), float64 dna_turb(time, zf), float64 dni_het(time, zf), float64 dni_sip(time, zf), float64 dni_hom(time, zf)
groups:
Generate 2D and 3D outputs#
2D#
# read list of requested variables
vars_2d_list = pd.read_excel('https://docs.google.com/spreadsheets/d/1Vl8jYGviet7EtXZuQiitrx4NSkV1x27aJAhxxjBb9zI/export?format=xlsx&gid=1756539842')
vars_2d_list
| standard_name | variable_id | units | dimensions | comment (10-min instantaneous) | |
|---|---|---|---|---|---|
| 0 | time | time | s | – | dimension, seconds since 2020-03-12 18:00:00 |
| 1 | x | x | m | – | dimension |
| 2 | y | y | m | – | dimension |
| 3 | surface_upward_sensible_heat_flux | hfss | W m-2 | time, x, y | – |
| 4 | surface_upward_latent_heat_flux | hfls | W m-2 | time, x, y | – |
| 5 | surface_friction_velocity | ustar | m s-1 | time, x, y | – |
| 6 | surface_eastward_wind | us | m s-1 | time, x, y | at 10-m for comparison with SAR satellite meas... |
| 7 | surface_northward_wind | vs | m s-1 | time, x, y | at 10-m for comparison with SAR satellite meas... |
| 8 | precipitation_flux_at_surface | pr | kg m-2 s-1 | time, x, y | all hydrometeors |
| 9 | atmosphere_mass_content_of_liquid_water | lwp | kg m-2 | time, x, y | all liquid hydrometeors |
| 10 | atmosphere_mass_content_of_ice | iwp | kg m-2 | time, x, y | all ice hydrometeors |
| 11 | atmosphere_optical_thickness | opt | 1 | time, x, y | all hydrometeors for comparison with VIIRS; mi... |
| 12 | pseudo-albedo | alb | 1 | time, x, y | opt/(opt+13), where opt = all hydrometeor opti... |
| 13 | toa_outgoing_longwave_flux_atmospheric_window | olr11 | K | time, x, y | flux at top-of-atmosphere in window containing... |
# Column (2D) variable conversions
# Dimensions
# time (s), seconds since 2020-03-12 18:00:00
# x, x (m)
# y, y (m)
def calc_cs(standard_name):
# Return -999 if no data
tmp=xr.DataArray(np.ones(les_cs['shf_bar'].shape))*-999.
#
# match to LES variable names and determine output
if standard_name=='surface_upward_sensible_heat_flux':
return les_cs['shf_bar'],''
elif standard_name=='surface_upward_latent_heat_flux':
return les_cs['lhf_bar'],''
elif standard_name=='surface_friction_velocity':
return les_cs['ustar'],''
elif standard_name=='surface_eastward_wind':
return les_cs['us'],''
elif standard_name=='surface_northward_wind':
return les_cs['vs'],''
elif standard_name=='precipitation_flux_at_surface':
# all hydrometeors
if sb:
return (les_cs['rprcp'])/alvl+(les_cs['iprcp']+les_cs['sprcp']+les_cs['gprcp'])/alvi,''
else:
return (les_cs['cprcp']+les_cs['rprcp'])/alvl+les_cs['sprcp']/alvi,''
elif standard_name=='atmosphere_mass_content_of_liquid_water':
# all liquid hydrometeors
return les_cs['lwp']+les_cs['rwp'],''
elif standard_name=='atmosphere_mass_content_of_ice':
if sb:
return les_cs['iwp']+les_cs['swp']+les_cs['gwp'],''
else:
return les_cs['swp'],''
elif standard_name=='atmosphere_optical_thickness':
# all hydrometeors for comparison with VIIRS; mid-visible, may assume geometric scatterers
return les_cs['tau_liq']+les_cs['tau_ice'],'No gases or vapors'
elif standard_name=='pseudo-albedo':
# opt/(opt+13), where opt = all hydrometeor optical depth as reported in domain-mean
return (les_cs['tau_liq']+les_cs['tau_ice'])/(les_cs['tau_liq']+les_cs['tau_ice']+13),''
elif standard_name=='toa_outgoing_longwave_flux_atmospheric_window':
# flux at top-of-atmosphere in window containing 10.8 µm, please specify wavelength range
return les_cs['toa_lwu'],'Total LW'
#
print('Undefined cs:'+standard_name+'!')
return tmp,'NA'
Generate outputs#
# Open source files
if not os.path.exists(input_filename_cs): raise Exception('Input file not found!')
les_cs = xr.open_dataset(input_filename_cs)
# create DEPHY output file
dephy_filename = my_gitdir + my_outfile + '_2D.nc'
if os.path.exists(dephy_filename):
try:
os.remove(dephy_filename)
print('The file '+dephy_filename+' has been deleted successfully')
except Exception as e:
print('Failed to delete '+dephy_filename+': ',e)
dephy_file = Dataset(dephy_filename,mode='w',format='NETCDF3_CLASSIC')
# create global attributes
dephy_file.title='UCLALES-SALSA results for COMBLE-MIP case'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Tomi Raatikainen (tomi.raatikainen@fmi.fi)'
dephy_file.source=exp_name
dephy_file.version=dt.datetime.now().strftime('%Y-%m-%d %H:%M:%S')
dephy_file.format_version='DEPHY SCM format version 1.6'
dephy_file.script='convert_UCLALES-SALSA_output_to_dephy_format.ipynb'
dephy_file.startDate='2020-03-12T22:00:00Z'
# Additional note about microphysics
dephy_file.note=micro
# create dimensions
# x, y
nx = les_cs.sizes['xt']
x = dephy_file.createDimension('x', nx)
x = dephy_file.createVariable('x', np.float64, ('x',))
x.units = 'm'
x.long_name = 'x'
x[:] = les_cs['xt'].data
ny = les_cs.sizes['yt']
y = dephy_file.createDimension('y', ny)
y = dephy_file.createVariable('y', np.float64, ('y',))
y.units = 'm'
y.long_name = 'y'
y[:] = les_cs['yt'].data
# time
t10min=True #False
if t10min:
# 10 min resolution
nt = les_cs.sizes['time']
time = dephy_file.createDimension('time', nt)
time = dephy_file.createVariable('time', np.float64, ('time',))
time.units = 'seconds since ' + dephy_file.startDate
time.long_name = 'time'
time[:] = 5.*np.floor(les_cs['time'].data/5.) # Model time may exceed target time by 1 s
else:
# 60 min resolution
nt = (les_cs.sizes['time']-1)/6+1
time = dephy_file.createDimension('time', nt)
time = dephy_file.createVariable('time', np.float64, ('time',))
time.units = 'seconds since ' + dephy_file.startDate
time.long_name = 'time'
time[:] = 5.*np.floor(les_cs['time'][::6].data/5.) # Model time may exceed target time by 1 s
preci = np.float64 # Default
#preci = np.float32 # To save disc space
# create variables
for index in vars_2d_list.index[3:]:
# Output
std_name = vars_2d_list.iat[index,0]
var_name = vars_2d_list.iat[index,1]
new = dephy_file.createVariable(var_name, preci, ('time','x','y'))
new.units = vars_2d_list.units.iat[index]
new.long_name = std_name
try:
dat,note=calc_cs(std_name)
# NA means no data
if note=='NA':
print('Ignoring '+std_name+' (no data)')
else:
if t10min:
new[:]=dat # 10 min resolution
else:
new[:]=dat[::6,:,:] # 60 min resolution
if len(note)>0: new.note=note
except Exception as e:
print('Failed: '+std_name,e)
print(dephy_file)
# Close all files
dephy_file.close()
les_cs.close()
Ignoring toa_outgoing_longwave_flux_atmospheric_window (no data)
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
title: UCLALES-SALSA results for COMBLE-MIP case
reference: https://github.com/ARM-Development/comble-mip
authors: Tomi Raatikainen (tomi.raatikainen@fmi.fi)
source: COMBLE_SB_FixN
version: 2024-03-15 09:59:49
format_version: DEPHY SCM format version 1.6
script: convert_UCLALES-SALSA_output_to_dephy_format.ipynb
startDate: 2020-03-12T22:00:00Z
note: Without SALSA, i.e., using the 2M bulk microphysics by Seifert and Beheng (SB). This version has diagnostic clouds (saturation adjustment + fixed CDNC) and prognostic 2M rain, ice, snow and graupel categories.
dimensions(sizes): x(256), y(256), time(121)
variables(dimensions): float64 x(x), float64 y(y), float64 time(time), float32 hfss(time, x, y), float32 hfls(time, x, y), float32 ustar(time, x, y), float32 us(time, x, y), float32 vs(time, x, y), float32 pr(time, x, y), float32 lwp(time, x, y), float32 iwp(time, x, y), float32 opt(time, x, y), float32 alb(time, x, y)
groups:
3D#
# read list of requested variables
vars_3d_list = pd.read_excel('https://docs.google.com/spreadsheets/d/1Vl8jYGviet7EtXZuQiitrx4NSkV1x27aJAhxxjBb9zI/export?format=xlsx&gid=1233994833')
vars_3d_list
| standard_name | variable_id | units | dimensions | comment (60-min instantaneous, red=required for EMC2) | |
|---|---|---|---|---|---|
| 0 | time | time | s | – | dimension, seconds since 2020-03-12 18:00:00 |
| 1 | x | x | m | – | dimension |
| 2 | y | y | m | – | dimension |
| 3 | height | zf | m | – | dimension, altitude of mid-level points above ... |
| 4 | air_pressure | pa | Pa | time, height | – |
| 5 | eastward_wind | ua | m s-1 | time, height, x, y | – |
| 6 | northward_wind | va | m s-1 | time, height, x, y | – |
| 7 | upward_air_velocity | wa | m s-1 | time, height, x, y | – |
| 8 | air_temperature | ta | K | time, height, x, y | – |
| 9 | air_dry_density | rhoa | kg m-3 | time, height, x, y | – |
| 10 | water_vapor_mixing_ratio | qv | kg kg-1 | time, height, x, y | per kg dry air |
| 11 | mass_mixing_ratio_of_cloud_liquid_water_in_air | qlc | kg kg-1 | time, height, x, y | per kg dry air; default breakdown to cloud wat... |
| 12 | mass_mixing_ratio_of_rain_water_in_air | qlr | kg kg-1 | time, height, x, y | – |
| 13 | mass_mixing_ratio_of_cloud_ice_in_air | qic | kg kg-1 | time, height, x, y | default breakdown to cloud_ice, snow and graup... |
| 14 | mass_mixing_ratio_of_snow_in_air | qis | kg kg-1 | time, height, x, y | – |
| 15 | mass_mixing_ratio_of_graupel_in_air | qig | kg kg-1 | time, height, x, y | – |
| 16 | number_of_dry_aerosol_in_air | na | kg-1 | time, height, x, y | per kg dry air; when using prognostic aerosol |
| 17 | number_of_cloud_droplets_in_air | nlc | kg-1 | time, height, x, y | when using prognostic aerosol |
| 18 | number_of_rain_droplets_in_air | nlr | kg-1 | time, height, x, y | – |
| 19 | number_of_cloud_ice_crystals_in_air | nic | kg-1 | time, height, x, y | default breakdown to cloud_ice, snow and graup... |
| 20 | number_of_snow_crystals_in_air | nis | kg-1 | time, height, x, y | – |
| 21 | number_of_graupel_crystals_in_air | nig | kg-1 | time, height, x, y | – |
| 22 | dissipation_rate_of_turbulent_kinetic_energy | eps | m2 s-3 | time, height, x, y | including numerical diffusion contributions if... |
| 23 | mass_weighted_fall_speed_of_liquid_cloud_water | vmlc | m s-1 | time, height, x, y | EMC2 uses mass-weighted fall-speed profiles fo... |
| 24 | mass_weighted_fall_speed_of_rain | vmlr | m s-1 | time, height, x, y | – |
| 25 | mass_weighted_fall_speed_of_ice_cloud | vmic | m s-1 | time, height, x, y | if other ice-phase categories are used, provid... |
| 26 | mass_weighted_fall_speed_of_snow | vmis | m s-1 | time, height, x, y | – |
| 27 | mass_weighted_fall_speed_of_graupel | vmig | m s-1 | time, height, x, y | – |
| 28 | tendency_of_condensation | tend_cond | s-1 | time, height, x, y | combined tendencies of cloud liquid and rain w... |
| 29 | tendency_of_deposition | tend_dep | s-1 | time, height, x, y | tendency of ice mixing ratio due to deposition... |
| 30 | tendency_of_riming | tend_rim | s-1 | time, height, x, y | tendency of ice mixing ratio due to riming; su... |
# Analysis (3D) variable conversions
# Dimensions
# time (s), seconds since 2020-03-12 18:00:00
# x, x (m)
# y, y (m)
# height, zf (m), altitude of mid-level points above sea surface
def calc_an(standard_name):
# Return -999 if no data
tmp=xr.DataArray(np.ones(les_an['p'].shape))*-999.
#
# match to LES variable names and determine output
if standard_name=='air_pressure':
return les_an['p'],''
elif standard_name=='eastward_wind':
# Note: need to interpolate for cell centers!
return les_an['u'],''
elif standard_name=='northward_wind':
# Note: need to interpolate for cell centers!
return les_an['v'],''
elif standard_name=='upward_air_velocity':
# Note: need to interpolate for cell centers!
return les_an['w'],''
elif standard_name=='air_temperature':
# temp=theta*(press/p00)^(Rd/cp)
return les_an['theta']*(les_an['p']/p00)**(Rd/cp),'Calculated'
elif standard_name=='air_dry_density':
# rho=press/(Rd*temp)
return les_an['p']/(Rd*les_an['theta']*(les_an['p']/p00)**(Rd/cp)),'Calculated'
elif standard_name=='water_vapor_mixing_ratio':
# per kg dry air
if sb:
return les_an['q']-les_an['l']-les_an['r']-les_an['i']-les_an['s']-les_an['g'],'Calculated'
else:
return les_an['q']-les_an['l']-les_an['r']-les_an['s'],'Calculated'
elif standard_name=='mass_mixing_ratio_of_cloud_liquid_water_in_air':
# per kg dry air; default breakdown to cloud water and rain; if other categories are used,
# provide additional standard_name and ql* variable_id; all ql* variables will be summed
# to obtain total liquid mixing ratio
if sb:
return les_an['l'],''
else:
return les_an['l'],'Includes aerosol water'
elif standard_name=='mass_mixing_ratio_of_rain_water_in_air':
return les_an['r'],''
elif standard_name=='mass_mixing_ratio_of_cloud_ice_in_air':
# default breakdown to cloud_ice, snow and graupel; if other categories are used,
# provide additional standard_name and qiX variable_id; all qi* variables will be summed
# to obtain total ice mixing ratio
if sb:
return les_an['i'],''
else:
# SALSA has just snow
return les_an['s'],'Called snow in SALSA'
elif standard_name=='mass_mixing_ratio_of_snow_in_air':
if sb:
return les_an['s'],''
else:
return tmp,'NA'
elif standard_name=='mass_mixing_ratio_of_graupel_in_air':
if sb:
return les_an['g'],''
else:
return tmp,'NA'
elif standard_name=='number_of_dry_aerosol_in_air':
# per kg dry air; when using prognostic aerosol
if sb:
return tmp,'NA'
else:
return les_an['C_Naa'],''
elif standard_name=='number_of_cloud_droplets_in_air':
# when using prognostic aerosol
if sb:
return tmp,'NA'
else:
return les_an['C_Nca'],''
elif standard_name=='number_of_rain_droplets_in_air':
if sb:
return les_an['n'],''
else:
return les_an['C_Nrt'],''
elif standard_name=='number_of_cloud_ice_crystals_in_air':
# default breakdown to cloud_ice, snow and graupel; if other categories are used,
# provide additional standard_name and niX variable_id; all ni* variables will be
# summed to obtain total ice number mixing ratio
if sb:
return les_an['ni'],''
else:
return les_an['C_Nst'],'Called snow in SALSA'
elif standard_name=='number_of_snow_crystals_in_air':
if sb:
return les_an['ns'],''
else:
return tmp,'NA'
elif standard_name=='number_of_graupel_crystals_in_air':
if sb:
return les_an['ng'],''
else:
return tmp,'NA'
elif standard_name=='dissipation_rate_of_turbulent_kinetic_energy':
# including numerical diffusion contributions if relevant; report as negative
return les_an['diss'],''
elif standard_name=='mass_weighted_fall_speed_of_liquid_cloud_water':
# Curretly not available
return tmp,'NA'
elif standard_name=='mass_weighted_fall_speed_of_rain':
return tmp,'NA'
elif standard_name=='mass_weighted_fall_speed_of_ice_cloud':
return tmp,'NA'
elif standard_name=='mass_weighted_fall_speed_of_snow':
return tmp,'NA'
elif standard_name=='mass_weighted_fall_speed_of_graupel':
return tmp,'NA'
elif standard_name=='tendency_of_condensation':
# combined tendencies of cloud liquid and rain water mixing ratios due to
# condensation/evaporation (+ for condensation, - for evaporation)
if sb:
return les_an['cond_rr'],'Rain only'
else:
return les_an['cond_ra']+les_an['cond_rc']+les_an['cond_rr'],'Aerosol, cloud and rain'
elif standard_name=='tendency_of_deposition':
# tendency of ice mixing ratio due to deposition/sublimation; sum over all
# ice hydrometeor classes (+ for deposition, - for sublimation)
if sb:
return les_an['cond_ri']+les_an['cond_rs']+les_an['cond_rg'],''
else:
return les_an['cond_rs'],''
elif standard_name=='tendency_of_riming':
# tendency of ice mixing ratio due to riming; sum over all ice hydrometeor classes
if sb:
return -les_an['rimi_rc']-les_an['rimi_rr'],'Cloud and rain'
else:
return les_an['coag_rs'],''
#
print('Undefined an:'+standard_name+'!')
return tmp,'NA'
Generate outputs#
# Open source files
if not os.path.exists(input_filename_an): raise Exception('Input file not found!')
les_an = xr.open_dataset(input_filename_an)
# create DEPHY output file
dephy_filename = my_gitdir + my_outfile + '_3D.nc'
if os.path.exists(dephy_filename):
try:
os.remove(dephy_filename)
print('The file '+dephy_filename+' has been deleted successfully')
except Exception as e:
print('Failed to delete '+dephy_filename+': ',e)
dephy_file = Dataset(dephy_filename,mode='w',format='NETCDF4')
# Note 1: NETCDF3_CLASSIC is not valid for 2+ GB files
# Note 2: even with float32 and zlib=True the output is at least 7 GB
# create global attributes
dephy_file.title='UCLALES-SALSA results for COMBLE-MIP case'
dephy_file.reference='https://github.com/ARM-Development/comble-mip'
dephy_file.authors='Tomi Raatikainen (tomi.raatikainen@fmi.fi)'
dephy_file.source=exp_name
dephy_file.version=dt.datetime.now().strftime('%Y-%m-%d %H:%M:%S')
dephy_file.format_version='DEPHY SCM format version 1.6'
dephy_file.script='convert_UCLALES-SALSA_output_to_dephy_format.ipynb'
dephy_file.startDate='2020-03-12T22:00:00Z'
# Additional note about microphysics
dephy_file.note=micro
# create dimensions
# x, y
nx = les_an.sizes['xt']
x = dephy_file.createDimension('x', nx)
x = dephy_file.createVariable('x', np.float64, ('x',))
x.units = 'm'
x.long_name = 'x'
x[:] = les_an['xt'].data
ny = les_an.sizes['yt']
y = dephy_file.createDimension('y', ny)
y = dephy_file.createVariable('y', np.float64, ('y',))
y.units = 'm'
y.long_name = 'y'
y[:] = les_an['yt'].data
# zf = dimension, altitude of mid-level points above sea surface = zf
nz = les_an.sizes['zt']-1
zf = dephy_file.createDimension('zf', nz)
zf = dephy_file.createVariable('zf', np.float64, ('zf',))
zf.units = 'm'
zf.long_name = 'altitude'
zf[:] = les_an['zt'][1:].data
# time
nt = les_an.sizes['time']
time = dephy_file.createDimension('time', nt)
time = dephy_file.createVariable('time', np.float64, ('time',))
time.units = 'seconds since ' + dephy_file.startDate
time.long_name = 'time'
time[:] = 5*np.floor(les_an['time'].data/5.) # Model time may exceed target time by 1 s
preci = np.float64 # Default
#preci = np.float32 # To save disc space
# create variables
for index in vars_3d_list.index[4:]:
# Output
std_name = vars_3d_list.iat[index,0]
var_name = vars_3d_list.iat[index,1]
new = dephy_file.createVariable(var_name, preci, ('time','x','y','zf'), zlib=True)
new.units = vars_3d_list.units.iat[index]
new.long_name = std_name
try:
dat,note=calc_an(std_name)
# NA means no data
if note=='NA':
print('Ignoring '+std_name+' (no data)')
else:
new[:]=dat[:,:,:,1:] # Without z<0
if len(note)>0: new.note=note
except Exception as e:
print('Failed: '+std_name,e)
print(dephy_file)
# Close all files
dephy_file.close()
les_an.close()
---------------------------------------------------------------------------
Exception Traceback (most recent call last)
Cell In[13], line 2
1 # Open source files
----> 2 if not os.path.exists(input_filename_an): raise Exception('File not found!')
3 les_an = xr.open_dataset(input_filename_an)
5 # create DEPHY output file
Exception: File not found!