diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/saber_blocks.yaml deleted file mode 100644 index 0a04b634e..000000000 --- a/parm/soca/berror/saber_blocks.yaml +++ /dev/null @@ -1,101 +0,0 @@ -covariance model: hybrid -components: -- covariance: - # This will setup B = K D C_v C_h C_h^{T} C_v^{T} D K^{T} - covariance model: SABER - saber central block: - saber block name: EXPLICIT_DIFFUSION - active variables: [tocn, socn, ssh, cicen] - geometry: - mom6_input_nml: mom_input.nml - fields metadata: ./fields_metadata.yaml - group mapping: - - name: ocean - variables: [tocn, socn, ssh] - - name: ice - variables: [cicen] - read: - groups: - - name: ocean - horizontal: - filename: hz_ocean.nc - vertical: - filename: vt_ocean.nc - - name: ice - horizontal: - filename: hz_ice.nc - - linear variable change: - input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] - output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] - linear variable changes: - - - linear variable change name: BkgErrFILT - ocean_depth_min: 500 # [m] - rescale_bkgerr: 1.0 - efold_z: 1500.0 # [m] - - - linear variable change name: BkgErrSOCA - read_from_file: 3 - basename: ./static_ens/ - ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' - ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' - date: '{{ATM_WINDOW_MIDDLE}}' - t_min: 0.1 - t_max: 5.0 - s_min: 0.1 - s_max: 1.0 - ssh_min: 0.0 # std ssh=0 => ssh balance applied as - ssh_max: 1.0 # strong constraint - cicen_min: 0.1 - cicen_max: 0.5 - hicen_min: 0.0 - hicen_max: 0.0 - standard deviation: true - - - linear variable change name: BalanceSOCA - ksshts: - nlayers: 10 - - weight: - value: 0.25 -#- covariance: -# covariance model: ensemble -# members from template: -# template: -# read_from_file: 1 -# date: '{{ATM_WINDOW_MIDDLE}}' -# basename: ./static_ens/ -# ocn_filename: 'ocn.pert.steric.%mem%.nc' -# ice_filename: 'ice.%mem%.nc' -# state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon] -# pattern: '%mem%' -# nmembers: ${ENS_SIZE} -# localization: -# localization method: SABER -# saber central block: -# saber block name: EXPLICIT_DIFFUSION -# active variables: [tocn, socn, ssh] -# geometry: -# mom6_input_nml: mom_input.nml -# fields metadata: ./fields_metadata.yaml -# group mapping: -# - name: ocean -# variables: [tocn, socn, ssh, uocn, vocn] -# - name: ice -# variables: [cicen, hicen, hsnon] -# read: -# groups: -# - name: ocean -# horizontal: -# filename: hz_ocean.nc -# - name: ice -# horizontal: -# filename: hz_ice.nc -# -# weight: -# read_from_file: 3 -# basename: ./ -# ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' -# ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' -# date: '{{ATM_WINDOW_MIDDLE}}' diff --git a/parm/soca/berror/soca_diagb.yaml.j2 b/parm/soca/berror/soca_diagb.yaml.j2 new file mode 100644 index 000000000..c57a4d82f --- /dev/null +++ b/parm/soca/berror/soca_diagb.yaml.j2 @@ -0,0 +1,59 @@ +geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + +date: '{{ ATM_WINDOW_MIDDLE }}' + +background: + date: '{{ ATM_WINDOW_MIDDLE }}' + basename: ./bkg/ + ocn_filename: '{{ RUN }}.ocean.t{{ gcyc }}z.inst.f009.nc' + ice_filename: '{{ RUN }}.agg_ice.t{{ gcyc }}z.inst.f009.nc' + read_from_file: 1 + +background error: + datadir: ./ + date: '{{ ATM_WINDOW_MIDDLE }}' + exp: bkgerr_stddev + type: incr + +variables: + name: [tocn, socn, uocn, vocn, hocn, ssh, cicen, hicen, hsnon, mom6_mld] + +rescale: 2.0 # rescales the filtered std. dev. by "rescale" +min sst: 0.0 # Added to sst bkg. err. +max ssh: 0.0 # Limits the amplitude of the unbalanced bkg err +number of halo points: 4 +number of neighbors: 16 +simple smoothing: + horizontal iterations: 2 + vertical iterations: 1 + +# TODO(G): Start using when the normalization is optional +#diffusion: +# saber block name: EXPLICIT_DIFFUSION +# active variables: [tocn, socn, ssh, cicen, hicen, hsnon] +# geometry: +# mom6_input_nml: mom_input.nml +# fields metadata: ./fields_metadata.yaml +# group mapping: +# - name: ocean +# variables: +# - tocn +# - socn +# - ssh +# - name: ice +# variables: +# - cicen +# - hicen +# - hsnon +# read: +# groups: +# - name: ocean +# horizontal: +# filename: hz_ocean.nc +# vertical: +# filename: vt_ocean.nc +# - name: ice +# horizontal: +# filename: hz_ice.nc diff --git a/parm/soca/berror/soca_hybrid_bmat.yaml b/parm/soca/berror/soca_hybrid_bmat.yaml new file mode 100644 index 000000000..6e3d35351 --- /dev/null +++ b/parm/soca/berror/soca_hybrid_bmat.yaml @@ -0,0 +1,92 @@ +covariance model: hybrid +components: +- covariance: + covariance model: SABER + saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh, cicen] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: + - tocn + - socn + - ssh + - name: ice + variables: + - cicen + read: + groups: + - name: ocean + horizontal: + filename: hz_ocean.nc + vertical: + filename: vt_ocean.nc + - name: ice + horizontal: + filename: hz_ice.nc + + saber outer blocks: + - saber block name: StdDev + read: + model file: + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./ + ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + read_from_file: 3 + + linear variable change: + input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + linear variable changes: + - linear variable change name: BalanceSOCA + + weight: + value: 1.00 + +- covariance: + covariance model: ensemble + members from template: + template: + read_from_file: 1 + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./static_ens/ + ocn_filename: 'ocn.pert.steric.%mem%.nc' + ice_filename: 'ice.%mem%.nc' + state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon] + pattern: '%mem%' + nmembers: ${ENS_SIZE} + localization: + localization method: SABER + saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: [tocn, socn, ssh, uocn, vocn] + - name: ice + variables: [cicen, hicen, hsnon] + read: + groups: + - name: ocean + multivariate strategy: duplicated + horizontal: + filename: hz_ocean.nc + vertical: + strategy: duplicated + - name: ice + horizontal: + filename: hz_ice.nc + + weight: + read_from_file: 3 + basename: ./ + ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc' + date: '{{ATM_WINDOW_MIDDLE}}' diff --git a/parm/soca/berror/soca_parameters_diffusion_vt.yaml b/parm/soca/berror/soca_parameters_diffusion_vt.yaml index 8597bfe15..98370be7d 100644 --- a/parm/soca/berror/soca_parameters_diffusion_vt.yaml +++ b/parm/soca/berror/soca_parameters_diffusion_vt.yaml @@ -25,6 +25,8 @@ background error: - name: vt_ocean vertical: as gaussian: true - fixed value: 5.0 + from file: + filename: vt_scales.nc + variable name: vt write: filename: vt_ocean.nc diff --git a/parm/soca/berror/soca_static_bmat.yaml b/parm/soca/berror/soca_static_bmat.yaml new file mode 100644 index 000000000..36e9274a7 --- /dev/null +++ b/parm/soca/berror/soca_static_bmat.yaml @@ -0,0 +1,42 @@ +covariance model: SABER +saber central block: + saber block name: EXPLICIT_DIFFUSION + active variables: [tocn, socn, ssh, cicen] + geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + group mapping: + - name: ocean + variables: + - tocn + - socn + - ssh + - name: ice + variables: + - cicen + read: + groups: + - name: ocean + horizontal: + filename: hz_ocean.nc + vertical: + filename: vt_ocean.nc + - name: ice + horizontal: + filename: hz_ice.nc + +saber outer blocks: +- saber block name: StdDev + read: + model file: + date: '{{ATM_WINDOW_MIDDLE}}' + basename: ./ + ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_MIDDLE}}.nc' + read_from_file: 3 + +linear variable change: + input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh] + linear variable changes: + - linear variable change name: BalanceSOCA diff --git a/parm/soca/berror/soca_vtscales.yaml.j2 b/parm/soca/berror/soca_vtscales.yaml.j2 new file mode 100644 index 000000000..bad7d623f --- /dev/null +++ b/parm/soca/berror/soca_vtscales.yaml.j2 @@ -0,0 +1,13 @@ +gridspec_filename: soca_gridspec.nc +restart_filename: ./INPUT/MOM.res.nc +mld_filename: 'ocn.bkgerr_stddev.incr.{{ ATM_WINDOW_MIDDLE }}.nc' +output_filename: ./vt_scales.nc +output_variable_vt: vt +output_variable_hz: hz + +VT_MIN: 5 +VT_MAX: 15 + +HZ_ROSSBY_MULT: 1.0 +HZ_MAX: 200e3 +HZ_MIN_GRID_MULT: 2.0 diff --git a/parm/soca/fields_metadata.yaml b/parm/soca/fields_metadata.yaml new file mode 100644 index 000000000..d239c4000 --- /dev/null +++ b/parm/soca/fields_metadata.yaml @@ -0,0 +1,180 @@ +# -------------------------------------------------------------------------------------------------- +# Field metadata for SOCA. Each field can contain the following information: +# +# name: Internal name used by soca code and config files +# grid: "h", "u", or "v" (Default: h) +# masked: use land mask if true (Default: true) +# levels: "1" or "full_ocn" (Default: 1) +# getval_name: variable name expected by GetValues (Default: ) +# getval_name_surface: GetValues variable name for 2D surface of a 3D field (Default: ) +# io_file: The restart file domain "ocn", "sfc", or "ice" (Default: ) +# io_name: The variable name used in the restart IO (Default: ) +# -------------------------------------------------------------------------------------------------- + + +# -------------------------------------------------------------------------------------------------- +# Ocean state variables +# -------------------------------------------------------------------------------------------------- +- name: tocn + levels: full_ocn + getval name: sea_water_potential_temperature + getval name surface: sea_surface_temperature + io file: ocn + io name: Temp + fill value: 0.0 + +- name: socn + levels: full_ocn + getval name: sea_water_salinity + getval name surface: sea_surface_salinity + io file: ocn + io name: Salt + property: positive_definite + fill value: 0.0 + +- name: uocn + grid: u + levels: full_ocn + getval name: eastward_sea_water_velocity + getval name surface: surface_eastward_sea_water_velocity + io file: ocn + io name: u + fill value: 0.0 + +- name: vocn + grid: v + levels: full_ocn + getval name: northward_sea_water_velocity + getval name surface: surface_northward_sea_water_velocity + io file: ocn + io name: v + fill value: 0.0 + +- name: hocn + levels: full_ocn + getval name: sea_water_cell_thickness + io file: ocn + io name: h + fill value: 0.001 + vert interp: false + +- name: ssh + getval name: sea_surface_height_above_geoid + io file: ocn + io name: ave_ssh + fill value: 0.0 + +- name: mom6_mld + io file: ocn + io name: MLD + fill value: 0.0 + +# -------------------------------------------------------------------------------------------------- +# ice state variables +# -------------------------------------------------------------------------------------------------- +- name: hicen + getval name: sea_ice_category_thickness + io file: ice + io name: hicen + property: positive_definite + fill value: 0.0 + +- name: cicen + getval name: sea_ice_category_area_fraction + getval name surface: sea_ice_area_fraction # note: not accurate, should be "sum" not "surface" + io file: ice + io name: aicen + fill value: 0.0 + +- name: hsnon + getval name: sea_ice_category_snow_thickness + io file: ice + io name: hsnon + property: positive_definite + fill value: 0.0 + +# -------------------------------------------------------------------------------------------------- +# wave state variables +# -------------------------------------------------------------------------------------------------- +- name: swh + getval name: sea_surface_wave_significant_height + io file: wav + io name: hs + property: positive_definite + +# -------------------------------------------------------------------------------------------------- +# sea surface variables +# -------------------------------------------------------------------------------------------------- +- name: sw + masked: false + getval name: net_downwelling_shortwave_radiation + io file: sfc + io name: sw_rad + +- name: lw + masked: false + getval name: net_downwelling_longwave_radiation + io file: sfc + io name: lw_rad + +- name: lhf + masked: false + getval name: upward_latent_heat_flux_in_air + io file: sfc + io name: latent_heat + +- name: shf + masked: false + getval name: upward_sensible_heat_flux_in_air + io file: sfc + io name: sens_heat + +- name: us + masked: false + getval name: friction_velocity_over_water + io file: sfc + io name: fric_vel + + +# -------------------------------------------------------------------------------------------------- +# BGC +# -------------------------------------------------------------------------------------------------- +- name: chl + levels: full_ocn + getval name: mass_concentration_of_chlorophyll_in_sea_water + getval name surface: sea_surface_chlorophyll + io file: ocn + io name: chl + property: positive_definite + +- name: biop + levels: full_ocn + getval name: molar_concentration_of_biomass_in_sea_water_in_p_units + getval name surface: sea_surface_biomass_in_p_units + io file: ocn + io name: biomass_p + property: positive_definite + +# -------------------------------------------------------------------------------------------------- +# built-in derived variables that you don't need to worry about +# -------------------------------------------------------------------------------------------------- +- name: distance_from_coast + masked: false + +- name: layer_depth + levels: full_ocn + vert interp: false + +- name: mesoscale_representation_error + +- name: mld + +- name: sea_floor_depth_below_sea_surface + +- name: sea_area_fraction + masked: false + +- name: surface_temperature_where_sea + +- name: sea_water_depth + levels: full_ocn diff --git a/parm/soca/obs/config/adt_rads_all.yaml b/parm/soca/obs/config/adt_rads_all.yaml index fe41a50da..74469902b 100644 --- a/parm/soca/obs/config/adt_rads_all.yaml +++ b/parm/soca/obs/config/adt_rads_all.yaml @@ -41,7 +41,7 @@ obs filters: variables: [GeoVaLs/mesoscale_representation_error, ObsError/absoluteDynamicTopography] coefs: [0.1, - 0.5] + 1.0] - filter: Domain Check where: - variable: { name: GeoVaLs/sea_ice_area_fraction} diff --git a/parm/soca/obs/config/sst_abi_g16_l3c.yaml b/parm/soca/obs/config/sst_abi_g16_l3c.yaml index 689207b14..b9d999843 100644 --- a/parm/soca/obs/config/sst_abi_g16_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g16_l3c.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_abi_g17_l3c.yaml b/parm/soca/obs/config/sst_abi_g17_l3c.yaml index b2c3cb692..a7ada5774 100644 --- a/parm/soca/obs/config/sst_abi_g17_l3c.yaml +++ b/parm/soca/obs/config/sst_abi_g17_l3c.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml index 5597d53c5..ed0efd721 100644 --- a/parm/soca/obs/config/sst_ahi_h08_l3c.yaml +++ b/parm/soca/obs/config/sst_ahi_h08_l3c.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml index a3fbfe6ca..73c048891 100644 --- a/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_ma_l3u.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml index 9c9f9cfbf..ccac3cc4e 100644 --- a/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mb_l3u.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml index 57e759017..9c34bcd56 100644 --- a/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml +++ b/parm/soca/obs/config/sst_avhrr_mc_l3u.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml index 36f155576..04021235f 100644 --- a/parm/soca/obs/config/sst_viirs_n20_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_n20_l3u.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml index c1647ada8..6df13028a 100644 --- a/parm/soca/obs/config/sst_viirs_npp_l3u.yaml +++ b/parm/soca/obs/config/sst_viirs_npp_l3u.yaml @@ -40,7 +40,17 @@ obs filters: where: - variable: {name: GeoVaLs/sea_surface_temperature} minvalue: -1.0 -- filter: Domain Check - where: - - variable: {name: GeoVaLs/distance_from_coast} - minvalue: 100e3 +#- filter: Domain Check +# where: +# - variable: {name: GeoVaLs/distance_from_coast} +# minvalue: 100e3 +- filter: Perform Action + action: + name: assign error + error function: + name: ObsFunction/LinearCombination + options: + variables: + - ObsError/seaSurfaceTemperature + coefs: + - 0.05 diff --git a/parm/soca/obsprep/obsprep_config.yaml b/parm/soca/obsprep/obsprep_config.yaml index 3e4d55c6a..cb7e4e9dc 100644 --- a/parm/soca/obsprep/obsprep_config.yaml +++ b/parm/soca/obsprep/obsprep_config.yaml @@ -22,8 +22,8 @@ observations: type: nc dmpdir regex: 'rads_adt_*.nc' window: - back: 8 # look back 8 six-hourly obs dumps - forward: 1 # look forward 1 six-hourly bin + back: 0 + forward: 0 error ratio: 0.4 # meters per day # Ice concentration diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index 1b7f7c215..9ae2b3646 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -56,54 +56,75 @@ else fi fi -################################################################################ -# Write ensemble weights for the hybrid envar -$APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err -fi ################################################################################ -# Ensemble perturbations for the EnVAR and diagonal of static B -# Static B: -# - Compute moments of original ensemble with the balanced ssh removed -# from the statistics -# -# Ensemble of perturbations: -# - apply the linear steric height balance to each members, using the -# deterministic for trajectory. -# - add the unblanced ssh to the steric ssh field -# Diagnostics: -# - variance explained by the steric heigh -# - moments +# Ensemble processing +if [ "${NMEM_ENS}" -gt 3 ]; then + ################################################################################ + # Write ensemble weights for the hybrid envar + $APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi -# Process static ensemble -clean_yaml soca_clim_ens_moments.yaml -$APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err + ################################################################################ + # Ensemble perturbations for the EnVAR + # Ensemble diagnostics: + # - Compute moments of original ensemble with the balanced ssh removed + # from the statistics + # + # Ensemble of perturbations: + # - apply the linear steric height balance to each members, using the + # deterministic for trajectory. + # - add the unblanced ssh to the steric ssh field + # Diagnostics: + # - variance explained by the steric heigh + # - moments + clean_yaml soca_clim_ens_moments.yaml + $APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi fi ################################################################################ -# Set decorrelation scales for the static B -$APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml +# Compute the parametric diag of B +$APRUN_OCNANAL $JEDI_BIN/gdas_soca_diagb.x soca_diagb.yaml export err=$?; err_chk if [ $err -gt 0 ]; then exit $err fi ################################################################################ -# Initialize diffusion blocks -clean_yaml soca_parameters_diffusion_hz.yaml -$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml -export err=$?; err_chk -if [ $err -gt 0 ]; then - exit $err +# Horizontal diffusion +if [ ! -f "ocn.cor_rh.incr.0001-01-01T00:00:00Z.nc" ]; then + # Set decorrelation scales for the static B + $APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi + # Initialize the horizontal diffusion block and normalize + clean_yaml soca_parameters_diffusion_hz.yaml + $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml + export err=$?; err_chk + if [ $err -gt 0 ]; then + exit $err + fi +else + echo "Diffusion weights already calculated, skipping the init of the horizontal diffusion" fi +################################################################################ +# Initialize the vertical diffusion block +# Compute the MLD dependent vertical scales +# TODO(G): Probably not the right way to execute this script +python ${HOMEgfs}/sorc/gdas.cd/sorc/soca/tools/calc_scales.py soca_vtscales.yaml clean_yaml soca_parameters_diffusion_vt.yaml + +# Initialize the vertical diffusion block and normalize $APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_vt.yaml export err=$?; err_chk if [ $err -gt 0 ]; then diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index e527509fd..54cce2363 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -53,6 +53,7 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): bdatedt = datetime.strptime(cdate, '%Y%m%d%H') - timedelta(hours=3) bdate = datetime.strftime(bdatedt, '%Y-%m-%dT%H:00:00Z') mdate = datetime.strftime(datetime.strptime(cdate, '%Y%m%d%H'), '%Y-%m-%dT%H:00:00Z') +nmem_ens = int(os.getenv('NMEM_ENS')) logger.info(f"---------------- Copy from RUNDIR to COMOUT") @@ -65,12 +66,13 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): domains = ['ocn', 'ice'] for domain in domains: # Copy of the diagonal of the background error for the cycle - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.bkgerr_stddev.incr.{bdate}.nc'), + post_file_list.append([os.path.join(anl_dir, f'{domain}.bkgerr_stddev.incr.{mdate}.nc'), os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.bkgerr_stddev.nc')]) # Copy the recentering error - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), - os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) + if nmem_ens > 2: + post_file_list.append([os.path.join(anl_dir, 'static_ens', f'{domain}.ssh_recentering_error.incr.{bdate}.nc'), + os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}.recentering_error.nc')]) # Copy the ice and ocean increments post_file_list.append([os.path.join(anl_dir, 'Data', f'{domain}.3dvarfgat_pseudo.incr.{mdate}.nc'), @@ -81,9 +83,10 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.{domain}ana.nc')]) # Copy of the ssh diagnostics -for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: - post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), - os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) +if nmem_ens > 2: + for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']: + post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'), + os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')]) # Copy DA grid (computed for the start of the window) post_file_list.append([os.path.join(anl_dir, 'soca_gridspec.nc'), diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 171ea8154..d80bfc834 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -10,7 +10,8 @@ from soca import bkg_utils from datetime import datetime, timedelta import pytz -from wxflow import (Logger, Template, TemplateConstants, YAMLFile, FileHandler) +from wxflow import (Logger, Template, TemplateConstants, + YAMLFile, FileHandler, AttrDict, parse_j2yaml) logger = Logger() @@ -52,6 +53,7 @@ def find_clim_ens(input_date): """ Find the clim. ens. that is the closest to the DA window """ + logger.info(f"$$$$$$$$$$$$$$$$$$$$$$$ {os.getenv('SOCA_INPUT_FIX_DIR')}") ens_clim_dir = os.path.join(os.getenv('SOCA_INPUT_FIX_DIR'), 'bkgerr', 'ens') dirs = glob.glob(os.path.join(ens_clim_dir, '*')) @@ -66,9 +68,10 @@ def find_clim_ens(input_date): comin_obs = os.getenv('COMIN_OBS') anl_dir = os.getenv('DATA') staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR') +nmem_ens = 0 +nmem_ens = int(os.getenv('NMEM_ENS')) if os.getenv('DOHYBVAR') == "YES": dohybvar = True - nmem_ens = int(os.getenv('NMEM_ENS')) else: dohybvar = False @@ -101,7 +104,9 @@ def find_clim_ens(input_date): envconfig = {'window_begin': f"{window_begin.strftime('%Y-%m-%dT%H:%M:%SZ')}", 'ATM_WINDOW_BEGIN': window_begin_iso, 'ATM_WINDOW_MIDDLE': window_middle_iso, - 'ATM_WINDOW_LENGTH': f"PT{os.getenv('assim_freq')}H"} + 'ATM_WINDOW_LENGTH': f"PT{os.getenv('assim_freq')}H", + 'gcyc': gcyc, + 'RUN': RUN} stage_cfg = YAMLFile(path=os.path.join(gdas_home, 'parm', 'templates', 'stage.yaml')) stage_cfg = Template.substitute_structure(stage_cfg, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) stage_cfg = Template.substitute_structure(stage_cfg, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) @@ -137,7 +142,6 @@ def find_clim_ens(input_date): # stage ensemble members if dohybvar: logger.info("---------------- Stage ensemble members") - nmem_ens = int(os.getenv('NMEM_ENS')) ens_member_list = [] for mem in range(1, nmem_ens+1): for domain in ['ocean', 'ice']: @@ -159,17 +163,20 @@ def find_clim_ens(input_date): cice_fname = os.path.realpath(os.path.join(static_ens, "ice."+str(mem)+".nc")) bkg_utils.cice_hist2fms(cice_fname, cice_fname) else: - logger.info("---------------- Stage offline ensemble members") - ens_member_list = [] - clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) - nmem_ens = len(glob.glob(os.path.join(clim_ens_dir, 'ocean.*.nc'))) - for domain in ['ocean', 'ice']: - for mem in range(1, nmem_ens+1): - fname = domain+"."+str(mem)+".nc" - fname_in = os.path.join(clim_ens_dir, fname) - fname_out = os.path.join(static_ens, fname) - ens_member_list.append([fname_in, fname_out]) - FileHandler({'copy': ens_member_list}).sync() + if nmem_ens >= 3: + logger.info("---------------- Stage offline ensemble members") + ens_member_list = [] + clim_ens_dir = find_clim_ens(pytz.utc.localize(window_begin, is_dst=None)) + list_of_members = glob.glob(os.path.join(clim_ens_dir, 'ocean.*.nc')) + nmem_ens = min(nmem_ens, len(list_of_members)) + for domain in ['ocean', 'ice']: + for mem in range(1, nmem_ens+1): + fname = domain+"."+str(mem)+".nc" + fname_in = os.path.join(clim_ens_dir, fname) + fname_out = os.path.join(static_ens, fname) + ens_member_list.append([fname_in, fname_out]) + FileHandler({'copy': ens_member_list}).sync() + os.environ['ENS_SIZE'] = str(nmem_ens) ################################################################################ @@ -189,6 +196,11 @@ def find_clim_ens(input_date): # generate the YAML file for the post processing of the clim. ens. B berror_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') +logger.info(f"---------------- generate soca_diagb.yaml") +conf = parse_j2yaml(path=os.path.join(berror_yaml_dir, 'soca_diagb.yaml.j2'), + data=envconfig) +conf.save(os.path.join(anl_dir, 'soca_diagb.yaml')) + logger.info(f"---------------- generate soca_ensb.yaml") berr_yaml = os.path.join(anl_dir, 'soca_ensb.yaml') berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_ensb.yaml') @@ -230,6 +242,11 @@ def find_clim_ens(input_date): config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) config.save(diffu_hz_yaml) +logger.info(f"---------------- generate soca_vtscales.yaml") +conf = parse_j2yaml(path=os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_vtscales.yaml.j2'), + data=envconfig) +conf.save(os.path.join(anl_dir, 'soca_vtscales.yaml')) + logger.info(f"---------------- generate soca_parameters_diffusion_vt.yaml") diffu_vt_yaml = os.path.join(anl_dir, 'soca_parameters_diffusion_vt.yaml') diffu_vt_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') @@ -254,13 +271,13 @@ def find_clim_ens(input_date): yaml_name='bkg_list.yaml') os.environ['BKG_LIST'] = 'bkg_list.yaml' -# select the SABER BLOCKS to use -if 'SABER_BLOCKS_YAML' in os.environ and os.environ['SABER_BLOCKS_YAML']: - saber_blocks_yaml = os.getenv('SABER_BLOCKS_YAML') - logger.info(f"using non-default SABER blocks yaml: {saber_blocks_yaml}") +# select the B-matrix to use +if nmem_ens > 3: + # Use a hybrid static-ensemble B + os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_hybrid_bmat.yaml') else: - logger.info(f"using default SABER blocks yaml") - os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'saber_blocks.yaml') + # Use a static B + os.environ['SABER_BLOCKS_YAML'] = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_static_bmat.yaml') # substitute templated variables in the var config logger.info(f"{config}") diff --git a/scripts/exgdas_global_marine_analysis_run.sh b/scripts/exgdas_global_marine_analysis_run.sh index 3e938df37..7ee52cd5f 100755 --- a/scripts/exgdas_global_marine_analysis_run.sh +++ b/scripts/exgdas_global_marine_analysis_run.sh @@ -38,7 +38,7 @@ function clean_yaml() } ################################################################################ -# run 3DVAR FGAT +# run the variational application cp var.yaml var_original.yaml clean_yaml var.yaml $APRUN_OCNANAL $JEDI_BIN/soca_var.x var.yaml diff --git a/test/soca/gw/static.sh b/test/soca/gw/static.sh index 05407ea92..30e6b25c9 100755 --- a/test/soca/gw/static.sh +++ b/test/soca/gw/static.sh @@ -14,8 +14,8 @@ lowres=${project_source_dir}/sorc/soca/test/Data cp -L ${lowres}/workdir/{diag_table,field_table} ${soca_static} cp -L ${project_source_dir}/test/soca/fix/MOM_input ${soca_static} -cp -L ${lowres}/{fields_metadata.yml,godas_sst_bgerr.nc,rossrad.dat} ${soca_static} -mv ${soca_static}/fields_metadata.yml ${soca_static}/fields_metadata.yaml +cp -L ${lowres}/rossrad.dat ${soca_static} +cp -L ${project_source_dir}/parm/soca/fields_metadata.yaml ${soca_static} cp -L ${project_source_dir}/sorc/soca/test/testinput/obsop_name_map.yml ${soca_static}/obsop_name_map.yaml cp -L ${lowres}/72x35x25/input.nml ${soca_static}/inputnml cp -L ${lowres}/72x35x25/INPUT/{hycom1_25.nc,ocean_mosaic.nc,grid_spec.nc,layer_coord25.nc,ocean_hgrid.nc,ocean_topog.nc} ${soca_static}/INPUT diff --git a/utils/soca/CMakeLists.txt b/utils/soca/CMakeLists.txt index 39da43a66..62996f2af 100644 --- a/utils/soca/CMakeLists.txt +++ b/utils/soca/CMakeLists.txt @@ -23,3 +23,9 @@ ecbuild_add_executable( TARGET gdassoca_obsstats.x SOURCES gdassoca_obsstats.cc ) target_compile_features( gdassoca_obsstats.x PUBLIC cxx_std_17) target_link_libraries( gdassoca_obsstats.x PUBLIC NetCDF::NetCDF_CXX oops ioda) + +# Diag B +ecbuild_add_executable( TARGET gdas_soca_diagb.x + SOURCES gdas_soca_diagb.cc ) +target_compile_features( gdas_soca_diagb.x PUBLIC cxx_std_17) +target_link_libraries( gdas_soca_diagb.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) diff --git a/utils/soca/gdas_soca_diagb.cc b/utils/soca/gdas_soca_diagb.cc new file mode 100644 index 000000000..7273f2e75 --- /dev/null +++ b/utils/soca/gdas_soca_diagb.cc @@ -0,0 +1,8 @@ +#include "gdas_soca_diagb.h" +#include "oops/runs/Run.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + gdasapp::SocaDiagB diagb; + return run.execute(diagb); +} diff --git a/utils/soca/gdas_soca_diagb.h b/utils/soca/gdas_soca_diagb.h new file mode 100644 index 000000000..62cff84f7 --- /dev/null +++ b/utils/soca/gdas_soca_diagb.h @@ -0,0 +1,465 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/mesh.h" +#include "atlas/mesh/actions/BuildEdges.h" +#include "atlas/mesh/actions/BuildHalo.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/util/Earth.h" +#include "atlas/util/Geometry.h" +#include "atlas/util/Point.h" + +#include "oops/base/FieldSet3D.h" +#include "oops/base/GeometryData.h" +#include "oops/mpi/mpi.h" +#include "oops/runs/Application.h" +#include "oops/util/DateTime.h" +#include "oops/util/Duration.h" +#include "oops/util/FieldSetHelpers.h" +#include "oops/util/FieldSetOperations.h" +#include "oops/util/Logger.h" + +#include "soca/ExplicitDiffusion/ExplicitDiffusion.h" +#include "soca/ExplicitDiffusion/ExplicitDiffusionParameters.h" +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" +#include "soca/State/State.h" + +namespace gdasapp { + /** + * SocaDiagB Class Implementation + * + * Implements variance partitioning within the GDAS for the ocean and sea-ice. It is used as a proxy for the + * diagonal of B. + * This class is designed to partition the variance of ocean and sea-ice fields by analyzing the variance within + * predefined geographical bins. + * This approach allows for an ensemble-free estimate of a flow-dependent background error. + */ + + // ----------------------------------------------------------------------------- + // Since we can't have useful data members in SocaDiagB because execute is "const" + // and the constructor does not know about the configuration, + // we're going around the issue with the use of a simple data structure + struct SocaDiagBConfig { + util::DateTime cycleDate; + oops::Variables socaVars; + double sshMax; // poor man's alternative to limiting the unbalanced ssh variance + double depthMin; // set the bkg error to 0 for depth < depthMin + bool diffusion; // apply explicit diffusion to the std. dev. fields + bool simpleSmoothing; // Simple spatial averaging + int niterHoriz; // number of iterations for the horizontal averaging + int niterVert; // number of iterations for the vertical averaging + double rescale; // inflation/deflation of the variance after filtering + }; + + // ----------------------------------------------------------------------------- + class SocaDiagB : public oops::Application { + public: + // ----------------------------------------------------------------------------- + explicit SocaDiagB(const eckit::mpi::Comm & comm = oops::mpi::world()) + : Application(comm) { + } + + // ----------------------------------------------------------------------------- + /** + * Sets up and returns a SocaDiagBConfig object configured according to the provided + * eckit::Configuration + */ + SocaDiagBConfig setup(const eckit::Configuration & fullConfig) const { + SocaDiagBConfig diagBConfig; + + /// Get the date + // ------------- + oops::Log::info() << "====================== date" << std::endl; + std::string strdt; + fullConfig.get("date", strdt); + diagBConfig.cycleDate = util::DateTime(strdt); + + /// Get the list of variables + // -------------------------- + oops::Log::info() << "====================== variables" << std::endl; + oops::Variables vars(fullConfig, "variables.name"); + diagBConfig.socaVars = vars; + + // Get the max std dev for ssh + diagBConfig.sshMax = 0.0; + fullConfig.get("max ssh", diagBConfig.sshMax); + + // Get the minimum depth for which to partition the 3D field's std. dev. + diagBConfig.depthMin = 0.0; + fullConfig.get("min depth", diagBConfig.depthMin); + + // Explicit diffusion + diagBConfig.diffusion = false; + if (fullConfig.has("diffusion")) { + diagBConfig.diffusion = true; + } + + // Simple smoothing parameters + diagBConfig.simpleSmoothing = false; + if (fullConfig.has("simple smoothing")) { + diagBConfig.simpleSmoothing = true; + fullConfig.get("simple smoothing.horizontal iterations", diagBConfig.niterHoriz); + fullConfig.get("simple smoothing.vertical iterations", diagBConfig.niterVert); + } + + // Variance rescaling + fullConfig.get("rescale", diagBConfig.rescale); + + return diagBConfig; + } + + // ----------------------------------------------------------------------------- + static const std::string classname() {return "gdasapp::SocaDiagB";} + + // ----------------------------------------------------------------------------- + /** + * Variance partitioning at jnode and level + */ + void stdDevFilt(const int jnode, + const int level, + const int nbz, + const double depthMin, + const std::vector neighbors, + const int nzMld, + const atlas::array::ArrayView h, + const atlas::array::ArrayView bkg, + const atlas::array::ArrayView bathy, + atlas::array::ArrayView& stdDevBkg, + bool doBathy = true, + int minn = 6) const { + // Early exit if too shallow + if ( doBathy & bathy(jnode, 0) < depthMin ) { + stdDevBkg(jnode, level) = 0.0; + return; + } + + int nbh = neighbors.size(); + std::vector local; + for (int nn = 0; nn < neighbors.size(); ++nn) { + int levelMin = std::max(0, level - nbz); + int levelMax = level + nbz; + // Do weird things withing the MLD + if (level < nzMld) { + // If in the MLD, compute the std. dev. throughout the MLD + levelMin = 0; + levelMax = 1; // nzMld; Projecting the surface variance down the water column for now + } + // 2D case + if (bkg.shape(1) == 1) { + levelMin = 0; + levelMax = 0; // nzMld; Projecting the surface variance down the water column for now + } + for (int ll = levelMin; ll <= levelMax; ++ll) { + if ( abs(h(neighbors[nn], ll)) <= 0.1 ) { + continue; + } + local.push_back(bkg(neighbors[nn], ll)); + } + } + if (local.size() >= minn) { + // Mean + double mean = std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + + // Standard deviation + double stdDev(0.0); + for (int nn = 0; nn < nbh; ++nn) { + stdDev += std::pow(local[nn] - mean, 2.0); + } + if (stdDev > 0.0 || local.size() > 2) { + stdDevBkg(jnode, level) = std::sqrt(stdDev / (local.size() - 1)); + } + } + } + + // ----------------------------------------------------------------------------- + /** + * Implementation of the virtual execute method from the Application parent class + */ + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + /// Initialize the paramters for D + SocaDiagBConfig configD = setup(fullConfig); + + /// Setup the soca geometry + oops::Log::info() << "====================== geometry" << std::endl; + const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); + const soca::Geometry geom(geomConfig, this->getComm()); + + /// Read the background + // -------------------- + oops::Log::info() << "====================== read bkg" << std::endl; + soca::State xb(geom, configD.socaVars, configD.cycleDate); + const eckit::LocalConfiguration bkgConfig(fullConfig, "background"); + xb.read(bkgConfig); + atlas::FieldSet xbFs; + xb.toFieldSet(xbFs); + oops::Log::info() << "Background:" << std::endl; + oops::Log::info() << xb << std::endl; + + /// Create the mesh connectivity (Copy/paste of Francois's stuff) + // -------------------------------------------------------------- + // Build edges, then connections between nodes and edges + oops::Log::info() << "====================== build mesh connectivity" << std::endl; + int nbHalos(2); + fullConfig.get("number of halo points", nbHalos); + int nbNeighbors(4); + fullConfig.get("number of neighbors", nbNeighbors); + atlas::functionspace::NodeColumns nodeColumns = geom.functionSpace(); + atlas::Mesh mesh = nodeColumns.mesh(); + atlas::mesh::actions::build_edges(mesh); + atlas::mesh::actions::build_node_to_edge_connectivity(mesh); + atlas::mesh::actions::build_halo(mesh, nbHalos); + const auto & node2edge = mesh.nodes().edge_connectivity(); + const auto & edge2node = mesh.edges().node_connectivity(); + + // Lambda function to get the neighbors of a node + // (Copy/paste from Francois's un-merged oops branch) + const auto get_neighbors_of_node = [&](const int node) { + std::vector neighbors{}; + neighbors.reserve(nbNeighbors); + if (node >= mesh.nodes().size()) { + return neighbors; + } + // Use node2edge and edge2node connectivities to find neighbors of node + const int nb_edges = node2edge.cols(node); + for (size_t ie = 0; ie < nb_edges; ++ie) { + const int edge = node2edge(node, ie); + ASSERT(edge2node.cols(edge) == 2); // sanity check, maybe not in production/release build + const int node0 = edge2node(edge, 0); + const int node1 = edge2node(edge, 1); + ASSERT(node == node0 || node == node1); // sanity check edge contains initial node + if (node != node0) { + neighbors.push_back(node0); + } else { + neighbors.push_back(node1); + } + } + return neighbors; + }; + + /// Compute local std. dev. as a proxy of the bkg error + // ---------------------------------------------------- + oops::Log::info() << "====================== start variance partitioning" << std::endl; + // Identify halo nodes + const auto ghostView = + atlas::array::make_view(geom.functionSpace().ghost()); + + // Create the background error fieldset + oops::Log::info() << "====================== allocate std. dev. field set" << std::endl; + soca::Increment bkgErr(geom, configD.socaVars, configD.cycleDate); + bkgErr.zero(); + atlas::FieldSet bkgErrFs; + bkgErr.toFieldSet(bkgErrFs); + + // Get the layer thicknesses and convert to layer depth + oops::Log::info() << "====================== calculate layer depth" << std::endl; + auto viewHocn = atlas::array::make_view(xbFs["hocn"]); + atlas::array::ArrayT depth(viewHocn.shape(0), viewHocn.shape(1)); + auto viewDepth = atlas::array::make_view(depth); + for (atlas::idx_t jnode = 0; jnode < depth.shape(0); ++jnode) { + viewDepth(jnode, 0) = 0.5 * viewHocn(jnode, 0); + for (atlas::idx_t level = 1; level < depth.shape(1); ++level) { + viewDepth(jnode, level) = viewDepth(jnode, level-1) + + 0.5 * (viewHocn(jnode, level-1) + viewHocn(jnode, level)); + } + } + + // Compute the bathymetry + oops::Log::info() << "====================== calculate bathymetry" << std::endl; + atlas::array::ArrayT bathy(viewHocn.shape(0), 1); + auto viewBathy = atlas::array::make_view(bathy); + for (atlas::idx_t jnode = 0; jnode < viewHocn.shape(0); ++jnode) { + viewBathy(jnode, 0) = 0.0; + for (atlas::idx_t level = 0; level < viewHocn.shape(1); ++level) { + viewBathy(jnode, 0) += viewHocn(jnode, level); + } + } + + // Get the mixed layer depth + auto mld = atlas::array::make_view(xbFs["mom6_mld"]); + atlas::array::ArrayT mldindex(mld.shape(0), mld.shape(1)); + auto viewMldindex = atlas::array::make_view(mldindex); + + for (atlas::idx_t jnode = 0; jnode < viewHocn.shape(0); ++jnode) { + std::vector testMld; + for (atlas::idx_t level = 0; level < viewDepth.shape(1); ++level) { + testMld.push_back(std::abs(viewDepth(jnode, level) - mld(jnode, 0))); + } + auto minVal = std::min_element(testMld.begin(), testMld.end()); + + viewMldindex(jnode, 0) = std::distance(testMld.begin(), minVal); + } + + // Update the layer thickness halo + nodeColumns.haloExchange(xbFs["hocn"]); + + // Loop through variables + for (auto & var : configD.socaVars.variables()) { + // Update the halo + nodeColumns.haloExchange(xbFs[var]); + + // Skip the layer thickness variable + if (var == "hocn") { + continue; + } + oops::Log::info() << "====================== std dev for " << var << std::endl; + auto bkg = atlas::array::make_view(xbFs[var]); + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + + // Loop through nodes + for (atlas::idx_t jnode = 0; jnode < xbFs[var].shape(0); ++jnode) { + // Early exit if thickness is 0 or on a ghost cell + if (ghostView(jnode) > 0 || abs(viewHocn(jnode, 0)) <= 0.1) { + continue; + } + + // get indices of the cell's neighbors + auto neighbors = get_neighbors_of_node(jnode); + + // 2D case + if (xbFs[var].shape(1) == 1) { + // Std. dev. of the partition + stdDevFilt(jnode, 0, 0, + configD.depthMin, neighbors, 0, + viewHocn, bkg, viewBathy, stdDevBkg, false, 4); + if (var == "ssh") { + // TODO(G): Extract the unbalanced ssh variance, in the mean time, do this: + stdDevBkg(jnode, 0) = std::min(configD.sshMax, stdDevBkg(jnode, 0)); + } + } else { + // 3D case + int nbz = 1; // Number of closest point in the vertical, above and below + int nzMld = std::max(10, viewMldindex(jnode, 0)); + for (atlas::idx_t level = 0; level < xbFs[var].shape(1) - nbz; ++level) { + // Std. dev. of the partition + stdDevFilt(jnode, level, nbz, + configD.depthMin, neighbors, nzMld, + viewHocn, bkg, viewBathy, stdDevBkg, true, 6); + } + } // end 3D case + } // end jnode + } // end var + + /// Smooth the fields + // ------------------ + if (configD.simpleSmoothing) { + for (auto & var : configD.socaVars.variables()) { + // Skip the layer thickness variable + if (var == "hocn") { + continue; + } + + // Horizontal averaging + for (int iter = 0; iter < configD.niterHoriz; ++iter) { + // Update the halo points + nodeColumns.haloExchange(bkgErrFs[var]); + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + + // Loops through nodes and levels + for (atlas::idx_t level = 0; level < xbFs[var].shape(1); ++level) { + for (atlas::idx_t jnode = 0; jnode < xbFs["tocn"].shape(0); ++jnode) { + // Early exit if on a ghost cell + if (ghostView(jnode) > 0) { + continue; + } + + // Ocean or ice node, do something + std::vector local; + auto neighbors = get_neighbors_of_node(jnode); + int nbh = neighbors.size(); + for (int nn = 0; nn < neighbors.size(); ++nn) { + int nbNode = neighbors[nn]; + if ( abs(viewHocn(nbNode, level)) <= 0.1 ) { + continue; + } + local.push_back(stdDevBkg(nbNode, level)); + } + + if (local.size() > 2) { + stdDevBkg(jnode, level) = + std::accumulate(local.begin(), local.end(), 0.0) / local.size(); + } + + // Reset to 0 over land + if (abs(viewHocn(jnode, level)) <= 0.1) { + stdDevBkg(jnode, level) = 0.0; + } + } + } + } + + // Vertical averaging + if (xbFs[var].shape(1) == 1) { + oops::Log::info() << "skipping vertical smoothing for " << var << std::endl; + } else { + auto stdDevBkg = atlas::array::make_view(bkgErrFs[var]); + auto tmpArray(stdDevBkg); + for (int iter = 0; iter < configD.niterVert; ++iter) { + for (atlas::idx_t jnode = 0; jnode < xbFs["tocn"].shape(0); ++jnode) { + for (atlas::idx_t level = 1; level < xbFs[var].shape(1)-1; ++level) { + stdDevBkg(jnode, level) = (tmpArray(jnode, level-1) + + tmpArray(jnode, level) + + tmpArray(jnode, level+1)) / 3.0; + } + stdDevBkg(jnode, 0) = stdDevBkg(jnode, 1); + } + } + } + } + } + + /// Use explicit diffusion to smooth the background error + // ------------------------------------------------------ + // TODO(G): Use this once Travis adds the option to skip the normalization. + // The output is currently in [0, ~1000] + // Initialize the diffusion central block + if (fullConfig.has("diffusion")) { + const eckit::LocalConfiguration diffusionConfig(fullConfig, "diffusion"); + soca::ExplicitDiffusionParameters params; + params.deserialize(diffusionConfig); + oops::GeometryData geometryData(geom.functionSpace(), + bkgErrFs["tocn"], true, this->getComm()); + const oops::FieldSet3D dumyXb(configD.cycleDate, this->getComm()); + soca::ExplicitDiffusion diffuse(geometryData, configD.socaVars, + diffusionConfig, params, dumyXb, dumyXb); + diffuse.read(); + + // Smooth the field + oops::FieldSet3D dx(configD.cycleDate, this->getComm()); + dx.deepCopy(bkgErrFs); + diffuse.multiply(dx); + bkgErrFs = dx.fieldSet(); + } + + // Rescale + util::multiplyFieldSet(bkgErrFs, configD.rescale); + + // We want to write with soca, not atlas: Syncronize with soca Increment + bkgErr.fromFieldSet(bkgErrFs); + + // Save the background error + const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error"); + bkgErr.write(bkgErrorConfig); + + return 0; + } + + // ----------------------------------------------------------------------------- + + private: + std::string appname() const { + return "gdasapp::SocaDiagB"; + } + + // ----------------------------------------------------------------------------- + }; +} // namespace gdasapp diff --git a/utils/test/testref/ghrsst2ioda.test b/utils/test/testref/ghrsst2ioda.test index 1413f9d04..ab0a0b561 100644 --- a/utils/test/testref/ghrsst2ioda.test +++ b/utils/test/testref/ghrsst2ioda.test @@ -21,6 +21,6 @@ latitude: Max: 77.95 Sum: 623.423 datetime: - Min: 1269445504 - Max: 1269445504 - Sum: 10155564032 + Min: 1269445446 + Max: 1269445446 + Sum: 10155563568