{ "cells": [ { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "import netCDF4 as nc\n", "import numpy as np\n", "import pandas as pd\n", "import os\n", "import matplotlib.pyplot as plt\n", "\n", "def utowgrid(ds,var): \n", " #-- ds has to be 3d or pr PALM output (netCDF) (1st dimension: time, 2nd dimension: z)\n", " #-- var: variable that's output on the u grid and should be converted to w grid (by averaging over the two closest u-levels)\n", " #-- averages output (u-grid) levels to match w-grid levels\n", " vals = np.squeeze(ds[var])\n", " # zu_new = np.zeros(zu.shape)\n", " vals_new = np.zeros((vals.shape))\n", " # zu_new[0] = zu[0]\n", " vals_new[:,0] = vals[:,0]\n", " # zu_new[-1] = np.nan\n", " vals_new[:,-1] = np.nan\n", " for i in np.arange(1,vals.shape[1]-1,1):\n", " # zu_new[i] = (zu[i]+zu[i+1])/2\n", " vals_new[:,i] = (vals[:,i]+vals[:,i+1])/2\n", " return vals_new\n", "\n", "def utowgrid_1d(ds,var): #-- ds has to be PALM output (netCDF) dimension (NO variable with time dimension!) / var: variable that's output on the u grid and should be converted to w grid (by averaging over the two closest u-levels)\n", " #-- averages output on u-grid levels to match w-grid levels\n", " vals = np.squeeze(ds[var])\n", " # zu_new = np.zeros(zu.shape)\n", " vals_new = np.zeros((vals.shape))\n", " # zu_new[0] = zu[0]\n", " vals_new[0] = vals[0]\n", " # zu_new[-1] = np.nan\n", " vals_new[-1] = np.nan\n", " for i in np.arange(1,vals.shape[0]-1,1):\n", " # zu_new[i] = (zu[i]+zu[i+1])/2\n", " vals_new[i] = (vals[i]+vals[i+1])/2\n", " return vals_new\n", "\n", "# def pr23d_time(t_3d, t_pr, pr_var): #-- average 15 min output (only one height) to 30 min output (3d)\n", "# res = np.zeros((t_3d.shape[0]))\n", "# for i in np.arange(t_3d.shape[0]):\n", "# indx = np.where((t_pr >= t_3d[i]) & (t_pr < t_3d[i]+1799))[0]\n", "# res[i] = np.average(pr_var[indx],axis=0)\n", "# return res\n", "\n", "def timestep3d(var_in, time_in, time_3d):\n", " if len(var_in.shape) == 1:\n", " out = np.zeros((time_3d.shape))\n", " elif len(var_in.shape) == 2:\n", " out = np.zeros((time_3d.shape[0],var_in.shape[1]))\n", " elif len(var_in.shape) == 3:\n", " out = np.zeros((time_3d.shape[0],var_in.shape[1],var_in.shape[2]))\n", " elif len(var_in.shape) == 4:\n", " out = np.zeros((time_3d.shape[0],var_in.shape[1],var_in.shape[2],var_in.shape[3]))\n", " for i in np.arange(len(time_3d)):\n", " if i == 0:\n", " out[i] = var_in[np.where((time_in > 0) & (time_in <= time_3d[i]))].mean(axis=0)\n", " else:\n", " out[i] = var_in[np.where((time_in > time_3d[i-1]) & (time_in <= time_3d[i]))].mean(axis=0)\n", " return out\n", "\n", "def find_nearest(array, value):\n", " array = np.asarray(array)\n", " idx = (np.abs(array - value)).argmin()\n", " return array[idx]\n", "\n", "def CanopyTopVals(arr1,arr2,lad): # arr1: surface values / arr2: atmospheric values\n", " # (1) copy surface value in every location (Y,X)\n", " out = arr1.copy()\n", " # (2) where lad[1] > 0, there are trees at least up to the first grid cell --> overwrite surface temperature with air temperature from that height at those locations\n", " # --> loop through all lad z levels\n", " for zl in np.arange(1,lad.shape[0]):\n", " out[:,np.where(lad[zl,:,:]>0)[0],np.where(lad[zl,:,:]>0)[1]] = arr2[:,zl][:,np.where(lad[zl,:,:]>0)[0],np.where(lad[zl,:,:]>0)[1]] \n", " return out\n", "\n", "#def c1pr_to_c23dav(c1_var, t_pr, t_3d, z_pr, z_3d): # input should be a profile variable from the c1 output, already interpolated to zw values (if it is a zu variable)\n", "# #-- replace first time step value by nan\n", "# c1_var[0] = np.nan \n", "# #-- average over 30 minutes (pr output is in 15 min intervals)\n", "# c1_var_hh = timestep3d(c1_var,t_pr,t_3d) \n", "# #-- create dataframe with z levels of child2 output (fine grid), fill with values where available\n", "# c1_var_fg = np.zeros((len(t_3d),len(z_3d)))\n", "# for i in np.arange(len(z_3d)):\n", "# if z_3d[i] in z_pr:\n", "# c1_var_fg[:,i] =c1_var_hh[:,np.where(z_pr == z_3d[i])[0][0]]\n", "# else:\n", "# c1_var_fg[:,i] = np.nan\n", "# #-- create pandas dataframe and interpolate\n", "# df = pd.DataFrame(c1_var_fg.T,columns = t_3d,index=z_3d)\n", "# df_int = df.interpolate(method='linear',axis=0)\n", "# out = df_int.to_numpy().T\n", "# return out" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [], "source": [ "#-- settings / simulation information / paths \n", "IOP = 'ches_IOP2'\n", "EM = 'ensemble.member.1'\n", "path = '/air/lwanner/'+IOP+'/'+EM+'/'\n", "path_static_file = '/bog/incoming/CHEESEHEAD/palm/realistic_runs/'+IOP+'/INPUT/'# copy here or find out where this is located on co2 'C:/Luise/PALM/CHEESEHEAD_simulations/'\n", "pathOUT = '/air/lwanner/palm_processed/'+IOP+'/'+EM+'/' #'C:/Luise/PALM/CHEESEHEAD_simulations/ches_'+IOP+'.'+sim+'/processed/'\n", "pathPLOT = '/air/lwanner/palm_plots/'+IOP+'/'+EM+'/'\n", "if not os.path.exists(pathOUT):\n", " os.makedirs(pathOUT)\n", "if not os.path.exists(pathPLOT):\n", " os.makedirs(pathPLOT)\n", "twr_sinfo = {'nw1': ['02', 32, 'pine'],\n", " 'nw2': ['03', 12, 'aspen'],\n", " 'nw3': ['04', 3, 'tussoc'],\n", " 'nw4': ['05', 32, 'lake'],\n", " 'ne1': ['06', 32, 'pine'],\n", " 'ne2': ['07', 32, 'pine'],\n", " 'ne3': ['08', 32, 'hardwood'],\n", " 'ne4': ['09', 32, 'maple'],\n", " 'sw1': ['10', 32, 'aspen'],\n", " 'sw2': ['11', 25, 'aspen'],\n", " 'sw3': ['12', 32, 'hardwood'],\n", " 'sw4': ['13', 32, 'hardwood'],\n", " 'se2': ['14', 32, 'hardwood'],\n", " 'se3': ['15', 32, 'aspen'],\n", " 'se4': ['16', 3, 'tussoc'],\n", " 'se5': ['17', 12, 'aspen'],\n", " 'se6': ['18', 32, 'pine']}\n", "\n", "#-- mask id in palm output\n", "msk_id = {'nw1':'02', 'nw2':'03', 'nw3':'04', 'nw4':'05',\n", " 'ne1':'06', 'ne2':'07', 'ne3':'08', 'ne4':'09',\n", " 'sw1':'10', 'sw2':'11', 'sw3':'12', 'sw4':'13',\n", " 'se2':'14', 'se3':'15', 'se4':'16', 'se5':'17', 'se6':'18'}\n", "#-- measurement height (in the field)\n", "h_m = {'nw1':32, 'nw2':12, 'nw3': 3, 'nw4':32,\n", " 'ne1':32, 'ne2':32, 'ne3':32, 'ne4':32,\n", " 'sw1':32, 'sw2':25, 'sw3':32, 'sw4':32,\n", " 'se2':32, 'se3':32, 'se4': 3, 'se5':12, 'se6':32}\n", "#-- canopy height (in the field)\n", "h_c = {'nw1':25.0, 'nw2': 3.0, 'nw3': .3, 'nw4':20.1,\n", " 'ne1':33.2, 'ne2':19.2, 'ne3':18.3, 'ne4':18.3,\n", " 'sw1':24.4, 'sw2':19.2, 'sw3':15.0, 'sw4':25.9,\n", " 'se2':24.4, 'se3':14.3, 'se4': .3, 'se5': 3.1, 'se6':21.6}\n", "\n", "\n", "origin = {'y_c2': 20070, #-- values are from the ches_IOP2_p3d file (see nesting section)\n", " 'x_c2': 18180, #-- only needs to be used when diffrent domains are combined\n", " 'y_c1': 10890,\n", " 'x_c1': 10710,\n", " 'y_p' : 0,\n", " 'x_p' : 0}\n", "\n", "cp = 1004.67 # Foken book p. 238\n", "#lv = 2501000 # Foken book p. 238" ] }, { "cell_type": "code", "execution_count": 113, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "output saved to /air/lwanner/palm_processed/ches_IOP2/ensemble.member.1/\n" ] } ], "source": [ "#-- first, determine canopy heights at each tower site\n", "static_c2 = nc.Dataset(path_static_file+IOP+'_static_N03') #-- static file of child 2 (contains fine grid lad information)\n", "twrs = ['nw1', 'nw2','ne1', 'ne2', 'ne3','ne4', 'sw1','sw2','sw3','sw4','se2', 'se3','se5']\n", "\n", "h_c_m = dict() #-- canopy height in model\n", "h_m_m = dict() #-- measurement height in model\n", "for t in np.arange(len(twrs)):\n", " twr = twrs[t]\n", " data_3d = nc.Dataset(path+'DATA_3D_AV_NETCDF_N03_M'+msk_id[twr]+'_'+IOP+'_'+EM)\n", " \n", " x_3d = np.squeeze(data_3d['x'])\n", " y_3d = np.squeeze(data_3d['y'])\n", " z_lad = np.squeeze(static_c2['zlad'])\n", " x_lad = np.squeeze(static_c2['x'])\n", " y_lad = np.squeeze(static_c2['y'])\n", " z_3d = np.squeeze(data_3d['zw_3d'])\n", " zu_3d = np.squeeze(data_3d['zu_3d'])\n", " #-- get lad at the tower site\n", " lad = np.average(np.squeeze(static_c2['lad'])[:, np.where((y_lad >= y_3d.min()) & (y_lad <= y_3d.max()))[0],:][:,:,np.where((x_lad >= x_3d.min()) & (x_lad <= x_3d.max()))[0]],axis=(1,2))\n", " if lad[1:].max() <= 0:\n", " h_c_m[twr] = 0 # height of canopy in model\n", " h_m_m[twr] = z_3d[5] # height of measurement in model\n", " else:\n", " h_c_m[twr] = z_lad[np.where(lad[1:] == 0)[0].min()-1]\n", " h_m_m[twr] = z_3d[int(np.where(zu_3d == z_lad[np.where(lad[1:] == 0)[0].min()])[0]+4)]\n", "#-- save to csv\n", "pd.DataFrame([h_c_m]).to_csv(pathOUT+'PALM_canopy_heights_virtual_sites_'+IOP+'.csv',index=False)\n", "pd.DataFrame([h_m_m]).to_csv(pathOUT+'PALM_measurement_heights_virtual_sites_'+IOP+'.csv',index=False)\n", "print('output saved to '+ pathOUT)" ] }, { "cell_type": "code", "execution_count": 114, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tcanopy created\n" ] } ], "source": [ "#-- get surface temperature map from child 1 3d output and child 1 static input\n", "data_3dc1 = nc.Dataset(path+'DATA_3D_AV_NETCDF_N02'+'_slice_'+IOP+'_'+EM)\n", "data_xyc1 = nc.Dataset(path+'DATA_2D_XY_AV_NETCDF_N02'+'slice_'+IOP+'_'+EM)\n", "static_c1 = nc.Dataset(path_static_file+IOP+'_static_N02') \n", "lad = np.squeeze(static_c1['lad'])\n", "Tsurf = np.squeeze(data_xyc1['tsurf*_xy'])\n", "Tair = np.squeeze(data_3dc1['theta'])\n", "\n", "#-- get canopy temperature without looping through each y,x location (takes long)\n", "Tcanopy = CanopyTopVals(Tsurf,Tair,lad)\n", "print('Tcanopy created')" ] }, { "cell_type": "code", "execution_count": 115, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "saved surface_temperature_ches_IOP2_ensemble.member.1_nw1.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_nw2.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_ne1.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_ne2.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_ne3.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_ne4.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_sw1.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_sw2.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_sw3.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_sw4.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_se2.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_se3.nc\n", "saved surface_temperature_ches_IOP2_ensemble.member.1_se5.nc\n", "output stored in /air/lwanner/palm_processed/ches_IOP2/ensemble.member.1/\n" ] } ], "source": [ "#-- extract patches for towers **TALL TOWER and smaller towers are MISSING**\n", "twrs = ['nw1', 'nw2','ne1', 'ne2', 'ne3','ne4', 'sw1','sw2','sw3','sw4','se2', 'se3','se5']\n", "# (1) get list of tower locations for x and y\n", "x_c1 = np.squeeze(data_xyc1['x'])+origin['x_c1']\n", "y_c1 = np.squeeze(data_xyc1['y'])+origin['y_c1']\n", "\n", "ylocs = []\n", "xlocs = []\n", "for i in np.arange(len(twrs)):\n", " twr = twrs[i]\n", " #-- x and y location of the tower\n", " data_3d = nc.Dataset(path+'DATA_3D_AV_NETCDF_N03_M'+msk_id[twr]+'_'+IOP+'_'+EM)\n", " yloc = np.median(np.squeeze(data_3d['y']))+origin['y_c2']\n", " xloc = np.median(np.squeeze(data_3d['x']))+origin['x_c2']\n", " \n", " ylocs.append(yloc)\n", " xlocs.append(xloc)\n", " idx_y = np.where(y_c1 == find_nearest(y_c1,yloc))[0][0]\n", " idx_x = np.where(x_c1 == find_nearest(x_c1,xloc))[0][0]\n", " \n", " dxy = y_c1[1]-y_c1[0]\n", " \n", " l = 6300 # length of the area #-- too small, but matches the size of the LES that was used to determine the scaling functions\n", " n = int(l/2/dxy)\n", " \n", " ymin = idx_y-n\n", " ymax = idx_y+n\n", " xmin = idx_x-n\n", " xmax = idx_x+n\n", " Tc_twr = Tcanopy[:,ymin:ymax, xmin:xmax].copy()\n", " \n", " #-- save as netcdf file\n", " ncout = nc.Dataset(pathOUT + 'surface_temperature_'+IOP+'_'+EM+'_'+twr+'.nc','w','NETCDF3')\n", " ncout.createDimension('time',Tc_twr.shape[0]);\n", " ncout.createDimension('y',Tc_twr.shape[1]);\n", " ncout.createDimension('x',Tc_twr.shape[2]);\n", " time = ncout.createVariable('time','float32',('time'));\n", " time.setncattr('units','s'); time[:] = np.squeeze(data_3d['time']);\n", " yvar = ncout.createVariable('y' ,'float32',('y' ));\n", " yvar.setncattr('units','m'); yvar[:] = y_c1[ymin:ymax];\n", " xvar = ncout.createVariable('x','float32',('x'));\n", " xvar.setncattr('units','m'); xvar[:] = x_c1[xmin:xmax];\n", " Ts = ncout.createVariable('T_surf','float32',('time','y','x'));\n", " Ts.setncattr('units','K'); Ts[:] = Tc_twr; \n", " ncout.close()\n", " print('saved '+ 'surface_temperature_'+IOP+'_'+EM+'_'+twr+'.nc')\n", "print('output stored in '+ pathOUT)" ] }, { "cell_type": "code", "execution_count": 116, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "saved Output_02_nw1_ches_IOP2_ensemble.member.1.csv\n", "saved Output_03_nw2_ches_IOP2_ensemble.member.1.csv\n", "saved Output_06_ne1_ches_IOP2_ensemble.member.1.csv\n", "saved Output_07_ne2_ches_IOP2_ensemble.member.1.csv\n", "saved Output_08_ne3_ches_IOP2_ensemble.member.1.csv\n", "saved Output_09_ne4_ches_IOP2_ensemble.member.1.csv\n", "saved Output_10_sw1_ches_IOP2_ensemble.member.1.csv\n", "saved Output_11_sw2_ches_IOP2_ensemble.member.1.csv\n", "saved Output_12_sw3_ches_IOP2_ensemble.member.1.csv\n", "saved Output_13_sw4_ches_IOP2_ensemble.member.1.csv\n", "saved Output_14_se2_ches_IOP2_ensemble.member.1.csv\n", "saved Output_15_se3_ches_IOP2_ensemble.member.1.csv\n", "saved Output_17_se5_ches_IOP2_ensemble.member.1.csv\n", "output stored in /air/lwanner/palm_processed/ches_IOP2/ensemble.member.1/\n" ] } ], "source": [ "# -- continue with calculating fluxes, use the canopy height information to choose a virtual measurement height\n", "data_pr = nc.Dataset(path+'DATA_1D_PR_NETCDF_N03'+'slice_'+IOP+'_'+EM)\n", "data_pr_c1 = nc.Dataset(path+'DATA_1D_PR_NETCDF_N02slice_'+IOP+'_'+EM)\n", "data_xy = nc.Dataset(path+'DATA_2D_XY_AV_NETCDF_N03'+'slice_'+IOP+'_'+EM)\n", "data_ts = nc.Dataset(path+'DATA_1D_TS_NETCDF_N02'+'slice_'+IOP+'_'+EM)\n", "\n", "#twrs = ['nw1', 'ne2', 'sw1', 'se3'] # picked out 4 different towers (forested, measurement height in the field: 32 m)\n", "twrs = ['nw1', 'nw2','ne1', 'ne2', 'ne3','ne4', 'sw1','sw2','sw3','sw4','se2', 'se3','se5']\n", "for i in np.arange(len(twrs)):\n", " twr = twrs[i]\n", " # data_ms = nc.Dataset(path+'DATA_MASK_NETCDF_N03_M'+twr_info[twr][0]+'slice_subsetLW')\n", " data_3d = nc.Dataset(path+'DATA_3D_AV_NETCDF_N03_M'+msk_id[twr]+'_'+IOP+'_'+EM)\n", "\n", " #-- coordinates\n", " z_3d = np.squeeze(data_3d['zw_3d'])\n", " zu_3d = np.squeeze(data_3d['zu_3d'])\n", " z_pr = np.squeeze(data_pr['zw\"theta\"'])\n", " zu_pr = np.squeeze(data_pr['ztheta'])\n", " x_3d = np.squeeze(data_3d['x'])\n", " y_3d = np.squeeze(data_3d['y'])\n", " x_xy = np.squeeze(data_xy['x'])\n", " y_xy = np.squeeze(data_xy['y'])\n", " \n", " xyIndx = [np.where(y_xy == y_3d[0])[0][0],np.where(y_xy == y_3d[-1])[0][0]+1,np.where(x_xy == x_3d[0])[0][0],np.where(x_xy == x_3d[-1])[0][0]+1]\n", "\n", " z_idx_c = int(np.where(zu_3d == h_c_m[twr])[0]) #-- \"surface\" height (at u grid) --> canopy top\n", " z_idx_m = int(np.where(z_3d == h_m_m[twr])[0]) #-- measurement height (at w grid) --> 10 m above canopy top\n", " time_3d = np.squeeze(data_3d['time'])\n", " time_pr = np.squeeze(data_pr['time'])\n", " \n", " #-- get variables to calculate temporal and spatial covariance\n", " rho_c1 = utowgrid(data_pr_c1,'rho')[1:]\n", " t_pr_c1 = np.squeeze(data_pr_c1['time'])[1:]\n", " z_rho_wgrid = utowgrid_1d(data_pr_c1,'zrho')\n", " rho_c1_3dt = timestep3d(rho_c1, t_pr_c1, time_3d)\n", " rho = np.zeros((len(time_3d),len(z_3d))) #-- interpolate to fine grid z levels\n", " for j in np.arange(len(time_3d)):\n", " rho[j,:] = np.interp(z_3d,z_rho_wgrid,rho_c1[j,:])\n", " # rho = 1.225 # Foken book p. 238\n", " corrH = rho*cp #--> use actual rho!!! \n", " w = np.squeeze(data_3d['w']).mean(axis=(2,3)) #-- horizontal average over 3x3 grid points\n", " t = utowgrid(data_3d, 'theta').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- output on u-grid --> converted\n", " q = utowgrid(data_3d, 'q').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- output on u-grid --> converted\n", " u = utowgrid(data_3d, 'u').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- output on u-grid --> converted\n", " v = utowgrid(data_3d, 'v').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- output on u-grid --> converted\n", " wt = np.squeeze(data_3d['wtheta']).mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- in K m s-1 --> convert to W m-2 later\n", " wq = np.squeeze(data_3d['wq']).mean(axis=(2,3)) #-- horizontal average over 3x3 grid points #-- in kg kg-1 m s-1 --> convert to W m-2 later\n", " wu = utowgrid(data_3d,'uw').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points \n", " wv = utowgrid(data_3d,'vw').mean(axis=(2,3)) #-- horizontal average over 3x3 grid points \n", " #-- calculate lv\n", " lv = (2.501 - 0.00237*(t-273.15)) *1000000 # uses theta instead of T! (Foken 2017, p. 331)\n", " corrLE = rho*lv\n", " #-- replace lowest grid level with surface fluxes from xy output\n", " # #-- 30-min averaged surface fluxes \n", " hf_s = np.squeeze(np.squeeze(data_xy['shf*_xy' ])[:,xyIndx[0]:xyIndx[1],xyIndx[2]:xyIndx[3]]).mean(axis=(1,2)) #-- surface sensible heat flux of last 30 minutes at tower location\n", " lf_s = np.squeeze(np.squeeze(data_xy['qsws*_xy'])[:,xyIndx[0]:xyIndx[1],xyIndx[2]:xyIndx[3]]).mean(axis=(1,2)) #-- surface latent heat flux of last 30 minutes\n", " \n", " wt[:,0] = hf_s.copy()\n", " wq[:,0] = lf_s.copy()\n", " # TAKE FLUX AT CANOPY TOP INSTEAD!!!\n", " hf_s_sgs = timestep3d(np.squeeze(data_pr['w\"theta\"'])[:,z_idx_c],time_pr,time_3d) #pr23d_time(time_3d,time_pr,np.squeeze(data_pr['w\"theta\"'])[:,z_idx_c]) #-- horizontally homogeneous sgs contribution\n", " lf_s_sgs = timestep3d(np.squeeze(data_pr['w\"q\"'])[:,z_idx_c],time_pr,time_3d) #-- horizontally homogeneous sgs contribution\n", " \n", " #-- compute resolved part of \"surface\" flux\n", " #hf_s_res = np.zeros((wt.shape[0]))\n", " #for tstp in np.arange(wt.shape[0]): #change #???\n", " # hf_s_res[tstp] = (wt - w*t)[tstp,z_idx_c] *rho[tstp,z_idx_c]*cp\n", " hf_s_res = (wt - w*t)[:,z_idx_c] * corrH[:,z_idx_c]\n", " lf_s_res = (wq - w*q)[:,z_idx_c] * corrLE[:,z_idx_c]\n", " #-- total \"surface flux\n", " hf_s = hf_s_res + hf_s_sgs\n", " lf_s = lf_s_res + lf_s_sgs\n", " \n", " #-- subgridscale contribution of turbulent fluxes ONLY at tower height\n", " hf_sgs = timestep3d(np.squeeze(data_pr['w\"theta\"'])[:,:len(z_3d)],time_pr,time_3d) #-- horizontally homogeneous sgs contribution\n", " lf_sgs = timestep3d(np.squeeze(data_pr['w\"q\"'])[:,:len(z_3d)],time_pr,time_3d) #-- horizontally homogeneous sgs contribution\n", " wu_sgs = timestep3d(np.squeeze(data_pr['w\"u\"'])[:,:len(z_3d)],time_pr,time_3d) #-- horizontally homogeneous sgs contribution\n", " wv_sgs = timestep3d(np.squeeze(data_pr['w\"v\"'])[:,:len(z_3d)],time_pr,time_3d) #-- horizontally homogeneous sgs contribution\n", " \n", " #-- compute resolved part of turbulent fluxes ONLY at tower location/height (see PALM documentation (https://palm.muk.uni-hannover.de/trac/wiki/doc/app/runtime_parameters#data_output))\n", " #hf_res = np.zeros((wt.shape)) #change #???\n", " #for tstp in np.arange(wt.shape[0]): #change #???\n", " # hf_res[tstp] = (wt-w*t)[tstp,:] * corrH[tstp,:] #change #???\n", " hf_res = (wt - w*t) * corrH\n", " lf_res = (wq - w*q) * corrLE\n", " wu_res = (wu - w*u)\n", " wv_res = (wv - w*v)\n", " #-- calculate total turbulent fluxes (subgrid scale + resolved part)\n", " hf_t = hf_res + hf_sgs\n", " lf_t = lf_res + lf_sgs\n", " wu_t = wu_res + wu_sgs\n", " wv_t = wv_res + wv_sgs\n", " \n", " #-- plot\n", "# for tstp in np.arange(wt.shape[0]):\n", "# fig, axs = plt.subplots(1,2,figsize=(15,10))\n", "# axs[0].plot(hf_t[tstp,1:],z_3d[1:],color='red',linestyle='solid',label='hf_tot')\n", "# axs[0].plot(hf_res[tstp,1:],z_3d[1:],color='red',linestyle='dashed',label='hf_res')\n", "# axs[0].plot(hf_sgs[tstp,1:],z_3d[1:],color='red',linestyle='dotted',label='hf_sgs')\n", "# axs[0].legend()\n", "# axs[0].set_ylabel('z (m)')\n", "# axs[0].set_xlabel('H (W m$^{-2}$)')\n", "# axs[1].plot(lf_t[tstp,1:],z_3d[1:],color='blue',linestyle='solid',label='lf_tot')\n", "# axs[1].plot(lf_res[tstp,1:],z_3d[1:],color='blue',linestyle='dashed',label='lf_res')\n", "# axs[1].plot(lf_sgs[tstp,1:],z_3d[1:],color='blue',linestyle='dotted',label='lf_sgs')\n", "# fig.suptitle(str(time_3d[tstp]/3600))\n", "# axs[1].legend()\n", "# axs[1].set_xlabel('LE (W m$^{-2}$)')\n", "# plt.savefig(pathPLOT+'FluxProfiles_'+IOP+'_'+EM+'_'+twr+'_'+str(time_3d[tstp]/3600)+'.png',dpi=200,bbox_inches='tight')\n", " # print('fluxplot saved: '+str(time_3d[time_3d[tstp]/3600]))\n", "# plt.close()\n", " \n", " hf_all = np.stack([hf_sgs,hf_res,hf_t])\n", " lf_all = np.stack([lf_sgs,lf_res,lf_t])\n", " np.save(pathOUT+'H_profiles_'+IOP+'_'+EM+'_'+twr+'_new.npy',hf_all)\n", " np.save(pathOUT+'LE_profiles_'+IOP+'_'+EM+'_'+twr+'_new.npy',lf_all)\n", " \n", " #-- extract values at virtual measurement height\n", " hf_sgs_m = hf_sgs[:,z_idx_m]\n", " hf_res_m = hf_res[:,z_idx_m]\n", " hf_tot_m = hf_t[:,z_idx_m]\n", " lf_sgs_m = lf_sgs[:,z_idx_m]\n", " lf_res_m = lf_res[:,z_idx_m]\n", " lf_tot_m = lf_t[:,z_idx_m]\n", " #-- extract values at surface height (canopy top)\n", " hf_sgs_c = hf_sgs[:,z_idx_c]\n", " hf_res_c = hf_res[:,z_idx_c]\n", " hf_tot_c = hf_t[:,z_idx_c]\n", " lf_sgs_c = lf_sgs[:,z_idx_c]\n", " lf_res_c = lf_res[:,z_idx_c]\n", " lf_tot_c = lf_t[:,z_idx_c]\n", " #-- calculate dispersive fluxes\n", " hf_d = hf_tot_c - hf_tot_m\n", " lf_d = lf_tot_c - lf_tot_m\n", " #-- get some other variables needed for further processing\n", " #-- zi\n", " zi = np.squeeze(data_ts['zi_wtheta'])\n", " time_ts = np.squeeze(data_ts['time'])\n", " #-- average zi over similar time intervals as 3d av data\n", " zi_av = timestep3d(zi,time_ts,time_3d)\n", " \n", " #-- Ug\n", " Ug = np.power(np.power(u[:,:],2)+np.power(v[:,:],2),1/2)[:,z_idx_m]\n", " \n", " #-- calculate ustar and wstar\n", " #-- calculate ustar\n", " # https://glossary.ametsoc.org/wiki/Friction_velocity\n", " ust = np.power((np.power(wu_t,2)+np.power(wv_t,2)),1/4)[:,z_idx_m]\n", " #-- calculate wstar\n", " # https://glossary.ametsoc.org/wiki/Deardorff_velocity\n", " g = 9.81 # m s-2\n", " T = t[:,z_idx_m] \n", " wst = np.power(np.abs(g/T * zi_av * (hf_tot_m/corrH[:,z_idx_m])),1/3)# (g = 9.81 m s-2)\n", " #-- calculate ustar/wstar\n", " ust_wst = ust/wst\n", "\n", " #-- get surface information from xy N02!!! #??? yes for grassland/water areas, otherwise: surface temperature at tree top --> 3d N02 entire domain but only very limited z extent, containing theta only (lad from static input file N02)\n", " ds_Ts = nc.Dataset(pathOUT+'surface_temperature_'+IOP+'_'+EM+'_'+twr+'.nc')\n", " Tsurf = np.squeeze(ds_Ts['T_surf'])\n", " Tsurf_mean = Tsurf.mean(axis=(1,2))\n", " dT = np.max(Tsurf,axis=(1,2)) -np.min(Tsurf,axis=(1,2))\n", " \n", " #-- gradients of mixing ratio and potential temperature --> use child1 because it covers entire boundary layer\n", " q_c1 = timestep3d(utowgrid_1d(data_pr_c1,'q')[1:], t_pr_c1, time_3d)\n", " t_c1 = timestep3d(utowgrid_1d(data_pr_c1,'theta')[1:], t_pr_c1, time_3d)\n", " z_c1 = timestep3d(utowgrid_1d(data_pr_c1,'ztheta')[1:], t_pr_c1, time_3d)\n", " \n", " dq = np.zeros((len(time_3d)))\n", " dt = np.zeros((len(time_3d)))\n", " for j in np.arange(len(time_3d)):\n", " #-- find the z-level closest to 1/2 zi\n", " z_index = np.where(np.abs(z_c1-zi_av[j]/2) == np.nanmin(np.abs(z_c1-zi_av[j]/2)))[0][-1]\n", " if q_c1[j,z_idx_m] == 0:\n", " dq[j] = np.nan\n", " else:\n", " dq[j] = (q_c1[j,z_idx_m]-q_c1[j,z_index])/q_c1[j,z_idx_m]\n", " dt[j] = (t_c1[j,z_idx_m]-t_c1[j,z_index])/t_c1[j,z_idx_m]\n", " \n", " #-- save profiles\n", " pr_all = np.stack([t_c1,q_c1])\n", " np.save(pathOUT+'t_q_profiles_c1wgrid'+IOP+'_'+EM+'_'+twr+'_new.npy',pr_all)\n", " \n", " #-- stack all fluxes to pandas dataframe and then save\n", " #-- later: read lacunarity file and add lacunarity value\n", " \n", " #-- fluxes, ustar, wstar, Ug, T, deltaT, zi, z, (lc) \n", " OUT = pd.DataFrame(np.stack((time_3d,hf_tot_c, hf_tot_m, hf_d, hf_res_m, hf_sgs_m,lf_tot_c, lf_tot_m, lf_d, lf_res_m, lf_sgs_m,zi_av,Ug,T,Tsurf_mean,dT,ust,wst,ust_wst,dt,dq)).T,columns = ['time','Hsurf','Hturb','Hdisp','Hres','Hsgs','LEsurf','LEturb','LEdisp','LEres','LEsgs','zi','Ug','Tair','Tsurf_mean','deltaTsurf','ustar','wstar','ust_wst','delta_T','delta_q'] )\n", " # '''\n", " # MISSING: lc\n", " # '''\n", " OUT.to_csv(pathOUT+'Output_'+msk_id[twr]+'_'+twr+'_'+IOP+'_'+EM+'_new.csv',index=False)\n", " print('saved ' +'Output_'+msk_id[twr]+'_'+twr+'_'+IOP+'_'+EM+'.csv')\n", " \n", " #-- add output with values used to calculate ustar\n", " ustar_out = np.stack([wu,wv,w,u,v,wu_sgs,wv_sgs])\n", " np.save(pathOUT+'ustar_output'+IOP+'_'+EM+'_'+twr+'_new.npy',ustar_out)\n", "print('output stored in '+pathOUT)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }