diff --git a/build.sh b/build.sh index a27da74e5..3a9e98fa8 100755 --- a/build.sh +++ b/build.sh @@ -71,7 +71,7 @@ while getopts "p:t:c:hvdfa" opt; do done case ${BUILD_TARGET} in - hera | orion | hercules) + hera | orion | hercules | wcoss2 | noaacloud) echo "Building GDASApp on $BUILD_TARGET" source $dir_root/ush/module-setup.sh module use $dir_root/modulefiles diff --git a/modulefiles/GDAS/noaacloud.intel.lua b/modulefiles/GDAS/noaacloud.intel.lua new file mode 100644 index 000000000..b976d75bc --- /dev/null +++ b/modulefiles/GDAS/noaacloud.intel.lua @@ -0,0 +1,93 @@ +help([[ +Load environment for running the GDAS application with Intel compilers and MPI. +]]) + +local pkgName = myModuleName() +local pkgVersion = myModuleVersion() +local pkgNameVer = myModuleFullName() + +prepend_path("MODULEPATH", "/contrib/spack-stack/spack-stack-1.6.0/envs/unified-env/install/modulefiles/Core") + +-- below two lines get us access to the spack-stack modules +load("stack-intel/2021.3.0") +load("stack-intel-oneapi-mpi/2021.3.0") +load("cmake/3.23.1") +load("gettext/0.19.8.1") +--load("libunistring/1.1") +--load("libidn2/2.3.4") +load("pcre2/10.42") +load("curl/8.4.0") +load("zlib/1.2.13") +load("git/1.8.3.1") +load("pkg-config/0.27.1") +load("hdf5/1.14.0") +load("parallel-netcdf/1.12.2") +load("netcdf-c/4.9.2") +load("nccmp/1.9.0.1") +load("netcdf-fortran/4.6.1") +load("nco/5.0.6") +load("parallelio/2.5.10") +load("wget/1.14") +load("boost/1.83.0") +load("bufr/12.0.1") +load("git-lfs/2.4.1") +load("ecbuild/3.7.2") +load("openjpeg/2.3.1") +load("eccodes/2.32.0") +load("eigen/3.4.0") +load("openblas/0.3.24") +load("eckit/1.24.5") +load("fftw/3.3.10") +load("fckit/0.11.0") +load("fiat/1.2.0") +load("ectrans/1.2.0") +load("atlas/0.35.1") +load("sp/2.5.0") +load("gsl-lite/0.37.0") +load("libjpeg/2.1.0") +load("krb5/1.15.1") +load("libtirpc/1.3.3") +load("hdf/4.2.15") +load("jedi-cmake/1.4.0") +load("libpng/1.6.37") +load("libxt/1.1.5") +load("libxmu/1.1.4") +load("libxpm/3.5.12") +load("libxaw/1.0.13") +load("udunits/2.2.28") +load("ncview/2.1.9") +load("netcdf-cxx4/4.3.1") +load("json/3.10.5") +load("crtm/2.4.0.1") +load("rocoto/1.3.3") +load("prod_util/2.1.1") + +load("py-jinja2/3.0.3") +load("py-netcdf4/1.5.8") +load("py-pybind11/2.11.0") +load("py-pycodestyle/2.11.0") +load("py-pyyaml/6.0") +load("py-scipy/1.11.3") +load("py-xarray/2023.7.0") +load("py-f90nml/1.4.3") +load("py-pip/23.1.2") + +-- hack for wxflow +--prepend_path("PYTHONPATH", "/contrib/Wei.Huang/data/glopara/data/gdasapp-wxflow-20240307/src/wxflow") + +setenv("CC","mpiicc") +setenv("FC","mpiifort") +setenv("CXX","mpiicpc") + +local mpiexec = '/apps/slurm/default/bin/srun' +local mpinproc = '-n' +setenv('MPIEXEC_EXEC', mpiexec) +setenv('MPIEXEC_NPROC', mpinproc) + +setenv("CRTM_FIX","/contrib/Wei.Huang/data/hack-orion/crtm/2.4.0_skylab_3.0") +setenv("GDASAPP_TESTDATA","/contrib/Wei.Huang/data/hack-orion/data") + +whatis("Name: ".. pkgName) +whatis("Version: ".. pkgVersion) +whatis("Category: GDASApp") +whatis("Description: Load all libraries needed for GDASApp") diff --git a/modulefiles/GDAS/wcoss2.intel.lua b/modulefiles/GDAS/wcoss2.intel.lua new file mode 100644 index 000000000..d21099b2e --- /dev/null +++ b/modulefiles/GDAS/wcoss2.intel.lua @@ -0,0 +1,51 @@ +help([[ +Load environment for running the GDAS application with Intel compilers and MPI. +]]) + +local pkgName = myModuleName() +local pkgVersion = myModuleVersion() +local pkgNameVer = myModuleFullName() + +prepend_path("MODULEPATH", "/apps/dev/lmodules/core") + +load("PrgEnv-intel/8.2.0") +load("cmake/3.20.2") +load("craype") +load("cray-pals") +load("git/2.29.0") +load("intel/19.1.3.304") +load("cray-mpich/8.1.12") +load("hdf5/1.12.2") +load("netcdf/4.7.4") +load("udunits/2.2.28") +load("eigen/3.4.0") +load("boost/1.79.0") +load("gsl-lite/v0.40.0") +load("sp/2.4.0") +load("python/3.8.6") +load("ecbuild/3.7.0") +load("qhull/2020.2") +load("eckit/1.24.4") +load("fckit/0.11.0") +load("atlas/0.35.0") +load("nccmp") + +-- hack for pybind11 +setenv("pybind11_ROOT", "/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/lib/python3.8/site-packages/pybind11/share/cmake/pybind11") + +-- hack for wxflow +--prepend_path("PYTHONPATH", "/scratch1/NCEPDEV/da/python/gdasapp/wxflow/20240307/src") + +local mpiexec = '/pe/intel/compilers_and_libraries_2020.4.304/linux/mpi/intel64/bin/mpirun' +local mpinproc = '-n' +setenv('MPIEXEC_EXEC', mpiexec) +setenv('MPIEXEC_NPROC', mpinproc) + +setenv("CRTM_FIX","/lfs/h2/emc/da/noscrub/emc.da/GDASApp/fix/crtm/2.4.0") +setenv("GDASAPP_TESTDATA","/lfs/h2/emc/da/noscrub/emc.da/GDASApp/data") +setenv("GDASAPP_UNIT_TEST_DATA_PATH", "/lfs/h2/emc/da/noscrub/emc.da/GDASApp/data/test") + +whatis("Name: ".. pkgName) +whatis("Version: ".. pkgVersion) +whatis("Category: GDASApp") +whatis("Description: Load all libraries needed for GDASApp") diff --git a/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 new file mode 100644 index 000000000..2a0c2258a --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_n20 + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_n20.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_n20_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 new file mode 100644 index 000000000..ce8a377a2 --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_npp + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_npp.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_npp_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json new file mode 100644 index 000000000..eab19c6b2 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json @@ -0,0 +1,16 @@ +{ + "data_format" : "bufr_d", + "data_type" : "satwnd", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "subsets" : [ "NC005091" ], + "data_description" : "NC005091 S-NPP/NOAA-20 SATWIND VIIRS IR(LW)", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "VIIRS", "sensor_full_name": "Visible Infrared Imaging Radiometer Suite", "sensor_id": 999 }, + "satellite_info" : [ + { "satellite_name": "npp", "satellite_full_name": "SUOMI-NPP", "satellite_id": 224, "launch time": "YYYYMMDD" }, + { "satellite_name": "n20", "satellite_full_name": "NOAA-20", "satellite_id": 225, "launch time": "YYYYMMDD" }, + ] +} diff --git a/parm/soca/berror/soca_diagb.yaml.j2 b/parm/soca/berror/soca_diagb.yaml.j2 index 3e6b9d2ec..fcfd7442f 100644 --- a/parm/soca/berror/soca_diagb.yaml.j2 +++ b/parm/soca/berror/soca_diagb.yaml.j2 @@ -7,8 +7,8 @@ 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' + ocn_filename: 'gdas.ocean.t{{ gcyc }}z.inst.f009.nc' + ice_filename: 'gdas.agg_ice.t{{ gcyc }}z.inst.f009.nc' read_from_file: 1 background error: diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 54cce2363..d207840f0 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -50,6 +50,7 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): cyc = str(os.getenv('cyc')).zfill(2) bcyc = str((int(cyc) - 3) % 24).zfill(2) gcyc = str((int(cyc) - 6) % 24).zfill(2) # previous cycle +cdatedt = datetime.strptime(cdate, '%Y%m%d%H') 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') @@ -93,8 +94,12 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{bcyc}z.ocngrid.nc')]) # Copy the CICE analysis restart -cdateice = pdy + '.' + cyc + '0000' -post_file_list.append([os.path.join(anl_dir, 'Data', f'{cdateice}.cice_model.res.nc'), +if os.getenv('DOIAU') == "YES": + cice_rst_date = bdatedt.strftime('%Y%m%d.%H%M%S') +else: + cice_rst_date = cdatedt.strftime('%Y%m%d.%H%M%S') + +post_file_list.append([os.path.join(anl_dir, 'Data', f'{cice_rst_date}.cice_model.res.nc'), os.path.join(com_ice_restart, f'{cdate}.cice_model_anl.res.nc')]) FileHandler({'copy': post_file_list}).sync() @@ -171,7 +176,7 @@ def create_obs_space(data): # TODO(GorA): this should be setup properly in the g-w once gdassoca_obsstats is in develop gdassoca_obsstats_exec = os.path.join(os.getenv('HOMEgfs'), 'sorc', 'gdas.cd', 'build', 'bin', 'gdassoca_obsstats.x') -command = f"{os.getenv('launcher')} {gdassoca_obsstats_exec} {stats_yaml}" +command = f"{os.getenv('launcher')} -n 1 {gdassoca_obsstats_exec} {stats_yaml}" logger.info(f"{command}") result = subprocess.run(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index d80bfc834..265526eb0 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -70,10 +70,6 @@ def find_clim_ens(input_date): 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 -else: - dohybvar = False # create analysis directories diags = os.path.join(anl_dir, 'diags') # output dir for soca DA obs space @@ -89,12 +85,25 @@ def find_clim_ens(input_date): window_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - half_assim_freq window_begin_iso = window_begin.strftime('%Y-%m-%dT%H:%M:%SZ') window_middle_iso = window_middle.strftime('%Y-%m-%dT%H:%M:%SZ') -fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') RUN = os.getenv('RUN') cyc = os.getenv('cyc') gcyc = os.getenv('gcyc') PDY = os.getenv('PDY') +# hybrid-envar switch +if os.getenv('DOHYBVAR') == "YES": + dohybvar = True +else: + dohybvar = False + +# switch for the cycling type +if os.getenv('DOIAU') == "YES": + # forecast initialized at the begining of the DA window + fcst_begin = window_begin +else: + # forecast initialized at the middle of the DA window + fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') + ################################################################################ # fetch observations @@ -327,7 +336,7 @@ def find_clim_ens(input_date): ################################################################################ # Copy initial condition -bkg_utils.stage_ic(bkg_dir, anl_dir, RUN, gcyc) +bkg_utils.stage_ic(bkg_dir, anl_dir, gcyc) ################################################################################ # prepare input.nml diff --git a/sorc/ioda b/sorc/ioda index 91eb2c764..206eba708 160000 --- a/sorc/ioda +++ b/sorc/ioda @@ -1 +1 @@ -Subproject commit 91eb2c7643b33c1af2470d289eefb7aec2667387 +Subproject commit 206eba7084c08a3bc9a1c6ccbbff3e63d7cb602e diff --git a/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py new file mode 100755 index 000000000..18f87db04 --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py @@ -0,0 +1,470 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import numpy.ma as ma +from pyiodaconv import bufr +import calendar +import json +import time +import math +import datetime +import os +from datetime import datetime +from pyioda import ioda_obs_space as ioda_ospace +from wxflow import Logger + +# ==================================================================== +# Satellite Winds (AMV) BUFR dump file for VIIRS/S-NPP,NOAA-20 +# ==================================================================== +# Subset | Spectral Band | Code (002023) | ObsType +# -------------------------------------------------------------------- +# NC005091 | IRLW (Freq < 5E+13) | Method 1 | 260 +# ==================================================================== + +# Define and initialize global variables +global float32_fill_value +global int32_fill_value +global int64_fill_value + +float32_fill_value = np.float32(0) +int32_fill_value = np.int32(0) +int64_fill_value = np.int64(0) + + +def Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd): + + uob = (-wspd * np.sin(np.radians(wdir))).astype(np.float32) + vob = (-wspd * np.cos(np.radians(wdir))).astype(np.float32) + + return uob, vob + + +def Get_ObsType(swcm, chanfreq): + + obstype = swcm.copy() + + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 260, obstype) # IRLW + + if not np.any(np.isin(obstype, [260])): + raise ValueError("Error: Unassigned ObsType found ... ") + + return obstype + + +def bufr_to_ioda(config, logger): + + subsets = config["subsets"] + logger.debug(f"Checking subsets = {subsets}") + + # Get parameters from configuration + subsets = config["subsets"] + data_format = config["data_format"] + data_type = config["data_type"] + data_description = config["data_description"] + data_provider = config["data_provider"] + cycle_type = config["cycle_type"] + dump_dir = config["dump_directory"] + ioda_dir = config["ioda_directory"] + cycle = config["cycle_datetime"] + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + + satellite_info_array = config["satellite_info"] + sensor_name = config["sensor_info"]["sensor_name"] + sensor_full_name = config["sensor_info"]["sensor_full_name"] + sensor_id = config["sensor_info"]["sensor_id"] + + # Get derived parameters + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + reference_time = datetime.strptime(cycle, "%Y%m%d%H") + reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") + + # General informaton + converter = 'BUFR to IODA Converter' + process_level = 'Level-2' + platform_description = 'TERRA/AQUA' + sensor_description = 'Moderate Resolution Imaging Spectroradiometer' + + logger.info(f'sensor_name = {sensor_name}') + logger.info(f'sensor_full_name = {sensor_full_name}') + logger.info(f'sensor_id = {sensor_id}') + logger.info(f'reference_time = {reference_time}') + + bufrfile = f"{cycle_type}.t{hh}z.{data_type}.tm00.{data_format}" + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), 'atmos', bufrfile) + if not os.path.isfile(DATA_PATH): + logger.info(f"DATA_PATH {DATA_PATH} does not exist") + return + logger.debug(f"The DATA_PATH is: {DATA_PATH}") + + # ============================================ + # Make the QuerySet for all the data we want + # ============================================ + start_time = time.time() + + logger.info('Making QuerySet') + q = bufr.QuerySet(subsets) + + # MetaData + q.add('latitude', '*/CLATH') + q.add('longitude', '*/CLONH') + q.add('satelliteId', '*/SAID') + q.add('year', '*/YEAR') + q.add('month', '*/MNTH') + q.add('day', '*/DAYS') + q.add('hour', '*/HOUR') + q.add('minute', '*/MINU') + q.add('second', '*/SECO') + q.add('satelliteZenithAngle', '*/SAZA') + q.add('sensorCentralFrequency', '*/SCCF') + q.add('pressure', '*/PRLC[1]') + + # Processing Center + q.add('dataProviderOrigin', '*/OGCE[1]') + q.add('windGeneratingApplication', '*/AMVQIC/GNAPS') + +# # Quality Infomation (Quality Indicator w/o forecast) + q.add('qualityInformationWithoutForecast', '*/AMVQIC/PCCF') + + # Wind Retrieval Method Information + q.add('windComputationMethod', '*/SWCM') + + # ObsValue + q.add('windDirection', '*/WDIR') + q.add('windSpeed', '*/WSPD') + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for making QuerySet : {running_time} seconds') + + # ============================================================== + # Open the BUFR file and execute the QuerySet to get ResultSet + # Use the ResultSet returned to get numpy arrays of the data + # ============================================================== + start_time = time.time() + + logger.info('Executing QuerySet to get ResultSet') + with bufr.File(DATA_PATH) as f: + try: + r = f.execute(q) + except Exception as err: + logger.info(f'Return with {err}') + return + + # MetaData + satid = r.get('satelliteId') + year = r.get('year') + month = r.get('month') + day = r.get('day') + hour = r.get('hour') + minute = r.get('minute') + second = r.get('second') + lat = r.get('latitude') + lon = r.get('longitude') + satzenang = r.get('satelliteZenithAngle') + pressure = r.get('pressure', type='float') + chanfreq = r.get('sensorCentralFrequency', type='float') + + # Processing Center + ogce = r.get('dataProviderOrigin') + ga = r.get('windGeneratingApplication') + + # Quality Information + qi = r.get('qualityInformationWithoutForecast', type='float') + # For TERRA/AQUA MODIS data, qi w/o forecast (qifn) is packaged in same + # vector of qi with ga = 5 (QI without forecast), and EE is + # packaged in same vector of qi with ga=7 (Estimated Error (EE) in m/s + # converted to a percent confidence)shape (4,nobs). Must conduct a + # search and extract the correct vector for gnap and qi + # 1. Find dimension-sizes of ga and qi (should be the same!) + gDim1, gDim2 = np.shape(ga) + qDim1, qDim2 = np.shape(qi) + logger.info(f'Generating Application and Quality Information SEARCH:') + logger.info(f'Dimension size of GNAP ({gDim1},{gDim2})') + logger.info(f'Dimension size of PCCF ({qDim1},{qDim2})') + # 2. Initialize gnap and qifn as None, and search for dimension of + # ga with values of 5. If the same column exists for qi, assign + # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no + # appropriate GNAP/PCCF combination was found + gnap = None + qifn = None + for i in range(gDim2): + if np.unique(ga[:, i].squeeze()) == 5: + if i <= qDim2: + logger.info(f'GNAP/PCCF found for column {i}') + gnap = ga[:, i].squeeze() + qifn = qi[:, i].squeeze() + else: + logger.info(f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') + if (gnap is None) & (qifn is None): + logger.info(f'ERROR: GNAP == 5 NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') + # If EE is needed, key search on np.unique(ga[:,i].squeeze()) == 7 instead + + # Wind Retrieval Method Information + swcm = r.get('windComputationMethod') + # ObsValue + # Wind direction and Speed + wdir = r.get('windDirection', type='float') + wspd = r.get('windSpeed') + + # DateTime: seconds since Epoch time + # IODA has no support for numpy datetime arrays dtype=datetime64[s] + timestamp = r.get_datetime('year', 'month', 'day', 'hour', 'minute', 'second').astype(np.int64) + + # Check BUFR variable generic dimension and type + + # Global variables declaration + # Set global fill values + float32_fill_value = satzenang.fill_value + int32_fill_value = satid.fill_value + int64_fill_value = timestamp.fill_value.astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + logger.info('Creating derived variables') + logger.debug('Creating derived variables - wind components (uob and vob)') + + uob, vob = Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd) + + logger.debug(f' uob min/max = {uob.min()} {uob.max()}') + logger.debug(f' vob min/max = {vob.min()} {vob.max()}') + + obstype = Get_ObsType(swcm, chanfreq) + + height = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + stnelev = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for creating derived variables : {running_time} seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Create IODA ObsSpace and Write IODA output based on satellite ID') + + # Find unique satellite identifiers in data to process + unique_satids = np.unique(satid) + logger.info(f'Number of Unique satellite identifiers: {len(unique_satids)}') + logger.info(f'Unique satellite identifiers: {unique_satids}') + + logger.debug(f'Loop through unique satellite identifier {unique_satids}') + total_ob_processed = 0 + for sat in unique_satids.tolist(): + start_time = time.time() + + matched = False + for satellite_info in satellite_info_array: + if (satellite_info["satellite_id"] == sat): + matched = True + satellite_id = satellite_info["satellite_id"] + satellite_name = satellite_info["satellite_name"] + satinst = sensor_name.lower()+'_'+satellite_name.lower() + logger.debug(f'Split data for {satinst} satid = {sat}') + + if matched: + + # Define a boolean mask to subset data from the original data object + mask = satid == sat + # MetaData + lon2 = lon[mask] + lat2 = lat[mask] + timestamp2 = timestamp[mask] + satid2 = satid[mask] + satzenang2 = satzenang[mask] + chanfreq2 = chanfreq[mask] + obstype2 = obstype[mask] + pressure2 = pressure[mask] + height2 = height[mask] + stnelev2 = stnelev[mask] + + # Processing Center + ogce2 = ogce[mask] + + # QC Info + qifn2 = qifn[mask] + + # Method + swcm2 = swcm[mask] + + # ObsValue + wdir2 = wdir[mask] + wspd2 = wspd[mask] + uob2 = uob[mask] + vob2 = vob[mask] + + # Timestamp Range + timestamp2_min = datetime.fromtimestamp(timestamp2.min()) + timestamp2_max = datetime.fromtimestamp(timestamp2.max()) + + # Check unique observation time + unique_timestamp2 = np.unique(timestamp2) + logger.debug(f'Processing output for satid {sat}') + + # Create the dimensions + dims = { + 'Location': np.arange(0, wdir2.shape[0]) + } + + # Create IODA ObsSpace + iodafile = f"{cycle_type}.t{hh}z.{data_type}.{satinst}.tm00.nc" + OUTPUT_PATH = os.path.join(ioda_dir, iodafile) + logger.info(f'Create output file : {OUTPUT_PATH}') + obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) + + # Create Global attributes + logger.debug('Write global attributes') + obsspace.write_attr('Converter', converter) + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('sensorLongDescription', sensor_description) + + # Create IODA variables + logger.debug('Write variables: name, type, units, and attributes') + # Longitude + obsspace.create_var('MetaData/longitude', dtype=lon2.dtype, fillval=lon2.fill_value) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon2) + + # Latitude + obsspace.create_var('MetaData/latitude', dtype=lat.dtype, fillval=lat2.fill_value) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat2) + + # Datetime + obsspace.create_var('MetaData/dateTime', dtype=np.int64, fillval=int64_fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Datetime') \ + .write_data(timestamp2) + + # Satellite Identifier + obsspace.create_var('MetaData/satelliteIdentifier', dtype=satid2.dtype, fillval=satid2.fill_value) \ + .write_attr('long_name', 'Satellite Identifier') \ + .write_data(satid2) + + # Sensor Zenith Angle + obsspace.create_var('MetaData/satelliteZenithAngle', dtype=satzenang2.dtype, fillval=satzenang2.fill_value) \ + .write_attr('units', 'degree') \ + .write_attr('valid_range', np.array([0, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Satellite Zenith Angle') \ + .write_data(satzenang2) + + # Sensor Centrall Frequency + obsspace.create_var('MetaData/sensorCentralFrequency', dtype=chanfreq2.dtype, fillval=chanfreq2.fill_value) \ + .write_attr('units', 'Hz') \ + .write_attr('long_name', 'Satellite Channel Center Frequency') \ + .write_data(chanfreq2) + + # Data Provider + obsspace.create_var('MetaData/dataProviderOrigin', dtype=ogce2.dtype, fillval=ogce2.fill_value) \ + .write_attr('long_name', 'Identification of Originating/Generating Center') \ + .write_data(ogce2) + + # Quality: Percent Confidence - Quality Information Without Forecast + obsspace.create_var('MetaData/qualityInformationWithoutForecast', dtype=qifn2.dtype, fillval=qifn2.fill_value) \ + .write_attr('long_name', 'Quality Information Without Forecast') \ + .write_data(qifn2) + + # Wind Computation Method + obsspace.create_var('MetaData/windComputationMethod', dtype=swcm2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Satellite-derived Wind Computation Method') \ + .write_data(swcm2) + + # Pressure + obsspace.create_var('MetaData/pressure', dtype=pressure2.dtype, fillval=pressure2.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Pressure') \ + .write_data(pressure2) + + # Height (mimic prepbufr) + obsspace.create_var('MetaData/height', dtype=height2.dtype, fillval=height2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Height of Observation') \ + .write_data(height2) + + # Station Elevation (mimic prepbufr) + obsspace.create_var('MetaData/stationElevation', dtype=stnelev2.dtype, fillval=stnelev2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Station Elevation') \ + .write_data(stnelev2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windEastward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windNorthward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # U-Wind Component + obsspace.create_var('ObsValue/windEastward', dtype=uob2.dtype, fillval=uob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Eastward Wind Component') \ + .write_data(uob2) + + # V-Wind Component + obsspace.create_var('ObsValue/windNorthward', dtype=vob2.dtype, fillval=vob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Northward Wind Component') \ + .write_data(vob2) + + end_time = time.time() + running_time = end_time - start_time + total_ob_processed += len(satid2) + logger.debug(f'Number of observation processed : {len(satid2)}') + logger.debug(f'Processing time for splitting and output IODA for {satinst} : {running_time} seconds') + + else: + logger.info(f"Do not find this satellite id in the configuration: satid = {sat}") + + logger.info("All Done!") + logger.info(f'Total number of observation processed : {total_ob_processed}') + + +if __name__ == '__main__': + + start_time = time.time() + + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config', type=str, help='Input JSON configuration', required=True) + parser.add_argument('-v', '--verbose', help='print debug logging information', + action='store_true') + args = parser.parse_args() + + log_level = 'DEBUG' if args.verbose else 'INFO' + logger = Logger('BUFR2IODA_satwind_amv_ahi.py', level=log_level, colored_log=True) + + with open(args.config, "r") as json_file: + config = json.load(json_file) + + bufr_to_ioda(config, logger) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f"Total running time: {running_time} seconds") diff --git a/ush/module-setup.sh b/ush/module-setup.sh index ef7f73699..6a2aabf73 100755 --- a/ush/module-setup.sh +++ b/ush/module-setup.sh @@ -105,6 +105,10 @@ elif [[ $MACHINE_ID = discover* ]]; then export PATH=$PATH:$SPACK_ROOT/bin . $SPACK_ROOT/share/spack/setup-env.sh +elif [[ $MACHINE_ID = noaacloud* ]]; then + # We are on NOAA Cloud + module purge + else echo WARNING: UNKNOWN PLATFORM 1>&2 fi diff --git a/ush/soca/bkg_utils.py b/ush/soca/bkg_utils.py index 08be3dd89..488fcf697 100755 --- a/ush/soca/bkg_utils.py +++ b/ush/soca/bkg_utils.py @@ -101,13 +101,12 @@ def gen_bkg_list(bkg_path, out_path, window_begin=' ', yaml_name='bkg.yaml', ice bkg_date = window_begin # Construct list of background file names - RUN = os.getenv('RUN') cyc = str(os.getenv('cyc')).zfill(2) gcyc = str((int(cyc) - 6) % 24).zfill(2) # previous cycle fcst_hrs = list(range(3, 10, dt_pseudo)) files = [] for fcst_hr in fcst_hrs: - files.append(os.path.join(bkg_path, f'{RUN}.ocean.t'+gcyc+'z.inst.f'+str(fcst_hr).zfill(3)+'.nc')) + files.append(os.path.join(bkg_path, f'gdas.ocean.t'+gcyc+'z.inst.f'+str(fcst_hr).zfill(3)+'.nc')) # Identify the ocean background that will be used for the vertical coordinate remapping ocn_filename_ic = os.path.splitext(os.path.basename(files[0]))[0]+'.nc' @@ -159,17 +158,17 @@ def gen_bkg_list(bkg_path, out_path, window_begin=' ', yaml_name='bkg.yaml', ice FileHandler({'copy': bkg_list_src_dst}).sync() -def stage_ic(bkg_dir, anl_dir, RUN, gcyc): +def stage_ic(bkg_dir, anl_dir, gcyc): # Copy and rename initial condition ics_list = [] # ocean IC's - mom_ic_src = os.path.join(bkg_dir, f'{RUN}.ocean.t{gcyc}z.inst.f003.nc') + mom_ic_src = os.path.join(bkg_dir, f'gdas.ocean.t{gcyc}z.inst.f003.nc') mom_ic_dst = os.path.join(anl_dir, 'INPUT', 'MOM.res.nc') ics_list.append([mom_ic_src, mom_ic_dst]) # seaice IC's - cice_ic_src = os.path.join(bkg_dir, f'{RUN}.agg_ice.t{gcyc}z.inst.f003.nc') + cice_ic_src = os.path.join(bkg_dir, f'gdas.agg_ice.t{gcyc}z.inst.f003.nc') cice_ic_dst = os.path.join(anl_dir, 'INPUT', 'cice.res.nc') ics_list.append([cice_ic_src, cice_ic_dst]) FileHandler({'copy': ics_list}).sync() diff --git a/ush/soca/marine_recenter.py b/ush/soca/marine_recenter.py index 4c49bb530..5bb689e71 100644 --- a/ush/soca/marine_recenter.py +++ b/ush/soca/marine_recenter.py @@ -128,7 +128,7 @@ def initialize(self): ################################################################################ # Copy initial condition - bkg_utils.stage_ic(self.config.bkg_dir, self.runtime_config.DATA, RUN, gcyc) + bkg_utils.stage_ic(self.config.bkg_dir, self.runtime_config.DATA, gcyc) ################################################################################ # stage ensemble members diff --git a/utils/obsproc/IcecAmsr2Ioda.h b/utils/obsproc/IcecAmsr2Ioda.h index 681bbcd8d..3241d3ae2 100644 --- a/utils/obsproc/IcecAmsr2Ioda.h +++ b/utils/obsproc/IcecAmsr2Ioda.h @@ -13,8 +13,10 @@ #include // NOLINT +#include "ioda/../../../../core/IodaUtils.h" #include "ioda/Group.h" #include "ioda/ObsGroup.h" + #include "oops/util/dateFunctions.h" #include "NetCDFToIodaConverter.h" @@ -83,22 +85,34 @@ namespace gdasapp { int minute = oneTmpdateTimeVal[i+4]; int second = static_cast(oneTmpdateTimeVal[i+5]); - // Replace Fillvalue -9999 to 0 to avoid crash in dateToJulian + // Avoid crash util in ioda::convertDtimeToTimeOffsets if (year == -9999 || month == -9999 || day == -9999 || hour == -9999 || minute == -9999 || second == -9999) { year = month = day = hour = minute = second = 0; } - // Convert a date to Julian date - uint64_t julianDate = util::datefunctions::dateToJulian(year, month, day); - - // Subtract Julian day from January 1, 1970 (convert to epoch) - int daysSinceEpoch = julianDate - 2440588; - - // Calculate seconds only from HHMMSS - int secondsOffset = util::datefunctions::hmsToSeconds(hour, minute, second); - - iodaVars.datetime_(index) = static_cast(daysSinceEpoch*86400.0f) + secondsOffset; + // Construct iso8601 string format for each dateTime + std::stringstream ss; + ss << std::setfill('0') + << std::setw(4) << year << '-' + << std::setw(2) << month << '-' + << std::setw(2) << day << 'T' + << std::setw(2) << hour << ':' + << std::setw(2) << minute << ':' + << std::setw(2) << second << 'Z'; + std::string formattedDateTime = ss.str(); + util::DateTime dateTime(formattedDateTime); + + // Set epoch time for AMSR2_ICEC + util::DateTime epochDtime("1970-01-01T00:00:00Z"); + + // Convert Obs DateTime objects to epoch time offsets in seconds + // 0000-00-00T00:00:00Z will be converterd to negative seconds + int64_t timeOffsets + = ioda::convertDtimeToTimeOffsets(epochDtime, {dateTime})[0]; + + // Update datetime Eigen Arrays + iodaVars.datetime_(index) = timeOffsets; index++; } @@ -106,7 +120,7 @@ namespace gdasapp { for (int i = 0; i < iodaVars.location_; i++) { iodaVars.longitude_(i) = oneDimLonVal[i]; iodaVars.latitude_(i) = oneDimLatVal[i]; - iodaVars.obsVal_(i) = static_cast(oneDimObsVal[i]*0.01f); + iodaVars.obsVal_(i) = static_cast(oneDimObsVal[i]*0.01); iodaVars.obsError_(i) = 0.1; // Do something for obs error iodaVars.preQc_(i) = oneDimFlagsVal[i]; // Store optional metadata, set ocean basins to -999 for now diff --git a/ush/soca/gdassoca_obsstats.py b/utils/soca/fig_gallery/gdassoca_obsstats.py similarity index 68% rename from ush/soca/gdassoca_obsstats.py rename to utils/soca/fig_gallery/gdassoca_obsstats.py index 4741687bb..3f92f799d 100644 --- a/ush/soca/gdassoca_obsstats.py +++ b/utils/soca/fig_gallery/gdassoca_obsstats.py @@ -10,6 +10,21 @@ import matplotlib.pyplot as plt import matplotlib.dates as mdates +colors = [ + "lightsteelblue", + "lightgreen", + "lightpink", + "lightsalmon", + "lightcoral", + "lightgoldenrodyellow", + "paleturquoise", + "palegreen", + "palegoldenrod", + "peachpuff", + "mistyrose", + "lavender" +] + class ObsStats: def __init__(self): @@ -19,14 +34,15 @@ def read_csv(self, filepaths): # Iterate through the list of file paths and append their data for filepath in filepaths: new_data = pd.read_csv(filepath) + # Convert date to datetime for easier plotting new_data['date'] = pd.to_datetime(new_data['date'], format='%Y%m%d%H') self.data = pd.concat([self.data, new_data], ignore_index=True) + self.data.sort_values('date', inplace=True) - def plot_timeseries(self, ocean, variable): + def plot_timeseries(self, ocean, variable, inst="", dirout=""): # Filter data for the given ocean and variable filtered_data = self.data[(self.data['Ocean'] == ocean) & (self.data['Variable'] == variable)] - if filtered_data.empty: print("No data available for the given ocean and variable combination.") return @@ -36,13 +52,16 @@ def plot_timeseries(self, ocean, variable): # Plot settings fig, axs = plt.subplots(3, 1, figsize=(10, 15), sharex=True) - fig.suptitle(f'{variable} statistics, {ocean} ocean', fontsize=16, fontweight='bold') + fig.suptitle(f'{inst} {variable} statistics, {ocean} ocean', fontsize=16, fontweight='bold') + exp_counter = 0 for exp in experiments: - exp_data = filtered_data[filtered_data['Exp'] == exp] + exp_data = self.data[(self.data['Ocean'] == ocean) & + (self.data['Variable'] == variable) & + (self.data['Exp'] == exp)] # Plot RMSE - axs[0].plot(exp_data['date'], exp_data['RMSE'], marker='o', linestyle='-', label=exp) + axs[0].plot(exp_data['date'], exp_data['RMSE'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[0].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[0].xaxis.set_major_locator(mdates.DayLocator()) axs[0].tick_params(labelbottom=False) @@ -51,42 +70,41 @@ def plot_timeseries(self, ocean, variable): axs[0].grid(True) # Plot Bias - axs[1].plot(exp_data['date'], exp_data['Bias'], marker='o', linestyle='-', label=exp) + axs[1].plot(exp_data['date'], exp_data['Bias'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[1].xaxis.set_major_locator(mdates.DayLocator()) axs[1].tick_params(labelbottom=False) axs[1].set_ylabel('Bias') - axs[1].legend() axs[1].grid(True) # Plot Count - axs[2].plot(exp_data['date'], exp_data['Count'], marker='o', linestyle='-', label=exp) + axs[2].plot(exp_data['date'], exp_data['Count'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[2].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[2].xaxis.set_major_locator(mdates.DayLocator()) axs[2].set_ylabel('Count') - axs[2].legend() axs[2].grid(True) # Improve layout and show plot plt.tight_layout(rect=[0, 0.03, 1, 0.95]) - plt.savefig(f'{ocean}_test.png') - + plt.savefig(f'{dirout}/{inst}_{variable}_{ocean}.png') if __name__ == "__main__": - epilog = ["Usage examples: ./gdassoca_obsstats.py --exp gdas.*.csv"] + epilog = ["Usage examples: ./gdassoca_obsstats.py --exps cp1/COMROOT/cp1 cp2/COMROOT/cp2 --inst sst_abi_g16_l3c --dirout cp1vscp2"] parser = argparse.ArgumentParser(description="Observation space RMSE's and BIAS's", formatter_class=argparse.RawDescriptionHelpFormatter, epilog=os.linesep.join(epilog)) parser.add_argument("--exps", nargs='+', required=True, help="Path to the experiment's COMROOT") - parser.add_argument("--variable", required=True, help="ombg_qc or ombg_noqc") - parser.add_argument("--figname", required=True, help="The name of the output figure") + parser.add_argument("--inst", required=True, help="The name of the instrument/platform (ex: sst_abi_g16_l3c)") + parser.add_argument("--dirout", required=True, help="Output directory") args = parser.parse_args() flist = [] + inst = args.inst + os.makedirs(args.dirout, exist_ok=True) + for exp in args.exps: - print(exp) - wc = exp + '/*.*/??/analysis/ocean/*.t??z.stats.csv' + wc = exp + f'/*.*/??/analysis/ocean/*{inst}*.stats.csv' flist.append(glob.glob(wc)) flist = sum(flist, []) @@ -94,4 +112,4 @@ def plot_timeseries(self, ocean, variable): obsStats.read_csv(flist) for var, ocean in product(['ombg_noqc', 'ombg_qc'], ['Atlantic', 'Pacific', 'Indian', 'Arctic', 'Southern']): - obsStats.plot_timeseries(ocean, var) + obsStats.plot_timeseries(ocean, var, inst=inst, dirout=args.dirout) diff --git a/utils/soca/fig_gallery/time_series_omb.html b/utils/soca/fig_gallery/time_series_omb.html new file mode 100644 index 000000000..4e06b5e56 --- /dev/null +++ b/utils/soca/fig_gallery/time_series_omb.html @@ -0,0 +1,120 @@ + + + + + + + + +
+ + + +
+
+ + + + + + +
+
+ + +
+ +
+ + + + diff --git a/utils/soca/gdassoca_obsstats.h b/utils/soca/gdassoca_obsstats.h index 98e70bf3b..5f95574a1 100644 --- a/utils/soca/gdassoca_obsstats.h +++ b/utils/soca/gdassoca_obsstats.h @@ -52,32 +52,18 @@ namespace gdasapp { std::vector obsSpaces; fullConfig.get("obs spaces", obsSpaces); - // check if the number of workers is correct. Only the serial and 1 diag file - // per pe works for now - ASSERT(getComm().size() <= 1 || getComm().size() >= obsSpaces.size()); - - // initialize pe's color and communicator's name - int color(0); - std::string commNameStr = "comm_idle"; - if (this->getComm().rank() < obsSpaces.size()) { - color = 1; - std::string commNameStr = "comm_work"; - } - - // Create the new communicator () - char const *commName = commNameStr.c_str(); - eckit::mpi::Comm & commObsSpace = this->getComm().split(color, commName); + // only the serial case works for now. + ASSERT(getComm().size() == 1); - // spread out the work over pes - if (color > 0) { + for (int i = 0; i < obsSpaces.size(); i++) { // get the obs space configuration - auto obsSpace = obsSpaces[commObsSpace.rank()]; + auto obsSpace = obsSpaces[i]; eckit::LocalConfiguration obsConfig(obsSpace, "obs space"); // get the obs diag file std::string obsFile; obsConfig.get("obsdatain.engine.obsfile", obsFile); - oops::Log::info() << "========= Processing" << obsFile + oops::Log::info() << "========= Processing " << obsFile << " date: " << extractDateFromFilename(obsFile) << std::endl; @@ -86,8 +72,7 @@ namespace gdasapp { obsSpace.get("variable", variable); // read the obs space - ioda::ObsSpace ospace(obsConfig, commObsSpace, timeWindow, - oops::mpi::myself()); + ioda::ObsSpace ospace(obsConfig, getComm(), timeWindow, getComm()); const size_t nlocs = ospace.nlocs(); oops::Log::info() << "nlocs =" << nlocs << std::endl; std::vector var(nlocs); @@ -125,7 +110,6 @@ namespace gdasapp { // Close the file outputFile.close(); } - getComm().barrier(); return EXIT_SUCCESS; }