From 791b3466021ef30b7da3ac762f89bbed5ab25369 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 8 Jan 2025 17:54:26 -0700 Subject: [PATCH 01/21] read_obs_sm_CYGNSS stub --- .../clsm_ensupd_read_obs.F90 | 115 +++++++++++++++--- 1 file changed, 100 insertions(+), 15 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 05eeaf6..03baa46 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -1830,19 +1830,19 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! open bufr file - call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it - open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') - call openbf(lnbufr, 'SEC3', lnbufr) - call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) - call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) +! call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it +! open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') +! call openbf(lnbufr, 'SEC3', lnbufr) +! call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) +! call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) - msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) +! msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) - loop_report: do while( ireadsb(lnbufr) == 0 ) +! loop_report: do while( ireadsb(lnbufr) == 0 ) ! columns of tmp_data: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 - call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') +! call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') N_obs = N_obs + 1 @@ -1853,11 +1853,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & tmp_data(N_obs,:) = tmp_vdata - end do loop_report - end do msg_report +! end do loop_report +! end do msg_report - call closbf(lnbufr) - close(lnbufr) +! call closbf(lnbufr) +! close(lnbufr) end do ! end file loop @@ -2121,7 +2121,72 @@ subroutine read_obs_sm_ASCAT_EUMET( & if (associated(tmp_jtime)) deallocate(tmp_jtime) end subroutine read_obs_sm_ASCAT_EUMET - + + ! **************************************************************************** + + subroutine read_obs_sm_CYGNSS( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, CYGNSS_sm, CYGNSS_sm_std, CYGNSS_lon, CYGNSS_lat, CYGNSS_time ) + +!--------------------------------------------------------------------- +! +! Routine to read in ASCAT surface degree of saturation (sfds) obs from +! BUFR files for both ascending and descending passes. +! +! ASCAT_sm and ASCAT_sm_std outputs from this subroutine are in wetness fraction (i.e., 0-1) units! +! +! Read in EUMETSAT level 2 soil moisture product 25 km (SMO), PPF software version 5.0. +! Soil moisture derived from re-sampled (spatially averaged) backscatter (sigma0) values +! on a 25-km orbit swath grid. Input data files are in BUFR file format. +! +! EUMETSAT BUFR files contain data for both ascending and descending half-orbits. +! The BUFR field DOMO ("Direction of motion of moving observing platform") could be used to +! separate Asc and Desc. (The BUFR files do not contain any explicit orbit indicator variable.) +! According to Pamela Schoebel-Pattiselanno, EUMETSAT User Services Helpdesk: +! "When the value (of DOMO) is between 180 and 270 degrees, it is the descending part +! of the orbit. When it is between 270 and 360 degrees, it is the ascending part." +! +! Q. Liu, Nov 2019 - based on read_obs_sm_ASCAT +! A. Fox, reichle, Sep 2023 - updated +! A. Fox, reichle, Feb 2024 - added ASCAT obs mask +! +! -------------------------------------------------------------------- + + use netcdf + implicit none + + ! inputs: + + type(date_time_type), intent(in) :: date_time + + integer, intent(in) :: dtstep_assim, N_catd + + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input + + type(grid_def_type), intent(in) :: tile_grid_d + + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & + N_tile_in_cell_ij + + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input + + type(obs_param_type), intent(in) :: this_obs_param + + ! outputs: + + logical, intent(out) :: found_obs + + real, intent(out), dimension(N_catd) :: CYGNSS_sm ! sfds obs [fraction] (i.e., 0-1) + real, intent(out), dimension(N_catd) :: CYGNSS_sm_std ! sfds obs err std [fraction] (i.e., 0-1) + real, intent(out), dimension(N_catd) :: CYGNSS_lon, CYGNSS_lat + real*8, intent(out), dimension(N_catd) :: CYGNSS_time ! J2000 seconds + + ! --------------- + + end subroutine read_obs_sm_CYGNSS + ! *************************************************************************** subroutine read_sm_ASCAT_bin( & @@ -2307,7 +2372,6 @@ subroutine read_sm_ASCAT_bin( & end subroutine read_sm_ASCAT_bin - ! ***************************************************************** subroutine read_obs_LaRC_Tskin( & @@ -8568,7 +8632,28 @@ subroutine read_obs( & date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & tmp_obs, tmp_std_obs ) - end if + end if + + case ('CYGNSS_SM','CYGNSS_SM_daily') + + call read_obs_sm_CYGNSS( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & + tmp_time) + + ! scale observations to model climatology + + if (this_obs_param%scale .and. found_obs) then + + scaled_obs = .true. + + call scale_obs_sfmc_zscore( N_catd, tile_coord, & + date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & + tmp_obs, tmp_std_obs ) + + end if case ('isccp_tskin_gswp2_v1') From 9ffc47ae93740a43b2b8e3bb392f23093dc05837 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Wed, 8 Jan 2025 18:06:33 -0700 Subject: [PATCH 02/21] edit description --- .../clsm_ensupd_read_obs.F90 | 30 +++++-------------- 1 file changed, 7 insertions(+), 23 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 03baa46..3609731 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2130,29 +2130,13 @@ subroutine read_obs_sm_CYGNSS( & this_obs_param, & found_obs, CYGNSS_sm, CYGNSS_sm_std, CYGNSS_lon, CYGNSS_lat, CYGNSS_time ) -!--------------------------------------------------------------------- -! -! Routine to read in ASCAT surface degree of saturation (sfds) obs from -! BUFR files for both ascending and descending passes. -! -! ASCAT_sm and ASCAT_sm_std outputs from this subroutine are in wetness fraction (i.e., 0-1) units! -! -! Read in EUMETSAT level 2 soil moisture product 25 km (SMO), PPF software version 5.0. -! Soil moisture derived from re-sampled (spatially averaged) backscatter (sigma0) values -! on a 25-km orbit swath grid. Input data files are in BUFR file format. -! -! EUMETSAT BUFR files contain data for both ascending and descending half-orbits. -! The BUFR field DOMO ("Direction of motion of moving observing platform") could be used to -! separate Asc and Desc. (The BUFR files do not contain any explicit orbit indicator variable.) -! According to Pamela Schoebel-Pattiselanno, EUMETSAT User Services Helpdesk: -! "When the value (of DOMO) is between 180 and 270 degrees, it is the descending part -! of the orbit. When it is between 270 and 360 degrees, it is the ascending part." -! -! Q. Liu, Nov 2019 - based on read_obs_sm_ASCAT -! A. Fox, reichle, Sep 2023 - updated -! A. Fox, reichle, Feb 2024 - added ASCAT obs mask -! -! -------------------------------------------------------------------- + !--------------------------------------------------------------------- + ! + ! Routine to read in CYGNSS soil moisture obs from netCDF files. + ! + ! A. Fox, Jan 2025 + ! + ! -------------------------------------------------------------------- use netcdf implicit none From 05eb461b5746e90da05222bf311c8acdd1bb76a0 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 9 Jan 2025 12:36:48 -0700 Subject: [PATCH 03/21] add a new type of obs --- .../clsm_ensupd_glob_param.F90 | 2 +- .../clsm_ensupd_read_obs.F90 | 40 ++++++++++++++++++- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 38 ++++++++++++++++++ 3 files changed, 77 insertions(+), 3 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 b/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 index 678a097..6d1ef60 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 @@ -48,7 +48,7 @@ module clsm_ensupd_glob_param ! total number of all obs species defined in "ensupd" namelist file ! (regardless of whether "assim" flag is true or false) - integer, parameter :: N_obs_species_nml = 53 + integer, parameter :: N_obs_species_nml = 54 ! ---------------------------------------------------------------------- ! diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 253528f..d79e279 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2186,13 +2186,49 @@ subroutine read_obs_sm_CYGNSS( & logical, intent(out) :: found_obs - real, intent(out), dimension(N_catd) :: CYGNSS_sm ! sfds obs [fraction] (i.e., 0-1) - real, intent(out), dimension(N_catd) :: CYGNSS_sm_std ! sfds obs err std [fraction] (i.e., 0-1) + real, intent(out), dimension(N_catd) :: CYGNSS_sm ! sm obs [cm3 cm-3] + real, intent(out), dimension(N_catd) :: CYGNSS_sm_std ! sm obs err std [cm3 cm-3] real, intent(out), dimension(N_catd) :: CYGNSS_lon, CYGNSS_lat real*8, intent(out), dimension(N_catd) :: CYGNSS_time ! J2000 seconds ! --------------- + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + + type(date_time_type) :: date_time_obs_beg, date_time_obs_end + type(date_time_type) :: date_time_up, date_time_low + + character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' + character(len=400) :: err_msg + + ! ------------------------------------------------------------------- + + ! initialize + + found_obs = .false. + + ! determine operating time range of sensor + + if (trim(this_obs_param%descr) == 'CYGNSS_SM') then + date_time_obs_beg = date_time_type(2018, 8, 1, 0, 0, 0,-9999,-9999) + date_time_obs_end = date_time_type(2100, 1, 1, 0, 0, 0,-9999,-9999) + else + err_msg = 'Unknown obs_param%descr: ' // trim(this_obs_param%descr) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! return if date_time falls outside operating time range + + if ( datetime_lt_refdatetime(date_time, date_time_obs_beg) .or. & + datetime_lt_refdatetime(date_time_obs_end, date_time) ) return + + ! cygnss sm subdaily observations are for 4 timeslices per day - 00-06, 06-12, 12-18, 18-24 UTC + ! we will assimilate these at 03, 09, 15, 21 UTC + + write (logunit,*) 'This is where we will read CYGNSS obs at YY:', date_time%year, ' DD:', date_time%day, ' HH:', date_time%hour + + return + end subroutine read_obs_sm_CYGNSS ! *************************************************************************** diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index 5085098..c6cf01d 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2363,6 +2363,44 @@ obs_param_nml(53)%xcorr = 0. obs_param_nml(53)%ycorr = 0. obs_param_nml(53)%adapt = 0 +! -------------------------------------------------------------------------------------------------------------- +! +! 54 = CYGNSS_SM (CYGNSS Level 3 soil moisture version 3.2) +! +! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 + +obs_param_nml(49)%descr = 'CYGNSS_SM' +obs_param_nml(49)%orbit = 3 +obs_param_nml(49)%pol = 0 +obs_param_nml(49)%N_ang = 0 +obs_param_nml(49)%freq = 0 +obs_param_nml(49)%FOV = 20. +obs_param_nml(49)%FOV_units = 'km' +obs_param_nml(49)%assim = .false. +obs_param_nml(49)%scale = .false. +obs_param_nml(49)%getinnov = .false. +obs_param_nml(49)%RTM_ID = 0 +obs_param_nml(49)%bias_Npar = 0 +obs_param_nml(49)%bias_trel = 864000 +obs_param_nml(49)%bias_tcut = 432000 +obs_param_nml(49)%nodata = -9999. +obs_param_nml(49)%varname = 'sfmc' +obs_param_nml(49)%units = 'cm3 cm-3' +obs_param_nml(49)%path = '/discover/nobackup/amfox/CYGNSS_obs/' +obs_param_nml(49)%name = '' +obs_param_nml(49)%maskpath = '' +obs_param_nml(49)%maskname = '' +obs_param_nml(49)%scalepath = '' +obs_param_nml(49)%scalename = '' +obs_param_nml(49)%flistpath = '' +obs_param_nml(49)%flistname = '' +obs_param_nml(49)%errstd = 9. +obs_param_nml(49)%std_normal_max = 2.5 +obs_param_nml(49)%zeromean = .true. +obs_param_nml(49)%coarsen_pert = .false. +obs_param_nml(49)%xcorr = 0.25 +obs_param_nml(49)%ycorr = 0.25 +obs_param_nml(49)%adapt = 0 ! -------------------------------------------------------------------- From 7fa8180134306f46187b59fa694424df86ddfa78 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 9 Jan 2025 14:39:37 -0700 Subject: [PATCH 04/21] fix number in ensupd.nml --- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 64 +++++++++---------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index c6cf01d..3724d02 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2369,38 +2369,38 @@ obs_param_nml(53)%adapt = 0 ! ! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 -obs_param_nml(49)%descr = 'CYGNSS_SM' -obs_param_nml(49)%orbit = 3 -obs_param_nml(49)%pol = 0 -obs_param_nml(49)%N_ang = 0 -obs_param_nml(49)%freq = 0 -obs_param_nml(49)%FOV = 20. -obs_param_nml(49)%FOV_units = 'km' -obs_param_nml(49)%assim = .false. -obs_param_nml(49)%scale = .false. -obs_param_nml(49)%getinnov = .false. -obs_param_nml(49)%RTM_ID = 0 -obs_param_nml(49)%bias_Npar = 0 -obs_param_nml(49)%bias_trel = 864000 -obs_param_nml(49)%bias_tcut = 432000 -obs_param_nml(49)%nodata = -9999. -obs_param_nml(49)%varname = 'sfmc' -obs_param_nml(49)%units = 'cm3 cm-3' -obs_param_nml(49)%path = '/discover/nobackup/amfox/CYGNSS_obs/' -obs_param_nml(49)%name = '' -obs_param_nml(49)%maskpath = '' -obs_param_nml(49)%maskname = '' -obs_param_nml(49)%scalepath = '' -obs_param_nml(49)%scalename = '' -obs_param_nml(49)%flistpath = '' -obs_param_nml(49)%flistname = '' -obs_param_nml(49)%errstd = 9. -obs_param_nml(49)%std_normal_max = 2.5 -obs_param_nml(49)%zeromean = .true. -obs_param_nml(49)%coarsen_pert = .false. -obs_param_nml(49)%xcorr = 0.25 -obs_param_nml(49)%ycorr = 0.25 -obs_param_nml(49)%adapt = 0 +obs_param_nml(54)%descr = 'CYGNSS_SM' +obs_param_nml(54)%orbit = 3 +obs_param_nml(54)%pol = 0 +obs_param_nml(54)%N_ang = 0 +obs_param_nml(54)%freq = 0 +obs_param_nml(54)%FOV = 20. +obs_param_nml(54)%FOV_units = 'km' +obs_param_nml(54)%assim = .false. +obs_param_nml(54)%scale = .false. +obs_param_nml(54)%getinnov = .false. +obs_param_nml(54)%RTM_ID = 0 +obs_param_nml(54)%bias_Npar = 0 +obs_param_nml(54)%bias_trel = 864000 +obs_param_nml(54)%bias_tcut = 432000 +obs_param_nml(54)%nodata = -9999. +obs_param_nml(54)%varname = 'sfmc' +obs_param_nml(54)%units = 'cm3 cm-3' +obs_param_nml(54)%path = '/discover/nobackup/amfox/CYGNSS_obs/' +obs_param_nml(54)%name = '' +obs_param_nml(54)%maskpath = '' +obs_param_nml(54)%maskname = '' +obs_param_nml(54)%scalepath = '' +obs_param_nml(54)%scalename = '' +obs_param_nml(54)%flistpath = '' +obs_param_nml(54)%flistname = '' +obs_param_nml(54)%errstd = 9. +obs_param_nml(54)%std_normal_max = 2.5 +obs_param_nml(54)%zeromean = .true. +obs_param_nml(54)%coarsen_pert = .false. +obs_param_nml(54)%xcorr = 0.25 +obs_param_nml(54)%ycorr = 0.25 +obs_param_nml(54)%adapt = 0 ! -------------------------------------------------------------------- From e6b1dc8d01d3316d9217af072ef1fba1ca20a10b Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 10 Jan 2025 12:41:21 -0700 Subject: [PATCH 05/21] open and read netcdf file --- .../clsm_ensupd_read_obs.F90 | 200 +++++++++++++++++- 1 file changed, 191 insertions(+), 9 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index d79e279..f664ce3 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2190,26 +2190,53 @@ subroutine read_obs_sm_CYGNSS( & real, intent(out), dimension(N_catd) :: CYGNSS_sm_std ! sm obs err std [cm3 cm-3] real, intent(out), dimension(N_catd) :: CYGNSS_lon, CYGNSS_lat real*8, intent(out), dimension(N_catd) :: CYGNSS_time ! J2000 seconds + + ! ------------------------------------------------------------------- + ! local variables: - ! --------------- - - character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + integer, parameter :: max_obs = 100000 ! max number of obs read by subroutine (expecting < 6 hr assim window) + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' type(date_time_type) :: date_time_obs_beg, date_time_obs_end type(date_time_type) :: date_time_up, date_time_low - character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' - character(len=400) :: err_msg - + character(400) :: err_msg + character(300) :: tmpfname, tmpname + character( 2) :: MM, DD, HH, MI + character( 4) :: YYYY + + integer :: pos, i, idx, N_obs, j, HHcnt + integer :: ierr, ncid + integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid + integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid + integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop + + logical :: file_exists + logical :: HH03, HH09, HH12, HH15, HH21 + + real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_time(:) + + real, allocatable :: sm_d(:,:,:), sm_subd(:,:,:), sigma_d(:,:,:), sigma_subd(:,:,:) + real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:) + real :: tmp1_lon(max_obs), tmp1_lat(max_obs), tmp1_obs(max_obs), tmp1_err(max_obs), tmp1_time(max_obs) + ! ------------------------------------------------------------------- ! initialize + nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_time ) + found_obs = .false. + HH03 = .false. + HH09 = .false. + HH12 = .false. + HH15 = .false. + HH21 = .false. ! determine operating time range of sensor - if (trim(this_obs_param%descr) == 'CYGNSS_SM') then + if (trim(this_obs_param%descr) == 'CYGNSS_SM' .or. trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then date_time_obs_beg = date_time_type(2018, 8, 1, 0, 0, 0,-9999,-9999) date_time_obs_end = date_time_type(2100, 1, 1, 0, 0, 0,-9999,-9999) else @@ -2222,12 +2249,167 @@ subroutine read_obs_sm_CYGNSS( & if ( datetime_lt_refdatetime(date_time, date_time_obs_beg) .or. & datetime_lt_refdatetime(date_time_obs_end, date_time) ) return + ! find the time bounds of the assimilation window + ! (date_time-dtstep_assim/2,date_time+dtstep_assim/2] + + date_time_low = date_time + call augment_date_time( -(dtstep_assim/2), date_time_low) + date_time_up = date_time + call augment_date_time( (dtstep_assim/2), date_time_up) + + ! establish if assimilation window encompasses 12:00 UTC (daily CYGNSS obs time) + + if (date_time_low%hour < 12 .and. date_time_up%hour >= 12) then + HH12 = .true. + end if + + ! return if assimilating daily obs and HH12 is false + + if (trim(this_obs_param%descr) == 'CYGNSS_SM_daily' .and. .not. HH12) return + ! cygnss sm subdaily observations are for 4 timeslices per day - 00-06, 06-12, 12-18, 18-24 UTC ! we will assimilate these at 03, 09, 15, 21 UTC + + HHcnt = 0 + + if (date_time_low%hour < 3 .and. date_time_up%hour >= 3) then + HH03 = .true. + HHcnt = HHcnt + 1 + end if + if (date_time_low%hour > 21 .and. date_time_up%hour >= 3) then ! wrap-around from previous day + HH03 = .true. + HHcnt = HHcnt + 1 + end if + if (date_time_low%hour < 9 .and. date_time_up%hour >= 9) then + HH09 = .true. + HHcnt = HHcnt + 1 + end if + if (date_time_low%hour < 15 .and. date_time_up%hour >= 15) then + HH15 = .true. + HHcnt = HHcnt + 1 + end if + if (date_time_low%hour < 21 .and. date_time_up%hour >= 21) then + HH21 = .true. + HHcnt = HHcnt + 1 + end if + + ! return if assimilating subdaily obs and no HH is true + + if (trim(this_obs_param%descr) == 'CYGNSS_SM' .and. .not. (HH03 .or. HH09 .or. HH15 .or. HH21)) return + + ! if assimilation window encompasses more than one time sllice, abort + + if (HHcnt > 1) then + err_msg = 'Can not assimilate subdaily CYGNSS obs from more than one time slice' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! determine filename for CYGNSS obs file + + write (YYYY,'(i4.4)') date_time%year + write (MM, '(i2.2)') date_time%month + write (DD, '(i2.2)') date_time%day + write (HH, '(i2.2)') date_time%hour + write (MI, '(i2.2)') date_time%min + + tmpname = trim(this_obs_param%name) + + ! Replace s20180801 with sYYYYMMDD + pos = index(tmpname, "s20180801") + if (pos > 0) then + tmpname = tmpname(1:pos-1) // "s" // YYYY // MM // DD // tmpname(pos+9:) + endif + + ! Replace e20180801 with eYYYYMMDD + pos = index(tmpname, "e20180801") + if (pos > 0) then + tmpname = tmpname(1:pos-1) // "e" // YYYY // MM // DD // tmpname(pos+9:) + endif + + tmpfname = trim(this_obs_param%path) // '/Y' // YYYY // '/M' // MM // '/' // trim(tmpname) // '.nc' + + if (logit) write (logunit, '(400A)') 'Reading CYGNSS soil moisture data from file: ', trim(tmpfname) + + ! Check if file exists + + inquire(file=tmpfname, exist=file_exists) + + if (.not. file_exists) then + err_msg = 'CYGNSS SM obs file not found!' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! Open the NetCDF observation file + ierr = nf90_open(trim(tmpfname), nf90_nowrite, ncid) + + ! get variable dimension IDs + ierr = nf90_inq_dimid(ncid, 'lon', lon_dimid) + ierr = nf90_inq_dimid(ncid, 'lat', lat_dimid) + ierr = nf90_inq_dimid(ncid, 'time', time_dimid) + ierr = nf90_inq_dimid(ncid, 'timeslices', timeslices_dimid) + ierr = nf90_inq_dimid(ncid, 'startstop', startstop_dimid) + + ! dimensions sizes + ierr = nf90_inquire_dimension(ncid, lon_dimid, len=N_lon) + ierr = nf90_inquire_dimension(ncid, lat_dimid, len=N_lat) + ierr = nf90_inquire_dimension(ncid, time_dimid, len=N_time) + ierr = nf90_inquire_dimension(ncid, timeslices_dimid, len=N_timeslices) + ierr = nf90_inquire_dimension(ncid, startstop_dimid, len=N_startstop) + + ! get variable IDs + ierr = nf90_inq_varid(ncid, 'SM_daily', sm_d_varid) + ierr = nf90_inq_varid(ncid, 'SM_subdaily', sm_subd_varid) + ierr = nf90_inq_varid(ncid, 'SIGMA_daily', sigma_d_varid) + ierr = nf90_inq_varid(ncid, 'SIGMA_subdaily', sigma_subd_varid) + ierr = nf90_inq_varid(ncid, 'timeintervals', timeintervals_varid) + ierr = nf90_inq_varid(ncid, 'latitude', lat_varid) + ierr = nf90_inq_varid(ncid, 'longitude', lon_varid) + + ! allocate memory for the variables + allocate(sm_d(N_lon, N_lat, N_time)) + allocate(sm_subd(N_lon, N_lat, N_timeslices)) + allocate(sigma_d(N_lon, N_lat, N_time)) + allocate(sigma_subd(N_lon, N_lat, N_timeslices)) + allocate(timeintervals(N_startstop, N_timeslices)) + allocate(latitudes(N_lon, N_lat)) + allocate(longitudes(N_lon, N_lat)) + + ! read the variables + ierr = nf90_get_var(ncid, sm_d_varid, sm_d) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading sm_d:', ierr + end if + ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading sm_subd:', ierr + end if + ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading sigma_d:', ierr + end if + ierr = nf90_get_var(ncid, sigma_subd_varid, sigma_subd) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading sigma_subd:', ierr + end if + ierr = nf90_get_var(ncid, timeintervals_varid, timeintervals) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading timeintervals:', ierr + end if + ierr = nf90_get_var(ncid, lat_varid, latitudes) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading latitudes:', ierr + end if + ierr = nf90_get_var(ncid, lon_varid, longitudes) + if (ierr /= nf90_noerr) then + write(*,*) 'Error reading longitudes:', ierr + end if + + ! close the file + ierr = nf90_close(ncid) - write (logunit,*) 'This is where we will read CYGNSS obs at YY:', date_time%year, ' DD:', date_time%day, ' HH:', date_time%hour + write (logunit,*) 'This is where we will read CYGNSS obs at YY:', date_time%year, ' DD:', date_time%day, ' HH:', date_time%hour - return + return end subroutine read_obs_sm_CYGNSS From f3b15a0a7cb9369d1ef5dcc65ef4b3560bd45d65 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 10 Jan 2025 18:56:20 -0700 Subject: [PATCH 06/21] read mask and fill temp arrays --- .../clsm_ensupd_read_obs.F90 | 182 +++++++++++++++--- 1 file changed, 150 insertions(+), 32 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index f664ce3..2721b02 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2194,31 +2194,36 @@ subroutine read_obs_sm_CYGNSS( & ! ------------------------------------------------------------------- ! local variables: - integer, parameter :: max_obs = 100000 ! max number of obs read by subroutine (expecting < 6 hr assim window) - character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + integer, parameter :: max_obs = 50000 ! max number of obs read by subroutine (expecting < 6 hr assim window) + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' type(date_time_type) :: date_time_obs_beg, date_time_obs_end type(date_time_type) :: date_time_up, date_time_low character(400) :: err_msg - character(300) :: tmpfname, tmpname + character(300) :: tmpfname, tmpname, tmpmaskname character( 2) :: MM, DD, HH, MI character( 4) :: YYYY integer :: pos, i, idx, N_obs, j, HHcnt integer :: ierr, ncid - integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid + integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid, lon_dimid_m, lat_dimid_m integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid - integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop + integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop, N_lon_m, N_lat_m + integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid, few_obs_varid, low_signal_varid + integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) logical :: file_exists logical :: HH03, HH09, HH12, HH15, HH21 - real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_time(:) + real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_time(:) real, allocatable :: sm_d(:,:,:), sm_subd(:,:,:), sigma_d(:,:,:), sigma_subd(:,:,:) - real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:) + real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:), tmp_sm(:,:), tmp_sigma(:,:) + real, allocatable :: time(:), timeslices(:,:) + real, allocatable :: latitudes_m(:,:), longitudes_m(:,:) + real :: tmp1_lon(max_obs), tmp1_lat(max_obs), tmp1_obs(max_obs), tmp1_err(max_obs), tmp1_time(max_obs) ! ------------------------------------------------------------------- @@ -2373,39 +2378,152 @@ subroutine read_obs_sm_CYGNSS( & allocate(timeintervals(N_startstop, N_timeslices)) allocate(latitudes(N_lon, N_lat)) allocate(longitudes(N_lon, N_lat)) + allocate(tmp_sm(N_lon, N_lat)) + allocate(tmp_sigma(N_lon, N_lat)) ! read the variables - ierr = nf90_get_var(ncid, sm_d_varid, sm_d) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading sm_d:', ierr - end if - ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading sm_subd:', ierr - end if - ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading sigma_d:', ierr - end if - ierr = nf90_get_var(ncid, sigma_subd_varid, sigma_subd) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading sigma_subd:', ierr - end if ierr = nf90_get_var(ncid, timeintervals_varid, timeintervals) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading timeintervals:', ierr - end if ierr = nf90_get_var(ncid, lat_varid, latitudes) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading latitudes:', ierr - end if ierr = nf90_get_var(ncid, lon_varid, longitudes) - if (ierr /= nf90_noerr) then - write(*,*) 'Error reading longitudes:', ierr + + ! read either subdaily or daily soil moisture and sigma variables + if (this_obs_param%descr == 'CYGNSS_SM' ) then + + ! subdaily observations required + ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) + ierr = nf90_get_var(ncid, sigma_subd_varid, sigma_subd) + + ! find the time slice that corresponds to the assimilation window + idx = -1 + do i = 1, N_timeslices + if ((date_time%hour == 3 .and. timeintervals(1,i) == 0.0) .or. & + (date_time%hour == 9 .and. timeintervals(1,i) == 0.25) .or. & + (date_time%hour == 15 .and. timeintervals(1,i) == 0.5) .or. & + (date_time%hour == 21 .and. timeintervals(1,i) == 0.75)) then + idx = i + exit + end if + end do + if (idx == -1) then + err_msg = 'Error: No matching time interval found for HH = ' // HH + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + else + write(*,*) 'Index of matching time interval: ', idx + end if + tmp_sm = sm_subd(:,:,idx) + tmp_sigma = sigma_subd(:,:,idx) + else + + ! daily observations required + ierr = nf90_get_var(ncid, sm_d_varid, sm_d) + ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) + tmp_sm = sm_d(:,:,1) + tmp_sigma = sigma_d(:,:,1) end if - ! close the file + ! close the obs file ierr = nf90_close(ncid) + + ! get name for CYGNSS mask file + + tmpmaskname = trim(this_obs_param%maskpath) // '/' // trim(this_obs_param%maskname) // '.nc' + + inquire(file=tmpfname, exist=file_exists) + + if (.not. file_exists) then + err_msg = 'CYGNSS mask file not found!' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! open the CYGNSS mask file + + ierr = nf90_open(trim(tmpmaskname), nf90_nowrite, ncid) + + ! get variable dimension IDs + ierr = nf90_inq_dimid(ncid, 'lon', lon_dimid) + ierr = nf90_inq_dimid(ncid, 'lat', lat_dimid) + + ! dimensions sizes + ierr = nf90_inquire_dimension(ncid, lon_dimid, len=N_lon_m) + ierr = nf90_inquire_dimension(ncid, lat_dimid, len=N_lat_m) + + ! get variable IDs + ierr = nf90_inq_varid(ncid, 'longitude', longitudes_m_varid) + ierr = nf90_inq_varid(ncid, 'latitude', latitudes_m_varid) + ierr = nf90_inq_varid(ncid, 'flag_small_SM_range', small_SM_range_varid) + ierr = nf90_inq_varid(ncid, 'flag_poor_SMAP', poor_SMAP_varid) + ierr = nf90_inq_varid(ncid, 'flag_high_ubrmsd', high_ubrmsd_varid) + ierr = nf90_inq_varid(ncid, 'flag_few_obs', few_obs_varid) + ierr = nf90_inq_varid(ncid, 'flag_low_signal', low_signal_varid) + + + ! allocate memory for the variables + allocate(latitudes_m(N_lon_m, N_lat_m)) + allocate(longitudes_m(N_lon_m, N_lat_m)) + allocate(small_SM_range_flag(N_lon_m, N_lat_m)) + allocate(poor_SMAP_flag(N_lon_m, N_lat_m)) + allocate(high_ubrmsd_flag(N_lon_m, N_lat_m)) + allocate(few_obs_flag(N_lon_m, N_lat_m)) + allocate(low_signal_flag(N_lon_m, N_lat_m)) + + ! read the variables + ierr = nf90_get_var(ncid, latitudes_m_varid, latitudes_m) + ierr = nf90_get_var(ncid, longitudes_m_varid, longitudes_m) + ierr = nf90_get_var(ncid, small_SM_range_varid, small_SM_range_flag) + ierr = nf90_get_var(ncid, poor_SMAP_varid, poor_SMAP_flag) + ierr = nf90_get_var(ncid, high_ubrmsd_varid, high_ubrmsd_flag) + ierr = nf90_get_var(ncid, few_obs_varid, few_obs_flag) + ierr = nf90_get_var(ncid, low_signal_varid, low_signal_flag) + + ! close the mask file + ierr = nf90_close(ncid) + + ! check the obs data and mask data are the same resolution + if (N_lon /= N_lon_m .or. N_lat /= N_lat_m) then + err_msg = 'The mask file ' // trim(this_obs_param%maskname) // ' does not match the obs resolution' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + + ! fill tmp arrays + N_obs = 0 + + do i = 1, N_lon + do j = 1, N_lat + if (tmp_sm(i,j) .ne. this_obs_param%nodata .and. & + small_SM_range_flag(i,j) .ne. 1 .and. & + poor_SMAP_flag(i,j) .ne. 1 .and. & + high_ubrmsd_flag(i,j) .ne. 1 .and. & + few_obs_flag(i,j) .ne. 1 .and. & + low_signal_flag(i,j) .ne. 1 ) then + + ! valid observation + N_obs = N_obs + 1 + if (N_obs > max_obs) then + err_msg = 'Attempting to read too many obs - how long is your assimilation window?' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + tmp1_lon(N_obs) = longitudes(i,j) + tmp1_lat(N_obs) = latitudes(i,j) + tmp1_obs(N_obs) = tmp_sm(i,j) + tmp1_err(N_obs) = tmp_sigma(i,j) + tmp1_time(N_obs) = date_time%year + end if + end do + end do + + write(*,*) 'Number of observations: ', N_obs + + allocate(tmp_lon(N_obs)) + allocate(tmp_lat(N_obs)) + allocate(tmp_obs(N_obs)) + allocate(tmp_err(N_obs)) + allocate(tmp_time(N_obs)) + + tmp_lon = tmp1_lon(1:N_obs) + tmp_lat = tmp1_lat(1:N_obs) + tmp_obs = tmp1_obs(1:N_obs) + tmp_err = tmp1_err(1:N_obs) + tmp_time = tmp1_time(1:N_obs) write (logunit,*) 'This is where we will read CYGNSS obs at YY:', date_time%year, ' DD:', date_time%day, ' HH:', date_time%hour From 3f9a03ceb36dd55bbdf041aac3de9b7b6084c83c Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 13 Jan 2025 11:51:31 -0700 Subject: [PATCH 07/21] add CYGNSS_SM_daily --- .../clsm_ensupd_glob_param.F90 | 2 +- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 40 +++++++++++++++++++ 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 b/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 index 6d1ef60..aebfb3e 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_glob_param.F90 @@ -48,7 +48,7 @@ module clsm_ensupd_glob_param ! total number of all obs species defined in "ensupd" namelist file ! (regardless of whether "assim" flag is true or false) - integer, parameter :: N_obs_species_nml = 54 + integer, parameter :: N_obs_species_nml = 55 ! ---------------------------------------------------------------------- ! diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index 3724d02..c6ac56b 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2402,8 +2402,48 @@ obs_param_nml(54)%xcorr = 0.25 obs_param_nml(54)%ycorr = 0.25 obs_param_nml(54)%adapt = 0 +! -------------------------------------------------------------------------------------------------------------- +! +! 55 = CYGNSS_SM_daily (CYGNSS Level 3 soil moisture version 3.2) +! +! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 + +obs_param_nml(55)%descr = 'CYGNSS_SM_daily' +obs_param_nml(55)%orbit = 3 +obs_param_nml(55)%pol = 0 +obs_param_nml(55)%N_ang = 0 +obs_param_nml(55)%freq = 0 +obs_param_nml(55)%FOV = 20. +obs_param_nml(55)%FOV_units = 'km' +obs_param_nml(55)%assim = .false. +obs_param_nml(55)%scale = .false. +obs_param_nml(55)%getinnov = .false. +obs_param_nml(55)%RTM_ID = 0 +obs_param_nml(55)%bias_Npar = 0 +obs_param_nml(55)%bias_trel = 864000 +obs_param_nml(55)%bias_tcut = 432000 +obs_param_nml(55)%nodata = -9999. +obs_param_nml(55)%varname = 'sfmc' +obs_param_nml(55)%units = 'cm3 cm-3' +obs_param_nml(55)%path = '/discover/nobackup/amfox/CYGNSS_obs/' +obs_param_nml(55)%name = '' +obs_param_nml(55)%maskpath = '' +obs_param_nml(55)%maskname = '' +obs_param_nml(55)%scalepath = '' +obs_param_nml(55)%scalename = '' +obs_param_nml(55)%flistpath = '' +obs_param_nml(55)%flistname = '' +obs_param_nml(55)%errstd = 9. +obs_param_nml(55)%std_normal_max = 2.5 +obs_param_nml(55)%zeromean = .true. +obs_param_nml(55)%coarsen_pert = .false. +obs_param_nml(55)%xcorr = 0.25 +obs_param_nml(55)%ycorr = 0.25 +obs_param_nml(55)%adapt = 0 + ! -------------------------------------------------------------------- + / ! =========================== EOF ======================================= From caf2d253d6560dec2a4baed6516374ae94cf5435 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 13 Jan 2025 16:54:17 -0700 Subject: [PATCH 08/21] tile space --- .../clsm_ensupd_read_obs.F90 | 146 +++++++++++++++--- 1 file changed, 126 insertions(+), 20 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 2721b02..c5dbb1b 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2206,31 +2206,33 @@ subroutine read_obs_sm_CYGNSS( & character( 2) :: MM, DD, HH, MI character( 4) :: YYYY - integer :: pos, i, idx, N_obs, j, HHcnt + integer :: pos, i, idx, N_obs, j, HHcnt, ind, ii integer :: ierr, ncid integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid, lon_dimid_m, lat_dimid_m integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop, N_lon_m, N_lat_m integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid, few_obs_varid, low_signal_varid - integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) + integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) logical :: file_exists logical :: HH03, HH09, HH12, HH15, HH21 - real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_time(:) - + real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_jtime(:) + integer, dimension(:), pointer :: tmp_tile_num + integer, dimension(N_catd) :: N_obs_in_tile + real, allocatable :: sm_d(:,:,:), sm_subd(:,:,:), sigma_d(:,:,:), sigma_subd(:,:,:) real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:), tmp_sm(:,:), tmp_sigma(:,:) real, allocatable :: time(:), timeslices(:,:) real, allocatable :: latitudes_m(:,:), longitudes_m(:,:) - real :: tmp1_lon(max_obs), tmp1_lat(max_obs), tmp1_obs(max_obs), tmp1_err(max_obs), tmp1_time(max_obs) + real :: tmp1_lon(max_obs), tmp1_lat(max_obs), tmp1_obs(max_obs), tmp1_err(max_obs), tmp1_jtime(max_obs) ! ------------------------------------------------------------------- ! initialize - nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_time ) + nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_jtime ) found_obs = .false. HH03 = .false. @@ -2502,32 +2504,136 @@ subroutine read_obs_sm_CYGNSS( & err_msg = 'Attempting to read too many obs - how long is your assimilation window?' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - tmp1_lon(N_obs) = longitudes(i,j) - tmp1_lat(N_obs) = latitudes(i,j) - tmp1_obs(N_obs) = tmp_sm(i,j) - tmp1_err(N_obs) = tmp_sigma(i,j) - tmp1_time(N_obs) = date_time%year + tmp1_lon( N_obs) = longitudes(i,j) + tmp1_lat( N_obs) = latitudes(i,j) + tmp1_obs( N_obs) = tmp_sm(i,j) + tmp1_err( N_obs) = tmp_sigma(i,j) + tmp1_jtime(N_obs) = datetime_to_J2000seconds( date_time, J2000_epoch_id ) end if end do end do - write(*,*) 'Number of observations: ', N_obs + if (logit) then + + write (logunit,*) 'read_obs_sm_CYGNSS: read ', N_obs, ' at date_time = ', date_time, ' from: ', tmpfname + write (logunit,*) '----------' + write (logunit,*) 'max(obs)=',maxval(tmp1_obs(1:N_obs)), ', min(obs)=',minval(tmp1_obs(1:N_obs)), & + ', avg(obs)=',sum(tmp1_obs(1:N_obs))/N_obs + + end if allocate(tmp_lon(N_obs)) allocate(tmp_lat(N_obs)) allocate(tmp_obs(N_obs)) allocate(tmp_err(N_obs)) - allocate(tmp_time(N_obs)) + allocate(tmp_jtime(N_obs)) + + tmp_lon = tmp1_lon(1:N_obs) + tmp_lat = tmp1_lat(1:N_obs) + tmp_obs = tmp1_obs(1:N_obs) + tmp_err = tmp1_err(1:N_obs) + tmp_jtime = tmp1_jtime(1:N_obs) + + ! ---------------------------------------------------------------- + ! + ! 2.) for each observation + ! a) determine grid cell that contains lat/lon + ! b) determine tile within grid cell that contains lat/lon + + if (N_obs>0) then + + allocate(tmp_tile_num(N_obs)) + + call get_tile_num_for_obs(N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + N_obs, tmp_lat, tmp_lon, & + this_obs_param, & + tmp_tile_num ) + + ! ---------------------------------------------------------------- + ! + ! 3.) compute super-obs for each tile from all obs w/in that tile + ! (also eliminate observations that are not in domain) + + CYGNSS_sm = 0. + CYGNSS_lon = 0. + CYGNSS_lat = 0. + CYGNSS_time = 0.0D0 + + N_obs_in_tile = 0 + + do ii=1,N_obs + + ind = tmp_tile_num(ii) ! 1<=tmp_tile_num<=N_catd (unless nodata) + + if (ind>0) then ! this step eliminates obs outside domain + + CYGNSS_sm( ind) = CYGNSS_sm( ind) + tmp_obs( ii) + CYGNSS_lon( ind) = CYGNSS_lon( ind) + tmp_lon( ii) + CYGNSS_lat( ind) = CYGNSS_lat( ind) + tmp_lat( ii) + CYGNSS_time(ind) = CYGNSS_time(ind) + tmp_jtime(ii) + + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 + + end if + + end do - tmp_lon = tmp1_lon(1:N_obs) - tmp_lat = tmp1_lat(1:N_obs) - tmp_obs = tmp1_obs(1:N_obs) - tmp_err = tmp1_err(1:N_obs) - tmp_time = tmp1_time(1:N_obs) + ! normalize and set obs error std-dev + + do ii=1,N_catd - write (logunit,*) 'This is where we will read CYGNSS obs at YY:', date_time%year, ' DD:', date_time%day, ' HH:', date_time%hour + ! set observation error standard deviation + + CYGNSS_sm_std(ii) = this_obs_param%errstd/100. ! change units from percent (0-100) to fraction (0-1) + + ! normalize - return + if (N_obs_in_tile(ii)>1) then + + CYGNSS_sm( ii) = CYGNSS_sm( ii)/real(N_obs_in_tile(ii)) + CYGNSS_lon( ii) = CYGNSS_lon( ii)/real(N_obs_in_tile(ii)) + CYGNSS_lat( ii) = CYGNSS_lat( ii)/real(N_obs_in_tile(ii)) + CYGNSS_time( ii) = CYGNSS_time(ii)/real(N_obs_in_tile(ii),kind(0.0D0)) + + elseif (N_obs_in_tile(ii)==0) then + + CYGNSS_sm( ii) = this_obs_param%nodata + CYGNSS_lon( ii) = this_obs_param%nodata + CYGNSS_lat( ii) = this_obs_param%nodata + CYGNSS_time( ii) = real(this_obs_param%nodata,kind(0.0D0)) + CYGNSS_sm_std(ii) = this_obs_param%nodata + + else + + ! nothing to do if N_obs_in_tile(ii)==1 (and assuming N_obs_in_tile is never negative) + + end if + + end do + + ! clean up + + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) + + if (any(N_obs_in_tile>0)) then + + found_obs = .true. + + else + + found_obs = .false. + + end if + + end if + + ! clean up + + if (associated(tmp_obs)) deallocate(tmp_obs) + if (associated(tmp_lon)) deallocate(tmp_lon) + if (associated(tmp_lat)) deallocate(tmp_lat) + if (associated(tmp_jtime)) deallocate(tmp_jtime) end subroutine read_obs_sm_CYGNSS From 2085488be49d4348630b18a597287d3da39a230d Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 17 Jan 2025 15:09:25 -0700 Subject: [PATCH 09/21] comments and error units --- .../clsm_ensupd_read_obs.F90 | 44 ++++++++++++------- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 4 +- 2 files changed, 29 insertions(+), 19 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index c5dbb1b..b0050a4 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -1854,19 +1854,19 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! open bufr file -! call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it -! open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') -! call openbf(lnbufr, 'SEC3', lnbufr) -! call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) -! call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) + call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it + open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') + call openbf(lnbufr, 'SEC3', lnbufr) + call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) + call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) -! msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) + msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) -! loop_report: do while( ireadsb(lnbufr) == 0 ) + loop_report: do while( ireadsb(lnbufr) == 0 ) ! columns of tmp_data: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 -! call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') + call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') N_obs = N_obs + 1 @@ -1877,11 +1877,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & tmp_data(N_obs,:) = tmp_vdata -! end do loop_report -! end do msg_report + end do loop_report + end do msg_report -! call closbf(lnbufr) -! close(lnbufr) + call closbf(lnbufr) + close(lnbufr) end do ! end file loop @@ -2148,7 +2148,7 @@ end subroutine read_obs_sm_ASCAT_EUMET ! **************************************************************************** - subroutine read_obs_sm_CYGNSS( & + subroutine read_obs_sm_CYGNSS( & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & @@ -2156,7 +2156,17 @@ subroutine read_obs_sm_CYGNSS( & !--------------------------------------------------------------------- ! - ! Routine to read in CYGNSS soil moisture obs from netCDF files. + ! Routine to read in CYGNSS soil moisture obs from netCDF files and apply mask. + ! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 + ! CYGNSS soil moisture obs are in volumetric units (cm3 cm-3) with a spatial + ! resolution of 0.3 Decimal Degrees x 0.37 Decimal Degrees. + ! They are daily files countaining daily (obs_param_nml(54)%descr = CYGNSS_SM_daily) + ! and subdaily SM estimates (obs_param_nml(53)%descr = CYGNSS_SM). The + ! subdaily values are for 6 hr time intervals, 0-6hr, 6-12hr, 12-18hr, 18-24hr. + ! The subdaily values are assimilated at 3z, 9z, 15z, 21z. If the assimilation + ! window is more than 6 hr the code will error out. The daily values are assimilated + ! at 12z. The assumption is that either the daily or subdaily values are assimilated, + ! but the code will assimilate both if in the nml. ! ! A. Fox, Jan 2025 ! @@ -2536,7 +2546,7 @@ subroutine read_obs_sm_CYGNSS( & ! ---------------------------------------------------------------- ! - ! 2.) for each observation + ! for each observation ! a) determine grid cell that contains lat/lon ! b) determine tile within grid cell that contains lat/lon @@ -2552,7 +2562,7 @@ subroutine read_obs_sm_CYGNSS( & ! ---------------------------------------------------------------- ! - ! 3.) compute super-obs for each tile from all obs w/in that tile + ! compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) CYGNSS_sm = 0. @@ -2585,7 +2595,7 @@ subroutine read_obs_sm_CYGNSS( & ! set observation error standard deviation - CYGNSS_sm_std(ii) = this_obs_param%errstd/100. ! change units from percent (0-100) to fraction (0-1) + CYGNSS_sm_std(ii) = this_obs_param%errstd ! fraction units (cm3 cm-3) ! normalize diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index c6ac56b..d42301d 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2394,7 +2394,7 @@ obs_param_nml(54)%scalepath = '' obs_param_nml(54)%scalename = '' obs_param_nml(54)%flistpath = '' obs_param_nml(54)%flistname = '' -obs_param_nml(54)%errstd = 9. +obs_param_nml(54)%errstd = 0.1 obs_param_nml(54)%std_normal_max = 2.5 obs_param_nml(54)%zeromean = .true. obs_param_nml(54)%coarsen_pert = .false. @@ -2433,7 +2433,7 @@ obs_param_nml(55)%scalepath = '' obs_param_nml(55)%scalename = '' obs_param_nml(55)%flistpath = '' obs_param_nml(55)%flistname = '' -obs_param_nml(55)%errstd = 9. +obs_param_nml(55)%errstd = 0.1 obs_param_nml(55)%std_normal_max = 2.5 obs_param_nml(55)%zeromean = .true. obs_param_nml(55)%coarsen_pert = .false. From cb49cb0b81ba65591d5009432fd25a84f44990a1 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 17 Jan 2025 15:16:48 -0700 Subject: [PATCH 10/21] todo: scaling --- GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index b0050a4..9421b7c 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -9149,9 +9149,7 @@ subroutine read_obs( & scaled_obs = .true. - call scale_obs_sfmc_zscore( N_catd, tile_coord, & - date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & - tmp_obs, tmp_std_obs ) + ! TODO: implement scaling for CYGNSS obs end if From 9ca0f0f4eff79f3731148ce8660a20dea4430cb8 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Mon, 20 Jan 2025 13:20:52 -0500 Subject: [PATCH 11/21] Fixed typos and alignment in new CYGNSS soil moisture obs reader (clsm_ensupd_read_obs.F90) --- .../clsm_ensupd_read_obs.F90 | 85 ++++++++++--------- 1 file changed, 44 insertions(+), 41 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 9421b7c..97fbd59 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2158,10 +2158,10 @@ subroutine read_obs_sm_CYGNSS( & ! ! Routine to read in CYGNSS soil moisture obs from netCDF files and apply mask. ! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 - ! CYGNSS soil moisture obs are in volumetric units (cm3 cm-3) with a spatial + ! CYGNSS soil moisture obs are in volumetric units (m3 m-3) with a spatial ! resolution of 0.3 Decimal Degrees x 0.37 Decimal Degrees. - ! They are daily files countaining daily (obs_param_nml(54)%descr = CYGNSS_SM_daily) - ! and subdaily SM estimates (obs_param_nml(53)%descr = CYGNSS_SM). The + ! They are daily files containing daily (obs_param_nml(.)%descr = CYGNSS_SM_daily) + ! and subdaily SM estimates (obs_param_nml(.)%descr = CYGNSS_SM). The ! subdaily values are for 6 hr time intervals, 0-6hr, 6-12hr, 12-18hr, 18-24hr. ! The subdaily values are assimilated at 3z, 9z, 15z, 21z. If the assimilation ! window is more than 6 hr the code will error out. The daily values are assimilated @@ -2245,6 +2245,7 @@ subroutine read_obs_sm_CYGNSS( & nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_jtime ) found_obs = .false. + HH03 = .false. HH09 = .false. HH12 = .false. @@ -2314,7 +2315,7 @@ subroutine read_obs_sm_CYGNSS( & if (trim(this_obs_param%descr) == 'CYGNSS_SM' .and. .not. (HH03 .or. HH09 .or. HH15 .or. HH21)) return - ! if assimilation window encompasses more than one time sllice, abort + ! if assimilation window encompasses more than one time slice, abort if (HHcnt > 1) then err_msg = 'Can not assimilate subdaily CYGNSS obs from more than one time slice' @@ -2383,15 +2384,15 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_inq_varid(ncid, 'longitude', lon_varid) ! allocate memory for the variables - allocate(sm_d(N_lon, N_lat, N_time)) - allocate(sm_subd(N_lon, N_lat, N_timeslices)) - allocate(sigma_d(N_lon, N_lat, N_time)) - allocate(sigma_subd(N_lon, N_lat, N_timeslices)) - allocate(timeintervals(N_startstop, N_timeslices)) - allocate(latitudes(N_lon, N_lat)) - allocate(longitudes(N_lon, N_lat)) - allocate(tmp_sm(N_lon, N_lat)) - allocate(tmp_sigma(N_lon, N_lat)) + allocate(sm_d( N_lon, N_lat, N_time )) + allocate(sm_subd( N_lon, N_lat, N_timeslices )) + allocate(sigma_d( N_lon, N_lat, N_time )) + allocate(sigma_subd( N_lon, N_lat, N_timeslices )) + allocate(timeintervals(N_startstop, N_timeslices )) + allocate(latitudes( N_lon, N_lat )) + allocate(longitudes( N_lon, N_lat )) + allocate(tmp_sm( N_lon, N_lat )) + allocate(tmp_sigma( N_lon, N_lat )) ! read the variables ierr = nf90_get_var(ncid, timeintervals_varid, timeintervals) @@ -2399,19 +2400,19 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_get_var(ncid, lon_varid, longitudes) ! read either subdaily or daily soil moisture and sigma variables - if (this_obs_param%descr == 'CYGNSS_SM' ) then + if ( trim(this_obs_param%descr) == 'CYGNSS_SM' ) then ! subdaily observations required - ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) + ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) ierr = nf90_get_var(ncid, sigma_subd_varid, sigma_subd) ! find the time slice that corresponds to the assimilation window idx = -1 do i = 1, N_timeslices - if ((date_time%hour == 3 .and. timeintervals(1,i) == 0.0) .or. & - (date_time%hour == 9 .and. timeintervals(1,i) == 0.25) .or. & - (date_time%hour == 15 .and. timeintervals(1,i) == 0.5) .or. & - (date_time%hour == 21 .and. timeintervals(1,i) == 0.75)) then + if ((date_time%hour == 3 .and. timeintervals(1,i) == 0.0 ) .or. & + (date_time%hour == 9 .and. timeintervals(1,i) == 0.25) .or. & + (date_time%hour == 15 .and. timeintervals(1,i) == 0.5 ) .or. & + (date_time%hour == 21 .and. timeintervals(1,i) == 0.75) ) then idx = i exit end if @@ -2422,15 +2423,17 @@ subroutine read_obs_sm_CYGNSS( & else write(*,*) 'Index of matching time interval: ', idx end if - tmp_sm = sm_subd(:,:,idx) + tmp_sm = sm_subd( :,:,idx) tmp_sigma = sigma_subd(:,:,idx) + else ! daily observations required - ierr = nf90_get_var(ncid, sm_d_varid, sm_d) + ierr = nf90_get_var(ncid, sm_d_varid, sm_d) ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) - tmp_sm = sm_d(:,:,1) + tmp_sm = sm_d( :,:,1) tmp_sigma = sigma_d(:,:,1) + end if ! close the obs file @@ -2470,13 +2473,13 @@ subroutine read_obs_sm_CYGNSS( & ! allocate memory for the variables - allocate(latitudes_m(N_lon_m, N_lat_m)) - allocate(longitudes_m(N_lon_m, N_lat_m)) + allocate(latitudes_m( N_lon_m, N_lat_m)) + allocate(longitudes_m( N_lon_m, N_lat_m)) allocate(small_SM_range_flag(N_lon_m, N_lat_m)) - allocate(poor_SMAP_flag(N_lon_m, N_lat_m)) - allocate(high_ubrmsd_flag(N_lon_m, N_lat_m)) - allocate(few_obs_flag(N_lon_m, N_lat_m)) - allocate(low_signal_flag(N_lon_m, N_lat_m)) + allocate(poor_SMAP_flag( N_lon_m, N_lat_m)) + allocate(high_ubrmsd_flag( N_lon_m, N_lat_m)) + allocate(few_obs_flag( N_lon_m, N_lat_m)) + allocate(low_signal_flag( N_lon_m, N_lat_m)) ! read the variables ierr = nf90_get_var(ncid, latitudes_m_varid, latitudes_m) @@ -2514,10 +2517,10 @@ subroutine read_obs_sm_CYGNSS( & err_msg = 'Attempting to read too many obs - how long is your assimilation window?' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - tmp1_lon( N_obs) = longitudes(i,j) - tmp1_lat( N_obs) = latitudes(i,j) - tmp1_obs( N_obs) = tmp_sm(i,j) - tmp1_err( N_obs) = tmp_sigma(i,j) + tmp1_lon( N_obs) = longitudes(i,j) + tmp1_lat( N_obs) = latitudes( i,j) + tmp1_obs( N_obs) = tmp_sm( i,j) + tmp1_err( N_obs) = tmp_sigma( i,j) tmp1_jtime(N_obs) = datetime_to_J2000seconds( date_time, J2000_epoch_id ) end if end do @@ -2532,16 +2535,16 @@ subroutine read_obs_sm_CYGNSS( & end if - allocate(tmp_lon(N_obs)) - allocate(tmp_lat(N_obs)) - allocate(tmp_obs(N_obs)) - allocate(tmp_err(N_obs)) + allocate(tmp_lon( N_obs)) + allocate(tmp_lat( N_obs)) + allocate(tmp_obs( N_obs)) + allocate(tmp_err( N_obs)) allocate(tmp_jtime(N_obs)) - tmp_lon = tmp1_lon(1:N_obs) - tmp_lat = tmp1_lat(1:N_obs) - tmp_obs = tmp1_obs(1:N_obs) - tmp_err = tmp1_err(1:N_obs) + tmp_lon = tmp1_lon( 1:N_obs) + tmp_lat = tmp1_lat( 1:N_obs) + tmp_obs = tmp1_obs( 1:N_obs) + tmp_err = tmp1_err( 1:N_obs) tmp_jtime = tmp1_jtime(1:N_obs) ! ---------------------------------------------------------------- @@ -2595,7 +2598,7 @@ subroutine read_obs_sm_CYGNSS( & ! set observation error standard deviation - CYGNSS_sm_std(ii) = this_obs_param%errstd ! fraction units (cm3 cm-3) + CYGNSS_sm_std(ii) = this_obs_param%errstd ! volumetric units (m3 m-3) ! normalize From 95775040705d67b04be11bcb9a004765a2021e55 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Mon, 20 Jan 2025 13:51:22 -0500 Subject: [PATCH 12/21] Changed units of CYGNSS soil moisture obs to SI standard (LDASsa_DEFAULT_inputs_ensupd.nml) --- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index d42301d..a4bfe27 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2385,7 +2385,7 @@ obs_param_nml(54)%bias_trel = 864000 obs_param_nml(54)%bias_tcut = 432000 obs_param_nml(54)%nodata = -9999. obs_param_nml(54)%varname = 'sfmc' -obs_param_nml(54)%units = 'cm3 cm-3' +obs_param_nml(54)%units = 'm3 m-3' obs_param_nml(54)%path = '/discover/nobackup/amfox/CYGNSS_obs/' obs_param_nml(54)%name = '' obs_param_nml(54)%maskpath = '' @@ -2424,7 +2424,7 @@ obs_param_nml(55)%bias_trel = 864000 obs_param_nml(55)%bias_tcut = 432000 obs_param_nml(55)%nodata = -9999. obs_param_nml(55)%varname = 'sfmc' -obs_param_nml(55)%units = 'cm3 cm-3' +obs_param_nml(55)%units = 'm3 m-3' obs_param_nml(55)%path = '/discover/nobackup/amfox/CYGNSS_obs/' obs_param_nml(55)%name = '' obs_param_nml(55)%maskpath = '' From 6fe11e9892ce24fe8f32eb96995bdd6f5b472c98 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 24 Jan 2025 16:38:51 -0700 Subject: [PATCH 13/21] first pass at addressing review comments --- CHANGELOG.md | 1 + .../clsm_ensupd_read_obs.F90 | 82 +++++++++++-------- 2 files changed, 50 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 390b789..c7bea22 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- CYGNSS soil moisture reader - New update_type for joint 3d soil moisture and 1d snow analysis (Tb+sfmc+sfds+SCF obs) ### Changed diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 97fbd59..31ed9a0 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2161,7 +2161,7 @@ subroutine read_obs_sm_CYGNSS( & ! CYGNSS soil moisture obs are in volumetric units (m3 m-3) with a spatial ! resolution of 0.3 Decimal Degrees x 0.37 Decimal Degrees. ! They are daily files containing daily (obs_param_nml(.)%descr = CYGNSS_SM_daily) - ! and subdaily SM estimates (obs_param_nml(.)%descr = CYGNSS_SM). The + ! and subdaily SM estimates (obs_param_nml(.)%descr = CYGNSS_SM_6hr). The ! subdaily values are for 6 hr time intervals, 0-6hr, 6-12hr, 12-18hr, 18-24hr. ! The subdaily values are assimilated at 3z, 9z, 15z, 21z. If the assimilation ! window is more than 6 hr the code will error out. The daily values are assimilated @@ -2204,8 +2204,9 @@ subroutine read_obs_sm_CYGNSS( & ! ------------------------------------------------------------------- ! local variables: - integer, parameter :: max_obs = 50000 ! max number of obs read by subroutine (expecting < 6 hr assim window) - character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 + + integer, parameter :: max_obs = 50000 ! max number of daily obs read by subroutine (expecting > 6 hr assim window) + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' type(date_time_type) :: date_time_obs_beg, date_time_obs_end @@ -2254,7 +2255,7 @@ subroutine read_obs_sm_CYGNSS( & ! determine operating time range of sensor - if (trim(this_obs_param%descr) == 'CYGNSS_SM' .or. trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then + if (trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' .or. trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then date_time_obs_beg = date_time_type(2018, 8, 1, 0, 0, 0,-9999,-9999) date_time_obs_end = date_time_type(2100, 1, 1, 0, 0, 0,-9999,-9999) else @@ -2294,7 +2295,7 @@ subroutine read_obs_sm_CYGNSS( & HH03 = .true. HHcnt = HHcnt + 1 end if - if (date_time_low%hour > 21 .and. date_time_up%hour >= 3) then ! wrap-around from previous day + if (date_time_low%hour > 21 .and. date_time_up%hour >= 3) then ! wrap-around from previous day for 3 hr < assim window <= 6 hr HH03 = .true. HHcnt = HHcnt + 1 end if @@ -2310,10 +2311,14 @@ subroutine read_obs_sm_CYGNSS( & HH21 = .true. HHcnt = HHcnt + 1 end if + if (date_time_low%hour < 21 .and. date_time_up%hour < 3) then ! wrap-around into next day for 3 hr < assim window <= 6 hr + HH21 = .true. + HHcnt = HHcnt + 1 + end if ! return if assimilating subdaily obs and no HH is true - if (trim(this_obs_param%descr) == 'CYGNSS_SM' .and. .not. (HH03 .or. HH09 .or. HH15 .or. HH21)) return + if (trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' .and. .not. (HH03 .or. HH09 .or. HH15 .or. HH21)) return ! if assimilation window encompasses more than one time slice, abort @@ -2330,25 +2335,15 @@ subroutine read_obs_sm_CYGNSS( & write (HH, '(i2.2)') date_time%hour write (MI, '(i2.2)') date_time%min - tmpname = trim(this_obs_param%name) - - ! Replace s20180801 with sYYYYMMDD - pos = index(tmpname, "s20180801") - if (pos > 0) then - tmpname = tmpname(1:pos-1) // "s" // YYYY // MM // DD // tmpname(pos+9:) - endif + ! construct the file name using hardwired name and date_time - ! Replace e20180801 with eYYYYMMDD - pos = index(tmpname, "e20180801") - if (pos > 0) then - tmpname = tmpname(1:pos-1) // "e" // YYYY // MM // DD // tmpname(pos+9:) - endif + tmpname = 'cyg.ddmi.s'// YYYY // MM // DD //'-030000-e'// YYYY // MM // DD //'-210000.l3.grid-soil-moisture-36km.a32.d33.nc' tmpfname = trim(this_obs_param%path) // '/Y' // YYYY // '/M' // MM // '/' // trim(tmpname) // '.nc' if (logit) write (logunit, '(400A)') 'Reading CYGNSS soil moisture data from file: ', trim(tmpfname) - ! Check if file exists + ! check if file exists inquire(file=tmpfname, exist=file_exists) @@ -2357,7 +2352,7 @@ subroutine read_obs_sm_CYGNSS( & call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - ! Open the NetCDF observation file + ! open the NetCDF observation file ierr = nf90_open(trim(tmpfname), nf90_nowrite, ncid) ! get variable dimension IDs @@ -2383,11 +2378,7 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_inq_varid(ncid, 'latitude', lat_varid) ierr = nf90_inq_varid(ncid, 'longitude', lon_varid) - ! allocate memory for the variables - allocate(sm_d( N_lon, N_lat, N_time )) - allocate(sm_subd( N_lon, N_lat, N_timeslices )) - allocate(sigma_d( N_lon, N_lat, N_time )) - allocate(sigma_subd( N_lon, N_lat, N_timeslices )) + ! allocate memory for universal variables allocate(timeintervals(N_startstop, N_timeslices )) allocate(latitudes( N_lon, N_lat )) allocate(longitudes( N_lon, N_lat )) @@ -2400,7 +2391,11 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_get_var(ncid, lon_varid, longitudes) ! read either subdaily or daily soil moisture and sigma variables - if ( trim(this_obs_param%descr) == 'CYGNSS_SM' ) then + if ( trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' ) then + + ! allocate memory for subdaily variables + allocate(sm_subd( N_lon, N_lat, N_timeslices )) + allocate(sigma_subd(N_lon, N_lat, N_timeslices )) ! subdaily observations required ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) @@ -2409,10 +2404,10 @@ subroutine read_obs_sm_CYGNSS( & ! find the time slice that corresponds to the assimilation window idx = -1 do i = 1, N_timeslices - if ((date_time%hour == 3 .and. timeintervals(1,i) == 0.0 ) .or. & - (date_time%hour == 9 .and. timeintervals(1,i) == 0.25) .or. & - (date_time%hour == 15 .and. timeintervals(1,i) == 0.5 ) .or. & - (date_time%hour == 21 .and. timeintervals(1,i) == 0.75) ) then + if ((HH03 .and. timeintervals(1,i) == 0.0 ) .or. & + (HH09 .and. timeintervals(1,i) == 0.25) .or. & + (HH15 .and. timeintervals(1,i) == 0.5 ) .or. & + (HH21 .and. timeintervals(1,i) == 0.75) ) then idx = i exit end if @@ -2428,6 +2423,10 @@ subroutine read_obs_sm_CYGNSS( & else + ! allocate memory for daily variables + allocate(sm_d( N_lon, N_lat, N_time )) + allocate(sigma_d(N_lon, N_lat, N_time )) + ! daily observations required ierr = nf90_get_var(ncid, sm_d_varid, sm_d) ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) @@ -2584,7 +2583,7 @@ subroutine read_obs_sm_CYGNSS( & CYGNSS_sm( ind) = CYGNSS_sm( ind) + tmp_obs( ii) CYGNSS_lon( ind) = CYGNSS_lon( ind) + tmp_lon( ii) CYGNSS_lat( ind) = CYGNSS_lat( ind) + tmp_lat( ii) - CYGNSS_time(ind) = CYGNSS_time(ind) + tmp_jtime(ii) + CYGNSS_time(ind) = tmp_jtime(ii) ! time is the same for all observations N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 @@ -2607,7 +2606,6 @@ subroutine read_obs_sm_CYGNSS( & CYGNSS_sm( ii) = CYGNSS_sm( ii)/real(N_obs_in_tile(ii)) CYGNSS_lon( ii) = CYGNSS_lon( ii)/real(N_obs_in_tile(ii)) CYGNSS_lat( ii) = CYGNSS_lat( ii)/real(N_obs_in_tile(ii)) - CYGNSS_time( ii) = CYGNSS_time(ii)/real(N_obs_in_tile(ii),kind(0.0D0)) elseif (N_obs_in_tile(ii)==0) then @@ -2642,6 +2640,24 @@ subroutine read_obs_sm_CYGNSS( & end if ! clean up + + deallocate(timeintervals) + deallocate(latitudes) + deallocate(longitudes) + deallocate(tmp_sm) + deallocate(tmp_sigma) + deallocate(latitudes_m) + deallocate(longitudes_m) + deallocate(small_SM_range_flag) + deallocate(poor_SMAP_flag) + deallocate(high_ubrmsd_flag) + deallocate(few_obs_flag) + deallocate(low_signal_flag) + + if (allocated(sm_d)) deallocate(sm_d) + if (allocated(sm_subd)) deallocate(sm_subd) + if (allocated(sigma_d)) deallocate(sigma_d) + if (allocated(sigma_subd)) deallocate(sigma_subd) if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) @@ -9137,7 +9153,7 @@ subroutine read_obs( & end if - case ('CYGNSS_SM','CYGNSS_SM_daily') + case ('CYGNSS_SM_6hr','CYGNSS_SM_daily') call read_obs_sm_CYGNSS( & date_time, dtstep_assim, N_catd, tile_coord, & From a27412edd816384f9c4902d790685c452ffd2fc5 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Fri, 24 Jan 2025 20:16:05 -0500 Subject: [PATCH 14/21] bugfix in fname --- GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 31ed9a0..4d0d087 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2339,7 +2339,7 @@ subroutine read_obs_sm_CYGNSS( & tmpname = 'cyg.ddmi.s'// YYYY // MM // DD //'-030000-e'// YYYY // MM // DD //'-210000.l3.grid-soil-moisture-36km.a32.d33.nc' - tmpfname = trim(this_obs_param%path) // '/Y' // YYYY // '/M' // MM // '/' // trim(tmpname) // '.nc' + tmpfname = trim(this_obs_param%path) // '/Y' // YYYY // '/M' // MM // '/' // trim(tmpname) if (logit) write (logunit, '(400A)') 'Reading CYGNSS soil moisture data from file: ', trim(tmpfname) From 3638e93452312e2539f95b98b9ee489c7cb44ea0 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Mon, 27 Jan 2025 17:03:24 -0700 Subject: [PATCH 15/21] further improvements after review --- .../clsm_ensupd_read_obs.F90 | 171 ++++++++---------- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 12 +- 2 files changed, 78 insertions(+), 105 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 4d0d087..9770e0b 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2205,7 +2205,7 @@ subroutine read_obs_sm_CYGNSS( & ! local variables: - integer, parameter :: max_obs = 50000 ! max number of daily obs read by subroutine (expecting > 6 hr assim window) + integer, parameter :: max_obs = 50000 ! max number of daily obs read by subroutine (assim window <= 24 hr) character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' @@ -2217,26 +2217,26 @@ subroutine read_obs_sm_CYGNSS( & character( 2) :: MM, DD, HH, MI character( 4) :: YYYY - integer :: pos, i, idx, N_obs, j, HHcnt, ind, ii + integer :: i, idx, N_obs, j, ind, ii, dt_obs, obs_hour integer :: ierr, ncid integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid, lon_dimid_m, lat_dimid_m integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop, N_lon_m, N_lat_m integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid, few_obs_varid, low_signal_varid + integer :: start(3), count(3) integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) logical :: file_exists - logical :: HH03, HH09, HH12, HH15, HH21 real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_jtime(:) integer, dimension(:), pointer :: tmp_tile_num integer, dimension(N_catd) :: N_obs_in_tile - real, allocatable :: sm_d(:,:,:), sm_subd(:,:,:), sigma_d(:,:,:), sigma_subd(:,:,:) real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:), tmp_sm(:,:), tmp_sigma(:,:) real, allocatable :: time(:), timeslices(:,:) real, allocatable :: latitudes_m(:,:), longitudes_m(:,:) + real :: date_time_low_hour, date_time_up_hour real :: tmp1_lon(max_obs), tmp1_lat(max_obs), tmp1_obs(max_obs), tmp1_err(max_obs), tmp1_jtime(max_obs) ! ------------------------------------------------------------------- @@ -2247,12 +2247,6 @@ subroutine read_obs_sm_CYGNSS( & found_obs = .false. - HH03 = .false. - HH09 = .false. - HH12 = .false. - HH15 = .false. - HH21 = .false. - ! determine operating time range of sensor if (trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' .or. trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then @@ -2268,64 +2262,57 @@ subroutine read_obs_sm_CYGNSS( & if ( datetime_lt_refdatetime(date_time, date_time_obs_beg) .or. & datetime_lt_refdatetime(date_time_obs_end, date_time) ) return + ! determine if reading daily or subdaily obs + + if (trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then + dt_obs = 24 ! obs frequency in hours + elseif (trim(this_obs_param%descr) == 'CYGNSS_SM_6hr') then + dt_obs = 6 + else + err_msg = 'Unknown obs_param%descr: ' // trim(this_obs_param%descr) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + ! find the time bounds of the assimilation window ! (date_time-dtstep_assim/2,date_time+dtstep_assim/2] - date_time_low = date_time - call augment_date_time( -(dtstep_assim/2), date_time_low) - date_time_up = date_time - call augment_date_time( (dtstep_assim/2), date_time_up) - - ! establish if assimilation window encompasses 12:00 UTC (daily CYGNSS obs time) + date_time_low = date_time + call augment_date_time( -(dtstep_assim/2), date_time_low) + date_time_up = date_time + call augment_date_time( (dtstep_assim/2), date_time_up) - if (date_time_low%hour < 12 .and. date_time_up%hour >= 12) then - HH12 = .true. - end if - - ! return if assimilating daily obs and HH12 is false - - if (trim(this_obs_param%descr) == 'CYGNSS_SM_daily' .and. .not. HH12) return + ! account for minutes and seconds in date_time_low and date_time_up - ! cygnss sm subdaily observations are for 4 timeslices per day - 00-06, 06-12, 12-18, 18-24 UTC - ! we will assimilate these at 03, 09, 15, 21 UTC + date_time_low_hour = date_time_low%hour + date_time_low%min/60.0 + date_time_low%sec/3600.0 + date_time_up_hour = date_time_up%hour + date_time_up%min/60.0 + date_time_up%sec/3600.0 - HHcnt = 0 + ! deal with assimilation window that spans 0z - if (date_time_low%hour < 3 .and. date_time_up%hour >= 3) then - HH03 = .true. - HHcnt = HHcnt + 1 - end if - if (date_time_low%hour > 21 .and. date_time_up%hour >= 3) then ! wrap-around from previous day for 3 hr < assim window <= 6 hr - HH03 = .true. - HHcnt = HHcnt + 1 - end if - if (date_time_low%hour < 9 .and. date_time_up%hour >= 9) then - HH09 = .true. - HHcnt = HHcnt + 1 - end if - if (date_time_low%hour < 15 .and. date_time_up%hour >= 15) then - HH15 = .true. - HHcnt = HHcnt + 1 - end if - if (date_time_low%hour < 21 .and. date_time_up%hour >= 21) then - HH21 = .true. - HHcnt = HHcnt + 1 - end if - if (date_time_low%hour < 21 .and. date_time_up%hour < 3) then ! wrap-around into next day for 3 hr < assim window <= 6 hr - HH21 = .true. - HHcnt = HHcnt + 1 + if (date_time_low_hour > date_time_up_hour) then + date_time_low_hour = date_time_low_hour - 24.0 end if - ! return if assimilating subdaily obs and no HH is true - - if (trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' .and. .not. (HH03 .or. HH09 .or. HH15 .or. HH21)) return + obs_hour = -9999 + + if (dt_obs == 6) then ! subdaily obs + if (date_time_low_hour < 3.0 .and. date_time_up_hour >= 3.0) then + obs_hour = 3 + elseif (date_time_low_hour < 9.0 .and. date_time_up_hour >= 9.0) then + obs_hour = 9 + elseif (date_time_low_hour < 15.0 .and. date_time_up_hour >= 15.0) then + obs_hour = 15 + elseif (date_time_low_hour < 21.0 .and. date_time_up_hour >= 21.0) then + obs_hour = 21 + end if + else ! daily obs + if (date_time_low_hour < 12.0 .and. date_time_up_hour >= 12.0) then + obs_hour = 12 + end if + end if - ! if assimilation window encompasses more than one time slice, abort + ! return if assimilation window does not encompass any obs - if (HHcnt > 1) then - err_msg = 'Can not assimilate subdaily CYGNSS obs from more than one time slice' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - end if + if (obs_hour < 0) return ! determine filename for CYGNSS obs file @@ -2385,55 +2372,46 @@ subroutine read_obs_sm_CYGNSS( & allocate(tmp_sm( N_lon, N_lat )) allocate(tmp_sigma( N_lon, N_lat )) + ! define the count array for the NetCDF read operations + count = (/ N_lon, N_lat, 1 /) + ! read the variables ierr = nf90_get_var(ncid, timeintervals_varid, timeintervals) ierr = nf90_get_var(ncid, lat_varid, latitudes) ierr = nf90_get_var(ncid, lon_varid, longitudes) - ! read either subdaily or daily soil moisture and sigma variables - if ( trim(this_obs_param%descr) == 'CYGNSS_SM_6hr' ) then + idx = -1 - ! allocate memory for subdaily variables - allocate(sm_subd( N_lon, N_lat, N_timeslices )) - allocate(sigma_subd(N_lon, N_lat, N_timeslices )) - - ! subdaily observations required - ierr = nf90_get_var(ncid, sm_subd_varid, sm_subd) - ierr = nf90_get_var(ncid, sigma_subd_varid, sigma_subd) - - ! find the time slice that corresponds to the assimilation window - idx = -1 + ! read either subdaily or daily soil moisture and sigma variables + + if (dt_obs == 6) then do i = 1, N_timeslices - if ((HH03 .and. timeintervals(1,i) == 0.0 ) .or. & - (HH09 .and. timeintervals(1,i) == 0.25) .or. & - (HH15 .and. timeintervals(1,i) == 0.5 ) .or. & - (HH21 .and. timeintervals(1,i) == 0.75) ) then - idx = i - exit - end if + if ((obs_hour == 3 .and. timeintervals(1,i) == 0.0 ) .or. & + (obs_hour == 9 .and. timeintervals(1,i) == 0.25) .or. & + (obs_hour == 15 .and. timeintervals(1,i) == 0.5 ) .or. & + (obs_hour == 21 .and. timeintervals(1,i) == 0.75) ) then + idx = i + exit + end if end do - if (idx == -1) then - err_msg = 'Error: No matching time interval found for HH = ' // HH - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - else - write(*,*) 'Index of matching time interval: ', idx - end if - tmp_sm = sm_subd( :,:,idx) - tmp_sigma = sigma_subd(:,:,idx) - else + if (idx == -1) then + write(err_msg, '(A,I2)') 'Error: No matching time interval found for obs_hour = ', obs_hour + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + endif - ! allocate memory for daily variables - allocate(sm_d( N_lon, N_lat, N_time )) - allocate(sigma_d(N_lon, N_lat, N_time )) + start = (/ 1, 1, idx /) + ierr = nf90_get_var(ncid, sm_subd_varid, tmp_sm, start=start, count=count) + ierr = nf90_get_var(ncid, sigma_subd_varid, tmp_sigma, start=start, count=count) + + else ! daily obs + idx = 1 - ! daily observations required - ierr = nf90_get_var(ncid, sm_d_varid, sm_d) - ierr = nf90_get_var(ncid, sigma_d_varid, sigma_d) - tmp_sm = sm_d( :,:,1) - tmp_sigma = sigma_d(:,:,1) + start = (/ 1, 1, idx /) + ierr = nf90_get_var(ncid, sm_d_varid, tmp_sm, start=start, count=count) + ierr = nf90_get_var(ncid, sigma_d_varid, tmp_sigma, start=start, count=count) - end if + endif ! close the obs file ierr = nf90_close(ncid) @@ -2653,11 +2631,6 @@ subroutine read_obs_sm_CYGNSS( & deallocate(high_ubrmsd_flag) deallocate(few_obs_flag) deallocate(low_signal_flag) - - if (allocated(sm_d)) deallocate(sm_d) - if (allocated(sm_subd)) deallocate(sm_subd) - if (allocated(sigma_d)) deallocate(sigma_d) - if (allocated(sigma_subd)) deallocate(sigma_subd) if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index a4bfe27..c8d4b4b 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2365,11 +2365,11 @@ obs_param_nml(53)%adapt = 0 ! -------------------------------------------------------------------------------------------------------------- ! -! 54 = CYGNSS_SM (CYGNSS Level 3 soil moisture version 3.2) +! 54 = CYGNSS_SM_6hr (CYGNSS Level 3 soil moisture version 3.2) ! ! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 -obs_param_nml(54)%descr = 'CYGNSS_SM' +obs_param_nml(54)%descr = 'CYGNSS_SM_6hr' obs_param_nml(54)%orbit = 3 obs_param_nml(54)%pol = 0 obs_param_nml(54)%N_ang = 0 @@ -2394,10 +2394,10 @@ obs_param_nml(54)%scalepath = '' obs_param_nml(54)%scalename = '' obs_param_nml(54)%flistpath = '' obs_param_nml(54)%flistname = '' -obs_param_nml(54)%errstd = 0.1 +obs_param_nml(54)%errstd = 0.04 obs_param_nml(54)%std_normal_max = 2.5 obs_param_nml(54)%zeromean = .true. -obs_param_nml(54)%coarsen_pert = .false. +obs_param_nml(54)%coarsen_pert = .true. obs_param_nml(54)%xcorr = 0.25 obs_param_nml(54)%ycorr = 0.25 obs_param_nml(54)%adapt = 0 @@ -2433,10 +2433,10 @@ obs_param_nml(55)%scalepath = '' obs_param_nml(55)%scalename = '' obs_param_nml(55)%flistpath = '' obs_param_nml(55)%flistname = '' -obs_param_nml(55)%errstd = 0.1 +obs_param_nml(55)%errstd = 0.04 obs_param_nml(55)%std_normal_max = 2.5 obs_param_nml(55)%zeromean = .true. -obs_param_nml(55)%coarsen_pert = .false. +obs_param_nml(55)%coarsen_pert = .true. obs_param_nml(55)%xcorr = 0.25 obs_param_nml(55)%ycorr = 0.25 obs_param_nml(55)%adapt = 0 From a676245a22d5d241e85306c5446e812cded95ab2 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Thu, 30 Jan 2025 09:09:34 -0500 Subject: [PATCH 16/21] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a99df3..6ee77c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- New update_type for joint 3d soil moisture and 1d snow analysis (Tb+sfmc+sfds+SCF obs). - Added CYGNSS soil moisture reader. - Added M21C surface met forcing. From 0599e292c0c20a074b1c2e5d062f35b3a9230aad Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 30 Jan 2025 09:35:41 -0500 Subject: [PATCH 17/21] minimal cleanup of CYGNSS obs reader (clsm_ensupd_read_obs.F90) --- .../clsm_ensupd_read_obs.F90 | 29 ++++++++++--------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 9770e0b..2e2fd6f 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2222,7 +2222,9 @@ subroutine read_obs_sm_CYGNSS( & integer :: lon_dimid, lat_dimid, time_dimid, timeslices_dimid, startstop_dimid, lon_dimid_m, lat_dimid_m integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop, N_lon_m, N_lat_m - integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid, few_obs_varid, low_signal_varid + integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid + integer :: few_obs_varid, low_signal_varid + integer :: start(3), count(3) integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) @@ -2230,7 +2232,8 @@ subroutine read_obs_sm_CYGNSS( & real, dimension(:), pointer :: tmp_lon(:), tmp_lat(:), tmp_obs(:), tmp_err(:), tmp_jtime(:) integer, dimension(:), pointer :: tmp_tile_num - integer, dimension(N_catd) :: N_obs_in_tile + + integer, dimension(N_catd) :: N_obs_in_tile real, allocatable :: timeintervals(:,:), latitudes(:,:), longitudes(:,:), tmp_sm(:,:), tmp_sigma(:,:) real, allocatable :: time(:), timeslices(:,:) @@ -2243,7 +2246,7 @@ subroutine read_obs_sm_CYGNSS( & ! initialize - nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_jtime ) + nullify( tmp_lon, tmp_lat, tmp_obs, tmp_err, tmp_jtime, tmp_tile_num ) found_obs = .false. @@ -2281,7 +2284,7 @@ subroutine read_obs_sm_CYGNSS( & date_time_up = date_time call augment_date_time( (dtstep_assim/2), date_time_up) - ! account for minutes and seconds in date_time_low and date_time_up + ! get fractional time-of-day (hours) date_time_low_hour = date_time_low%hour + date_time_low%min/60.0 + date_time_low%sec/3600.0 date_time_up_hour = date_time_up%hour + date_time_up%min/60.0 + date_time_up%sec/3600.0 @@ -2295,17 +2298,17 @@ subroutine read_obs_sm_CYGNSS( & obs_hour = -9999 if (dt_obs == 6) then ! subdaily obs - if (date_time_low_hour < 3.0 .and. date_time_up_hour >= 3.0) then - obs_hour = 3 - elseif (date_time_low_hour < 9.0 .and. date_time_up_hour >= 9.0) then - obs_hour = 9 + if (date_time_low_hour < 3.0 .and. date_time_up_hour >= 3.0) then + obs_hour = 3 + elseif (date_time_low_hour < 9.0 .and. date_time_up_hour >= 9.0) then + obs_hour = 9 elseif (date_time_low_hour < 15.0 .and. date_time_up_hour >= 15.0) then obs_hour = 15 elseif (date_time_low_hour < 21.0 .and. date_time_up_hour >= 21.0) then obs_hour = 21 end if else ! daily obs - if (date_time_low_hour < 12.0 .and. date_time_up_hour >= 12.0) then + if (date_time_low_hour < 12.0 .and. date_time_up_hour >= 12.0) then obs_hour = 12 end if end if @@ -2405,11 +2408,10 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_get_var(ncid, sigma_subd_varid, tmp_sigma, start=start, count=count) else ! daily obs - idx = 1 - start = (/ 1, 1, idx /) - ierr = nf90_get_var(ncid, sm_d_varid, tmp_sm, start=start, count=count) - ierr = nf90_get_var(ncid, sigma_d_varid, tmp_sigma, start=start, count=count) + start = (/ 1, 1, 1 /) + ierr = nf90_get_var(ncid, sm_d_varid, tmp_sm, start=start, count=count) + ierr = nf90_get_var(ncid, sigma_d_varid, tmp_sigma, start=start, count=count) endif @@ -2448,7 +2450,6 @@ subroutine read_obs_sm_CYGNSS( & ierr = nf90_inq_varid(ncid, 'flag_few_obs', few_obs_varid) ierr = nf90_inq_varid(ncid, 'flag_low_signal', low_signal_varid) - ! allocate memory for the variables allocate(latitudes_m( N_lon_m, N_lat_m)) allocate(longitudes_m( N_lon_m, N_lat_m)) From d8665db54a0af61d1e8fa5171c2e8c3e7388608a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 30 Jan 2025 10:11:56 -0500 Subject: [PATCH 18/21] check that CYGNSS daily and subdaily sm obs are not both assimilated; check that dtstep_assim<=6h for CYGNSS sm obs (clsm_ensupd_read_obs.F90, clsm_ensupd_upd_routines.F90) --- .../clsm_ensupd_read_obs.F90 | 32 ++++++++++++------- .../clsm_ensupd_upd_routines.F90 | 23 +++++++++++++ 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 2e2fd6f..23675b3 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2156,19 +2156,20 @@ subroutine read_obs_sm_CYGNSS( & !--------------------------------------------------------------------- ! - ! Routine to read in CYGNSS soil moisture obs from netCDF files and apply mask. + ! Routine to read CYGNSS soil moisture obs from netCDF file and apply mask. ! https://podaac.jpl.nasa.gov/dataset/CYGNSS_L3_SOIL_MOISTURE_V3.2 - ! CYGNSS soil moisture obs are in volumetric units (m3 m-3) with a spatial - ! resolution of 0.3 Decimal Degrees x 0.37 Decimal Degrees. - ! They are daily files containing daily (obs_param_nml(.)%descr = CYGNSS_SM_daily) - ! and subdaily SM estimates (obs_param_nml(.)%descr = CYGNSS_SM_6hr). The - ! subdaily values are for 6 hr time intervals, 0-6hr, 6-12hr, 12-18hr, 18-24hr. - ! The subdaily values are assimilated at 3z, 9z, 15z, 21z. If the assimilation - ! window is more than 6 hr the code will error out. The daily values are assimilated - ! at 12z. The assumption is that either the daily or subdaily values are assimilated, - ! but the code will assimilate both if in the nml. - ! - ! A. Fox, Jan 2025 + ! CYGNSS soil moisture obs are in volumetric units (m3 m-3) on the + ! 36-km EASEv2 global cylindrical grid. + ! Obs are in files containing daily (obs_param_nml(.)%descr = CYGNSS_SM_daily) + ! and subdaily SM estimates (obs_param_nml(.)%descr = CYGNSS_SM_6hr). + ! Daily values are for 0-24z and assimilated at 12z. + ! Subdaily values are for 6 hr time intervals, 0-6z, 6-12z, 12-18z, and 18-24z + ! and assimilated at 3z, 9z, 15z, 21z. + ! Assimilation window must not be longer than 6 hours. + ! Subroutine read_ens_upd_inputs() ensures that only daily *or* subdaily + ! obs are assimilated. + ! + ! A. Fox, reichle, Jan 2025 ! ! -------------------------------------------------------------------- @@ -2265,6 +2266,13 @@ subroutine read_obs_sm_CYGNSS( & if ( datetime_lt_refdatetime(date_time, date_time_obs_beg) .or. & datetime_lt_refdatetime(date_time_obs_end, date_time) ) return + ! make sure assimilation window is no longer than 6 hours + + if (dtstep_assim>6*60*60) then + err_msg = 'dtstep_assim too long for CYGNSS soil moisture obs' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if + ! determine if reading daily or subdaily obs if (trim(this_obs_param%descr) == 'CYGNSS_SM_daily') then diff --git a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 index 0ab23b8..d078df4 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_upd_routines.F90 @@ -617,7 +617,30 @@ subroutine read_ens_upd_inputs( & end if + ! make sure only sub-daily or daily CYGNSS sm obs are assimilated + k=0 + + do i=1,N_obs_param + + if (obs_param(i)%assim) then + + select case (trim(obs_param(i)%descr)) + + case('CYGNSS_SM_daily','CYGNSS_SM_6hr'); k=k+1 + + end select + + end if + + end do + + if (k>1) then + + err_msg = 'cannot assimilate CYGNSS_SM_daily *and* CYGNSS_SM_6hr' + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + + end if ! ------------------------------------------------------------- ! From 220de4ec0b30308376d9d42c558a389c815f20df Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 30 Jan 2025 11:27:16 -0700 Subject: [PATCH 19/21] update obs file location --- GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml index c8d4b4b..b1c1764 100644 --- a/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml +++ b/GEOSldas_App/LDASsa_DEFAULT_inputs_ensupd.nml @@ -2386,7 +2386,7 @@ obs_param_nml(54)%bias_tcut = 432000 obs_param_nml(54)%nodata = -9999. obs_param_nml(54)%varname = 'sfmc' obs_param_nml(54)%units = 'm3 m-3' -obs_param_nml(54)%path = '/discover/nobackup/amfox/CYGNSS_obs/' +obs_param_nml(54)%path = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/CYGNSS/' obs_param_nml(54)%name = '' obs_param_nml(54)%maskpath = '' obs_param_nml(54)%maskname = '' @@ -2425,7 +2425,7 @@ obs_param_nml(55)%bias_tcut = 432000 obs_param_nml(55)%nodata = -9999. obs_param_nml(55)%varname = 'sfmc' obs_param_nml(55)%units = 'm3 m-3' -obs_param_nml(55)%path = '/discover/nobackup/amfox/CYGNSS_obs/' +obs_param_nml(55)%path = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/CYGNSS/' obs_param_nml(55)%name = '' obs_param_nml(55)%maskpath = '' obs_param_nml(55)%maskname = '' From e2747183fc6ab31675c121e8bd2ca48c374c9a50 Mon Sep 17 00:00:00 2001 From: amfox37 Date: Thu, 30 Jan 2025 11:29:15 -0700 Subject: [PATCH 20/21] update flag values and tolerance in float compare --- .../clsm_ensupd_read_obs.F90 | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 23675b3..98f9c74 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -2209,6 +2209,7 @@ subroutine read_obs_sm_CYGNSS( & integer, parameter :: max_obs = 50000 ! max number of daily obs read by subroutine (assim window <= 24 hr) character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 character(len=*), parameter :: Iam = 'read_obs_sm_CYGNSS' + real, parameter :: tolerance = 1.0e-6 ! tolerance for floating point comparisons type(date_time_type) :: date_time_obs_beg, date_time_obs_end type(date_time_type) :: date_time_up, date_time_low @@ -2224,7 +2225,7 @@ subroutine read_obs_sm_CYGNSS( & integer :: sm_d_varid, sm_subd_varid, sigma_d_varid, sigma_subd_varid, timeintervals_varid, lat_varid, lon_varid integer :: N_lon, N_lat, N_time, N_timeslices, N_startstop, N_lon_m, N_lat_m integer :: longitudes_m_varid, latitudes_m_varid, small_SM_range_varid, poor_SMAP_varid, high_ubrmsd_varid - integer :: few_obs_varid, low_signal_varid + integer :: few_obs_varid, low_signal_varid, good_flag_value integer :: start(3), count(3) integer, allocatable :: small_SM_range_flag(:,:), poor_SMAP_flag(:,:), high_ubrmsd_flag(:,:), few_obs_flag(:,:), low_signal_flag(:,:) @@ -2397,10 +2398,10 @@ subroutine read_obs_sm_CYGNSS( & if (dt_obs == 6) then do i = 1, N_timeslices - if ((obs_hour == 3 .and. timeintervals(1,i) == 0.0 ) .or. & - (obs_hour == 9 .and. timeintervals(1,i) == 0.25) .or. & - (obs_hour == 15 .and. timeintervals(1,i) == 0.5 ) .or. & - (obs_hour == 21 .and. timeintervals(1,i) == 0.75) ) then + if ((obs_hour == 3 .and. (abs(timeintervals(1,i) - 0.0 ) < tolerance)) .or. & + (obs_hour == 9 .and. (abs(timeintervals(1,i) - 0.25) < tolerance)) .or. & + (obs_hour == 15 .and. (abs(timeintervals(1,i) - 0.5 ) < tolerance)) .or. & + (obs_hour == 21 .and. (abs(timeintervals(1,i) - 0.75) < tolerance)) ) then idx = i exit end if @@ -2485,17 +2486,19 @@ subroutine read_obs_sm_CYGNSS( & call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if + good_flag_value = 255 ! should really be 0 but is 255 because of unsigned v. signed byte issues + ! fill tmp arrays N_obs = 0 do i = 1, N_lon do j = 1, N_lat - if (tmp_sm(i,j) .ne. this_obs_param%nodata .and. & - small_SM_range_flag(i,j) .ne. 1 .and. & - poor_SMAP_flag(i,j) .ne. 1 .and. & - high_ubrmsd_flag(i,j) .ne. 1 .and. & - few_obs_flag(i,j) .ne. 1 .and. & - low_signal_flag(i,j) .ne. 1 ) then + if (tmp_sm(i,j) .ne. this_obs_param%nodata .and. & + small_SM_range_flag(i,j) == good_flag_value .and. & + poor_SMAP_flag(i,j) == good_flag_value .and. & + high_ubrmsd_flag(i,j) == good_flag_value .and. & + few_obs_flag(i,j) == good_flag_value .and. & + low_signal_flag(i,j) == good_flag_value ) then ! valid observation N_obs = N_obs + 1 From 9a1d6edff9c3eca9bc31e5f2b56dcc8bc6e5050b Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 31 Jan 2025 11:31:23 -0500 Subject: [PATCH 21/21] for CYGNSS sm obs, added call to scale_obs_sfmc_zscore() (clsm_ensupd_read_obs.F90) --- .../clsm_ensupd_read_obs.F90 | 116 +++++++++--------- 1 file changed, 59 insertions(+), 57 deletions(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index 98f9c74..22ab792 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -9098,65 +9098,67 @@ subroutine read_obs( & case ('ASCAT_SM_A', 'ASCAT_SM_D' ) - - call read_obs_sm_ASCAT( & - work_path, exp_id, & - date_time, dtstep_assim, N_catd, tile_coord, & - tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & - this_obs_param, & - found_obs, tmp_obs, tmp_std_obs ) - - ! scale observations to model climatology - - if (this_obs_param%scale .and. found_obs) then - - scaled_obs = .true. - - call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & - tmp_obs, tmp_std_obs ) - - end if - + + call read_obs_sm_ASCAT( & + work_path, exp_id, & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, tmp_obs, tmp_std_obs ) + + ! scale observations to model climatology + + if (this_obs_param%scale .and. found_obs) then + + scaled_obs = .true. + + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & + tmp_obs, tmp_std_obs ) + + end if + case ('ASCAT_META_SM','ASCAT_METB_SM','ASCAT_METC_SM' ) - - call read_obs_sm_ASCAT_EUMET( & - date_time, dtstep_assim, N_catd, tile_coord, & - tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & - this_obs_param, & - found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & - tmp_time) - - ! scale observations to model climatology - - if (this_obs_param%scale .and. found_obs) then - - scaled_obs = .true. - - call scale_obs_sfmc_zscore( N_catd, tile_coord, & - date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & - tmp_obs, tmp_std_obs ) - - end if - + + call read_obs_sm_ASCAT_EUMET( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & + tmp_time) + + ! scale observations to model climatology + + if (this_obs_param%scale .and. found_obs) then + + scaled_obs = .true. + + call scale_obs_sfmc_zscore( N_catd, tile_coord, & + date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & + tmp_obs, tmp_std_obs ) + + end if + case ('CYGNSS_SM_6hr','CYGNSS_SM_daily') - - call read_obs_sm_CYGNSS( & - date_time, dtstep_assim, N_catd, tile_coord, & - tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & - this_obs_param, & - found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & - tmp_time) - - ! scale observations to model climatology - - if (this_obs_param%scale .and. found_obs) then - - scaled_obs = .true. - - ! TODO: implement scaling for CYGNSS obs - - end if - + + call read_obs_sm_CYGNSS( & + date_time, dtstep_assim, N_catd, tile_coord, & + tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & + this_obs_param, & + found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & + tmp_time) + + ! scale observations to model climatology + + if (this_obs_param%scale .and. found_obs) then + + scaled_obs = .true. + + call scale_obs_sfmc_zscore( N_catd, tile_coord, & + date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & + tmp_obs, tmp_std_obs ) + + end if + case ('isccp_tskin_gswp2_v1') call read_obs_isccp_tskin_gswp2_v1( &