Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revise location_index for IODA file in Trajectory sampler #3279

Open
wants to merge 13 commits into
base: develop
Choose a base branch
from
Open

Revise location_index for IODA file in Trajectory sampler #3279

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
bfd3607
step 1: print out obs(k)%nobs_epoch, obs(k)%count_location_in_matchin…
metdyn Jan 2, 2025
19f7ead
step 1.1 print: aircraft. missing pts in obervation
metdyn Jan 3, 2025
bf55de1
Add if (traj%use_NWP_1_file == 1) option
metdyn Jan 3, 2025
9379440
NWP_1_file complete: # of obs points recover, no missing points
metdyn Jan 4, 2025
80185d6
Introduce use_NWP_1_file=1 option to trajectory sampler, in order tha…
metdyn Jan 4, 2025
5308922
Add item in CHANGELOG.md
metdyn Jan 4, 2025
3ad76f7
code cleanup; add clear notes
metdyn Jan 6, 2025
127075c
Add integer parameter: use_NWP_1_file_param
metdyn Jan 6, 2025
8b633dc
use one line for _ASSERT()
metdyn Jan 6, 2025
f51b43a
Merge branch 'develop' into feature/ygyu/NPW_1_file
metdyn Jan 6, 2025
46d77c4
change variable use_NWP_1_file from integer to logical
metdyn Jan 7, 2025
8569d56
Merge branch 'develop' into feature/ygyu/NPW_1_file
metdyn Jan 7, 2025
63ac2cd
Merge branch 'develop' into feature/ygyu/NPW_1_file
metdyn Jan 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 2 additions & 0 deletions base/MAPL_ObsUtil.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 0 additions & 3 deletions gridcomps/History/MAPL_HistoryGridComp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.)
Expand Down
7 changes: 2 additions & 5 deletions gridcomps/History/Sampler/MAPL_TrajectoryMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,11 +48,6 @@ module HistoryTrajectoryMod
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
Expand All @@ -62,6 +57,8 @@ 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, dimension(2), parameter :: use_NWP_1_file_param = [0, 1]
metdyn marked this conversation as resolved.
Show resolved Hide resolved
integer :: epoch ! unit: second
integer(kind=ESMF_KIND_I8) :: epoch_index(2)
real(kind=ESMF_KIND_R8), pointer:: obsTime(:)
Expand Down
102 changes: 78 additions & 24 deletions gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand All @@ -71,6 +72,16 @@
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
_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
end if

call ESMF_ConfigGetAttribute(config, value=STR1, default="", &
label=trim(string) // 'obs_file_begin:', _RC)
Expand Down Expand Up @@ -540,8 +551,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)
Expand All @@ -553,7 +566,7 @@
type(ESMF_VM) :: vm
integer :: mypet, petcount, mpic

integer :: i, j, k, L, ii, jj
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
Expand All @@ -564,6 +577,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
Expand Down Expand Up @@ -612,17 +626,29 @@
trim(this%var_name_lat),trim(this%var_name_time))

L=0
fid_s=this%obsfile_Ts_index
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 ! 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 ! index 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)

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, &
Expand All @@ -632,6 +658,12 @@
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==this%obsfile_Ts_index) then
this%obs(k)%count_location_until_matching_file = jj2
elseif (j==this%obsfile_Ts_index+1) then
this%obs(k)%count_location_in_matching_file = num_times
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))
Expand All @@ -641,6 +673,9 @@
enddo
arr(1)=len

!
! __ 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)
Expand All @@ -652,6 +687,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, &
Expand All @@ -671,7 +707,9 @@
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
! for each obs type use the correct starting point
location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file
end do
len = len + num_times
end if
Expand All @@ -681,7 +719,6 @@
end if
end if


call ESMF_VMAllFullReduce(vm, sendData=arr, recvData=nx_sum, &
count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc)
if (nx_sum == 0) then
Expand All @@ -708,10 +745,15 @@
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(:)
!! use index to sort togehter multiple arrays
allocate(location_index_ioda_full_aux, source=location_index_ioda_full, _STAT)
do jj = 1, nx_sum
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
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)
Expand All @@ -721,7 +763,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:
! 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:
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)
Expand Down Expand Up @@ -753,15 +803,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
Expand Down Expand Up @@ -826,10 +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 %i12', &
'obs(', k, ')%nobs_epoch', this%obs(k)%nobs_epoch )
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)
Expand Down Expand Up @@ -1199,7 +1248,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
Expand Down Expand Up @@ -1232,11 +1281,17 @@
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


!
! __ I designed a method to return from regridding if no valid points exist
Expand Down Expand Up @@ -1427,8 +1482,7 @@
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

x_subset(1) = index1+1
x_subset(2) = index2
Expand Down
Loading