Add dimensionless vertical coordinate to criteria for import vertical…

… flip Some dimensionless vertical coordinates are proxies for pressure and should therefore have the direction of their values used to determine if the import should be vertically flipped. Signed-off-by: Lizzie Lundgren <[email protected]>
geoschem · Sep 19, 2024 · 8236f7a · 8236f7a
1 parent 231d53c
commit 8236f7a
Show file tree

Hide file tree

Showing 5 changed files with 60 additions and 8 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Add debug logger call to print out the name of the container being read from `ExtData.rc`
+- Added dimensionless vertical coordinates to criteria for vertically flipping imports in ExtData
 
 ### Changed
 

diff --git a/base/FileMetadataUtilities.F90 b/base/FileMetadataUtilities.F90
@@ -490,20 +490,23 @@ function get_variable_attribute(this,var_name,attr_name,rc) result(units)
 
    end function get_variable_attribute
 
-   subroutine get_coordinate_info(this,coordinate_name,coordSize,coordUnits,coords,rc)
+   subroutine get_coordinate_info(this,coordinate_name,coordSize,coordUnits,coordStandardName,coords,rc)
       class (FileMetadataUtils), intent(inout) :: this
       character(len=*), intent(in) :: coordinate_name
       integer, optional, intent(out) :: coordSize
       character(len=*), optional, intent(out) :: coordUnits
+      character(len=*), optional, intent(out) :: coordStandardName
       real, allocatable, optional,  intent(inout) :: coords(:)
       integer, optional, intent(out) :: rc
 
       integer :: status
+      logical :: isPresent
       character(:), allocatable :: fname
       class(CoordinateVariable), pointer :: var
       type(Attribute), pointer :: attr
       character(len=:), pointer :: vdim
       class(*), pointer :: coordUnitPtr
+      class(*), pointer :: coordStandardNamePtr
       class(*), pointer :: ptr(:)
 
       fname = this%get_file_name(_RC)
@@ -526,6 +529,23 @@ subroutine get_coordinate_info(this,coordinate_name,coordSize,coordUnits,coords,
          end select
       end if 
 
+      if (present(coordStandardName)) then
+         ! Allow for missing standard_name
+         isPresent = var%is_attribute_present('standard_name')
+         if ( isPresent) then
+            attr => var%get_attribute('standard_name')
+            coordStandardNamePtr => attr%get_value()
+            select type(coordStandardNamePtr)
+            type is (character(*))
+               coordStandardName = trim(coordStandardNamePtr)
+            class default
+               _FAIL('coordinate standard name must be string in '//fname)
+            end select
+         else
+            coordStandardName = 'missing'
+         endif
+      end if
+
       if (present(coords)) then
          ptr => var%get_coordinate_data()
          _ASSERT(associated(ptr),"coord variable coordinate data not found in "//fname)

diff --git a/gridcomps/ExtData/ExtDataGridCompMod.F90 b/gridcomps/ExtData/ExtDataGridCompMod.F90
@@ -154,6 +154,7 @@ MODULE MAPL_ExtDataGridCompMod
      character(len=4)             :: importVDir = "down"
      character(len=4)             :: fileVDir = "down"
      character(len=ESMF_MAXSTR)   :: levUnit
+     character(len=ESMF_MAXSTR)   :: levStandardName
      logical                      :: havePressure = .false.
   end type PrimaryExport
 
@@ -2087,6 +2088,7 @@ subroutine GetLevs(item, time, state, allowExtrap, rc)
         type(ESMF_Field)           :: field
         real, allocatable          :: levFile(:)
         character(len=ESMF_MAXSTR) :: buff,levunits,tlevunits,temp_name
+        character(len=ESMF_MAXSTR) :: levstandardname,tlevstandardname
         logical                    :: found,lFound,intOK
         integer                    :: maxOffset
         character(len=:), allocatable :: levname
@@ -2194,15 +2196,29 @@ subroutine GetLevs(item, time, state, allowExtrap, rc)
         levName = metadata%get_level_name(rc=status)
         _VERIFY(status)
         if (trim(levName) /='') then
-           call metadata%get_coordinate_info(levName,coordSize=item%lm,coordUnits=tLevUnits,coords=levFile,__RC__)
+           call metadata%get_coordinate_info(levName,coordSize=item%lm,coordUnits=tLevUnits, &
+                coordStandardName=tLevStandardName,coords=levFile,__RC__)
            ! GCHP: use trim to avoid gfortran issue
            !levUnits=MAPL_TrimString(tlevUnits)
            levUnits=trim(tlevUnits)
-           ! check if pressure
+           ! check if vertical coordinate is pressure or dimensionless pressure proxy
            item%levUnit = ESMF_UtilStringLowerCase(levUnits)
-           if (trim(item%levUnit) == 'hpa' .or. trim(item%levUnit)=='pa') then
+           if ( ( trim(item%levUnit) == 'hpa' ) .or. &
+                ( trim(item%levUnit) =='pa' ) .or. &
+                ( trim(item%levUnit) == 'sigma_level' ) ) then
               item%havePressure = .true.
-           end if
+           elseif ( trim(levStandardName) /= 'missing' ) then
+              ! check standard_name attribute to catch cases where lev values are
+              ! dimensionless vertical coordinates as proxy for pressure
+              levStandardName= trim(tlevStandardName)
+              item%levStandardName = ESMF_UtilStringLowerCase(levStandardName)
+              if ( ( index(trim(item%levStandardName),"atmosphere_hybrid_sigma_pressure_coordinate") > 0 ) .or. &
+                   ( index(trim(item%levStandardName),"atmosphere_sigma_coordinate") > 0 ) .or. &
+                   ( index(trim(item%levStandardName),"atmosphere_ln_pressure_coordinate") > 0 ) ) then
+                 item%havePressure = .true.
+              endif
+           endif
+
            if (item%havePressure) then
               if (levFile(1)>levFile(size(levFile))) item%fileVDir="up"
            else

diff --git a/gridcomps/ExtData2G/ExtDataGridCompNG.F90 b/gridcomps/ExtData2G/ExtDataGridCompNG.F90
@@ -924,6 +924,7 @@ subroutine GetLevs(item, rc)
 
         real, allocatable          :: levFile(:)
         character(len=ESMF_MAXSTR) :: levunits,tlevunits
+        character(len=ESMF_MAXSTR) :: levstandardname,tlevstandardname
         character(len=:), allocatable :: levname
         character(len=:), pointer :: positive
         type(Variable), pointer :: var
@@ -945,12 +946,25 @@ subroutine GetLevs(item, rc)
         levName = item%file_metadata%get_level_name(rc=status)
         _VERIFY(status)
         if (trim(levName) /='') then
-           call item%file_metadata%get_coordinate_info(levName,coordSize=item%lm,coordUnits=tLevUnits,coords=levFile,_RC)
+           call item%file_metadata%get_coordinate_info(levName,coordSize=item%lm,coordUnits=tLevUnits, &
+                coordStandardName=tLevStandardName,coords=levFile,_RC)
            levUnits=MAPL_TrimString(tlevUnits)
-           ! check if pressure
+           ! check if vertical coordinate is pressure or dimensionless pressure proxy
            item%levUnit = ESMF_UtilStringLowerCase(levUnits)
-           if (trim(item%levUnit) == 'hpa' .or. trim(item%levUnit)=='pa') then
+           if ( ( trim(item%levUnit) == 'hpa' ) .or. &
+                ( trim(item%levUnit) == 'pa' ) .or. &
+                ( trim(item%levUnit) == 'sigma_level' ) ) then
               item%havePressure = .true.
+           elseif ( trim(levStandardName) /= 'missing' ) then
+              ! check standard_name attribute to catch cases where lev values are
+              ! dimensionless vertical coordinates as proxy for pressure
+              levStandardName= trim(tlevStandardName)
+              item%levStandardName = ESMF_UtilStringLowerCase(levStandardName)
+              if ( ( index(trim(item%levStandardName),"atmosphere_hybrid_sigma_pressure_coordinate") > 0 ) .or. &
+                   ( index(trim(item%levStandardName),"atmosphere_sigma_coordinate") > 0 ) .or. &
+                   ( index(trim(item%levStandardName),"atmosphere_ln_pressure_coordinate") > 0 ) ) then
+                 item%havePressure = .true.
+              endif
            end if
            if (item%havePressure) then
               if (levFile(1)>levFile(size(levFile))) item%fileVDir="up"

diff --git a/gridcomps/ExtData2G/ExtDataTypeDef.F90 b/gridcomps/ExtData2G/ExtDataTypeDef.F90
@@ -62,6 +62,7 @@ module MAPL_ExtDataTypeDef
      character(len=4)             :: importVDir = "down"
      character(len=4)             :: fileVDir = "down"
      character(len=ESMF_MAXSTR)   :: levUnit
+     character(len=ESMF_MAXSTR)   :: levStandardName
      logical                      :: havePressure = .false.
      type(ExtDataPointerUpdate) :: update_freq