diff --git a/.gitignore b/.gitignore index 2d99d7c84..6f0f7e656 100644 --- a/.gitignore +++ b/.gitignore @@ -49,6 +49,7 @@ rtcoef* advance_time closest_member_tool +cmems_sst_to_obs compare_states compute_error create_fixed_network_seq diff --git a/guide/available-observation-converters.rst b/guide/available-observation-converters.rst index 1623246c6..a074b8d6e 100644 --- a/guide/available-observation-converters.rst +++ b/guide/available-observation-converters.rst @@ -21,6 +21,7 @@ Each directory has at least one converter: - ``BATS``: :ref:`bats` - ``CHAMP``: :ref:`champ` - ``cice``: :ref:`cice_to_obs` +- ``cmems_sst_l3s``: :ref:`cmems_sst_to_obs` - ``CNOFS``: See ``DART/observations/obs_converters/CNOFS`` - ``CONAGUA``: :ref:`conagua` - ``COSMOS``: :ref:`cosmos` diff --git a/index.rst b/index.rst index abcc57774..4cac87f4e 100644 --- a/index.rst +++ b/index.rst @@ -335,6 +335,7 @@ References observations/obs_converters/CHAMP/README observations/obs_converters/BATS/readme observations/obs_converters/cice/cice_to_obs + observations/obs_converters/cmems_sst_l3s/readme observations/obs_converters/CONAGUA/README observations/obs_converters/COSMOS/COSMOS_to_obs observations/obs_converters/COSMOS/COSMOS_development diff --git a/observations/obs_converters/cmems_sst_l3s/cmems_sst_to_obs.f90 b/observations/obs_converters/cmems_sst_l3s/cmems_sst_to_obs.f90 new file mode 100644 index 000000000..98a3532e8 --- /dev/null +++ b/observations/obs_converters/cmems_sst_l3s/cmems_sst_to_obs.f90 @@ -0,0 +1,309 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download + +! The following is an obs converter for Copernicus Sea Surface +! Temperature L3S product. The data is based on a blended global +! high resolution ODYSSEA SST Multi-sensor L3 composite. + +! The converter reads in the adjusted SST (bias corrected) +! and its associated error sd from the incoming netcdf file. +! A lower bound on the SD is applied to avoid over-fitting +! a relatively raw data product. Observations from multiple +! files (at different times) are collected and added to the +! observation sequence file. + +! The observed temperature is converted from K to C. + + +program cmems_sst_to_obs + +use types_mod, only : i4, i8, r8, t_kelvin, MISSING_R8, digits12 +use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, get_time, & + set_date, set_time, print_date, operator(+), operator(-) +use utilities_mod, only : initialize_utilities, find_namelist_in_file, E_MSG, & + nmlfileunit, error_handler, do_nml_term, E_ERR, & + finalize_utilities, do_nml_file, get_next_filename, & + find_textfile_dims, file_exist, check_namelist_read +use location_mod, only : VERTISSURFACE +use obs_sequence_mod, only : obs_type, obs_sequence_type, init_obs, get_num_obs, & + static_init_obs_sequence, init_obs_sequence, & + set_copy_meta_data, set_qc_meta_data, write_obs_seq, & + destroy_obs_sequence, destroy_obs +use obs_utilities_mod, only : create_3d_obs, add_obs_to_seq +use obs_kind_mod, only : SATELLITE_BLENDED_SST +use netcdf_utilities_mod, only : nc_check, nc_open_file_readonly, nc_close_file, & + nc_get_variable, nc_get_attribute_from_variable, & + nc_get_dimension_size +use netcdf + +implicit none + +character(len=*), parameter :: source = 'cmems_sst_to_obs' + +! Lower bound for obs_error_sd +real(r8), parameter :: OBS_ERROR_SD_MIN = 0.05_r8 ! Add this to the namelist? +integer, parameter :: OBS_QC_LOW_QUALITY = 3 + +! File variables +character(len=512) :: string1 +character(len=256) :: next_infile +integer :: io, iunit, filenum +integer :: num_new_obs, nfiles +logical :: first_obs = .true. + +! Obs sequence +type(obs_sequence_type) :: obs_seq +type(obs_type) :: obs, prev_obs +type(time_type) :: obs_time, prev_time +integer :: num_copies = 1, & ! number of copies in sequence + num_qc = 1 ! number of QC entries +real(r8) :: obs_qc = 0.0_r8 + +! Data arrays +real(r8), allocatable :: lon(:), lat(:) ! Location +real(r8), allocatable :: sst(:,:), err(:,:) ! Obs and associated error +real(r8), allocatable :: iqc(:,:) ! Incoming data QC + +!------------------------------------------------------------------------ +! Declare namelist parameters +character(len=256) :: file_list = 'sst_file_list.txt' +character(len=256) :: file_out = 'obs_seq.sst' +integer :: avg_obs_per_file = 500000 +logical :: debug = .true. + + +namelist /cmems_sst_to_obs_nml/ file_list, & + file_out, & + avg_obs_per_file, & + debug + +! Start Converter +call initialize_utilities() + +! Read the namelist options +call find_namelist_in_file('input.nml', 'cmems_sst_to_obs_nml', iunit) +read(iunit, nml = cmems_sst_to_obs_nml, iostat = io) +call check_namelist_read(iunit, io, 'cmems_sst_to_obs_nml') + +if (do_nml_file()) write(nmlfileunit, nml=cmems_sst_to_obs_nml) +if (do_nml_term()) write( * , nml=cmems_sst_to_obs_nml) + +! Set the calendar kind +call set_calendar_type(GREGORIAN) + +! Get number of observations +call find_textfile_dims(file_list, nfiles) +num_new_obs = avg_obs_per_file * nfiles + +! Initialize obs seq file +call static_init_obs_sequence() +call init_obs(obs, num_copies, num_qc) +call init_obs(prev_obs, num_copies, num_qc) + +if (file_exist(file_out)) then + string1 = 'Output file: '//trim(adjustl(file_out))//' exists. Replacing it ...' + call error_handler(E_MSG, source, string1) +else + string1 = 'Creating "'//trim(adjustl(file_out))//'" file.' + call error_handler(E_MSG, source, string1) +endif + +call init_obs_sequence(obs_seq, num_copies, num_qc, num_new_obs) +call set_copy_meta_data(obs_seq, num_copies, 'observation') +call set_qc_meta_data(obs_seq, num_qc, 'QC') + +! Loop over the obs files +do filenum = 1, nfiles + next_infile = get_next_filename(file_list, filenum) + + if (debug) write(*, '(/, A, i0, X, A)') 'Input file: #', filenum, trim(next_infile) + + ! Collect data from netcdf file + call collect_data(next_infile, obs, prev_obs, prev_time) + call cleanup() +enddo + +call finish_obs() +call finalize_utilities() + +contains + +!--------------------------------------------------------------- +! Read SST and associated metadata and add them to the obs +! sequence file +subroutine collect_data(datafile, obs, prev_obs, prev_time) + +character(len=*), parameter :: routine = 'collect_data' + +character(len=*), intent(in) :: datafile +type(obs_type), intent(inout) :: obs, prev_obs +type(time_type), intent(inout) :: prev_time + +integer :: ncid, nlon, nlat, ntime +integer :: ilon, ilat +real(r8) :: missing_sst + +integer :: year, month, day, hour, minute, second, ios +real(digits12) :: time_since_init +character(len=64) :: datestr +type(time_type) :: base_date +integer(i8) :: big_integer +integer :: some_seconds, some_days + +integer :: osec, oday + +ncid = nc_open_file_readonly(datafile) +ntime = nc_get_dimension_size(ncid, 'time', routine) + +if (ntime /= 1) then + string1 = '"time" dimension > 1 not supported.' + call error_handler(E_ERR, routine, string1, source) +endif + +! Dimensions +nlat = nc_get_dimension_size(ncid, 'latitude', routine) +nlon = nc_get_dimension_size(ncid, 'longitude', routine) + +! Memory +allocate(lat(nlat), lon(nlon)) +allocate(sst(nlon, nlat), iqc(nlon, nlat), err(nlon, nlat)) + +! Location data +call nc_get_variable(ncid, 'latitude', lat, routine) +call nc_get_variable(ncid, 'longitude', lon, routine) +where(lon < 0.0_r8) lon = lon + 360.0_r8 + +! SST +call nc_get_variable(ncid, 'adjusted_sea_surface_temperature', sst, routine) +call nc_get_attribute_from_variable(ncid, 'adjusted_sea_surface_temperature', & + '_FillValue', missing_sst, routine) + +! Incoming QC +! 0: no_data, 1: bad_data, 2: worst_quality, 3: low_quality, +! 4: acceptable_quality, 5: best_quality +call nc_get_variable(ncid, 'quality_level', iqc, routine) + +! Errors +call nc_get_variable(ncid, 'sses_standard_deviation', err, routine) + +do ilat = 1, nlat + do ilon = 1, nlon + if (err(ilon, ilat) /= err(ilon, ilat)) err(ilon, ilat) = OBS_ERROR_SD_MIN + + ! Clamp the errors + if (sst(ilon, ilat) /= missing_sst .and. sst(ilon, ilat) == sst(ilon, ilat)) then + err(ilon, ilat) = max(err(ilon, ilat), OBS_ERROR_SD_MIN) + endif + enddo +enddo + +! Time +call nc_get_variable(ncid, 'time', time_since_init) +call get_time_units(ncid, datestr) + +!call nc_get_attribute_from_variable(ncid, 'time', 'units', datestr, routine) + +read(datestr, '(14x, i4, 5(1x,i2))', iostat=ios) year, month, day, hour, minute, second +if (ios /= 0) then + write(string1, *) 'Unable to read time variable units. Error status was ', ios + call error_handler(E_ERR, routine, string1, source) +endif +big_integer = int(time_since_init, i8) + +base_date = set_date(year, month, day, hour, minute, second) +some_days = big_integer / 86400 +big_integer = big_integer - (some_days * 86400) +some_seconds = int(big_integer, i4) + +obs_time = base_date + set_time(some_seconds, some_days) +if (debug) call print_date(obs_time, str= ' * Reading SST; date ') + +call nc_close_file(ncid, source) + +! Start adding to the sequence +call get_time(obs_time, osec, oday) + +do ilat = 1, nlat + do ilon = 1, nlon + if (sst(ilon, ilat) == missing_sst .or. & + sst(ilon, ilat) /= sst(ilon, ilat)) cycle + if (iqc(ilon, ilat) <= OBS_QC_LOW_QUALITY) cycle + + if (debug) then + write(*, '(4X, 4(A, F10.6))') 'lon: ', lon(ilon), ', lat: ', lat(ilat), & + ', sst: ', sst(ilon, ilat), ', QC: ', iqc(ilon, ilat) + endif + + call create_3d_obs(lat(ilat), lon(ilon), 0.0_r8, VERTISSURFACE, & + sst(ilon, ilat)-t_kelvin, SATELLITE_BLENDED_SST, & + err(ilon, ilat), oday, osec, obs_qc, obs) + call add_obs_to_seq(obs_seq, obs, obs_time, prev_obs, prev_time, first_obs) + enddo +enddo + +end subroutine collect_data + +!------------------------------------------------------ +! Get the reference time. This could be using +! NF90_STRING (and not NF90_CHAR). If that's the +! case, just fall back to what we expect the reference +! to be. +subroutine get_time_units(ncid, datestr) + +integer, intent(in) :: ncid +character(len=*), intent(out) :: datestr + +integer :: status, varid, xtype, attlen + +! fallback default +datestr = 'seconds since 1970-01-01 00:00:00' + +status = nf90_inq_varid(ncid, 'time', varid) +if (status /= nf90_noerr) return + +status = nf90_inquire_attribute(ncid, varid, 'units', xtype=xtype, len=attlen) +if (status /= nf90_noerr) return + +if (xtype == nf90_char) then + status = nf90_get_att(ncid, varid, 'units', datestr) + if (status /= nf90_noerr) then + ! keep fallback + endif +endif + +end subroutine get_time_units + + +!------------------------------------------------------------ +! All collected obs (if any) are written in the seq file. +subroutine finish_obs() + +integer :: obs_num + +obs_num = get_num_obs(obs_seq) + +if (obs_num > 0) then + if (debug) write(*, '(/, A, i0, A)') '>>>> Ready to write ', obs_num, ' observations:' + + call write_obs_seq(obs_seq, file_out) + call destroy_obs(obs) + call destroy_obs_sequence(obs_seq) +else + string1 = 'No obs were converted.' + call error_handler(E_MSG, source, string1) +endif + +call error_handler(E_MSG, source, 'Finished successfully.') + +end subroutine finish_obs + +!------------------------------------------------------------ +! Clear up memory +subroutine cleanup() + +deallocate(lon, lat, sst, iqc, err) + +end subroutine cleanup + +end program cmems_sst_to_obs diff --git a/observations/obs_converters/cmems_sst_l3s/readme.rst b/observations/obs_converters/cmems_sst_l3s/readme.rst new file mode 100644 index 000000000..1a984c6fe --- /dev/null +++ b/observations/obs_converters/cmems_sst_l3s/readme.rst @@ -0,0 +1,113 @@ +.. _cmems_sst_to_obs: + +CMEMS L3S SST observations +========================== + +.. contents:: + :depth: 3 + :local: + +Overview +-------- +This converter reads Sea Surface Temperature (SST) observations from the +Copernicus Marine Environment Monitoring Service (CMEMS) Level-3S (L3S) +blended SST product and converts them into a DART ``obs_seq`` file. + +The input data consist of bias-corrected, multi-sensor SST observations +from the ODYSSEA system. These products are distributed in netCDF format +and include per-pixel uncertainty estimates and quality control flags. + +The converter extracts: + +- SST observations (``adjusted_sea_surface_temperature``) +- Observation uncertainties (``sses_standard_deviation``) +- Quality flags (``quality_level``) +- Observation time and location + +and writes them into a DART observation sequence file using the +``SATELLITE_BLENDED_SST`` observation type. The unit of +temperature is C. + +For more information about the SST product, visit the +Copernicus Marine Service: https://marine.copernicus.eu + +Input data +---------- +The converter expects netCDF files from the CMEMS SST L3S product, e.g., ``SST_GLO_PHY_L3S_*.nc`` + +Required variables: + +- ``adjusted_sea_surface_temperature`` (Kelvin) +- ``sses_standard_deviation`` (Kelvin) +- ``quality_level`` (0–5) +- ``latitude`` +- ``longitude`` +- ``time`` + +Only observations with quality level greater than 3 are retained. + +Obs Error handling +------------------ +The input uncertainty field ``sses_standard_deviation`` is used as the +observation error. + +To ensure numerical stability and prevent unrealistically small errors +from dominating the assimilation, a lower bound is applied: + +:: + + err = max(sses_standard_deviation, OBS_ERROR_SD_MIN) + +with default: + +:: + + OBS_ERROR_SD_MIN = 0.05 °C + +NaN values in the uncertainty field are replaced with this minimum value. +This lower bound represents a practical safeguard against over-weighting +individual observations and does not replace more sophisticated error +modeling or inflation during assimilation. + +Namelist +-------- + +The converter uses the following namelist: + +:: + + &cmems_sst_to_obs_nml + file_list = 'sst_file_list.txt' ! text file listing input netCDF files + file_out = 'obs_seq.sst' ! output obs_seq filename + avg_obs_per_file = 500000 ! pre-allocation hint + debug = .true. ! print diagnostic output + / + +Usage +----- +1. Create a file list: + +:: + + ls *.nc > sst_file_list.txt + +2. Edit ``input.nml`` to configure options + +3. Build the converter: + +:: + + cd work + ./quickbuild.sh + +4. Run: + +:: + + ./cmems_sst_to_obs + +5. Output: + +:: + + obs_seq.sst diff --git a/observations/obs_converters/cmems_sst_l3s/work/input.nml b/observations/obs_converters/cmems_sst_l3s/work/input.nml new file mode 100644 index 000000000..67a3076ef --- /dev/null +++ b/observations/obs_converters/cmems_sst_l3s/work/input.nml @@ -0,0 +1,40 @@ +&preprocess_nml + overwrite_output = .true. + input_obs_qty_mod_file = '../../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../../assimilation_code/modules/observations/obs_kind_mod.f90' + input_obs_def_mod_file = '../../../forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../forward_operators/obs_def_mod.f90' + obs_type_files = '../../../forward_operators/obs_def_ocean_mod.f90' + quantity_files = '../../../../assimilation_code/modules/observations/ocean_quantities_mod.f90' + / + +&cmems_sst_to_obs_nml + file_list = 'sst_file_list.txt' + file_out = 'obs_seq.sst' + debug = .true. + / + +&utilities_nml + module_details = .false. + / + +&obs_kind_nml + / + +&location_nml + / + +&obs_sequence_nml + write_binary_obs_sequence = .false. + / + +&obs_sequence_tool_nml + filename_seq = '' + filename_out = '' + print_only = .false. + gregorian_cal = .true. + first_obs_days = -1 + first_obs_seconds = -1 + last_obs_days = -1 + last_obs_seconds = -1 + / diff --git a/observations/obs_converters/cmems_sst_l3s/work/quickbuild.sh b/observations/obs_converters/cmems_sst_l3s/work/quickbuild.sh new file mode 100755 index 000000000..91e4fe9be --- /dev/null +++ b/observations/obs_converters/cmems_sst_l3s/work/quickbuild.sh @@ -0,0 +1,39 @@ +#!/usr/bin/env bash + +# DART software - Copyright UCAR. This open source software is provided +# by UCAR, "as is", without charge, subject to all terms of use at +# http://www.image.ucar.edu/DAReS/DART/DART_download + +main() { + +export DART=$(git rev-parse --show-toplevel) +source "$DART"/build_templates/buildconvfunctions.sh + +CONVERTER=cmems_sst_l3s +LOCATION=threed_sphere + +programs=( +advance_time +cmems_sst_to_obs +obs_seq_to_netcdf +obs_sequence_tool +) + +# build arguments +arguments "$@" + +# clean the directory +\rm -f -- *.o *.mod Makefile .cppdefs + +# build and run preprocess before making any other DART executables +buildpreprocess + +# build +buildconv + +# clean up +\rm -f -- *.o *.mod + +} + +main "$@"