Trajectory sampler for IODA file: step-1

CHANGELOG.md

@@ -6,7 +6,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
-
+- sampling IODA file with trajectory sampler (step-1): make it run
 ### Added
 
 - Convert ExtData to use ESMF HConfig for YAML parsing rather than YaFYAML

gridcomps/History/MAPL_HistoryCollection.F90

@@ -101,7 +101,7 @@ module MAPL_HistoryCollectionMod
      character(len=ESMF_MAXSTR)         :: output_grid_label
      type(GriddedIOItemVector)          :: items
      character(len=ESMF_MAXSTR)         :: currentFile
-     character(len=ESMF_MAXPATHLEN)     :: trackFile
+     character(len=ESMF_MAXPATHLEN)     :: obsFile
      character(len=ESMF_MAXPATHLEN)     :: stationIdFile
      logical                            :: splitField
      logical                            :: regex

gridcomps/History/MAPL_HistoryGridComp.F90

@@ -873,9 +873,9 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
             label=trim(string) // 'station_id_file:', _RC)
 
 ! Get an optional file containing a 1-D track for the output
-       call ESMF_ConfigGetAttribute(cfg, value=list(n)%trackFile, default="", &
+       call ESMF_ConfigGetAttribute(cfg, value=list(n)%obsFile, default="", &
                                     label=trim(string) // 'track_file:', _RC)
-       if (trim(list(n)%trackfile) /= '') list(n)%timeseries_output = .true.
+       if (trim(list(n)%obsFile) /= '') list(n)%timeseries_output = .true.
        call ESMF_ConfigGetAttribute(cfg, value=list(n)%recycle_track, default=.false., &
                                     label=trim(string) // 'recycle_track:', _RC)
 
@@ -2338,6 +2338,9 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
 
     do n=1,nlist
        if (list(n)%disabled) cycle
+       string = trim( list(n)%collection ) // '.'
+       cfg = ESMF_ConfigCreate(_RC)
+       call ESMF_ConfigLoadFile(cfg, filename = trim(string)//'rcx', _RC)
        if (list(n)%format == 'CFIOasync') then
           list(n)%format = 'CFIO'
           if (mapl_am_i_root()) write(*,*)'Chose CFIOasync setting to CFIO, update your History.rc file'
@@ -2367,8 +2370,9 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
              list(n)%timeInfo = TimeData(clock,tm,MAPL_nsecf(list(n)%frequency),IntState%stampoffset(n),integer_time=intstate%integer_time)
           end if
           if (list(n)%timeseries_output) then
-             list(n)%trajectory = HistoryTrajectory(trim(list(n)%trackfile),_RC)
+             list(n)%trajectory = HistoryTrajectory(cfg,string,_RC)
              call list(n)%trajectory%initialize(list(n)%items,list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,recycle_track=list(n)%recycle_track,_RC)
+
           elseif (list(n)%sampler_spec == 'station') then
              list(n)%station_sampler = StationSampler (trim(list(n)%stationIdFile),_RC)
              call list(n)%station_sampler%add_metadata_route_handle(list(n)%bundle,list(n)%timeInfo,vdata=list(n)%vdata,_RC)
@@ -2384,6 +2388,7 @@ subroutine Initialize ( gc, import, dumexport, clock, rc )
              call list(n)%mGriddedIO%set_param(write_collection_id=collection_id)
           end if
        end if
+       call ESMF_ConfigDestroy(cfg, _RC)
    end do
 
 ! Echo History List Data Structure

gridcomps/History/MAPL_HistoryTrajectoryMod.F90

@@ -15,6 +15,8 @@ module HistoryTrajectoryMod
    use MAPL_VerticalDataMod
    use MAPL_BaseMod
    use MAPL_CommsMod
+   use MAPL_SortMod
+   use MAPL_NetCDF
    use MAPL_LocstreamRegridderMod
    use, intrinsic :: iso_fortran_env, only: REAL32
    use, intrinsic :: iso_fortran_env, only: REAL64
@@ -28,6 +30,7 @@ module HistoryTrajectoryMod
       type(ESMF_LocStream) :: root_locstream,dist_locstream
       type(LocStreamFactory) :: locstream_factory
       type(ESMF_Time), allocatable :: times(:)
+      real(kind=REAL64), allocatable :: times_R8(:)
       real(kind=REAL64), allocatable :: lons(:),lats(:)
       type(ESMF_FieldBundle) :: bundle
       type(ESMF_FieldBundle) :: output_bundle
@@ -43,6 +46,15 @@ module HistoryTrajectoryMod
       character(LEN=ESMF_MAXPATHLEN) :: file_name
       type(TimeData) :: time_info
       logical :: recycle_track
+      character(len=ESMF_MAXSTR)     :: obsFile
+      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)     :: var_name_time
+      character(len=ESMF_MAXSTR)     :: var_name_lat
+      character(len=ESMF_MAXSTR)     :: var_name_lon
+      character(len=ESMF_MAXSTR)     :: datetime_units
       contains
          procedure :: initialize
          procedure :: create_variable
@@ -55,47 +67,102 @@ module HistoryTrajectoryMod
          procedure :: get_file_start_time
          procedure :: get
          procedure :: reset_times_to_current_day
+         procedure :: sort_arrays_by_time
+         procedure :: time_real_to_ESMF
+
    end type
 
    interface HistoryTrajectory
-      module procedure HistoryTrajectory_from_file
+      module procedure HistoryTrajectory_from_config
    end interface HistoryTrajectory
 
    contains
 
-      function HistoryTrajectory_from_file(filename,unusable,rc) result(trajectory)
-         type(HistoryTrajectory) :: trajectory
-         character(len=*), intent(in) :: filename
+      function HistoryTrajectory_from_config(config,string,unusable,rc) result(traj)
+         use pflogger, only : Logger, logging
+         type(HistoryTrajectory) :: traj
+         type(ESMF_Config), intent(inout) :: config
+         character(len=*), intent(in)     :: string
          class (KeywordEnforcer), optional, intent(in) :: unusable
          integer, optional, intent(out) :: rc
          integer :: status
 
+         character(len=ESMF_MAXSTR) :: filename
          type(NetCDF4_FileFormatter) :: formatter
-         type(FileMetadataUtils) :: metadata
+         type(FileMetadataUtils) :: metadata_utils
          type(FileMetadata) :: basic_metadata
-         integer :: num_times
+         integer(ESMF_KIND_I8) :: num_times
+
+         integer :: ncid, grpid, ncid0
+         integer :: dimid(10),  dimlen(10)
+         integer :: len
+         integer :: i
+         character(len=ESMF_MAXSTR) :: grp_name
+         character(len=ESMF_MAXSTR) :: dim_name(10)
+         character(len=ESMF_MAXSTR) :: var_name_lon
+         character(len=ESMF_MAXSTR) :: var_name_lat
+         character(len=ESMF_MAXSTR) :: var_name_time
+         type(Logger), pointer :: lgr
 
          _UNUSED_DUMMY(unusable)
 
+         call ESMF_ConfigGetAttribute(config, value=traj%obsFile, default="", &
+              label=trim(string) // 'track_file:', _RC)
+         call ESMF_ConfigGetAttribute(config, value=traj%nc_index, default="", &
+              label=trim(string) // 'nc_Index:', _RC)
+         call ESMF_ConfigGetAttribute(config, value=traj%nc_time, default="", &
+              label=trim(string) // 'nc_Time:', _RC)
+         call ESMF_ConfigGetAttribute(config, value=traj%nc_longitude, default="", &
+              label=trim(string) // 'nc_Longitude:', _RC)
+         call ESMF_ConfigGetAttribute(config, value=traj%nc_latitude, default="", &
+              label=trim(string) // 'nc_Latitude:', _RC)
+
+         traj%datetime_units = "seconds since 1970-01-01 00:00:00"
+
+         filename=trim(traj%obsFile)
          call formatter%open(trim(filename),pFIO_READ,_RC)
-         basic_metadata = formatter%read(_RC)
-         call metadata%create(basic_metadata,trim(filename))
-         num_times = metadata%get_dimension("time",_RC)
-         allocate(trajectory%lons(num_times),trajectory%lats(num_times),_STAT)
-         if (metadata%is_var_present("longitude")) then
-            call formatter%get_var("longitude",trajectory%lons,_RC)
-         end if
-         if (metadata%is_var_present("latitude")) then
-            call formatter%get_var("latitude",trajectory%lats,_RC)
-         end if
+         if (traj%nc_index == '') then
+            basic_metadata = formatter%read(_RC)
+            call metadata_utils%create(basic_metadata,trim(filename))
+            num_times = metadata_utils%get_dimension("time",_RC)
+            allocate(traj%lons(num_times),traj%lats(num_times),_STAT)
+            if (metadata_utils%is_var_present("longitude")) then
+               call formatter%get_var("longitude",traj%lons,_RC)
+            end if
+            if (metadata_utils%is_var_present("latitude")) then
+               call formatter%get_var("latitude",traj%lats,_RC)
+            end if
+            call metadata_utils%get_time_info(timeVector=traj%times,_RC)
+         else
+            i=index(traj%nc_longitude, '/')
+            _ASSERT (i>0, 'group name not found')
+            grp_name = traj%nc_longitude(1:i-1)
+            traj%var_name_lat = traj%nc_latitude(i+1:)
+            traj%var_name_lon = traj%nc_longitude(i+1:)
+            traj%var_name_time= traj%nc_time(i+1:)
+
+            call formatter%open(trim(filename),pFIO_READ,_RC)
+            basic_metadata = formatter%read(_RC)
+            call metadata_utils%create(basic_metadata,trim(filename))
+            num_times = basic_metadata%get_dimension(trim(traj%nc_index),_RC)
+            len = num_times
+
+            allocate(traj%lons(len),traj%lats(len),_STAT)
+            allocate(traj%times_R8(len),traj%times(len),_STAT)
+            call formatter%get_var(traj%var_name_lon,  traj%lons, group_name=grp_name, count=[len], rc=status)
+            call formatter%get_var(traj%var_name_lat,  traj%lats, group_name=grp_name, count=[len], rc=status)
+            call formatter%get_var(traj%var_name_time, traj%times_R8, group_name=grp_name, count=[len], rc=status)
+
+            call traj%sort_arrays_by_time(_RC)
+            call traj%time_real_to_ESMF(_RC)
+         endif
 
-         call metadata%get_time_info(timeVector=trajectory%times,_RC)
-         trajectory%locstream_factory = LocStreamFactory(trajectory%lons,trajectory%lats,_RC)
-         trajectory%root_locstream = trajectory%locstream_factory%create_locstream(_RC)
+         traj%locstream_factory = LocStreamFactory(traj%lons,traj%lats,_RC)
+         traj%root_locstream = traj%locstream_factory%create_locstream(_RC)
 
          _RETURN(_SUCCESS)
 
-      end function HistoryTrajectory_from_file
+      end function HistoryTrajectory_from_config
 
       subroutine initialize(this,items,bundle,timeInfo,unusable,vdata,recycle_track,rc)
          class(HistoryTrajectory), intent(inout) :: this
@@ -563,4 +630,65 @@ subroutine reset_times_to_current_day(this,rc)
 
       end subroutine reset_times_to_current_day
 
+
+      subroutine sort_arrays_by_time(this,rc)
+         class(HistoryTrajectory), intent(inout) :: this
+         integer, optional, intent(out) :: rc
+         integer :: status
+
+         integer :: i, len
+         integer, allocatable :: IA(:)
+         real(ESMF_KIND_R8), allocatable :: X(:), Y(:)
+         integer(ESMF_KIND_I8), allocatable :: IX(:)
+
+         len = size (this%times_R8)
+         allocate (IA(len), IX(len), X(len))
+         do i=1, len
+            IX(i)=this%times_R8(i)
+            IA(i)=i
+         enddo
+         call MAPL_Sort(IX,IA)
+
+         X = this%lons
+         do i=1, len
+            this%lons(i) = X(IA(i))
+         enddo
+         X = this%lats
+         do i=1, len
+            this%lats(i) = X(IA(i))
+         enddo
+         X = this%times_R8
+         do i=1, len
+            this%times_R8(i) = X(IA(i))
+         enddo
+
+         _RETURN(_SUCCESS)
+       end subroutine sort_arrays_by_time
+
+
+       subroutine time_real_to_ESMF (this,rc)
+         class(HistoryTrajectory), intent(inout) :: this
+         integer, optional, intent(out) :: rc
+         integer :: status
+
+         integer :: i, len
+         integer :: int_time
+         type(ESMF_TimeInterval) :: interval
+         type(ESMF_Time) :: time0
+         type(ESMF_Time) :: time1
+         character(len=:), allocatable :: tunit
+         character(len=ESMF_MAXSTR) :: datetime_units
+
+         datetime_units = this%datetime_units
+         len = size (this%times_R8)
+
+         do i=1, len
+            int_time = this%times_R8(i)
+            call convert_NetCDF_DateTime_to_ESMF(int_time, datetime_units, interval, time0, time1=time1, tunit=tunit, _RC)
+            this%times(i) = time1
+         enddo
+
+         _RETURN(_SUCCESS)
+       end subroutine time_real_to_ESMF
+
 end module HistoryTrajectoryMod