Adding ASCAT SM reader and matlab util for clim_stats

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

config/NCEP_Shared.sparse

@@ -0,0 +1,8 @@
+!/*
+/NCEP_sp
+/NCEP_w3
+/NCEP_bufr
+/NCEP_bacio
+/NCEP_sfcio
+/NCEP_sigio
+/CMakeLists.txt

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
+
 ! --------------------------------------------------------------------
 
 /

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
 

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

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 =========================================================================