diff --git a/components.yaml b/components.yaml index 9ccb6dad..b5fee251 100644 --- a/components.yaml +++ b/components.yaml @@ -49,6 +49,6 @@ MAPL: GEOSgcm_GridComp: local: ./src/Components/GEOSldas_GridComp/@GEOSgcm_GridComp remote: ../GEOSgcm_GridComp.git - branch: bugfix/rreichle/snow_excs + branch: develop sparse: ./config/GEOSgcm_GridComp_ldas.sparse develop: develop diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index 2d0b55b1..1075229e 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -13,7 +13,7 @@ ecbuild_add_executable ( TARGET tile_bin2nc4.x SOURCES tile_bin2nc4.F90 LIBS MAPL) - + ecbuild_add_executable ( TARGET mwrtm_bin2nc4.x SOURCES util/inputs/mwRTM_params/mwrtm_bin2nc4.F90 @@ -45,3 +45,5 @@ install( FILES ${rc_files} ${nml_files} lenkf.j.template DESTINATION etc ) + +esma_add_subdirectories(util/inputs/ASCAT_sm_mask) diff --git a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml index 0234887a..ece206c4 100644 --- a/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2177,7 +2177,7 @@ obs_param_nml(49)%bias_tcut = 432000 obs_param_nml(49)%nodata = -9999. obs_param_nml(49)%varname = 'sfds' obs_param_nml(49)%units = '%' -obs_param_nml(49)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_A/' +obs_param_nml(49)%path = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/ASCAT_EUMETSAT/Metop_A/' obs_param_nml(49)%name = 'M02-ASCA-ASCSMO02' obs_param_nml(49)%maskpath = '' obs_param_nml(49)%maskname = '' @@ -2216,7 +2216,7 @@ obs_param_nml(50)%bias_tcut = 432000 obs_param_nml(50)%nodata = -9999. obs_param_nml(50)%varname = 'sfds' obs_param_nml(50)%units = '%' -obs_param_nml(50)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_B/' +obs_param_nml(50)%path = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/ASCAT_EUMETSAT/Metop_B/' obs_param_nml(50)%name = 'M01-ASCA-ASCSMO02' obs_param_nml(50)%maskpath = '' obs_param_nml(50)%maskname = '' @@ -2255,7 +2255,7 @@ obs_param_nml(51)%bias_tcut = 432000 obs_param_nml(51)%nodata = -9999. obs_param_nml(51)%varname = 'sfds' obs_param_nml(51)%units = '%' -obs_param_nml(51)%path = '/discover/nobackup/projects/gmao/merra/iau/merra_land/DATA/ASCAT_EUMETSAT/Metop_C/' +obs_param_nml(51)%path = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/ASCAT_EUMETSAT/Metop_C/' obs_param_nml(51)%name = 'M03-ASCA-ASCSMO02' obs_param_nml(51)%maskpath = '' obs_param_nml(51)%maskname = '' diff --git a/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt new file mode 100644 index 00000000..11114b07 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt @@ -0,0 +1,4 @@ +# build without installation + +add_executable(ascat_mask_maker.x ascat_mask_maker.F90) +target_link_libraries(ascat_mask_maker.x MAPL) diff --git a/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 new file mode 100644 index 00000000..e9b4bbd0 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 @@ -0,0 +1,242 @@ +! This program produces a combined mask for use with the assimilation of ASCAT soil moisture retrievals in GEOSldas. +! The combined mask is based on component masks from: +! +! Lindorfer, R., Wagner, W., Hahn, S., Kim, H., Vreugdenhil, M., Gruber, A., Fischer, M., & Trnka, M. (2023). +! Global Scale Maps of Subsurface Scattering Signals Impacting ASCAT Soil Moisture Retrievals (1.0.0) [Data set]. +! TU Wien. https://doi.org/10.48436/9a2y9-e5z14 +! +! The program provides the possibility to combine different masks (default is combination of subsurface and wetland masks) +! and interpolates onto a regular grid with a (hardwired) 0.1 degree lat/lon spacing and -90/-180 degree lower left +! corner used for quick indexing in the ASCAT observation reader QC routine. +! +! Author: AM Fox, March, 2024 + +program ascat_mask_maker + + use netcdf + + implicit none + + integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2), N_lon, N_lat + integer :: i, j, closest_index, mask_mode + integer, dimension(:), allocatable :: cold_mask, wet_mask, veg_mask, subsurface_mask, combined_mask + + integer(kind=1), dimension(:,:), allocatable :: mask_out + integer(kind=1) :: missing_value + + real, dimension(:), allocatable :: asc_lon, asc_lat + real, dimension(:), allocatable :: lon, lat, distances + real :: d_lon, d_lat, ll_lon, ll_lat + + character(200) :: fname_in, mask_description, fname_out + + ! -------------------------------------------------------------------------------- + ! + ! hardwired variables + + ! ASCAT soil moisture mask file from Lindorfer et al 2023 + + fname_in = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/ASCAT_EUMETSAT/Mask/subsurface_scattering_ASCAT_ERA5_Land.nc' + + ! Specification of how to combine the masks + ! Mask_mode = 1 (default) combines subsurface and wetland masks + ! Mask_mode = 2 uses only the subsurface mask + ! Mask_mode = 3 uses only the wetland mask + ! Mask_mode = 4 combines subsurface, wetland and vegetation masks + + mask_mode = 1 + + ! Specification of output grid and missing value + + d_lon = 0.1 + d_lat = 0.1 + ll_lon = -180.0 ! longitude of boundary of lower-left grid cell (longitude of lower left corner of grid) + ll_lat = -90.0 ! latitude of boundary of lower-left grid cell (latitude of lower left corner of grid) + + missing_value = -128 + + ! Specify the NetCDF file name for the output mask + + fname_out = 'ascat_combined_mask_p1.nc' + + ! ------------------------------- + + ! Open the NetCDF input file + ierr = nf90_open(fname_in, nf90_nowrite, ncid) + if (ierr /= nf90_noerr) stop 'Error opening file: ' // trim(fname_in) + + ! Data in original mask file are on the 12.5 km fixed Earth grid used for ASCAT (WARP5 grid) and + ! stored in the NetCDF file as 1-dimensional arrays of length N_gpi (over land only). + + ! Get the dimension ID + ierr = nf90_inq_dimid(ncid, 'gpi', dimid) + if (ierr /= nf90_noerr) stop 'Error getting dimension ID' + + ! Get the length of the dimension + ierr = nf90_inquire_dimension(ncid, dimid, len = N_gpi) + if (ierr /= nf90_noerr) stop 'Error inquiring dimension' + + print*, 'N_gpi = ', N_gpi + + ! Allocate the arrays + allocate(asc_lon( N_gpi)) + allocate(asc_lat( N_gpi)) + allocate(cold_mask( N_gpi)) + allocate(wet_mask( N_gpi)) + allocate(veg_mask( N_gpi)) + allocate(subsurface_mask(N_gpi)) + allocate(combined_mask( N_gpi)) + allocate(distances( N_gpi)) + + ! Get the variable IDs and read the variables + ierr = nf90_inq_varid(ncid, 'lon', varid) + ierr = nf90_get_var(ncid, varid, asc_lon) + ierr = nf90_inq_varid(ncid, 'lat', varid) + ierr = nf90_get_var(ncid, varid, asc_lat) + ierr = nf90_inq_varid(ncid, 'cold_mask', varid) + ierr = nf90_get_var(ncid, varid, cold_mask) + ierr = nf90_inq_varid(ncid, 'wet_mask', varid) + ierr = nf90_get_var(ncid, varid, wet_mask) + ierr = nf90_inq_varid(ncid, 'veg_mask', varid) + ierr = nf90_get_var(ncid, varid, veg_mask) + ierr = nf90_inq_varid(ncid, 'subsurface_mask', varid) + ierr = nf90_get_var(ncid, varid, subsurface_mask) + + ! Close the NetCDF file + ierr = nf90_close(ncid) + if (ierr /= nf90_noerr) stop 'Error closing file' + + ! Combine the masks (1-dim arrays) + select case (mask_mode) + case (1) + ! Combine wet_mask and subsurface_mask + mask_description = 'Combined subsurface and wetland mask' + where (wet_mask /= 0) + combined_mask = 1 + elsewhere + combined_mask = subsurface_mask + end where + case (2) + ! Use only subsurface_mask + mask_description = 'Used only subsurface mask' + combined_mask = subsurface_mask + case (3) + ! Use only wet_mask + mask_description = 'Used only wetland mask' + combined_mask = wet_mask + case (4) + ! Combine subsurface_mask, wet_mask, and veg_mask + mask_description = 'Combined subsurface, wetland, and vegetation mask' + where (wet_mask /= 0 .or. veg_mask /= 0) + combined_mask = 1 + elsewhere + combined_mask = subsurface_mask + end where + end select + + ! Re-map "combined_mask" from WARP5 input grid to regular lat/lon output grid (2-dim array) + + N_lon = nint( 360.0 / d_lon ) + N_lat = nint( 180.0 / d_lat ) + + allocate(lon(N_lon)) + allocate(lat(N_lat)) + + lon = [((ll_lon + (d_lon / 2)) + i * d_lon, i = 0, N_lon - 1)] ! NB using grid cell centers for nearest neighbor search + lat = [((ll_lat + (d_lat / 2)) + i * d_lat, i = 0, N_lat - 1)] + + allocate(mask_out( N_lon, N_lat)) + + do i = 1, N_lon + print*, lon(i) + do j = 1, N_lat + distances = (asc_lon - lon(i))**2 + (asc_lat - lat(j))**2 + closest_index = minloc(distances, dim = 1) + if (distances(closest_index) > 0.14**2) then + mask_out(i, j) = missing_value ! Note: ASCAT EUMETSAT reader masks everything .ne. 0 + else + mask_out(i, j) = combined_mask(closest_index) + end if + end do + end do + + ! Write out the mask to netcdf + ierr = nf90_create(fname_out, nf90_clobber, ncid) + if (ierr /= nf90_noerr) stop 'Error creating file: ' // fname_out + + ! Define the dimensions + ierr = nf90_def_dim(ncid, 'lon', N_lon, dimids(1)) + ierr = nf90_def_dim(ncid, 'lat', N_lat, dimids(2)) + + ! Define the global attributes + ierr = nf90_put_att(ncid, nf90_global, 'title', 'ASCAT combined mask') + ierr = nf90_put_att(ncid, nf90_global, 'source', 'Lindorfer et al 2023 doi:10.48436/9a2y9-e5z14') + ierr = nf90_put_att(ncid, nf90_global, 'description', mask_description) + + ! Define the variables + ierr = nf90_def_var(ncid, 'lat', nf90_real, dimids(2), varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'latitude') + ierr = nf90_put_att(ncid, varid, 'long_name', 'grid cell center latitude') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees_north') + ierr = nf90_put_att(ncid, varid, 'axis', 'Y') + + ierr = nf90_def_var(ncid, 'lon', nf90_real, dimids(1), varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'longitude') + ierr = nf90_put_att(ncid, varid, 'long_name', 'grid cell center longitude') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees_east') + ierr = nf90_put_att(ncid, varid, 'axis', 'X') + + ierr = nf90_def_var(ncid, 'mask', nf90_byte, dimids, varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'combined_mask') + ierr = nf90_put_att(ncid, varid, 'long_name', 'combined mask') + ierr = nf90_put_att(ncid, varid, 'units', 'boolean') + ierr = nf90_put_att(ncid, varid, '_FillValue', missing_value) + + ierr = nf90_def_var(ncid, 'll_lon', nf90_real, varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'longitude of lower left corner') + ierr = nf90_put_att(ncid, varid, 'long_name', 'longitude of lower left corner') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees_east') + ierr = nf90_put_att(ncid, varid, 'axis', 'X') + + ierr = nf90_def_var(ncid, 'll_lat', nf90_real, varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'latitude of lower left corner') + ierr = nf90_put_att(ncid, varid, 'long_name', 'latitude of lower left corner') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees_north') + ierr = nf90_put_att(ncid, varid, 'axis', 'Y') + + ierr = nf90_def_var(ncid, 'd_lon', nf90_real, varid) + ierr = nf90_put_att(ncid, varid, 'standard_name', 'longitude grid spacing') + ierr = nf90_put_att(ncid, varid, 'long_name', 'longitude grid spacing') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees') + ierr = nf90_put_att(ncid, varid, 'axis', 'X') + + ierr = nf90_def_var(ncid, 'd_lat', nf90_real, varid) + ierr = nf90_put_att(ncid, varid, 'long_name', 'latitude grid spacing') + ierr = nf90_put_att(ncid, varid, 'units', 'degrees') + ierr = nf90_put_att(ncid, varid, 'axis', 'Y') + + ! End define mode + ierr = nf90_enddef(ncid) + + ! Write the variables + ierr = nf90_inq_varid(ncid, 'lat', varid) + ierr = nf90_put_var(ncid, varid, lat) + ierr = nf90_inq_varid(ncid, 'lon', varid) + ierr = nf90_put_var(ncid, varid, lon) + ierr = nf90_inq_varid(ncid, 'mask', varid) + ierr = nf90_put_var(ncid, varid, mask_out) + if (ierr /= nf90_noerr) stop 'Error writing variable' + ierr = nf90_inq_varid(ncid, 'll_lon', varid) + ierr = nf90_put_var(ncid, varid, ll_lon) + ierr = nf90_inq_varid(ncid, 'll_lat', varid) + ierr = nf90_put_var(ncid, varid, ll_lat) + ierr = nf90_inq_varid(ncid, 'd_lon', varid) + ierr = nf90_put_var(ncid, varid, d_lon) + ierr = nf90_inq_varid(ncid, 'd_lat', varid) + ierr = nf90_put_var(ncid, varid, d_lat) + + ! Close the NetCDF file + ierr = nf90_close(ncid) + if (ierr /= nf90_noerr) stop 'Error closing file' + +end program ascat_mask_maker