From 45bdba58d66d7a8f816f52709a543f20744e38fe Mon Sep 17 00:00:00 2001 From: Yonggang Yu Date: Thu, 16 Jan 2025 09:02:01 -0700 Subject: [PATCH] Trajectory sampler: use a single obs file at each Epoch time and restore obs location for output (#3326) --- CHANGELOG.md | 2 + base/MAPL_ObsUtil.F90 | 14 +- base/Plain_netCDF_Time.F90 | 49 ++++ gridcomps/History/MAPL_HistoryGridComp.F90 | 34 +-- .../History/Sampler/MAPL_TrajectoryMod.F90 | 8 +- .../Sampler/MAPL_TrajectoryMod_smod.F90 | 269 ++++++++++++++---- 6 files changed, 299 insertions(+), 77 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d2fdf3ce41b7..a6b309c9ffa0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added utility to prepare inputs for ExtDatDriver.x so that ExtData can simulate a real GEOS run - 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 +- The trajectory sampler netCDF output variable `location_index_in_iodafile` can be turned off, after we add two control variables: `use_NWP_1_file` and `restore_2_obs_vector` for users. When set to true, the two options will select only one obs file at each Epoch interval, and will rotate the output field index back to the location vector inthe obs file before generating netCDF output. +- Support `splitfield: 1` in HISTORY.rc for trajectory sampler ### Changed diff --git a/base/MAPL_ObsUtil.F90 b/base/MAPL_ObsUtil.F90 index 1b96fad92ae8..aed4c9adcd59 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 @@ -38,6 +40,7 @@ module MAPL_ObsUtilMod real(kind=REAL64), allocatable :: lats(:) real(kind=REAL64), allocatable :: times_R8(:) integer, allocatable :: location_index_ioda(:) + integer, allocatable :: restore_index(:) real(kind=REAL32), allocatable :: p2d(:) real(kind=REAL32), allocatable :: p3d(:,:) end type obs_unit @@ -50,7 +53,7 @@ module MAPL_ObsUtilMod character (len=ESMF_MAXSTR) :: var_name_time='' character (len=ESMF_MAXSTR) :: file_name_template='' integer :: ngeoval=0 - integer :: nentry_name=0 + integer :: nfield_name_mx=12 character (len=ESMF_MAXSTR), allocatable :: field_name(:,:) !character (len=ESMF_MAXSTR), allocatable :: field_name(:) end type obs_platform @@ -768,7 +771,7 @@ function copy_platform_nckeys(a, rc) copy_platform_nckeys%var_name_lon = a%var_name_lon copy_platform_nckeys%var_name_lat = a%var_name_lat copy_platform_nckeys%var_name_time = a%var_name_time - copy_platform_nckeys%nentry_name = a%nentry_name + copy_platform_nckeys%nfield_name_mx = a%nfield_name_mx _RETURN(_SUCCESS) end function copy_platform_nckeys @@ -781,7 +784,7 @@ function union_platform(a, b, rc) integer, optional, intent(out) :: rc character (len=ESMF_MAXSTR), allocatable :: field_name_loc(:,:) - integer :: nfield, nentry_name + integer :: nfield, nfield_name_mx integer, allocatable :: tag(:) integer :: i, j, k integer :: status @@ -802,9 +805,9 @@ function union_platform(a, b, rc) enddo union_platform%ngeoval=k nfield=k - nentry_name=union_platform%nentry_name + nfield_name_mx=union_platform%nfield_name_mx if ( allocated (union_platform%field_name) ) deallocate(union_platform%field_name) - allocate(union_platform%field_name(nentry_name, nfield)) + allocate(union_platform%field_name(nfield_name_mx, nfield)) do i=1, a%ngeoval union_platform%field_name(:,i) = a%field_name(:,i) enddo @@ -950,6 +953,7 @@ subroutine fglob(search_name, filename, rc) ! give the last name if (lenmax < slen) then if (MAPL_AM_I_ROOT()) write(6,*) 'pathlen vs filename_max_char_len: ', slen, lenmax _FAIL ('PATHLEN is greater than filename_max_char_len') + STOP 'lenmax < slen' end if if (slen>0) filename(1:slen)=c_filename(1:slen) diff --git a/base/Plain_netCDF_Time.F90 b/base/Plain_netCDF_Time.F90 index b9c163816647..924ba9323978 100644 --- a/base/Plain_netCDF_Time.F90 +++ b/base/Plain_netCDF_Time.F90 @@ -840,6 +840,13 @@ subroutine split_string_by_space (string_in, length_mx, & wc=0 ios=0 string = trim( adjustl(string_in) ) + str_piece(:)='' + i = index (string, mark) + if (i==0) then + nseg=1 + str_piece(1)=string + return + end if do while (ios==0) i = index (string, mark) !!print*, 'index=', i @@ -858,5 +865,47 @@ subroutine split_string_by_space (string_in, length_mx, & end subroutine split_string_by_space + subroutine split_string_by_seperator (string_in, length_mx, seperator_in, & + mxseg, nseg, str_piece, jstatus) + character (len=length_mx), intent (in) :: string_in + integer, intent (in) :: length_mx + character (len=length_mx), intent (in) :: seperator_in + integer, intent (in) :: mxseg + integer, intent (out):: nseg + character (len=length_mx), intent (out):: str_piece(mxseg) + integer, intent (out):: jstatus + character (len=length_mx) :: string_sc, string_oper, string_aux + character (len=1) :: mark, CH + integer :: ios + integer :: wc + integer :: len1, len2 + ! + ! "xxxx; yy; zz; uu, vv," + ! seperator = ";," + ! + + !__ s1. replace seperator by space + ! + string_sc = trim( adjustl(string_in) ) + string_oper = trim( adjustl(seperator_in) ) + len1 = len_trim(string_sc) + len2 = len_trim(string_oper) + string_aux=string_sc + do i = 1, len1 + CH = string_sc(i:i) + do j = 1, len2 + mark = string_oper(j:j) + if (CH==mark) then + string_aux(i:i)=' ' + end if + end do + end do + + !__ s2. split by space + call split_string_by_space (string_aux, length_mx, & + mxseg, nseg, str_piece, jstatus) + + return + end subroutine split_string_by_seperator end module MAPL_scan_pattern_in_file diff --git a/gridcomps/History/MAPL_HistoryGridComp.F90 b/gridcomps/History/MAPL_HistoryGridComp.F90 index c0ee057d4dd9..a6e95338cd1e 100644 --- a/gridcomps/History/MAPL_HistoryGridComp.F90 +++ b/gridcomps/History/MAPL_HistoryGridComp.F90 @@ -5433,8 +5433,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 @@ -5466,7 +5464,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) integer :: nseg integer :: nseg_ub integer :: nfield, nplatform - integer :: nentry_name + integer :: nfield_name_max logical :: obs_flag integer, allocatable :: map(:) type(Logger), pointer :: lgr @@ -5493,7 +5491,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.) @@ -5557,7 +5554,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) - ! __ s2.1 scan fields: only determine ngeoval / nentry_name = nword + ! __ s2.1 scan fields: only determine ngeoval / nfield_name_max = nword allocate (str_piece(mxseg)) rewind(unitr) do k=1, nplf @@ -5581,10 +5578,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end if enddo PLFS(k)%ngeoval = ngeoval - PLFS(k)%nentry_name = nseg_ub + nseg_ub = PLFS(k)%nfield_name_mx allocate ( PLFS(k)%field_name (nseg_ub, ngeoval) ) PLFS(k)%field_name = '' - !! print*, 'k, ngeoval, nentry_name', k, ngeoval, nseg_ub + !! print*, 'k, ngeoval, nfield_name_max', k, ngeoval, nseg_ub end do @@ -5682,12 +5679,12 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) call split_string_by_space (line, length_mx, mxseg, & nplatform, str_piece, status) - !! to do: add debug - !write(6,*) 'line for obsplatforms=', trim(line) - !write(6,*) 'split string, nplatform=', nplatform - !write(6,*) 'nplf=', nplf - !write(6,*) 'str_piece=', str_piece(1:nplatform) - + call lgr%debug('%a %a', 'line for obsplatforms=', trim(line)) + call lgr%debug('%a %i6', 'split string, nplatform=', nplatform) + call lgr%debug('%a %i6', 'nplf=', nplf) + !if (mapl_am_i_root()) then + ! write(6,*) ' str_piece=', str_piece(1:nplatform) + !end if ! ! a) union the platform @@ -5703,7 +5700,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do end do deallocate(str_piece) - !! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !if (mapl_am_i_root()) then + ! write(6,*) 'collection n=',n, 'map(:)=', map(:) + !end if + ! __ write common nc_index,time,lon,lat k=map(1) ! plat form # 1 @@ -5722,10 +5722,10 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) end do nfield = p1%ngeoval - nentry_name = p1%nentry_name + nfield_name_max = p1%nfield_name_mx do j=1, nfield line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(p1%field_name(i,j)) enddo if (j==1) then @@ -5744,7 +5744,7 @@ subroutine regen_rcx_for_obs_platform (config, nlist, list, rc) write(unitw, '(a)') trim(adjustl(PLFS(k)%file_name_template)) do j=1, PLFS(k)%ngeoval line='' - do i=1, nentry_name + do i=1, nfield_name_max line = trim(line)//' '//trim(adjustl(PLFS(k)%field_name(i,j))) enddo write(unitw, '(a)') trim(adjustl(line)) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 index 49c2eff96b38..7cfccd4b6d40 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod.F90 @@ -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 @@ -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 + logical :: use_NWP_1_file = .false. + logical :: restore_2_obs_vector = .false. integer :: epoch ! unit: second integer(kind=ESMF_KIND_I8) :: epoch_index(2) real(kind=ESMF_KIND_R8), pointer:: obsTime(:) @@ -74,6 +71,7 @@ module HistoryTrajectoryMod integer :: obsfile_Te_index logical :: active ! case: when no obs. exist logical :: level_by_level = .false. + ! ! note ! for MPI_GATHERV of 3D data in procedure :: append_file ! we have choice LEVEL_BY_LEVEL or ALL_AT_ONCE (timing in sec below for extdata) diff --git a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 index 4ae3193970c9..d751474c3483 100644 --- a/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 +++ b/gridcomps/History/Sampler/MAPL_TrajectoryMod_smod.F90 @@ -25,33 +25,36 @@ use, intrinsic :: iso_fortran_env, only: REAL64 use, intrinsic :: iso_fortran_env, only: INT64 implicit none - contains module procedure HistoryTrajectory_from_config use BinIOMod + use MAPL_scan_pattern_in_file use pflogger, only : Logger, logging type(ESMF_Time) :: currTime type(ESMF_TimeInterval) :: epoch_frequency type(ESMF_TimeInterval) :: obs_time_span integer :: time_integer, second integer :: status - character(len=ESMF_MAXSTR) :: STR1, line + character(len=ESMF_MAXSTR) :: STR1, line, splitter character(len=ESMF_MAXSTR) :: symd, shms integer :: nline, col integer, allocatable :: ncol(:) character(len=ESMF_MAXSTR), allocatable :: word(:) + character(len=ESMF_MAXSTR), allocatable :: str_piece(:) integer :: nobs, head, jvar logical :: tend integer :: i, j, k, k2, M integer :: count, idx integer :: unitr, unitw + integer :: length_mx, mxseg, nseg type(GriddedIOitem) :: item type(Logger), pointer :: lgr 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 +74,22 @@ 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=.false., & + label=trim(string)//'use_NWP_1_file:', _RC) + call ESMF_ConfigGetAttribute(config, value=traj%restore_2_obs_vector, default=.false., & + label=trim(string)//'restore_2_obs_vector:', _RC) + if (mapl_am_I_root()) then + if (traj%use_NWP_1_file) then + write(6,105) 'WARNING: Traj sampler: use_NWP_1_file is true' + write(6,105) 'WARNING: USER needs to check if observation file is fetched correctly' + end if + if (traj%restore_2_obs_vector) then + write(6,105) 'WARNING: Traj sampler: restore_2_obs_vector is true' + end if + end if + if (.NOT. traj%use_NWP_1_file .AND. traj%restore_2_obs_vector) then + _FAIL('use_NWP_1_file=.false. and restore_2_obs_vector=.true. is not allowed') + end if call ESMF_ConfigGetAttribute(config, value=STR1, default="", & label=trim(string) // 'obs_file_begin:', _RC) @@ -147,6 +166,9 @@ ! __ s3. retrieve template and geoval, set metadata file_handle lgr => logging%get_logger('HISTORY.sampler') + length_mx = ESMF_MAXSTR + mxseg = 100 + allocate (str_piece(mxseg)) if ( nobs == 0 ) then ! ! treatment-1: @@ -192,11 +214,18 @@ ! 1-item case: file template or one-var ! 2-item : var1 , 'root' case STR1=trim(word(1)) + elseif ( count == 3 ) then + ! the Alias case + the splitField case + ! 3-item : var1 , 'root', var1_alias case + ! 3-item : var1 , 'root', 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870', + ! 3-item : 'u;v' vector interpolation is not handled + STR1=trim(word(3)) else - ! 3-item : var1 , 'root', var1_alias case STR1=trim(word(3)) + call lgr%debug('%a %i8', 'WARNING: there are more than 3 field_names in platform rcx' ) end if deallocate(word, _STAT) + if ( index(trim(STR1), '-----') == 0 ) then if (head==1 .AND. trim(STR1)/='') then nobs=nobs+1 @@ -205,13 +234,29 @@ head=0 else if (trim(STR1)/='') then - jvar=jvar+1 - idx = index(STR1,";") - if (idx==0) then - traj%obs(nobs)%geoval_xname(jvar) = STR1 + splitter=';,' + call split_string_by_seperator (STR1, length_mx, splitter, mxseg, & + nseg, str_piece, status) + if (count < 3) then + ! case + ! 'var1' + ! 'var1' , 'ROOT' + ! 'u;v' , 'ROOT' + jvar=jvar+1 + if (nseg==1) then + traj%obs(nobs)%geoval_xname(jvar) = STR1 + else + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(1)) + traj%obs(nobs)%geoval_yname(jvar) = trim(str_piece(2)) + end if else - traj%obs(nobs)%geoval_xname(jvar) = trim(STR1(1:idx-1)) - traj%obs(nobs)%geoval_yname(jvar) = trim(STR1(idx+1:)) + ! case + ! 'var1' , 'ROOT' , alias + ! 'var1' , 'ROOT' , 'TOTEXTTAU470;TOTEXTTAU550;TOTEXTTAU870,' split_field + do j=1, nseg + jvar=jvar+1 + traj%obs(nobs)%geoval_xname(jvar) = trim(str_piece(j)) + end do end if end if end if @@ -223,6 +268,15 @@ enddo end if + !!if (mapl_am_i_root()) then + !! do k=1, nobs + !! do j=1, traj%obs(k)%ngeoval + !! write(6, '(2x,a,2x,2i10,2x,a)') & + !! 'traj%obs(k)%geoval_xname(j), k, j, xname ', k, j, trim(traj%obs(k)%geoval_xname(j)) + !! end do + !! end do + !!end if + do k=1, traj%nobs_type allocate (traj%obs(k)%metadata, _STAT) @@ -335,10 +389,12 @@ call v%add_attribute('long_name', 'dateTime') call this%obs(k)%metadata%add_variable(this%var_name_time,v) - v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) - call v%add_attribute('units', '1') - call v%add_attribute('long_name', 'Location index in corresponding IODA file') - call this%obs(k)%metadata%add_variable(this%location_index_name,v) + if (.NOT. this%restore_2_obs_vector) then + v = Variable(type=PFIO_INT32,dimensions=this%index_name_x) + call v%add_attribute('units', '1') + call v%add_attribute('long_name', 'Location index in corresponding IODA file') + call this%obs(k)%metadata%add_variable(this%location_index_name,v) + end if v = variable(type=PFIO_REAL64,dimensions=this%index_name_x) call v%add_attribute('units','degrees_east') @@ -436,6 +492,7 @@ iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'create new bundle, this%items%xname= ', trim(item%xname) if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%bundle,trim(item%xname),field=src_field,_RC) call ESMF_FieldGet(src_field,rank=rank,_RC) @@ -486,8 +543,7 @@ do k=1, this%nobs_type if (this%obs(k)%nobs_epoch > 0) then filename=trim(this%obs(k)%name)//trim(filename_suffix) - call lgr%debug('%a %a', & - "Sampling to new file : ",trim(filename)) + write(6,'(1x,a,2x,a)') "Sampling to new file :",trim(filename) call this%obs(k)%file_handle%create(trim(filename),_RC) call this%obs(k)%file_handle%write(this%obs(k)%metadata,_RC) end if @@ -540,8 +596,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 +611,7 @@ type(ESMF_VM) :: vm integer :: mypet, petcount, mpic - integer :: i, j, k, L, ii, jj + integer :: i, j, k, L, ii, jj, jj2 integer :: fid_s, fid_e integer(kind=ESMF_KIND_I8) :: j0, j1 integer(kind=ESMF_KIND_I8) :: jt1, jt2 @@ -612,17 +670,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) 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, & @@ -632,6 +702,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)) @@ -641,6 +717,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) @@ -652,6 +731,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 +751,14 @@ 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 + ! if: use_nwp_1_file: index_ioda = [ 1, Nobs ] : restore_2_obs_vector is exact + ! else index_ioda = [ -M, 0 ] + [1, Nobs1] + [Nob1+1, Nobs2] : restore_2_obs_vector may fail + ! File1 File(center) File3 + ! why: bc we have no restriction on observation file name vs content, hence unexpected things can happen + ! use use_nwp_1_file + restore_2_obs_vector only when filename and content are systematic + location_index_ioda_full(len+jj) = jj2 - this%obs(k)%count_location_until_matching_file end do len = len + num_times end if @@ -681,7 +768,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 @@ -708,10 +794,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) @@ -721,7 +812,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) 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) @@ -753,15 +852,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 @@ -803,6 +896,9 @@ allocate (this%obs(k)%lats(nx2), _STAT) allocate (this%obs(k)%times_R8(nx2), _STAT) allocate (this%obs(k)%location_index_ioda(nx2), _STAT) + if (this%use_NWP_1_file) then + allocate (this%obs(k)%restore_index(nx2), _STAT) + end if enddo allocate(ix(this%nobs_type), _STAT) @@ -815,6 +911,11 @@ this%obs(k)%lats(ix(k)) = lats_full(j) this%obs(k)%times_R8(ix(k)) = times_R8_full(j) this%obs(k)%location_index_ioda(ix(k)) = location_index_ioda_full(j) + if (this%use_NWP_1_file) then + ! only this case, we have exact obs in 1_file <-> sampling match + this%obs(k)%restore_index(location_index_ioda_full(j)) = ix(k) + end if + ! !if (mod(k,10**8)==1) then ! write(6,*) 'this%obs(k)%times_R8(ix(k))', this%obs(k)%times_R8(ix(k)) !endif @@ -826,10 +927,14 @@ 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: the time boundary issue can appear when we use python convention [T1, T2) but obs files donot + ! ioda file split [1/2data 15Z : 1/2data 21Z ] [ 1/2data 21Z : 1/2data 3Z] (aircraft) + ! ___x x x x x ___ ---------------------------------- o --o ---o -- o -- + ! debug: negative index (extra) at Tmin missing Tmax + ! debug shows: overcount at Tmin and missing points at Tmax + ! use_NPW_1_file=.true. solves this issue, due to special treatment at time boundaries + ! end if else allocate(this%lons(0),this%lats(0),_STAT) @@ -941,6 +1046,9 @@ real(REAL32), pointer :: p_src(:,:),p_dst(:,:), p_dst_t(:,:) ! _t: transpose real(REAL32), pointer :: p_dst_rt(:,:), p_acc_rt_3d(:,:) real(REAL32), pointer :: pt1(:), pt2(:) + real(REAL64), allocatable :: aux_R8(:) + real(REAL64), allocatable :: aux_R4(:) + integer, allocatable :: vec(:) type(ESMF_Field) :: acc_field_2d_chunk, acc_field_3d_chunk, chunk_field real(REAL32), pointer :: p_acc_chunk_3d(:,:),p_acc_chunk_2d(:) @@ -949,7 +1057,7 @@ integer :: lm integer :: rank integer :: status - integer :: j, k, ig + integer :: j, j2, k, kz, ig integer, allocatable :: ix(:) type(ESMF_VM) :: vm integer :: mypet, petcount, mpic, iroot @@ -977,14 +1085,31 @@ nx=this%obs(k)%nobs_epoch if (nx >0) then if (mapl_am_i_root()) then - call this%obs(k)%file_handle%put_var(this%var_name_time, real(this%obs(k)%times_R8), & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%var_name_lon, this%obs(k)%lons, & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%var_name_lat, this%obs(k)%lats, & - start=[is], count=[nx], _RC) - call this%obs(k)%file_handle%put_var(this%location_index_name, this%obs(k)%location_index_ioda, & - start=[is], count=[nx], _RC) + if (this%restore_2_obs_vector) then + ! restore back to obs vector + allocate (aux_R8(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + aux_R8(1:nx) = this%obs(k)%times_R8(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_time, aux_R8, & + start=[is], count=[nx], _RC) + aux_R8(1:nx) = this%obs(k)%lons(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_lon, aux_R8, & + start=[is], count=[nx], _RC) + aux_R8(1:nx) = this%obs(k)%lats(vec(1:nx)) + call this%obs(k)%file_handle%put_var(this%var_name_lat, aux_R8, & + start=[is], count=[nx], _RC) + deallocate (aux_R8, vec) + else + ! default: location in time sequence + call this%obs(k)%file_handle%put_var(this%var_name_time, this%obs(k)%times_R8, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%var_name_lon, this%obs(k)%lons, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%var_name_lat, this%obs(k)%lats, & + start=[is], count=[nx], _RC) + call this%obs(k)%file_handle%put_var(this%location_index_name, this%obs(k)%location_index_ioda, & + start=[is], count=[nx], _RC) + end if end if end if enddo @@ -1053,9 +1178,12 @@ allocate ( p_dst_rt(lm, 1) ) end if + iter = this%items%begin() do while (iter /= this%items%end()) item => iter%get() + !!if (mapl_am_I_root()) print*, 'regrid: item%xname= ', trim(item%xname) + if (item%itemType == ItemTypeScalar) then call ESMF_FieldBundleGet(this%acc_bundle,trim(item%xname),field=acc_field,_RC) call ESMF_FieldGet(acc_field,rank=rank,_RC) @@ -1095,16 +1223,33 @@ endif enddo deallocate(ix, _STAT) + + ! rotate 2d field + if (this%restore_2_obs_vector) then + do k=1, this%nobs_type + nx = this%obs(k)%nobs_epoch + if (nx>0) then + allocate (aux_R4(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + aux_R4(1:nx) = this%obs(k)%p2d(vec(1:nx)) + this%obs(k)%p2d(1:nx) = aux_R4(1:nx) + deallocate (aux_R4, vec) + end if + end do + end if + do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch if (nx>0) then do ig = 1, this%obs(k)%ngeoval + !! print*, 'this%obs(k)%geoval_xname(ig)= ', this%obs(k)%geoval_xname(ig) + if (trim(item%xname) == trim(this%obs(k)%geoval_xname(ig))) then call this%obs(k)%file_handle%put_var(trim(item%xname), this%obs(k)%p2d(1:nx), & start=[is],count=[nx]) end if - enddo + end do endif enddo do k=1, this%nobs_type @@ -1145,7 +1290,6 @@ call ESMF_FieldDestroy(dst_field,noGarbage=.true.,_RC) call ESMF_FieldDestroy(src_field,noGarbage=.true.,_RC) - if (mapl_am_i_root()) then ! !-- pack fields to obs(k)%p3d and put_var @@ -1164,6 +1308,23 @@ this%obs(k)%p3d(ix(k),:) = p_acc_rt_3d(j,:) enddo deallocate(ix, _STAT) + + ! rotate 3d field + if (this%restore_2_obs_vector) then + do k=1, this%nobs_type + nx = this%obs(k)%nobs_epoch + if (nx>0) then + allocate (aux_R4(nx), vec(nx)) + vec(1:nx) = this%obs(k)%restore_index(1:nx) + do kz=1, lm + aux_R4(1:nx) = this%obs(k)%p3d(vec(1:nx), kz) + this%obs(k)%p3d(1:nx, kz) = aux_R4(1:nx) + end do + deallocate (aux_R4, vec) + end if + end do + end if + do k=1, this%nobs_type is = 1 nx = this%obs(k)%nobs_epoch @@ -1199,7 +1360,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 @@ -1232,11 +1393,17 @@ call ESMF_ClockGet(this%clock,timeStep=dur, _RC ) timeset(1) = current_time - dur timeset(2) = current_time + if (this%use_NWP_1_file) 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 @@ -1351,6 +1518,9 @@ if (allocated (this%obs(k)%location_index_ioda)) then deallocate (this%obs(k)%location_index_ioda) end if + if (allocated (this%obs(k)%restore_index)) then + deallocate (this%obs(k)%restore_index) + end if if (allocated(this%obs(k)%p2d)) then deallocate (this%obs(k)%p2d) endif @@ -1427,8 +1597,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