Skip to content

Commit

Permalink
Addressing review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@abishekg7
abishekg7 committed Oct 22, 2024
1 parent 53d9332 commit 0c4e0cf
Showing 1 changed file with 35 additions and 31 deletions.
66 changes: 35 additions & 31 deletions src/core_init_atmosphere/mpas_init_atm_gwd.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ module mpas_init_atm_gwd
! A derived type to hold contents of a tile (both topo and landuse)
type :: mpas_gwd_tile_type

real (kind=R4KIND), dimension(:,:), pointer :: topo_array
integer (kind=I1KIND), dimension(:,:), pointer :: land_array
real (kind=R4KIND), dimension(:,:), pointer :: topo_array => null()
integer (kind=I1KIND), dimension(:,:), pointer :: landuse_array => null()
! coordinates of the tile to be read.
! NB: tile_start_x can be used as is to read landuse tiles, but need an
! adjustment to account for the shifting of topo array start_lon to -180.0.
integer :: tile_start_x, tile_start_y
integer :: tile_start_x = -1, tile_start_y = -1
! linked list next pointer
type (mpas_gwd_tile_type), pointer :: next => null()

Expand All @@ -51,14 +51,14 @@ subroutine read_geogrid(fname, rarray, nx, ny, nz, isigned, endian, &
end subroutine read_geogrid
end interface

real (kind=RKIND), parameter :: Re = 6371229.0_RKIND ! Earth radius in MPAS-Atmosphere
real (kind=RKIND), parameter :: Re = 6371229.0_RKIND ! Earth radius in MPAS-Atmosphere
real (kind=RKIND), parameter :: Pi = 2.0_RKIND * asin(1.0_RKIND)
real (kind=RKIND), parameter :: rad2deg = 180.0_RKIND / Pi

integer, parameter :: topo_x = 43200 ! x-dimension of global 30-arc-second topography array
integer, parameter :: topo_y = 21600 ! y-dimension of global 30-arc-second topography array
integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second topography
integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second topography
integer, parameter :: topo_x = 43200 ! x-dimension of global 30-arc-second topography array
integer, parameter :: topo_y = 21600 ! y-dimension of global 30-arc-second topography array
integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second topography
integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second topography

real (kind=RKIND), parameter :: pts_per_degree = real(topo_x,RKIND) / 360.0_RKIND

Expand Down Expand Up @@ -170,6 +170,7 @@ function compute_gwd_fields(domain) result(iErr)
case('GMTED2010')
call mpas_log_write('--- Using GMTED2010 terrain dataset for GWDO static fields')
geog_sub_path = 'topo_gmted2010_30s/'

! NB: the GMTED2010 data on disk actually has start_lon = 0.0, but the read_global_30s_topo()
! routine will shift the dataset when writing to the topo array so that the start_lon seen
! by the rest of this code is -180.0.
Expand Down Expand Up @@ -271,7 +272,7 @@ function compute_gwd_fields(domain) result(iErr)
call get_box_size_from_lat_lon(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc, nx, ny)

call get_box(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny, &
& geog_data_path, geog_sub_path, tilesHead, box, box_landuse, dxm, box_mean)
geog_data_path, geog_sub_path, tilesHead, box, box_landuse, dxm, box_mean)

!
! With a box of 30-arc-second data for the current grid cell, call
Expand Down Expand Up @@ -364,7 +365,6 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny)
!
ny = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re))


end subroutine get_box_size_from_lat_lon


Expand All @@ -387,12 +387,13 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result

implicit none

type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead
integer, intent(in) :: box_x, box_y
type(mpas_gwd_tile_type), pointer, intent(in) :: tilesHead
type(mpas_gwd_tile_type), pointer :: thisTile
character(len=*), intent(in) :: path
character(len=*), intent(in) :: sub_path

type(mpas_gwd_tile_type), pointer, intent(out) :: thisTile

integer :: tile_start_x, tile_start_y, tile_start_x_topo

! Need special handling for the x-coordinates of topo tiles, due to the shift by topo_shift
Expand Down Expand Up @@ -431,27 +432,30 @@ end function get_tile_from_box_point
!
! function add_tile
!
!> \brief Routine to read in a new topo and landmask tile
!> \brief Routine to read in a new topo and landuse tile, and add
!> these tiles to the head of the linked list tilesHead
!> \author Abishek Gopal
!> \date 05 Sep 2024
!> \details
!> Routine to read in a new topo and landmask tile, given the tile
!> Routine to read in a new topo and landuse tile, given the tile
! coordinates
!-----------------------------------------------------------------------
function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) result(newTile)

implicit none

integer :: iErr
type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead
integer, intent(in) :: tile_start_x, tile_start_y, tile_start_x_topo
type(mpas_gwd_tile_type), pointer :: tilesHead
character(len=*), intent(in) :: path
character(len=*), intent(in) :: sub_path
type(mpas_gwd_tile_type), pointer :: newTile

type(mpas_gwd_tile_type), pointer, intent(out) :: newTile

integer :: iErr

allocate(newTile)
allocate(newTile%topo_array(tile_x,tile_y))
allocate(newTile%land_array(tile_x,tile_y))
allocate(newTile%landuse_array(tile_x,tile_y))
newTile%tile_start_x = tile_start_x
newTile%tile_start_y = tile_start_y
newTile%next => tilesHead
Expand All @@ -462,7 +466,7 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path
return
end if

iErr = read_30s_landuse_tile(path, sub_path, newTile%land_array, newTile%tile_start_x, newTile%tile_start_y)
iErr = read_30s_landuse_tile(path, sub_path, newTile%landuse_array, newTile%tile_start_x, newTile%tile_start_y)
if (iErr /= 0) then
call mpas_log_write('Error reading global 30-arc-sec landuse for GWD statistics', messageType=MPAS_LOG_ERR)
return
Expand All @@ -488,19 +492,20 @@ function free_tile_list(tilesHead) result(iErr)

implicit none

integer :: iErr
type(mpas_gwd_tile_type), pointer :: tilesHead
type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead

integer :: iErr

type(mpas_gwd_tile_type), pointer :: thisTile

! loop over tiles
do while (associated(tilesHead))
thisTile => tilesHead
tilesHead => thisTile % next
deallocate(thisTile%topo_array)
deallocate(thisTile%land_array)
deallocate(thisTile%landuse_array)
deallocate(thisTile)
enddo ! associated(thisTile)
end do ! associated(thisTile)

iErr = 0

Expand All @@ -524,7 +529,7 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re

character(len=*), intent(in) :: path
character(len=*), intent(in) :: sub_path
real(kind=R4KIND), dimension(:,:), pointer :: topo
real(kind=R4KIND), dimension(:,:), pointer, intent(inout) :: topo
integer, intent(in) :: tile_start_x
integer, intent(in) :: tile_start_y

Expand Down Expand Up @@ -555,10 +560,8 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re
iy = tile_start_y
ix = tile_start_x

!write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//trim(sub_path), tile_start_x, '-', (tile_start_x+tile_x-1), '.', &
!tile_start_y, '-', (tile_start_y+tile_y-1)
write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//trim(sub_path), ix, '-', (ix+tile_x-1), '.', &
iy, '-', (iy+tile_y-1)
iy, '-', (iy+tile_y-1)
call mpas_f_to_c_string(filename, c_filename)
call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, &
wordsize, istatus)
Expand Down Expand Up @@ -595,7 +598,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start

character(len=*), intent(in) :: path
character(len=*), intent(in) :: sub_path
integer (kind=I1KIND), dimension(:,:), pointer :: landuse
integer (kind=I1KIND), dimension(:,:), pointer, intent(inout) :: landuse
integer, intent(in) :: tile_start_x
integer, intent(in) :: tile_start_y

Expand All @@ -621,7 +624,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start
nz = 1

write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//'/landuse_30s/', tile_start_x, '-', (tile_start_x+tile_x-1), '.', &
tile_start_y, '-', (tile_start_y+tile_y-1)
tile_start_y, '-', (tile_start_y+tile_y-1)
call mpas_f_to_c_string(filename, c_filename)
call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, &
wordsize, istatus)
Expand Down Expand Up @@ -709,7 +712,7 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse
integer, intent(in) :: nx, ny
character(len=*), intent(in) :: path
character(len=*), intent(in) :: sub_path
type(mpas_gwd_tile_type), pointer, intent(in) :: tilesHead
type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead
real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell
integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell
real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell
Expand Down Expand Up @@ -770,14 +773,15 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse
ix = (ii - thisTile%tile_start_x) + 1
jx = (jj - thisTile%tile_start_y) + 1
box(i,j) = thisTile%topo_array(ix, jx)
box_landuse(i,j) = thisTile%land_array(ix, jx)
box_landuse(i,j) = thisTile%landuse_array(ix, jx)
sg_lat = (start_lat + (real(jj-1,RKIND) + 0.5) / pts_per_degree) / rad2deg ! Add 0.5 for cell center
dxm(i,j) = sg_delta * cos(sg_lat)
box_mean = box_mean + box(i,j)

end do
end do


!
! Compute mean topography in the extracted box
!
Expand Down

0 comments on commit 0c4e0cf

Please sign in to comment.