diff --git a/components.yaml b/components.yaml index 642ab88a..664ce5a3 100644 --- a/components.yaml +++ b/components.yaml @@ -19,6 +19,13 @@ ecbuild: remote: ../ecbuild.git tag: geos/v1.3.0 +NCEP_Shared: + local: ./src/Shared/@NCEP_Shared + remote: ../NCEP_Shared.git + tag: v1.3.0 + sparse: ./config/NCEP_Shared.sparse + develop: main + GMAO_Shared: local: ./src/Shared/@GMAO_Shared remote: ../GMAO_Shared.git diff --git a/config/NCEP_Shared.sparse b/config/NCEP_Shared.sparse new file mode 100644 index 00000000..b0801a63 --- /dev/null +++ b/config/NCEP_Shared.sparse @@ -0,0 +1,8 @@ +!/* +/NCEP_sp +/NCEP_w3 +/NCEP_bufr +/NCEP_bacio +/NCEP_sfcio +/NCEP_sigio +/CMakeLists.txt diff --git a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml index 259eeb52..b2266730 100644 --- a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -83,7 +83,7 @@ fcsterr_inflation_fac = -9999. ! 0 = n/a [eg., in situ obs] ! 1 = ascending ! 2 = descending -! 3 = ascending and descending +! 3 = ascending or descending ! 4 = geostationary ! %pol = polarization ! 0 = n/a [eg., multi-pol. retrieval] @@ -2064,6 +2064,117 @@ obs_param_nml(48)%xcorr = 0.1875 obs_param_nml(48)%ycorr = 0.1875 obs_param_nml(48)%adapt = 0 +! -------------------------------------------------------------------- +! +! 49 = ASCAT_META_SM (ASCAT soil moisture ascending and descending orbits) +! +! https://navigator.eumetsat.int/product/EO:EUM:DAT:METOP:SOMO25 + +obs_param_nml(49)%descr = 'ASCAT_META_SM' +obs_param_nml(49)%orbit = 3 +obs_param_nml(49)%pol = 0 +obs_param_nml(49)%N_ang = 0 +obs_param_nml(49)%freq = 0 +obs_param_nml(49)%FOV = 20. +obs_param_nml(49)%FOV_units = 'km' +obs_param_nml(49)%assim = .false. +obs_param_nml(49)%scale = .false. +obs_param_nml(49)%getinnov = .false. +obs_param_nml(49)%RTM_ID = 0 +obs_param_nml(49)%bias_Npar = 0 +obs_param_nml(49)%bias_trel = 864000 +obs_param_nml(49)%bias_tcut = 432000 +obs_param_nml(49)%nodata = -9999. +obs_param_nml(49)%varname = 'sfds' +obs_param_nml(49)%units = '%' +obs_param_nml(49)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_A/' +obs_param_nml(49)%name = 'M02-ASCA-ASCSMO02' +obs_param_nml(49)%scalepath = '' +obs_param_nml(49)%scalename = '' +obs_param_nml(49)%flistpath = '' +obs_param_nml(49)%flistname = '' +obs_param_nml(49)%errstd = 9. +obs_param_nml(49)%std_normal_max = 2.5 +obs_param_nml(49)%zeromean = .true. +obs_param_nml(49)%coarsen_pert = .false. +obs_param_nml(49)%xcorr = 0.25 +obs_param_nml(49)%ycorr = 0.25 +obs_param_nml(49)%adapt = 0 + +! -------------------------------------------------------------------- +! +! 50 = ASCAT_METB_SM (ASCAT soil moisture ascending and descending orbits) +! +! https://navigator.eumetsat.int/product/EO:EUM:DAT:METOP:SOMO25 + +obs_param_nml(50)%descr = 'ASCAT_METB_SM' +obs_param_nml(50)%orbit = 3 +obs_param_nml(50)%pol = 0 +obs_param_nml(50)%N_ang = 0 +obs_param_nml(50)%freq = 0 +obs_param_nml(50)%FOV = 20. +obs_param_nml(50)%FOV_units = 'km' +obs_param_nml(50)%assim = .false. +obs_param_nml(50)%scale = .false. +obs_param_nml(50)%getinnov = .false. +obs_param_nml(50)%RTM_ID = 0 +obs_param_nml(50)%bias_Npar = 0 +obs_param_nml(50)%bias_trel = 864000 +obs_param_nml(50)%bias_tcut = 432000 +obs_param_nml(50)%nodata = -9999. +obs_param_nml(50)%varname = 'sfds' +obs_param_nml(50)%units = '%' +obs_param_nml(50)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_B/' +obs_param_nml(50)%name = 'M01-ASCA-ASCSMO02' +obs_param_nml(50)%scalepath = '' +obs_param_nml(50)%scalename = '' +obs_param_nml(50)%flistpath = '' +obs_param_nml(50)%flistname = '' +obs_param_nml(50)%errstd = 9. +obs_param_nml(50)%std_normal_max = 2.5 +obs_param_nml(50)%zeromean = .true. +obs_param_nml(50)%coarsen_pert = .false. +obs_param_nml(50)%xcorr = 0.25 +obs_param_nml(50)%ycorr = 0.25 +obs_param_nml(50)%adapt = 0 + +! -------------------------------------------------------------------- +! +! 51 = ASCAT_METC_SM (ASCAT soil moisture ascending and descending orbits) +! +! https://navigator.eumetsat.int/product/EO:EUM:DAT:METOP:SOMO25 + +obs_param_nml(51)%descr = 'ASCAT_METC_SM' +obs_param_nml(51)%orbit = 3 +obs_param_nml(51)%pol = 0 +obs_param_nml(51)%N_ang = 0 +obs_param_nml(51)%freq = 0 +obs_param_nml(51)%FOV = 20. +obs_param_nml(51)%FOV_units = 'km' +obs_param_nml(51)%assim = .false. +obs_param_nml(51)%scale = .false. +obs_param_nml(51)%getinnov = .false. +obs_param_nml(51)%RTM_ID = 0 +obs_param_nml(51)%bias_Npar = 0 +obs_param_nml(51)%bias_trel = 864000 +obs_param_nml(51)%bias_tcut = 432000 +obs_param_nml(51)%nodata = -9999. +obs_param_nml(51)%varname = 'sfds' +obs_param_nml(51)%units = '%' +obs_param_nml(51)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_C/' +obs_param_nml(51)%name = 'M03-ASCA-ASCSMO02' +obs_param_nml(51)%scalepath = '' +obs_param_nml(51)%scalename = '' +obs_param_nml(51)%flistpath = '' +obs_param_nml(51)%flistname = '' +obs_param_nml(51)%errstd = 9. +obs_param_nml(51)%std_normal_max = 2.5 +obs_param_nml(51)%zeromean = .true. +obs_param_nml(51)%coarsen_pert = .false. +obs_param_nml(51)%xcorr = 0.25 +obs_param_nml(51)%ycorr = 0.25 +obs_param_nml(51)%adapt = 0 + ! -------------------------------------------------------------------- / diff --git a/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m b/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m index 3c243b27..13272941 100644 --- a/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m +++ b/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m @@ -241,7 +241,7 @@ header = [y1 m1 d1 0 0 0 y2 m2 d2 0 0 0 tc.N_tile 1]; - tile_data = nanmean(L2_tau_tile([nidx_pre nidx nidx_nxt],:),1); + tile_data = mean(L2_tau_tile([nidx_pre nidx nidx_nxt],:),1,"omitnan"); tile_data(isnan(tile_data)) = 1.e15; fwrite( ifp, 14*4, int_precision ); % fortran_tag @@ -260,7 +260,7 @@ end % ======================== -% The final vegopacity.bin file contains data averaged (nanmean) of Asc +% The final vegopacity.bin file contains data averaged (mean) of Asc % and Des tau climatology. VOD is climatology,the other 2 parameters are % time constant with maximum spatial coverage data_clim_tile = NaN + ones(48, tc.N_tile,2); @@ -295,7 +295,7 @@ data_clim_tile(data_clim_tile < 0.) = 0.; % set small negative values to 0 % averaging A, D values -tile_data = nanmean(data_clim_tile,3); +tile_data = mean(data_clim_tile,3,"omitnan"); fname_out = strrep(fname, L2_Ascdes,'_AD_'); if fill_small_gaps diff --git a/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m b/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m index c717b9be..3125e892 100644 --- a/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m +++ b/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m @@ -4,7 +4,7 @@ L2_version,start_time, end_time) % function to compute the 2 RTM constant variables: Albedo and Roughness(H) based on the -% preprocessed data. Constants are taken as the long term temporal nanmean +% preprocessed data. Constants are taken as the long term temporal mean % with maximum coverage across Asc/Desc passes. % Q. Liu 18 Jul 2022 @@ -212,7 +212,7 @@ data_clim_tile(data_clim_tile < 0.) = 0.; % set small negative values to 0 % averaging A, D values - tile_data = nanmean(data_clim_tile,2); + tile_data = mean(data_clim_tile,2,"omitnan"); eval(['out_mwRTM.',out_Para{iPara},'=tile_data;']) end diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/Run_get_model_and_obs_clim_stats_latlon_grid.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/Run_get_model_and_obs_clim_stats_latlon_grid.m new file mode 100644 index 00000000..015f6762 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/Run_get_model_and_obs_clim_stats_latlon_grid.m @@ -0,0 +1,115 @@ +clear + +% ------------------------------------------------------------------- +% Begin user-defined inputs +% ------------------------------------------------------------------- + +% addpath('../../shared/matlab/'); +addpath('/discover/nobackup/amfox/current_GEOSldas/GEOSldas/src/Applications/LDAS_App/util/shared/matlab') + +% Define the Open Loop experiment path, run name, domain, and output prefix + +exp_path = '/discover/nobackup/amfox/Experiments/OLv7_M36_ascat'; +exp_run = {'OLv7_M36_ascat'}; +domain = 'SMAP_EASEv2_M36_GLOBAL'; +prefix_out = 'M36_zscore_stats_'; + +% Define the Open Loop experiment start and end dates + +start_month = 4; +start_year = 2015; +end_month = 3; +end_year = 2021; + +% Define the species names + +species_names = {'ASCAT_META_SM','ASCAT_METB_SM','ASCAT_METC_SM'}; + +% Define whether to combine species + +combine_species_stats = 1; % 1 to combine all species into single set of statistics + +% Define the grid resolution (degrees) + +grid_resolution = 0.25; + +% Define moving window size over which statistics are calculated, +% and minimum number of data points required to calculate statistics + +w_days = 75; +Ndata_min = 5; + +% Define the assimilation time step and initial time + +dt_assim = 3*60*60; +t0_assim = 0; + +% Define print intervals + +print_each_DOY = 1; +print_each_pentad = 0; +print_all_pentads = 1; + +% Define output directory (takes form "domain"/stats/"out_dir") +out_dir = 'z_score_clim_quarter_degree'; + +% Define the months to run over, 1:12, plus a number of months required to complete the window +run_months = [1:12 1:ceil(w_days/30)]; + +% ------------------------------------------------------------------- +% End user-defined inputs +% ------------------------------------------------------------------- + +% Calculate the earliest and latest years for each month in the experiment +earliest_year = zeros(length(run_months),1); +latest_year = zeros(length(run_months),1); + +cnt = 0; +for month = run_months + % Initialize the earliest and latest year variables + cnt = cnt + 1; + + % Check if the current year/month combination is earlier than the earliest + if datenum(start_year, month, 1) < datenum(start_year, start_month, 1) + earliest_year(cnt) = start_year+1; + else + earliest_year(cnt) = start_year; + end + + % Check if the current year/month combination is later than the latest + if datenum(end_year, month, 1) > datenum(end_year, end_month,1) + latest_year(cnt) = end_year-1; + else + latest_year(cnt) = end_year; + end +end + +% assume "ldas_obsparam" file is available at 0z on first day of start_month/start_year + +YYYY = num2str( start_year, '%4.4d' ); +MM = num2str( start_month, '%2.2d' ); + +obs_param_fname = [exp_path, '/', exp_run{1}, '/output/', domain, '/rc_out/Y', YYYY, ... + '/M', MM, '/',exp_run{1}, '.ldas_obsparam.', YYYY, MM, '01_0000z.txt']; + +[N_obs_param, obs_param ] = read_obsparam(obs_param_fname); + +species =[]; + +for i = 1:length(species_names) + add_species = obs_param(strcmp(species_names(i),{obs_param.descr})).species; + species = union(species,add_species); +end + +if combine_species_stats + disp('Calculating stats by combining multiple species'); +end + +% Calculate the climatology statistics + +get_model_and_obs_clim_stats_latlon_grid( species_names, run_months, exp_path, exp_run{1}, domain, earliest_year, ... + latest_year, dt_assim, t0_assim, species, combine_species_stats, ... + grid_resolution, w_days, Ndata_min, prefix_out, print_each_DOY, ... + print_each_pentad, print_all_pentads, out_dir ); + +% ================= EOF ========================================================================= diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m index 492839a2..8d8b242e 100644 --- a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m @@ -1,8 +1,8 @@ -function [] = get_model_and_obs_clim_stats( varname, ... +function [] = get_model_and_obs_clim_stats( varname, ... run_months, exp_path, exp_run, domain, start_year, end_year, ... - dt_assim, t0_assim, species, obs_param, ... - hscale, inc_angle, int_Asc, w_days, Ndata_min, prefix, ... + dt_assim, t0_assim, species, obs_param, ... + hscale, inc_angle, int_Asc, w_days, Ndata_min, prefix, ... convert_grid , time_of_day_in_hours ) % @@ -107,19 +107,20 @@ disp('ASSUMING EASEv2 M36 observations'); if ~isempty(strfind(domain,'M36')) && isempty(strfind(obs_param(species(1)).descr, '_E')) - tol = 1E-3; + tol = 1E-3; else - tol = 2; + tol = 2; end % output specs -overwrite = 1; +overwrite = 1; + +Nf = 5; %5 fields per polarization -Nf = 5; %5 fields per polarization N_out_fields = 2*Nf+4; %14; -write_ind_latlon = 'latlon_id'; %'latlon'; +write_ind_latlon = 'latlon_id'; %'latlon'; N_angle = length(inc_angle); N_pol = 2; @@ -159,46 +160,46 @@ D(1) = 1; P(1) = 1; if mi_m > 1 - D(1) = sum(days_in_month( 2014, [1:mi_m-1]))+1; - P(1) = ceil(D(1)/5); + D(1) = sum(days_in_month( 2014, [1:mi_m-1]))+1; + P(1) = ceil(D(1)/5); end if ma_m > 1 - D(2) = sum(days_in_month( 2014, [1:ma_m])); + D(2) = sum(days_in_month( 2014, [1:ma_m])); else - D(2) = 1; + D(2) = 1; end P(2) = floor(D(2)/5); if run_months(1) ~= run_months(end) && run_months(2) ~= run_months(end) - disp('WARNING: incomplete pentad-windows; loop through additional months to get complete pentads'); + disp('WARNING: incomplete pentad-windows; loop through additional months to get complete pentads'); end -fname_out_base = [ outpath, '/', prefix, ... - num2str(min(start_year)),'_doy',num2str(D(1)),'_',... - num2str(max(end_year)), '_doy',num2str(D(2)),... - '_hscale_', num2str(hscale,'%2.2f'), '_', ... - 'W_', num2str(w_days),'d_Nmin_', num2str(Ndata_min)]; +fname_out_base = [ outpath, '/', prefix, ... + num2str(min(start_year)),'_doy',num2str(D(1)),'_', ... + num2str(max(end_year)), '_doy',num2str(D(2)), ... + '_hscale_', num2str(hscale,'%2.2f'), '_', ... + 'W_', num2str(w_days),'d_Nmin_', num2str(Ndata_min)]; -fname_out_base_p = [ outpath, '/', prefix, ... - num2str(min(start_year)),'_p',num2str(P(1)),'_',... - num2str(max(end_year)), '_p',num2str(P(2)),... - '_hscale_', num2str(hscale,'%2.2f'), '_', ... - 'W_', num2str(round(w_days/5)),'p_Nmin_', num2str(Ndata_min)]; +fname_out_base_p = [ outpath, '/', prefix, ... + num2str(min(start_year)),'_p',num2str(P(1)),'_', ... + num2str(max(end_year)), '_p',num2str(P(2)), ... + '_hscale_', num2str(hscale,'%2.2f'), '_', ... + 'W_', num2str(round(w_days/5)),'p_Nmin_', num2str(Ndata_min)]; %fname_out_base = [fname_out_base, spec_tag]; if (int_Asc == 1) - Orbit_tag = '_A'; %'_Asc'; + Orbit_tag = '_A'; %'_Asc'; else - Orbit_tag = '_D'; %'_Desc'; + Orbit_tag = '_D'; %'_Desc'; end -fname_out_base = [fname_out_base, Orbit_tag]; +fname_out_base = [fname_out_base, Orbit_tag]; fname_out_base_p = [fname_out_base_p, Orbit_tag]; - + if exist( 'time_of_day_in_hours', 'var') - fname_out_base = [fname_out_base, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z']; + fname_out_base = [fname_out_base, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z']; fname_out_base_p = [fname_out_base_p, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z']; end @@ -210,7 +211,7 @@ fname = [inpath, '/rc_out/', exp_run, '.ldas_tilecoord.bin']; fnameg= [inpath, '/rc_out/', exp_run, '.ldas_tilegrids.bin']; -[ tile_coord ] = read_tilecoord( fname ); +[ tile_coord ] = read_tilecoord( fname ); [ tile_grid ] = read_tilegrids( fnameg ); N_tile = length(tile_coord.tile_id); @@ -226,101 +227,101 @@ tile_coord_tile_id = tile_coord.tile_id; if (exist('convert_grid')) - - %1) convert to M36 EASE indices - %2) convert back to lat/lon at center of obs - if (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv2'))) - gridid = 'M36'; - [central_row,central_col] = EASEv2_latlon2ind(central_lat,central_lon,gridid,1); - [central_lat,central_lon] = EASEv2_ind2latlon(central_row,central_col,gridid); - elseif (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv1'))) - error('Must provide smapeasev1_latlon2ind() and smapeasev1_ind2latlon()!') - gridid = 'M36'; - [central_row,central_col] = smapeasev1_latlon2ind(central_lat,central_lon,gridid); - [central_lat,central_lon] = smapeasev1_ind2latlon(central_row,central_col,gridid); - else - error(['Unable to convert to ',convert_grid]) - end - - row_col_tmp = [central_row central_col]; - [unique_rc, ia, ic] = unique(row_col_tmp,'rows'); - - max_Hx_c = length(find(mode(ic)==ic)); - - %know which exact M09 tiles are actually administering the obs - %------------------- - tmp_lon = central_lon(ia)+tmp_shift_lon; - tmp_lat = central_lat(ia)+tmp_shift_lat; - - [N_tile_in_cell_ij, tile_num_in_cell_ij] = get_tile_num_in_cell_ij(... - tile_coord, tile_grid); + + %1) convert to M36 EASE indices + %2) convert back to lat/lon at center of obs + if (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv2'))) + gridid = 'M36'; + [central_row,central_col] = EASEv2_latlon2ind(central_lat,central_lon,gridid,1); + [central_lat,central_lon] = EASEv2_ind2latlon(central_row,central_col,gridid); + elseif (~isempty(strfind(convert_grid, 'M36')) && ~isempty(strfind(convert_grid, 'EASEv1'))) + error('Must provide smapeasev1_latlon2ind() and smapeasev1_ind2latlon()!') + gridid = 'M36'; + [central_row,central_col] = smapeasev1_latlon2ind(central_lat,central_lon,gridid); + [central_lat,central_lon] = smapeasev1_ind2latlon(central_row,central_col,gridid); + else + error(['Unable to convert to ',convert_grid]) + end + + row_col_tmp = [central_row central_col]; + [unique_rc, ia, ic] = unique(row_col_tmp,'rows'); + + max_Hx_c = length(find(mode(ic)==ic)); + + %know which exact M09 tiles are actually administering the obs + %------------------- + tmp_lon = central_lon(ia)+tmp_shift_lon; + tmp_lat = central_lat(ia)+tmp_shift_lat; + + [N_tile_in_cell_ij, tile_num_in_cell_ij] = get_tile_num_in_cell_ij( ... + tile_coord, tile_grid); + + this_FOV = 20; + option = 'FOV_in_km'; + %overwrite ia with actual administering tile number + [ia] = get_tile_num_for_obs( tile_coord, tile_grid, ... + N_tile_in_cell_ij, tile_num_in_cell_ij, ... + option, this_FOV, tmp_lat, tmp_lon); + + ia = ia(ia>0 & ~isnan(ia)); + + obsnum = NaN+zeros(length(ic),1); + obsnum(ia) = [1:length(ia)]; + + N_tile_obs = length(ia); + + %------------------- + + if store_all_M09inM36 - this_FOV = 20; - option = 'FOV_in_km'; - %overwrite ia with actual administering tile number - [ia] = get_tile_num_for_obs(tile_coord, tile_grid,... - N_tile_in_cell_ij, tile_num_in_cell_ij,... - option, this_FOV, tmp_lat, tmp_lon); - - ia = ia(ia>0 & ~isnan(ia)); - - obsnum = NaN+zeros(length(ic),1); - obsnum(ia) = [1:length(ia)]; - - N_tile_obs = length(ia); - - %------------------- - - if store_all_M09inM36 - - %Not maintained/elaborated - tile_coord_tile_id = zeros(N_tile_obs,max_Hx_c); - - disp(['centralizing obs on ',convert_grid,' grid before doing stats: max ',num2str(max_Hx_c),'tiles per obs cell']) - - for i=1:N_tile_obs - - tmp_ind = find(row_col_tmp(:,1) == unique_rc(i,1) & row_col_tmp(:,2) == unique_rc(i,2)); - tile_coord_tile_id(i,1:length(tmp_ind)) = tile_coord.tile_id(tmp_ind); - - end - - else - - tile_coord_tile_id = tile_coord.tile_id(ia); - - end - + %Not maintained/elaborated + tile_coord_tile_id = zeros(N_tile_obs,max_Hx_c); + + disp(['centralizing obs on ',convert_grid,' grid before doing stats: max ',num2str(max_Hx_c),'tiles per obs cell']) + + for i=1:N_tile_obs + + tmp_ind = find(row_col_tmp(:,1) == unique_rc(i,1) & row_col_tmp(:,2) == unique_rc(i,2)); + + tile_coord_tile_id(i,1:length(tmp_ind)) = tile_coord.tile_id(tmp_ind); + + end + + else + + tile_coord_tile_id = tile_coord.tile_id(ia); + + end + else - - N_tile_obs = N_tile; - ia = 1:N_tile; - ic = 1:N_tile; - obsnum = 1:N_tile; - + + N_tile_obs = N_tile; + ia = 1:N_tile; + ic = 1:N_tile; + obsnum = 1:N_tile; + end -lon_out = tile_coord.com_lon(ia); %NaN+zeros(N_tile,1); -lat_out = tile_coord.com_lat(ia); %NaN+zeros(N_tile,1); - +lon_out = tile_coord.com_lon(ia); %NaN+zeros(N_tile,1); +lat_out = tile_coord.com_lat(ia); %NaN+zeros(N_tile,1); if hscale>0 - - for i=1:N_tile_obs - - this_lat = lat_out(i); - this_lon = lon_out(i); - - tmp_sq_distance = ... + + for i=1:N_tile_obs + + this_lat = lat_out(i); + this_lon = lon_out(i); + + tmp_sq_distance = ... (central_lon - this_lon).^2 + ... (central_lat - this_lat).^2; - - hscale_ind{i} = find( tmp_sq_distance <= hscale^2 ); - end - + + hscale_ind{i} = find( tmp_sq_distance <= hscale^2 ); + end + else - - hscale_ind = num2cell(ia); + + hscale_ind = num2cell(ia); end @@ -330,13 +331,13 @@ % N_pol and N_angle be specified here. Then subsample specifically % when the files are written out. -o_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); -m_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); -o_data2 = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); -m_data2 = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); -N_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); +o_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); +m_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); +o_data2 = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); +m_data2 = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); +N_data = NaN+zeros(N_pol,N_tile_obs,N_angle,w_days); -data_out = NaN+zeros(N_out_fields,N_tile_obs,N_angle); +data_out = NaN+zeros(N_out_fields,N_tile_obs,N_angle); % ------------------------------------------------------------- @@ -347,441 +348,381 @@ count = 0; for imonth = 1:length(run_months) - month = run_months(imonth); - -for day = 1:days_in_month( 2014, month) %2014 = random non-leap year + month = run_months(imonth); + + for day = 1:days_in_month( 2014, month) %2014 = random non-leap year + if count < w_days - count = count + 1; - else - count = w_days; - end - - for seconds_in_day = t0_assim:dt_assim:(86400-1) - - hour = floor(seconds_in_day/3600); - - % check if diurnal stats are needed - - if exist('time_of_day_in_hours','var') - tmp_hour = time_of_day_in_hours; + count = count + 1; else - tmp_hour = hour; % all hours of day will be included + count = w_days; end - - if hour==tmp_hour - - minute = floor( (seconds_in_day-hour*3600)/60 ); - - seconds = seconds_in_day-hour*3600-minute*60; - - if (seconds~=0) - input('something is wrong! Ctrl-c now') + + for seconds_in_day = t0_assim:dt_assim:(86400-1) + + hour = floor(seconds_in_day/3600); + + % check if diurnal stats are needed + + if exist('time_of_day_in_hours','var') + tmp_hour = time_of_day_in_hours; + else + tmp_hour = hour; % all hours of day will be included end - - - for year = start_year(imonth):end_year(imonth) - + + if hour==tmp_hour + + minute = floor( (seconds_in_day-hour*3600)/60 ); + + seconds = seconds_in_day-hour*3600-minute*60; + + if (seconds~=0) + input('something is wrong! Ctrl-c now') + end + + for year = start_year(imonth):end_year(imonth) + YYYYMMDD = [ num2str(year, '%4.4d'), ... - num2str(month, '%2.2d'), ... - num2str(day, '%2.2d') ]; - + num2str(month, '%2.2d'), ... + num2str(day, '%2.2d') ]; + HHMM = [ num2str(hour, '%2.2d'), ... - num2str(minute, '%2.2d') ]; - + num2str(minute, '%2.2d') ]; + % read innov files - - fname = [ inpath, '/ana/ens_avg/', ... - 'Y', YYYYMMDD(1:4), '/', ... - 'M', YYYYMMDD(5:6), '/', ... - exp_run, '.ens_avg.ldas_ObsFcstAna.', ... - YYYYMMDD, '_', HHMM, 'z.bin' ]; - - ifp = fopen( fname, 'r', 'l' ); - - if (ifp > 0) %Proceed only if file exists (e.g. irregular SMOS swaths!) - - fclose(ifp); - - [date_time, ... - obs_assim, ... - obs_species, ... - obs_tilenum, ... - obs_lon, ... - obs_lat, ... - obs_obs, ... - obs_obsvar, ... - obs_fcst, ... - obs_fcstvar, ... - obs_ana, ... - obs_anavar ... - ] = ... - read_ObsFcstAna( fname ); - - % remove tiles when there is no obs_fcst (obs_fcst == 0 in innov output when - % missing) - idx = find(obs_fcst == 0); - obs_assim(idx) = []; - obs_species(idx) = []; - obs_tilenum(idx) =[]; - obs_lon(idx) =[]; - obs_lat(idx) = []; - obs_obs(idx) = []; - obs_obsvar(idx) = []; - obs_fcst(idx) = []; - obs_fcstvar(idx) = []; - obs_ana(idx) = []; - obs_anavar(idx) = []; - - % extract species of interest + fname = [ inpath, '/ana/ens_avg/', ... + 'Y', YYYYMMDD(1:4), '/', ... + 'M', YYYYMMDD(5:6), '/', ... + exp_run, '.ens_avg.ldas_ObsFcstAna.', ... + YYYYMMDD, '_', HHMM, 'z.bin' ]; - ind = []; - - for this_species = species - - ind = find( obs_species == this_species); - - if (~isempty(ind)) - + ifp = fopen( fname, 'r', 'l' ); + + if (ifp > 0) % Proceed only if file exists (e.g. irregular SMOS swaths!) + + fclose(ifp); + + [ date_time, ... + obs_assim, ... + obs_species, ... + obs_tilenum, ... + obs_lon, ... + obs_lat, ... + obs_obs, ... + obs_obsvar, ... + obs_fcst, ... + obs_fcstvar, ... + obs_ana, ... + obs_anavar ... + ] = ... + read_ObsFcstAna( fname ); + + % remove tiles where obs_fcst is no-data (note: read_ObsFcstAna() returns NaN) + + idx = isnan(obs_fcst); + + obs_assim( idx) = []; + obs_species(idx) = []; + obs_tilenum(idx) = []; + obs_lon( idx) = []; + obs_lat( idx) = []; + obs_obs( idx) = []; + obs_obsvar( idx) = []; + obs_fcst( idx) = []; + obs_fcstvar(idx) = []; + obs_ana( idx) = []; + obs_anavar( idx) = []; + + % extract species of interest + + ind = []; + + for this_species = species + + ind = find( obs_species == this_species); + + if (~isempty(ind)) + obs_tilenum_i = obs_tilenum(ind); obs_obs_i = obs_obs(ind); obs_fcst_i = obs_fcst(ind); obs_lon_i = obs_lon(ind); obs_lat_i = obs_lat(ind); - + % Check if any location receives more than 1 obs (or 1 species) - + tmp = sort(obs_tilenum_i); same_tile = find(diff(tmp)==0); - + if (~isempty(same_tile)) - error('multiple obs of the same species at one location? - only last one in line is used'); + error('multiple obs of the same species at one location? - only last one in line is used'); end - + % Organize the data in a big matrix - + angle = obs_param(this_species == [obs_param.species]).ang; pol = obs_param(this_species == [obs_param.species]).pol; - - %pol intrinsically gives an index - %now find the index for the angle + + % pol intrinsically gives an index + % now find the index for the angle angle_i = find(angle(1) == inc_angle); - + % Only writes lat-lon at exact obs locations, but with % hscale>0, these obs are spread outside their exact % location. This allows to calculate stats at lan-lons % where no obs are available. - + %lon_out(obs_tilenum_i) = obs_lon_i; %lat_out(obs_tilenum_i) = obs_lat_i; - - %obs_lat/lon are the actual M36 lat/lons, *not* the - %administering tiles, so the lat/lons for the obs and those - %in the tile_coord would not be identical. - %Still, they should be in the - %neighbourhood, so check here if that is true. + + % obs_lat/lon are the actual M36 lat/lons, *not* the + % administering tiles, so the lat/lons for the obs and those + % in the tile_coord would not be identical. + % Still, they should be in the + % neighbourhood, so check here if that is true. if (any(abs(tile_coord.com_lat(obs_tilenum_i)-obs_lat_i) > tol) || ... any(abs(tile_coord.com_lon(obs_tilenum_i)-obs_lon_i) > tol) ) - error('Something wrong with tile_lat/lon') + error('Something wrong with tile_lat/lon') end - - %map model tiles (e.g. all M09) to observation administering - %tiles (could be a reduced subset of all M09) - %-------------------------------------------------------- + + % map model tiles (e.g. all M09) to observation administering + % tiles (could be a reduced subset of all M09) + % -------------------------------------------------------- obs_i = obsnum(obs_tilenum_i); - %-------------------------------------------------------- + % -------------------------------------------------------- if (hscale == 0) - - %11 May 2015: sum the obs and fcst within each day; - %and across years! - %some obs can be found at multiple hours within a day - %e.g. at the poles. - %**nansum of NaN's** result in zero, this need to be - %taken care of - o_data(pol(1),obs_i,angle_i,count) = nansum([o_data(pol(1),obs_i,angle_i,count); obs_obs_i' ]); - m_data(pol(1),obs_i,angle_i,count) = nansum([m_data(pol(1),obs_i,angle_i,count); obs_fcst_i']); - - %X^2 - o_data2(pol(1),obs_i,angle_i,count) = nansum([o_data2(pol(1),obs_i,angle_i,count); obs_obs_i'.^2 ]); - m_data2(pol(1),obs_i,angle_i,count) = nansum([m_data2(pol(1),obs_i,angle_i,count); obs_fcst_i'.^2]); - - %Sum of obs or model elements at each location - N_data(pol(1),obs_i,angle_i,count) = nansum([N_data(pol(1),obs_i,angle_i,count); ~isnan([obs_obs_i])']); - + + % 11 May 2015: sum the obs and fcst within each day; + % and across years! + % some obs can be found at multiple hours within a day + % e.g. at the poles. + % **sum(...,"omitnan") of NaNs** results in zero, this need to be + % taken care of + o_data( pol(1),obs_i,angle_i,count) = sum([o_data( pol(1),obs_i,angle_i,count); obs_obs_i' ], "omitnan"); + m_data( pol(1),obs_i,angle_i,count) = sum([m_data( pol(1),obs_i,angle_i,count); obs_fcst_i' ], "omitnan"); + + % X^2 + o_data2(pol(1),obs_i,angle_i,count) = sum([o_data2(pol(1),obs_i,angle_i,count); obs_obs_i'.^2 ], "omitnan"); + m_data2(pol(1),obs_i,angle_i,count) = sum([m_data2(pol(1),obs_i,angle_i,count); obs_fcst_i'.^2 ], "omitnan"); + + % Sum of obs or model elements at each location + N_data(pol(1), obs_i,angle_i,count) = sum([N_data( pol(1),obs_i,angle_i,count); ~isnan([obs_obs_i])'], "omitnan"); + else - - for i_ind = 1:length(obs_obs_i) - - %introduce a spatial effect of each observation on - %neighbouring statistics (through hscale) - s_eff = unique(hscale_ind{obs_i(i_ind)}); - %hscale_ind =[obs space] % - - %Sum of X - o_data(pol(1),s_eff,angle_i,count) = ... - nansum([o_data(pol(1),s_eff,angle_i,count); repmat(obs_obs_i(i_ind),1,length(s_eff))]); - m_data(pol(1),s_eff,angle_i,count) = ... - nansum([m_data(pol(1),s_eff,angle_i,count); repmat(obs_fcst_i(i_ind),1,length(s_eff))]); - - %Sum of X^2 - o_data2(pol(1),s_eff,angle_i,count) = ... - nansum([o_data2(pol(1),s_eff,angle_i,count); repmat(obs_obs_i(i_ind).^2,1,length(s_eff))]); - m_data2(pol(1),s_eff,angle_i,count) = ... - nansum([m_data2(pol(1),s_eff,angle_i,count); repmat(obs_fcst_i(i_ind).^2,1,length(s_eff))]); - - %Sum of obs or model elements at each location - N_data(pol(1),s_eff,angle_i,count) = ... - nansum([N_data(pol(1),s_eff,angle_i,count); repmat(~isnan([obs_obs_i(i_ind)]),1,length(s_eff)) ]); - - end - - end - - end - - end - + + for i_ind = 1:length(obs_obs_i) + + % introduce a spatial effect of each observation on + % neighbouring statistics (through hscale) + s_eff = unique(hscale_ind{obs_i(i_ind)}); + %hscale_ind =[obs space] % + + % Sum of X + o_data(pol(1),s_eff,angle_i,count) = ... + sum([o_data( pol(1),s_eff,angle_i,count); repmat( obs_obs_i( i_ind), 1,length(s_eff))], "omitnan"); + m_data(pol(1),s_eff,angle_i,count) = ... + sum([m_data( pol(1),s_eff,angle_i,count); repmat( obs_fcst_i(i_ind), 1,length(s_eff))], "omitnan"); + + % Sum of X^2 + o_data2(pol(1),s_eff,angle_i,count) = ... + sum([o_data2(pol(1),s_eff,angle_i,count); repmat( obs_obs_i( i_ind).^2,1,length(s_eff))], "omitnan"); + m_data2(pol(1),s_eff,angle_i,count) = ... + sum([m_data2(pol(1),s_eff,angle_i,count); repmat( obs_fcst_i(i_ind).^2,1,length(s_eff))], "omitnan"); + + % Sum of obs or model elements at each location + N_data(pol(1),s_eff,angle_i,count) = ... + sum([N_data( pol(1),s_eff,angle_i,count); repmat(~isnan([obs_obs_i(i_ind)]), 1,length(s_eff))], "omitnan"); + + end + + end % (hscale == 0) + + end % ~isempty(ind) + + end % species + end % if file present - - end % loop over multiple years - - end % time_of_day_in_hours - - end % seconds_in_day - - %count = count+1; - - if count >= w_days %wait initially until enough data is built up - - end_time.year = 2014; - end_time.month = month; - end_time.day = day; - end_time.hour = hour; - end_time.min = minute; - end_time.sec = seconds; - - start_time = augment_date_time( -floor(w_days*(24*60*60)), end_time ); - - % At the end of each day, collect the obs and fcst of the last - % w_day period, and write out a statistics-file at [w_day - floor(w_day/2)] - - o_data(abs(o_data - nodata) <= nodata_tol) = NaN; - m_data(abs(o_data - nodata) <= nodata_tol) = NaN; - - % data_out = zeros(N_out_fields,1:N_tiles,N_angle); - - for pol=[0 1] - - pp = pol*Nf; - - N_hscale_window = nansum(N_data(1+pol,:,:,1:w_days),4); - - if w_days == 95 - N_hscale_inner_window = nansum(N_data(1+pol,:,:,((w_days+1)/2-15):((w_days+1)/2+15)),4); - end + + end % loop over multiple years + + end % hour == tmp_hour (time_of_day_in_hours) - % OBSERVATIONS - %---------------- - %o_data is a sum over neighbouring obs above; - %here then take a sum over the time steps in the window - data_out(1+pp,:,:) = nansum(o_data(1+pol,:,:,1:w_days),4); - - %then make the average, by dividing over the sum of the number of - %timesteps and influencing obs at each location - data_out(1+pp,:,:) = data_out(1+pp,:,:)./N_hscale_window; - - %stdv_H = sqrt(E[X^2] - E[X]^2) - data_out(2+pp,:,:) = nansum(o_data2(1+pol,:,:,1:w_days),4); - data_out(2+pp,:,:) = data_out(2+pp,:,:)./N_hscale_window; - data_out(2+pp,:,:) = sqrt( data_out(2+pp,:,:) - data_out(1+pp,:,:).^2); - - % MODEL - %---------------- - data_out(3+pp,:,:) = nansum(m_data(1+pol,:,:,1:w_days),4); - data_out(3+pp,:,:) = data_out(3+pp,:,:)./N_hscale_window; - - data_out(4+pp,:,:) = nansum(m_data2(1+pol,:,:,1:w_days),4); - data_out(4+pp,:,:) = data_out(4+pp,:,:)./N_hscale_window; - data_out(4+pp,:,:) = sqrt( data_out(4+pp,:,:) - data_out(3+pp,:,:).^2); - - data_out(5+pp,:,:) = N_hscale_window; - - % Toss out stats that are based on too little data - - data_out([1:5]+pp,N_hscale_window= w_days %wait initially until enough data is built up + + end_time.year = 2014; + end_time.month = month; + end_time.day = day; + end_time.hour = hour; + end_time.min = minute; + end_time.sec = seconds; + + start_time = augment_date_time( -floor(w_days*(24*60*60)), end_time ); + + % At the end of each day, collect the obs and fcst of the last + % w_day period, and write out a statistics-file at [w_day - floor(w_day/2)] - if w_days == 95 + o_data(abs(o_data - nodata) <= nodata_tol) = NaN; + m_data(abs(o_data - nodata) <= nodata_tol) = NaN; + + % data_out = zeros(N_out_fields,1:N_tiles,N_angle); + + for pol=[0 1] + + pp = pol*Nf; + + N_hscale_window = sum(N_data(1+pol,:,:,1:w_days), 4,"omitnan"); + + if w_days == 95 + N_hscale_inner_window = sum(N_data(1+pol,:,:,((w_days+1)/2-15):((w_days+1)/2+15)),4,"omitnan"); + end + + % OBSERVATIONS + %---------------- + % o_data is a sum over neighbouring obs above; + % here then take a sum over the time steps in the window + data_out(1+pp,:,:) = sum( o_data( 1+pol,:,:,1:w_days),4,"omitnan"); + + % then make the average, by dividing over the sum of the number of + % timesteps and influencing obs at each location + data_out(1+pp,:,:) = data_out( 1+pp, :,:)./N_hscale_window; + + %stdv_H = sqrt(E[X^2] - E[X]^2) + data_out(2+pp,:,:) = sum( o_data2( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(2+pp,:,:) = data_out( 2+pp, :,:)./N_hscale_window; + data_out(2+pp,:,:) = sqrt( data_out( 2+pp, :,:) - data_out(1+pp,:,:).^2); + + % MODEL + %---------------- + data_out(3+pp,:,:) = sum( m_data( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(3+pp,:,:) = data_out( 3+pp, :,:)./N_hscale_window; + + data_out(4+pp,:,:) = sum( m_data2( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(4+pp,:,:) = data_out( 4+pp, :,:)./N_hscale_window; + data_out(4+pp,:,:) = sqrt( data_out( 4+pp, :,:) - data_out(3+pp,:,:).^2); + + data_out(5+pp,:,:) = N_hscale_window; + + % Toss out stats that are based on too little data + + data_out( [1:5]+pp,N_hscale_window < Ndata_min ) = NaN; + + if w_days == 95 data_out([1:5]+pp,N_hscale_inner_window < (Ndata_min/2.5)) = NaN; + end + end - end - -% fill_M09basedonM36 = 0; -% %X -% if ~exist('convert_grid') && fill_M09basedonM36 -% -% for kk = 1:size(data_out,1) -% -% idx = find(~isnan(data_out(kk,:))); -% for dd = 1:length(idx) -% -% this_lon = lon_out(idx(dd)); -% this_lat = lat_out(idx(dd)); -% this_data = data_out(kk, idx(dd)); -% -% -% -% [M36_row, M36_col] = smapeasev2_latlon2ind(this_lat,this_lon,'M36'); -% -% M09_row = (M36_row*4):((M36_row+1)*4-1); -% M09_col = (M36_col*4):((M36_col+1)*4-1); -% -% [M09_col, M09_row] = meshgrid(M09_col, M09_row); -% -% [M09_lat, M09_lon] = smapeasev2_ind2latlon(M09_row(:), M09_col(:), 'M09'); -% -% clear s_eff -% for ii = 1:length(M09_lat) -% tmp_sq_distance = ... -% (lon_out - M09_lon(ii)).^2 + ... -% (lat_out - M09_lat(ii)).^2; -% -% [tmp,s_eff(ii)] = min(tmp_sq_distance); -% -% -% %s_eff(ii) = intersect(find(lon_out>M09_lon(ii)-1e-4 & lon_outM09_lat(ii)-1e-4 & lat_out 1) -% disp('something is wrong, M36 might be overwritten') -% pause -% end -% -% data_out(kk,s_eff) = repmat(this_data, 1, length(s_eff)); -% -% end -% end -% end - - - % Get the actual obs/model at the center point (for debugging only!!) - - data_out(11,:,:) = o_data(1,:,:,w_days-floor(w_days/2.0))./N_data(1,:,:,w_days-floor(w_days/2.0)); - data_out(12,:,:) = m_data(1,:,:,w_days-floor(w_days/2.0))./N_data(1,:,:,w_days-floor(w_days/2.0)); - data_out(13,:,:) = o_data(2,:,:,w_days-floor(w_days/2.0))./N_data(2,:,:,w_days-floor(w_days/2.0)); - data_out(14,:,:) = m_data(2,:,:,w_days-floor(w_days/2.0))./N_data(2,:,:,w_days-floor(w_days/2.0)); - - % Get rid of NaN before writing a file - - data_out(isnan(data_out)) = nodata; - %lon_out(isnan(lon_out)) = nodata; - %lat_out(isnan(lat_out)) = nodata; - - % write output file - - date_time = end_time; - date_time = augment_date_time( -floor(w_days*(24*60*60)/2.0), date_time ); - - % always 365 files - - DOY = date_time.dofyr; - - if(is_leap_year(date_time.year) && DOY>=59) - - DOY = DOY-1; - - error('This code should never hit a leap year'); + % Get the actual obs/model at the center point (for debugging only!!) - end - - - fname_out = [fname_out_base, '_DOY', num2str(DOY,'%3.3d'), '.bin']; - - % check whether output file exists - - if (exist(fname_out)==2 && overwrite) - - disp(['output file exists. overwriting', fname_out]) - - elseif (exist(fname_out)==2 && ~overwrite) - - disp(['output file exists. not overwriting. returning']) - disp(['writing ', fname_out]) - return - - else - - disp(['creating ', fname_out]) - - end - - % compress data before writing in file. - - %idx_keep = find(any(abs(data_out -nodata) > nodata_tol,1)); - %lon_out_write = lon_out(idx_keep); - %lat_out_write = lat_out(idx_keep); - %data_out_write = data_out(:,idx_keep); - %tile_coord_tile_id_write = tile_coord_tile_id(idx_keep); - - - % write output for each DOY, sorted by all tiles - - if print_each_DOY + data_out(11,:,:) = o_data(1,:,:,w_days-floor(w_days/2.0))./N_data(1,:,:,w_days-floor(w_days/2.0)); + data_out(12,:,:) = m_data(1,:,:,w_days-floor(w_days/2.0))./N_data(1,:,:,w_days-floor(w_days/2.0)); + data_out(13,:,:) = o_data(2,:,:,w_days-floor(w_days/2.0))./N_data(2,:,:,w_days-floor(w_days/2.0)); + data_out(14,:,:) = m_data(2,:,:,w_days-floor(w_days/2.0))./N_data(2,:,:,w_days-floor(w_days/2.0)); + + % Get rid of NaN before writing a file + + data_out(isnan(data_out)) = nodata; + %lon_out(isnan(lon_out)) = nodata; + %lat_out(isnan(lat_out)) = nodata; + + % write output file + + date_time = end_time; + date_time = augment_date_time( -floor(w_days*(24*60*60)/2.0), date_time ); + + % always 365 files + + DOY = date_time.dofyr; + + if(is_leap_year(date_time.year) && DOY>=59) - write_seqbin_file(fname_out, lon_out, lat_out, ... - inc_angle, data_out(:,:,:), int_Asc, 0, ... %instead of writing the version#, write Ndata_min=0 - start_time, end_time, overwrite, ... - N_out_fields, write_ind_latlon, 'scaling',... - tile_coord_tile_id) - else + DOY = DOY-1; + + error('This code should never hit a leap year'); + + end + + + fname_out = [fname_out_base, '_DOY', num2str(DOY,'%3.3d'), '.bin']; - % if DOY is at middle of pentad, then copy the DOY to a pentad file - % DOY = pentad*5 - 2; ==> pentad = (DOY + 2)/5; - - pentad = (DOY + 2)/5; + % check whether output file exists + + if (exist(fname_out)==2 && overwrite) + + disp(['output file exists. overwriting', fname_out]) + + elseif (exist(fname_out)==2 && ~overwrite) + + disp(['output file exists. not overwriting. returning']) + disp(['writing ', fname_out]) + return + + else + + disp(['creating ', fname_out]) + + end - if mod((DOY + 2),5) == 0 + % write output for each DOY, sorted by all tiles - write_seqbin_file(fname_out, lon_out, lat_out, ... - inc_angle, data_out(:,:,:), int_Asc, 0, ... - start_time, end_time, overwrite, ... - N_out_fields, write_ind_latlon, 'scaling',... - tile_coord_tile_id) + if print_each_DOY + + write_seqbin_file(fname_out, lon_out, lat_out, ... + inc_angle, data_out(:,:,:), int_Asc, 0, ... % instead of writing the version#, write Ndata_min=0 + start_time, end_time, overwrite, ... + N_out_fields, write_ind_latlon, 'scaling', ... + tile_coord_tile_id) + else + + % if DOY is at middle of pentad, then copy the DOY to a pentad file + % DOY = pentad*5 - 2; ==> pentad = (DOY + 2)/5; + + pentad = (DOY + 2)/5; + + if mod((DOY + 2),5) == 0 + + write_seqbin_file(fname_out, lon_out, lat_out, ... + inc_angle, data_out(:,:,:), int_Asc, 0, ... + start_time, end_time, overwrite, ... + N_out_fields, write_ind_latlon, 'scaling', ... + tile_coord_tile_id) + + fname_out_p = [fname_out_base_p, '_p', num2str(pentad,'%2.2d'), '.bin']; + + copyfile(fname_out,fname_out_p); + + end - fname_out_p = [fname_out_base_p, '_p', num2str(pentad,'%2.2d'), '.bin']; - - copyfile(fname_out,fname_out_p); - end - end - - %clear idx_keep lon_out_write lat_out_write data_out_write tile_coord_tile_id_write - - % shift the window by one day and make room for the next day at the end - - o_data(:,:,:,1:w_days-1) = o_data(:,:,:,2:w_days); - m_data(:,:,:,1:w_days-1) = m_data(:,:,:,2:w_days); - o_data2(:,:,:,1:w_days-1) = o_data2(:,:,:,2:w_days); - m_data2(:,:,:,1:w_days-1) = m_data2(:,:,:,2:w_days); - N_data(:,:,:,1:w_days-1) = N_data(:,:,:,2:w_days); - - o_data(:,:,:,w_days) = NaN; - m_data(:,:,:,w_days) = NaN; - o_data2(:,:,:,w_days) = NaN; - m_data2(:,:,:,w_days) = NaN; - N_data(:,:,:,w_days) = NaN; - - data_out = NaN+0.0.*data_out; - - end - -end % day + %clear idx_keep lon_out_write lat_out_write data_out_write tile_coord_tile_id_write + + % shift the window by one day and make room for the next day at the end + + o_data( :,:,:,1:w_days-1) = o_data( :,:,:,2:w_days); + m_data( :,:,:,1:w_days-1) = m_data( :,:,:,2:w_days); + o_data2(:,:,:,1:w_days-1) = o_data2(:,:,:,2:w_days); + m_data2(:,:,:,1:w_days-1) = m_data2(:,:,:,2:w_days); + N_data( :,:,:,1:w_days-1) = N_data( :,:,:,2:w_days); + + o_data( :,:,:,w_days) = NaN; + m_data( :,:,:,w_days) = NaN; + o_data2(:,:,:,w_days) = NaN; + m_data2(:,:,:,w_days) = NaN; + N_data( :,:,:,w_days) = NaN; + + data_out = NaN+0.0.*data_out; + + end + + end % day end % month diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats_latlon_grid.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats_latlon_grid.m new file mode 100644 index 00000000..1d3a0c74 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats_latlon_grid.m @@ -0,0 +1,352 @@ + +function [] = get_model_and_obs_clim_stats_latlon_grid( species_names, ... + run_months, exp_path, exp_run, domain, start_year, end_year, ... + dt_assim, t0_assim, species, combine_species_stats, ... + resol, w_days, Ndata_min, prefix, print_each_DOY, ... + print_each_pentad, print_all_pentads,out_dir ) +% +% Adapted from get_model_and_obs_clim_stats.m +% +% Compute mean, stdv of model and observations from tile-based +% "innov" files for a selection of species on an Earth-fixed global +% lat/lon grid with resolution "resol". +% +% The main purpose of this function is to aggregate the information +% from the "innov" files so that the climatology statistics can +% be used in scaling of the observations before assimilation. +% +% One file with statistics is generated for every DOY (1,...,365). +% The temporal smoothing/averaging window (w_days) is given in days. +% +% We calculate the bias correction factors (scaling parameters) and +% write on an Earth-fixed lat/lon grid as there is no Earth-fixed +% regular grid for ASCAT observations. +% +% A. M. Fox - 27 Oct 2023 +% +% ------------------------------------------------------------------- +% begin user-defined inputs +% ------------------------------------------------------------------- + +nodata = -9999; +nodata_tol = 1e-4; +overwrite = 1; +Nf = 7; +N_pentads = 73; + +disp('ASSUMING obs are not on Earth-fixed regular grid (e.g., ASCAT)'); +disp(['Calculating scaling parameters on grid with resolution = ', num2str(resol) , ' degrees']); + +if combine_species_stats + N_species = 1; +else + N_species = length(species); +end + +inpath = [ exp_path, '/', exp_run, '/output/', domain ]; + +outpath = [ inpath, '/stats/', out_dir ]; + +% create outpath if it doesn't exist +if ~exist(outpath, 'dir') + mkdir(outpath); +end + +% assemble output file name +ind = start_year == min(start_year); +mi_m = min(run_months(ind)); +ind = end_year == max(end_year); +ma_m = max(run_months(ind)); + +D(1) = 1; +P(1) = 1; +if mi_m > 1 + D(1) = sum(days_in_month(2014, 1:mi_m-1)) + 1; + P(1) = ceil(D(1) / 5); +end +D(2) = sum(days_in_month(2014, 1:ma_m)); +P(2) = floor(D(2) / 5); + +fname_out_base_d = [outpath, '/', prefix, ... + num2str(min(start_year)), '_doy', num2str(D(1)), '_', ... + num2str(max(end_year)), '_doy', num2str(D(2)), ... + '_W_', num2str(w_days), 'd_Nmin_', num2str(Ndata_min)]; + +fname_out_base_p = [outpath, '/', prefix, ... + num2str(min(start_year)), '_p', num2str(P(1)), '_', ... + num2str(max(end_year)), '_p', num2str(P(2)), ... + '_W_', num2str(round(w_days/5)), 'p_Nmin_', num2str(Ndata_min)]; + +%====================================================== + +% Define lat/lon grid (dateline-on-edge, pole-on-edge) +% Define lower-left corner coordinates and grid cell size +ll_lon = -180; +ll_lat = -90; + +d_lon = resol; +d_lat = resol; + +% Calculate number of longitude and latitude grid cells +n_lon = round(360 / d_lon); +n_lat = round(180 / d_lat); + +% Calculate longitude and latitude values for the grid +ll_lons = linspace(ll_lon, ll_lon + (n_lon-1)*d_lon, n_lon); +ll_lats = linspace(ll_lat, ll_lat + (n_lat-1)*d_lat, n_lat); + +% Create grid index +obsnum = (1:n_lon*n_lat)'; +[i_out, j_out] = ind2sub([n_lon, n_lat], obsnum); +lon_out = ll_lons(i_out)'; +lat_out = ll_lats(j_out)'; +N_gridcells = length(obsnum); + +% initialize output statistics +o_data_sum = NaN(N_species, N_gridcells, w_days); +m_data_sum = NaN(N_species, N_gridcells, w_days); +o_data_sum2 = NaN(N_species, N_gridcells, w_days); +m_data_sum2 = NaN(N_species, N_gridcells, w_days); +m_data_min = NaN(N_species, N_gridcells, w_days); +m_data_max = NaN(N_species, N_gridcells, w_days); +N_data = NaN(N_species, N_gridcells, w_days); + +data_out = NaN(N_species, Nf, N_gridcells, N_pentads); +data2D = NaN(Nf, N_gridcells); + +% ------------------------------------------------------------- + +% make sure t0_assim is *first* analysis time in a day + +t0_assim = mod( t0_assim, dt_assim ); + +count = 0; + +for imonth = 1:length(run_months) + + month = run_months(imonth); + + for day = 1:days_in_month( 2014, month) %2014 = random non-leap year + + if count < w_days + count = count + 1; + else + count = w_days; + end + + for seconds_in_day = t0_assim:dt_assim:(86400-1) + + hour = floor( seconds_in_day/3600); + minute = floor((seconds_in_day-hour*3600)/60); + seconds = seconds_in_day-hour*3600-minute*60; + + if (seconds ~= 0) + input('something is wrong! Ctrl-c now') + end + + for year = start_year(imonth):end_year(imonth) + + YYYYMMDD = [num2str(year, '%4.4d'), num2str(month, '%2.2d'), num2str(day, '%2.2d')]; + HHMM = [num2str(hour, '%2.2d'), num2str(minute, '%2.2d')]; + + % read innov files + fname = [inpath, '/ana/ens_avg/', 'Y', YYYYMMDD(1:4), '/', 'M', YYYYMMDD(5:6), '/', exp_run, '.ens_avg.ldas_ObsFcstAna.', YYYYMMDD, '_', HHMM, 'z.bin']; + ifp = fopen(fname, 'r', 'l'); + + if (ifp > 0) % Proceed only if file exists + fclose(ifp); + [date_time, obs_assim, obs_species, obs_tilenum, obs_lon, obs_lat, obs_obs, obs_obsvar, obs_fcst, obs_fcstvar, obs_ana, obs_anavar] = read_ObsFcstAna(fname); + + % remove tiles where obs_fcst is no-data (note: read_ObsFcstAna() returns NaN) + + idx = isnan(obs_fcst); + + obs_assim( idx) = []; + obs_species(idx) = []; + obs_tilenum(idx) = []; + obs_lon( idx) = []; + obs_lat( idx) = []; + obs_obs( idx) = []; + obs_obsvar( idx) = []; + obs_fcst( idx) = []; + obs_fcstvar(idx) = []; + obs_ana( idx) = []; + obs_anavar( idx) = []; + + % extract species of interest + ind = []; + for scnt = 1:N_species + + if combine_species_stats + ind = find(ismember(obs_species, species)); + else + this_species = species(scnt); + ind = find(obs_species == this_species); + end + + if ~isempty(ind) + obs_tilenum_i = obs_tilenum(ind); + obs_obs_i = obs_obs( ind); + obs_fcst_i = obs_fcst( ind); + obs_lon_i = obs_lon( ind); + obs_lat_i = obs_lat( ind); + + % Check if any location receives more than 1 obs (or 1 species) + tmp = sort(obs_tilenum_i); + same_tile = find(diff(tmp) == 0, 1); + if ~isempty(same_tile) && ~combine_species_stats + error('multiple obs of the same species at one location? - only last one in line is used'); + end + + % Put obs lat/lon on our grid and figure out obsnum/grid index + i_idx = floor((obs_lon_i - ll_lon) / d_lon) + 1; + j_idx = floor((obs_lat_i - ll_lat) / d_lat) + 1; + [~, obs_idx] = ismember([i_idx, j_idx], [i_out, j_out], 'rows'); + obs_i = obsnum(obs_idx); + + o_data_sum( scnt, obs_i, count) = sum([o_data_sum( scnt, obs_i, count); obs_obs_i' ], "omitnan"); + m_data_sum( scnt, obs_i, count) = sum([m_data_sum( scnt, obs_i, count); obs_fcst_i' ], "omitnan"); + o_data_sum2(scnt, obs_i, count) = sum([o_data_sum2(scnt, obs_i, count); obs_obs_i'.^2 ], "omitnan"); + m_data_sum2(scnt, obs_i, count) = sum([m_data_sum2(scnt, obs_i, count); obs_fcst_i'.^2 ], "omitnan"); + m_data_min( scnt, obs_i, count) = min([m_data_min( scnt, obs_i, count); obs_fcst_i' ] ); + m_data_max( scnt, obs_i, count) = max([m_data_max( scnt, obs_i, count); obs_fcst_i' ] ); + N_data( scnt, obs_i, count) = sum([N_data( scnt, obs_i, count); ~isnan(obs_obs_i)' ], "omitnan"); + end + end + end + end + end + + if count >= w_days %wait initially until enough data is built up + end_time.year = 2014; + end_time.month = month; + end_time.day = day; + end_time.hour = hour; + end_time.min = minute; + end_time.sec = seconds; + + start_time = augment_date_time( -floor(w_days*(24*60*60)), end_time ); + + % At the end of each day, collect the obs and fcst of the last + % w_day period, and write out a statistics-file at [w_day - floor(w_day/2)] + o_data_sum(abs(o_data_sum - nodata) <= nodata_tol) = NaN; + m_data_sum(abs(m_data_sum - nodata) <= nodata_tol) = NaN; + + for i = 1:N_species + + N_window = sum(N_data( i,:,1:w_days), 3,"omitnan"); + + data2D(1,:) = sum(o_data_sum( i,:,1:w_days), 3,"omitnan")./N_window; + data2D(2,:) = sqrt(sum(o_data_sum2(i,:,1:w_days), 3,"omitnan")./N_window - data2D(1,:).^2); + data2D(3,:) = sum(m_data_sum( i,:,1:w_days), 3,"omitnan")./N_window; + data2D(4,:) = sqrt(sum(m_data_sum2(i,:,1:w_days), 3,"omitnan")./N_window - data2D(3,:).^2); + data2D(5,:) = N_window; + data2D(6,:) = min(m_data_min( i,:,1:w_days),[],3); % Want to use minimum mean daily value + data2D(7,:) = max(m_data_max( i,:,1:w_days),[],3); % Want to use maximum mean daily value + + % Set NaNs where there is not enough data + data2D([1:Nf],N_window=59) + error('This code should never hit a leap year'); + end + + if print_each_DOY + pentad = floor((DOY + 2)/5); + if combine_species_stats + fname_out = [fname_out_base_d, '_sp_ALL_DOY', num2str(DOY,'%3.3d'), '.nc4']; + else + fname_out = [fname_out_base_d,'_sp_', char(species_names(i)),'_DOY', num2str(DOY,'%3.3d'), '.nc4']; + end + if (exist(fname_out)==2 && overwrite) + disp(['Output file exists. overwriting', fname_out]) + elseif (exist(fname_out)==2 && ~overwrite) + disp(['Output file exists. not overwriting. returning']) + disp(['Writing ', fname_out]) + return + else + disp(['Creating ', fname_out]) + end + + % Write out the data + write_netcdf_latlon_grid( fname_out, i_out, j_out, ll_lons, ll_lats, data2D, pentad, ... + start_time, end_time, overwrite, Nf, ll_lon, ll_lat, d_lon, d_lat) + end + + if mod((DOY + 2),5) == 0 + pentad = (DOY + 2)/5; + data_out(i,:,:,pentad) = data2D; + start_time_p(pentad) = start_time; + end_time_p(pentad) = end_time; + if print_each_pentad + if combine_species_stats + fname_out = [fname_out_base_p, '_sp_ALL_p', num2str(pentad,'%2.2d'), '.nc4']; + else + fname_out = [fname_out_base_p, '_sp_', char(species_names(i)),'_p', num2str(pentad,'%2.2d'), '.nc4']; + end + if (exist(fname_out)==2 && overwrite) + disp(['Output file exists. overwriting', fname_out]) + elseif (exist(fname_out)==2 && ~overwrite) + disp(['Output file exists. not overwriting. returning']) + disp(['Writing ', fname_out]) + return + else + disp(['Creating ', fname_out]) + end + + % Write out the data + write_netcdf_latlon_grid( fname_out, i_out, j_out, ll_lons, ll_lats, data2D, pentad, ... + start_time, end_time, overwrite, Nf, ll_lon, ll_lat, d_lon, d_lat ) + end + end + + % Shift the data in the window to make room for next day + o_data_sum( i,:,1:w_days-1) = o_data_sum( i,:,2:w_days); + m_data_sum( i,:,1:w_days-1) = m_data_sum( i,:,2:w_days); + o_data_sum2(i,:,1:w_days-1) = o_data_sum2(i,:,2:w_days); + m_data_sum2(i,:,1:w_days-1) = m_data_sum2(i,:,2:w_days); + m_data_min( i,:,1:w_days-1) = m_data_min( i,:,2:w_days); + m_data_max( i,:,1:w_days-1) = m_data_max( i,:,2:w_days); + N_data( i,:,1:w_days-1) = N_data( i,:,2:w_days); + o_data_sum( i,:,w_days ) = NaN; + m_data_sum( i,:,w_days ) = NaN; + o_data_sum2(i,:,w_days ) = NaN; + m_data_sum2(i,:,w_days ) = NaN; + m_data_min( i,:,w_days ) = NaN; + m_data_max( i,:,w_days ) = NaN; + N_data( i,:,w_days ) = NaN; + + data2D = NaN+0.0.*data2D; + end + end % count >= w_days + end % day +end % month + +if print_all_pentads + for i = 1:N_species + data_o = squeeze(data_out(i,:,:,:)); + + if combine_species_stats + fname_out = [fname_out_base_d, '_sp_ALL_all_pentads.nc4']; + else + fname_out = [fname_out_base_d,'_sp_', char(species_names(i)),'_all_pentads.nc4']; + end + + if (exist(fname_out)==2 && overwrite) + disp(['Output file exists. overwriting', fname_out]) + elseif (exist(fname_out)==2 && ~overwrite) + disp(['Output file exists. not overwriting. returning']) + disp(['Writing ', fname_out]) + return + else + disp(['Creating ', fname_out]) + end + + write_netcdf_latlon_grid( fname_out, i_out, j_out, ll_lons, ll_lats, data_o, [1:73], ... + start_time_p, end_time_p, overwrite, Nf, ll_lon, ll_lat, d_lon, d_lat ) + end +end + + +% ==================== EOF ============================================== diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_tile_num_for_obs.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_tile_num_for_obs.m index ed1e9b8a..559853cf 100644 --- a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_tile_num_for_obs.m +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_tile_num_for_obs.m @@ -66,8 +66,10 @@ % map from i_ind, j_ind to tile_num - if ( ~isempty(strfind(tile_grid.gridtype, 'EASE_M')) || ... - ~isempty(strfind(tile_grid.gridtype, 'EASEv2_M')) ) + if ( ~isempty(strfind(tile_grid.gridtype, 'EASE_M')) || ... + ~isempty(strfind(tile_grid.gridtype, 'EASE-M')) || ... + ~isempty(strfind(tile_grid.gridtype, 'EASEv2-M')) || ... + ~isempty(strfind(tile_grid.gridtype, 'EASEv2_M')) ) % ASSUMPTION: tiles match EASE or EASEv2 grid cells exactly % (unless "outside" the domain, eg. water surface) diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/write_netcdf_latlon_grid.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/write_netcdf_latlon_grid.m new file mode 100644 index 00000000..49969ad5 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/write_netcdf_latlon_grid.m @@ -0,0 +1,232 @@ +function [] = write_netcdf_latlon_grid( fname, colind, rowind, ll_lons, ll_lats, ... + data, pentad, start_time, end_time, overwrite, N_out_fields, ll_lon, ll_lat, d_lon, d_lat ) + + int_precision = 'NC_INT'; % precision of fortran tag + float_precision = 'NC_DOUBLE'; % precision of data in input file + + % Define the compression level (0-9, where 0 is no compression and 9 is maximum compression) + compression_level = 5; + + version = 0; + + % check dimensions + if size(data,1)~=N_out_fields + error('ERROR: size of data incompatible with N_out_fields') + end + + % check for presence of optional input "overwrite" + if ~exist('overwrite','var') + overwrite = 0; % default: do NOT overwrite existing files + end + + % check if file exists + if exist(fname,'file') + if overwrite==0 + disp(['RETURNING!!! -- NOT OVERWRITING EXISTING FILE ', fname]) + return + else + disp(['OVERWRITING ', fname]) + end + else + disp(['writing ', fname]) + end + + % Convert the cell arrays to matrices using cell2mat (to deal with all_pentad case) + year_mat = cell2mat({start_time.year}'); + month_mat = cell2mat({start_time.month}'); + day_mat = cell2mat({start_time.day}'); + hour_mat = cell2mat({start_time.hour}'); + min_mat = cell2mat({start_time.min}'); + sec_mat = cell2mat({start_time.sec}'); + % Use the matrices as input to the datetime function + d = datetime(year_mat, month_mat, day_mat, hour_mat, min_mat, sec_mat); + % Convert to serial date number + serialNum = datenum(d); + % Subtract serial date number of January 1, 1950 + daysSince1950 = serialNum - datenum('January 1, 1950'); + tmp_start_time = daysSince1950; + + % Convert the cell arrays to matrices using cell2mat (to deal with all_pentad case) + year_mat = cell2mat({end_time.year}'); + month_mat = cell2mat({end_time.month}'); + day_mat = cell2mat({end_time.day}'); + hour_mat = cell2mat({end_time.hour}'); + min_mat = cell2mat({end_time.min}'); + sec_mat = cell2mat({end_time.sec}'); + % Use the matrices as input to the datetime function + d = datetime(year_mat, month_mat, day_mat, hour_mat, min_mat, sec_mat); + % Convert to serial date number + serialNum = datenum(d); + % Subtract serial date number of January 1, 1950 + daysSince1950 = serialNum - datenum('January 1, 1950'); + tmp_end_time = daysSince1950; + + N_lon = length(ll_lons); + N_lat = length(ll_lats); + + % Have we got multiple pentads + if ismatrix(data) + N_pentad = 1; + else + N_pentad = size(data,3); + end + + % create netCDF file + netcdf.setDefaultFormat('FORMAT_NETCDF4'); + ncid = netcdf.create(fname, 'NETCDF4'); + + % define dimensions + dimid_pentad = netcdf.defDim(ncid, 'pentad', N_pentad); + dimid_lon = netcdf.defDim(ncid, 'lon', N_lon); + dimid_lat = netcdf.defDim(ncid, 'lat', N_lat); + + % define variables + + varid_version = netcdf.defVar(ncid, 'version', int_precision, []); + + varid_ll_lon = netcdf.defVar(ncid, 'll_lon', float_precision, []); + netcdf.putAtt(ncid, varid_ll_lon, 'standard_name', 'longitude of lower left corner'); + netcdf.putAtt(ncid, varid_ll_lon, 'long_name', 'longitude of lower left corner'); + netcdf.putAtt(ncid, varid_ll_lon, 'units', 'degrees_east'); + netcdf.putAtt(ncid, varid_ll_lon, 'axis', 'X'); + + varid_ll_lat = netcdf.defVar(ncid, 'll_lat', float_precision, []); + netcdf.putAtt(ncid, varid_ll_lat, 'standard_name', 'latitude of lower left corner'); + netcdf.putAtt(ncid, varid_ll_lat, 'long_name', 'latitude of lower left corner'); + netcdf.putAtt(ncid, varid_ll_lat, 'units', 'degrees_north'); + netcdf.putAtt(ncid, varid_ll_lat, 'axis', 'Y'); + + varid_d_lon = netcdf.defVar(ncid, 'd_lon', float_precision, []); + netcdf.putAtt(ncid, varid_d_lon, 'standard_name', 'longitude grid spacing'); + netcdf.putAtt(ncid, varid_d_lon, 'long_name', 'longitude grid spacing'); + netcdf.putAtt(ncid, varid_d_lon, 'units', 'degrees'); + netcdf.putAtt(ncid, varid_d_lon, 'axis', 'X'); + + varid_d_lat = netcdf.defVar(ncid, 'd_lat', float_precision, []); + netcdf.putAtt(ncid, varid_d_lat, 'standard_name', 'latitude grid spacing'); + netcdf.putAtt(ncid, varid_d_lat, 'long_name', 'latitude grid spacing'); + netcdf.putAtt(ncid, varid_d_lat, 'units', 'degrees'); + netcdf.putAtt(ncid, varid_d_lat, 'axis', 'Y'); + + varid_pentad = netcdf.defVar(ncid, 'pentad', int_precision, [dimid_pentad]); + netcdf.putAtt(ncid, varid_pentad, 'standard_name', 'pentad'); + netcdf.putAtt(ncid, varid_pentad, 'long_name', 'pentad'); + netcdf.putAtt(ncid, varid_pentad, 'units', '1'); + netcdf.putAtt(ncid, varid_pentad, 'axis', 'T'); + + varid_start_time = netcdf.defVar(ncid, 'start_time', float_precision, [dimid_pentad]); + netcdf.putAtt(ncid, varid_start_time, 'standard_name', 'start time'); + netcdf.putAtt(ncid, varid_start_time, 'long_name', 'start time'); + netcdf.putAtt(ncid, varid_start_time, 'axis', 'T'); + netcdf.putAtt(ncid, varid_start_time, 'units', 'days since 1950-01-01 00:00:00.0 +0000'); + + varid_end_time = netcdf.defVar(ncid, 'end_time', float_precision, [dimid_pentad]); + netcdf.putAtt(ncid, varid_end_time, 'standard_name', 'end time'); + netcdf.putAtt(ncid, varid_end_time, 'long_name', 'end time'); + netcdf.putAtt(ncid, varid_end_time, 'axis', 'T'); + netcdf.putAtt(ncid, varid_end_time, 'units', 'days since 1950-01-01 00:00:00.0 +0000'); + + varid_om = netcdf.defVar(ncid, 'o_mean', float_precision, [dimid_lat dimid_lon dimid_pentad]); + netcdf.defVarDeflate(ncid,varid_om,true,true,compression_level); + netcdf.putAtt(ncid, varid_om, 'standard_name', 'observation mean'); + netcdf.putAtt(ncid, varid_om, 'long_name', 'Observation mean for pentad calculated over all years for window length'); + netcdf.putAtt(ncid, varid_om, 'units', 'Degree of saturation (0-1)'); + + varid_ov = netcdf.defVar(ncid, 'o_std', float_precision, [dimid_lat dimid_lon dimid_pentad]); + netcdf.defVarDeflate(ncid,varid_ov,true,true,compression_level); + netcdf.putAtt(ncid, varid_ov, 'standard_name', 'observation standard deviation'); + netcdf.putAtt(ncid, varid_ov, 'long_name', 'Observation standard deviation for pentad calculated over all years for window length'); + netcdf.putAtt(ncid, varid_ov, 'units', 'Degree of saturation (0-1)'); + + varid_mm = netcdf.defVar(ncid, 'm_mean', float_precision, [dimid_lat dimid_lon dimid_pentad]); + netcdf.defVarDeflate(ncid,varid_mm,true,true,compression_level); + netcdf.putAtt(ncid, varid_mm, 'standard_name', 'model mean'); + netcdf.putAtt(ncid, varid_mm, 'long_name', 'Model mean for pentad calculated over all years for window length'); + netcdf.putAtt(ncid, varid_mm, 'units', 'Surface soil moisture (m^3 m^-3)'); + + varid_mv = netcdf.defVar(ncid, 'm_std', float_precision, [dimid_lat dimid_lon dimid_pentad]); + netcdf.defVarDeflate(ncid,varid_mv,true,true,compression_level); + netcdf.putAtt(ncid, varid_mv, 'standard_name', 'model standard deviation'); + netcdf.putAtt(ncid, varid_mv, 'long_name', 'Model standard deviation for pentad calculated over all years for window length'); + netcdf.putAtt(ncid, varid_mv, 'units', 'Surface soil moisture (m^3 m^-3)'); + + varid_mi = netcdf.defVar(ncid, 'm_min', float_precision, [dimid_lat dimid_lon]); + netcdf.defVarDeflate(ncid,varid_mi,true,true,compression_level); + netcdf.putAtt(ncid, varid_mi, 'standard_name', 'model minimum'); + netcdf.putAtt(ncid, varid_mi, 'long_name', 'Model minimum calculated over all years'); + netcdf.putAtt(ncid, varid_mi, 'units', 'Surface soil moisture (m^3 m^-3)'); + + varid_ma = netcdf.defVar(ncid, 'm_max', float_precision, [dimid_lat dimid_lon]); + netcdf.defVarDeflate(ncid,varid_ma,true,true,compression_level); + netcdf.putAtt(ncid, varid_ma, 'standard_name', 'model maximum'); + netcdf.putAtt(ncid, varid_ma, 'long_name', 'Model maximum calculated over all years'); + netcdf.putAtt(ncid, varid_ma, 'units', 'Surface soil moisture (m^3 m^-3)'); + + varid_ndata = netcdf.defVar(ncid, 'n_data', float_precision, [dimid_lat dimid_lon dimid_pentad]); + netcdf.defVarDeflate(ncid,varid_ndata,true,true,compression_level); + netcdf.putAtt(ncid, varid_ndata, 'standard_name', 'number of data points'); + netcdf.putAtt(ncid, varid_ndata, 'long_name', 'Number of data points for pentad calculated over all years for window length'); + netcdf.putAtt(ncid, varid_ndata, 'units', '1'); + + % end define mode + netcdf.endDef(ncid); + + % write data + netcdf.putVar(ncid, varid_pentad, pentad); + netcdf.putVar(ncid, varid_start_time, tmp_start_time); + netcdf.putVar(ncid, varid_end_time, tmp_end_time); + + netcdf.putVar(ncid, varid_ll_lon, ll_lon); + netcdf.putVar(ncid, varid_ll_lat, ll_lat); + netcdf.putVar(ncid, varid_d_lon, d_lon); + netcdf.putVar(ncid, varid_d_lat, d_lat); + + if N_pentad ==1 + + data_out = ones(N_out_fields,N_lat,N_lon ) * -999.0; + + for n = 1:N_out_fields + for i = 1:length(colind) + data_out(n,rowind(i),colind(i)) = data(n,i); + end + end + + netcdf.putVar(ncid,varid_om, data_out(1,:,:) ); + netcdf.putVar(ncid,varid_ov, data_out(2,:,:) ); + netcdf.putVar(ncid,varid_mm, data_out(3,:,:) ); + netcdf.putVar(ncid,varid_mv, data_out(4,:,:) ); + netcdf.putVar(ncid,varid_mi, data_out(6,:,:) ); + netcdf.putVar(ncid,varid_ma, data_out(7,:,:) ); + netcdf.putVar(ncid,varid_ndata, data_out(5,:,:) ); + + else + + data_out = ones(N_out_fields,N_lat,N_lon,N_pentad) * -999.0; + + for n = 1:N_out_fields + for i = 1:length(colind) + data_out(n,rowind(i),colind(i),:) = data(n,i,:); + end + end + + netcdf.putVar(ncid,varid_om, data_out(1,:,:,:) ); + netcdf.putVar(ncid,varid_ov, data_out(2,:,:,:) ); + netcdf.putVar(ncid,varid_mm, data_out(3,:,:,:) ); + netcdf.putVar(ncid,varid_mv, data_out(4,:,:,:) ); + netcdf.putVar(ncid,varid_ndata, data_out(5,:,:,:) ); + netcdf.putVar(ncid,varid_ma, max(data_out(7,:,:,:),[],4)); % Max over all pentads, always only 2D + + min_data = squeeze(data_out(6, :, :,:)); + min_data(min_data < -9998) = NaN; % Switch current missing value to NaN before calculating min + min_data_out = min(min_data,[],3); + min_data_out(isnan(min_data_out)) = -9999.; + + netcdf.putVar(ncid,varid_mi, min_data_out); % Min over all pentads, always only 2D + end + + % close netCDF file + netcdf.close(ncid); + +end + +% ================ EOF ======================================================== diff --git a/src/Applications/LDAS_App/util/postproc/climatology/get_model_clim_stats.m b/src/Applications/LDAS_App/util/postproc/climatology/get_model_clim_stats.m index fef20cc8..7082a6c3 100644 --- a/src/Applications/LDAS_App/util/postproc/climatology/get_model_clim_stats.m +++ b/src/Applications/LDAS_App/util/postproc/climatology/get_model_clim_stats.m @@ -358,11 +358,11 @@ tmp = reshape(squeeze(hist_data(tile, s, :)),1,[]); - data_out(s,tile,1) = nanmean(tmp); % mean - data_out(s,tile,2) = nanstd(tmp); % stdv - data_out(s,tile,3) = min(tmp); % min - data_out(s,tile,4) = max(tmp); % max - data_out(s,tile,5) = sum(~isnan(tmp)); % N_data + data_out(s,tile,1) = mean( tmp,"omitnan"); % mean + data_out(s,tile,2) = std( tmp,"omitnan"); % stdv + data_out(s,tile,3) = min( tmp ); % min + data_out(s,tile,4) = max( tmp ); % max + data_out(s,tile,5) = sum(~isnan(tmp) ); % N_data % determine the CDF-parameters, or the edges for each % percentile diff --git a/src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m b/src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m index 7b14e671..7f510450 100644 --- a/src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m +++ b/src/Applications/LDAS_App/util/postproc/write_smapL4SMqa.m @@ -1258,7 +1258,6 @@ % % GDL, 16 Jan 2014 % GDL, 30 Jan 2014: remove NaN prior to the calculation of stats -% and avoid using nanmean, nanstdv, etc.. % GDL, 1 Feb 2014: - weighted statistics for mean and stdv % - add 6th statistic (frac): % the actual used fraction of N_data (based on weights) diff --git a/src/Applications/LDAS_App/util/shared/matlab/pentad_of_year.m b/src/Applications/LDAS_App/util/shared/matlab/pentad_of_year.m new file mode 100644 index 00000000..cc827841 --- /dev/null +++ b/src/Applications/LDAS_App/util/shared/matlab/pentad_of_year.m @@ -0,0 +1,14 @@ + +function pentad = pentad_of_year(day_of_year, year) + +if (is_leap_year(year) & day_of_year>=59) + + pentad = floor((day_of_year-2)/5)+1; + +else + + pentad = floor((day_of_year-1)/5)+1; + +end + +% ======================= EOF ================================== diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt index dc2312bc..13ca535e 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/CMakeLists.txt @@ -16,7 +16,7 @@ find_package(HDF5 REQUIRED COMPONENTS Fortran) esma_add_library (${this} SRCS ${SRCS} SUBCOMPONENTS ${alldirs} - DEPENDENCIES GEOS_LdasShared GEOSens_GridComp GEOSlandpert_GridComp GEOSland_GridComp makebcs MAPL GMAO_gfio_r4 hdf5hl_fortran hdf5_fortran ${NETCDF_LIBRARIES} + DEPENDENCIES GEOS_LdasShared GEOSens_GridComp GEOSlandpert_GridComp GEOSland_GridComp makebcs MAPL NCEP_bufr_r4i4 GMAO_gfio_r4 hdf5hl_fortran hdf5_fortran ${NETCDF_LIBRARIES} INCLUDES ${INC_ESMF} ${INC_HDF5}) target_compile_definitions (${this} PRIVATE LDAS_MPI) diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 index f8ba2682..f24124fd 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/GEOS_LandAssimGridComp.F90 @@ -1907,7 +1907,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) allocate(mwRTM_param(0)) - endif + endif call get_enkf_increments( & date_time_new, & diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 index d3cdcf16..b1a22716 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 @@ -44,7 +44,7 @@ module clsm_ensupd_glob_param ! total number of all obs species defined in namelist file ! (regardless of whether "assim" flag is true or false) - integer, parameter :: N_obs_species_nml = 48 + integer, parameter :: N_obs_species_nml = 51 ! ---------------------------------------------------------------------- ! diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 50e0fbf8..e7647e4c 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -1545,7 +1545,467 @@ subroutine read_obs_sm_ASCAT( & if (associated(tmp_lat)) deallocate(tmp_lat) end subroutine read_obs_sm_ASCAT + + ! **************************************************************************** + + subroutine read_obs_sm_ASCAT_EUMET( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, ASCAT_sm, ASCAT_sm_std, ASCAT_lon, ASCAT_lat, ASCAT_time ) + + !--------------------------------------------------------------------- + ! + ! Routine to read in ASCAT surface degree of saturation (sfds) obs from + ! BUFR files for both ascending and descending passes. + ! + ! ASCAT_sm and ASCAT_sm_std outputs from this subroutine are in wetness fraction (i.e., 0-1) units! + ! + ! Read in EUMETSAT level 2 soil moisture product 25 km (SMO), PPF software version 5.0. + ! Soil moisture derived from re-sampled (spatially averaged) backscatter (sigma0) values + ! on a 25-km orbit swath grid. Input data files are in BUFR file format. + ! + ! EUMETSAT BUFR files contain data for both ascending and descending half-orbits. + ! The BUFR field DOMO ("Direction of motion of moving observing platform") could be used to + ! separate Asc and Desc. (The BUFR files do not contain any explicit orbit indicator variable.) + ! According to Pamela Schoebel-Pattiselanno, EUMETSAT User Services Helpdesk: + ! "When the value (of DOMO) is between 180 and 270 degrees, it is the descending part + ! of the orbit. When it is between 270 and 360 degrees, it is the ascending part." + ! + ! Q. Liu, Nov 2019 - based on read_obs_sm_ASCAT + ! A. Fox, reichle, Sep 2023 - updated + ! + ! -------------------------------------------------------------------- + + implicit none + + ! inputs: + + type(date_time_type), intent(in) :: date_time + + integer, intent(in) :: dtstep_assim, N_catd + + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input + + type(grid_def_type), intent(in) :: tile_grid_d + + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & + N_tile_in_cell_ij + + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input + + type(obs_param_type), intent(in) :: this_obs_param + + ! outputs: + + logical, intent(out) :: found_obs + + real, intent(out), dimension(N_catd) :: ASCAT_sm ! sfds obs [fraction] (i.e., 0-1) + real, intent(out), dimension(N_catd) :: ASCAT_sm_std ! sfds obs err std [fraction] (i.e., 0-1) + real, intent(out), dimension(N_catd) :: ASCAT_lon, ASCAT_lat + real*8, intent(out), dimension(N_catd) :: ASCAT_time ! J2000 seconds + + ! --------------- + + ! Each obs file contains ~100-110 minutes (1 full orbit?) of observations (dt_ASCAT_obsfile). + ! File name indicates start time of swath. + + integer, parameter :: dt_ASCAT_obsfile = 110*60 ! seconds + + integer, parameter :: N_fnames_max = 15 ! max number of files per day + + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + + character( 15) :: str_date_time + character( 80) :: fname_of_fname_list + character(300) :: tmpfname + + type(date_time_type) :: date_time_tmp + type(date_time_type) :: date_time_low, date_time_low_fname + type(date_time_type) :: date_time_up + + integer :: ii, ind, N_tmp, N_files, kk, N_obs, N_fnames, N_fnames_tmp, obs_dir_hier + + character(300), dimension(:), allocatable :: fnames, tmpfnames + + ! -------------------- + ! + ! variables for BUFR read + + real*8, dimension(14) :: tmp_vdata + + integer, parameter :: lnbufr = 50 ! BUFR file unit number + integer, parameter :: max_rec = 50000 ! max number of obs after QC (expecting < 6 hr assim window) + integer, parameter :: max_obs = 250000 ! max number of obs read by subroutine (expecting < 6 hr assim window) + integer :: idate, iret, ireadmg, ireadsb + + character(8) :: subset + + ! -------------------- + + character(100), dimension(2*N_fnames_max) :: fname_list ! max 2 days of files + + real, dimension(:), allocatable :: tmp1_obs, tmp1_lat, tmp1_lon + real*8, dimension(:), allocatable :: tmp1_jtime + real*8, dimension(:,:), allocatable :: tmp_data + + real, dimension(:), pointer :: tmp_obs, tmp_lat, tmp_lon + real*8, dimension(:), pointer :: tmp_jtime + + integer, dimension(:), pointer :: tmp_tile_num + + integer, dimension(N_catd) :: N_obs_in_tile + + character(len=*), parameter :: Iam = 'read_obs_sm_ASCAT_EUMET' + character(len=400) :: err_msg + + ! ------------------------------------------------------------------- + + nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num, tmp_jtime ) + + ! --------------- + + ! initialize + + found_obs = .false. + + ! find files that are within half-open interval + ! (date_time-dtstep_assim/2,date_time+dtstep_assim/2] + + date_time_low = date_time + call augment_date_time( -(dtstep_assim/2), date_time_low) + date_time_up = date_time + call augment_date_time( (dtstep_assim/2), date_time_up) + + ! calculate "extra" date_time_low to catch files w/ swath start time stamps before window + ! but containing relevant obs + + date_time_low_fname = date_time_low + call augment_date_time( -dt_ASCAT_obsfile, date_time_low_fname) + + ! ---------------------------------------------------------------- + ! + ! read file with list of ASCAT file names for first day + + fname_of_fname_list = 'dummy' ! make sure it is properly defined in obs_param nml + + obs_dir_hier = 1 + + call read_obs_fnames( date_time_low_fname, this_obs_param, & + fname_of_fname_list, N_fnames_max, & + N_fnames, fname_list(1:N_fnames_max), obs_dir_hier ) + + ! if needed, read file with list of ASCAT file names for second day and add + ! file names into "fname_list" + + if (date_time_low_fname%day /= date_time_up%day) then + + call read_obs_fnames( date_time_up, this_obs_param, & + fname_of_fname_list, N_fnames_max, & + N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_fnames_max)), obs_dir_hier ) + + N_fnames = N_fnames + N_fnames_tmp + + end if + + tmpfnames = fname_list(1:N_fnames) + + ! ---------------------------------------------------------------- + ! + ! find files that have obs within assimilation window + + N_tmp = 0 + + do kk = 1,N_fnames + + tmpfname = fname_list(kk) + + ! Are we in the required assimilation window? + ! + ! e.g. Y2019/M07/M01-ASCA-ASCSMO02-NA-5.0-20190702075700.000000000Z-20190702084627-1350204.bfr + ! + ! 12345678901234567890123456789012345678901234567890123456789012345678901234567890 + ! 1 2 3 4 5 6 7 + + str_date_time = tmpfname(36:49) + + read(str_date_time( 1: 4), *) date_time_tmp%year + read(str_date_time( 5: 6), *) date_time_tmp%month + read(str_date_time( 7: 8), *) date_time_tmp%day + read(str_date_time( 9:10), *) date_time_tmp%hour + read(str_date_time(11:12), *) date_time_tmp%min + read(str_date_time(13:14), *) date_time_tmp%sec + + if ( datetime_lt_refdatetime( date_time_low_fname, date_time_tmp ) .and. & + datetime_le_refdatetime( date_time_tmp, date_time_up ) ) then + + N_tmp = N_tmp + 1 + + tmpfnames(N_tmp) = trim(this_obs_param%path) // '/' // trim(tmpfname) + + end if + + end do + + fnames = tmpfnames(1:N_tmp) + N_files = N_tmp + + ! ---------------------------------------------------------------- + ! + ! loop through files and read obs + metadata into tmp_data + + if (N_files>0) then + + ! read and process data if files are found + + allocate(tmp1_lon( max_rec )) + allocate(tmp1_lat( max_rec )) + allocate(tmp1_obs( max_rec )) + allocate(tmp1_jtime(max_rec )) + + allocate(tmp_data( max_obs,14)) + + N_obs = 0 + + do kk = 1,N_files + + ! open bufr file + + call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it + open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') + call openbf(lnbufr, 'SEC3', lnbufr) + call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) + call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) + + msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) + + loop_report: do while( ireadsb(lnbufr) == 0 ) + + ! columns of tmp_data: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 + + call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') + + N_obs = N_obs + 1 + + if (N_obs > max_obs) then + err_msg = 'Attempting to read too many obs - how long is your assimilation window?' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + tmp_data(N_obs,:) = tmp_vdata + + end do loop_report + end do msg_report + + call closbf(lnbufr) + close(lnbufr) + + end do ! end file loop + + ! ---------------------------------------------------------------- + ! + ! select obs within assimilation window and from desired orbit direction; apply basic QC based on obs info + + N_tmp = 0 + + do kk = 1,N_obs + + date_time_tmp%year = int(tmp_data(kk, 1)) + date_time_tmp%month = int(tmp_data(kk, 2)) + date_time_tmp%day = int(tmp_data(kk, 3)) + date_time_tmp%hour = int(tmp_data(kk, 4)) + date_time_tmp%min = int(tmp_data(kk, 5)) + date_time_tmp%sec = int(tmp_data(kk, 6)) + + ! skip if record outside of current assim window + if ( datetime_lt_refdatetime( date_time_tmp, date_time_low ) .and. & + datetime_le_refdatetime( date_time_up, date_time_tmp ) ) cycle + + ! skip if record contains invalid soil moisture value + if ( tmp_data(kk, 7) > 100. .or. tmp_data(kk, 7) < 0. ) cycle + + ! to distinguish orbit directions, must read "DOMO" from BUFR file + ! + ! 180 <= DOMO < 270 : descending + ! 270 < DOMO <= 360 : ascending + ! + ! if (index(this_obs_param%descr,'_A') /=0 .and. (tmp_data(kk, 8) < 270 .or. tmp_data(kk, 8) > 360)) cycle + ! if (index(this_obs_param%descr,'_D') /=0 .and. (tmp_data(kk, 8) < 180 .or. tmp_data(kk, 8) >= 270)) cycle + + ! skip if processing flag is set + if(int(tmp_data(kk, 8)) /= 0) cycle + + ! skip if correction flag is set ("good" values are 0 and 4) + if (.not. ( (int(tmp_data(kk, 9)) == 0) .or. (int(tmp_data(kk, 9)) == 4)) ) cycle + + ! skip if land fraction is missing or < 0.9 + if(tmp_data(kk, 10) > 1. .or. tmp_data(kk, 10) < 0.9 ) cycle + + ! skip if topographic complexity > 10% + if(tmp_data(kk, 11) > 10.) cycle + + ! skip if inundation and wetland fraction > 10% + if(tmp_data(kk, 12) > 10.) cycle + + N_tmp = N_tmp + 1 ! passed all QC + + if (N_tmp > max_rec) then + err_msg = 'Too many obs have passed QC - how long is your assimilation window?' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + tmp1_jtime(N_tmp) = datetime_to_J2000seconds( date_time_tmp, J2000_epoch_id ) + + tmp1_lat( N_tmp) = tmp_data(kk, 13) + tmp1_lon( N_tmp) = tmp_data(kk, 14) + + tmp1_obs( N_tmp) = tmp_data(kk, 7)/100. ! change units from percent (0-100) to fraction (0-1) + + end do + + if (logit) then + + write (logunit,*) 'read_obs_sm_ASCAT_EUMET: read ', N_tmp, ' at date_time = ', date_time, ' from:' + do ii=1,N_files + write (logunit,*) trim(fnames(ii)) + end do + write (logunit,*) '----------' + write (logunit,*) 'max(obs)=',maxval(tmp1_obs(1:N_tmp)), ', min(obs)=',minval(tmp1_obs(1:N_tmp)), & + ', avg(obs)=',sum(tmp1_obs(1:N_tmp))/N_tmp + + end if + + deallocate(fnames) + + ! copy "good" obs with lat/lon/time into tmp_* (pointers) + + allocate(tmp_jtime(N_tmp)) + allocate(tmp_lon( N_tmp)) + allocate(tmp_lat( N_tmp)) + allocate(tmp_obs( N_tmp)) + + tmp_jtime = tmp1_jtime(1:N_tmp) + tmp_lon = tmp1_lon( 1:N_tmp) + tmp_lat = tmp1_lat( 1:N_tmp) + tmp_obs = tmp1_obs( 1:N_tmp) + + deallocate(tmp1_jtime) + deallocate(tmp1_lon) + deallocate(tmp1_lat) + deallocate(tmp1_obs) + deallocate(tmp_data) + + else + + N_tmp = 0 + + end if + + ! ---------------------------------------------------------------- + ! + ! 2.) for each observation + ! a) determine grid cell that contains lat/lon + ! b) determine tile within grid cell that contains lat/lon + + if (N_tmp>0) then + + allocate(tmp_tile_num(N_tmp)) + + call get_tile_num_for_obs(N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + N_tmp, tmp_lat, tmp_lon, & + this_obs_param, & + tmp_tile_num ) + + + ! ---------------------------------------------------------------- + ! + ! 3.) compute super-obs for each tile from all obs w/in that tile + ! (also eliminate observations that are not in domain) + + ASCAT_sm = 0. + ASCAT_lon = 0. + ASCAT_lat = 0. + ASCAT_time = 0.0D0 + + N_obs_in_tile = 0 + + do ii=1,N_tmp + + ind = tmp_tile_num(ii) ! 1<=tmp_tile_num<=N_catd (unless nodata) + + if (ind>0) then ! this step eliminates obs outside domain + + ASCAT_sm( ind) = ASCAT_sm( ind) + tmp_obs( ii) + ASCAT_lon( ind) = ASCAT_lon( ind) + tmp_lon( ii) + ASCAT_lat( ind) = ASCAT_lat( ind) + tmp_lat( ii) + ASCAT_time(ind) = ASCAT_time(ind) + tmp_jtime(ii) + + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 + + end if + + end do + + ! normalize and set obs error std-dev + + do ii=1,N_catd + + ! set observation error standard deviation + + ASCAT_sm_std(ii) = this_obs_param%errstd/100. ! change units from percent (0-100) to fraction (0-1) + + ! normalize + + if (N_obs_in_tile(ii)>1) then + + ASCAT_sm( ii) = ASCAT_sm( ii)/real(N_obs_in_tile(ii)) + ASCAT_lon( ii) = ASCAT_lon( ii)/real(N_obs_in_tile(ii)) + ASCAT_lat( ii) = ASCAT_lat( ii)/real(N_obs_in_tile(ii)) + ASCAT_time( ii) = ASCAT_time(ii)/real(N_obs_in_tile(ii),kind(0.0D0)) + + elseif (N_obs_in_tile(ii)==0) then + + ASCAT_sm( ii) = this_obs_param%nodata + ASCAT_lon( ii) = this_obs_param%nodata + ASCAT_lat( ii) = this_obs_param%nodata + ASCAT_time( ii) = real(this_obs_param%nodata,kind(0.0D0)) + ASCAT_sm_std(ii) = this_obs_param%nodata + + else + + ! nothing to do if N_obs_in_tile(ii)==1 (and assuming N_obs_in_tile is never negative) + + end if + + end do + + ! clean up + + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) + + if (any(N_obs_in_tile>0)) then + + found_obs = .true. + + else + + found_obs = .false. + + end if + + end if + + ! clean up + + if (associated(tmp_obs)) deallocate(tmp_obs) + if (associated(tmp_lon)) deallocate(tmp_lon) + if (associated(tmp_lat)) deallocate(tmp_lat) + if (associated(tmp_jtime)) deallocate(tmp_jtime) + + end subroutine read_obs_sm_ASCAT_EUMET + ! *************************************************************************** subroutine read_sm_ASCAT_bin( & @@ -4811,7 +5271,7 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! read file with list of SMAP file names for first day - call read_obs_SMAP_fnames( date_time_low_fname, this_obs_param, & + call read_obs_fnames( date_time_low_fname, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames, fname_list(1:N_halforbits_max) ) @@ -4820,7 +5280,7 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & if (date_time_low_fname%day /= date_time_upp%day) then - call read_obs_SMAP_fnames( date_time_upp, this_obs_param, & + call read_obs_fnames( date_time_upp, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_halforbits_max)) ) @@ -5647,7 +6107,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! read file with list of SMAP file names for first day - call read_obs_SMAP_fnames( date_time_low_fname, this_obs_param, & + call read_obs_fnames( date_time_low_fname, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames, fname_list(1:N_halforbits_max) ) @@ -5656,7 +6116,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & if (date_time_low_fname%day /= date_time_upp%day) then - call read_obs_SMAP_fnames( date_time_upp, this_obs_param, & + call read_obs_fnames( date_time_upp, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_halforbits_max)) ) @@ -6683,14 +7143,24 @@ end subroutine turn_off_assim_SMAP_L1CTb ! ***************************************************************** - subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & - fname_of_fname_list, N_max, N_fnames, fname_list ) + subroutine read_obs_fnames( date_time, this_obs_param, & + fname_of_fname_list, N_max, N_fnames, fname_list, & + obs_dir_hier ) - ! read the file within a SMAP Yyyyy/Mmm/Ddd directory that lists - ! the SMAP h5 file names; preface file names with "Yyyyy/Mmm/Ddd" + ! read the file within an obs Yyyyy/Mmm/Ddd directory that lists + ! the obs file names, preface file names with "Yyyyy/Mmm/Ddd", + ! and return in "fname_list" + ! + ! optional input argument: + ! obs_dir_hier==1 : preface file names with "Yyyyy/Mmm" instead + ! + ! this subroutine is needed when obs file names cannot be predicted + ! and must be provided in a short text file that lists the file names + ! (e.g., SMAP Tb or soil moisture h5 files, ASCAT soil moisture BUFR files) ! ! reichle, 3 Jan 2014 ! reichle, 8 Jun 2017: Use "%flistpath" and "%flistname" from "obs_param_type". + ! A M Fox, reichle, 22 Sep 2023: added optional argument obs_dir_hier ! ! --------------------------------------------------------------------------------- @@ -6700,7 +7170,7 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & type(obs_param_type), intent(in) :: this_obs_param - character( *), intent(in) :: fname_of_fname_list + character( *), intent(in) :: fname_of_fname_list integer, intent(in) :: N_max @@ -6708,6 +7178,8 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & character(100), dimension(N_max), intent(out) :: fname_list + integer, optional, intent(in) :: obs_dir_hier + ! local variables character(300) :: fname @@ -6716,12 +7188,13 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & character( 80) :: tmpstr80 character( 14) :: YYYYMMDDdir + character( 10) :: YYYYMMdir character( 4) :: YYYY character( 2) :: MM, DD integer :: ii, istat - character(len=*), parameter :: Iam = 'read_obs_SMAP_fnames' + character(len=*), parameter :: Iam = 'read_obs_fnames' character(len=400) :: err_msg ! --------------------------------------------------------------------- @@ -6731,6 +7204,7 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & write (DD ,'(i2.2)') date_time%day YYYYMMDDdir = 'Y' // YYYY // '/M' // MM // '/D' // DD // '/' + YYYYMMdir = 'Y' // YYYY // '/M' // MM // '/' ! initialize default values @@ -6745,7 +7219,7 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & fname = trim(fpath_tmp) // '/' // YYYYMMDDdir // trim(fname_tmp) open( 10, file=fname, form='formatted', action='read', iostat=istat) - + if (istat/=0) then err_msg = 'cannot open file ' // trim(fname) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) @@ -6771,14 +7245,30 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - ! preface file names with "Yyyyy/Mmm/Ddd" - + ! preface file names with "Yyyyy/Mmm/Ddd" (default) + fname_list(ii) = YYYYMMDDdir // trim(tmpstr80) + if (present(obs_dir_hier)) then + + if (obs_dir_hier == 1) then + + ! preface file names with "Yyyyy/Mmm" + + fname_list(ii) = YYYYMMdir // trim(tmpstr80) + + else + + err_msg = 'Unrecognized value of optional argument obs_dir_hier' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end if + end if + else exit - + end if end do @@ -6787,7 +7277,7 @@ subroutine read_obs_SMAP_fnames( date_time, this_obs_param, & N_fnames = ii - end subroutine read_obs_SMAP_fnames + end subroutine read_obs_fnames ! ***************************************************************** @@ -7032,6 +7522,27 @@ subroutine read_obs( & tmp_obs, tmp_std_obs ) end if + + case ('ASCAT_META_SM','ASCAT_METB_SM','ASCAT_METC_SM' ) + + call read_obs_sm_ASCAT_EUMET( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & + tmp_time) + + ! scale observations to model climatology + + if (this_obs_param%scale .and. found_obs) then + + scaled_obs = .true. + + call scale_obs_sfmc_zscore( N_catd, tile_coord, & + date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & + tmp_obs, tmp_std_obs ) + + end if case ('isccp_tskin_gswp2_v1') @@ -7205,7 +7716,7 @@ subroutine read_obs( & 'SMAP_L1C_Tbh_E27_D', 'SMAP_L1C_Tbv_E27_D', & 'SMAP_L2AP_Tbh_A', 'SMAP_L2AP_Tbv_A', & 'SMAP_L2AP_Tbh_D', 'SMAP_L2AP_Tbv_D' ) - + call read_obs_SMAP_halforbit_Tb( & date_time, N_catd, this_obs_param, & dtstep_assim, tile_coord, tile_grid_d, & @@ -7222,7 +7733,7 @@ subroutine read_obs( & this_obs_param, tmp_obs, tmp_std_obs, tmp_assim ) end if - + case('LaRC_tskin-GOESW', 'LaRC_tskin-GOESE', 'LaRC_tskin-MET09', & 'LaRC_tskin-FY2E-', 'LaRC_tskin-MTST2') @@ -7645,6 +8156,251 @@ subroutine scale_obs_tskin_zscore( N_catd, tile_coord, & end subroutine scale_obs_tskin_zscore + + ! ***************************************************************** + + subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & + date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & + tmp_obs, tmp_std_obs ) + + ! scale sfmc or sfds obs to model climatology via standard-normal-deviate (zscore) + ! scaling using seasonally varying (pentad) stats + ! Grid information is read from a NetCDF file + ! + ! Scaling parameters are read from a NetCDF file that contains the following: + ! variables: + ! int version ; + ! double ll_lon ; + ! ll_lon:standard_name = "longitude of lower left corner" ; + ! double ll_lat ; + ! ll_lat:standard_name = "latitude of lower left corner" ; + ! double d_lon ; + ! d_lon:standard_name = "longitude grid spacing" ; + ! double d_lat ; + ! d_lat:standard_name = "latitude grid spacing" ; + ! int pentad(pentad) ; + ! pentad:standard_name = "pentad" ; + ! double start_time(pentad) ; + ! start_time:standard_name = "start time" ; + ! double end_time(pentad) ; + ! end_time:standard_name = "end time" ; + ! double o_mean(pentad, lon, lat) ; + ! o_mean:standard_name = "observation mean" ; + ! double o_std(pentad, lon, lat) ; + ! o_std:standard_name = "observation standard deviation" ; + ! double m_mean(pentad, lon, lat) ; + ! m_mean:standard_name = "model mean" ; + ! double m_std(pentad, lon, lat) ; + ! m_std:standard_name = "model standard deviation" ; + ! double m_min(lon, lat) ; + ! m_min:standard_name = "model minimum" ; + ! double m_max(lon, lat) ; + ! m_max:standard_name = "model maximum" ; + ! double n_data(pentad, lon, lat) ; + ! n_data:standard_name = "number of data points" ; + ! + ! A M Fox, reichle, April 2023 + + use netcdf + implicit none + + integer, intent(in) :: N_catd + + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input + + type(date_time_type), intent(in) :: date_time + + type(obs_param_type), intent(in) :: this_obs_param + + real, intent(in), dimension(N_catd) :: tmp_lon, tmp_lat + real*8, intent(in), dimension(N_catd) :: tmp_time ! J2000 seconds + + ! inout + + real, intent(inout), dimension(N_catd) :: tmp_obs + real, intent(inout), dimension(N_catd) :: tmp_std_obs + + ! ------------------- + + character(300) :: fname + + integer :: i, ind, pp, j_ind, i_ind + integer :: ncid, varid, ierr, ierr2 + integer :: pentad_dimid, lon_dimid, lat_dimid + integer :: N_pentad, N_lon, N_lat + integer :: pentad_varid, lon_varid, lat_varid + integer :: o_mean_varid, o_std_varid, m_mean_varid, m_std_varid, m_min_varid, m_max_varid + integer :: ll_lon_varid, ll_lat_varid, dlon_varid, dlat_varid + integer, dimension(3) :: start, icount + + logical :: file_exists + + real :: tmpreal, this_lon, this_lat, ll_lon, ll_lat, dlon, dlat + + integer, dimension(:), allocatable :: sclprm_tile_id + integer, dimension(:), allocatable :: pentads + + real, dimension(:,:), allocatable :: sclprm_mean_obs, sclprm_std_obs + real, dimension(:,:), allocatable :: sclprm_mean_mod, sclprm_std_mod + real, dimension(:,:), allocatable :: sclprm_min_mod, sclprm_max_mod + + character(len=*), parameter :: Iam = ' scale_obs_sfmc_zscore' + character(len=400) :: err_msg + + ! ------------------------------------------------------------------ + + ! Read scaling parameters from file + + fname = trim(this_obs_param%scalepath) // '/' // & + trim(this_obs_param%scalename) // '.nc4' + + if (logit) write (logunit,*) 'scaling obs species ', this_obs_param%species, ':' + if (logit) write (logunit,'(400A)') ' reading ', trim(fname) + + ! Check if file exists + + inquire(file=fname, exist=file_exists) + + if (.not. file_exists) then + + err_msg = 'Scaling parameter file not found' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end if + + ! Determine pentad to use + + pp = date_time%pentad + + ! Open the NetCDF file + + ierr = nf90_open(fname, nf90_nowrite, ncid) + + ! Get the dimension and variable IDs + + ierr = nf90_inq_varid(ncid, 'll_lon', ll_lon_varid) + ierr = nf90_inq_varid(ncid, 'll_lat', ll_lat_varid) + ierr = nf90_inq_varid(ncid, 'd_lon', dlon_varid) + ierr = nf90_inq_varid(ncid, 'd_lat', dlat_varid) + + ! Get the dimension sizes + + ierr = nf90_inq_dimid(ncid, 'pentad', pentad_dimid) + ierr = nf90_inq_dimid(ncid, 'lon', lon_dimid) + ierr = nf90_inq_dimid(ncid, 'lat', lat_dimid) + + ierr = nf90_inquire_dimension(ncid, pentad_dimid, len = N_pentad) + ierr = nf90_inquire_dimension(ncid, lon_dimid, len = N_lon) + ierr = nf90_inquire_dimension(ncid, lat_dimid, len = N_lat) + + ! Get the variable IDs + + ierr = nf90_inq_varid(ncid, 'o_mean', o_mean_varid) + ierr = nf90_inq_varid(ncid, 'o_std', o_std_varid) + ierr = nf90_inq_varid(ncid, 'm_mean', m_mean_varid) + ierr = nf90_inq_varid(ncid, 'm_std', m_std_varid) + ierr = nf90_inq_varid(ncid, 'm_min', m_min_varid) + ierr = nf90_inq_varid(ncid, 'm_max', m_max_varid) + + ! Read grid variables + + ierr = nf90_get_var(ncid, ll_lon_varid, ll_lon) + ierr = nf90_get_var(ncid, ll_lat_varid, ll_lat) + ierr = nf90_get_var(ncid, dlon_varid, dlon) + ierr = nf90_get_var(ncid, dlat_varid, dlat) + + start = [1, 1, pp] + icount = [N_lat, N_lon, 1 ] + + ! Read mean and std variables + + allocate(sclprm_mean_obs(N_lat, N_lon), sclprm_std_obs(N_lat, N_lon)) + allocate(sclprm_mean_mod(N_lat, N_lon), sclprm_std_mod(N_lat, N_lon)) + allocate(sclprm_min_mod( N_lat, N_lon), sclprm_max_mod(N_lat, N_lon)) + + ierr = nf90_get_var(ncid, o_mean_varid, sclprm_mean_obs, start, icount) + ierr = nf90_get_var(ncid, o_std_varid, sclprm_std_obs, start, icount) + ierr = nf90_get_var(ncid, m_mean_varid, sclprm_mean_mod, start, icount) + ierr = nf90_get_var(ncid, m_std_varid, sclprm_std_mod, start, icount) + ierr = nf90_get_var(ncid, m_min_varid, sclprm_min_mod) + ierr = nf90_get_var(ncid, m_max_varid, sclprm_max_mod) + + ! Close the netcdf file + + ierr = nf90_close(ncid) + + ! -------------------------------------------------------------- + + ! Scale observations (at this point all obs are of same type because + ! of the way the subroutine is called from subroutine read_obs() + + do i=1,N_catd + + ! Check for no-data-values in observation (any neg value is no_data) + + if (tmp_obs(i)>=0.) then + + ! ll_lon and ll_lat refer to lower left corner of grid cell + ! (as opposed to the grid point in the center of the grid cell) + + this_lon = tmp_lon(i) + this_lat = tmp_lat(i) + + ! Find indices for current tile lat and lon on scaling parameter grid + + i_ind = ceiling((this_lon - ll_lon)/dlon) + j_ind = ceiling((this_lat - ll_lat)/dlat) + + ! Check for no-data-values in observation and fit parameters + ! (any negative number could be no-data-value for observations) + + if ( sclprm_mean_obs(j_ind, i_ind)>0. .and. & + sclprm_mean_mod(j_ind, i_ind)>0. .and. & + sclprm_std_obs(j_ind, i_ind)>=0. .and. & + sclprm_std_mod(j_ind, i_ind)>=0. ) then + + ! Scale via standard normal deviates + + tmpreal = sclprm_std_mod(j_ind, i_ind)/sclprm_std_obs(j_ind, i_ind) + + tmp_obs(i) = sclprm_mean_mod(j_ind, i_ind) & + + tmpreal*(tmp_obs(i)-sclprm_mean_obs(j_ind, i_ind)) + + ! Check of tmp_obs is within range of model climatology + + if (tmp_obs(i)sclprm_max_mod(j_ind, i_ind)) then + + tmp_obs(i) = sclprm_max_mod(j_ind, i_ind) + + end if + + ! Scale observation error std + + tmp_std_obs(i) = tmpreal*tmp_std_obs(i) + + else + + tmp_obs(i) = this_obs_param%nodata + + end if + + end if + + end do + + deallocate(sclprm_mean_obs) + deallocate(sclprm_std_obs) + deallocate(sclprm_mean_mod) + deallocate(sclprm_std_mod) + deallocate(sclprm_min_mod) + deallocate(sclprm_max_mod) + + end subroutine scale_obs_sfmc_zscore + ! ******************************************************************************** subroutine scale_obs_Tb_zscore( N_catd, tile_coord, date_time, this_obs_param, & @@ -8083,7 +8839,7 @@ subroutine collect_obs( implicit none - character(*), intent(in) :: work_path + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id type(date_time_type), intent(in) :: date_time @@ -8691,4 +9447,3 @@ end program test ! ******* EOF ************************************************************* - diff --git a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 64b29673..2b18be6a 100644 --- a/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/src/Components/GEOSldas_GridComp/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -1124,12 +1124,12 @@ subroutine get_obs_pred( & select case (trim(obs_param(i)%varname)) - case ('sfmc') + case ('sfmc', 'sfds') get_sfmc_l = .true. get_sfmc_lH = .true. get_tsurf_l = .true. ! needed for model-based QC - + case ('rzmc') get_rzmc_l = .true. @@ -1641,7 +1641,7 @@ subroutine get_obs_pred( & select case (trim(obs_param(this_species)%varname)) - case ('sfmc') + case ('sfmc', 'sfds') tmp_data(1:N_tmp) = sfmc_lH( ind_tmp(1:N_tmp), n_e ) @@ -3595,15 +3595,16 @@ subroutine cat_enkf_increments( & select_update_type: select case (update_type) - case (1) select_update_type ! 1d soil moisture analysis; sfmc obs + case (1) select_update_type ! 1d soil moisture analysis; sfmc and/or sfds obs ! this 1d update requires that obs are on same tile space as model if (logit) write (logunit,*) 'get 1d soil moisture increments; sfmc obs' - N_select_varnames = 1 + N_select_varnames = 2 select_varnames(1) = 'sfmc' + select_varnames(2) = 'sfds' call get_select_species( & N_select_varnames, select_varnames(1:N_select_varnames), & @@ -3662,16 +3663,17 @@ subroutine cat_enkf_increments( & end do - case (2) select_update_type ! 3d soil moisture analysis; sfmc obs + case (2) select_update_type ! 3d soil moisture analysis; sfmc and/or sfds obs ! update each tile separately using all observations within ! the customized halo around each tile if (logit) write (logunit,*) 'get 3d soil moisture increments; sfmc obs' - N_select_varnames = 1 + N_select_varnames = 2 select_varnames(1) = 'sfmc' + select_varnames(2) = 'sfds' call get_select_species( & N_select_varnames, select_varnames(1:N_select_varnames), & diff --git a/src/Components/GEOSldas_GridComp/Shared/enkf_types.F90 b/src/Components/GEOSldas_GridComp/Shared/enkf_types.F90 index 46e5ef5b..c36f7127 100644 --- a/src/Components/GEOSldas_GridComp/Shared/enkf_types.F90 +++ b/src/Components/GEOSldas_GridComp/Shared/enkf_types.F90 @@ -53,7 +53,7 @@ module enkf_types ! added OminusA information for adaptive filtering - reichle, 1 Feb 2007 ! added "varname" field to "obs_param_type" - reichle 14 Jun 2011 ! major revisions to "obs_type" fields - reichle 16 Jun 2011 - ! added "units" field to "obs_param_tpye" - reichle 22 Nov 2011 + ! added "units" field to "obs_param_type" - reichle 22 Nov 2011 type :: obs_type @@ -96,6 +96,11 @@ module enkf_types integer :: species ! identifier for type of measurement integer :: orbit ! type of (half-)orbit + ! 0 = n/a [eg., in situ obs] + ! 1 = ascending + ! 2 = descending + ! 3 = ascending or descending + ! 4 = geostationary integer :: pol ! polarization ! 0 = n/a [eg., multi-pol. retrieval] diff --git a/src/Shared/.gitignore b/src/Shared/.gitignore index d5e49170..dba5156e 100644 --- a/src/Shared/.gitignore +++ b/src/Shared/.gitignore @@ -4,3 +4,6 @@ /@MAPL /MAPL /MAPL@ +/@NCEP_Shared +/NCEP_Shared +/NCEP_Shared@ diff --git a/src/Shared/CMakeLists.txt b/src/Shared/CMakeLists.txt index e7ace2f3..822331ee 100644 --- a/src/Shared/CMakeLists.txt +++ b/src/Shared/CMakeLists.txt @@ -1,3 +1,4 @@ esma_add_subdirectory (MAPL) esma_add_subdirectory (GMAO_Shared) esma_add_subdirectory (GEOS_Util) +esma_add_subdirectory (NCEP_Shared)