Simulation with ice part 1: FixN

import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import netCDF4
import datetime
from netCDF4 import Dataset
import re

native_outputs_path='/proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_NEW_SIMULATIONS/'
converted_outputs_path='/proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED/'
converted_outputs_path_time_mod='/proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED_time_mod/'

# 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#.iloc[:, :5]
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 22: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... 
# 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#.iloc[:, :5]
standard_name variable_id units dimensions comment (10-min instantaneous)
0 time time s – dimension, seconds since 2020-03-12 22:00:00
1 x x m – dimension, grid cell center
2 y y m – dimension, grid cell center
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... 
vars_3d_list = pd.read_excel('https://docs.google.com/spreadsheets/d/1Vl8jYGviet7EtXZuQiitrx4NSkV1x27aJAhxxjBb9zI/export?format=xlsx&gid=1233994833')
vars_3d_list#.iloc[:, :5]
standard_name variable_id units dimensions comment (60-min instantaneous, red=required for EMC2)
0 time time s – dimension, seconds since 2020-03-12 22:00:00
1 x x m – dimension, grid cell center
2 y y m – dimension, grid cell center
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 tendency_of_condensation tend_cond s-1 time, height, x, y combined tendencies of cloud liquid and rain w...
24 tendency_of_deposition tend_dep s-1 time, height, x, y tendency of ice mixing ratio due to deposition...
25 tendency_of_riming tend_rim s-1 time, height, x, y tendency of ice mixing ratio due to riming; su... 
Output_Mean = {row['variable_id']: {
        'std_name': str(row['standard_name']) if pd.notna(row['standard_name']) else '',
        'var_id': str(row['variable_id']),
        'units': str(row['units']) if pd.notna(row['units']) else '',
        'dim': str(row['dimensions']) if pd.notna(row['dimensions']) else ''
    }
    for _, row in vars_mean_list.iterrows()
    if pd.notna(row['variable_id'])
}

Output_2D = {row['variable_id']: {
        'std_name': str(row['standard_name']) if pd.notna(row['standard_name']) else '',
        'var_id': str(row['variable_id']),
        'units': str(row['units']) if pd.notna(row['units']) else '',
        'dim': str(row['dimensions']) if pd.notna(row['dimensions']) else ''
    }
    for _, row in vars_2d_list.iterrows()
    if pd.notna(row['variable_id'])
}

Output_3D = {row['variable_id']: {
        'std_name': str(row['standard_name']) if pd.notna(row['standard_name']) else '',
        'var_id': str(row['variable_id']),
        'units': str(row['units']) if pd.notna(row['units']) else '',
        'dim': str(row['dimensions']) if pd.notna(row['dimensions']) else ''
    }
    for _, row in vars_3d_list.iterrows()
    if pd.notna(row['variable_id'])
}

def Output_add_attributes(output_data,output_name,mixed_phase):
    
    if mixed_phase:

        Output_Mean_attr={'title': 'Mean output for '+output_name,
                            'reference': 'COMBLE',
                            'author': 'Alejandro Baro Perez (alejandro.baro.perez@misu.su.se)',
                            'version': 'Created on 2026-06-18',
                            'format_version': 'DEPHY SCM format version 1',}
    
    else:
        
        Output_Mean_attr={'case': 'COMBLE_CAO','title': 'Mean output for '+output_name,
                            'reference': 'COMBLE',
                            'author': 'Alejandro Baro Perez (alejandro.baro.perez@misu.su.se)',
                            'version': 'Created on 2026-06-18',
                            'format_version': 'DEPHY SCM format version 1',}        
        
    
    output_data.attrs=Output_Mean_attr
    
    return(output_data)


def convert_time_new(output_xarray,yyear,mmonth,dday,hhour,time_format):
    
    time_arr=output_xarray.time.values
    
    if(time_format=='days'):
        t_origin = datetime.datetime(year=yyear, month=mmonth, day=dday, hour=hhour)
        date_list = [datetime.timedelta(days=x) + t_origin for x in time_arr[:].astype(np.float64)]
    if(time_format=='seconds'):
        t_origin = datetime.datetime(year=yyear, month=mmonth, day=dday, hour=hhour)
        date_list = [datetime.timedelta(seconds=x) + t_origin for x in time_arr[:].astype(np.float64)]
    
    # Set seconds and microseconds to zero
    date_list = [dt.replace(second=0, microsecond=0) for dt in date_list]
    
    date_list=np.asarray(date_list)
    
    output_xarray=output_xarray.assign(time=np.asarray(date_list))

    return(output_xarray)




#=============================================================================================================

# OUTPUT STATISTIC
# domain-mean quantities and statistics every 10 minutes (time-averages reported at the end of each period)

def create_stat_output(MIMICA_prof,MIMICA_prof_cl,MIMICA_surf,MIMICA_TS,output_name,mixed_phase):
    
    #Ensuring that the time values are the same for prof,prof_cl,surf and TS outputs
    MIMICA_TS=MIMICA_TS.assign_coords(time=MIMICA_surf.time.values)
    MIMICA_prof=MIMICA_prof.assign_coords(time=MIMICA_surf.time.values)
    MIMICA_prof_cl=MIMICA_prof_cl.assign_coords(time=MIMICA_surf.time.values)

    Statistical_output = xr.Dataset()
    
    #Adding the surface pressure (it is constant in MIMICA)
    psurf=99544.5 #[Pa]
    ps_array=np.full(MIMICA_surf.time.shape,psurf)
    ps_dataset = xr.Dataset({"ps": ("time",ps_array,)},coords={"time": MIMICA_surf.time.values},)   
    
    Statistical_output['ps']=ps_dataset.ps
    Statistical_output.ps.attrs['standard_name']='surface_pressure'
    Statistical_output.ps.attrs['units']='Pa'        
    
    # Friction velocity, surface fluxes
    DEPHY_name=['ts','ustar', 'z0','z0h','hfss','hfls', 'ol','rlut', 'rlus'  , 'rlds']
    MIMICA_name=['SST','U_star', 'Z0_m','z0_h','SHF','LHF','L_MO','OLR', 'Surface_LW_up','Surface_LW_down']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        Statistical_output[dephy]=MIMICA_TS[mimica]
        Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name']   
        Statistical_output[dephy].attrs['units']=Output_Mean[dephy]['units'] 
 
     # Water paths
    DEPHY_name=['lwpc','lwpr','iwp']
    MIMICA_name=['LWPc','LWPr', 'IWP']  
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        Statistical_output[dephy]=MIMICA_TS[mimica]
        Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name']  
        Statistical_output[dephy].attrs['units']=Output_Mean[dephy]['units'] 
 
    #Cloud area fraction
    Statistical_output['clt']=MIMICA_TS['C_COL']
    Statistical_output['clt'].attrs['standard_name']= Output_Mean['clt']['std_name'] 
    Statistical_output['clt'].attrs['units']=Output_Mean['clt']['units']     
         
    #Optical depths
    Statistical_output['od'] =(MIMICA_surf['OPTC']+MIMICA_surf['OPTI']+MIMICA_surf['OPTR']).mean(dim=('Y','X'))
    Statistical_output['od'].attrs['standard_name']= Output_Mean['od']['std_name'] 
    Statistical_output['od'].attrs['units']=Output_Mean['od']['units']
    
    Statistical_output['odlc'] =MIMICA_surf['OPTC'].mean(dim=('Y','X'))
    Statistical_output['odlc'].attrs['standard_name']= Output_Mean['odlc']['std_name']        
    Statistical_output['odlc'].attrs['units']=Output_Mean['odlc']['units']
    
    # Surface flux of aerosols MIMICA
    if 'Surf_naer_flux' in MIMICA_surf:
        Statistical_output['ssaf']=MIMICA_surf['Surf_naer_flux'].mean(dim=('Y','X'))
        Statistical_output['ssaf'].attrs['standard_name']= Output_Mean['ssaf']['std_name']                                                                       
        Statistical_output['ssaf'].attrs['units']= Output_Mean['ssaf']['units']
    
    
    # U,V,temperature,potential temperature,water vapor mixing fraction(ratio)
    DEPHY_name=['pa', 'ua', 'va','ta','theta','tke_res','tke_sgs']
    MIMICA_name=['P_tot','U', 'V', 'T','PT'  , 'TKE' ,  'TKE_sgs']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        Statistical_output[dephy]=MIMICA_prof[mimica]
        Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name']  
        Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units']    

    # Dry air density is calculated from MIMICA full density and the total water mixing ratio :
    RHO_dry= MIMICA_prof['RHO']/(1+MIMICA_prof['Qt'])
    Statistical_output['rhoa']=RHO_dry
    Statistical_output['rhoa'].attrs['standard_name']= Output_Mean['rhoa']['std_name']  
    Statistical_output['rhoa'].attrs['units']= Output_Mean['rhoa']['units']

    # Vertical velocity variance
    Statistical_output['w2']=MIMICA_prof['W_var']
    Statistical_output['w2'].attrs['standard_name']= Output_Mean['w2']['std_name']  
    Statistical_output['w2'].attrs['units']= Output_Mean['w2']['units']
        
    ## Precipitation flux at surface. Change from precipitation rate (MIMICA) to precipitation flux
    #Considering that
    water_factor= 1./86400
    ice_factor= 1./86400
    

    DEPHY_name=['pr','pri']
    MIMICA_name=['R_rate','R_ice_rate']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mixed_phase:
            #Assuming that the liquid phase precipitation in negligible in the mixed phase case
            Statistical_output['pri']=ice_factor*MIMICA_TS['R_ice_rate']
            Statistical_output['pri'].attrs['standard_name']= Output_Mean[dephy]['std_name']
            Statistical_output['pri'].attrs['units']= Output_Mean[dephy]['units']

        Statistical_output['pr']=water_factor*MIMICA_TS['R_rate']                    
        Statistical_output['pr'].attrs['standard_name']= Output_Mean[dephy]['std_name']
        Statistical_output['pr'].attrs['units']= Output_Mean[dephy]['units']        

    #Precipitation flux in air in ice phase (considering the fall speed of the hydrometeors and the mixing ratio of the hydrometeors in air)
    if mixed_phase:
        Statistical_output['prfi']=MIMICA_prof['Qi']*MIMICA_prof['Vpi']+MIMICA_prof['Qs']*MIMICA_prof['Vps']+MIMICA_prof['Qg']*MIMICA_prof['Vpg']
        Statistical_output['prfi'].attrs['standard_name']= Output_Mean['prfi']['std_name'] 
        Statistical_output['prfi'].attrs['units']= Output_Mean['prfi']['units']

        #Precipitation flux in air (considering the fall speed of the hydrometeors and the mixing ratio of the hydrometeors in air) 
        Statistical_output['prf']=Statistical_output['prfi']+ MIMICA_prof['Vpc']*MIMICA_prof['Qc']+MIMICA_prof['Qr']*MIMICA_prof['Vpr']   
    else:
        Statistical_output['prf']=MIMICA_prof['Vpc']*MIMICA_prof['Qc']+MIMICA_prof['Qr']*MIMICA_prof['Vpr']

    Statistical_output['prf'].attrs['standard_name']= Output_Mean['prf']['std_name'] 
    Statistical_output['prf'].attrs['units']= Output_Mean['prf']['units'] 


    # Relative humidity, and relative humidity over ice:
    DEPHY_name=['hur','huri']
    MIMICA_name=['RH','RHi']    
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        Statistical_output[dephy]=1e-2*MIMICA_prof[mimica]
        Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name']
        Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units'] 

    # Mixing fraction from MIMICA are directly mass mixing ratio in DEPHY
    DEPHY_name=[ 'qv','qlc','qlr','qic','qis','qig']
    MIMICA_name=['Qv','Qc', 'Qr','Qi', 'Qs', 'Qg']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mimica in MIMICA_prof:
            Statistical_output[dephy]=MIMICA_prof[mimica]/RHO_dry
            Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name']
            Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units'] 


    # Hydrometeors number density: 
    DEPHY_name=['nlc','nlr','nic','nis','nig']
    MIMICA_name=['Nc','Nr','Ni','Ns','Ng']    
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mimica in MIMICA_prof:
            Statistical_output[dephy]=MIMICA_prof[mimica]/RHO_dry
            Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name'] 
            Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units'] 
      
    # Hydrometeors number density in cloud:
    DEPHY_name=['nlcic','niic']
    MIMICA_name=['Nc','Ni']    
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mimica in MIMICA_prof_cl:
            #Statistical_output[dephy]=MIMICA_prof_cl[mimica]/MIMICA_prof_cl['RHO']
            RHO_dry_cl=MIMICA_prof_cl['RHO']/(1+MIMICA_prof_cl['Qt'])
            Statistical_output[dephy]=save_division(MIMICA_prof_cl[mimica],RHO_dry_cl)
            Statistical_output[dephy].attrs['standard_name']= Output_Mean[dephy]['std_name'] 
            Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units']       

    # Number of total aerosol in air #For later   !activated + unactivated???
    if 'N_AERO1' in MIMICA_prof:
        DEPHY_name=['na1','na2','na3']
        MIMICA_name=['N_AERO1','N_AERO2','N_AERO3']
        for dephy,mimica in zip(DEPHY_name,MIMICA_name):
            Statistical_output[dephy]=MIMICA_prof[mimica]
            Statistical_output[dephy].attrs['units']= Output_Mean[dephy]['units']

    # Renaming vertical coordinate   
    Statistical_output=Statistical_output.rename({'Z': 'zf'})
    Statistical_output.zf.attrs['standard_name']='height'
    Statistical_output.zf.attrs['comment']='altitude of mid level points above sea surface'

    # Adding the global attributes to the file
    Statistical_output=Output_add_attributes(Statistical_output,output_name,mixed_phase)

    # Changing time format
    Statistical_output=convert_time_new(Statistical_output,2020,3,12,22,'seconds')
    
    return Statistical_output

def create_2D_output(MIMICA_surf,output_name,mixed_phase):
    
    D2_output = xr.Dataset()
    
    # Friction velocity, surface fluxes
    DEPHY_name=['ustar','hfss','hfls']
    MIMICA_name=['U_star','SHF','LHF']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        D2_output[dephy]=MIMICA_surf[mimica]
        D2_output[dephy].attrs['standard_name']= Output_2D[dephy]['std_name']  
        D2_output[dephy].attrs['units']= Output_2D[dephy]['units']         
        
    # Liquid water path, ice water path
    DEPHY_name=['lwp','iwp']
    MIMICA_name=['LWP','IWP']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        D2_output[dephy]=MIMICA_surf[mimica]
        D2_output[dephy].attrs['standard_name']= Output_2D[dephy]['std_name']     
        D2_output[dephy].attrs['units']= Output_2D[dephy]['units']       
        
    # Precipitation flux. Change from precipitation rate (MIMICA) to precipitation flux
    #water_factor= 0.997/86400
    #ice_factor= 0.917/86400

    water_factor= 1./86400
    ice_factor= 1./86400
    
    if mixed_phase:
        D2_output['pr']=ice_factor*MIMICA_surf['R_rate']
    else:
        D2_output['pr']=water_factor*MIMICA_surf['R_rate']        
            
    D2_output['pr'].attrs['standard_name']= Output_2D['pr']['std_name']
    D2_output['pr'].attrs['units']= Output_2D['pr']['units']     
    
    # Atmosphere optical thickness
    OPT= MIMICA_surf['OPTC']+MIMICA_surf['OPTI']+MIMICA_surf['OPTR']
    D2_output['opt'] = OPT  
    D2_output['opt'].attrs['standard_name']= Output_2D['opt']['std_name'] 
    D2_output['opt'].attrs['units']= Output_2D['opt']['units'] 
    
    #Pseudo-albedo
    D2_output['alb'] = OPT/(OPT+13.)
    D2_output['alb'].attrs['standard_name']= Output_2D['alb']['std_name']  
    D2_output['alb'].attrs['units']= Output_2D['alb']['units']      
    
    # Outgoing longwave radiation in atmospheric window  
    D2_output['olr11'] = MIMICA_surf['OLR_11mum']
    D2_output['olr11'].attrs['starndard_name']= Output_2D['olr11']['std_name']
    D2_output['olr11'].attrs['units']= Output_2D['olr11']['units']  
    D2_output['olr11'].attrs['description']= 'Outgoing longwave flux in window containing 10.8 micrometers (wavelength range 10-12.5 micrometers)'
    
    # Surface winds u,v
    D2_output['us'] = MIMICA_surf['U_zlevel_1']
    D2_output['us'].attrs['starndard_name']= Output_2D['us']['std_name']   
    D2_output['us'].attrs['units']= Output_2D['us']['units']
    
    D2_output['vs'] = MIMICA_surf['V_zlevel_1']
    D2_output['vs'].attrs['starndard_name']= Output_2D['vs']['std_name']       
    D2_output['vs'].attrs['units']= Output_2D['vs']['units']
    
    # Changing time format
    D2_output=convert_time_new(D2_output,2020,3,12,22,'seconds')
    
    # Adding the global attributes to the file
    D2_output=Output_add_attributes(D2_output,output_name,mixed_phase)
      
          
    return D2_output


def create_3D_output(MIMICA_3D,output_name,mixed_phase):
    
    D3_output = xr.Dataset()

    # U,V,temperature,potential temperature,water vapor mixing fraction(ratio)
    DEPHY_name=[ 'pa','ua','va','wa','ta']
    MIMICA_name=['P_tot','U', 'V','W', 'T']
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        D3_output[dephy]=MIMICA_3D[mimica]
        D3_output[dephy].attrs['standard_name']= Output_3D[dephy]['std_name']  
        D3_output[dephy].attrs['units']= Output_3D[dephy]['units'] 

    RHO_dry= MIMICA_3D['RHO']/(1+MIMICA_3D['Qt'])
    D3_output['rhoa']=RHO_dry
    D3_output['rhoa'].attrs['standard_name']= Output_3D['rhoa']['std_name']  
    D3_output['rhoa'].attrs['units']= Output_3D['rhoa']['units']
    
    # Mixing fraction from MIMICA are directly mass mixing ratio in DEPHY
    DEPHY_name=[ 'qv','qlc','qlr','qic','qis','qig']
    MIMICA_name=['Qv','Qc', 'Qr','Qi', 'Qs', 'Qg' ]
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mimica in MIMICA_3D:
            D3_output[dephy]=MIMICA_3D[mimica]/RHO_dry
            D3_output[dephy].attrs['standard_name']= Output_3D[dephy]['std_name']
            D3_output[dephy].attrs['units']= Output_3D[dephy]['units'] 
                        
    # Hydrometeors number density:   #Question here
    DEPHY_name=['nlc','nlr','nic','nis','nig']
    MIMICA_name=['Nc','Nr','Ni','Ns','Ng']    
    for dephy,mimica in zip(DEPHY_name,MIMICA_name):
        if mimica in MIMICA_3D:
            D3_output[dephy]=MIMICA_3D[mimica]/RHO_dry
            D3_output[dephy].attrs['standard_name']= Output_3D[dephy]['std_name']        
            D3_output[dephy].attrs['units']= Output_3D[dephy]['units'] 
                
    #Check these variables. Are they mass weighted?       
    #DEPHY_name=[ 'vmlc','vmlr','vmic','vmis','vmig']
    #MIMICA_name=[ 'Vpc','Vpr','Vpi','Vps','Vpg']
    #for dephy,mimica in zip(DEPHY_name,MIMICA_name):
    #    if mimica in MIMICA_3D:
    #        D3_output[dephy]=MIMICA_3D[mimica]
    #        D3_output[dephy].attrs['standard_name']= Output_3D[dephy]['std_name']
    #        D3_output[dephy].attrs['units']= Output_3D[dephy]['units']      
    
            
    # Changing time format
    D3_output=convert_time_new(D3_output,2020,3,12,22,'seconds')
    
    # Adding the global attributes to the file
    D3_output=Output_add_attributes(D3_output,output_name,mixed_phase)
            

    
    return D3_output
    
    
    
def save_division(numerator,denominator):

    # Replace zeros in the denominator with NaN to avoid division by zero
    safe_denominator = denominator.where(denominator != 0, np.nan)

    # Perform the division
    resulting_variable = numerator / safe_denominator

    # Replace NaN values in the result with zero
    resulting_variable = resulting_variable.fillna(0)

    # Assign the resulting variable back to the dataset if needed
    #COMBLE_FixN_256_new_prof_cl['resulting_variable'] = resulting_variable
    
    return resulting_variable   


def convert_time_from_minutes_to_seconds_new(file, file_time_modified):

    # Remove existing file if it exists
    if os.path.exists(file_time_modified):
        os.remove(file_time_modified)

    ds = xr.open_dataset(file)    

    # Ensure all units attributes are strings to avoid xarray CF decoding errors later
    for var in ds.variables:
        if 'units' in ds[var].attrs and not isinstance(ds[var].attrs['units'], str):
            ds[var].attrs['units'] = str(ds[var].attrs['units'])

    # Force correct time encoding in seconds
    ds.time.encoding["units"] = "seconds since 2020-03-12 22:00:00"
    ds.time.encoding["calendar"] = "proleptic_gregorian"

    # Save modified dataset
    ds.to_netcdf(file_time_modified)

    return ds

Simulation with ice part 1: FixN#

#Reading MIMICA native outputs
COMBLE_FixN_256_prof=   xr.open_dataset(native_outputs_path+'COMBLE_FixN_256x256/OUTPUT/profiles_tot.nc')
COMBLE_FixN_256_prof_cl=xr.open_dataset(native_outputs_path+'COMBLE_FixN_256x256/OUTPUT/profiles_cl.nc')
COMBLE_FixN_256_surf=   xr.open_dataset(native_outputs_path+'COMBLE_FixN_256x256/OUTPUT/slice_surf.nc')
COMBLE_FixN_256_3D=     xr.open_dataset(native_outputs_path+'COMBLE_FixN_256x256/OUTPUT/COMBLE_MPI.nc')
COMBLE_FixN_256_TS=     xr.open_dataset(native_outputs_path+'COMBLE_FixN_256x256/OUTPUT/T_S.nc')

Stat_output_FixN_256=create_stat_output(COMBLE_FixN_256_prof,COMBLE_FixN_256_prof_cl,COMBLE_FixN_256_surf,
                                         COMBLE_FixN_256_TS,'MIMICA_Lx25_dx100_FixN',True)

OUTPUT_2D_FixN_256=create_2D_output(COMBLE_FixN_256_surf,'MIMICA_Lx25_dx100_FixN_2D',True)

#OUTPUT_3D_FixN_256=create_3D_output(COMBLE_FixN_256_3D,'MIMICA_Lx25_dx100_FixN_3D_alt',True)

Create netcdf files for FixN#

Stat_output_FixN_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN'+'.nc',mode='w',engine="netcdf4")

OUTPUT_2D_FixN_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN_2D'+'.nc',mode='w',engine="netcdf4")

#OUTPUT_3D_FixN_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN_3D'+'.nc',mode='w',engine="netcdf4")

Converting time from minutes to seconds#

#Doing this only because Tim Juliano could not read properly the netcdf files after the conversion
file_names=[ 'MIMICA_Lx25_dx100_FixN.nc','MIMICA_Lx25_dx100_FixN_2D.nc']#,  'MIMICA_Lx25_dx100_FixN_3D.nc', 

for ii_file in file_names:
    print("Converting time in file ",ii_file, "saving in",converted_outputs_path_time_mod+ii_file)
    convert_time_from_minutes_to_seconds_new(converted_outputs_path+ii_file,converted_outputs_path_time_mod+ii_file)  

Converting time in file  MIMICA_Lx25_dx100_FixN.nc saving in /proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED_time_mod/MIMICA_Lx25_dx100_FixN.nc
Converting time in file  MIMICA_Lx25_dx100_FixN_2D.nc saving in /proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED_time_mod/MIMICA_Lx25_dx100_FixN_2D.nc

Simulation without ice part 1#

#Reading MIMICA native outputs
COMBLE_FixN_noice_256_prof=   xr.open_dataset(native_outputs_path+'COMBLE_FixN_noice_256x256/OUTPUT/profiles_tot.nc')
COMBLE_FixN_noice_256_prof_cl=xr.open_dataset(native_outputs_path+'COMBLE_FixN_noice_256x256/OUTPUT/profiles_cl.nc')
COMBLE_FixN_noice_256_surf=   xr.open_dataset(native_outputs_path+'COMBLE_FixN_noice_256x256/OUTPUT/slice_surf.nc')
COMBLE_FixN_noice_256_3D=     xr.open_dataset(native_outputs_path+'COMBLE_FixN_noice_256x256/OUTPUT/COMBLE_MPI.nc')
COMBLE_FixN_noice_256_TS=     xr.open_dataset(native_outputs_path+'COMBLE_FixN_noice_256x256/OUTPUT/T_S.nc')

Stat_output_FixN_noice_256=create_stat_output(COMBLE_FixN_noice_256_prof,COMBLE_FixN_noice_256_prof_cl,COMBLE_FixN_noice_256_surf,
                                         COMBLE_FixN_noice_256_TS,'MIMICA_Lx25_dx100_FixN_noice',False)

OUTPUT_2D_FixN_noice_256=create_2D_output(COMBLE_FixN_noice_256_surf,'MIMICA_Lx25_dx100_FixN_noice_2D',False)

#OUTPUT_3D_FixN_noice_256=create_3D_output(COMBLE_FixN_noice_256_3D,'MIMICA_Lx25_dx100_FixN_noice_3D',False)

Create Netcdf files for FixN noice#

Stat_output_FixN_noice_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN_noice'+'.nc',mode='w',engine="netcdf4")

OUTPUT_2D_FixN_noice_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN_noice_2D'+'.nc',mode='w',engine="netcdf4")

#OUTPUT_3D_FixN_noice_256.to_netcdf(path=converted_outputs_path+'MIMICA_Lx25_dx100_FixN_noice_3D'+'.nc',mode='w',engine="netcdf4")

Converting time from minutes to seconds#

#Doing this only because Tim Juliano could not read properly the netcdf files after the conversion
file_names= ['MIMICA_Lx25_dx100_FixN_noice.nc','MIMICA_Lx25_dx100_FixN_noice_2D.nc'] #  'MIMICA_Lx25_dx100_FixN_noice_3D.nc' ]

for ii_file in file_names:
    print("Converting time in file ",ii_file, "saving in",converted_outputs_path_time_mod+ii_file)
    convert_time_from_minutes_to_seconds_new(converted_outputs_path+ii_file,converted_outputs_path_time_mod+ii_file) 

Converting time in file  MIMICA_Lx25_dx100_FixN_noice.nc saving in /proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED_time_mod/MIMICA_Lx25_dx100_FixN_noice.nc
Converting time in file  MIMICA_Lx25_dx100_FixN_noice_2D.nc saving in /proj/bolinc/users/x_abaro/MIMICA_radiation/COMBLE_CONVERTED_time_mod/MIMICA_Lx25_dx100_FixN_noice_2D.nc
Contents