diff --git a/base/MAPL_TripolarGridFactory.F90 b/base/MAPL_TripolarGridFactory.F90 index b629e33299ea..9024ac102bd1 100644 --- a/base/MAPL_TripolarGridFactory.F90 +++ b/base/MAPL_TripolarGridFactory.F90 @@ -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 @@ -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 @@ -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 diff --git a/base/MAPL_VerticalMethods.F90 b/base/MAPL_VerticalMethods.F90 index 6759ab8e3bef..21cf2a1e31fe 100644 --- a/base/MAPL_VerticalMethods.F90 +++ b/base/MAPL_VerticalMethods.F90 @@ -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 @@ -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 @@ -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