Skip to content

Commit

Permalink
Adding ASCAT SM reader and matlab util for clim_stats (#656)
Browse files Browse the repository at this point in the history
  • Loading branch information
gmao-rreichle committed Dec 24, 2023
2 parents 61d59b7 + 544a119 commit 8459e1b
Show file tree
Hide file tree
Showing 21 changed files with 2,096 additions and 549 deletions.
7 changes: 7 additions & 0 deletions components.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 8 additions & 0 deletions config/NCEP_Shared.sparse
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
!/*
/NCEP_sp
/NCEP_w3
/NCEP_bufr
/NCEP_bacio
/NCEP_sfcio
/NCEP_sigio
/CMakeLists.txt
113 changes: 112 additions & 1 deletion src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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

! --------------------------------------------------------------------

/
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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);
Expand Down Expand Up @@ -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

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Original file line number Diff line number Diff line change
@@ -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 =========================================================================
Loading

0 comments on commit 8459e1b

Please sign in to comment.