Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixes tripolar metadata ouput #528

Merged
merged 2 commits into from
Aug 26, 2020
Merged

Fixes tripolar metadata ouput #528

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6c16d9d
add 2D lons/lats to tripolar output, fix so you can specify a file wi…
bena-nasa Aug 26, 2020
9a75525
removed what appears to be used routines. Can't find any reference to…
bena-nasa Aug 26, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 12 additions & 121 deletions base/MAPL_TripolarGridFactory.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,6 @@ module MAPL_TripolarGridFactoryMod

procedure :: equals

procedure :: read_grid_dimensions
procedure :: read_grid_coordinates
procedure :: check_and_fill_consistency
procedure :: generate_grid_name
procedure :: to_string
Expand Down Expand Up @@ -556,125 +554,6 @@ function generate_grid_name(this) result(name)

end function generate_grid_name


subroutine read_grid_coordinates(this, longitudes, latitudes, unusable, rc)
class (TripolarGridFactory), intent(in) :: this
real(kind=REAL64), allocatable, intent(out) :: longitudes(:,:)
real(kind=REAL64), allocatable, intent(out) :: latitudes(:,:)
class (KeywordEnforcer), optional, intent(out) :: unusable
integer, optional, intent(out) :: rc

include 'netcdf.inc'

integer :: status
character(len=*), parameter :: Iam = MOD_NAME // 'read_grid_coordinates()'

integer :: xid, yid
integer :: start(2), counts(2)
integer :: pet, ndes
logical :: i_am_root
integer :: ncid
type (ESMF_VM) :: vm

real(kind=REAL64), allocatable :: lons(:,:), lats(:,:)

_UNUSED_DUMMY(unusable)

call ESMF_VMGetCurrent(vm, rc=status)
_VERIFY(status)

call ESMF_VMGet(vm, localpet=pet, petCount=ndes, rc=status)
_VERIFY(status)

i_am_root = (pet == 0)

if (i_am_root) then
allocate(longitudes(this%im_world, this%jm_world), stat=status)
_VERIFY(status)
allocate(latitudes(this%im_world, this%jm_world), stat=status)
_VERIFY(status)

ncid = ncopn(this%grid_file_name, NCNOWRIT, status)
_VERIFY(status)

xid = ncvid(ncid, 'x_T', status)
_VERIFY(status)

yid = ncvid(ncid, 'y_T', status)
_VERIFY(status)

call ncvgt(ncid, xid, start, counts, lons, status)
_VERIFY(status)
call ncvgt(ncid, yid, start, counts, lats, status)
_VERIFY(status)

call ncclos(ncid, status)
_VERIFY(status)
else
allocate(longitudes(0,0))
allocate(latitudes(0,0))
end if

end subroutine read_grid_coordinates

subroutine read_grid_dimensions(this, unusable, rc)
use MAPL_CommsMod
class (TripolarGridFactory), intent(inout) :: this
class (KeywordEnforcer), optional, intent(out) :: unusable
integer, optional, intent(out) :: rc

include 'netcdf.inc'

integer :: status
character(len=*), parameter :: Iam = MOD_NAME // 'read_grid_dimensions()'

integer :: xid
character(len=128) :: name
integer :: type
integer :: n
integer :: dims(MAXVDIMS)
integer :: natt
integer :: pet, ndes
logical :: i_am_root
integer :: ncid
type (ESMF_VM) :: vm

_UNUSED_DUMMY(unusable)

call ESMF_VMGetCurrent(vm, rc=status)
_VERIFY(status)

call ESMF_VMGet(vm, localpet=pet, petCount=ndes, rc=status)
_VERIFY(status)

i_am_root = (pet == 0)

if (i_am_root) then
ncid = ncopn(this%grid_file_name, NCNOWRIT, status)
_VERIFY(status)

xid = ncvid(ncid, 'x_T', status)
_VERIFY(status)

call ncvinq (ncid, xid, name, type, n, dims, natt, status)
_VERIFY(status)

associate (im => this%im_world, jm => this%jm_world)
call ncdinq(ncid, dims(1), name, im, status)
_VERIFY(status)
call ncdinq(ncid, dims(2), name, jm, status)
_VERIFY(status)
end associate
end if

call MAPL_CommsBCast(vm, this%im_world, 1, root=0, rc=status)
_VERIFY(status)
call MAPL_CommsBCast(vm, this%jm_world, 1, root=0, rc=status)
_VERIFY(status)


end subroutine read_grid_dimensions

subroutine init_halo(this, unusable, rc)
class (TripolarGridFactory), intent(inout) :: this
class (KeywordEnforcer), optional, intent(in) :: unusable
Expand Down Expand Up @@ -901,9 +780,21 @@ subroutine append_metadata(this, metadata)
class (TripolarGridFactory), intent(inout) :: this
type (FileMetadata), intent(inout) :: metadata

type (Variable) :: v

call metadata%add_dimension('lon', this%im_world)
call metadata%add_dimension('lat', this%jm_world)

v = Variable(PFIO_REAL32, dimensions='lon,lat')
call v%add_attribute('long_name','longitude')
call v%add_attribute('units','degrees_east')
call metadata%add_variable('lons',v)

v = Variable(PFIO_REAL32, dimensions='lon,lat')
call v%add_attribute('long_name','latitude')
call v%add_attribute('units','degrees_north')
call metadata%add_variable('lats',v)


end subroutine append_metadata

Expand Down
72 changes: 43 additions & 29 deletions base/MAPL_VerticalMethods.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ function newVerticalData(levels,vcoord,vscale,vunit,rc) result(vdata)
else
vdata%interp_levels = vdata%scaled_levels
endif
else
else
vdata%regrid_type = VERTICAL_METHOD_SELECT
end if
end function newVerticalData
Expand Down Expand Up @@ -291,7 +291,7 @@ subroutine append_vertical_metadata(this,metadata,bundle,rc)
integer :: ungrdsize
real, allocatable :: ungridded_coord(:)
real, allocatable :: ungridded_coords(:,:)
logical :: unGrdNameCheck, unGrdUnitCheck, unGrdCoordCheck, have_ungrd
logical :: unGrdNameCheck, unGrdUnitCheck, unGrdCoordCheck, have_ungrd, found_mixed_ce

integer :: status
type(Variable) :: v
Expand Down Expand Up @@ -368,45 +368,59 @@ subroutine append_vertical_metadata(this,metadata,bundle,rc)
end do
end if

found_mixed_ce=.false.
lm=1
haveVert = any(varDims/=1)
if (haveVert) then
if (this%regrid_type == VERTICAL_METHOD_NONE ) then
vlast=1
do i=1,numVars
if (varDims(i)/=1) then
if (vlast/=1) then
if (vlast /= varDims(i)) then
_ASSERT(.false.,"have mixed level/edge")
end if
else
vlast=varDims(i)
vloc=location(i)
vlast=1
do i=1,numVars
if (varDims(i)/=1) then
if (vlast/=1) then
if (vlast /= varDims(i)) then
found_mixed_ce=.true.
end if
else
vlast=varDims(i)
vloc=location(i)
end if
end do
lm=vlast
if (vloc == MAPL_VLocationCenter) then
vlb = 1
else if (vloc == MAPL_VlocationEdge) then
vlb = 0
else
vlb = 1
end if
else if (this%regrid_type == VERTICAL_METHOD_ETA2LEV .or. &
this%regrid_type == VERTICAL_METHOD_SELECT) then
lm = size(this%levs)
end if
end do
lm=vlast
if (vloc == MAPL_VLocationCenter) then
vlb = 1
else if (vloc == MAPL_VlocationEdge) then
vlb = 0
else
vlb = 1
end if
end if

if (this%regrid_type == VERTICAL_METHOD_ETA2LEV) then
lm = size(this%levs)
vlb = 1
end if
if (this%regrid_type == VERTICAL_METHOD_SELECT) then
if (lm == size(this%levs)) then
this%regrid_type = VERTICAL_METHOD_NONE
else
lm = size(this%levs)
vlb=1
end if
end if
if (this%regrid_type == VERTICAL_METHOD_NONE) then
_ASSERT(.not.(found_mixed_ce),'have mixed level/edge')
end if


if (haveVert) then
this%lm=lm
if (this%regrid_type == VERTICAL_METHOD_NONE) then
allocate(this%levs(lm))
do i=1,lm
this%levs(i)=vlb+i-1
enddo
if (.not.allocated(this%levs)) then
allocate(this%levs(lm))
do i=1,lm
this%levs(i)=vlb+i-1
enddo
end if
if (have_ungrd) then
if (allocated(ungridded_coord)) then
this%levs=ungridded_coord
Expand Down