Skip to content

Commit

Permalink
Merge pull request #8 from GEOS-ESM/bugfix/yvikhlya/tripolar_grid_fix
Browse files Browse the repository at this point in the history 
This fixes a bug which caused a crash when creating MOM grid in coupl…
  • Loading branch information
@tclune
tclune authored Aug 1, 2019
2 parents 07de57d + 7623b78 commit 32ed6b9
Show file tree
Hide file tree
Showing 2 changed files with 142 additions and 113 deletions.
3 changes: 3 additions & 0 deletions MAPL_Base/MAPL_GridManager.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,7 @@ end subroutine add_prototype

function make_clone(this, grid_type, unusable, rc) result(factory)
use MAPL_LatLonGridFactoryMod, only: LatLonGridFactory
use MAPL_TripolarGridFactoryMod, only: TripolarGridFactory
class (AbstractGridFactory), allocatable :: factory
class (GridManager), intent(inout) :: this
character(len=*), intent(in) :: grid_type
Expand All @@ -102,11 +103,13 @@ function make_clone(this, grid_type, unusable, rc) result(factory)
!---------------
logical, save :: initialized = .false.
type (LatLonGridFactory) :: latlon_factory
type (TripolarGridFactory) :: tripolar_factory

_UNUSED_DUMMY(unusable)

if (.not. initialized) then
call this%prototypes%insert('LatLon', latlon_factory)
call this%prototypes%insert('Tripolar',tripolar_factory)
initialized = .true.
end if

Expand Down
252 changes: 139 additions & 113 deletions MAPL_Base/MAPL_TripolarGridFactory.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ module MAPL_TripolarGridFactoryMod
character(len=*), parameter :: UNDEFINED_CHAR = '**'

character(len=*), parameter :: GRID_NAME_DEFAULT = 'UNKNOWN'
character(len=*), parameter :: GRID_FILE_NAME_DEFAULT = 'UNKNOWN'

type, extends(AbstractGridFactory) :: TripolarGridFactory
private
Expand Down Expand Up @@ -63,7 +64,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 All @@ -85,13 +85,16 @@ module MAPL_TripolarGridFactoryMod


function TripolarGridFactory_from_parameters(unusable, grid_file_name, grid_name, &
& nx, ny, rc) result(factory)
& im_world,jm_world,lm,nx, ny, rc) result(factory)
type (TripolarGridFactory) :: factory
class (KeywordEnforcer), optional, intent(in) :: unusable

! grid details:
character(len=*), intent(in) :: grid_file_name ! required
character(len=*), optional, intent(in) :: grid_name
integer, optional, intent(in) :: im_world
integer, optional, intent(in) :: jm_world
integer, optional, intent(in) :: lm

! decomposition:
integer, optional, intent(in) :: nx
Expand All @@ -100,20 +103,17 @@ function TripolarGridFactory_from_parameters(unusable, grid_file_name, grid_name

integer :: status
character(len=*), parameter :: Iam = MOD_NAME // 'TripolarGridFactory_from_parameters'
logical :: exists

if (present(unusable)) print*,shape(unusable)

call set_with_default(factory%grid_name, grid_name, GRID_NAME_DEFAULT)
call set_with_default(factory%grid_file_name, grid_file_name, GRID_FILE_NAME_DEFAULT)

call set_with_default(factory%ny, nx, UNDEFINED_INTEGER)
call set_with_default(factory%nx, ny, UNDEFINED_INTEGER)

factory%grid_file_name = grid_file_name
inquire(file=grid_file_name, exist=exists)
_ASSERT(exists)

call factory%read_grid_dimensions()
call set_with_default(factory%im_world, im_world, UNDEFINED_INTEGER)
call set_with_default(factory%jm_world, jm_world, UNDEFINED_INTEGER)
call set_with_default(factory%lm, lm, UNDEFINED_INTEGER)

call factory%check_and_fill_consistency(rc=status)
_VERIFY(status)
Expand Down Expand Up @@ -157,19 +157,20 @@ function create_basic_grid(this, unusable, rc) result(grid)

_UNUSED_DUMMY(unusable)

grid = ESMF_GridCreate( &
& name = this%grid_name, &
& countsPerDEDim1=this%ims, &
& countsPerDEDim2=this%jms, &
& indexFlag=ESMF_INDEX_DELOCAL, &
& gridEdgeLWidth=[0,0], &
& gridEdgeUWidth=[0,0], &
& coordDep1=[1,2], &
& coordDep2=[1,2], &
& rc=status)
grid = ESMF_GridCreate1PeriDim( &
name=trim(this%grid_name) ,&
countsPerDEDim1=this%ims, &
countsPerDEDim2=this%jms, &
indexFlag=ESMF_INDEX_DELOCAL, &
gridEdgeLWidth=[0,0], &
gridEdgeUWidth=[0,1], &
coordDep1=[1,2], &
coordDep2=[1,2], &
poleKindFlag=[ESMF_POLEKIND_MONOPOLE,ESMF_POLEKIND_BIPOLE], &
coordSys=ESMF_COORDSYS_SPH_RAD, rc=status)
_VERIFY(status)

! Allocate coords at default stagger location
!Allocate coords at default stagger location
call ESMF_GridAddCoord(grid, rc=status)
_VERIFY(status)

Expand All @@ -185,44 +186,87 @@ function create_basic_grid(this, unusable, rc) result(grid)
end function create_basic_grid

subroutine add_horz_coordinates(this, grid, unusable, rc)
use MAPL_BaseMod, only: MAPL_grid_interior
use MAPL_BaseMod, only: MAPL_grid_interior, MAPL_gridget
use MAPL_CommsMod
use MAPL_IOMod
use MAPL_ConstantsMod
class (TripolarGridFactory), intent(in) :: this
type (ESMF_Grid), intent(inout) :: grid
class (KeywordEnforcer), optional, intent(in) :: unusable
integer, optional, intent(out) :: rc

integer :: i_1, i_n, j_1, j_n ! regional array bounds
real(kind=ESMF_KIND_R8), pointer :: centers(:,:)
real(kind=ESMF_KIND_R8), allocatable :: longitudes(:,:)
real(kind=ESMF_KIND_R8), allocatable :: latitudes(:,:)
integer :: status
character(len=*), parameter :: Iam = MOD_NAME // 'add_horz_coordinates'

integer :: UNIT
integer :: IM, JM
integer :: IMSTART, JMSTART
integer :: IM_WORLD, JM_WORLD
integer :: DUMMYI, DUMMYJ

integer :: COUNTS(3), DIMS(3)
type(ESMF_DELayout) :: LAYOUT
type(ESMF_DistGrid) :: DISTGRID
real(ESMF_KIND_R8), allocatable :: x(:,:), y(:,:)
real(ESMF_KIND_R8), pointer :: gridx(:,:), gridy(:,:)

_UNUSED_DUMMY(unusable)
! get IM, JM and IM_WORLD, JM_WORLD
call MAPL_GridGet(GRID, localCellCountPerDim=COUNTS, globalCellCountPerDim=DIMS, RC=STATUS)
_VERIFY(STATUS)

IM = COUNTS(1)
JM = COUNTS(2)
IM_WORLD = DIMS(1)
JM_WORLD = DIMS(2)

! get global index of the lower left corner
!------------------------------------------
call MAPL_GRID_INTERIOR(GRID,IMSTART,DUMMYI,JMSTART,DUMMYJ)

call ESMF_GridGetCoord(grid, localDE=0, coordDim=1, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=gridx, rc=status)
_VERIFY(STATUS)

call ESMF_GridGetCoord(grid, localDE=0, coordDim=2, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=gridy, rc=status)
_VERIFY(STATUS)

allocate(x(IM_WORLD, JM_WORLD), stat=status)
_VERIFY(STATUS)
allocate(y(IM_WORLD, JM_WORLD), stat=status)
_VERIFY(STATUS)

call ESMF_GridGet (GRID, distGrid=distGrid, rc=STATUS)
_VERIFY(STATUS)
call ESMF_DistGridGet(distGRID, delayout=layout, rc=STATUS)
_VERIFY(STATUS)

UNIT = GETFILE(this%grid_file_name, form="formatted", rc=status)
call READ_PARALLEL(LAYOUT, X, unit=UNIT)
call READ_PARALLEL(LAYOUT, Y, unit=UNIT)
call FREE_FILE(UNIT)


! Make sure the longitudes are between -180 and 180 degrees
!ALT disable this for AR5 X = mod(X + 72180._8,360._8) - 180._8 ! -180<= lon0 <180
! Convert to radians
X = X * (MAPL_PI_R8)/180._8
Y = Y * (MAPL_PI_R8)/180._8


! Modify grid coordinates
!------------------------
GRIDX = X(IMSTART:IMSTART+IM-1,JMSTART:JMSTART+JM-1)
GRIDY = Y(IMSTART:IMSTART+IM-1,JMSTART:JMSTART+JM-1)

! Clean-up
!---------
deallocate(y)
deallocate(x)

call this%read_grid_coordinates(longitudes, latitudes)

call MAPL_grid_interior(grid, i_1, i_n, j_1, j_n)

! First we handle longitudes:
call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=centers, rc=status)
_VERIFY(status)

call ArrayScatter(centers, longitudes, grid, rc=status)
_VERIFY(status)

! Now latitudes
call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=centers, rc=status)
_VERIFY(status)
call ArrayScatter(centers, latitudes, grid, rc=status)
_VERIFY(status)

deallocate(longitudes, latitudes)
_RETURN(_SUCCESS)

end subroutine add_horz_coordinates
Expand All @@ -245,14 +289,15 @@ subroutine initialize_from_config_with_prefix(this, config, prefix, unusable, rc
call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRIDNAME:', default=GRID_NAME_DEFAULT)
this%grid_name = trim(tmp)

call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRID_FILE_NAME:', rc=status)
call ESMF_ConfigGetAttribute(config, tmp, label=prefix//'GRIDSPEC:', rc=status)
_VERIFY(status)
this%grid_file_name = trim(tmp)
call this%read_grid_dimensions()
!call this%read_grid_dimensions()

call ESMF_ConfigGetAttribute(config, this%nx, label=prefix//'NX:', default=UNDEFINED_INTEGER)
call ESMF_ConfigGetAttribute(config, this%ny, label=prefix//'NY:', default=UNDEFINED_INTEGER)

call ESMF_ConfigGetAttribute(config, this%im_world, label=prefix//'IM_WORLD:', default=UNDEFINED_INTEGER)
call ESMF_ConfigGetAttribute(config, this%jm_world, label=prefix//'JM_WORLD:', default=UNDEFINED_INTEGER)
call ESMF_ConfigGetAttribute(config, this%lm, label=prefix//'LM:', default=UNDEFINED_INTEGER)

call this%check_and_fill_consistency(rc=status)
Expand Down Expand Up @@ -324,7 +369,7 @@ subroutine check_and_fill_consistency(this, unusable, rc)
integer, optional, intent(out) :: rc

character(len=*), parameter :: Iam = MOD_NAME // 'check_and_fill_consistency'

integer :: status
_UNUSED_DUMMY(unusable)

if (.not. allocated(this%grid_name)) then
Expand All @@ -337,10 +382,51 @@ subroutine check_and_fill_consistency(this, unusable, rc)
_ASSERT(mod(this%jm_world, this%ny) == 0)

! local extents
this%ims = spread(this%im_world / this%nx, 1, this%nx)
this%jms = spread(this%jm_world / this%ny, 1, this%ny)
call verify(this%nx, this%im_world, this%ims, rc=status)
call verify(this%ny, this%jm_world, this%jms, rc=status)
!this%ims = spread(this%im_world / this%nx, 1, this%nx)
!this%jms = spread(this%jm_world / this%ny, 1, this%ny)

_RETURN(_SUCCESS)

contains

subroutine verify(n, m_world, ms, rc)
integer, intent(inout) :: n
integer, intent(inout) :: m_world
integer, allocatable, intent(inout) :: ms(:)
integer, optional, intent(out) :: rc

integer :: status

if (allocated(ms)) then
_ASSERT(size(ms) > 0)

if (n == UNDEFINED_INTEGER) then
n = size(ms)
else
_ASSERT(n == size(ms))
end if

if (m_world == UNDEFINED_INTEGER) then
m_world = sum(ms)
else
_ASSERT(m_world == sum(ms))
end if

else

_ASSERT(n /= UNDEFINED_INTEGER)
_ASSERT(m_world /= UNDEFINED_INTEGER)
allocate(ms(n), stat=status)
_VERIFY(status)
call MAPL_DecomposeDim(m_world, ms, n)

end if

_RETURN(_SUCCESS)

end subroutine verify

end subroutine check_and_fill_consistency

Expand Down Expand Up @@ -444,66 +530,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
Expand Down

0 comments on commit 32ed6b9

Please sign in to comment.