From 8236f7aac525ceac9bf3e29e804708b4c773abcb Mon Sep 17 00:00:00 2001 From: Lizzie Lundgren Date: Thu, 19 Sep 2024 10:32:03 -0400 Subject: [PATCH] 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 --- CHANGELOG.md | 1 + base/FileMetadataUtilities.F90 | 22 ++++++++++++++++++++- gridcomps/ExtData/ExtDataGridCompMod.F90 | 24 +++++++++++++++++++---- gridcomps/ExtData2G/ExtDataGridCompNG.F90 | 20 ++++++++++++++++--- gridcomps/ExtData2G/ExtDataTypeDef.F90 | 1 + 5 files changed, 60 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a5ceb79d3bf8..4356b8703ccd 100644 --- 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 index 695bffe89e5b..0f0bc749d54d 100644 --- 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 index 333605aec3ef..eaa705eba66f 100644 --- 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 index c1984c538e4c..7af3c01cc610 100644 --- 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 index f7f7ec75ded3..e3ad51654e24 100644 --- 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