From bfd360773ccf8e87627e6cc02ae2cda605b96061 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 2 Jan 2025 15:09:38 -0700 Subject: [PATCH 01/10] step 1: print out obs(k)%nobs_epoch, obs(k)%count_location_in_matching_file, diff --- base/MAPL_ObsUtil.F90 | 2 + .../Sampler/MAPL_TrajectoryMod_smod.F90 | 53 ++++++++++++++++--- 2 files changed, 48 insertions(+), 7 deletions(-) diff --git a/base/MAPL_ObsUtil.F90 b/base/MAPL_ObsUtil.F90 index 1b96fad92ae8..59af1aec46a2 100644 --- a/base/MAPL_ObsUtil.F90 +++ b/base/MAPL_ObsUtil.F90 @@ -26,6 +26,8 @@ module MAPL_ObsUtilMod type :: obs_unit integer :: nobs_epoch integer :: ngeoval + integer :: count_location_until_matching_file + integer :: count_location_in_matching_file logical :: export_all_geoval type(FileMetadata), allocatable :: metadata type(NetCDF4_FileFormatter), allocatable :: file_handle diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 4ae3193970c9..f12bd06c0983 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -540,8 +540,10 @@ real(REAL64) :: t_shift integer, allocatable :: obstype_id_full(:) integer, allocatable :: location_index_ioda_full(:) + integer, allocatable :: location_index_ioda_full_aux(:) integer, allocatable :: IA_full(:) + real(ESMF_KIND_R8), pointer :: ptAT(:) type(ESMF_routehandle) :: RH type(ESMF_Time) :: timeset(2) @@ -553,7 +555,7 @@ type(ESMF_VM) :: vm integer :: mypet, petcount, mpic - integer :: i, j, k, L, ii, jj + integer :: i, j, k, L, ii, jj, jj2, jj3 integer :: fid_s, fid_e integer(kind=ESMF_KIND_I8) :: j0, j1 integer(kind=ESMF_KIND_I8) :: jt1, jt2 @@ -612,7 +614,7 @@ trim(this%var_name_lat),trim(this%var_name_time)) L=0 - fid_s=this%obsfile_Ts_index + fid_s=this%obsfile_Ts_index ! this is downshifted by 1 in MAPL_ObsUtil.F90 fid_e=this%obsfile_Te_index call lgr%debug('%a %i10 %i10', & @@ -620,9 +622,14 @@ arr(1)=0 ! len_full if (mapl_am_I_root()) then + ! + ! __ s1. scan 1st time: determine len len = 0 do k=1, this%nobs_type j = max (fid_s, L) + jj2 = 0 ! obs(k) count location + this%obs(k)%count_location_until_matching_file = 0 ! init + this%obs(k)%count_location_in_matching_file = 0 ! init do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & @@ -632,6 +639,14 @@ call lgr%debug('%a %a', 'exist: true filename :', trim(filename)) call get_ncfile_dimension(filename, tdim=num_times, key_time=this%index_name_x, _RC) len = len + num_times + jj2 = jj2 + num_times + if (j==fid_s) then + this%obs(k)%count_location_until_matching_file = jj2 + write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_until_matching_file', this%obs(k)%count_location_until_matching_file + elseif (j==fid_s+1) then + this%obs(k)%count_location_in_matching_file = num_times + write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_in_matching_file', this%obs(k)%count_location_in_matching_file + end if else call lgr%debug('%a %i10', 'non-exist: filename fid j :', j) call lgr%debug('%a %a', 'non-exist: missing filename:', trim(filename)) @@ -641,6 +656,9 @@ enddo arr(1)=len + + ! + ! __ s2. scan 2nd time: read time ect. into array, set location_index if (len>0) then allocate(lons_full(len),lats_full(len),_STAT) allocate(times_R8_full(len),_STAT) @@ -652,6 +670,7 @@ ii = 0 do k=1, this%nobs_type j = max (fid_s, L) + jj2 = 0 ! obs(k) count location do while (j<=fid_e) filename = get_filename_from_template_use_index( & this%obsfile_start_time, this%obsfile_interval, & @@ -671,7 +690,10 @@ times_R8_full(len+1:len+num_times) = times_R8_full(len+1:len+num_times) + t_shift obstype_id_full(len+1:len+num_times) = k do jj = 1, num_times - location_index_ioda_full(len+jj) = jj + jj2 = jj2 + 1 + location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file +!! ygyu debug +!! write(6,'(2x,a,2x,i10)') 'location_index_ioda_full(len+jj) = ', location_index_ioda_full(len+jj) end do len = len + num_times end if @@ -681,6 +703,8 @@ end if end if + !call MPI_Barrier(mpic,ierr) + !stop 'nail 1' call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) @@ -708,10 +732,16 @@ call MAPL_CommsBcast(vm, this%datetime_units, N=ESMF_MAXSTR, ROOT=MAPL_Root, _RC) - if (mapl_am_I_root()) then call sort_index (times_R8_full, IA_full, _RC) - location_index_ioda_full(:) = IA_full(:) + !! ygyu: this is a complete reshuffle, but output the index of the sorted array + !! I need the content of location_index_ioda_full + allocate(location_index_ioda_full_aux, source=location_index_ioda_full, _STAT) + do jj = 1, nx_sum + location_index_ioda_full(jj) = location_index_ioda_full_aux(IA_full(jj)) + end do + deallocate(location_index_ioda_full_aux) + !! location_index_ioda_full(:) = IA_full(:) ! NVHPC dies with NVFORTRAN-S-0155-Could not resolve generic procedure sort_multi_arrays_by_time call sort_four_arrays_by_time(lons_full, lats_full, times_R8_full, obstype_id_full, _RC) call ESMF_ClockGet(this%clock,currTime=current_time,_RC) @@ -827,8 +857,17 @@ 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) do k=1, this%nobs_type - call lgr%debug('%a %i4 %a %i12', & - 'obs(', k, ')%nobs_epoch', this%obs(k)%nobs_epoch ) + call lgr%debug('%a %i4 %a %i12 %i12 %i12', & + 'obs(', k, ')%nobs_epoch, %count_location_in_matching_file, diff', & + this%obs(k)%nobs_epoch, this%obs(k)%count_location_in_matching_file, & + this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file ) + j = this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file + _ASSERT(j==0, 'count_location_in_matching_file diff from cutted nobs_epoch') + j = minval(this%obs(k)%location_index_ioda(:)) + _ASSERT(j==1, 'this%obs(k)%location_index_ioda(:) did not start from 1') + j = maxval(this%obs(k)%location_index_ioda(:)) + _ASSERT(j==this%obs(k)%count_location_in_matching_file, & + 'this%obs(k)%location_index_ioda(:) did not end at obs file max pts') enddo end if else From 19f7ead48a83a6d883b8c2ecb121e254b1a74661 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 3 Jan 2025 09:07:03 -0700 Subject: [PATCH 02/10] step 1.1 print: aircraft. missing pts in obervation --- .../Sampler/MAPL_TrajectoryMod_smod.F90 | 60 ++++++++++++++----- 1 file changed, 46 insertions(+), 14 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index f12bd06c0983..d2968d11329b 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -555,7 +555,7 @@ type(ESMF_VM) :: vm integer :: mypet, petcount, mpic - integer :: i, j, k, L, ii, jj, jj2, jj3 + integer :: i, j, k, L, ii, jj, jj2, kk integer :: fid_s, fid_e integer(kind=ESMF_KIND_I8) :: j0, j1 integer(kind=ESMF_KIND_I8) :: jt1, jt2 @@ -566,6 +566,7 @@ integer :: arr(1) integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch + integer, allocatable :: marker(:) integer :: nx2 logical :: EX ! file logical :: zero_obs @@ -783,15 +784,9 @@ arr(1)=nx else !! doulbe check - ! (x1, x2] design in bisect + ! (x1, x2] design in bisect : y(n) < x <= y(n+1), n is intercept index this%epoch_index(1)= jt1 + 1 -!! ! (x1, x2] design in bisect -!! if (jt1==0) then -!! this%epoch_index(1)= 1 -!! else -!! this%epoch_index(1)= jt1 -!! endif _ASSERT(jt2<=len, 'bisect index for this%epoch_index(2) failed') if (jt2==0) then this%epoch_index(2)= 1 @@ -857,17 +852,54 @@ 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) do k=1, this%nobs_type - call lgr%debug('%a %i4 %a %i12 %i12 %i12', & - 'obs(', k, ')%nobs_epoch, %count_location_in_matching_file, diff', & + call lgr%debug('%a %i4 %a80 %i12 %i12 %i12', & + 'obs(', k, ') '//trim(this%obs(k)%name)//' %nobs_epoch, %count_location_in_matching_file, diff', & this%obs(k)%nobs_epoch, this%obs(k)%count_location_in_matching_file, & this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file ) + !! ygyu: ioda file split [1/2 15Z : 1/2 21Z ] [ 1/2 21Z : 1/2 3Z] + !! now we know + !! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- + !! negative index (extra) at Tmin missing Tmax + + ! count non-positive index on the left Tmin (extra) + ! + jj = 0 + do j = 1, this%obs(k)%nobs_epoch + if ( this%obs(k)%location_index_ioda(j) <= 0 ) then + jj = jj + 1 + end if + end do + + kk = this%obs(k)%count_location_in_matching_file + allocate ( marker(kk), _STAT ) + marker(:) = 0 + do j = 1, this%obs(k)%nobs_epoch + kk = this%obs(k)%location_index_ioda(j) + if ( kk > 0 .and. kk <= this%obs(k)%nobs_epoch ) then + marker(kk) = 1 ! exist + end if + end do + do kk = 1, this%obs(k)%count_location_in_matching_file + if ( marker(kk) == 0 ) then + write(6, '(2x,a,2x,10i12)') trim(this%obs(k)%name)//' missing pts in obervation ', kk + end if + end do + deallocate(marker) + j = this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file - _ASSERT(j==0, 'count_location_in_matching_file diff from cutted nobs_epoch') + if (j/=0) & + write(6, *) trim(this%obs(k)%name)//'count_location_in_matching_file diff from cutted nobs_epoch' + !!_ASSERT(j==0, trim(this%obs(k)%name)//'count_location_in_matching_file diff from cutted nobs_epoch') j = minval(this%obs(k)%location_index_ioda(:)) - _ASSERT(j==1, 'this%obs(k)%location_index_ioda(:) did not start from 1') + if (j/=1) & + write(6,*) trim(this%obs(k)%name)//' min =', j, & + ' this%obs(k)%location_index_ioda(:) start diff from 1' + !!_ASSERT(j==1, trim(this%obs(k)%name)//'this%obs(k)%location_index_ioda(:) did not start from 1') j = maxval(this%obs(k)%location_index_ioda(:)) - _ASSERT(j==this%obs(k)%count_location_in_matching_file, & - 'this%obs(k)%location_index_ioda(:) did not end at obs file max pts') + if (j/=this%obs(k)%count_location_in_matching_file) & + write(6,*) trim(this%obs(k)%name)//' max =', j, & + ' this%obs(k)%location_index_ioda(:) end diff from obs file max pts' + !!_ASSERT(j==this%obs(k)%count_location_in_matching_file, trim(this%obs(k)%name)//'this%obs(k)%location_index_ioda(:) did not end at obs file max pts') enddo end if else From bf55de118026ff0372b908c3828e9a48de3836f0 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 3 Jan 2025 14:37:50 -0700 Subject: [PATCH 03/10] Add if (traj%use_NWP_1_file == 1) option --- gridcomps/History/MAPL_HistoryGridComp.F90 | 3 - .../History/Sampler/MAPL_TrajectoryMod.F90 | 3 + .../Sampler/MAPL_TrajectoryMod_smod.F90 | 60 +++++++++++++++---- 3 files changed, 51 insertions(+), 15 deletions(-) diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index 17deb2db1299..be6e2a361e91 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -5384,8 +5384,6 @@ function get_acc_offset(current_time,ref_time,rc) result(acc_offset) ! __ for each collection: find union fields, write to collection.rcx ! __ note: this subroutine is called by MPI root only ! - ! __ note: this subroutine is called by MPI root only - ! subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) use MAPL_scan_pattern_in_file use MAPL_ObsUtilMod, only : obs_platform, union_platform @@ -5444,7 +5442,6 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) length_mx = ESMF_MAXSTR2 mxseg = 100 - ! __ s1. scan get platform name + index_name_x var_name_lat/lon/time do k=1, nplf call scan_begin(unitr, 'PLATFORM.', .false.) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 49c2eff96b38..f66ecf3747b7 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -40,6 +40,7 @@ module HistoryTrajectoryMod type(VerticalData) :: vdata logical :: do_vertical_regrid + type(LocstreamRegridder) :: regridder type(TimeData) :: time_info type(ESMF_Clock) :: clock @@ -47,6 +48,7 @@ module HistoryTrajectoryMod type(ESMF_Time) :: RingTime type(ESMF_TimeInterval), public :: epoch_frequency + integer :: nobs_type ! character(len=ESMF_MAXSTR) :: nc_index ! character(len=ESMF_MAXSTR) :: nc_time @@ -62,6 +64,7 @@ module HistoryTrajectoryMod character(len=ESMF_MAXSTR) :: var_name_lon_full character(len=ESMF_MAXSTR) :: datetime_units character(len=ESMF_MAXSTR) :: Location_index_name + integer :: use_NWP_1_file integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index d2968d11329b..9f7c0d6dd82a 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -52,6 +52,7 @@ traj%clock=clock if (present(GENSTATE)) traj%GENSTATE => GENSTATE + call ESMF_ClockGet ( clock, CurrTime=currTime, _RC ) call ESMF_ConfigGetAttribute(config, value=time_integer, label=trim(string)//'Epoch:', default=0, _RC) _ASSERT(time_integer /= 0, 'Epoch value in config wrong') @@ -71,6 +72,15 @@ label=trim(string) // 'var_name_lat:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_time_full, default="", & label=trim(string) // 'var_name_time:', _RC) + call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=0, & + label=trim(string)//'use_NWP_1_file:', _RC) + if (mapl_am_I_root()) then + if (traj%use_NWP_1_file == 1) then + write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is ON' + write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' + end if + _ASSERT ( (traj%use_NWP_1_file==1).OR.(traj%use_NWP_1_file==0), 'use_NWP_1_file: wrong input') + end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_begin:', _RC) @@ -615,8 +625,15 @@ trim(this%var_name_lat),trim(this%var_name_time)) L=0 - fid_s=this%obsfile_Ts_index ! this is downshifted by 1 in MAPL_ObsUtil.F90 - fid_e=this%obsfile_Te_index + if (this%use_NWP_1_file == 1) then + ! NWP IODA 1 file case + fid_s=this%obsfile_Ts_index+1 ! this is downshifted by 1 in MAPL_ObsUtil.F90 + fid_e=fid_s + else + ! regular case for any trajectory + fid_s=this%obsfile_Ts_index ! this is downshifted by 1 in MAPL_ObsUtil.F90 + fid_e=this%obsfile_Te_index + end if call lgr%debug('%a %i10 %i10', & 'fid_s, fid_e', fid_s, fid_e) @@ -641,10 +658,10 @@ call get_ncfile_dimension(filename, tdim=num_times, key_time=this%index_name_x, _RC) len = len + num_times jj2 = jj2 + num_times - if (j==fid_s) then + if (j==this%obsfile_Ts_index) then this%obs(k)%count_location_until_matching_file = jj2 write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_until_matching_file', this%obs(k)%count_location_until_matching_file - elseif (j==fid_s+1) then + elseif (j==this%obsfile_Ts_index+1) then this%obs(k)%count_location_in_matching_file = num_times write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_in_matching_file', this%obs(k)%count_location_in_matching_file end if @@ -704,8 +721,9 @@ end if end if - !call MPI_Barrier(mpic,ierr) - !stop 'nail 1' +!! ygyu debug +! call MPI_Barrier(mpic,ierr) +! stop 'nail 1' call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) @@ -752,7 +770,15 @@ sec = hms_2_s(this%Epoch) j1 = j0 + int(sec, kind=ESMF_KIND_I8) jx0 = real ( j0, kind=ESMF_KIND_R8) - jx1 = real ( j1, kind=ESMF_KIND_R8) + if (this%use_NWP_1_file == 1) then + ! IODA case: Epoch + 1 second : for ioda location index recover to match observation files + ! this works because all IODA files has same end <= Epoch + ! and we use 1 file only + jx1 = real ( j1 + 1, kind=ESMF_KIND_R8) + else + ! normal case: strict design with cutoff time at Epoch + jx1 = real ( j1, kind=ESMF_KIND_R8) + end if nstart=1; nend=size(times_R8_full) call bisect( times_R8_full, jx0, jt1, n_LB=int(nstart, ESMF_KIND_I8), n_UB=int(nend, ESMF_KIND_I8), rc=rc) @@ -992,6 +1018,8 @@ ! defer destroy fieldB at regen_grid step ! + ! debug + !write(6,*) 'ip, npt=', ip, size(this%obsTime) _RETURN(_SUCCESS) end procedure create_grid @@ -1270,7 +1298,7 @@ integer :: x_subset(2) type(ESMF_Time) :: timeset(2) type(ESMF_Time) :: current_time - type(ESMF_TimeInterval) :: dur + type(ESMF_TimeInterval) :: dur, delT type(GriddedIOitemVectorIterator) :: iter type(GriddedIOitem), pointer :: item type(ESMF_Field) :: src_field,dst_field,acc_field @@ -1303,11 +1331,19 @@ call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) timeset(1) = current_time - dur timeset(2) = current_time + if (this%use_NWP_1_file == 1) then + ! + ! change UB to Epoch + 1 s to be inclusive for IODA + if ( ESMF_AlarmIsRinging (this%alarm) ) then + call ESMF_TimeIntervalSet(delT, s=1, _RC) + timeset(2) = current_time + delT + end if + end if call this%get_x_subset(timeset, x_subset, _RC) is=x_subset(1) ie=x_subset(2) - !! write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie - + ! debug + !write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie ! ! __ I designed a method to return from regridding if no valid points exist @@ -1498,8 +1534,8 @@ call bisect( this%obstime, rT1, index1, n_LB=lb, n_UB=ub, rc=rc) call bisect( this%obstime, rT2, index2, n_LB=lb, n_UB=ub, rc=rc) - ! (x1, x2] design in bisect - ! simple version + ! (x1, x2] design in bisect: y(n) < x <= y(n+1), n is output index + ! simple version x_subset(1) = index1+1 x_subset(2) = index2 From 93794403a628d6bf6d7c7c17e2d77422f16ae7fb Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 3 Jan 2025 20:42:11 -0700 Subject: [PATCH 04/10] NWP_1_file complete: # of obs points recover, no missing points --- gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 9f7c0d6dd82a..045135f192bb 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -886,7 +886,7 @@ !! now we know !! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- !! negative index (extra) at Tmin missing Tmax - + ! ! count non-positive index on the left Tmin (extra) ! jj = 0 @@ -915,17 +915,14 @@ j = this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file if (j/=0) & write(6, *) trim(this%obs(k)%name)//'count_location_in_matching_file diff from cutted nobs_epoch' - !!_ASSERT(j==0, trim(this%obs(k)%name)//'count_location_in_matching_file diff from cutted nobs_epoch') j = minval(this%obs(k)%location_index_ioda(:)) if (j/=1) & write(6,*) trim(this%obs(k)%name)//' min =', j, & ' this%obs(k)%location_index_ioda(:) start diff from 1' - !!_ASSERT(j==1, trim(this%obs(k)%name)//'this%obs(k)%location_index_ioda(:) did not start from 1') j = maxval(this%obs(k)%location_index_ioda(:)) if (j/=this%obs(k)%count_location_in_matching_file) & write(6,*) trim(this%obs(k)%name)//' max =', j, & ' this%obs(k)%location_index_ioda(:) end diff from obs file max pts' - !!_ASSERT(j==this%obs(k)%count_location_in_matching_file, trim(this%obs(k)%name)//'this%obs(k)%location_index_ioda(:) did not end at obs file max pts') enddo end if else From 80185d6a1dba104cf619e26e886c4a3d47cc7afa Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 3 Jan 2025 21:06:23 -0700 Subject: [PATCH 05/10] Introduce use_NWP_1_file=1 option to trajectory sampler, in order that location_index for IODA file exact match the single input NPW observation file default use_NWP_1_file=0 recovers the strict Time interval cut stratgy as python [T1, T1+Epoch). --- .../History/Sampler/MAPL_TrajectoryMod.F90 | 7 ----- .../Sampler/MAPL_TrajectoryMod_smod.F90 | 31 ++++++++----------- 2 files changed, 13 insertions(+), 25 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index f66ecf3747b7..451a0bb6be74 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -40,7 +40,6 @@ module HistoryTrajectoryMod type(VerticalData) :: vdata logical :: do_vertical_regrid - type(LocstreamRegridder) :: regridder type(TimeData) :: time_info type(ESMF_Clock) :: clock @@ -48,13 +47,7 @@ module HistoryTrajectoryMod type(ESMF_Time) :: RingTime type(ESMF_TimeInterval), public :: epoch_frequency - integer :: nobs_type -! character(len=ESMF_MAXSTR) :: nc_index -! character(len=ESMF_MAXSTR) :: nc_time -! character(len=ESMF_MAXSTR) :: nc_latitude -! character(len=ESMF_MAXSTR) :: nc_longitude - character(len=ESMF_MAXSTR) :: index_name_x character(len=ESMF_MAXSTR) :: var_name_time character(len=ESMF_MAXSTR) :: var_name_lat diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 045135f192bb..1737a03d8426 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -710,8 +710,7 @@ do jj = 1, num_times jj2 = jj2 + 1 location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file -!! ygyu debug -!! write(6,'(2x,a,2x,i10)') 'location_index_ioda_full(len+jj) = ', location_index_ioda_full(len+jj) + !! write(6,'(2x,a,2x,i10)') 'location_index_ioda_full(len+jj) = ', location_index_ioda_full(len+jj) end do len = len + num_times end if @@ -721,10 +720,6 @@ end if end if -!! ygyu debug -! call MPI_Barrier(mpic,ierr) -! stop 'nail 1' - call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, & count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) if (nx_sum == 0) then @@ -753,14 +748,12 @@ if (mapl_am_I_root()) then call sort_index (times_R8_full, IA_full, _RC) - !! ygyu: this is a complete reshuffle, but output the index of the sorted array - !! I need the content of location_index_ioda_full + !! use index to sort togehter multiple arrays allocate(location_index_ioda_full_aux, source=location_index_ioda_full, _STAT) do jj = 1, nx_sum location_index_ioda_full(jj) = location_index_ioda_full_aux(IA_full(jj)) end do deallocate(location_index_ioda_full_aux) - !! location_index_ioda_full(:) = IA_full(:) ! NVHPC dies with NVFORTRAN-S-0155-Could not resolve generic procedure sort_multi_arrays_by_time call sort_four_arrays_by_time(lons_full, lats_full, times_R8_full, obstype_id_full, _RC) call ESMF_ClockGet(this%clock,currTime=current_time,_RC) @@ -878,16 +871,18 @@ 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) do k=1, this%nobs_type - call lgr%debug('%a %i4 %a80 %i12 %i12 %i12', & - 'obs(', k, ') '//trim(this%obs(k)%name)//' %nobs_epoch, %count_location_in_matching_file, diff', & + call lgr%debug('%a %i4 %a %a80 %i12 %i12 %i12', & + 'obs(', k, ') ', trim(this%obs(k)%name)//' %nobs_epoch, %count_location_in_matching_file, diff', & this%obs(k)%nobs_epoch, this%obs(k)%count_location_in_matching_file, & this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file ) - !! ygyu: ioda file split [1/2 15Z : 1/2 21Z ] [ 1/2 21Z : 1/2 3Z] - !! now we know - !! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- - !! negative index (extra) at Tmin missing Tmax ! - ! count non-positive index on the left Tmin (extra) + ! diagnostic print + ! now we know + ! ioda file split [1/2 15Z : 1/2 21Z ] [ 1/2 21Z : 1/2 3Z] + ! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- + ! negative index (extra) at Tmin missing Tmax + ! + ! count non-positive index on the left Tmin (extra) ! jj = 0 do j = 1, this%obs(k)%nobs_epoch @@ -1016,7 +1011,7 @@ ! ! debug - !write(6,*) 'ip, npt=', ip, size(this%obsTime) + ! write(6,*) 'ip, npt=', ip, size(this%obsTime) _RETURN(_SUCCESS) end procedure create_grid @@ -1340,7 +1335,7 @@ is=x_subset(1) ie=x_subset(2) ! debug - !write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie + ! write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie ! ! __ I designed a method to return from regridding if no valid points exist From 5308922dd2540136190d7c6cd7bdb1f528b9e823 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Fri, 3 Jan 2025 21:45:34 -0700 Subject: [PATCH 06/10] Add item in CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1de3b647a206..730f1b2f3996 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 loggers when writing or reading weight files - Added new option to AGCM.rc `overwrite_checkpoint` to allow checkpoint files to be overwritten. By default still will not overwrite checkpoints +- Added use_NWP_1_file=1 option in Trajectory sampler (HISTORY.rc) to generate the correct location_index_ioda array, mapping back to location points in the single IODA observation input file ### Changed From 3ad76f7bde25779f41b70c0f8211c326e4b14d32 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Mon, 6 Jan 2025 12:34:22 -0700 Subject: [PATCH 07/10] code cleanup; add clear notes --- .../History/Sampler/MAPL_TrajectoryMod.F90 | 1 + .../Sampler/MAPL_TrajectoryMod_smod.F90 | 89 +++++-------------- 2 files changed, 23 insertions(+), 67 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 451a0bb6be74..6157dd17fbf4 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -58,6 +58,7 @@ module HistoryTrajectoryMod character(len=ESMF_MAXSTR) :: datetime_units character(len=ESMF_MAXSTR) :: Location_index_name integer :: use_NWP_1_file + integer, dimension(2), parameter :: use_NWP_1_file_param = [0, 1] integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 1737a03d8426..7e8cb58d172b 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -75,11 +75,12 @@ call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=0, & label=trim(string)//'use_NWP_1_file:', _RC) if (mapl_am_I_root()) then + _ASSERT ( ANY(traj%use_NPW_1_file_param == traj%use_NWP_1_file), & + 'use_NWP_1_file: wrong input') if (traj%use_NWP_1_file == 1) then write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is ON' write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' end if - _ASSERT ( (traj%use_NWP_1_file==1).OR.(traj%use_NWP_1_file==0), 'use_NWP_1_file: wrong input') end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & @@ -627,11 +628,11 @@ L=0 if (this%use_NWP_1_file == 1) then ! NWP IODA 1 file case - fid_s=this%obsfile_Ts_index+1 ! this is downshifted by 1 in MAPL_ObsUtil.F90 + fid_s=this%obsfile_Ts_index+1 ! index is downshifted by 1 in MAPL_ObsUtil.F90 fid_e=fid_s else ! regular case for any trajectory - fid_s=this%obsfile_Ts_index ! this is downshifted by 1 in MAPL_ObsUtil.F90 + fid_s=this%obsfile_Ts_index ! index is downshifted by 1 in MAPL_ObsUtil.F90 fid_e=this%obsfile_Te_index end if @@ -660,10 +661,8 @@ jj2 = jj2 + num_times if (j==this%obsfile_Ts_index) then this%obs(k)%count_location_until_matching_file = jj2 - write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_until_matching_file', this%obs(k)%count_location_until_matching_file elseif (j==this%obsfile_Ts_index+1) then this%obs(k)%count_location_in_matching_file = num_times - write(6,'(2x,a,2x,i10)') 'this%obs(k)%count_location_in_matching_file', this%obs(k)%count_location_in_matching_file end if else call lgr%debug('%a %i10', 'non-exist: filename fid j :', j) @@ -674,9 +673,9 @@ enddo arr(1)=len - ! - ! __ s2. scan 2nd time: read time ect. into array, set location_index + ! __ s2. scan 2nd time: read time ect. into array, + ! set location_index starting with the 1st element of the closest matched obs file if (len>0) then allocate(lons_full(len),lats_full(len),_STAT) allocate(times_R8_full(len),_STAT) @@ -709,8 +708,8 @@ obstype_id_full(len+1:len+num_times) = k do jj = 1, num_times jj2 = jj2 + 1 + ! for each obs type use the correct starting point location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file - !! write(6,'(2x,a,2x,i10)') 'location_index_ioda_full(len+jj) = ', location_index_ioda_full(len+jj) end do len = len + num_times end if @@ -751,7 +750,8 @@ !! use index to sort togehter multiple arrays allocate(location_index_ioda_full_aux, source=location_index_ioda_full, _STAT) do jj = 1, nx_sum - location_index_ioda_full(jj) = location_index_ioda_full_aux(IA_full(jj)) + ii = IA_full(jj) + location_index_ioda_full(jj) = location_index_ioda_full_aux(ii) end do deallocate(location_index_ioda_full_aux) ! NVHPC dies with NVFORTRAN-S-0155-Could not resolve generic procedure sort_multi_arrays_by_time @@ -764,12 +764,12 @@ j1 = j0 + int(sec, kind=ESMF_KIND_I8) jx0 = real ( j0, kind=ESMF_KIND_R8) if (this%use_NWP_1_file == 1) then - ! IODA case: Epoch + 1 second : for ioda location index recover to match observation files - ! this works because all IODA files has same end <= Epoch - ! and we use 1 file only + ! IODA case: + ! Upper bound time is set at 'Epoch + 1 second' to get the right index from bisect + ! jx1 = real ( j1 + 1, kind=ESMF_KIND_R8) else - ! normal case: strict design with cutoff time at Epoch + ! normal case: jx1 = real ( j1, kind=ESMF_KIND_R8) end if @@ -870,55 +870,15 @@ call lgr%debug('%a %i12 %i12 %i12', & 'epoch_index(1:2), nx', this%epoch_index(1), & this%epoch_index(2), this%nobs_epoch) - do k=1, this%nobs_type - call lgr%debug('%a %i4 %a %a80 %i12 %i12 %i12', & - 'obs(', k, ') ', trim(this%obs(k)%name)//' %nobs_epoch, %count_location_in_matching_file, diff', & - this%obs(k)%nobs_epoch, this%obs(k)%count_location_in_matching_file, & - this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file ) - ! - ! diagnostic print - ! now we know - ! ioda file split [1/2 15Z : 1/2 21Z ] [ 1/2 21Z : 1/2 3Z] - ! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- - ! negative index (extra) at Tmin missing Tmax - ! - ! count non-positive index on the left Tmin (extra) - ! - jj = 0 - do j = 1, this%obs(k)%nobs_epoch - if ( this%obs(k)%location_index_ioda(j) <= 0 ) then - jj = jj + 1 - end if - end do - - kk = this%obs(k)%count_location_in_matching_file - allocate ( marker(kk), _STAT ) - marker(:) = 0 - do j = 1, this%obs(k)%nobs_epoch - kk = this%obs(k)%location_index_ioda(j) - if ( kk > 0 .and. kk <= this%obs(k)%nobs_epoch ) then - marker(kk) = 1 ! exist - end if - end do - do kk = 1, this%obs(k)%count_location_in_matching_file - if ( marker(kk) == 0 ) then - write(6, '(2x,a,2x,10i12)') trim(this%obs(k)%name)//' missing pts in obervation ', kk - end if - end do - deallocate(marker) - - j = this%obs(k)%nobs_epoch - this%obs(k)%count_location_in_matching_file - if (j/=0) & - write(6, *) trim(this%obs(k)%name)//'count_location_in_matching_file diff from cutted nobs_epoch' - j = minval(this%obs(k)%location_index_ioda(:)) - if (j/=1) & - write(6,*) trim(this%obs(k)%name)//' min =', j, & - ' this%obs(k)%location_index_ioda(:) start diff from 1' - j = maxval(this%obs(k)%location_index_ioda(:)) - if (j/=this%obs(k)%count_location_in_matching_file) & - write(6,*) trim(this%obs(k)%name)//' max =', j, & - ' this%obs(k)%location_index_ioda(:) end diff from obs file max pts' - enddo + ! + ! Note: For IODA files, the default NPW_1_file=0 case shall reveal + ! the non-python array behavior in obs time sequence from observation files: for example: + ! ioda file split [1/2 15Z : 1/2 21Z ] [ 1/2 21Z : 1/2 3Z] (aircraft) + ! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- + ! negative index (extra) at Tmin missing Tmax + ! our trajectory ioda_index will show: overcount at Tmin and missing points at Tmax + ! NPW_1_file=1 solves this issue. + ! end if else allocate(this%lons(0),this%lats(0),_STAT) @@ -1010,8 +970,6 @@ ! defer destroy fieldB at regen_grid step ! - ! debug - ! write(6,*) 'ip, npt=', ip, size(this%obsTime) _RETURN(_SUCCESS) end procedure create_grid @@ -1334,8 +1292,6 @@ call this%get_x_subset(timeset, x_subset, _RC) is=x_subset(1) ie=x_subset(2) - ! debug - ! write(6,'(2x,a,4i10)') 'in regrid_accumulate is, ie=', is, ie ! ! __ I designed a method to return from regridding if no valid points exist @@ -1527,7 +1483,6 @@ call bisect( this%obstime, rT2, index2, n_LB=lb, n_UB=ub, rc=rc) ! (x1, x2] design in bisect: y(n) < x <= y(n+1), n is output index - ! simple version x_subset(1) = index1+1 x_subset(2) = index2 From 127075c62f170320b404a52f3177edd58a6b7cbe Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Mon, 6 Jan 2025 13:24:41 -0700 Subject: [PATCH 08/10] Add integer parameter: use_NWP_1_file_param --- gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 | 1 - gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 6157dd17fbf4..451a0bb6be74 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -58,7 +58,6 @@ module HistoryTrajectoryMod character(len=ESMF_MAXSTR) :: datetime_units character(len=ESMF_MAXSTR) :: Location_index_name integer :: use_NWP_1_file - integer, dimension(2), parameter :: use_NWP_1_file_param = [0, 1] integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 7e8cb58d172b..cfcf05f4f0ee 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -25,7 +25,7 @@ use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_fortran_env, only: INT64 implicit none - + integer, parameter :: use_NWP_1_file_param(2) = [0, 1] contains module procedure HistoryTrajectory_from_config @@ -75,7 +75,7 @@ call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=0, & label=trim(string)//'use_NWP_1_file:', _RC) if (mapl_am_I_root()) then - _ASSERT ( ANY(traj%use_NPW_1_file_param == traj%use_NWP_1_file), & + _ASSERT (ANY(use_NWP_1_file_param == traj%use_NWP_1_file), & 'use_NWP_1_file: wrong input') if (traj%use_NWP_1_file == 1) then write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is ON' From 8b633dc50c5563adb07144919d67fb207da83643 Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Mon, 6 Jan 2025 14:32:37 -0700 Subject: [PATCH 09/10] use one line for _ASSERT() --- gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index cfcf05f4f0ee..d0a818ba79b7 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -75,8 +75,7 @@ call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=0, & label=trim(string)//'use_NWP_1_file:', _RC) if (mapl_am_I_root()) then - _ASSERT (ANY(use_NWP_1_file_param == traj%use_NWP_1_file), & - 'use_NWP_1_file: wrong input') + _ASSERT (ANY(use_NWP_1_file_param == traj%use_NWP_1_file), 'use_NWP_1_file: wrong input') if (traj%use_NWP_1_file == 1) then write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is ON' write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' From 46d77c4dab72113ee294102338e34289358479be Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Tue, 7 Jan 2025 10:38:12 -0700 Subject: [PATCH 10/10] change variable use_NWP_1_file from integer to logical --- gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 | 2 +- .../History/Sampler/MAPL_TrajectoryMod_smod.F90 | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 451a0bb6be74..67eee93f452a 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -57,7 +57,7 @@ module HistoryTrajectoryMod character(len=ESMF_MAXSTR) :: var_name_lon_full character(len=ESMF_MAXSTR) :: datetime_units character(len=ESMF_MAXSTR) :: Location_index_name - integer :: use_NWP_1_file + logical :: use_NWP_1_file = .false. integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index d0a818ba79b7..924ac99cc7ac 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -25,7 +25,6 @@ use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_fortran_env, only: INT64 implicit none - integer, parameter :: use_NWP_1_file_param(2) = [0, 1] contains module procedure HistoryTrajectory_from_config @@ -72,11 +71,10 @@ label=trim(string) // 'var_name_lat:', _RC) call ESMF_ConfigGetAttribute(config, value=traj%var_name_time_full, default="", & label=trim(string) // 'var_name_time:', _RC) - call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=0, & + call ESMF_ConfigGetAttribute(config, value=traj%use_NWP_1_file, default=.false., & label=trim(string)//'use_NWP_1_file:', _RC) if (mapl_am_I_root()) then - _ASSERT (ANY(use_NWP_1_file_param == traj%use_NWP_1_file), 'use_NWP_1_file: wrong input') - if (traj%use_NWP_1_file == 1) then + if (traj%use_NWP_1_file) then write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is ON' write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' end if @@ -576,7 +574,7 @@ integer :: arr(1) integer :: sec integer, allocatable :: ix(:) ! counter for each obs(k)%nobs_epoch - integer, allocatable :: marker(:) + integer, allocatable :: marker(:) integer :: nx2 logical :: EX ! file logical :: zero_obs @@ -625,7 +623,7 @@ trim(this%var_name_lat),trim(this%var_name_time)) L=0 - if (this%use_NWP_1_file == 1) then + if (this%use_NWP_1_file) then ! NWP IODA 1 file case fid_s=this%obsfile_Ts_index+1 ! index is downshifted by 1 in MAPL_ObsUtil.F90 fid_e=fid_s @@ -762,10 +760,10 @@ sec = hms_2_s(this%Epoch) j1 = j0 + int(sec, kind=ESMF_KIND_I8) jx0 = real ( j0, kind=ESMF_KIND_R8) - if (this%use_NWP_1_file == 1) then + if (this%use_NWP_1_file) then ! IODA case: ! Upper bound time is set at 'Epoch + 1 second' to get the right index from bisect - ! + ! jx1 = real ( j1 + 1, kind=ESMF_KIND_R8) else ! normal case: @@ -1280,7 +1278,7 @@ call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) timeset(1) = current_time - dur timeset(2) = current_time - if (this%use_NWP_1_file == 1) then + if (this%use_NWP_1_file) then ! ! change UB to Epoch + 1 s to be inclusive for IODA if ( ESMF_AlarmIsRinging (this%alarm) ) then