diff --git a/scm/etc/scripts/scm_read_obs.py b/scm/etc/scripts/scm_read_obs.py index ac6b80100..ba9dc712d 100644 --- a/scm/etc/scripts/scm_read_obs.py +++ b/scm/etc/scripts/scm_read_obs.py @@ -7,249 +7,240 @@ import forcing_file_common as ffc def read_MOSAiC_obs(obs_file, time_slices, date): - obs_time_slice_indices = [] + obs_time_slice_indices = [] - obs_fid = Dataset(obs_file, 'r') + obs_fid = Dataset(obs_file, 'r') - obs_year = obs_fid.variables['year'][:] - obs_month = obs_fid.variables['month'][:] - obs_day = obs_fid.variables['day'][:] - obs_hour = obs_fid.variables['hour'][:] - obs_time = obs_fid.variables['time_offset'][:] + obs_year = obs_fid.variables['year'][:] + obs_month = obs_fid.variables['month'][:] + obs_day = obs_fid.variables['day'][:] + obs_hour = obs_fid.variables['hour'][:] + obs_time = obs_fid.variables['time_offset'][:] - obs_date = [] - for i in range(obs_hour.size): - obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) - obs_date = np.array(obs_date) + obs_date = [] + for i in range(obs_hour.size): + obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) + obs_date = np.array(obs_date) - for time_slice in time_slices: - start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) - end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) - start_date_index = np.where(obs_date == start_date)[0][0] - end_date_index = np.where(obs_date == end_date)[0][0] - obs_time_slice_indices.append([start_date_index, end_date_index]) + for time_slice in time_slices: + start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) + end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) + start_date_index = np.where(obs_date == start_date)[0][0] + end_date_index = np.where(obs_date == end_date)[0][0] + obs_time_slice_indices.append([start_date_index, end_date_index]) -#print(start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index]) + #print(start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index]) - #find the index corresponding to the start of the simulations - obs_start_index = np.where(obs_date == date[0][0])[0] - obs_time = obs_time - obs_time[obs_start_index] + #find the index corresponding to the start of the simulations + obs_start_index = np.where(obs_date == date[0][0])[0] + obs_time = obs_time - obs_time[obs_start_index] - obs_pres_l = obs_fid.variables['levels'][:]*100.0 #pressure levels in mb + obs_pres_l = obs_fid.variables['levels'][:]*100.0 #pressure levels in mb + + obs_T = obs_fid.variables['T'][:] + obs_q = obs_fid.variables['q'][:] + obs_qi= obs_fid.variables['qi'][:] + obs_ql= obs_fid.variables['ql'][:] + obs_u = obs_fid.variables['u'][:] + obs_v = obs_fid.variables['v'][:] + obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] + obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] + obs_tsk = obs_fid.variables['T_skin'][:] + obs_shf = obs_fid.variables['SH'][:] + obs_lhf = obs_fid.variables['LH'][:] + obs_t2m = obs_fid.variables['T_srf'][:] + obs_q2m = obs_fid.variables['q_srf'][:] + + obs_time_h = obs_time/3600.0 + + Rd = 287.0 + Rv = 461.0 + + e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) + e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity + obs_rh = np.clip(e/e_s, 0.0, 1.0) + + return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, + 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, + 'pres_l': obs_pres_l, 'T': obs_T, 'q': obs_q, 'rh': obs_rh, 'u': obs_u, 'v': obs_v, 'shf': obs_shf, + 'lhf': obs_lhf, 't2m': obs_t2m, 'q2m': obs_q2m, 'time_h': obs_time_h, 'tsfc': obs_tsk, + 'qv': obs_q, 'qi': obs_qi, 'ql': obs_ql, 'rad_net_srf': obs_rad_net_srf, 'sfc_dwn_lw': obs_lw_dn_srf} +# 'lwp': obs_lwp, 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} - obs_T = obs_fid.variables['T'][:] - obs_q = obs_fid.variables['q'][:] - obs_qi= obs_fid.variables['qi'][:] - obs_ql= obs_fid.variables['ql'][:] - obs_u = obs_fid.variables['u'][:] - obs_v = obs_fid.variables['v'][:] - obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] - obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] - obs_tsk = obs_fid.variables['T_skin'][:] - obs_shf = obs_fid.variables['SH'][:] - obs_lhf = obs_fid.variables['LH'][:] - obs_t2m = obs_fid.variables['T_srf'][:] - obs_q2m = obs_fid.variables['q_srf'][:] + obs_fid.close() - obs_time_h = obs_time/3600.0 + return return_dict - Rd = 287.0 - Rv = 461.0 +def read_twpice_obs(obs_file, time_slices, date): + obs_time_slice_indices = [] - e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) - e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity - obs_rh = np.clip(e/e_s, 0.0, 1.0) + obs_fid = Dataset(obs_file, 'r') - return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, - 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, - 'pres_l': obs_pres_l, 'T': obs_T, 'q': obs_q, 'rh': obs_rh, 'u': obs_u, 'v': obs_v, 'shf': obs_shf, - 'lhf': obs_lhf, 't2m': obs_t2m, 'q2m': obs_q2m, 'time_h': obs_time_h, 'tsfc': obs_tsk, - 'qv': obs_q, 'qi': obs_qi, 'ql': obs_ql, 'rad_net_srf': obs_rad_net_srf, 'sfc_dwn_lw': obs_lw_dn_srf} -# 'lwp': obs_lwp, 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} + obs_year = obs_fid.variables['year'][:] + obs_month = obs_fid.variables['month'][:] + obs_day = obs_fid.variables['day'][:] + obs_hour = obs_fid.variables['hour'][:] + obs_time = obs_fid.variables['time_offset'][:] - obs_fid.close() + obs_date = [] + for i in range(obs_hour.size): + obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) + obs_date = np.array(obs_date) - return return_dict + for time_slice in time_slices: + start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) + end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) + start_date_index = np.where(obs_date == start_date)[0][0] + end_date_index = np.where(obs_date == end_date)[0][0] + obs_time_slice_indices.append([start_date_index, end_date_index]) -def read_twpice_obs(obs_file, time_slices, date): - obs_time_slice_indices = [] + #find the index corresponding to the start of the simulations + obs_start_index = np.where(obs_date == date[0][0])[0] + obs_time = obs_time - obs_time[obs_start_index] - obs_fid = Dataset(obs_file, 'r') + obs_pres_l = obs_fid.variables['lev'][:]*100.0 #pressure levels in mb + + obs_cld = obs_fid.variables['cld'][:]/100.0 + obs_T = obs_fid.variables['T'][:] + obs_q = obs_fid.variables['q'][:]/1000.0 + obs_u = obs_fid.variables['u'][:] + obs_v = obs_fid.variables['v'][:] + obs_precip = obs_fid.variables['prec_srf'][:]/3.6E7 #convert from mm/hr to m/s + obs_shf = obs_fid.variables['SH'][:] + obs_lhf = obs_fid.variables['LH'][:] + obs_pwat = obs_fid.variables['PW'][:]*10.0 #convert from cm to kg/m2 + obs_lw_net_toa = obs_fid.variables['lw_net_toa'][:] + obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] + obs_sw_dn_toa = obs_fid.variables['sw_dn_toa'][:] + obs_sw_dn_srf = obs_fid.variables['sw_dn_srf'][:] + obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] + obs_lwp = obs_fid.variables['LWP'][:]*10.0 #convert from cm to kg/m2 + #obs_T_forcing = obs_fid.variables['dTdt'][:]*24.0 #convert from K/hour to K/day + #obs_q_forcing = obs_fid.variables['dqdt'][:]*24.0 #convert from g/kg/hour to g/kg/day + obs_h_advec_T = obs_fid.variables['T_adv_h'][:]*24.0 + obs_h_advec_q = obs_fid.variables['q_adv_h'][:]*24.0 + obs_v_advec_T = obs_fid.variables['T_adv_v'][:]*24.0 + obs_v_advec_q = obs_fid.variables['q_adv_v'][:]*24.0 + + obs_T_forcing = obs_h_advec_T + obs_v_advec_T + obs_q_forcing = obs_h_advec_q + obs_v_advec_q + + obs_time_h = obs_time/3600.0 + + Rd = 287.0 + Rv = 461.0 + + e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) + e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity + obs_rh = np.clip(e/e_s, 0.0, 1.0) + + obs_rh_500 = np.zeros(obs_rh.shape[0]) + index_500 = np.where(obs_pres_l[:]*0.01 < 500.0)[0][0] + lifrac = (obs_pres_l[index_500-1] - 50000.0)/(obs_pres_l[index_500-1] - obs_pres_l[index_500]) + for j in range(obs_rh.shape[0]): #loop over times + obs_rh_500[j] = obs_rh[j,index_500-1] + lifrac*(obs_rh[j,index_500] - obs_rh[j,index_500-1]) + #print index_500, pres_l[-1][j,index_500,k], pres_l[-1][j,index_500-1,k], rh_500_kj, rh[-1][j,index_500,k], rh[-1][j,index_500-1,k] + + return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, + 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, + 'pres_l': obs_pres_l, 'cld': obs_cld, 'T': obs_T, 'q': obs_q, 'u': obs_u, 'v': obs_v, + 'shf': obs_shf, 'lhf': obs_lhf, 'pwat': obs_pwat, 'time_h': obs_time_h, + 'tprcp_rate_accum': obs_precip, 'qv': obs_q, 'rh': obs_rh, 'rh_500': obs_rh_500, + 'lw_up_TOA_tot': obs_lw_net_toa, 'rad_net_srf': obs_rad_net_srf, 'sw_dn_TOA_tot': obs_sw_dn_toa, + 'lw_dn_sfc_tot': obs_lw_dn_srf, 'sw_dn_sfc_tot': obs_sw_dn_srf, 'lwp': obs_lwp, + 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} - obs_year = obs_fid.variables['year'][:] - obs_month = obs_fid.variables['month'][:] - obs_day = obs_fid.variables['day'][:] - obs_hour = obs_fid.variables['hour'][:] - obs_time = obs_fid.variables['time_offset'][:] + obs_fid.close() - obs_date = [] - for i in range(obs_hour.size): - obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) - obs_date = np.array(obs_date) + return return_dict - for time_slice in time_slices: +def read_arm_sgp_summer_1997_obs(obs_file, time_slices, date): + obs_time_slice_indices = [] + + obs_fid = Dataset(obs_file, 'r') + + obs_year = obs_fid.variables['Year'][:] + obs_month = obs_fid.variables['Month'][:] + obs_day = obs_fid.variables['Day'][:] + #obs_hour = obs_fid.variables['hour'][:] + obs_time = obs_fid.variables['time_offset'][:] + + #this file doesn't have the hour variable - calculate from the time offset (seconds from 00Z on 6/18/1997) + obs_hour = (((obs_time - 3)/3600.0)%24).astype(int) + + obs_date = [] + for i in range(obs_hour.size): + obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) + obs_date = np.array(obs_date) + + for time_slice in time_slices: start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) start_date_index = np.where(obs_date == start_date)[0][0] end_date_index = np.where(obs_date == end_date)[0][0] obs_time_slice_indices.append([start_date_index, end_date_index]) + #print start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index] + + #find the index corresponding to the start of the simulations + obs_start_index = np.where(obs_date == date[0][0])[0] + obs_time = obs_time - obs_time[obs_start_index] - #find the index corresponding to the start of the simulations - obs_start_index = np.where(obs_date == date[0][0])[0] - obs_time = obs_time - obs_time[obs_start_index] - - obs_pres_l = obs_fid.variables['lev'][:]*100.0 #pressure levels in mb - - obs_cld = obs_fid.variables['cld'][:]/100.0 - obs_T = obs_fid.variables['T'][:] - obs_q = obs_fid.variables['q'][:]/1000.0 - obs_u = obs_fid.variables['u'][:] - obs_v = obs_fid.variables['v'][:] - obs_precip = obs_fid.variables['prec_srf'][:]/3.6E7 #convert from mm/hr to m/s - obs_shf = obs_fid.variables['SH'][:] - obs_lhf = obs_fid.variables['LH'][:] - obs_pwat = obs_fid.variables['PW'][:]*10.0 #convert from cm to kg/m2 - obs_lw_net_toa = obs_fid.variables['lw_net_toa'][:] - obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] - obs_sw_dn_toa = obs_fid.variables['sw_dn_toa'][:] - obs_sw_dn_srf = obs_fid.variables['sw_dn_srf'][:] - obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] - obs_lwp = obs_fid.variables['LWP'][:]*10.0 #convert from cm to kg/m2 - #obs_T_forcing = obs_fid.variables['dTdt'][:]*24.0 #convert from K/hour to K/day - #obs_q_forcing = obs_fid.variables['dqdt'][:]*24.0 #convert from g/kg/hour to g/kg/day - obs_h_advec_T = obs_fid.variables['T_adv_h'][:]*24.0 - obs_h_advec_q = obs_fid.variables['q_adv_h'][:]*24.0 - obs_v_advec_T = obs_fid.variables['T_adv_v'][:]*24.0 - obs_v_advec_q = obs_fid.variables['q_adv_v'][:]*24.0 - - obs_T_forcing = obs_h_advec_T + obs_v_advec_T - obs_q_forcing = obs_h_advec_q + obs_v_advec_q - - obs_time_h = obs_time/3600.0 - - Rd = 287.0 - Rv = 461.0 - - e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) - e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity - obs_rh = np.clip(e/e_s, 0.0, 1.0) - - obs_rh_500 = np.zeros(obs_rh.shape[0]) - index_500 = np.where(obs_pres_l[:]*0.01 < 500.0)[0][0] - lifrac = (obs_pres_l[index_500-1] - 50000.0)/(obs_pres_l[index_500-1] - obs_pres_l[index_500]) - for j in range(obs_rh.shape[0]): #loop over times - obs_rh_500[j] = obs_rh[j,index_500-1] + lifrac*(obs_rh[j,index_500] - obs_rh[j,index_500-1]) - #print index_500, pres_l[-1][j,index_500,k], pres_l[-1][j,index_500-1,k], rh_500_kj, rh[-1][j,index_500,k], rh[-1][j,index_500-1,k] - - return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, - 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, - 'pres_l': obs_pres_l, 'cld': obs_cld, 'T': obs_T, 'q': obs_q, 'u': obs_u, 'v': obs_v, - 'shf': obs_shf, 'lhf': obs_lhf, 'pwat': obs_pwat, 'time_h': obs_time_h, - 'tprcp_rate_accum': obs_precip, 'qv': obs_q, 'rh': obs_rh, 'rh_500': obs_rh_500, - 'lw_up_TOA_tot': obs_lw_net_toa, 'rad_net_srf': obs_rad_net_srf, 'sw_dn_TOA_tot': obs_sw_dn_toa, - 'lw_dn_sfc_tot': obs_lw_dn_srf, 'sw_dn_sfc_tot': obs_sw_dn_srf, 'lwp': obs_lwp, - 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} - - obs_fid.close() - - return return_dict -def read_arm_sgp_summer_1997_obs(obs_file, time_slices, date): - obs_time_slice_indices = [] - - obs_fid = Dataset(obs_file, 'r') - - obs_year = obs_fid.variables['Year'][:] - obs_month = obs_fid.variables['Month'][:] - obs_day = obs_fid.variables['Day'][:] - #obs_hour = obs_fid.variables['hour'][:] - obs_time = obs_fid.variables['time_offset'][:] - - #this file doesn't have the hour variable - calculate from the time offset (seconds from 00Z on 6/18/1997) - obs_hour = (((obs_time - 3)/3600.0)%24).astype(int) - - obs_date = [] - for i in range(obs_hour.size): - obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0)) - obs_date = np.array(obs_date) - - for time_slice in time_slices: - start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) - end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) - start_date_index = np.where(obs_date == start_date)[0][0] - end_date_index = np.where(obs_date == end_date)[0][0] - obs_time_slice_indices.append([start_date_index, end_date_index]) - #print start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index] - - #find the index corresponding to the start of the simulations - obs_start_index = np.where(obs_date == date[0][0])[0] - obs_time = obs_time - obs_time[obs_start_index] - - - obs_pres_l = np.flipud(obs_fid.variables['lev'][:])*100.0 #pressure levels in mb - - obs_cld = np.fliplr(obs_fid.variables['ARSCL_Cld'][:,:,0,0])/100.0 - obs_T = np.fliplr(obs_fid.variables['Temp'][:,:,0,0]) - obs_q = np.fliplr(obs_fid.variables['H2O_Mixing_Ratio'][:,:,0,0]/1000.0) - obs_u = np.fliplr(obs_fid.variables['u_wind'][:,:,0,0]) - obs_v = np.fliplr(obs_fid.variables['v_wind'][:,:,0,0]) - obs_precip = obs_fid.variables['Prec'][:,0,0] - # obs_shf = obs_fid.variables['SH'][:] - # obs_lhf = obs_fid.variables['LH'][:] - # obs_pwat = obs_fid.variables['PW'][:] - # obs_lw_net_toa = obs_fid.variables['lw_net_toa'][:] - # obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] - # obs_sw_dn_toa = obs_fid.variables['sw_dn_toa'][:] - # obs_sw_dn_srf = obs_fid.variables['sw_dn_srf'][:] - # obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] - # obs_lwp = obs_fid.variables['LWP'][:]*10.0 #convert from cm to kg/m2 - # #obs_T_forcing = obs_fid.variables['dTdt'][:]*24.0 #convert from K/hour to K/day - # #obs_q_forcing = obs_fid.variables['dqdt'][:]*24.0 #convert from g/kg/hour to g/kg/day - # obs_h_advec_T = obs_fid.variables['T_adv_h'][:]*24.0 - # obs_h_advec_q = obs_fid.variables['q_adv_h'][:]*24.0 - # obs_v_advec_T = obs_fid.variables['T_adv_v'][:]*24.0 - # obs_v_advec_q = obs_fid.variables['q_adv_v'][:]*24.0 - # - # obs_T_forcing = obs_h_advec_T + obs_v_advec_T - # obs_q_forcing = obs_h_advec_q + obs_v_advec_q - # - # obs_time_h = obs_time/3600.0 - # - # Rd = 287.0 - # Rv = 461.0 - # - # e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) - # e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity - # obs_rh = np.clip(e/e_s, 0.0, 1.0) - # - # obs_rh_500 = np.zeros(obs_rh.shape[0]) - # index_500 = np.where(obs_pres_l[:]*0.01 < 500.0)[0][0] - # lifrac = (obs_pres_l[index_500-1] - 50000.0)/(obs_pres_l[index_500-1] - obs_pres_l[index_500]) - # for j in range(obs_rh.shape[0]): #loop over times - # obs_rh_500[j] = obs_rh[j,index_500-1] + lifrac*(obs_rh[j,index_500] - obs_rh[j,index_500-1]) - # #print index_500, pres_l[-1][j,index_500,k], pres_l[-1][j,index_500-1,k], rh_500_kj, rh[-1][j,index_500,k], rh[-1][j,index_500-1,k] - - return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, - 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, - 'pres_l': obs_pres_l, 'cld': obs_cld, 'T': obs_T, 'qv': obs_q, 'u': obs_u, 'v': obs_v, - 'precip': obs_precip}#, 'shf': obs_shf, 'lhf': obs_lhf, 'pwat': obs_pwat, 'time_h': obs_time_h, - # 'rain': obs_precip, 'rainc': obs_precip, 'qv': obs_q, 'rh': obs_rh, 'rh_500': obs_rh_500, - # 'lw_up_TOA_tot': obs_lw_net_toa, 'rad_net_srf': obs_rad_net_srf, 'sw_dn_TOA_tot': obs_sw_dn_toa, - # 'lw_dn_sfc_tot': obs_lw_dn_srf, 'sw_dn_sfc_tot': obs_sw_dn_srf, 'lwp': obs_lwp, - # 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} - - # return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, - # 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, - # 'pres_l': obs_pres_l, 'cld': obs_cld, 'T': obs_T, 'q': obs_q, 'u': obs_u, 'v': obs_v, - # 'precip': obs_precip, 'shf': obs_shf, 'lhf': obs_lhf, 'pwat': obs_pwat, 'time_h': obs_time_h, - # 'rain': obs_precip, 'rainc': obs_precip, 'qv': obs_q, 'rh': obs_rh, 'rh_500': obs_rh_500, - # 'lw_up_TOA_tot': obs_lw_net_toa, 'rad_net_srf': obs_rad_net_srf, 'sw_dn_TOA_tot': obs_sw_dn_toa, - # 'lw_dn_sfc_tot': obs_lw_dn_srf, 'sw_dn_sfc_tot': obs_sw_dn_srf, 'lwp': obs_lwp, - # 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} - - obs_fid.close() - - return return_dict + obs_pres_l = np.flipud(obs_fid.variables['lev'][:])*100.0 #pressure levels in mb + + obs_cld = np.fliplr(obs_fid.variables['ARSCL_Cld'][:,:,0,0])/100.0 + obs_T = np.fliplr(obs_fid.variables['Temp'][:,:,0,0]) + obs_q = np.fliplr(obs_fid.variables['H2O_Mixing_Ratio'][:,:,0,0]/1000.0) + obs_u = np.fliplr(obs_fid.variables['u_wind'][:,:,0,0]) + obs_v = np.fliplr(obs_fid.variables['v_wind'][:,:,0,0]) + obs_precip = obs_fid.variables['Prec'][:,0,0] + # obs_shf = obs_fid.variables['SH'][:] + # obs_lhf = obs_fid.variables['LH'][:] + # obs_pwat = obs_fid.variables['PW'][:] + # obs_lw_net_toa = obs_fid.variables['lw_net_toa'][:] + # obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:] + # obs_sw_dn_toa = obs_fid.variables['sw_dn_toa'][:] + # obs_sw_dn_srf = obs_fid.variables['sw_dn_srf'][:] + # obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:] + # obs_lwp = obs_fid.variables['LWP'][:]*10.0 #convert from cm to kg/m2 + # #obs_T_forcing = obs_fid.variables['dTdt'][:]*24.0 #convert from K/hour to K/day + # #obs_q_forcing = obs_fid.variables['dqdt'][:]*24.0 #convert from g/kg/hour to g/kg/day + # obs_h_advec_T = obs_fid.variables['T_adv_h'][:]*24.0 + # obs_h_advec_q = obs_fid.variables['q_adv_h'][:]*24.0 + # obs_v_advec_T = obs_fid.variables['T_adv_v'][:]*24.0 + # obs_v_advec_q = obs_fid.variables['q_adv_v'][:]*24.0 + # + # obs_T_forcing = obs_h_advec_T + obs_v_advec_T + # obs_q_forcing = obs_h_advec_q + obs_v_advec_q + # + # obs_time_h = obs_time/3600.0 + # + # Rd = 287.0 + # Rv = 461.0 + # + # e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa) + # e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity + # obs_rh = np.clip(e/e_s, 0.0, 1.0) + # + # obs_rh_500 = np.zeros(obs_rh.shape[0]) + # index_500 = np.where(obs_pres_l[:]*0.01 < 500.0)[0][0] + # lifrac = (obs_pres_l[index_500-1] - 50000.0)/(obs_pres_l[index_500-1] - obs_pres_l[index_500]) + # for j in range(obs_rh.shape[0]): #loop over times + # obs_rh_500[j] = obs_rh[j,index_500-1] + lifrac*(obs_rh[j,index_500] - obs_rh[j,index_500-1]) + # #print index_500, pres_l[-1][j,index_500,k], pres_l[-1][j,index_500-1,k], rh_500_kj, rh[-1][j,index_500,k], rh[-1][j,index_500-1,k] + + return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour, + 'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, + 'pres_l': obs_pres_l, 'cld': obs_cld, 'T': obs_T, 'qv': obs_q, 'u': obs_u, 'v': obs_v, + 'precip': obs_precip}#, 'shf': obs_shf, 'lhf': obs_lhf, 'pwat': obs_pwat, 'time_h': obs_time_h, + # 'rain': obs_precip, 'rainc': obs_precip, 'qv': obs_q, 'rh': obs_rh, 'rh_500': obs_rh_500, + # 'lw_up_TOA_tot': obs_lw_net_toa, 'rad_net_srf': obs_rad_net_srf, 'sw_dn_TOA_tot': obs_sw_dn_toa, + # 'lw_dn_sfc_tot': obs_lw_dn_srf, 'sw_dn_sfc_tot': obs_sw_dn_srf, 'lwp': obs_lwp, + # 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing} + + obs_fid.close() + + return return_dict def read_LASSO_obs(obs_file, time_slices, date): obs_time_slice_indices = []