diff --git a/CHANGELOG.md b/CHANGELOG.md index 09f0ead..77caf5d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- Revisions for handling of Nens and special nml and mwtrm path/files in coupled land-atm DAS. + ### Fixed - Fixed error from MAPL's ApplicationSupport.F90 to init UDUNITS diff --git a/GEOSldas_App/CMakeLists.txt b/GEOSldas_App/CMakeLists.txt index f13954e..0f4024a 100644 --- a/GEOSldas_App/CMakeLists.txt +++ b/GEOSldas_App/CMakeLists.txt @@ -39,7 +39,7 @@ configure_file(${file} ${file} @ONLY) install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/${file} DESTINATION bin) file(GLOB rc_files GEOSldas_*rc) -file(GLOB nml_files LDASsa_DEFAULT*nml) +file(GLOB nml_files LDASsa_DEFAULT*nml LandAtmDAS_nml/LDASsa_SPECIAL*nml ) install( FILES ${rc_files} ${nml_files} diff --git a/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensprop.nml b/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensprop.nml new file mode 100644 index 0000000..d76fef4 --- /dev/null +++ b/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensprop.nml @@ -0,0 +1,193 @@ +! +! namelist of EnKF inputs +! +! CUSTOM version for coupled land-atm DAS with soil moisture analysis based on SMAP Tb obs +! +! -------------------------------------------------------------------- + +&ens_prop_inputs + +progn_pert_dtstep = 10800 ! time step for generation of AR(1) prognostics perts [s] + +force_pert_dtstep = 10800 ! time step for generation of AR(1) forcing perts [s] + +! --------------------------------------------------------------------- +! +! forcing error (or forcing perturbation) parameters +! + +descr_force_pert%pcp = 'pcp' +descr_force_pert%sw = 'sw' +descr_force_pert%lw = 'lw' +descr_force_pert%tair = 'tair' +descr_force_pert%qair = 'qair' +descr_force_pert%wind = 'wind' + +! specify whether forcing perturbations are additive or multiplicative +! +! additive: typ = 0. +! multiplicative and lognormal: typ = 1. +! + +typ_force_pert%pcp = 1. +typ_force_pert%sw = 1. +typ_force_pert%lw = 0. + + +std_force_pert%pcp = 0.5 ! units if additive: [kg/m2/s], if multiplicative: [fraction] +std_force_pert%sw = 0.3 ! units if additive: [W/m^2] , if multiplicative: [fraction] +std_force_pert%lw = 20. ! units if additive: [W/m^2] , if multiplicative: [fraction] + +! read error std-dev from file? (if .false., default values above apply) + +stdfromfile_force_pert%pcp = .false. +stdfromfile_force_pert%sw = .false. +stdfromfile_force_pert%lw = .false. + +! specify file name (with full path) that contains std-dev values + +stdfilename_force_pert = '' + +! enforce zero (sample) mean across ensemble? + +zeromean_force_pert%pcp = .true. +zeromean_force_pert%sw = .true. +zeromean_force_pert%lw = .true. + +! allow perturbations to be computed on coarsened grid? + +coarsen_force_pert%pcp = .true. +coarsen_force_pert%sw = .true. +coarsen_force_pert%lw = .true. + +! max perturbation relative to standard normal +! (limits on range of random numbers: specify max absolute value +! allowed to be drawn from a standard normal distribution) + +std_normal_max_force_pert%pcp = 3. +std_normal_max_force_pert%sw = 3. +std_normal_max_force_pert%lw = 2.5 + +! spatial correlation of forcing perturbations + +xcorr_force_pert%pcp = 0.5 ! [deg] +xcorr_force_pert%sw = 0.5 ! [deg] +xcorr_force_pert%lw = 0.5 ! [deg] + +ycorr_force_pert%pcp = 0.5 ! [deg] +ycorr_force_pert%sw = 0.5 ! [deg] +ycorr_force_pert%lw = 0.5 ! [deg] + +! temporal correlation of forcing perturbations + +tcorr_force_pert%pcp = 86400. ! [s] +tcorr_force_pert%sw = 86400. ! [s] +tcorr_force_pert%lw = 86400. ! [s] + +! correlation coefficients -1 <= rho <= 1 +! +! specify only essential information, the other side of off-diagonals and +! the diagonal will be filled in later (subroutines read_ens_prop_inputs +! and get_force_pert_inputs) + +ccorr_force_pert%pcp%sw = -.8 +ccorr_force_pert%pcp%lw = .5 + +ccorr_force_pert%sw%lw = -.5 + + +! --------------------------------------------------------------------- +! +! model error (or progn_pert) parameters +! +! the mean is computed according to "typ" for unbiased perturbations +! and not specified here + +! string that describes the prognostics to be perturbed +! (see subroutine apply_progn_pert() for details) + +descr_progn_pert%catdef = 'catdef' +descr_progn_pert%srfexc = 'srfexc' + + +! specify whether model error is additive or multiplicative +! +! additive: typ = 0. +! multiplicative and lognormal: typ = 1. +! +! real numbers are used so that "assemble_state()" can +! be used to assemble the model error parameters + +typ_progn_pert%catdef = 0. +typ_progn_pert%srfexc = 0. + + +! The perturbation (or error) std-dev can be specified as a spatially constant +! (default) value. Alternatively, perturbation std-dev values can be read from +! a netcdf-4 input file (where they may be spatially constant or distributed). +! See subroutines get_progn_pert_param() and get_force_pert_param(). +! +! Turn off all perturbations by setting std-dev values to zero and +! "stdfromfile" to false. +! +! Default, spatially homogeneous perturbations std-dev +! (used unless std-devs are read from file, see below) + +std_progn_pert%catdef = 0.24 ! units if additive: [kg/m2/HOUR], if multiplicative: [fraction/HOUR] +std_progn_pert%srfexc = 0.16 ! units if additive: [kg/m2/HOUR], if multiplicative: [fraction/HOUR] + + +! read error std-dev from file? (if .false., default values above apply) + +stdfromfile_progn_pert%catdef = .false. +stdfromfile_progn_pert%srfexc = .false. + +! specify file name (with full path) that contains std-dev values + +stdfilename_progn_pert = '' + +! enforce zero (sample) mean across ensemble? + +zeromean_progn_pert%catdef = .true. +zeromean_progn_pert%srfexc = .true. + +! allow perturbations to be computed on coarsened grid? + +coarsen_progn_pert%catdef = .false. +coarsen_progn_pert%srfexc = .false. + +! max perturbation relative to standard normal +! (limits on range of random numbers: specify max absolute value +! allowed to be drawn from a standard normal distribution) + +std_normal_max_progn_pert%catdef = 2.5 +std_normal_max_progn_pert%srfexc = 2.5 + +! model error spatial correlation [deg] +! (x runs east-west, y runs north-south) + +xcorr_progn_pert%catdef = 0.3 ! [deg] +xcorr_progn_pert%srfexc = 0.3 ! [deg] + +ycorr_progn_pert%catdef = 0.3 ! [deg] +ycorr_progn_pert%srfexc = 0.3 ! [deg] + +! model error temporal correlation [s] + +tcorr_progn_pert%catdef = 10800. ! [s] +tcorr_progn_pert%srfexc = 10800. ! [s] + +! correlation coefficients -1 <= rho <= 1 +! +! specify only essential information, the other side of off-diagonals and +! the diagonal will be filled in later (subroutines read_ens_prop_inputs +! and get_force_pert_inputs) +! +! (the default input list below was put together with matlab +! script create_ccorr_cat_progn_default.m) + +ccorr_progn_pert%catdef%srfexc = 0.0 + +/ + +! =========================== EOF ======================================= diff --git a/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensupd.nml b/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensupd.nml new file mode 100644 index 0000000..e8ecb12 --- /dev/null +++ b/GEOSldas_App/LandAtmDAS_nml/LDASsa_SPECIAL_inputs_ensupd.nml @@ -0,0 +1,404 @@ +! +! namelist of EnKF inputs for land EnKF update +! +! CUSTOM version for coupled land-atm DAS with soil moisture analysis based on SMAP Tb obs +! +! ---------------------------------------------------------------------- + +&ens_upd_inputs + +! ---------------------------------------------------------------------- +! +! update type - for details see subroutine cat_enkf_update() +! (note: all 3d updates use compact support) +! +! local = "1d", regional = "3d" +! +! # = no longer supported +! +! update_type | analysis state vector | assimilated obs +! ------------------------------------------------------------------------------------------------- +! 0 | NO assimilation, NO bias correction | n/a +! # 1 | 1d soil moisture | sfmc +! # 2 | 3d soil moisture | sfmc +! 3 | 1d Tskin (assim incr NOT applied, use w/ bias corr) | Tskin +! 4 | 1d Tskin/ght1 (assim incr applied, use w/ or w/o bias corr) | Tskin +! 5 | 1d Tskin/ght1 (assim incr NOT applied, use w/ bias corr) | Tskin +! 6 | 1d soil moisture/Tskin/ght(1) | Tb +! 7 | 3d Tskin/ght1 update | Tskin +! 8 | 3d soil moisture/Tskin/ght(1) | Tb +! 9 | 1d Tskin/ght1 update | FT +! 10 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | Tb +! 11 | 1d snow analysis (Toure et al. 2018 empirical gain) | SCF +! 12 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | sfmc/sfds, Tb, SCF +! | & 1d snow analysis (Toure et al. 2018 empirical gain) | +! 13 | 3d soil moisture/Tskin/ght(1) excl. catdef unless PEATCLSM tile | sfmc/sfds, Tb + +update_type = 10 + +out_obslog = .true. +out_ObsFcstAna = .true. +out_smapL4SMaup = .false. + +! --------------------------------------------------------------------- +! +! Compact support parameters - for 3d updates +! +! All correlations vanish outside of an ellipse with semi-axes xcompact +! and ycompact + +xcompact = 1.25 ! [deg] longitude +ycompact = 1.25 ! [deg] latitude + +! --------------------------------------------------------------------- +! +! forecast error covariance inflaction factor +! +! - assigns more weight to observations in analysis by inflating forecast error covariance +! - works on std-dev, i.e., var_inflated = var * inflation_fac**2 +! - typical values: 1 <= inflation_fac <= 1.5 +! - to turn off, set to any negative real number + +fcsterr_inflation_fac = -9999. + +! --------------------------------------------------------------------- +! +! Definition of measurement species and parameters +! +! NOTE: When additional types of measurements are included here, +! at least the following parameters and subroutines must be adapted: +! +! - N_obs_species_nml in clsm_ensupd_glob_param.f90 +! [- read_ens_upd_inputs()] +! [- collect_obs()] +! - read_obs() +! - get_obs_pred() +! - cat_enkf_update() +! +! +! Definition of obs_param_nml fields (see also enkf_types.F90): +! +! %descr = description +! %species = identifier for type of measurement +! %orbit = type of (half-)orbit +! 0 = n/a [eg., in situ obs] +! 1 = ascending +! 2 = descending +! 3 = ascending or descending +! 4 = geostationary +! %pol = polarization +! 0 = n/a [eg., multi-pol. retrieval] +! 1 = horizontal +! 2 = vertical +! 3 = ... +! %N_ang = # satellite viewing angles in species (radiance obs only) +! %ang = vector of satellite viewing angles +! %freq = frequency [Hz] +! %FOV = field-of-view *radius*, see NOTES below +! (if FOV==0. equate obs footprint w/ tile) +! %FOV_units = field-of-view units ('km' or 'deg'), see NOTES below +! %assim = Should this obs type be assimilated (state update)? (logical) +! %scale = Should this obs be scaled? (logical) +! %getinnov = Should innov be computed for this obs type (logical) +! (innovations are always computed if assim==.true.) +! %RTM_ID = ID of radiative transfer model to use for Tb forward modeling +! (subroutine get_obs_pred()) +! 0 = none +! 1 = L-band tau-omega model as in De Lannoy et al. 2013 (doi:10.1175/JHM-D-12-092.1) (SMOS) +! 2 = same as 1 but without Pellarin atm corr (SMAP) +! 3 = same as 1 but with Mironov and SMAP L2_SM pol mixing (SMOS) +! 4 = same as 3 but without Pellarin atm corr (targeted for SMAP L4_SM Version 8) +! %bias_Npar = number of obs bias states tracked per day (integer) +! %bias_trel = e-folding time scale of obs bias memory [s] +! %bias_tcut = cutoff time for confident obs bias estimate [s] +! %nodata = no-data-value +! %varname = equivalent model variable name (for "Obs_pred") +! %units = units (eg., 'K' or 'm3/m3') +! %path = path to measurement files +! %name = name identifier for file containing measurements +! %maskpath = path to obs mask file +! %maskname = filename for obs mask +! %scalepath = path to file(s) with scaling parameters +! %scalename = filename for scaling parameters +! %flistpath = path to file with list of obs file names +! %flistname = name of file with list of obs file names +! %errstd = default obs error std +! %std_normal_max = maximum allowed perturbation (relative to N(0,1)) +! %zeromean = enforce zero mean across ensemble +! %coarsen_pert = generate obs perturbations on coarser grid (see pert_param_type%coarsen) +! %xcorr = correlation length (deg) in longitude direction +! %ycorr = correlation length (deg) in latitude direction +! +! For observation perturbations, always use: +! +! tcorr = 0. (never temporally correlated) +! typ = 0 (always additive) +! ccorr = 0. (never cross-correlated) +! +! (these are specified in get_obs_pert_inputs() and not here) +! +! +! NOTES: +! +! Field-of-view (FOV) can be specified in units of [km] or [deg] lat/lon. +! Note the special case of FOV=0. below. +! If FOV is specified in units of [km], the FOV in units of [deg] lat/lon that +! is used to compute observation predictions will depend on latitude. +! If FOV is specified in units of [deg] lat/lon, its value remains constant and +! is independent of latitude. +! The choice of units also determines the shape function that is used to +! compute the observation predictions. +! Units of [km] are meant for observations that are based on relatively +! coarse-scale measurements (such as microwave data). The resolution of such obs +! in units of [km] is approximately constant across the globe and independent +! of latitude. Observation predictions are computed by averaging tile-based +! model forecasts out to a distance of fac_search_FOV_km*FOV using a Gaussian kernel, +! where fac_search_FOV_km=2.0 as of 28 March 2015. +! Specifically, the normalized square distance is defined as +! +! ndist2 = dx^2/FOV_x^2 + dy^2/FOV_y^2 +! +! where FOV_x and dx are the meridional FOV and the meridional distance between the obs +! and the tile (in units of deg lat/lon), with FOV_x proportional to 1/cos(lat). +! FOV_y and dy are the corresponding zonal values. +! The weights are then proportional to +! +! exp( -0.5*ndist2 ) +! +! The averaging is therefore over an ellipse in lat/lon space, with weights +! decreasing away from the center of the observation. +! A 2.0*FOV averaging footprint encapsulates about 91% of the power. A 1.0*FOV +! averaging footprint would encapsulate about 47% of the power. These numbers +! are meant to be approximately consistent with FOV numbers for microwave radiometers +! (see 3 Dec 2014 email from Ed Kim reproduced below). +! Note that weights are further adjusted based on tile area. +! Units of [deg] lat/lon are meant for observations that are based on +! relatively high-resolution measurements (such as infrared data). Such +! observations are often available on a lat/lon grid that is much coarser than +! the footprint of the underlying observations. The assimilated data product +! therefore has a resolution that varies with latitude. Observation predictions are +! computed by averaging over a constant kernel out to a distance of FOV. +! The averaging is therefore over a circle in lat/lon space, with weights that do not +! depend on the distance from the center of the observation. +! (Note that weights are further adjusted based on tile area.) +! If FOV=0., observation predictions are computed by assigning the model forecast +! associated with the tile to which the observation is formally assigned. +! This is useful if the resolution of the assimilated observations is higher +! than that of the model tile space. This might be the case for snow-cover-fraction +! observations. FOV=0 can also be useful for tile-based synthetic observations. +! +! +! ------------------------------------------------------------------------ +! +! Date: Wed, 3 Dec 2014 11:21:30 -0600 +! From: +! To: , +! Subject: FW: [SMAP] antenna pattern question +! +! Hi Rolf & Gabrielle, +! +! First, a little terminology: the weighted integral is what Level 1 folks call +! "beam efficiency". So, apparently, Steven is assuming the "-3dB beam efficiency" +! is ~50%. The calculated [SMAP] beam efficiency within the -3dB contour is +! 53.40% (v-pol), 53.83% (h-pol). +! If you draw the -3dB contour on the Earth's surface, for h-pol, 53.83% of the energy +! comes from inside the contour, and 100-53.83 = 46.17% comes from outside the contour. +! The accuracy of the 1/10 and 1/100 digits is questionable, anyway. +! So, if you used 53% for v-pol and 54% for h-pol, you should be fine. +! I guess this means Steven was not far off, if he is using "50%." +! This -3dB beam efficiency means we have significant energy coming from outside +! the 3dB footprint, which is the footprint we use to come up with the "40 km" footprint +! size number. +! And, this is why many folks who use microwave instruments prefer to use a contour that +! encloses a higher % of the beam energy as a better measure of the footprint size. +! One such measure is the "main beam efficiency (MBE)." This beamwidth is usually taken +! to be 2.5 times the 3dB beamwidth. The corresponding footprint size is then +! 2.5x 40km = 100km. +! The last calculation put the MBE at 89.23 for V-pol and 89.33 for H-pol. +! So, for h-pol, 89.33% of the energy comes from inside a 100km footprint, +! and 100-89.33 = 10.67% from outside. +! - Ed +! +! ------------------------------------------------------------------------ +! +! IMPORTANT: The number of measurement species defined below must *match* +! global parameter "N_obs_species_nml" +! +! Multi-angular observations (eg., SMOS) are defined as a single +! species here (in the nml file) and are later split into +! multiple species, each having a unique incidence angle +! (see subroutine read_ens_upd_inputs()) +! +! +! ------------------------------------------------------------------------ + + +! -------------------------------------------------------------------- +! +! SMAP L1C_TB brightness temperature (36 km EASE grid) +! +! "A" = ascending (6pm *SMAP* overpass) +! "D" = descending (6am *SMAP* overpass) +! +! "Tbh" = h-pol Tb +! "Tbv" = v-pol Tb +! +! ------------------- +! +! 31 = SMAP_L1C_Tbh_A + +obs_param_nml(31)%descr = 'SMAP_L1C_Tbh_A' +obs_param_nml(31)%orbit = 1 +obs_param_nml(31)%pol = 1 +obs_param_nml(31)%N_ang = 1 +obs_param_nml(31)%ang(1) = 40. +obs_param_nml(31)%freq = 1.41e9 +obs_param_nml(31)%FOV = 20. +obs_param_nml(31)%FOV_units = 'km' +obs_param_nml(31)%assim = .true. +obs_param_nml(31)%scale = .true. +obs_param_nml(31)%getinnov = .true. +obs_param_nml(31)%RTM_ID = 4 +obs_param_nml(31)%bias_Npar = 0 +obs_param_nml(31)%bias_trel = 864000 +obs_param_nml(31)%bias_tcut = 432000 +obs_param_nml(31)%nodata = -9999. +obs_param_nml(31)%varname = 'Tb' +obs_param_nml(31)%units = 'K' +obs_param_nml(31)%path = '/discover/nobackup/projects/gmao/smap/SMAP_L4/SMAP/OPS/L1C_TB/' +obs_param_nml(31)%name = '' +obs_param_nml(31)%maskpath = '' +obs_param_nml(31)%maskname = '' +obs_param_nml(31)%scalepath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/scaling/SMAP_L1C_Tb_pentad/LADAS_v000/x6C_GLOBAL/' +obs_param_nml(31)%scalename = 'ScMO_SMAP__e24_zscore_stats_2015_p19_2020_p18_hscale_0.00_W_9p_Nmin_20' +obs_param_nml(31)%flistpath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/obs/SMAP/L1C_TB_flist/SPL4SM_OL8000/' +obs_param_nml(31)%flistname = 'SMAP_L1C_TB_A_list.txt' +obs_param_nml(31)%errstd = 4. +obs_param_nml(31)%std_normal_max = 2.5 +obs_param_nml(31)%zeromean = .true. +obs_param_nml(31)%coarsen_pert = .true. +obs_param_nml(31)%xcorr = 0.25 +obs_param_nml(31)%ycorr = 0.25 +obs_param_nml(31)%adapt = 0 + +! ------------------- +! +! 32 = SMAP_L1C_Tbh_D + +obs_param_nml(32)%descr = 'SMAP_L1C_Tbh_D' +obs_param_nml(32)%orbit = 2 +obs_param_nml(32)%pol = 1 +obs_param_nml(32)%N_ang = 1 +obs_param_nml(32)%ang(1) = 40. +obs_param_nml(32)%freq = 1.41e9 +obs_param_nml(32)%FOV = 20. +obs_param_nml(32)%FOV_units = 'km' +obs_param_nml(32)%assim = .true. +obs_param_nml(32)%scale = .true. +obs_param_nml(32)%getinnov = .true. +obs_param_nml(32)%RTM_ID = 4 +obs_param_nml(32)%bias_Npar = 0 +obs_param_nml(32)%bias_trel = 864000 +obs_param_nml(32)%bias_tcut = 432000 +obs_param_nml(32)%nodata = -9999. +obs_param_nml(32)%varname = 'Tb' +obs_param_nml(32)%units = 'K' +obs_param_nml(32)%path = '/discover/nobackup/projects/gmao/smap/SMAP_L4/SMAP/OPS/L1C_TB/' +obs_param_nml(32)%name = '' +obs_param_nml(32)%maskpath = '' +obs_param_nml(32)%maskname = '' +obs_param_nml(32)%scalepath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/scaling/SMAP_L1C_Tb_pentad/LADAS_v000/x6C_GLOBAL/' +obs_param_nml(32)%scalename = 'ScMO_SMAP__e24_zscore_stats_2015_p19_2020_p18_hscale_0.00_W_9p_Nmin_20' +obs_param_nml(32)%flistpath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/obs/SMAP/L1C_TB_flist/SPL4SM_OL8000/' +obs_param_nml(32)%flistname = 'SMAP_L1C_TB_D_list.txt' +obs_param_nml(32)%errstd = 4. +obs_param_nml(32)%std_normal_max = 2.5 +obs_param_nml(32)%zeromean = .true. +obs_param_nml(32)%coarsen_pert = .true. +obs_param_nml(32)%xcorr = 0.25 +obs_param_nml(32)%ycorr = 0.25 +obs_param_nml(32)%adapt = 0 + +! ------------------- +! +! 33 = SMAP_L1C_Tbv_A + +obs_param_nml(33)%descr = 'SMAP_L1C_Tbv_A' +obs_param_nml(33)%orbit = 1 +obs_param_nml(33)%pol = 2 +obs_param_nml(33)%N_ang = 1 +obs_param_nml(33)%ang(1) = 40. +obs_param_nml(33)%freq = 1.41e9 +obs_param_nml(33)%FOV = 20. +obs_param_nml(33)%FOV_units = 'km' +obs_param_nml(33)%assim = .true. +obs_param_nml(33)%scale = .true. +obs_param_nml(33)%getinnov = .true. +obs_param_nml(33)%RTM_ID = 4 +obs_param_nml(33)%bias_Npar = 0 +obs_param_nml(33)%bias_trel = 864000 +obs_param_nml(33)%bias_tcut = 432000 +obs_param_nml(33)%nodata = -9999. +obs_param_nml(33)%varname = 'Tb' +obs_param_nml(33)%units = 'K' +obs_param_nml(33)%path = '/discover/nobackup/projects/gmao/smap/SMAP_L4/SMAP/OPS/L1C_TB/' +obs_param_nml(33)%name = '' +obs_param_nml(33)%maskpath = '' +obs_param_nml(33)%maskname = '' +obs_param_nml(33)%scalepath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/scaling/SMAP_L1C_Tb_pentad/LADAS_v000/x6C_GLOBAL/' +obs_param_nml(33)%scalename = 'ScMO_SMAP__e24_zscore_stats_2015_p19_2020_p18_hscale_0.00_W_9p_Nmin_20' +obs_param_nml(33)%flistpath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/obs/SMAP/L1C_TB_flist/SPL4SM_OL8000/' +obs_param_nml(33)%flistname = 'SMAP_L1C_TB_A_list.txt' +obs_param_nml(33)%errstd = 4. +obs_param_nml(33)%std_normal_max = 2.5 +obs_param_nml(33)%zeromean = .true. +obs_param_nml(33)%coarsen_pert = .true. +obs_param_nml(33)%xcorr = 0.25 +obs_param_nml(33)%ycorr = 0.25 +obs_param_nml(33)%adapt = 0 + +! ------------------- +! +! 34 = SMAP_L1C_Tbv_D + +obs_param_nml(34)%descr = 'SMAP_L1C_Tbv_D' +obs_param_nml(34)%orbit = 2 +obs_param_nml(34)%pol = 2 +obs_param_nml(34)%N_ang = 1 +obs_param_nml(34)%ang(1) = 40. +obs_param_nml(34)%freq = 1.41e9 +obs_param_nml(34)%FOV = 20. +obs_param_nml(34)%FOV_units = 'km' +obs_param_nml(34)%assim = .true. +obs_param_nml(34)%scale = .true. +obs_param_nml(34)%getinnov = .true. +obs_param_nml(34)%RTM_ID = 4 +obs_param_nml(34)%bias_Npar = 0 +obs_param_nml(34)%bias_trel = 864000 +obs_param_nml(34)%bias_tcut = 432000 +obs_param_nml(34)%nodata = -9999. +obs_param_nml(34)%varname = 'Tb' +obs_param_nml(34)%units = 'K' +obs_param_nml(34)%path = '/discover/nobackup/projects/gmao/smap/SMAP_L4/SMAP/OPS/L1C_TB/' +obs_param_nml(34)%name = '' +obs_param_nml(34)%maskpath = '' +obs_param_nml(34)%maskname = '' +obs_param_nml(34)%scalepath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/scaling/SMAP_L1C_Tb_pentad/LADAS_v000/x6C_GLOBAL/' +obs_param_nml(34)%scalename = 'ScMO_SMAP__e24_zscore_stats_2015_p19_2020_p18_hscale_0.00_W_9p_Nmin_20' +obs_param_nml(34)%flistpath = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/obs/SMAP/L1C_TB_flist/SPL4SM_OL8000/' +obs_param_nml(34)%flistname = 'SMAP_L1C_TB_D_list.txt' +obs_param_nml(34)%errstd = 4. +obs_param_nml(34)%std_normal_max = 2.5 +obs_param_nml(34)%zeromean = .true. +obs_param_nml(34)%coarsen_pert = .true. +obs_param_nml(34)%xcorr = 0.25 +obs_param_nml(34)%ycorr = 0.25 +obs_param_nml(34)%adapt = 0 + + +! -------------------------------------------------------------------- + + +/ + +! =========================== EOF ======================================= diff --git a/GEOSldas_App/ldas_setup b/GEOSldas_App/ldas_setup index d1f332c..76ba840 100755 --- a/GEOSldas_App/ldas_setup +++ b/GEOSldas_App/ldas_setup @@ -111,9 +111,8 @@ class LDASsetup: self.agcm_res = cmdLineArgs['agcm_res'] self.bcs_version = cmdLineArgs['bcs_version'] self.rstloc = cmdLineArgs['rstloc'] - self.varwindow = cmdLineArgs['varwindow'] - self.LDAS_nml_path = cmdLineArgs['LDAS_nml_path'] - self.mwrtm_path = cmdLineArgs['mwrtm_path'] + self.varwindow = cmdLineArgs['varwindow'] + self.nens = cmdLineArgs['nens'] # obsolete command line args self.runmodel = cmdLineArgs['runmodel'] @@ -190,81 +189,79 @@ class LDASsetup: ### check if ldas is coupled to adas; if so, set/overwrite input parameters accordingly if self.ladas_cpl > 0 : - # make sure all necessary command line arguments were supplied - assert self.nymdb is not None, "Error. Must have command line arg nymdb for coupled land-atm DAS.\n" - assert self.nhmsb is not None, "Error. Must have command line arg nhmsb for coupled land-atm DAS.\n" - assert self.agcm_res is not None, "Error. Must have command line arg agcm_res for coupled land-atm DAS.\n" - assert self.bcs_version is not None, "Error. Must have command line arg bcs_version for coupled land-atm DAS.\n" - assert self.rstloc is not None, "Error. Must have command line arg rstloc for coupled land-atm DAS.\n" - assert self.varwindow is not None, "Error. Must have command line arg varwindow for coupled land-atm DAS.\n" - assert self.LDAS_nml_path is not None, "Error. Must have command line arg LDAS_nml_path for coupled land-atm DAS.\n" - assert self.mwrtm_path is not None, "Error. Must have command line arg mwrtm_path for coupled land-atm DAS.\n" + # make sure all necessary command line arguments were supplied + assert self.nymdb is not None, "Error. Must have command line arg nymdb for coupled land-atm DAS.\n" + assert self.nhmsb is not None, "Error. Must have command line arg nhmsb for coupled land-atm DAS.\n" + assert self.agcm_res is not None, "Error. Must have command line arg agcm_res for coupled land-atm DAS.\n" + assert self.bcs_version is not None, "Error. Must have command line arg bcs_version for coupled land-atm DAS.\n" + assert self.rstloc is not None, "Error. Must have command line arg rstloc for coupled land-atm DAS.\n" + assert self.varwindow is not None, "Error. Must have command line arg varwindow for coupled land-atm DAS.\n" + assert self.nens is not None, "Error. Must have command line arg nens for coupled land-atmensDAS.\n" - self.rqdExeInp['BEG_DATE'] = f"{self.nymdb} {self.nhmsb}" - rstloc_ = self.rstloc.rstrip('/') # remove trailing '/' - assert os.path.isdir(rstloc_) # make sure rstloc_ is a valid directory - self.rstloc = os.path.abspath(rstloc_) - self.rqdExeInp['RESTART_PATH'] = os.path.dirname( self.rstloc) - self.rqdExeInp['RESTART_ID'] = os.path.basename(self.rstloc) - - self.adas_expdir = os.path.dirname( self.exphome) - self.rqdExeInp['ADAS_EXPDIR'] = self.adas_expdir - self.adas_expid = os.path.basename(self.adas_expdir) - self.rqdExeInp['MET_TAG'] = self.adas_expid + '__bkg' - if self.ladas_cpl == 1 : - # ldas coupled with determistic component of adas - self.rqdExeInp['EXP_ID'] = self.adas_expid + '_LDAS' - self.rqdExeInp['MET_PATH'] = self.adas_expdir + '/recycle/holdpredout' - self.rqdExeInp['ENSEMBLE_FORCING'] = 'NO' - self.nens = 24 # for now, hardwire nens for atm det - elif self.ladas_cpl == 2 : - # ldas coupled with ensemble component of adas - self.rqdExeInp['EXP_ID'] = self.adas_expid + '_LDAS4ens' - self.rqdExeInp['MET_PATH'] = self.adas_expdir + '/atmens/mem' - self.rqdExeInp['ENSEMBLE_FORCING'] = 'YES' - self.nens = 32 # for now, hardwire nens for atm ens - else : - exit("Error. Unknown value of self.ladas_cpl.\n") - self.rqdExeInp['NUM_LDAS_ENSEMBLE'] = self.nens - self.first_ens_id = 1 - self.rqdExeInp['FIRST_ENS_ID'] = self.first_ens_id - self.rqdExeInp['EXP_DOMAIN'] = 'CF'+ self.agcm_res +'x6C_GLOBAL' + self.rqdExeInp['BEG_DATE'] = f"{self.nymdb} {self.nhmsb}" + rstloc_ = self.rstloc.rstrip('/') # remove trailing '/' + assert os.path.isdir(rstloc_) # make sure rstloc_ is a valid directory + self.rstloc = os.path.abspath(rstloc_) + self.rqdExeInp['RESTART_PATH'] = os.path.dirname( self.rstloc) + self.rqdExeInp['RESTART_ID'] = os.path.basename(self.rstloc) + self.adas_expdir = os.path.dirname( self.exphome) + self.rqdExeInp['ADAS_EXPDIR'] = self.adas_expdir + self.adas_expid = os.path.basename(self.adas_expdir) + self.rqdExeInp['MET_TAG'] = self.adas_expid + '__bkg' + + if self.ladas_cpl == 1 : + # ldas coupled with determistic component of ADAS + self.rqdExeInp['EXP_ID'] = self.adas_expid + '_LDAS' + self.rqdExeInp['MET_PATH'] = self.adas_expdir + '/recycle/holdpredout' + self.rqdExeInp['ENSEMBLE_FORCING'] = 'NO' + elif self.ladas_cpl == 2 : + # ldas coupled with ensemble component of ADAS + self.rqdExeInp['EXP_ID'] = self.adas_expid + '_LDAS4ens' + self.rqdExeInp['MET_PATH'] = self.adas_expdir + '/atmens/mem' + self.rqdExeInp['ENSEMBLE_FORCING'] = 'YES' + else : + exit("Error. Unknown value of self.ladas_cpl.\n") - # when coupled to ADAS, "BCS_PATH" EXCLUDE bcs version info - # hard-wired BCS_PATH for now - self.rqdExeInp['BCS_PATH'] = "/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles" - self.rqdExeInp['BCS_PATH'] = self.rqdExeInp['BCS_PATH'].rstrip('/') + '/' + self.bcs_version - if self.bcs_version == "Icarus-NLv3" : - self.rqdExeInp['BCS_PATH'] = self.rqdExeInp['BCS_PATH'] + '_new_layout' - self.rqdExeInp['BCS_RESOLUTION'] = 'CF'+ self.agcm_res +'x6C_CF' + self.agcm_res +'x6C' - self.rqdExeInp['RESTART_DOMAIN'] = 'CF'+ self.agcm_res +'x6C_GLOBAL' - - # the following are not in default rqdExeInp list; hardwire for now - self.rqdExeInp['LAND_ASSIM'] = "YES" - self.rqdExeInp['MET_HINTERP'] = 0 - self.landassim_dt = 10800 # seconds - # make sure ADAS analysis window [minutes] is multiple of LANDASSIM_DT [seconds] - if int(self.varwindow) % (self.landassim_dt/60) == 0 : - self.rqdExeInp['LANDASSIM_DT'] = self.landassim_dt - else : - exit("Error. LANDASSIM_DT is inconsistent with ADAS analysis window.\n") - self.rqdExeInp['LANDASSIM_T0'] = "013000" # HHMMSS - jsgmt1 = "00000000" - jsgmt2 = hours_to_hhmmss(int(self.varwindow)/60) # convert minutes to HHMMSS - self.rqdExeInp['JOB_SGMT'] = f"{jsgmt1} {jsgmt2}" - self.rqdExeInp['NUM_SGMT'] = 1 - self.rqdExeInp['FORCE_DTSTEP'] = 3600 - - # determine END_DATE = BEG_DATE + TIME_STEP_OF_ADAS_CYCLE - _beg_date = datetime.strptime( self.rqdExeInp['BEG_DATE'], "%Y%m%d %H%M%S") - _hours = int(self.rqdExeInp['JOB_SGMT'][ 9:11]) - _end_date = _beg_date + timedelta(hours=int(self.varwindow)/60) - - self.rqdExeInp['END_DATE'] = _end_date.strftime("%Y%m%d %H%M%S") - self.rqdExeInp['NML_INPUT_PATH'] = self.LDAS_nml_path - self.rqdExeInp['MWRTM_PATH'] = self.mwrtm_path + self.rqdExeInp['NUM_LDAS_ENSEMBLE'] = self.nens # fvsetup finds Nens by counting restart files + self.first_ens_id = 1 # match ADAS convention + self.rqdExeInp['FIRST_ENS_ID'] = self.first_ens_id + + self.agcm_res = 'CF' + self.agcm_res # change format to "CFnnnn" + self.rqdExeInp['EXP_DOMAIN'] = self.agcm_res +'x6C_GLOBAL' + + # when coupled to ADAS, "BCS_PATH" EXCLUDE bcs version info + # hard-wired BCS_PATH for now + self.rqdExeInp['BCS_PATH'] = "/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles" + self.rqdExeInp['BCS_PATH'] = self.rqdExeInp['BCS_PATH'].rstrip('/') + '/' + self.bcs_version + if self.bcs_version == "Icarus-NLv3" : + self.rqdExeInp['BCS_PATH'] = self.rqdExeInp['BCS_PATH'] + '_new_layout' + self.rqdExeInp['BCS_RESOLUTION'] = self.agcm_res +'x6C_' + self.agcm_res +'x6C' + self.rqdExeInp['RESTART_DOMAIN'] = self.agcm_res +'x6C_GLOBAL' + + # the following are not in default rqdExeInp list; hardwire for now + self.rqdExeInp['MWRTM_PATH'] = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/RTM_params/RTMParam_SMAP_L4SM_v006/' + self.rqdExeInp['LAND_ASSIM'] = "YES" + self.rqdExeInp['MET_HINTERP'] = 0 + self.landassim_dt = 10800 # seconds + # make sure ADAS analysis window [minutes] is multiple of LANDASSIM_DT [seconds] + if int(self.varwindow) % (self.landassim_dt/60) == 0 : + self.rqdExeInp['LANDASSIM_DT'] = self.landassim_dt + else : + exit("Error. LANDASSIM_DT is inconsistent with ADAS analysis window.\n") + self.rqdExeInp['LANDASSIM_T0'] = "013000" # HHMMSS + jsgmt1 = "00000000" + jsgmt2 = hours_to_hhmmss(int(self.varwindow)/60) # convert minutes to HHMMSS + self.rqdExeInp['JOB_SGMT'] = f"{jsgmt1} {jsgmt2}" + self.rqdExeInp['NUM_SGMT'] = 1 + self.rqdExeInp['FORCE_DTSTEP'] = 3600 + + # determine END_DATE = BEG_DATE + TIME_STEP_OF_ADAS_CYCLE + _beg_date = datetime.strptime( self.rqdExeInp['BEG_DATE'], "%Y%m%d %H%M%S") + _hours = int(self.rqdExeInp['JOB_SGMT'][ 9:11]) + _end_date = _beg_date + timedelta(hours=int(self.varwindow)/60) + self.rqdExeInp['END_DATE'] = _end_date.strftime("%Y%m%d %H%M%S") - # end if self.ladas_cpl > 0 + # end if self.ladas_cpl > 0 ----------------------------------------------------------------------------------------- for key in rqdExeInpKeys : assert key in self.rqdExeInp,' "%s" is required in the input file %s' % (key,self.exeinpfile) @@ -632,22 +629,22 @@ class LDASsetup: # exefyl # ------ - self.bindir = os.path.dirname(os.path.realpath(__file__)) - self.blddir = self.bindir.rsplit('/',1)[0] - exefyl = '/bin/GEOSldas.x' - tmp_execfyl= self.blddir+exefyl + self.bindir = os.path.dirname(os.path.realpath(__file__)) + self.blddir = self.bindir.rsplit('/',1)[0] + exefyl = '/bin/GEOSldas.x' + tmp_execfyl = self.blddir + exefyl assert os.path.isfile(tmp_execfyl),\ 'Executable [%s] does not exist!' % tmp_execfyl - self.expdir = self.exphome + '/' + self.rqdExeInp['EXP_ID'] - self.rundir = self.expdir + '/run' - self.inpdir = self.expdir + '/input' - self.outdir = self.expdir + '/output' - self.scratchdir = self.expdir + '/scratch' - self.blddirLn = self.expdir + '/build' - self.out_path = self.outdir+'/'+self.rqdExeInp['EXP_DOMAIN'] - self.bcsdir = self.outdir+'/'+self.rqdExeInp['EXP_DOMAIN']+'/rc_out/' - self.rstdir = self.outdir+'/'+self.rqdExeInp['EXP_DOMAIN']+'/rs/' - self.exefyl = self.blddirLn+exefyl + self.expdir = self.exphome + '/' + self.rqdExeInp['EXP_ID'] + self.rundir = self.expdir + '/run' + self.inpdir = self.expdir + '/input' + self.outdir = self.expdir + '/output' + self.scratchdir = self.expdir + '/scratch' + self.blddirLn = self.expdir + '/build' + self.out_path = self.outdir + '/'+self.rqdExeInp['EXP_DOMAIN'] + self.bcsdir = self.outdir + '/'+self.rqdExeInp['EXP_DOMAIN']+'/rc_out/' + self.rstdir = self.outdir + '/'+self.rqdExeInp['EXP_DOMAIN']+'/rs/' + self.exefyl = self.blddirLn + exefyl my_ntasks_per_node = int(self.rqdRmInp['ntasks-per-node']) @@ -1238,17 +1235,22 @@ class LDASsetup: for nmlfile in default_nml: shortfile=self.rundir+'/'+nmlfile.split('/')[-1] shutil.copy2(nmlfile, shortfile) - # special nml + # special nml special_nml=[] - if 'NML_INPUT_PATH' in self.rqdExeInp : - special_nml = glob.glob(self.rqdExeInp['NML_INPUT_PATH']+'/LDASsa_SPECIAL_inputs_*.nml') - if self.ladas_cpl > 0: - special_nml = glob.glob(self.rqdExeInp['NML_INPUT_PATH']+'/'+ self.rqdExeInp['EXP_DOMAIN'] + '/LDASsa_SPECIAL_inputs_*.nml') - + if self.ladas_cpl > 0: + special_nml= glob.glob(etcdir+'/LDASsa_SPECIAL_inputs_*.nml') + else : + if 'NML_INPUT_PATH' in self.rqdExeInp : + special_nml = glob.glob(self.rqdExeInp['NML_INPUT_PATH']+'/LDASsa_SPECIAL_inputs_*.nml') + for nmlfile in special_nml: - shortfile=nmlfile.split('/')[-1] - shutil.copy2(nmlfile, self.rundir+'/'+shortfile) + shortfile=self.rundir+'/'+nmlfile.split('/')[-1] + shutil.copy2(nmlfile, shortfile) + if self.ladas_cpl > 0: + # edit resolution info in ensupd nml file + sp.run(['sed', '-i', 's//'+self.agcm_res+'/g', self.rundir+'/LDASsa_SPECIAL_inputs_ensupd.nml']) + # get optimzed NX and IMS optimized_distribution_file = tempfile.NamedTemporaryFile(delete=False) print ("Optimizing... decomposition of processes.... \n") @@ -1960,7 +1962,7 @@ def parseCmdLine(): ) p_setup.add_argument( '--agcm_res', - help='AGCM resolution associated with boundary conditions: xxxx (4 digits). Required when exeinp is dummy', + help='AGCM resolution associated with boundary conditions: nnnn (4 digits). Required when exeinp is dummy', type=str ) p_setup.add_argument( @@ -1979,16 +1981,10 @@ def parseCmdLine(): type=str ) p_setup.add_argument( - '--LDAS_nml_path', - help='path to LDAS (special) nml input files. Required when exeinp is dummy', + '--nens', + help='number of ensemble members. Required when exeinp is dummy', type=str - ) - p_setup.add_argument( - '--mwrtm_path', - help='path to LDAS L-band microwave radiative transfer model (mwrtm) parameters. Required when exeinp is dummy', - type=str - ) - + ) # obsolete command line args p_setup.add_argument(