Skip to content

Commit

Permalink
Merge branch 'develop' into feature/rreichle/GEOSgcm_GridComp_develop
Browse files Browse the repository at this point in the history
  • Loading branch information
gmao-rreichle committed Mar 19, 2024
2 parents a1b5a3b + cea441d commit fe99646
Show file tree
Hide file tree
Showing 4 changed files with 252 additions and 4 deletions.
4 changes: 3 additions & 1 deletion src/Applications/LDAS_App/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -45,3 +45,5 @@ install(
FILES ${rc_files} ${nml_files} lenkf.j.template
DESTINATION etc
)

esma_add_subdirectories(util/inputs/ASCAT_sm_mask)
6 changes: 3 additions & 3 deletions src/Applications/LDAS_App/LDASsa_DEFAULT_inputs_ensupd.nml
Original file line number Diff line number Diff line change
Expand Up @@ -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 = ''
Expand Down Expand Up @@ -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 = ''
Expand Down Expand Up @@ -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 = ''
Expand Down
Original file line number Diff line number Diff line change
@@ -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)
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit fe99646

Please sign in to comment.