Feature/bmauer/identify vertical coord extdata2g

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added experimental capabiltiy to regrid from constituents in mass mixing or emissions units from one set of hybrid sigma levels to model levels in ExtData2G
 - Added subdirectory GetHorzIJIndex in the Tests directory for testing subroutine MAPL_GetHorzIJIndex to ensure the subroutine fails is teh cubed-sphere assumption is violated and also if the returned I and J indices are not correct.
 
 ### Changed

CMakeLists.txt

@@ -230,6 +230,7 @@ add_subdirectory (base)
 add_subdirectory (MAPL)
 add_subdirectory (gridcomps)
 add_subdirectory (griddedio)
+add_subdirectory (vertical)
 if (BUILD_WITH_FARGPARSE)
    add_subdirectory (docs)
    add_subdirectory (benchmarks)

Tests/ExtDataDriverGridComp.F90

@@ -145,6 +145,7 @@ subroutine initialize_gc(gc, import_state, export_state, clock, rc)
     class(BaseProfiler), pointer :: t_p
     logical :: use_extdata2g
     integer :: useShmem
+    type(ESMF_StateItem_Flag) :: item_type
 
     _UNUSED_DUMMY(import_state)
     _UNUSED_DUMMY(export_state)
@@ -403,6 +404,12 @@ subroutine initialize_gc(gc, import_state, export_state, clock, rc)
              _VERIFY(status)
           end if
        END DO
+       ! now add PS from Root if available
+       call ESMF_StateGet(cap%exports(cap%root_id), 'PLE', item_type, _RC)
+       if (item_type == ESMF_STATEITEM_FIELD) then
+          call ESMF_StateGet(cap%exports(cap%root_id), 'PLE', field, _RC)
+          call MAPL_StateAdd(cap%imports(cap%extdata_id), field, _RC)
+       end if
        deallocate(itemtypes)
        deallocate(itemnames)
 

Tests/ExtDataRoot_GridComp.F90

@@ -132,6 +132,18 @@ subroutine SetServices ( GC, RC )
                units = 'na', &
                dims = vloc, &
                vlocation = MAPL_VLocationNone, _RC)
+         call MAPL_AddInternalSpec(GC,&
+               short_name='numlev', &
+               long_name='level number' , &
+               units = 'na', &
+               dims = MAPL_DimsHorzVert, &
+               vlocation = MAPL_VLocationCenter, _RC)
+         call MAPL_AddInternalSpec(GC,&
+               short_name='PLE', &
+               long_name='PLE' , &
+               units = 'Pa', &
+               dims = MAPL_DimsHorzVert, &
+               vlocation = MAPL_VLocationEdge, _RC)
          call MAPL_AddExportSpec(GC, &
                short_name='test_bundle', &
                long_name='test', &
@@ -167,6 +179,7 @@ SUBROUTINE Initialize_ ( GC, IMPORT, EXPORT, CLOCK, rc )
          character(len=ESMF_MaxStr) :: key, keyVal
          type(MAPL_MetaComp), pointer :: MAPL
          logical :: isPresent, fill_bundle
+         type(ESMF_State) :: internal
 
          call ESMF_GridCompGet( GC, name=comp_name, config=CF, _RC )
          call MAPL_GetObjectFromGC ( GC, MAPL, _RC )
@@ -211,9 +224,13 @@ SUBROUTINE Initialize_ ( GC, IMPORT, EXPORT, CLOCK, rc )
             call FillBundle(Export,_RC)
          end if
 
+         call MAPL_Get ( MAPL, internal_esmf_state=internal, cf=cf, _RC )
+         call FillState(internal,export,currTime,grid,synth,cf,_RC)
+
          _RETURN(ESMF_SUCCESS)
       contains
-
+
+
             subroutine set_locstream(rc)
 
             integer, optional, intent(out) :: rc
@@ -299,14 +316,14 @@ SUBROUTINE Run_ ( GC, IMPORT, EXPORT, CLOCK, rc )
 
          case(RunModeGenerateExports)
 
-            call FillState(internal,export,currTime,grid,synth,_RC)
+            call FillState(internal,export,currTime,grid,synth,cf,_RC)
 
          case(RunModeGenerateImports)
 
-            call FillState(internal,import,currTime,grid,synth,_RC)
+            call FillState(internal,import,currTime,grid,synth,cf,_RC)
 
          case(runModecompareImports)
-            call FillState(internal,export,currTime,grid,synth,_RC)
+            call FillState(internal,export,currTime,grid,synth,cf,_RC)
             call CompareState(import,export,0.001,_RC)
 
          case(runModeFillImport)
@@ -539,21 +556,22 @@ subroutine CopyState(inState,outState,rc)
 
    end subroutine CopyState
 
-   subroutine FillState(inState,outState,time,grid,Synth,rc)
+   subroutine FillState(inState,outState,time,grid,Synth,cf,rc)
 
       type(ESMF_State), intent(inout) :: inState
       type(ESMF_State), intent(inout) :: outState
       type(ESMF_Time),  intent(Inout) :: time
       type(ESMF_Grid),  intent(inout) :: grid
       type(SyntheticFieldSupport) :: synth
+      type(ESMF_Config), intent(inout) :: cf
       integer, optional, intent(out) :: rc
 
       integer :: status
-      real, pointer                       :: Exptr2(:,:), Exptr1(:)
+      real, pointer                       :: Exptr3(:,:,:), Exptr2(:,:), Exptr1(:)
       integer :: itemcount
       character(len=ESMF_MAXSTR), allocatable :: outNameList(:)
       type(ESMF_StateItem_Flag), allocatable :: item_type(:)
-      type(ESMF_Field) :: expf,farray(7)
+      type(ESMF_Field) :: expf,farray(9)
       type(ESMF_State) :: pstate
       character(len=:), pointer :: fexpr
       integer :: i1,in,j1,jn,ldims(3),i,j,seed_size,mypet
@@ -592,6 +610,12 @@ subroutine FillState(inState,outState,time,grid,Synth,rc)
                exPtr2(i,j)=j1+j-1
             enddo
          enddo
+         call MAPL_GetPointer(inState,exPtr3,'numlev',_RC)
+         do i=1,size(exPtr3,3)
+            exPtr3(:,:,i)=i
+         enddo
+         call MAPL_GetPointer(inState,exPtr3,'PLE',_RC)
+         call Fill_PLE(exPtr3, cf, _RC)
       end if
 
       if (synth%on_tiles) then
@@ -623,19 +647,71 @@ subroutine FillState(inState,outState,time,grid,Synth,rc)
       call ESMF_StateGet(inState,'j_index',farray(5),_RC)
       call ESMF_StateGet(inState,'doy',farray(6),_RC)
       call ESMF_StateGet(inState,'rand',farray(7),_RC)
+      call ESMF_StateGet(inState,'numlev',farray(8),_RC)
+      call ESMF_StateGet(inState,'PLE',farray(9),_RC)
       pstate = ESMF_StateCreate(_RC)
       call ESMF_StateAdd(pstate,farray,_RC)
 
       do i=1,itemCount
          if (item_type(i) == ESMF_STATEITEM_FIELD) then
             call ESMF_StateGet(outState,trim(outNameList(i)),expf,_RC)
             fexpr => synth%fillDefs%at(trim(outNameList(i)))
-            call MAPL_StateEval(pstate,fexpr,expf,_RC)
+            if (associated(fexpr)) then
+               call MAPL_StateEval(pstate,fexpr,expf,_RC)
+            end if
          end if
       enddo
 
       _RETURN(ESMF_SUCCESS)
 
+      contains 
+
+      subroutine fill_ple(ple_ptr, cf, rc)
+         real, pointer, intent(in) :: ple_ptr(:,:,:)
+         type(ESMF_Config), intent(inout) :: cf
+         integer, intent(out), optional :: rc
+
+         integer :: status
+
+         type(ESMF_HConfig) :: hconfig, akbk
+         character(len=3) :: km_str
+         integer :: km, i
+         character(len=ESMF_MAXPATHLEN) :: akbk_file, ps_file
+         real, allocatable :: ak(:), bk(:), ps(:,:)
+         real :: ps_val
+         logical :: has_akbk, has_psval
+         type(NetCDF4_FileFormatter) :: fs
+
+         call ESMF_ConfigFindLabel(cf,"akbk_file:",isPresent=has_akbk,_RC)
+         !call ESMF_ConfigFindLabel(cf,"ps_val:",isPresent=has_psval,_RC)
+         call ESMF_ConfigFindLabel(cf,"ps_file:",isPresent=has_psval,_RC)
+         if (has_akbk .and. has_psval) then
+            call ESMF_ConfigGetAttribute(cf,label="akbk_file:",value=akbk_file,_RC)
+            !call ESMF_ConfigGetAttribute(cf,label="ps_val:",value=ps_val,_RC)
+            call ESMF_ConfigGetAttribute(cf,label="ps_file:",value=ps_file,_RC)
+         else
+            _RETURN(_SUCCESS)
+         end if
+
+         allocate(ps(size(ple_ptr,1),size(ple_ptr,2))) 
+         call fs%open(trim(ps_file),pFIO_READ,_RC)
+         call fs%get_var("PS",ps, _RC)
+
+         hconfig = ESMF_HConfigCreate(filename=trim(akbk_file), _RC)
+         km = size(ple_ptr,3) - 1
+         write(km_str,'(i3.3)') km
+         akbk = ESMF_HConfigCreateAt(hconfig, keyString='L'//trim(km_str), _RC) 
+         ak = ESMF_HConfigAsR4Seq(akbk, keyString='ak', _RC)
+         bk = ESMF_HConfigAsR4Seq(akbk, keyString='bk', _RC)
+
+         do i=1,km+1
+            !ple_ptr(:,:,i-1) = ak(i)+ps_val*bk(i)
+            ple_ptr(:,:,i-1) = ak(i)+ps(:,:)*bk(i)
+         enddo
+
+         _RETURN(_SUCCESS)
+      end subroutine fill_ple
+
    end subroutine FillState
 
    subroutine FillBundle(inState,rc)
@@ -692,6 +768,7 @@ subroutine CompareState(State1,State2,tol,rc)
          _VERIFY(status)
          call ESMF_StateGet(State1,itemNameList=NameList,_RC)
          do i=1,itemCount
+            if (trim(nameList(i)) == 'PLE') cycle
             call ESMF_StateGet(State1,trim(nameList(i)),field1,_RC)
             call ESMF_StateGet(State2,trim(nameList(i)),field2,_RC)
             call ESMF_FieldGet(field1,rank=rank1,_RC)

Tests/ExtData_Testing_Framework/test_cases/case27/CAP.rc

@@ -1,4 +1,6 @@
 CASES::
+CAP1.rc
+CAP2.rc
 CAP3.rc
 CAP4.rc
 ::

Tests/ExtData_Testing_Framework/test_cases/case28/extdata.yaml

@@ -1,4 +1,4 @@
 Collections:
-   fstream1: {template: case1.%y4%m2%d2.nc4, valid_range: "2004-01-01/2005-01-01" }
+   fstream1: {template: case1.%y4%m2%d2.nc4}
 Exports:
    VAR2D: {variable: VAR2D, collection: fstream1}

Tests/ExtData_Testing_Framework/test_cases/case29/extdata.yaml

@@ -1,5 +1,5 @@
 Collections:
-   fstream1: {template: case1.%y4%m2%d2.nc4, valid_range: "2004-01-01/2005-01-01" }
+   fstream1: {template: case1.%y4%m2%d2.nc4}
 Samplings:
    s1: {update_frequency: PT24H, update_offset: PT12H, update_reference_time: '0'}
 Exports:

Tests/ExtData_Testing_Framework/test_cases/case37/CAP1.rc

@@ -2,7 +2,7 @@ ROOT_NAME: Root
 ROOT_CF: AGCM1.rc
 HIST_CF: HISTORY1.rc
 
-BEG_DATE:     20070926 230000
+BEG_DATE:     20070921 230000
 
-JOB_SGMT:     00000004 000000
+JOB_SGMT:     00000009 000000
 HEARTBEAT_DT:   3600

Tests/ExtData_Testing_Framework/test_cases/case37/extdata.yaml

@@ -1,12 +1,9 @@
 Collections:
-   fstream1: {template: "case1.%y4%m2%d2.nc4", valid_range: "2007-09-22/2007-09-30"}
-   fstream2: {template: "case2.%y4%m2%d2.nc4", valid_range: "2007-10-01/2007-10-08"}
+   fstream1: {template: "case1.%y4%m2%d2.nc4", valid_range: "2007-09-22/2007-10-08"}
 
 Samplings:
   s1:
     time_interpolation: False
 
 Exports:
-   VAR1:
-     - {starting: "2007-09-22", variable: VAR1, collection: fstream1, sample: s1, fail_on_missing_file: false}
-     - {starting: "2007-10-01", variable: VAR1, collection: fstream2, sample: s1, fail_on_missing_file: false}
+   VAR1:  {variable: VAR1, collection: fstream1, sample: s1, fail_on_missing_file: false}

Tests/ExtData_Testing_Framework/test_cases/test_case_descriptions.md

@@ -42,3 +42,4 @@ path_to_script/run_extdatadriver_cases.py --builddir path_to_geos_install/bin --
 35. Case1 with zstandard compression (only enabled if netcdf built with zstandard support)
 36. Case1 with zstandard compression and NetCDF granular_bitround quantization (only enabled if netcdf built with quantization
     support and zstandard support)
+37. Test that when time interpolation is off, missing files is allowed, that the last read value is persisted during the period of the missing file

base/ApplicationSupport.F90

@@ -5,6 +5,7 @@ module MAPL_ApplicationSupport
  use MAPL_KeywordEnforcerMod
  use pflogger, only: logging
  use pflogger, only: Logger
+ use udunits2f, initialize_udunits => initialize, finalize_udunits => finalize
  use MAPL_Profiler, initialize_profiler =>initialize, finalize_profiler =>finalize
 
  implicit none
@@ -43,6 +44,7 @@ subroutine MAPL_Initialize(unusable,comm,logging_config,rc)
       call initialize_profiler(comm=comm_world)
       call start_global_time_profiler(rc=status)
       _VERIFY(status)
+      call initialize_udunits(_RC)
       _RETURN(_SUCCESS)
 
    end subroutine MAPL_Initialize
@@ -56,6 +58,7 @@ subroutine MAPL_Finalize(unusable,comm,rc)
 
       _UNUSED_DUMMY(unusable)
 
+      call finalize_udunits()
       if (present(comm)) then
          comm_world = comm
       else

base/Base/Base_Base_implementation.F90

@@ -1270,6 +1270,7 @@ module function MAPL_FieldCreateNewgrid(FIELD, GRID, LM, NEWNAME, RC) RESULT(F)
     integer                 :: lbnds(ESMF_MAXDIM), ubnds(ESMF_MAXDIM)
     character(len=ESMF_MAXSTR) :: newName_
     character(len=ESMF_MAXSTR), parameter :: Iam='MAPL_FieldCreateNewgrid'
+    real, pointer :: ptr1d(:)
 
     call ESMF_FieldGet(FIELD, grid=fgrid, _RC)
 
@@ -1355,6 +1356,9 @@ module function MAPL_FieldCreateNewgrid(FIELD, GRID, LM, NEWNAME, RC) RESULT(F)
     ! otherwise we will overwrite it
     call ESMF_AttributeSet(F, NAME='DIMS', VALUE=DIMS, _RC)
 
+    call assign_fptr(f, ptr1d, _RC)
+    ptr1d = 0.0    
+
     _RETURN(ESMF_SUCCESS)
   end function MAPL_FieldCreateNewgrid
 

base/CMakeLists.txt

@@ -58,7 +58,7 @@ endif()
 
 esma_add_library(
   ${this} SRCS ${srcs}
-  DEPENDENCIES MAPL.shared MAPL.constants MAPL.profiler MAPL.pfio MAPL_cfio_r4 MAPL.field_utils PFLOGGER::pflogger
+  DEPENDENCIES MAPL.shared MAPL.constants MAPL.profiler MAPL.pfio MAPL_cfio_r4 MAPL.field_utils udunits2f PFLOGGER::pflogger
                GFTL_SHARED::gftl-shared-v2 GFTL_SHARED::gftl-shared-v1  GFTL::gftl-v2 GFTL::gftl-v1
                ESMF::ESMF NetCDF::NetCDF_Fortran MPI::MPI_Fortran
   TYPE ${MAPL_LIBRARY_TYPE})

gridcomps/ExtData2G/CMakeLists.txt

@@ -25,7 +25,7 @@ set (srcs
     )
 
 
-esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL.shared MAPL.base MAPL.generic MAPL.griddedio TYPE ${MAPL_LIBRARY_TYPE})
+esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL.shared MAPL.base MAPL.generic MAPL.griddedio udunits2f MAPL.vertical TYPE ${MAPL_LIBRARY_TYPE})
 target_link_libraries (${this} PUBLIC GFTL::gftl GFTL_SHARED::gftl-shared ESMF::ESMF NetCDF::NetCDF_Fortran
                                PRIVATE MPI::MPI_Fortran)
 target_include_directories (${this} PUBLIC $<BUILD_INTERFACE:${MAPL_SOURCE_DIR}/include>)

gridcomps/ExtData2G/ExtDataAbstractFileHandler.F90

@@ -19,6 +19,7 @@ module MAPL_ExtdataAbstractFileHandler
    private
    public :: ExtDataAbstractFileHandler
 
+   integer, parameter :: MAX_TRIALS = 10
    type, abstract :: ExtDataAbstractFileHandler
       character(:), allocatable :: file_template
       type(ESMF_TimeInterval) :: frequency
@@ -30,6 +31,7 @@ module MAPL_ExtdataAbstractFileHandler
          procedure :: initialize
          procedure :: make_metadata
          procedure :: get_time_on_file
+         procedure :: find_any_file
          procedure(get_file_bracket), deferred :: get_file_bracket
    end type
 
@@ -166,4 +168,39 @@ subroutine make_metadata(this,file,metadata,rc)
 
   end subroutine make_metadata
 
+  function find_any_file(this, current_time, rc) result(filename)
+     character(len=:), allocatable :: filename
+     class(ExtDataAbstractFileHandler), intent(inout) :: this
+     type(ESMF_Time), intent(in) :: current_time
+     integer, optional, intent(out) :: rc
+
+     integer :: status, i
+     type(ESMF_Time) :: useable_time
+     character(len=ESMF_MAXPATHLEN) :: trial_file
+     logical :: file_found
+
+     useable_time = current_time
+     if (allocated(this%valid_range)) then
+        useable_time = this%valid_range(1)
+     end if
+     call fill_grads_template(trial_file, this%file_template, time=useable_time, _RC)
+     inquire(file=trim(trial_file),exist=file_found)
+
+     if (file_found) then
+        filename = trial_file
+        _RETURN(_SUCCESS)
+     end if
+     do i=1, MAX_TRIALS
+        useable_time = useable_time + this%frequency
+        call fill_grads_template(trial_file, this%file_template, time=useable_time, _RC)
+        inquire(file=trim(trial_file),exist=file_found)
+        if (file_found) exit
+     enddo
+
+     _ASSERT(file_found,"Could not find any file to open to determine metadata after multiple trials")
+     filename = trial_file
+     _RETURN(_SUCCESS) 
+
+   end function find_any_file
+
 end module MAPL_ExtdataAbstractFileHandler

gridcomps/ExtData2G/ExtDataBracket.F90

@@ -184,7 +184,6 @@ subroutine interpolate_to_time(this,field,time,rc)
       real, pointer              :: var1d_right(:)   => null()
       integer :: status
       logical :: right_node_set, left_node_set
-      character(len=ESMF_MAXPATHLEN) :: left_file, right_file
 
       right_node_set = this%right_node%check_if_initialized(_RC)
       left_node_set = this%left_node%check_if_initialized(_RC)