From 4c3062bd4d6d81076fed280b96e39ba233943c7c Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 11 Mar 2024 14:56:57 -0600 Subject: [PATCH 01/12] add ascat mask maker --- src/Applications/LDAS_App/CMakeLists.txt | 6 + .../inputs/ascat_mask/ascat_mask_maker.F90 | 168 ++++++++++++++++++ 2 files changed, 174 insertions(+) create mode 100644 src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index 2d0b55b1..1c5b4dc9 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -13,6 +13,12 @@ ecbuild_add_executable ( TARGET tile_bin2nc4.x SOURCES tile_bin2nc4.F90 LIBS MAPL) + +ecbuild_add_executable ( + TARGET ascat_mask_maker.x + SOURCES util/inputs/ascat_mask/ascat_mask_maker.F90 + LIBS MAPL +) ecbuild_add_executable ( TARGET mwrtm_bin2nc4.x diff --git a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 new file mode 100644 index 00000000..5b348de2 --- /dev/null +++ b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 @@ -0,0 +1,168 @@ +program ascat_mask_maker + + use netcdf + + implicit none + + integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2) + integer :: i, j, closest_index + + integer(kind=1), dimension(:,:), allocatable :: mask_out + integer(kind=1) :: missing_value + + real, dimension(:), allocatable :: asc_lon, asc_lat, cold_mask, wet_mask, veg_mask, subsurface_mask, combined_mask + real, dimension(:), allocatable :: lon, lat, distances + real :: d_lon, d_lat, ll_lon, ll_lat + + ! Open the NetCDF file + ierr = nf90_open('/Users/amfox/Desktop/GEOSldas_diagnostics/test_data/clsm/subsurface_scattering_ASCAT_ERA5_Land.nc', nf90_nowrite, ncid) + if (ierr /= nf90_noerr) stop 'Error opening file' + + ! 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)) + + ! 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 + where (wet_mask == 1) + combined_mask = 1 + elsewhere + combined_mask = subsurface_mask + end where + + d_lon = 0.1 + d_lat = 0.1 + ll_lon = -180.0 + ll_lat = -90.0 + missing_value = -128 + + allocate(lon(int(360.0 / d_lon))) + allocate(lat(int(180.0 / d_lat))) + + lon = [(ll_lon + i * d_lon, i = 0, size(lon) - 1)] + lat = [(ll_lat + i * d_lat, i = 0, size(lat) - 1)] + + allocate(mask_out(size(lon), size(lat))) + + do i = 1, size(lon) + print*, lon(i) + do j = 1, size(lat) + allocate(distances(size(asc_lon))) + distances = sqrt((asc_lon - lon(i))**2 + (asc_lat - lat(j))**2) + closest_index = minloc(distances, dim = 1) + if (distances(closest_index) > 0.14) then + mask_out(i, j) = missing_value + else + mask_out(i, j) = combined_mask(closest_index) + end if + deallocate(distances) + end do + end do + + ! Write out the mask to netcdf + ierr = nf90_create('ascat_combined_mask_p1_f90.nc', nf90_clobber, ncid) + if (ierr /= nf90_noerr) stop 'Error creating file' + + ! Define the dimensions + ierr = nf90_def_dim(ncid, 'lon', size(lon), dimids(1)) + ierr = nf90_def_dim(ncid, 'lat', size(lat), dimids(2)) + + ! 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', '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', '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', 'subsurface_mask') + ierr = nf90_put_att(ncid, varid, 'long_name', 'Mask accounting for subsurface scattering') + 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 \ No newline at end of file From 53065bba3a558b15bc6d4046345625cbe68d14d3 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 12 Mar 2024 13:49:43 -0400 Subject: [PATCH 02/12] discover filepath --- .../LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 index 5b348de2..6034206d 100644 --- a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 +++ b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 @@ -15,7 +15,7 @@ program ascat_mask_maker real :: d_lon, d_lat, ll_lon, ll_lat ! Open the NetCDF file - ierr = nf90_open('/Users/amfox/Desktop/GEOSldas_diagnostics/test_data/clsm/subsurface_scattering_ASCAT_ERA5_Land.nc', nf90_nowrite, ncid) + ierr = nf90_open('/discover/nobackup/amfox/subsurface_scattering_ASCAT_ERA5_Land.nc', nf90_nowrite, ncid) if (ierr /= nf90_noerr) stop 'Error opening file' ! Get the dimension ID @@ -165,4 +165,4 @@ program ascat_mask_maker ierr = nf90_close(ncid) if (ierr /= nf90_noerr) stop 'Error closing file' -end program ascat_mask_maker \ No newline at end of file +end program ascat_mask_maker From ff241e71c0e8105b54bb82469f2d79e6ac8fa084 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 12 Mar 2024 14:04:44 -0600 Subject: [PATCH 03/12] added header comment --- .../util/inputs/ascat_mask/ascat_mask_maker.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 index 6034206d..a3f1944f 100644 --- a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 +++ b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 @@ -1,3 +1,14 @@ +! This program reads in a NetCDF file containing ASCAT soil moisture masks available 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 + +! It provides the possibility to combine different masks (default case is combination of subsurface and wetland masks) +! and interpolates onto a regular grid with a (hardwired) 0.1 degree lat/lon resolution 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 From ae97130435b698c5c163993f4e2422ba874e1cf8 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 12 Mar 2024 14:07:36 -0600 Subject: [PATCH 04/12] change default file write to current directory --- .../LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 index a3f1944f..06570135 100644 --- a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 +++ b/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 @@ -103,7 +103,7 @@ program ascat_mask_maker end do ! Write out the mask to netcdf - ierr = nf90_create('ascat_combined_mask_p1_f90.nc', nf90_clobber, ncid) + ierr = nf90_create('ascat_combined_mask_p1.nc', nf90_clobber, ncid) if (ierr /= nf90_noerr) stop 'Error creating file' ! Define the dimensions From 199c730d1e35acda8617a0510476fedc23d77d03 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 12 Mar 2024 17:50:43 -0400 Subject: [PATCH 05/12] minor cleanup and documentation (ascat_mask_maker.F90, CMakeLists.txt) --- src/Applications/LDAS_App/CMakeLists.txt | 2 +- .../ascat_mask_maker.F90 | 60 ++++++++++++------- 2 files changed, 41 insertions(+), 21 deletions(-) rename src/Applications/LDAS_App/util/inputs/{ascat_mask => ASCAT_sm_mask}/ascat_mask_maker.F90 (82%) diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index 1c5b4dc9..9f43cdd9 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -16,7 +16,7 @@ ecbuild_add_executable ( ecbuild_add_executable ( TARGET ascat_mask_maker.x - SOURCES util/inputs/ascat_mask/ascat_mask_maker.F90 + SOURCES util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 LIBS MAPL ) diff --git a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 similarity index 82% rename from src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 rename to src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 index 06570135..b1849a7d 100644 --- a/src/Applications/LDAS_App/util/inputs/ascat_mask/ascat_mask_maker.F90 +++ b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 @@ -2,11 +2,11 @@ ! 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 - +! ! It provides the possibility to combine different masks (default case is combination of subsurface and wetland masks) -! and interpolates onto a regular grid with a (hardwired) 0.1 degree lat/lon resolution and -90/-180 degree lower left +! 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 @@ -25,9 +25,33 @@ program ascat_mask_maker real, dimension(:), allocatable :: lon, lat, distances real :: d_lon, d_lat, ll_lon, ll_lat + character(400) :: fname_in + + ! -------------------------------------------------------------------------------- + ! + ! hardwired variables + + ! ASCAT soil moisture mask file from Lindorfer et al 2023 + + fname_in = '/discover/nobackup/amfox/subsurface_scattering_ASCAT_ERA5_Land.nc' + + ! Specification of output grid and missing value + + d_lon = 0.1 + d_lat = 0.1 + ll_lon = -180.0 + ll_lat = -90.0 + + missing_value = -128 + + ! ------------------------------- + ! Open the NetCDF file - ierr = nf90_open('/discover/nobackup/amfox/subsurface_scattering_ASCAT_ERA5_Land.nc', nf90_nowrite, ncid) + ierr = nf90_open(fname_in, nf90_nowrite, ncid) if (ierr /= nf90_noerr) stop 'Error opening file' + + ! 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) @@ -40,13 +64,13 @@ program ascat_mask_maker 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(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(combined_mask( N_gpi)) ! Get the variable IDs and read the variables ierr = nf90_inq_varid(ncid, 'lon', varid) @@ -66,18 +90,14 @@ program ascat_mask_maker ierr = nf90_close(ncid) if (ierr /= nf90_noerr) stop 'Error closing file' - ! Combine the masks + ! Combine the masks (1-dim arrays) where (wet_mask == 1) combined_mask = 1 elsewhere combined_mask = subsurface_mask end where - d_lon = 0.1 - d_lat = 0.1 - ll_lon = -180.0 - ll_lat = -90.0 - missing_value = -128 + ! Re-map "combined_mask" from WARP5 input grid to regular lat/lon output grid (2-dim array) allocate(lon(int(360.0 / d_lon))) allocate(lat(int(180.0 / d_lat))) @@ -91,9 +111,9 @@ program ascat_mask_maker print*, lon(i) do j = 1, size(lat) allocate(distances(size(asc_lon))) - distances = sqrt((asc_lon - lon(i))**2 + (asc_lat - lat(j))**2) + distances = (asc_lon - lon(i))**2 + (asc_lat - lat(j))**2 closest_index = minloc(distances, dim = 1) - if (distances(closest_index) > 0.14) then + if (distances(closest_index) > 0.14**2) then mask_out(i, j) = missing_value else mask_out(i, j) = combined_mask(closest_index) @@ -124,8 +144,8 @@ program ascat_mask_maker 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', 'subsurface_mask') - ierr = nf90_put_att(ncid, varid, 'long_name', 'Mask accounting for subsurface scattering') + 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) From 9aa8ec7b2c8c812d404cde1e878a9408dfbe32c9 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 13 Mar 2024 09:26:44 -0600 Subject: [PATCH 06/12] declare mask variables as integers --- .../util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) 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 index b1849a7d..4d909411 100644 --- 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 @@ -15,13 +15,14 @@ program ascat_mask_maker implicit none - integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2) - integer :: i, j, closest_index + integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2) + integer :: i, j, closest_index + 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, cold_mask, wet_mask, veg_mask, subsurface_mask, combined_mask + real, dimension(:), allocatable :: asc_lon, asc_lat real, dimension(:), allocatable :: lon, lat, distances real :: d_lon, d_lat, ll_lon, ll_lat From a1b5a3b9bd9ca787924ce41c66f45d8552e83e2f Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Thu, 14 Mar 2024 13:16:45 -0400 Subject: [PATCH 07/12] reverting to "develop" branch for GEOSgcm_GridComp (components.yaml) --- components.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 409cd2cb668477b57ce8df96dfcd4a2a1bf313d2 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Sun, 17 Mar 2024 15:05:28 -0400 Subject: [PATCH 08/12] move lat/lon to grid cell center --- .../LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 index 4d909411..3f0e3c60 100644 --- 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 @@ -103,8 +103,8 @@ program ascat_mask_maker allocate(lon(int(360.0 / d_lon))) allocate(lat(int(180.0 / d_lat))) - lon = [(ll_lon + i * d_lon, i = 0, size(lon) - 1)] - lat = [(ll_lat + i * d_lat, i = 0, size(lat) - 1)] + lon = [((ll_lon + (d_lon / 2)) + i * d_lon, i = 0, size(lon) - 1)] + lat = [((ll_lat + (d_lat / 2)) + i * d_lat, i = 0, size(lat) - 1)] allocate(mask_out(size(lon), size(lat))) From e0f3655c4de9e80aded7659b7649a14053ff8934 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 18 Mar 2024 16:15:58 -0600 Subject: [PATCH 09/12] add combination choices --- .../inputs/ASCAT_sm_mask/ascat_mask_maker.F90 | 61 +++++++++++++++---- 1 file changed, 49 insertions(+), 12 deletions(-) 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 index 3f0e3c60..a67e1803 100644 --- 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 @@ -16,7 +16,7 @@ program ascat_mask_maker implicit none integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2) - integer :: i, j, closest_index + 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 @@ -26,7 +26,7 @@ program ascat_mask_maker real, dimension(:), allocatable :: lon, lat, distances real :: d_lon, d_lat, ll_lon, ll_lat - character(400) :: fname_in + character(200) :: fname_in, mask_description, fname_out ! -------------------------------------------------------------------------------- ! @@ -34,8 +34,16 @@ program ascat_mask_maker ! ASCAT soil moisture mask file from Lindorfer et al 2023 - fname_in = '/discover/nobackup/amfox/subsurface_scattering_ASCAT_ERA5_Land.nc' + fname_in = '/discover/nobackup/amfox/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 @@ -44,6 +52,9 @@ program ascat_mask_maker ll_lat = -90.0 missing_value = -128 + + ! Specify the NetCDF file name for the output mask + fname_out = 'ascat_combined_mask_p1.nc' ! ------------------------------- @@ -92,18 +103,39 @@ program ascat_mask_maker if (ierr /= nf90_noerr) stop 'Error closing file' ! Combine the masks (1-dim arrays) - where (wet_mask == 1) - combined_mask = 1 - elsewhere - combined_mask = subsurface_mask - end where + select case (mask_mode) + case (1) + ! Combine wet_mask and subsurface_mask + mask_description = 'Combined subsurface and wetland mask' + where (wet_mask == 1) + 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 == 1 .or. veg_mask == 1) + 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) allocate(lon(int(360.0 / d_lon))) allocate(lat(int(180.0 / d_lat))) - lon = [((ll_lon + (d_lon / 2)) + i * d_lon, i = 0, size(lon) - 1)] + lon = [((ll_lon + (d_lon / 2)) + i * d_lon, i = 0, size(lon) - 1)] ! NB using grid cell centers for nearest neighbor search lat = [((ll_lat + (d_lat / 2)) + i * d_lat, i = 0, size(lat) - 1)] allocate(mask_out(size(lon), size(lat))) @@ -124,23 +156,28 @@ program ascat_mask_maker end do ! Write out the mask to netcdf - ierr = nf90_create('ascat_combined_mask_p1.nc', nf90_clobber, ncid) + ierr = nf90_create(fname_out, nf90_clobber, ncid) if (ierr /= nf90_noerr) stop 'Error creating file' ! Define the dimensions ierr = nf90_def_dim(ncid, 'lon', size(lon), dimids(1)) ierr = nf90_def_dim(ncid, 'lat', size(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') + 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', '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', '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') From 5805a55f9b74e299880bd7deff97fddfa958a3c7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 19 Mar 2024 11:44:59 -0400 Subject: [PATCH 10/12] additional cleanup and documentation of ascat_mask_maker.F90: - clarified comments and error messages - change "== 1" to "/= 0" for consistency with ASCAT EUMETSAT reader - avoid repeated allocate/deallocate within i,j loop - use nint() to avoid risky integer division when calculating N_lon, N_lat --- .../inputs/ASCAT_sm_mask/ascat_mask_maker.F90 | 53 ++++++++++--------- 1 file changed, 29 insertions(+), 24 deletions(-) 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 index a67e1803..83760b65 100644 --- 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 @@ -1,9 +1,11 @@ -! This program reads in a NetCDF file containing ASCAT soil moisture masks available from: +! 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 ! -! It provides the possibility to combine different masks (default case is combination of subsurface and wetland masks) +! 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. ! @@ -15,7 +17,7 @@ program ascat_mask_maker implicit none - integer :: ncid, varid, dimid, ierr, len, N_gpi, dimids(2) + 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 @@ -48,19 +50,20 @@ program ascat_mask_maker d_lon = 0.1 d_lat = 0.1 - ll_lon = -180.0 - ll_lat = -90.0 + 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 file + ! Open the NetCDF input file ierr = nf90_open(fname_in, nf90_nowrite, ncid) - if (ierr /= nf90_noerr) stop 'Error opening file' + 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). @@ -83,6 +86,7 @@ program ascat_mask_maker 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) @@ -107,7 +111,7 @@ program ascat_mask_maker case (1) ! Combine wet_mask and subsurface_mask mask_description = 'Combined subsurface and wetland mask' - where (wet_mask == 1) + where (wet_mask /= 0) combined_mask = 1 elsewhere combined_mask = subsurface_mask @@ -123,7 +127,7 @@ program ascat_mask_maker case (4) ! Combine subsurface_mask, wet_mask, and veg_mask mask_description = 'Combined subsurface, wetland, and vegetation mask' - where (wet_mask == 1 .or. veg_mask == 1) + where (wet_mask /= 0 .or. veg_mask /= 0) combined_mask = 1 elsewhere combined_mask = subsurface_mask @@ -132,40 +136,41 @@ program ascat_mask_maker ! Re-map "combined_mask" from WARP5 input grid to regular lat/lon output grid (2-dim array) - allocate(lon(int(360.0 / d_lon))) - allocate(lat(int(180.0 / d_lat))) + 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, size(lon) - 1)] ! NB using grid cell centers for nearest neighbor search - lat = [((ll_lat + (d_lat / 2)) + i * d_lat, i = 0, size(lat) - 1)] + 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(size(lon), size(lat))) + allocate(mask_out( N_lon, N_lat)) - do i = 1, size(lon) + do i = 1, N_lon print*, lon(i) - do j = 1, size(lat) - allocate(distances(size(asc_lon))) + 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 + 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 - deallocate(distances) 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' + if (ierr /= nf90_noerr) stop 'Error creating file: ' // fname_out ! Define the dimensions - ierr = nf90_def_dim(ncid, 'lon', size(lon), dimids(1)) - ierr = nf90_def_dim(ncid, 'lat', size(lat), dimids(2)) + 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') + 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 @@ -183,7 +188,7 @@ program ascat_mask_maker 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, 'long_name', 'combined mask') ierr = nf90_put_att(ncid, varid, 'units', 'boolean') ierr = nf90_put_att(ncid, varid, '_FillValue', missing_value) From f7d0750067af8bb263fdfffdb6a6119641ecd87b Mon Sep 17 00:00:00 2001 From: amfox37 Date: Tue, 19 Mar 2024 13:16:11 -0600 Subject: [PATCH 11/12] CMakeLists to build without install --- src/Applications/LDAS_App/CMakeLists.txt | 8 ++------ .../LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt | 4 ++++ 2 files changed, 6 insertions(+), 6 deletions(-) create mode 100644 src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt diff --git a/src/Applications/LDAS_App/CMakeLists.txt b/src/Applications/LDAS_App/CMakeLists.txt index 9f43cdd9..1075229e 100644 --- a/src/Applications/LDAS_App/CMakeLists.txt +++ b/src/Applications/LDAS_App/CMakeLists.txt @@ -14,12 +14,6 @@ ecbuild_add_executable ( SOURCES tile_bin2nc4.F90 LIBS MAPL) -ecbuild_add_executable ( - TARGET ascat_mask_maker.x - SOURCES util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 - LIBS MAPL -) - ecbuild_add_executable ( TARGET mwrtm_bin2nc4.x SOURCES util/inputs/mwRTM_params/mwrtm_bin2nc4.F90 @@ -51,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/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..ac0fbd08 --- /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) \ No newline at end of file From 48c5ece1894ca7a7817c6a924ad628ebaf0f5f79 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 19 Mar 2024 16:31:33 -0400 Subject: [PATCH 12/12] updated paths to ASCAT EUMETSAT soil moisture retrieval data and mask file (LDASsa_DEFAULT_inputs_ensupd.nml, ascat_mask_maker.F90) --- src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml | 6 +++--- .../LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt | 2 +- .../LDAS_App/util/inputs/ASCAT_sm_mask/ascat_mask_maker.F90 | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) 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 index ac0fbd08..11114b07 100644 --- a/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt +++ b/src/Applications/LDAS_App/util/inputs/ASCAT_sm_mask/CMakeLists.txt @@ -1,4 +1,4 @@ # build without installation add_executable(ascat_mask_maker.x ascat_mask_maker.F90) -target_link_libraries(ascat_mask_maker.x MAPL) \ No newline at end of file +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 index 83760b65..e9b4bbd0 100644 --- 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 @@ -36,7 +36,7 @@ program ascat_mask_maker ! ASCAT soil moisture mask file from Lindorfer et al 2023 - fname_in = '/discover/nobackup/amfox/subsurface_scattering_ASCAT_ERA5_Land.nc' + 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