Trajectory sampler: use a single obs file at each Epoch time and rest…

…ore obs location for output (#3326)

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
 

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)
 

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

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))

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)