From 6de7c5c4588abba9c045d1ce8432eeb2bca828c4 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 23 Aug 2024 16:11:50 +0000 Subject: [PATCH 01/27] First working version of parallel topo read to reduce mem usage. Some statistics are checked to be correct --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 577 +++++++++++++++++-- 1 file changed, 522 insertions(+), 55 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 25ef93c8c6..5889dabbb9 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -43,6 +43,9 @@ end subroutine read_geogrid 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 ! The following are set at the beginning of the compute_gwd_fields routine depending @@ -54,10 +57,12 @@ end subroutine read_geogrid real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) real (kind=R4KIND), dimension(:,:), pointer :: topo ! Global 30-arc-second topography + real (kind=R4KIND), dimension(:), pointer :: topo_new ! Global 30-arc-second topography real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer :: box_new ! Subset of topography covering a grid cell real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell - real (kind=RKIND) :: box_mean ! Mean value of topography in box - integer :: nx, ny ! Dimensions of box covering grid cell + real (kind=RKIND) :: box_mean, box_mean_new ! Mean value of topography in box + !integer :: nx, ny ! Dimensions of box covering grid cell integer (kind=I1KIND), dimension(:,:), pointer :: landuse ! Global 30-arc-second landuse integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell @@ -68,6 +73,12 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height + integer, dimension(:), pointer :: local_tile_x ! + integer, dimension(:), pointer :: local_tile_y ! + integer, dimension(:), pointer :: local_box_x ! + integer, dimension(:), pointer :: local_box_y ! + integer:: max_local_tiles, cum_local_tiles, num_local_box_points + contains @@ -116,12 +127,13 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND), pointer :: config_topo_data character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - + real (kind=RKIND) :: var2d_tmp, con_tmp!, oa1, oa2, oa3, oa4, ol1, ol2, ol3, ol4 ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse real (kind=RKIND) :: sum_var - + integer :: nx, ny + allocate(topo(topo_x,topo_y)) allocate(landuse(topo_x,topo_y)) @@ -192,13 +204,29 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'gamma', hgamma) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) - allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells + ! + ! Get number of points to extract in the zonal direction + ! + ! if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + ! nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) + ! else + ! nx = topo_x / 2 + ! end if - iErr = read_global_30s_topo(geog_data_path, geog_sub_path) - if (iErr /= 0) then - call mpas_log_write('Error reading global 30-arc-sec topography for GWD statistics', messageType=MPAS_LOG_ERR) - return - end if + ! ! + ! ! Get number of points to extract in the meridional direction + ! ! + ! ny = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re)) + + nx = 22000 + ny = 300 + + allocate(local_box_x(nx*ny)) + allocate(local_box_y(nx*ny)) + + allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells + + iErr = read_global_30s_landuse(geog_data_path) if (iErr /= 0) then @@ -225,11 +253,56 @@ function compute_gwd_fields(domain) result(iErr) call mpas_log_write('in the computation of GWD static fields.') end if + !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius) + + max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) + call mpas_log_write(' hi from compute_gwd_fields -- before allocate max_local_tiles:$i', intArgs=(/max_local_tiles/)) + allocate(local_tile_x(max_local_tiles)) + allocate(local_tile_y(max_local_tiles)) + ! ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. ! + cum_local_tiles=0 + do iCell=1,nCells + !do iCell=5,65 + + ! + ! First, get an estimate of the mean diameter (in meters) of the grid + ! cell by averaging the distances to each of the neighboring cells + ! + dc = 0.0 + do i=1,nEdgesOnCell(iCell) + dc = dc + dcEdge(edgesOnCell(i,iCell)) + end do + dc = dc / real(nEdgesOnCell(iCell),RKIND) + if (onUnitSphere) then + dc = dc * sphere_radius + end if + dc = dc * config_gwd_cell_scaling + + ! + call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) + + !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg/), intArgs=(/nx,ny/)) + + call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) + + !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius, dc) + iErr = determine_complete_tile_set(nx, ny) + + if (cum_local_tiles >= max_local_tiles) exit + end do + + iErr = read_global_30s_topo(geog_data_path, geog_sub_path) + if (iErr /= 0) then + call mpas_log_write('Error reading global 30-arc-sec topography for GWD statistics', messageType=MPAS_LOG_ERR) + return + end if + do iCell=1,nCells + !do iCell=5,65 ! ! First, get an estimate of the mean diameter (in meters) of the grid @@ -245,6 +318,13 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling + + !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg/), intArgs=(/nx,ny/)) + + !call get_box_points(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc) + call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) + !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius, dc) + !iErr = determine_complete_tile_set() ! ! Cut out a rectangular piece of the global 30-arc-second topography ! data that is centered at the lat/lon of the current cell being @@ -255,30 +335,37 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc) + call get_box(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) + + call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) + call get_box_new(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) ! ! With a box of 30-arc-second data for the current grid cell, call ! subroutines to compute each sub-grid orography statistic ! - var2d(iCell) = get_var() - con(iCell) = get_con() - oa1(iCell) = get_oa1() - oa2(iCell) = get_oa2() - oa3(iCell) = get_oa3() - oa4(iCell) = get_oa4() + var2d(iCell) = get_var(nx, ny) + var2d_tmp = get_var_new(nx, ny) + call mpas_log_write('iCell = $i .... var2d: $r var2d_tmp: $r', intArgs=(/iCell/), realArgs=(/var2d(iCell), var2d_tmp/)) + con(iCell) = get_con(nx, ny) + con_tmp = get_con_new(nx, ny) + call mpas_log_write('iCell = $i .... con: $r con_tmp: $r', intArgs=(/iCell/), realArgs=(/con(iCell), con_tmp/)) + oa1(iCell) = get_oa1(nx, ny) + oa2(iCell) = get_oa2(nx, ny) + oa3(iCell) = get_oa3(nx, ny) + oa4(iCell) = get_oa4(nx, ny) ! Critical height, to be used in OL computation ! See Appendix of Kim, Y-J, 1996: Representation of Sub-Grid Scale Orographic Effects ! in a General Circulation Model. J. Climate, 9, 2698-2717. hc = 1116.2_RKIND - 0.878_RKIND * var2d(iCell) - ol1(iCell) = get_ol1() - ol2(iCell) = get_ol2() - ol3(iCell) = get_ol3() - ol4(iCell) = get_ol4() + ol1(iCell) = get_ol1(nx, ny) + ol2(iCell) = get_ol2(nx, ny) + ol3(iCell) = get_ol3(nx, ny) + ol4(iCell) = get_ol4(nx, ny) - hlanduse(iCell) = get_dom_landmask() ! get dominant land mask in cell + hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell ! elvmax(iCell) = get_elvmax() ! htheta(iCell) = get_htheta() @@ -287,6 +374,7 @@ function compute_gwd_fields(domain) result(iErr) end do + ! Smooth variance at isolated points do iCell = 1,nCells sum_landuse = 0_I1KIND @@ -314,6 +402,140 @@ function compute_gwd_fields(domain) result(iErr) end function compute_gwd_fields + subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) + + implicit none + + real (kind=RKIND), intent(in) :: lat + real (kind=RKIND), intent(in) :: lon + real (kind=RKIND), intent(in) :: dx + + integer, intent(out) :: nx + integer, intent(out) :: ny + + ! + ! Get number of points to extract in the zonal direction + ! + if (cos(lat) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat))) + else + nx = topo_x / 2 + end if + + ! + ! Get number of points to extract in the meridional direction + ! + ny = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re)) + + + end subroutine get_box_size_from_lat_lon + + + + function unique_tile(ix, jx, len) result (is_unique) + + implicit none + + integer, intent(in) :: ix + integer, intent(in) :: jx + integer, intent(in) :: len + + integer :: i + logical :: is_unique + + do i = 1, len + if(local_tile_x(i)==ix .and. local_tile_y(i)==jx) then + is_unique = .false. + return + end if + end do + + is_unique = .true. + + end function unique_tile + + function locate_tile_in_set(tile_index_i, tile_index_j) result (tile_index) + + implicit none + + integer, intent(in) :: tile_index_i + integer, intent(in) :: tile_index_j + + integer :: itile, tile_index + + do itile = 1, cum_local_tiles + if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then + tile_index = itile + return + end if + end do + + tile_index = -1 + + end function locate_tile_in_set + + function determine_complete_tile_set(nx, ny) result(iErr) + + implicit none + !integer, intent(in) :: nCells + !REAL(kind=RKIND), DIMENSION(:), POINTER, intent(in) :: latCell + !REAL(kind=RKIND), DIMENSION(:), POINTER, intent(in) :: lonCell + integer :: iErr + integer, intent(in) :: nx, ny + + !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 :: ix, jx, ii, l,i, iCell, itile, jtile + + !max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) + !call mpas_log_write(' hi from determine_complete_tile_set -- before allocate max_local_tiles:$i', intArgs=(/max_local_tiles/)) + !allocate(local_tile_x(max_local_tiles)) + !allocate(local_tile_y(max_local_tiles)) + + + !call mpas_log_write(' hi from determine_complete_tile_set -- after allocate') + !do iCell=1,nCells + + ! + ! Cut out a rectangular piece of the global 30-arc-second topography + ! data that is centered at the lat/lon of the current cell being + ! processed and that is just large enough to cover the cell. The + ! rectangular array of topography data is stored in the module + ! variable 'box', and the dimensions of this array are given by the + ! module variables 'nx' and 'ny'. The get_box() routine also + ! computes the mean elevation in the array and stores that value in + ! the module variable 'box_mean'. + ! + !call get_box_points(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc) + + !ix = nint((lonCell(iCell)*rad2deg) * pts_per_degree) + !jx = nint((latCell(iCell)*rad2deg - start_lat) * pts_per_degree) + ii = cum_local_tiles + do l=1,nx*ny + itile = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 + jtile = floor( real(local_box_y(l) - 1) / real(tile_y)) * tile_y + 1 + !itile = (ix / tile_x) * tile_x + 1 + !jtile = (jx / tile_y) * tile_y + 1 + !call mpas_log_write(' Checking if point ($i, $i) with lon lat ($i, $i) is already present', intArgs=(/ix, jx, nint(lonCell(iCell)*rad2deg), nint(latCell(iCell)*rad2deg)/)) + !call mpas_log_write(' Checking if local box point ($i, $i) in tile ($i, $i) is already present', intArgs=(/local_box_x(l), local_box_y(l), itile, jtile/)) + if(unique_tile(itile, jtile, ii)) then + ii = ii + 1 + local_tile_x(ii) = itile + local_tile_y(ii) = jtile + !call mpas_log_write('---- Adding tile to local tile list ($i, $i)', intArgs=(/local_tile_x(ii),local_tile_y(ii)/)) + !call mpas_log_write('---- local box ($i, $i)', intArgs=(/local_box_x(l),local_box_y(l)/)) + end if + end do + + cum_local_tiles = ii + + !end do + + call mpas_log_write(' hi from end of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/ii, max_local_tiles/)) + ! + end function determine_complete_tile_set + !*********************************************************************** ! @@ -341,12 +563,14 @@ function read_global_30s_topo(path, sub_path) result(iErr) integer, parameter :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile integer (c_int) :: istatus - integer :: ix, iy, ishift, ix_shift + integer :: ix, iy, ishift, ix_shift, il + integer(kind=I8KIND):: ix_shift_long integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile type (c_ptr) :: tile_ptr character(len=StrKIND) :: filename + character(len=StrKIND):: message character(kind=c_char), dimension(StrKIND+1) :: c_filename allocate(tile(tile_x+2*tile_bdr,tile_y+2*tile_bdr,1)) @@ -362,6 +586,9 @@ function read_global_30s_topo(path, sub_path) result(iErr) ishift = 0 + allocate(topo_new(tile_x*tile_y*cum_local_tiles)) + call mpas_log_write('---- topo_x:$i topo_y: $i', intArgs=(/topo_x,topo_y/)) + call mpas_log_write('---- tile_x:$i tile_y: $i', intArgs=(/tile_x,tile_y/)) ! ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 @@ -370,6 +597,31 @@ function read_global_30s_topo(path, sub_path) result(iErr) ishift = topo_x / 2 end if + do il=1,cum_local_tiles + ix = local_tile_x(il) + iy = local_tile_y(il) + 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) + call mpas_f_to_c_string(filename, c_filename) + call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & + wordsize, istatus) + tile(:,:,:) = tile(:,:,:) * scalefactor + if (istatus /= 0) then + call mpas_log_write('Error reading topography tile '//trim(filename), messageType=MPAS_LOG_ERR) + iErr = 1 + return + end if + + !ix_shift = mod((ix-1) + ishift, topo_x) + 1 + ix_shift = (il-1)*tile_x*tile_y + 1 + call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,ishift/)) + !call mpas_log_write(' topo_x,y ($i,$i) ', intArgs=(/topo_x, topo_y/)) + !write(message,fmt='(A,i18,i18,A,i18,A,i18)') 'Reading tile (',il,cum_local_tiles,') topo_new start=', ix_shift_long,' end=',ix_shift_long+tile_x*tile_y-1 + !call mpas_log_write(message) + topo_new(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + + end do + do iy=1,topo_y,tile_y do ix=1,topo_x,tile_x write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//trim(sub_path), ix, '-', (ix+tile_x-1), '.', & @@ -476,9 +728,10 @@ end function read_global_30s_landuse !> \details 1 = land, 0 = water ! !----------------------------------------------------------------------- - integer (kind=I1KIND) function get_dom_landmask( ) + integer (kind=I1KIND) function get_dom_landmask(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j real (kind=RKIND) :: xland @@ -503,6 +756,65 @@ integer (kind=I1KIND) function get_dom_landmask( ) end function get_dom_landmask + subroutine get_box_points(lat, lon, nx, ny) + + implicit none + real (kind=RKIND), intent(in) :: lat, lon + + integer, intent(in) :: nx, ny + + integer :: i, j, ii, jj, ic, jc, l + real (kind=RKIND) :: sg_lat + + !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/lon,lat/), intArgs=(/nx,ny/)) + ! + ! + ! Find coordinates in global topography array of the box center + ! + ic = nint((lon) * pts_per_degree) + 1 + jc = nint((lat - start_lat) * pts_per_degree) + 1 + + if (ic <= 0) ic = ic + topo_x + if (ic > topo_x) ic = ic - topo_x + + ! + ! Extract sub-array (box) from global array; must properly account for + ! the periodicity in the longitude coordinate, as well as the poles + ! + l = 0 + do j=1,ny + do i=1,nx + + ii = i - nx/2 + ic + jj = j - ny/2 + jc + + if (jj <= 0) then + jj = -jj + 1 + ii = ii + topo_y + end if + if (jj > topo_y) then + jj = topo_y - (jj - topo_y - 1) + ii = ii + topo_y + end if + do while (ii <= 0) + ii = ii + topo_x + end do + do while (ii > topo_x) + ii = ii - topo_x + end do + l = l + 1 + local_box_x(l) = ii + local_box_y(l) = jj + + + end do + end do + num_local_box_points = l + + end subroutine get_box_points + + + !*********************************************************************** ! ! subroutine get_box @@ -523,30 +835,17 @@ end function get_dom_landmask !> this subroutine and stored in the module variable 'box_mean'. ! !----------------------------------------------------------------------- - subroutine get_box(lat, lon, dx) + subroutine get_box(lat, lon, nx, ny) implicit none - real (kind=RKIND), intent(in) :: lat, lon, dx + real (kind=RKIND), intent(in) :: lat, lon + integer, intent(in) :: nx, ny integer :: i, j, ii, jj, ic, jc real (kind=RKIND) :: sg_lat - - ! - ! Get number of points to extract in the zonal direction - ! - if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) - else - nx = topo_x / 2 - end if - - ! - ! Get number of points to extract in the meridional direction - ! - ny = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re)) - + !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/lon,lat/), intArgs=(/nx,ny/)) ! ! Find coordinates in global topography array of the box center ! @@ -610,6 +909,73 @@ subroutine get_box(lat, lon, dx) end subroutine get_box + subroutine get_box_new(lat, lon, nx, ny) + + implicit none + + real (kind=RKIND), intent(in) :: lat, lon + integer, intent(in) :: nx, ny + integer :: itile, ix_shift, ix, jx, ibox, tile_index_i, tile_index_j + integer :: tile_base, tile_offset, i, j + + + if (associated(box_new)) deallocate(box_new) + allocate(box_new(nx, ny)) + + + ! if (associated(box_landuse)) deallocate(box_landuse) + ! allocate(box_landuse(nx,ny)) + + ! if (associated(dxm)) deallocate(dxm) + ! allocate(dxm(nx,ny)) + + ! + ! Extract sub-array (box) from global array; must properly account for + ! the periodicity in the longitude coordinate, as well as the poles + ! + box_mean_new = 0.0 + ibox = 1 + do j=1,ny + do i=1,nx + + !j = 25 + !i = 16 + ibox = (j-1)*nx + i + ! pixels on the global topo grid + ix = local_box_x(ibox) + jx = local_box_y(ibox) + + + tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 + tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 + + + + itile = locate_tile_in_set(tile_index_i, tile_index_j) + + !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) + + !ix_shift = (itile-1)*tile_x*tile_y + 1 + tile_base = (itile-1)*tile_x*tile_y + 1 + !tile_base = local_tile_x(tile_index_i) + local_tile_y(tile_index_j)* + + tile_offset = (ix - local_tile_x(itile)) + (jx - local_tile_y(itile)) * tile_x + !tile_offset = mod(ix, local_tile_x(tile_index_i)) + mod(jx, local_tile_y(tile_index_j)) * tile_x + + box_new(i,j) = topo_new(tile_base + tile_offset) + box_mean_new = box_mean_new + box_new(i,j) + !call mpas_log_write('In box_new for (i,j):($i,$i), (tile_base,tile_offset):($i,$i)', intArgs=(/i,j,tile_base,tile_offset/)) + if (box_new(i,j) .ne. box(i,j)) call mpas_log_write('In box_new for (i,j):($i,$i), box_new(i,j)= $r box(i,j)= $r ', intArgs=(/i,j/), realArgs=(/box_new(i,j),box(i,j)/)) + ibox = ibox + 1 + end do + end do + ! + ! Compute mean topography in the extracted box + ! + box_mean_new = box_mean_new / real(nx*ny, RKIND) + + end subroutine get_box_new + + !*********************************************************************** ! ! function get_var @@ -620,11 +986,12 @@ end subroutine get_box !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_var() + real (kind=RKIND) function get_var(nx, ny) implicit none integer :: i, j + integer, intent(in) :: nx, ny real (kind=RKIND) :: s2 s2 = 0.0 @@ -640,6 +1007,27 @@ real (kind=RKIND) function get_var() end function get_var + real (kind=RKIND) function get_var_new(nx, ny) + + implicit none + + integer :: i, j + integer, intent(in) :: nx, ny + real (kind=RKIND) :: s2 + + s2 = 0.0 + + do j=1,ny + do i=1,nx + s2 = s2 + (box_new(i,j) - box_mean_new)**2 + end do + end do + + get_var_new = sqrt(s2 / real(nx*ny,RKIND)) + + end function get_var_new + + !*********************************************************************** ! ! function get_con @@ -650,11 +1038,12 @@ end function get_var !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_con() + real (kind=RKIND) function get_con(nx, ny) implicit none integer :: i, j + integer, intent(in) :: nx, ny real (kind=RKIND) :: s2, s4, var, xland, mean_land, mean_water, oro s2 = 0.0 @@ -715,6 +1104,72 @@ real (kind=RKIND) function get_con() end function get_con + real (kind=RKIND) function get_con_new(nx, ny) + + implicit none + + integer :: i, j + integer, intent(in) :: nx, ny + real (kind=RKIND) :: s2, s4, var, xland, mean_land, mean_water, oro + + s2 = 0.0 + s4 = 0.0 + mean_land = 0.0 + mean_water = 0.0 + xland = 0.0 + + ! + ! Compute grid-box mean + ! + do j=1,ny + do i=1,nx + if (box_landuse(i,j) /= WATER) then + xland = xland + 1.0 + mean_land = mean_land + box_new(i,j) + else + mean_water = mean_water + box_new(i,j) + end if + end do + end do + if (xland > 0.0) then + mean_land = mean_land / xland + end if + if (xland < real(nx*ny,kind=RKIND)) then + mean_water = mean_water / (real(nx*ny,kind=RKIND) - xland) + end if + xland = xland / real(nx*ny,kind=RKIND) + + if (xland >= 0.5_RKIND) then + oro = mean_land + else + oro = mean_water + end if + + do j=1,ny + do i=1,nx + s2 = s2 + (box_new(i,j) - box_mean_new)**2 + s4 = s4 + (box_new(i,j) - oro)**4 + end do + end do + + var = s2 / real(nx*ny,RKIND) + + if (sqrt(var) < 1.0) then + get_con_new = 0.0 + else + get_con_new = s4 / (var**2 * real(nx*ny,RKIND)) + end if + + ! + ! Zero-ing all convexity statistics over dominantly water points. + ! + if (xland < 0.5_RKIND) then + get_con_new = 0.0 + end if + + end function get_con_new + + !*********************************************************************** ! ! function get_oa1 @@ -727,11 +1182,12 @@ end function get_con !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa1() + real (kind=RKIND) function get_oa1(nx, ny) implicit none integer :: i, j + integer, intent(in) :: nx, ny integer :: nu, nd nu = 0 @@ -766,9 +1222,10 @@ end function get_oa1 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa2() + real (kind=RKIND) function get_oa2(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nu, nd @@ -807,9 +1264,10 @@ end function get_oa2 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa3() + real (kind=RKIND) function get_oa3(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nu, nd @@ -849,9 +1307,10 @@ end function get_oa3 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa4() + real (kind=RKIND) function get_oa4(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nu, nd @@ -889,9 +1348,10 @@ end function get_oa4 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol1() + real (kind=RKIND) function get_ol1(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nw @@ -922,9 +1382,10 @@ end function get_ol1 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol2() + real (kind=RKIND) function get_ol2(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nw @@ -955,9 +1416,10 @@ end function get_ol2 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol3() + real (kind=RKIND) function get_ol3(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nw @@ -994,9 +1456,10 @@ end function get_ol3 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol4() + real (kind=RKIND) function get_ol4(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j integer :: nw @@ -1033,9 +1496,10 @@ end function get_ol4 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_elvmax() + real (kind=RKIND) function get_elvmax(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j @@ -1062,9 +1526,10 @@ end function get_elvmax !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_htheta() + real (kind=RKIND) function get_htheta(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j real (kind=RKIND) :: dx, dy @@ -1110,9 +1575,10 @@ end function get_htheta !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_hgamma() + real (kind=RKIND) function get_hgamma(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j real (kind=RKIND) :: dx, dy @@ -1163,9 +1629,10 @@ end function get_hgamma !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_hsigma() + real (kind=RKIND) function get_hsigma(nx, ny) implicit none + integer, intent(in) :: nx, ny integer :: i, j real (kind=RKIND) :: dx, dy From a5d69c6956af519a300e0ddc0246df1cf28b3f9f Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 27 Aug 2024 20:50:27 +0000 Subject: [PATCH 02/27] all statistics seem correct --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 409 +++++++++++++++++-- 1 file changed, 379 insertions(+), 30 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 5889dabbb9..3018a82a25 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -53,6 +53,8 @@ end subroutine read_geogrid real (kind=RKIND) :: start_lat real (kind=RKIND) :: start_lon + integer :: topo_shift ! special handling + ! Nominal delta-x (in meters) for sub-grid topography cells real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) @@ -64,7 +66,9 @@ end subroutine read_geogrid real (kind=RKIND) :: box_mean, box_mean_new ! Mean value of topography in box !integer :: nx, ny ! Dimensions of box covering grid cell integer (kind=I1KIND), dimension(:,:), pointer :: landuse ! Global 30-arc-second landuse + integer (kind=I1KIND), dimension(:), pointer :: landuse_new ! Global 30-arc-second landuse integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell + integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse_new ! Subset of landuse covering a grid cell ! NB: At present, only the USGS GLCC land cover dataset is supported, so we can assume 16 == water ! See the read_global_30s_landuse function @@ -73,7 +77,7 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height - integer, dimension(:), pointer :: local_tile_x ! + integer, dimension(:), pointer :: local_tile_x, local_tile_x_land ! integer, dimension(:), pointer :: local_tile_y ! integer, dimension(:), pointer :: local_box_x ! integer, dimension(:), pointer :: local_box_y ! @@ -127,13 +131,13 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND), pointer :: config_topo_data character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - real (kind=RKIND) :: var2d_tmp, con_tmp!, oa1, oa2, oa3, oa4, ol1, ol2, ol3, ol4 + real (kind=RKIND) :: var2d_tmp, con_tmp, oa1_tmp, oa2_tmp, oa3_tmp, oa4_tmp, ol1_tmp, ol2_tmp, ol3_tmp, ol4_tmp, hlanduse_tmp ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse real (kind=RKIND) :: sum_var integer :: nx, ny - + character(len=StrKIND):: message allocate(topo(topo_x,topo_y)) allocate(landuse(topo_x,topo_y)) @@ -228,11 +232,7 @@ function compute_gwd_fields(domain) result(iErr) - iErr = read_global_30s_landuse(geog_data_path) - if (iErr /= 0) then - call mpas_log_write('Error reading global 30-arc-sec landuse for GWD statistics', messageType=MPAS_LOG_ERR) - return - end if + ! ! It is possible that this code is called before the mesh fields have been scaled @@ -258,15 +258,17 @@ function compute_gwd_fields(domain) result(iErr) max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) call mpas_log_write(' hi from compute_gwd_fields -- before allocate max_local_tiles:$i', intArgs=(/max_local_tiles/)) allocate(local_tile_x(max_local_tiles)) + allocate(local_tile_x_land(max_local_tiles)) allocate(local_tile_y(max_local_tiles)) + topo_shift = topo_x / 2 ! ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. ! cum_local_tiles=0 do iCell=1,nCells - !do iCell=5,65 + !do iCell=5,25 ! ! First, get an estimate of the mean diameter (in meters) of the grid @@ -301,8 +303,14 @@ function compute_gwd_fields(domain) result(iErr) return end if + iErr = read_global_30s_landuse(geog_data_path) + if (iErr /= 0) then + call mpas_log_write('Error reading global 30-arc-sec landuse for GWD statistics', messageType=MPAS_LOG_ERR) + return + end if + do iCell=1,nCells - !do iCell=5,65 + !do iCell=5,25 ! ! First, get an estimate of the mean diameter (in meters) of the grid @@ -346,14 +354,22 @@ function compute_gwd_fields(domain) result(iErr) ! var2d(iCell) = get_var(nx, ny) var2d_tmp = get_var_new(nx, ny) - call mpas_log_write('iCell = $i .... var2d: $r var2d_tmp: $r', intArgs=(/iCell/), realArgs=(/var2d(iCell), var2d_tmp/)) + if ( var2d(iCell) .ne. var2d_tmp) call mpas_log_write('iCell = $i .... var2d: $r var2d_tmp: $r', intArgs=(/iCell/), realArgs=(/var2d(iCell), var2d_tmp/)) con(iCell) = get_con(nx, ny) con_tmp = get_con_new(nx, ny) - call mpas_log_write('iCell = $i .... con: $r con_tmp: $r', intArgs=(/iCell/), realArgs=(/con(iCell), con_tmp/)) + if ( con(iCell) .ne. con_tmp) call mpas_log_write('iCell = $i .... con: $r con_tmp: $r', intArgs=(/iCell/), realArgs=(/con(iCell), con_tmp/)) oa1(iCell) = get_oa1(nx, ny) + oa1_tmp = get_oa1_new(nx, ny) + if ( oa1(iCell) .ne. oa1_tmp) call mpas_log_write('iCell = $i .... oa1: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa1(iCell), oa1_tmp/)) oa2(iCell) = get_oa2(nx, ny) + oa2_tmp = get_oa2_new(nx, ny) + if ( oa2(iCell) .ne. oa2_tmp) call mpas_log_write('iCell = $i .... oa2: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa2(iCell), oa2_tmp/)) oa3(iCell) = get_oa3(nx, ny) + oa3_tmp = get_oa3_new(nx, ny) + if ( oa3(iCell) .ne. oa3_tmp) call mpas_log_write('iCell = $i .... oa3: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa3(iCell), oa4_tmp/)) oa4(iCell) = get_oa4(nx, ny) + oa4_tmp = get_oa4_new(nx, ny) + if ( oa4(iCell) .ne. oa4_tmp) call mpas_log_write('iCell = $i .... oa4: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa4(iCell), oa4_tmp/)) ! Critical height, to be used in OL computation ! See Appendix of Kim, Y-J, 1996: Representation of Sub-Grid Scale Orographic Effects @@ -361,11 +377,24 @@ function compute_gwd_fields(domain) result(iErr) hc = 1116.2_RKIND - 0.878_RKIND * var2d(iCell) ol1(iCell) = get_ol1(nx, ny) + ol1_tmp = get_ol1_new(nx, ny) + if ( ol1(iCell) .ne. ol1_tmp) call mpas_log_write('iCell = $i .... ol1: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol1(iCell), ol1_tmp/)) ol2(iCell) = get_ol2(nx, ny) + ol2_tmp = get_ol2_new(nx, ny) + if ( ol2(iCell) .ne. ol2_tmp) call mpas_log_write('iCell = $i .... ol2: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol2(iCell), ol2_tmp/)) ol3(iCell) = get_ol3(nx, ny) + ol3_tmp = get_ol3_new(nx, ny) + if ( ol3(iCell) .ne. ol3_tmp) call mpas_log_write('iCell = $i .... ol3: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol3(iCell), ol3_tmp/)) ol4(iCell) = get_ol4(nx, ny) + ol4_tmp = get_ol4_new(nx, ny) + if ( ol4(iCell) .ne. ol4_tmp) call mpas_log_write('iCell = $i .... ol4: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol4(iCell), ol4_tmp/)) hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell + hlanduse_tmp = get_dom_landmask_new(nx, ny) ! get dominant land mask in cell + if ( hlanduse(iCell) .ne. hlanduse_tmp) then + write(message,fmt='(A,i10, A, i18, A, i18)') 'iCell = ',iCell,' .... hlanduse: ',hlanduse(iCell),' _tmp:',hlanduse_tmp + call mpas_log_write(message) + end if ! elvmax(iCell) = get_elvmax() ! htheta(iCell) = get_htheta() @@ -454,14 +483,19 @@ function unique_tile(ix, jx, len) result (is_unique) end function unique_tile - function locate_tile_in_set(tile_index_i, tile_index_j) result (tile_index) + function locate_tile_in_set(tile_index_i, tile_index_j, is_topo) result (tile_index) implicit none integer, intent(in) :: tile_index_i integer, intent(in) :: tile_index_j + logical, intent(in) :: is_topo - integer :: itile, tile_index + integer :: itile, tile_index, tile_index_i_local + + tile_index_i_local = tile_index_i + !if (is_topo) tile_index_i_local = tile_index_i + topo_shift + do itile = 1, cum_local_tiles if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then @@ -513,7 +547,12 @@ function determine_complete_tile_set(nx, ny) result(iErr) !jx = nint((latCell(iCell)*rad2deg - start_lat) * pts_per_degree) ii = cum_local_tiles do l=1,nx*ny - itile = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 + + if (local_box_x(l) > topo_shift) then + itile = floor( real(local_box_x(l) - topo_shift - 1) / real(tile_x)) * tile_x + 1 + else + itile = floor( real(local_box_x(l) + topo_shift - 1) / real(tile_x)) * tile_x + 1 + end if jtile = floor( real(local_box_y(l) - 1) / real(tile_y)) * tile_y + 1 !itile = (ix / tile_x) * tile_x + 1 !jtile = (jx / tile_y) * tile_y + 1 @@ -522,9 +561,10 @@ function determine_complete_tile_set(nx, ny) result(iErr) if(unique_tile(itile, jtile, ii)) then ii = ii + 1 local_tile_x(ii) = itile + local_tile_x_land(ii) = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 local_tile_y(ii) = jtile - !call mpas_log_write('---- Adding tile to local tile list ($i, $i)', intArgs=(/local_tile_x(ii),local_tile_y(ii)/)) - !call mpas_log_write('---- local box ($i, $i)', intArgs=(/local_box_x(l),local_box_y(l)/)) + call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) + !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) end if end do @@ -563,7 +603,7 @@ function read_global_30s_topo(path, sub_path) result(iErr) integer, parameter :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile integer (c_int) :: istatus - integer :: ix, iy, ishift, ix_shift, il + integer :: ix, iy, ix_shift, il integer(kind=I8KIND):: ix_shift_long integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor @@ -584,9 +624,10 @@ function read_global_30s_topo(path, sub_path) result(iErr) ny = tile_y + 2*tile_bdr nz = 1 - ishift = 0 + topo_shift = 0 allocate(topo_new(tile_x*tile_y*cum_local_tiles)) + allocate(landuse_new(tile_x*tile_y*cum_local_tiles)) call mpas_log_write('---- topo_x:$i topo_y: $i', intArgs=(/topo_x,topo_y/)) call mpas_log_write('---- tile_x:$i tile_y: $i', intArgs=(/tile_x,tile_y/)) ! @@ -594,12 +635,16 @@ function read_global_30s_topo(path, sub_path) result(iErr) ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 ! if (trim(sub_path) == 'topo_gmted2010_30s/') then - ishift = topo_x / 2 + topo_shift = topo_x / 2 end if do il=1,cum_local_tiles ix = local_tile_x(il) + + !if (ix<=topo_shift) ix = ix + topo_shift + iy = local_tile_y(il) + call mpas_log_write(' --- Try to read tile with ix,iy ($i,$i) ', intArgs=(/ix, iy/)) 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) call mpas_f_to_c_string(filename, c_filename) @@ -614,8 +659,8 @@ function read_global_30s_topo(path, sub_path) result(iErr) !ix_shift = mod((ix-1) + ishift, topo_x) + 1 ix_shift = (il-1)*tile_x*tile_y + 1 - call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,ishift/)) - !call mpas_log_write(' topo_x,y ($i,$i) ', intArgs=(/topo_x, topo_y/)) + call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) + call mpas_log_write(' --- with ix,iy ($i,$i) ', intArgs=(/ix, iy/)) !write(message,fmt='(A,i18,i18,A,i18,A,i18)') 'Reading tile (',il,cum_local_tiles,') topo_new start=', ix_shift_long,' end=',ix_shift_long+tile_x*tile_y-1 !call mpas_log_write(message) topo_new(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) @@ -636,7 +681,7 @@ function read_global_30s_topo(path, sub_path) result(iErr) return end if - ix_shift = mod((ix-1) + ishift, topo_x) + 1 + ix_shift = mod((ix-1) + topo_shift, topo_x) + 1 topo(ix_shift:(ix_shift+tile_x-1),iy:(iy+tile_y-1)) = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) end do @@ -673,13 +718,14 @@ function read_global_30s_landuse(path) result(iErr) integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second landuse integer (c_int) :: istatus - integer :: ix, iy + integer :: ix, iy, ix_shift, il integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile type (c_ptr) :: tile_ptr character(len=StrKIND) :: filename character(kind=c_char), dimension(StrKIND+1) :: c_filename + character(len=StrKIND):: message allocate(tile(tile_x,tile_y,1)) tile_ptr = c_loc(tile) @@ -692,6 +738,36 @@ function read_global_30s_landuse(path) result(iErr) ny = tile_y nz = 1 + + do il=1,cum_local_tiles + ix = local_tile_x_land(il) + iy = local_tile_y(il) + + write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//'/landuse_30s/', ix, '-', (ix+tile_x-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) + tile(:,:,:) = tile(:,:,:) * scalefactor + if (istatus /= 0) then + call mpas_log_write('Error reading landuse tile '//trim(filename)) + iErr = 1 + return + end if + + !ix_shift = mod((ix-1) + ishift, topo_x) + 1 + ix_shift = (il-1)*tile_x*tile_y + 1 + call mpas_log_write(' Reading tile ($i/$i) into landuse_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) + !call mpas_log_write(' topo_x,y ($i,$i) ', intArgs=(/topo_x, topo_y/)) + !write(message,fmt='(A,i18,i18,A,i18,A,i18)') 'Reading tile (',il,cum_local_tiles,') topo_new start=', ix_shift_long,' end=',ix_shift_long+tile_x*tile_y-1 + !call mpas_log_write(message) + landuse_new(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) + !write(message,fmt='(A,i18)') 'At 25,16 for landuse_new for (i,j):',landuse_new(ix_shift+25+16*tile_x) + !call mpas_log_write(message) + end do + + + do iy=1,topo_y,tile_y do ix=1,topo_x,tile_x write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//'/landuse_30s/', ix, '-', (ix+tile_x-1), '.', & @@ -756,6 +832,34 @@ integer (kind=I1KIND) function get_dom_landmask(nx, ny) end function get_dom_landmask + integer (kind=I1KIND) function get_dom_landmask_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + real (kind=RKIND) :: xland + xland = 0.0_RKIND + + ! Get dominant land/water mask in the box + do j=1,ny + do i=1,nx + if (box_landuse_new(i,j) /= WATER) then + xland = xland + 1.0_RKIND + end if + end do + end do + xland = xland / real(nx*ny,kind=RKIND) + + if (xland >= 0.5_RKIND) then + get_dom_landmask_new = 1_I1KIND + else + get_dom_landmask_new = 0_I1KIND + end if + + end function get_dom_landmask_new + + subroutine get_box_points(lat, lon, nx, ny) implicit none @@ -771,7 +875,7 @@ subroutine get_box_points(lat, lon, nx, ny) ! ! Find coordinates in global topography array of the box center ! - ic = nint((lon) * pts_per_degree) + 1 + ic = nint((lon - start_lon) * pts_per_degree) + 1 jc = nint((lat - start_lat) * pts_per_degree) + 1 if (ic <= 0) ic = ic + topo_x @@ -805,6 +909,9 @@ subroutine get_box_points(lat, lon, nx, ny) l = l + 1 local_box_x(l) = ii local_box_y(l) = jj + + !if (i == 25 .and. j ==16) call mpas_log_write('In get box for 25,16 (i,j):($i,$i)', intArgs=(/ii,jj/)) + end do @@ -915,13 +1022,18 @@ subroutine get_box_new(lat, lon, nx, ny) real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny - integer :: itile, ix_shift, ix, jx, ibox, tile_index_i, tile_index_j + integer :: itile, ix_shift, ix, ix_land, ix_topo, jx, ibox, tile_index_i, tile_index_j integer :: tile_base, tile_offset, i, j + real (kind=RKIND) :: sg_lat + character(len=StrKIND):: message if (associated(box_new)) deallocate(box_new) allocate(box_new(nx, ny)) + if (associated(box_landuse_new)) deallocate(box_landuse_new) + allocate(box_landuse_new(nx, ny)) + ! if (associated(box_landuse)) deallocate(box_landuse) ! allocate(box_landuse(nx,ny)) @@ -941,17 +1053,23 @@ subroutine get_box_new(lat, lon, nx, ny) !j = 25 !i = 16 ibox = (j-1)*nx + i - ! pixels on the global topo grid - ix = local_box_x(ibox) + ! pixels on the global topo grid + ix = local_box_x(ibox) !+ topo_shift + ix_topo = local_box_x(ibox) jx = local_box_y(ibox) + if (ix - topo_shift > 0 ) then + ix = ix - topo_shift + else + ix = ix + topo_shift + end if tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 - itile = locate_tile_in_set(tile_index_i, tile_index_j) + itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) @@ -963,8 +1081,19 @@ subroutine get_box_new(lat, lon, nx, ny) box_new(i,j) = topo_new(tile_base + tile_offset) box_mean_new = box_mean_new + box_new(i,j) + + box_landuse_new(i,j) = landuse_new(tile_base + tile_offset) + sg_lat = (start_lat + (real(jx-1,RKIND) + 0.5) / pts_per_degree) / rad2deg ! Add 0.5 for cell center + dxm(i,j) = sg_delta * cos(sg_lat) + !call mpas_log_write('In box_new for (i,j):($i,$i), (tile_base,tile_offset):($i,$i)', intArgs=(/i,j,tile_base,tile_offset/)) if (box_new(i,j) .ne. box(i,j)) call mpas_log_write('In box_new for (i,j):($i,$i), box_new(i,j)= $r box(i,j)= $r ', intArgs=(/i,j/), realArgs=(/box_new(i,j),box(i,j)/)) + + if (box_landuse_new(i,j) .ne. landuse(ix_topo,jx)) then + write(message,fmt='(A,i10,A,i10,A,i18,A,i18)') 'In box_landuse_new for (i,j):(',i,',',j,'), box_landuse_new(i,j)=', box_landuse_new(i,j),' box_landuse(i,j)=',box_landuse(i,j) + call mpas_log_write(message) + !call mpas_log_write('In box_landuse_new for (i,j):($i,$i), box_landuse_new(i,j)= $i box_landuse(i,j)= $i ', intArgs=(/i,j,box_landuse_new(i,j),box_landuse(i,j)/)) + end if ibox = ibox + 1 end do end do @@ -1123,7 +1252,7 @@ real (kind=RKIND) function get_con_new(nx, ny) ! do j=1,ny do i=1,nx - if (box_landuse(i,j) /= WATER) then + if (box_landuse_new(i,j) /= WATER) then xland = xland + 1.0 mean_land = mean_land + box_new(i,j) else @@ -1209,6 +1338,32 @@ real (kind=RKIND) function get_oa1(nx, ny) end function get_oa1 + real (kind=RKIND) function get_oa1_new(nx, ny) + + implicit none + + integer :: i, j + integer, intent(in) :: nx, ny + integer :: nu, nd + + nu = 0 + nd = 0 + do j=1,ny + do i=1,nx/2 + if (box_new(i,j) > box_mean_new) nu = nu + 1 + end do + do i=nx/2+1,nx + if (box_new(i,j) > box_mean_new) nd = nd + 1 + end do + end do + + if (nu + nd > 0) then + get_oa1_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) + else + get_oa1_new = 0.0 + end if + + end function get_oa1_new !*********************************************************************** ! @@ -1251,6 +1406,34 @@ real (kind=RKIND) function get_oa2(nx, ny) end function get_oa2 + real (kind=RKIND) function get_oa2_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nu, nd + + nu = 0 + nd = 0 + do j=1,ny/2 + do i=1,nx + if (box_new(i,j) > box_mean_new) nu = nu + 1 + end do + end do + do j=ny/2+1,ny + do i=1,nx + if (box_new(i,j) > box_mean_new) nd = nd + 1 + end do + end do + + if (nu + nd > 0) then + get_oa2_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) + else + get_oa2_new = 0.0 + end if + + end function get_oa2_new !*********************************************************************** ! @@ -1295,6 +1478,36 @@ real (kind=RKIND) function get_oa3(nx, ny) end function get_oa3 + real (kind=RKIND) function get_oa3_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nu, nd + real (kind=RKIND) :: ratio + + nu = 0 + nd = 0 + ratio = real(ny,RKIND)/real(nx,RKIND) + do j=1,ny + do i=1,nx + if (nint(real(i,RKIND) * ratio) < (ny - j)) then + if (box_new(i,j) > box_mean_new) nu = nu + 1 + else + if (box_new(i,j) > box_mean_new) nd = nd + 1 + end if + end do + end do + + if (nu + nd > 0) then + get_oa3_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) + else + get_oa3_new = 0.0 + end if + + end function get_oa3_new + !*********************************************************************** ! ! function get_oa4 @@ -1338,6 +1551,35 @@ real (kind=RKIND) function get_oa4(nx, ny) end function get_oa4 + real (kind=RKIND) function get_oa4_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nu, nd + real (kind=RKIND) :: ratio + + nu = 0 + nd = 0 + ratio = real(ny,RKIND)/real(nx,RKIND) + do j=1,ny + do i=1,nx + if (nint(real(i,RKIND) * ratio) < j) then + if (box_new(i,j) > box_mean_new) nu = nu + 1 + else + if (box_new(i,j) > box_mean_new) nd = nd + 1 + end if + end do + end do + + if (nu + nd > 0) then + get_oa4_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) + else + get_oa4_new = 0.0 + end if + +end function get_oa4_new !*********************************************************************** ! ! function get_ol1 @@ -1372,6 +1614,30 @@ real (kind=RKIND) function get_ol1(nx, ny) end function get_ol1 + real (kind=RKIND) function get_ol1_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nw + integer :: nt + + nw = 0 + nt = 0 + + do j=ny/4,3*ny/4 + do i=1,nx + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + + get_ol1_new = real(nw,RKIND) / real(nt,RKIND) + + end function get_ol1_new + + !*********************************************************************** ! ! function get_ol2 @@ -1406,6 +1672,31 @@ real (kind=RKIND) function get_ol2(nx, ny) end function get_ol2 + + real (kind=RKIND) function get_ol2_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nw + integer :: nt + + nw = 0 + nt = 0 + + do j=1,ny + do i=nx/4,3*nx/4 + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + + get_ol2_new = real(nw,RKIND) / real(nt,RKIND) + + end function get_ol2_new + + !*********************************************************************** ! ! function get_ol3 @@ -1445,6 +1736,35 @@ real (kind=RKIND) function get_ol3(nx, ny) end function get_ol3 + real (kind=RKIND) function get_ol3_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nw + integer :: nt + + nw = 0 + nt = 0 + + do j=1,ny/2 + do i=1,nx/2 + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + do j=ny/2+1,ny + do i=nx/2+1,nx + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + + get_ol3_new = real(nw,RKIND) / real(nt,RKIND) + + end function get_ol3_new + !*********************************************************************** ! @@ -1486,6 +1806,35 @@ real (kind=RKIND) function get_ol4(nx, ny) end function get_ol4 + real (kind=RKIND) function get_ol4_new(nx, ny) + + implicit none + integer, intent(in) :: nx, ny + + integer :: i, j + integer :: nw + integer :: nt + + nw = 0 + nt = 0 + + do j=ny/2+1,ny + do i=1,nx/2 + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + do j=1,ny/2 + do i=nx/2+1,nx + if (box_new(i,j) > hc) nw = nw + 1 + nt = nt + 1 + end do + end do + + get_ol4_new = real(nw,RKIND) / real(nt,RKIND) + + end function get_ol4_new + !*********************************************************************** ! ! function get_elvmax From 882d6b8a149e403043f155b209d852ae8ad16b07 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 27 Aug 2024 22:58:07 +0000 Subject: [PATCH 03/27] cleanup --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 696 +++---------------- 1 file changed, 80 insertions(+), 616 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 3018a82a25..0244b7a292 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -58,17 +58,13 @@ end subroutine read_geogrid ! Nominal delta-x (in meters) for sub-grid topography cells real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) - real (kind=R4KIND), dimension(:,:), pointer :: topo ! Global 30-arc-second topography - real (kind=R4KIND), dimension(:), pointer :: topo_new ! Global 30-arc-second topography + real (kind=R4KIND), dimension(:), pointer :: topo ! Global 30-arc-second topography real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell - real (kind=RKIND), dimension(:,:), pointer :: box_new ! Subset of topography covering a grid cell real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell - real (kind=RKIND) :: box_mean, box_mean_new ! Mean value of topography in box + real (kind=RKIND) :: box_mean ! Mean value of topography in box !integer :: nx, ny ! Dimensions of box covering grid cell - integer (kind=I1KIND), dimension(:,:), pointer :: landuse ! Global 30-arc-second landuse - integer (kind=I1KIND), dimension(:), pointer :: landuse_new ! Global 30-arc-second landuse + integer (kind=I1KIND), dimension(:), pointer :: landuse ! Global 30-arc-second landuse integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell - integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse_new ! Subset of landuse covering a grid cell ! NB: At present, only the USGS GLCC land cover dataset is supported, so we can assume 16 == water ! See the read_global_30s_landuse function @@ -81,7 +77,7 @@ end subroutine read_geogrid integer, dimension(:), pointer :: local_tile_y ! integer, dimension(:), pointer :: local_box_x ! integer, dimension(:), pointer :: local_box_y ! - integer:: max_local_tiles, cum_local_tiles, num_local_box_points + integer:: max_local_tiles, cum_local_tiles contains @@ -132,6 +128,8 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash real (kind=RKIND) :: var2d_tmp, con_tmp, oa1_tmp, oa2_tmp, oa3_tmp, oa4_tmp, ol1_tmp, ol2_tmp, ol3_tmp, ol4_tmp, hlanduse_tmp + real (kind=RKIND) :: elvmax_1, htheta_1, hgamma_1, hsigma_1, elvmax_tmp, htheta_tmp, hgamma_tmp, hsigma_tmp + ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse @@ -139,8 +137,6 @@ function compute_gwd_fields(domain) result(iErr) integer :: nx, ny character(len=StrKIND):: message - allocate(topo(topo_x,topo_y)) - allocate(landuse(topo_x,topo_y)) call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh) call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state) @@ -286,16 +282,15 @@ function compute_gwd_fields(domain) result(iErr) ! call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) - - !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg/), intArgs=(/nx,ny/)) call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) - - !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius, dc) + iErr = determine_complete_tile_set(nx, ny) if (cum_local_tiles >= max_local_tiles) exit - end do + end do + + call mpas_log_write(' End of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/cum_local_tiles, max_local_tiles/)) iErr = read_global_30s_topo(geog_data_path, geog_sub_path) if (iErr /= 0) then @@ -326,13 +321,7 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - - !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg/), intArgs=(/nx,ny/)) - - !call get_box_points(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc) call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) - !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius, dc) - !iErr = determine_complete_tile_set() ! ! Cut out a rectangular piece of the global 30-arc-second topography ! data that is centered at the lat/lon of the current cell being @@ -343,33 +332,19 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) - call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) - call get_box_new(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) + call get_box(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) ! ! With a box of 30-arc-second data for the current grid cell, call ! subroutines to compute each sub-grid orography statistic ! var2d(iCell) = get_var(nx, ny) - var2d_tmp = get_var_new(nx, ny) - if ( var2d(iCell) .ne. var2d_tmp) call mpas_log_write('iCell = $i .... var2d: $r var2d_tmp: $r', intArgs=(/iCell/), realArgs=(/var2d(iCell), var2d_tmp/)) con(iCell) = get_con(nx, ny) - con_tmp = get_con_new(nx, ny) - if ( con(iCell) .ne. con_tmp) call mpas_log_write('iCell = $i .... con: $r con_tmp: $r', intArgs=(/iCell/), realArgs=(/con(iCell), con_tmp/)) oa1(iCell) = get_oa1(nx, ny) - oa1_tmp = get_oa1_new(nx, ny) - if ( oa1(iCell) .ne. oa1_tmp) call mpas_log_write('iCell = $i .... oa1: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa1(iCell), oa1_tmp/)) oa2(iCell) = get_oa2(nx, ny) - oa2_tmp = get_oa2_new(nx, ny) - if ( oa2(iCell) .ne. oa2_tmp) call mpas_log_write('iCell = $i .... oa2: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa2(iCell), oa2_tmp/)) oa3(iCell) = get_oa3(nx, ny) - oa3_tmp = get_oa3_new(nx, ny) - if ( oa3(iCell) .ne. oa3_tmp) call mpas_log_write('iCell = $i .... oa3: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa3(iCell), oa4_tmp/)) oa4(iCell) = get_oa4(nx, ny) - oa4_tmp = get_oa4_new(nx, ny) - if ( oa4(iCell) .ne. oa4_tmp) call mpas_log_write('iCell = $i .... oa4: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/oa4(iCell), oa4_tmp/)) ! Critical height, to be used in OL computation ! See Appendix of Kim, Y-J, 1996: Representation of Sub-Grid Scale Orographic Effects @@ -377,29 +352,17 @@ function compute_gwd_fields(domain) result(iErr) hc = 1116.2_RKIND - 0.878_RKIND * var2d(iCell) ol1(iCell) = get_ol1(nx, ny) - ol1_tmp = get_ol1_new(nx, ny) - if ( ol1(iCell) .ne. ol1_tmp) call mpas_log_write('iCell = $i .... ol1: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol1(iCell), ol1_tmp/)) ol2(iCell) = get_ol2(nx, ny) - ol2_tmp = get_ol2_new(nx, ny) - if ( ol2(iCell) .ne. ol2_tmp) call mpas_log_write('iCell = $i .... ol2: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol2(iCell), ol2_tmp/)) ol3(iCell) = get_ol3(nx, ny) - ol3_tmp = get_ol3_new(nx, ny) - if ( ol3(iCell) .ne. ol3_tmp) call mpas_log_write('iCell = $i .... ol3: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol3(iCell), ol3_tmp/)) ol4(iCell) = get_ol4(nx, ny) - ol4_tmp = get_ol4_new(nx, ny) - if ( ol4(iCell) .ne. ol4_tmp) call mpas_log_write('iCell = $i .... ol4: $r _tmp: $r', intArgs=(/iCell/), realArgs=(/ol4(iCell), ol4_tmp/)) hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell - hlanduse_tmp = get_dom_landmask_new(nx, ny) ! get dominant land mask in cell - if ( hlanduse(iCell) .ne. hlanduse_tmp) then - write(message,fmt='(A,i10, A, i18, A, i18)') 'iCell = ',iCell,' .... hlanduse: ',hlanduse(iCell),' _tmp:',hlanduse_tmp - call mpas_log_write(message) - end if -! elvmax(iCell) = get_elvmax() -! htheta(iCell) = get_htheta() -! hgamma(iCell) = get_hgamma() -! hsigma(iCell) = get_hsigma() + ! elvmax_1 = get_elvmax(nx,ny) + ! htheta_1 = get_htheta(nx,ny) + ! hgamma_1 = get_hgamma(nx,ny) + ! hsigma_1 = get_hsigma(nx,ny) + end do @@ -491,11 +454,8 @@ function locate_tile_in_set(tile_index_i, tile_index_j, is_topo) result (tile_in integer, intent(in) :: tile_index_j logical, intent(in) :: is_topo - integer :: itile, tile_index, tile_index_i_local + integer :: itile, tile_index - tile_index_i_local = tile_index_i - !if (is_topo) tile_index_i_local = tile_index_i + topo_shift - do itile = 1, cum_local_tiles if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then @@ -511,9 +471,6 @@ end function locate_tile_in_set function determine_complete_tile_set(nx, ny) result(iErr) implicit none - !integer, intent(in) :: nCells - !REAL(kind=RKIND), DIMENSION(:), POINTER, intent(in) :: latCell - !REAL(kind=RKIND), DIMENSION(:), POINTER, intent(in) :: lonCell integer :: iErr integer, intent(in) :: nx, ny @@ -528,8 +485,7 @@ function determine_complete_tile_set(nx, ny) result(iErr) !allocate(local_tile_y(max_local_tiles)) - !call mpas_log_write(' hi from determine_complete_tile_set -- after allocate') - !do iCell=1,nCells + ! ! Cut out a rectangular piece of the global 30-arc-second topography @@ -541,38 +497,28 @@ function determine_complete_tile_set(nx, ny) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - !call get_box_points(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc) - - !ix = nint((lonCell(iCell)*rad2deg) * pts_per_degree) - !jx = nint((latCell(iCell)*rad2deg - start_lat) * pts_per_degree) - ii = cum_local_tiles - do l=1,nx*ny - - if (local_box_x(l) > topo_shift) then - itile = floor( real(local_box_x(l) - topo_shift - 1) / real(tile_x)) * tile_x + 1 - else - itile = floor( real(local_box_x(l) + topo_shift - 1) / real(tile_x)) * tile_x + 1 - end if - jtile = floor( real(local_box_y(l) - 1) / real(tile_y)) * tile_y + 1 - !itile = (ix / tile_x) * tile_x + 1 - !jtile = (jx / tile_y) * tile_y + 1 - !call mpas_log_write(' Checking if point ($i, $i) with lon lat ($i, $i) is already present', intArgs=(/ix, jx, nint(lonCell(iCell)*rad2deg), nint(latCell(iCell)*rad2deg)/)) - !call mpas_log_write(' Checking if local box point ($i, $i) in tile ($i, $i) is already present', intArgs=(/local_box_x(l), local_box_y(l), itile, jtile/)) - if(unique_tile(itile, jtile, ii)) then - ii = ii + 1 - local_tile_x(ii) = itile - local_tile_x_land(ii) = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 - local_tile_y(ii) = jtile - call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) - !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) - end if - end do - - cum_local_tiles = ii + ii = cum_local_tiles + do l=1,nx*ny + + if (local_box_x(l) > topo_shift) then + itile = floor( real(local_box_x(l) - topo_shift - 1) / real(tile_x)) * tile_x + 1 + else + itile = floor( real(local_box_x(l) + topo_shift - 1) / real(tile_x)) * tile_x + 1 + end if + jtile = floor( real(local_box_y(l) - 1) / real(tile_y)) * tile_y + 1 + + if(unique_tile(itile, jtile, ii)) then + ii = ii + 1 + local_tile_x(ii) = itile + local_tile_x_land(ii) = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 + local_tile_y(ii) = jtile + !call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) + !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) + end if + end do - !end do + cum_local_tiles = ii - call mpas_log_write(' hi from end of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/ii, max_local_tiles/)) ! end function determine_complete_tile_set @@ -626,8 +572,8 @@ function read_global_30s_topo(path, sub_path) result(iErr) topo_shift = 0 - allocate(topo_new(tile_x*tile_y*cum_local_tiles)) - allocate(landuse_new(tile_x*tile_y*cum_local_tiles)) + allocate(topo(tile_x*tile_y*cum_local_tiles)) + call mpas_log_write('---- topo_x:$i topo_y: $i', intArgs=(/topo_x,topo_y/)) call mpas_log_write('---- tile_x:$i tile_y: $i', intArgs=(/tile_x,tile_y/)) ! @@ -640,11 +586,8 @@ function read_global_30s_topo(path, sub_path) result(iErr) do il=1,cum_local_tiles ix = local_tile_x(il) - - !if (ix<=topo_shift) ix = ix + topo_shift - iy = local_tile_y(il) - call mpas_log_write(' --- Try to read tile with ix,iy ($i,$i) ', intArgs=(/ix, iy/)) + 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) call mpas_f_to_c_string(filename, c_filename) @@ -657,35 +600,13 @@ function read_global_30s_topo(path, sub_path) result(iErr) return end if - !ix_shift = mod((ix-1) + ishift, topo_x) + 1 ix_shift = (il-1)*tile_x*tile_y + 1 - call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) - call mpas_log_write(' --- with ix,iy ($i,$i) ', intArgs=(/ix, iy/)) - !write(message,fmt='(A,i18,i18,A,i18,A,i18)') 'Reading tile (',il,cum_local_tiles,') topo_new start=', ix_shift_long,' end=',ix_shift_long+tile_x*tile_y-1 - !call mpas_log_write(message) - topo_new(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + !call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) + + topo(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) end do - do iy=1,topo_y,tile_y - do ix=1,topo_x,tile_x - 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) - call mpas_f_to_c_string(filename, c_filename) - call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & - wordsize, istatus) - tile(:,:,:) = tile(:,:,:) * scalefactor - if (istatus /= 0) then - call mpas_log_write('Error reading topography tile '//trim(filename), messageType=MPAS_LOG_ERR) - iErr = 1 - return - end if - - ix_shift = mod((ix-1) + topo_shift, topo_x) + 1 - topo(ix_shift:(ix_shift+tile_x-1),iy:(iy+tile_y-1)) = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) - - end do - end do deallocate(tile) @@ -738,6 +659,7 @@ function read_global_30s_landuse(path) result(iErr) ny = tile_y nz = 1 + allocate(landuse(tile_x*tile_y*cum_local_tiles)) do il=1,cum_local_tiles ix = local_tile_x_land(il) @@ -755,38 +677,14 @@ function read_global_30s_landuse(path) result(iErr) return end if - !ix_shift = mod((ix-1) + ishift, topo_x) + 1 ix_shift = (il-1)*tile_x*tile_y + 1 - call mpas_log_write(' Reading tile ($i/$i) into landuse_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) - !call mpas_log_write(' topo_x,y ($i,$i) ', intArgs=(/topo_x, topo_y/)) - !write(message,fmt='(A,i18,i18,A,i18,A,i18)') 'Reading tile (',il,cum_local_tiles,') topo_new start=', ix_shift_long,' end=',ix_shift_long+tile_x*tile_y-1 - !call mpas_log_write(message) - landuse_new(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) - !write(message,fmt='(A,i18)') 'At 25,16 for landuse_new for (i,j):',landuse_new(ix_shift+25+16*tile_x) - !call mpas_log_write(message) + !call mpas_log_write(' Reading tile ($i/$i) into landuse_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) + + landuse(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) + end do - - do iy=1,topo_y,tile_y - do ix=1,topo_x,tile_x - write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//'/landuse_30s/', ix, '-', (ix+tile_x-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) - tile(:,:,:) = tile(:,:,:) * scalefactor - if (istatus /= 0) then - call mpas_log_write('Error reading landuse tile '//trim(filename)) - iErr = 1 - return - end if - - landuse(ix:(ix+tile_x-1),iy:(iy+tile_y-1)) = int(tile(1:tile_x,1:tile_y,1), kind=I1KIND) - - end do - end do - deallocate(tile) iErr = 0 @@ -832,33 +730,6 @@ integer (kind=I1KIND) function get_dom_landmask(nx, ny) end function get_dom_landmask - integer (kind=I1KIND) function get_dom_landmask_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - real (kind=RKIND) :: xland - xland = 0.0_RKIND - - ! Get dominant land/water mask in the box - do j=1,ny - do i=1,nx - if (box_landuse_new(i,j) /= WATER) then - xland = xland + 1.0_RKIND - end if - end do - end do - xland = xland / real(nx*ny,kind=RKIND) - - if (xland >= 0.5_RKIND) then - get_dom_landmask_new = 1_I1KIND - else - get_dom_landmask_new = 0_I1KIND - end if - - end function get_dom_landmask_new - subroutine get_box_points(lat, lon, nx, ny) @@ -870,7 +741,6 @@ subroutine get_box_points(lat, lon, nx, ny) integer :: i, j, ii, jj, ic, jc, l real (kind=RKIND) :: sg_lat - !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/lon,lat/), intArgs=(/nx,ny/)) ! ! ! Find coordinates in global topography array of the box center @@ -909,14 +779,9 @@ subroutine get_box_points(lat, lon, nx, ny) l = l + 1 local_box_x(l) = ii local_box_y(l) = jj - - !if (i == 25 .and. j ==16) call mpas_log_write('In get box for 25,16 (i,j):($i,$i)', intArgs=(/ii,jj/)) - - end do end do - num_local_box_points = l end subroutine get_box_points @@ -948,24 +813,15 @@ subroutine get_box(lat, lon, nx, ny) real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny - - integer :: i, j, ii, jj, ic, jc + integer :: itile, ix_shift, ix, ix_land, ix_topo, jx, ibox, tile_index_i, tile_index_j + integer :: tile_base, tile_offset, i, j real (kind=RKIND) :: sg_lat - - !call mpas_log_write('---- get box lon, lat ($10.5f, $10.5f) nx, ny ($i, $i)', realArgs=(/lon,lat/), intArgs=(/nx,ny/)) - ! - ! Find coordinates in global topography array of the box center - ! - ic = nint((lon - start_lon) * pts_per_degree) + 1 - jc = nint((lat - start_lat) * pts_per_degree) + 1 - - if (ic <= 0) ic = ic + topo_x - if (ic > topo_x) ic = ic - topo_x - + character(len=StrKIND):: message + if (associated(box)) deallocate(box) - allocate(box(nx,ny)) - + allocate(box(nx, ny)) + if (associated(box_landuse)) deallocate(box_landuse) allocate(box_landuse(nx,ny)) @@ -977,132 +833,47 @@ subroutine get_box(lat, lon, nx, ny) ! the periodicity in the longitude coordinate, as well as the poles ! box_mean = 0.0 - do j=1,ny - do i=1,nx - - ii = i - nx/2 + ic - jj = j - ny/2 + jc - - if (jj <= 0) then - jj = -jj + 1 - ii = ii + topo_y - end if - if (jj > topo_y) then - jj = topo_y - (jj - topo_y - 1) - ii = ii + topo_y - end if - do while (ii <= 0) - ii = ii + topo_x - end do - do while (ii > topo_x) - ii = ii - topo_x - end do - - box(i,j) = topo(ii,jj) - box_landuse(i,j) = landuse(ii,jj) - 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 - ! - box_mean = box_mean / real(nx*ny, RKIND) - - end subroutine get_box - - - subroutine get_box_new(lat, lon, nx, ny) - - implicit none - - real (kind=RKIND), intent(in) :: lat, lon - integer, intent(in) :: nx, ny - integer :: itile, ix_shift, ix, ix_land, ix_topo, jx, ibox, tile_index_i, tile_index_j - integer :: tile_base, tile_offset, i, j - real (kind=RKIND) :: sg_lat - character(len=StrKIND):: message - - - if (associated(box_new)) deallocate(box_new) - allocate(box_new(nx, ny)) - - if (associated(box_landuse_new)) deallocate(box_landuse_new) - allocate(box_landuse_new(nx, ny)) - - - ! if (associated(box_landuse)) deallocate(box_landuse) - ! allocate(box_landuse(nx,ny)) - - ! if (associated(dxm)) deallocate(dxm) - ! allocate(dxm(nx,ny)) - - ! - ! Extract sub-array (box) from global array; must properly account for - ! the periodicity in the longitude coordinate, as well as the poles - ! - box_mean_new = 0.0 ibox = 1 do j=1,ny do i=1,nx - !j = 25 - !i = 16 - ibox = (j-1)*nx + i - ! pixels on the global topo grid - ix = local_box_x(ibox) !+ topo_shift - ix_topo = local_box_x(ibox) - jx = local_box_y(ibox) - - if (ix - topo_shift > 0 ) then - ix = ix - topo_shift - else - ix = ix + topo_shift - end if - - tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 - tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 - + ibox = (j-1)*nx + i + ! pixels on the global topo grid + ix = local_box_x(ibox) !+ topo_shift + ix_topo = local_box_x(ibox) + jx = local_box_y(ibox) + if (ix - topo_shift > 0 ) then + ix = ix - topo_shift + else + ix = ix + topo_shift + end if - itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) + tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 + tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 - !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) + itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) - !ix_shift = (itile-1)*tile_x*tile_y + 1 - tile_base = (itile-1)*tile_x*tile_y + 1 - !tile_base = local_tile_x(tile_index_i) + local_tile_y(tile_index_j)* + - tile_offset = (ix - local_tile_x(itile)) + (jx - local_tile_y(itile)) * tile_x - !tile_offset = mod(ix, local_tile_x(tile_index_i)) + mod(jx, local_tile_y(tile_index_j)) * tile_x - - box_new(i,j) = topo_new(tile_base + tile_offset) - box_mean_new = box_mean_new + box_new(i,j) + !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) - box_landuse_new(i,j) = landuse_new(tile_base + tile_offset) - sg_lat = (start_lat + (real(jx-1,RKIND) + 0.5) / pts_per_degree) / rad2deg ! Add 0.5 for cell center - dxm(i,j) = sg_delta * cos(sg_lat) + tile_base = (itile-1)*tile_x*tile_y + 1 + tile_offset = (ix - local_tile_x(itile)) + (jx - local_tile_y(itile)) * tile_x + + box(i,j) = topo(tile_base + tile_offset) + box_landuse(i,j) = landuse(tile_base + tile_offset) + sg_lat = (start_lat + (real(jx-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) - !call mpas_log_write('In box_new for (i,j):($i,$i), (tile_base,tile_offset):($i,$i)', intArgs=(/i,j,tile_base,tile_offset/)) - if (box_new(i,j) .ne. box(i,j)) call mpas_log_write('In box_new for (i,j):($i,$i), box_new(i,j)= $r box(i,j)= $r ', intArgs=(/i,j/), realArgs=(/box_new(i,j),box(i,j)/)) - - if (box_landuse_new(i,j) .ne. landuse(ix_topo,jx)) then - write(message,fmt='(A,i10,A,i10,A,i18,A,i18)') 'In box_landuse_new for (i,j):(',i,',',j,'), box_landuse_new(i,j)=', box_landuse_new(i,j),' box_landuse(i,j)=',box_landuse(i,j) - call mpas_log_write(message) - !call mpas_log_write('In box_landuse_new for (i,j):($i,$i), box_landuse_new(i,j)= $i box_landuse(i,j)= $i ', intArgs=(/i,j,box_landuse_new(i,j),box_landuse(i,j)/)) - end if - ibox = ibox + 1 + ibox = ibox + 1 end do end do ! ! Compute mean topography in the extracted box ! - box_mean_new = box_mean_new / real(nx*ny, RKIND) + box_mean = box_mean / real(nx*ny, RKIND) - end subroutine get_box_new + end subroutine get_box !*********************************************************************** @@ -1136,26 +907,6 @@ real (kind=RKIND) function get_var(nx, ny) end function get_var - real (kind=RKIND) function get_var_new(nx, ny) - - implicit none - - integer :: i, j - integer, intent(in) :: nx, ny - real (kind=RKIND) :: s2 - - s2 = 0.0 - - do j=1,ny - do i=1,nx - s2 = s2 + (box_new(i,j) - box_mean_new)**2 - end do - end do - - get_var_new = sqrt(s2 / real(nx*ny,RKIND)) - - end function get_var_new - !*********************************************************************** ! @@ -1233,71 +984,6 @@ real (kind=RKIND) function get_con(nx, ny) end function get_con - real (kind=RKIND) function get_con_new(nx, ny) - - implicit none - - integer :: i, j - integer, intent(in) :: nx, ny - real (kind=RKIND) :: s2, s4, var, xland, mean_land, mean_water, oro - - s2 = 0.0 - s4 = 0.0 - mean_land = 0.0 - mean_water = 0.0 - xland = 0.0 - - ! - ! Compute grid-box mean - ! - do j=1,ny - do i=1,nx - if (box_landuse_new(i,j) /= WATER) then - xland = xland + 1.0 - mean_land = mean_land + box_new(i,j) - else - mean_water = mean_water + box_new(i,j) - end if - end do - end do - if (xland > 0.0) then - mean_land = mean_land / xland - end if - if (xland < real(nx*ny,kind=RKIND)) then - mean_water = mean_water / (real(nx*ny,kind=RKIND) - xland) - end if - xland = xland / real(nx*ny,kind=RKIND) - - if (xland >= 0.5_RKIND) then - oro = mean_land - else - oro = mean_water - end if - - do j=1,ny - do i=1,nx - s2 = s2 + (box_new(i,j) - box_mean_new)**2 - s4 = s4 + (box_new(i,j) - oro)**4 - end do - end do - - var = s2 / real(nx*ny,RKIND) - - if (sqrt(var) < 1.0) then - get_con_new = 0.0 - else - get_con_new = s4 / (var**2 * real(nx*ny,RKIND)) - end if - - ! - ! Zero-ing all convexity statistics over dominantly water points. - ! - if (xland < 0.5_RKIND) then - get_con_new = 0.0 - end if - - end function get_con_new - !*********************************************************************** ! @@ -1338,32 +1024,6 @@ real (kind=RKIND) function get_oa1(nx, ny) end function get_oa1 - real (kind=RKIND) function get_oa1_new(nx, ny) - - implicit none - - integer :: i, j - integer, intent(in) :: nx, ny - integer :: nu, nd - - nu = 0 - nd = 0 - do j=1,ny - do i=1,nx/2 - if (box_new(i,j) > box_mean_new) nu = nu + 1 - end do - do i=nx/2+1,nx - if (box_new(i,j) > box_mean_new) nd = nd + 1 - end do - end do - - if (nu + nd > 0) then - get_oa1_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) - else - get_oa1_new = 0.0 - end if - - end function get_oa1_new !*********************************************************************** ! @@ -1406,34 +1066,6 @@ real (kind=RKIND) function get_oa2(nx, ny) end function get_oa2 - real (kind=RKIND) function get_oa2_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nu, nd - - nu = 0 - nd = 0 - do j=1,ny/2 - do i=1,nx - if (box_new(i,j) > box_mean_new) nu = nu + 1 - end do - end do - do j=ny/2+1,ny - do i=1,nx - if (box_new(i,j) > box_mean_new) nd = nd + 1 - end do - end do - - if (nu + nd > 0) then - get_oa2_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) - else - get_oa2_new = 0.0 - end if - - end function get_oa2_new !*********************************************************************** ! @@ -1478,36 +1110,6 @@ real (kind=RKIND) function get_oa3(nx, ny) end function get_oa3 - real (kind=RKIND) function get_oa3_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nu, nd - real (kind=RKIND) :: ratio - - nu = 0 - nd = 0 - ratio = real(ny,RKIND)/real(nx,RKIND) - do j=1,ny - do i=1,nx - if (nint(real(i,RKIND) * ratio) < (ny - j)) then - if (box_new(i,j) > box_mean_new) nu = nu + 1 - else - if (box_new(i,j) > box_mean_new) nd = nd + 1 - end if - end do - end do - - if (nu + nd > 0) then - get_oa3_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) - else - get_oa3_new = 0.0 - end if - - end function get_oa3_new - !*********************************************************************** ! ! function get_oa4 @@ -1551,35 +1153,6 @@ real (kind=RKIND) function get_oa4(nx, ny) end function get_oa4 - real (kind=RKIND) function get_oa4_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nu, nd - real (kind=RKIND) :: ratio - - nu = 0 - nd = 0 - ratio = real(ny,RKIND)/real(nx,RKIND) - do j=1,ny - do i=1,nx - if (nint(real(i,RKIND) * ratio) < j) then - if (box_new(i,j) > box_mean_new) nu = nu + 1 - else - if (box_new(i,j) > box_mean_new) nd = nd + 1 - end if - end do - end do - - if (nu + nd > 0) then - get_oa4_new = real((nu - nd),RKIND) / real((nu + nd),RKIND) - else - get_oa4_new = 0.0 - end if - -end function get_oa4_new !*********************************************************************** ! ! function get_ol1 @@ -1613,31 +1186,6 @@ real (kind=RKIND) function get_ol1(nx, ny) end function get_ol1 - - real (kind=RKIND) function get_ol1_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nw - integer :: nt - - nw = 0 - nt = 0 - - do j=ny/4,3*ny/4 - do i=1,nx - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - - get_ol1_new = real(nw,RKIND) / real(nt,RKIND) - - end function get_ol1_new - - !*********************************************************************** ! ! function get_ol2 @@ -1672,31 +1220,6 @@ real (kind=RKIND) function get_ol2(nx, ny) end function get_ol2 - - real (kind=RKIND) function get_ol2_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nw - integer :: nt - - nw = 0 - nt = 0 - - do j=1,ny - do i=nx/4,3*nx/4 - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - - get_ol2_new = real(nw,RKIND) / real(nt,RKIND) - - end function get_ol2_new - - !*********************************************************************** ! ! function get_ol3 @@ -1736,35 +1259,6 @@ real (kind=RKIND) function get_ol3(nx, ny) end function get_ol3 - real (kind=RKIND) function get_ol3_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nw - integer :: nt - - nw = 0 - nt = 0 - - do j=1,ny/2 - do i=1,nx/2 - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - do j=ny/2+1,ny - do i=nx/2+1,nx - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - - get_ol3_new = real(nw,RKIND) / real(nt,RKIND) - - end function get_ol3_new - !*********************************************************************** ! @@ -1806,35 +1300,6 @@ real (kind=RKIND) function get_ol4(nx, ny) end function get_ol4 - real (kind=RKIND) function get_ol4_new(nx, ny) - - implicit none - integer, intent(in) :: nx, ny - - integer :: i, j - integer :: nw - integer :: nt - - nw = 0 - nt = 0 - - do j=ny/2+1,ny - do i=1,nx/2 - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - do j=1,ny/2 - do i=nx/2+1,nx - if (box_new(i,j) > hc) nw = nw + 1 - nt = nt + 1 - end do - end do - - get_ol4_new = real(nw,RKIND) / real(nt,RKIND) - - end function get_ol4_new - !*********************************************************************** ! ! function get_elvmax @@ -1967,7 +1432,6 @@ real (kind=RKIND) function get_hgamma(nx, ny) end function get_hgamma - !*********************************************************************** ! ! function get_hsigma From cb39d6858f3b2c8aa717c39cbca10539e5a41179 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 27 Aug 2024 23:42:44 +0000 Subject: [PATCH 04/27] more cleanups --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 101 +++++-------------- 1 file changed, 23 insertions(+), 78 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0244b7a292..af1cc50420 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -127,8 +127,6 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND), pointer :: config_topo_data character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - real (kind=RKIND) :: var2d_tmp, con_tmp, oa1_tmp, oa2_tmp, oa3_tmp, oa4_tmp, ol1_tmp, ol2_tmp, ol3_tmp, ol4_tmp, hlanduse_tmp - real (kind=RKIND) :: elvmax_1, htheta_1, hgamma_1, hsigma_1, elvmax_tmp, htheta_tmp, hgamma_tmp, hsigma_tmp ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell @@ -204,20 +202,6 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'gamma', hgamma) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) - ! - ! Get number of points to extract in the zonal direction - ! - ! if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - ! nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) - ! else - ! nx = topo_x / 2 - ! end if - - ! ! - ! ! Get number of points to extract in the meridional direction - ! ! - ! ny = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re)) - nx = 22000 ny = 300 @@ -225,10 +209,7 @@ function compute_gwd_fields(domain) result(iErr) allocate(local_box_y(nx*ny)) allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells - - - - + ! ! It is possible that this code is called before the mesh fields have been scaled @@ -249,10 +230,8 @@ function compute_gwd_fields(domain) result(iErr) call mpas_log_write('in the computation of GWD static fields.') end if - !iErr = determine_complete_tile_set(nCells, latCell, lonCell, config_gwd_cell_scaling, onUnitSphere, sphere_radius) max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) - call mpas_log_write(' hi from compute_gwd_fields -- before allocate max_local_tiles:$i', intArgs=(/max_local_tiles/)) allocate(local_tile_x(max_local_tiles)) allocate(local_tile_x_land(max_local_tiles)) allocate(local_tile_y(max_local_tiles)) @@ -264,30 +243,26 @@ function compute_gwd_fields(domain) result(iErr) ! cum_local_tiles=0 do iCell=1,nCells - !do iCell=5,25 - - ! - ! First, get an estimate of the mean diameter (in meters) of the grid - ! cell by averaging the distances to each of the neighboring cells - ! - dc = 0.0 - do i=1,nEdgesOnCell(iCell) - dc = dc + dcEdge(edgesOnCell(i,iCell)) - end do - dc = dc / real(nEdgesOnCell(iCell),RKIND) - if (onUnitSphere) then - dc = dc * sphere_radius - end if - dc = dc * config_gwd_cell_scaling - ! - call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) - - call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) + ! + ! First, get an estimate of the mean diameter (in meters) of the grid + ! cell by averaging the distances to each of the neighboring cells + ! + dc = 0.0 + do i=1,nEdgesOnCell(iCell) + dc = dc + dcEdge(edgesOnCell(i,iCell)) + end do + dc = dc / real(nEdgesOnCell(iCell),RKIND) + if (onUnitSphere) then + dc = dc * sphere_radius + end if + dc = dc * config_gwd_cell_scaling - iErr = determine_complete_tile_set(nx, ny) + call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) + call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) - if (cum_local_tiles >= max_local_tiles) exit + iErr = determine_complete_tile_set(nx, ny) + if (cum_local_tiles >= max_local_tiles) exit end do call mpas_log_write(' End of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/cum_local_tiles, max_local_tiles/)) @@ -305,7 +280,6 @@ function compute_gwd_fields(domain) result(iErr) end if do iCell=1,nCells - !do iCell=5,25 ! ! First, get an estimate of the mean diameter (in meters) of the grid @@ -366,7 +340,6 @@ function compute_gwd_fields(domain) result(iErr) end do - ! Smooth variance at isolated points do iCell = 1,nCells sum_landuse = 0_I1KIND @@ -474,29 +447,8 @@ function determine_complete_tile_set(nx, ny) result(iErr) integer :: iErr integer, intent(in) :: nx, ny - !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 :: ix, jx, ii, l,i, iCell, itile, jtile - !max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) - !call mpas_log_write(' hi from determine_complete_tile_set -- before allocate max_local_tiles:$i', intArgs=(/max_local_tiles/)) - !allocate(local_tile_x(max_local_tiles)) - !allocate(local_tile_y(max_local_tiles)) - - - - - ! - ! Cut out a rectangular piece of the global 30-arc-second topography - ! data that is centered at the lat/lon of the current cell being - ! processed and that is just large enough to cover the cell. The - ! rectangular array of topography data is stored in the module - ! variable 'box', and the dimensions of this array are given by the - ! module variables 'nx' and 'ny'. The get_box() routine also - ! computes the mean elevation in the array and stores that value in - ! the module variable 'box_mean'. - ! ii = cum_local_tiles do l=1,nx*ny @@ -550,7 +502,6 @@ function read_global_30s_topo(path, sub_path) result(iErr) integer (c_int) :: istatus integer :: ix, iy, ix_shift, il - integer(kind=I8KIND):: ix_shift_long integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile @@ -574,8 +525,6 @@ function read_global_30s_topo(path, sub_path) result(iErr) allocate(topo(tile_x*tile_y*cum_local_tiles)) - call mpas_log_write('---- topo_x:$i topo_y: $i', intArgs=(/topo_x,topo_y/)) - call mpas_log_write('---- tile_x:$i tile_y: $i', intArgs=(/tile_x,tile_y/)) ! ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 @@ -601,13 +550,10 @@ function read_global_30s_topo(path, sub_path) result(iErr) end if ix_shift = (il-1)*tile_x*tile_y + 1 - !call mpas_log_write(' Reading tile ($i/$i) into topo_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) - topo(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) end do - deallocate(tile) iErr = 0 @@ -678,13 +624,10 @@ function read_global_30s_landuse(path) result(iErr) end if ix_shift = (il-1)*tile_x*tile_y + 1 - !call mpas_log_write(' Reading tile ($i/$i) into landuse_new ($i,$i) ishift: $i)', intArgs=(/il,cum_local_tiles,ix_shift,ix_shift+tile_x*tile_y-1,topo_shift/)) - landuse(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) end do - deallocate(tile) iErr = 0 @@ -820,7 +763,7 @@ subroutine get_box(lat, lon, nx, ny) if (associated(box)) deallocate(box) - allocate(box(nx, ny)) + allocate(box(nx,ny)) if (associated(box_landuse)) deallocate(box_landuse) allocate(box_landuse(nx,ny)) @@ -868,6 +811,8 @@ subroutine get_box(lat, lon, nx, ny) ibox = ibox + 1 end do end do + + ! ! Compute mean topography in the extracted box ! @@ -907,7 +852,6 @@ real (kind=RKIND) function get_var(nx, ny) end function get_var - !*********************************************************************** ! ! function get_con @@ -984,7 +928,6 @@ real (kind=RKIND) function get_con(nx, ny) end function get_con - !*********************************************************************** ! ! function get_oa1 @@ -1186,6 +1129,7 @@ real (kind=RKIND) function get_ol1(nx, ny) end function get_ol1 + !*********************************************************************** ! ! function get_ol2 @@ -1432,6 +1376,7 @@ real (kind=RKIND) function get_hgamma(nx, ny) end function get_hgamma + !*********************************************************************** ! ! function get_hsigma From 4938b7fb096c8141e1a8b97ad95c3f83f2804d2c Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 29 Aug 2024 04:07:00 +0000 Subject: [PATCH 05/27] some debug --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index af1cc50420..0bbe2ad250 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -202,8 +202,8 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'gamma', hgamma) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) - nx = 22000 - ny = 300 + nx = topo_x + ny = topo_y allocate(local_box_x(nx*ny)) allocate(local_box_y(nx*ny)) @@ -296,6 +296,9 @@ function compute_gwd_fields(domain) result(iErr) dc = dc * config_gwd_cell_scaling call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) + + if (iCell==2621501) call mpas_log_write(' iCell is selected nx, ny = $i $i, lat lon $r $r', intArgs=(/nx, ny/), realArgs=(/latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg/)) + ! ! Cut out a rectangular piece of the global 30-arc-second topography ! data that is centered at the lat/lon of the current cell being From e939c677594d1b42e24899ebaf92744c7ee44e68 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 29 Aug 2024 22:49:10 +0000 Subject: [PATCH 06/27] some test --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0bbe2ad250..0cc2b0f883 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -258,7 +258,7 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) + call get_box_size_from_lat_lon(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg,dc,nx,ny) call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) iErr = determine_complete_tile_set(nx, ny) @@ -384,8 +384,8 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) ! ! Get number of points to extract in the zonal direction ! - if (cos(lat) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat))) + if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) else nx = topo_x / 2 end if From 30a433eec39bfb3bb906219b459967de9973fce9 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 29 Aug 2024 22:56:57 +0000 Subject: [PATCH 07/27] some test --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0cc2b0f883..2779b1496d 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -295,7 +295,7 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) + call get_box_size_from_lat_lon(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc, nx, ny) if (iCell==2621501) call mpas_log_write(' iCell is selected nx, ny = $i $i, lat lon $r $r', intArgs=(/nx, ny/), realArgs=(/latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg/)) From d01117b135cd578cccf719eed929e429fc1a3f7c Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 30 Aug 2024 20:32:04 +0000 Subject: [PATCH 08/27] reverting to new scheme rad2deg --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 23 +++++++++----------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 2779b1496d..49f74fda56 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -258,8 +258,8 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - call get_box_size_from_lat_lon(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg,dc,nx,ny) - call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) + call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) + call get_box_points(latCell(iCell),lonCell(iCell), nx, ny) iErr = determine_complete_tile_set(nx, ny) if (cum_local_tiles >= max_local_tiles) exit @@ -295,9 +295,7 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - call get_box_size_from_lat_lon(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc, nx, ny) - - if (iCell==2621501) call mpas_log_write(' iCell is selected nx, ny = $i $i, lat lon $r $r', intArgs=(/nx, ny/), realArgs=(/latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg/)) + call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) ! ! Cut out a rectangular piece of the global 30-arc-second topography @@ -309,8 +307,8 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box_points(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny) - call get_box(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, nx, ny) + call get_box_points(latCell(iCell),lonCell(iCell), nx, ny) + call get_box(nx, ny) ! ! With a box of 30-arc-second data for the current grid cell, call @@ -384,8 +382,8 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) ! ! Get number of points to extract in the zonal direction ! - if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) + if (cos(lat) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat))) else nx = topo_x / 2 end if @@ -691,8 +689,8 @@ subroutine get_box_points(lat, lon, nx, ny) ! ! Find coordinates in global topography array of the box center ! - ic = nint((lon - start_lon) * pts_per_degree) + 1 - jc = nint((lat - start_lat) * pts_per_degree) + 1 + ic = nint((lon*rad2deg - start_lon) * pts_per_degree) + 1 + jc = nint((lat*rad2deg - start_lat) * pts_per_degree) + 1 if (ic <= 0) ic = ic + topo_x if (ic > topo_x) ic = ic - topo_x @@ -753,11 +751,10 @@ end subroutine get_box_points !> this subroutine and stored in the module variable 'box_mean'. ! !----------------------------------------------------------------------- - subroutine get_box(lat, lon, nx, ny) + subroutine get_box(nx, ny) implicit none - real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny integer :: itile, ix_shift, ix, ix_land, ix_topo, jx, ibox, tile_index_i, tile_index_j integer :: tile_base, tile_offset, i, j From da1741154e8dd44e193d8c6f1826ac25399e6ecd Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 30 Aug 2024 23:12:26 +0000 Subject: [PATCH 09/27] wip --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 264 ++++++++++++++++--- 1 file changed, 225 insertions(+), 39 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 49f74fda56..ff77cce0af 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -19,6 +19,15 @@ module mpas_init_atm_gwd private + type :: mpas_gwd_tile_type + + real (kind=R4KIND), dimension(:), pointer :: array + integer :: tile_start_x, tile_start_y + ! linked list next pointer + type (mpas_gwd_tile_type), pointer :: next => null() + + end type mpas_gwd_tile_type + interface subroutine read_geogrid(fname, rarray, nx, ny, nz, isigned, endian, & wordsize, status) bind(C) @@ -242,42 +251,8 @@ function compute_gwd_fields(domain) result(iErr) ! grid cell in the mesh. ! cum_local_tiles=0 - do iCell=1,nCells - - ! - ! First, get an estimate of the mean diameter (in meters) of the grid - ! cell by averaging the distances to each of the neighboring cells - ! - dc = 0.0 - do i=1,nEdgesOnCell(iCell) - dc = dc + dcEdge(edgesOnCell(i,iCell)) - end do - dc = dc / real(nEdgesOnCell(iCell),RKIND) - if (onUnitSphere) then - dc = dc * sphere_radius - end if - dc = dc * config_gwd_cell_scaling - - call get_box_size_from_lat_lon(latCell(iCell),lonCell(iCell),dc,nx,ny) - call get_box_points(latCell(iCell),lonCell(iCell), nx, ny) - - iErr = determine_complete_tile_set(nx, ny) - if (cum_local_tiles >= max_local_tiles) exit - end do call mpas_log_write(' End of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/cum_local_tiles, max_local_tiles/)) - - iErr = read_global_30s_topo(geog_data_path, geog_sub_path) - if (iErr /= 0) then - call mpas_log_write('Error reading global 30-arc-sec topography for GWD statistics', messageType=MPAS_LOG_ERR) - return - end if - - iErr = read_global_30s_landuse(geog_data_path) - if (iErr /= 0) then - call mpas_log_write('Error reading global 30-arc-sec landuse for GWD statistics', messageType=MPAS_LOG_ERR) - return - end if do iCell=1,nCells @@ -307,8 +282,8 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box_points(latCell(iCell),lonCell(iCell), nx, ny) - call get_box(nx, ny) + call get_box_points(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path) + call get_box(nx, ny, geog_data_path, geog_sub_path) ! ! With a box of 30-arc-second data for the current grid cell, call @@ -431,6 +406,46 @@ function locate_tile_in_set(tile_index_i, tile_index_j, is_topo) result (tile_in integer :: itile, tile_index + do itile = 1, cum_local_tiles + if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then + tile_index = itile + return + end if + end do + + tile_index = -1 + + end function locate_tile_in_set + + + function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_index) + + implicit none + type(mpas_gwd_tile_type), pointer :: tilesHead + type(mpas_gwd_tile_type), pointer :: thisTile + integer, intent(in) :: tile_index_i + integer, intent(in) :: tile_index_j + logical, intent(in) :: is_topo + + integer :: itile, tile_index + + ! loop over tiles + thisTile => tilesHead + + do while (associated(thisTile)) + + if(thisTile=>tile_start_x==tile_index_i .and. thisTile=>tile_start_y==tile_index_j) then + tile_index = itile + return + end if + + thisTile => thisTile % next + + enddo ! associated(thisTile) + + + + do itile = 1, cum_local_tiles if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then tile_index = itile @@ -476,6 +491,151 @@ function determine_complete_tile_set(nx, ny) result(iErr) end function determine_complete_tile_set + function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny) result(iErr) + + implicit none + integer :: iErr + integer, intent(in) :: box_x, box_y, nx, ny + type(mpas_gwd_tile_type), pointer, intent(in) :: tilesHead + type(mpas_gwd_tile_type), pointer :: thisTile + + integer :: ix, jx, ii, l,i, iCell, tile_start_x, tile_start_y + + ii = cum_local_tiles + + + if (box_x > topo_shift) then + tile_start_x = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 + else + tile_start_x = floor( real(box_x + topo_shift - 1) / real(tile_x)) * tile_x + 1 + end if + tile_start_y = floor( real(box_y - 1) / real(tile_y)) * tile_y + 1 + + + ! loop over tiles + thisTile => tilesHead + + do while (associated(thisTile)) + + if(thisTile%tile_start_x==tile_start_x .and. thisTile%tile_start_y==tile_start_y) then + exit + end if + + thisTile => thisTile % next + + enddo ! associated(thisTile) + + if (.not. associated(thisTile)) then + thisTile = add_tile(tilesHead, tile_start_x, tile_start_y) + end if + + + if(unique_tile(itile, jtile, ii)) then + ii = ii + 1 + local_tile_x(ii) = itile + local_tile_x_land(ii) = floor( real(box_y - 1) / real(tile_x)) * tile_x + 1 + local_tile_y(ii) = jtile + !call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) + !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) + end if + + + cum_local_tiles = ii + + ! + end function determine_complete_tile_set + + + function add_tile(tilesHead, tile_start_x, tile_start_y) result(tilesHead) + + implicit none + integer :: iErr + integer, intent(in) :: tile_start_x, tile_start_y + type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead + type(mpas_gwd_tile_type), pointer :: newTile + + allocate(newTile) + allocate(newTile%array(topo_x*topo_y)) + newTile%tile_start_x = tile_start_x + newTile%tile_start_y = tile_start_y + newTile%next => tilesHead + + read_30s_topo_tile(path, sub_path, newTile%array, newTile%tile_start_x, newTile%tile_start_y) + tilesHead => newTile + + cum_local_tiles = cum_local_tiles + 1 + + end function add_tile + + function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) result(iErr) + + implicit none + + character(len=*), intent(in) :: path + character(len=*), intent(in) :: sub_path + + integer :: iErr + + 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 :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile + + integer (c_int) :: istatus + integer :: ix, iy, ix_shift, il + integer (c_int) :: isigned, endian, wordsize, nx, ny, nz + real (c_float) :: scalefactor + real (c_float), dimension(:,:,:), pointer, contiguous :: tile + type (c_ptr) :: tile_ptr + character(len=StrKIND) :: filename + character(len=StrKIND):: message + character(kind=c_char), dimension(StrKIND+1) :: c_filename + + allocate(tile(tile_x+2*tile_bdr,tile_y+2*tile_bdr,1)) + tile_ptr = c_loc(tile) + + isigned = 1 + endian = 0 + wordsize = 2 + scalefactor = 1.0 + nx = tile_x + 2*tile_bdr + ny = tile_y + 2*tile_bdr + nz = 1 + + topo_shift = 0 + + allocate(topo(tile_x*tile_y)) + + ! + ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location + ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 + ! + if (trim(sub_path) == 'topo_gmted2010_30s/') then + topo_shift = topo_x / 2 + end if + + 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) + call mpas_f_to_c_string(filename, c_filename) + call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & + wordsize, istatus) + tile(:,:,:) = tile(:,:,:) * scalefactor + if (istatus /= 0) then + call mpas_log_write('Error reading topography tile '//trim(filename), messageType=MPAS_LOG_ERR) + iErr = 1 + return + end if + + ix_shift = (il-1)*tile_x*tile_y + 1 + topo(1:(tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + + deallocate(tile) + + iErr = 0 + + end function read_30s_topo_tile + + + !*********************************************************************** ! ! function read_global_30s_topo @@ -675,9 +835,11 @@ end function get_dom_landmask - subroutine get_box_points(lat, lon, nx, ny) + subroutine get_box_points(lat, lon, nx, ny, path, sub_path) implicit none + character(len=*), intent(in) :: path + character(len=*), intent(in) :: sub_path real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny @@ -721,8 +883,32 @@ subroutine get_box_points(lat, lon, nx, ny) ii = ii - topo_x end do l = l + 1 - local_box_x(l) = ii - local_box_y(l) = jj + + thisTile = get_tile_from_box_point(tilesHead, ii, jj, nx, ny) + + ! ibox = (j-1)*nx + i + + ! if (ix - topo_shift > 0 ) then + ! ix = ix - topo_shift + ! else + ! ix = ix + topo_shift + ! end if + + ! tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 + ! tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 + + ! itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) + + !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) + + ! tile_base = (itile-1)*tile_x*tile_y + 1 + tile_offset = (ii - thisTile%tile_start_x) + (jx - thisTile%tile_start_y) * tile_x + + box(i,j) = thisTile%topo(tile_offset) + box_landuse(i,j) = landuse(tile_offset) + sg_lat = (start_lat + (real(jx-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 From bd42f771130e4e22bccdf49ae1c8b60c5c9420d9 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Wed, 4 Sep 2024 00:28:34 +0000 Subject: [PATCH 10/27] topo is correct --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 136 ++++++++++--------- 1 file changed, 71 insertions(+), 65 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index ff77cce0af..579fd3290c 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -137,6 +137,8 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash + TYPE(mpas_gwd_tile_type), POINTER :: tilesHead + ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse @@ -239,6 +241,7 @@ function compute_gwd_fields(domain) result(iErr) call mpas_log_write('in the computation of GWD static fields.') end if + tilesHead => null() max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) allocate(local_tile_x(max_local_tiles)) @@ -282,15 +285,15 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box_points(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path) - call get_box(nx, ny, geog_data_path, geog_sub_path) + call get_box_points(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path, tilesHead) + !call get_box(nx, ny, geog_data_path, geog_sub_path) ! ! With a box of 30-arc-second data for the current grid cell, call ! subroutines to compute each sub-grid orography statistic ! var2d(iCell) = get_var(nx, ny) - con(iCell) = get_con(nx, ny) + !con(iCell) = get_con(nx, ny) oa1(iCell) = get_oa1(nx, ny) oa2(iCell) = get_oa2(nx, ny) oa3(iCell) = get_oa3(nx, ny) @@ -306,7 +309,7 @@ function compute_gwd_fields(domain) result(iErr) ol3(iCell) = get_ol3(nx, ny) ol4(iCell) = get_ol4(nx, ny) - hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell + !hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell ! elvmax_1 = get_elvmax(nx,ny) ! htheta_1 = get_htheta(nx,ny) @@ -335,8 +338,8 @@ function compute_gwd_fields(domain) result(iErr) end do - deallocate(topo) - deallocate(landuse) + !deallocate(topo) + !deallocate(landuse) deallocate(hlanduse) iErr = 0 @@ -434,7 +437,7 @@ function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_i do while (associated(thisTile)) - if(thisTile=>tile_start_x==tile_index_i .and. thisTile=>tile_start_y==tile_index_j) then + if(thisTile%tile_start_x==tile_index_i .and. thisTile%tile_start_y==tile_index_j) then tile_index = itile return end if @@ -443,19 +446,7 @@ function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_i enddo ! associated(thisTile) - - - - do itile = 1, cum_local_tiles - if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then - tile_index = itile - return - end if - end do - - tile_index = -1 - - end function locate_tile_in_set + end function locate_tile_in_list function determine_complete_tile_set(nx, ny) result(iErr) @@ -491,26 +482,27 @@ function determine_complete_tile_set(nx, ny) result(iErr) end function determine_complete_tile_set - function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny) result(iErr) + function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path) result(thisTile) implicit none integer :: iErr integer, intent(in) :: box_x, box_y, nx, ny 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 integer :: ix, jx, ii, l,i, iCell, tile_start_x, tile_start_y - - ii = cum_local_tiles - if (box_x > topo_shift) then - tile_start_x = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 - else - tile_start_x = floor( real(box_x + topo_shift - 1) / real(tile_x)) * tile_x + 1 - end if + ! if (box_x > topo_shift) then + ! tile_start_x = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 + ! else + ! tile_start_x = floor( real(box_x + topo_shift - 1) / real(tile_x)) * tile_x + 1 + ! end if + tile_start_x = floor( real(box_x - 1) / real(tile_x)) * tile_x + 1 tile_start_y = floor( real(box_y - 1) / real(tile_y)) * tile_y + 1 - + ! loop over tiles thisTile => tilesHead @@ -526,33 +518,22 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny) result(iErr) enddo ! associated(thisTile) if (.not. associated(thisTile)) then - thisTile = add_tile(tilesHead, tile_start_x, tile_start_y) - end if - - - if(unique_tile(itile, jtile, ii)) then - ii = ii + 1 - local_tile_x(ii) = itile - local_tile_x_land(ii) = floor( real(box_y - 1) / real(tile_x)) * tile_x + 1 - local_tile_y(ii) = jtile - !call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) - !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) + thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, path, sub_path) end if - - - cum_local_tiles = ii ! - end function determine_complete_tile_set + end function get_tile_from_box_point - function add_tile(tilesHead, tile_start_x, tile_start_y) result(tilesHead) + function add_tile(tilesHead, tile_start_x, tile_start_y, path, sub_path) result(newTile) implicit none integer :: iErr integer, intent(in) :: tile_start_x, tile_start_y - type(mpas_gwd_tile_type), pointer, intent(inout) :: tilesHead + type(mpas_gwd_tile_type), pointer :: tilesHead type(mpas_gwd_tile_type), pointer :: newTile + character(len=*), intent(in) :: path + character(len=*), intent(in) :: sub_path allocate(newTile) allocate(newTile%array(topo_x*topo_y)) @@ -560,11 +541,13 @@ function add_tile(tilesHead, tile_start_x, tile_start_y) result(tilesHead) newTile%tile_start_y = tile_start_y newTile%next => tilesHead - read_30s_topo_tile(path, sub_path, newTile%array, newTile%tile_start_x, newTile%tile_start_y) + iErr = read_30s_topo_tile(path, sub_path, newTile%array, newTile%tile_start_x, newTile%tile_start_y) tilesHead => newTile cum_local_tiles = cum_local_tiles + 1 + call mpas_log_write('In add_tile for tile_start:($i,$i), cum_local_tiles:$i ', intArgs=(/tile_start_x, tile_start_y,cum_local_tiles /)) + end function add_tile function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) result(iErr) @@ -573,7 +556,9 @@ 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, intent(in) :: topo + integer, intent(in) :: tile_start_x + integer, intent(in) :: tile_start_y integer :: iErr integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second topography @@ -603,7 +588,7 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re topo_shift = 0 - allocate(topo(tile_x*tile_y)) + !allocate(topo(tile_x*tile_y)) ! ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location @@ -626,7 +611,8 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re end if ix_shift = (il-1)*tile_x*tile_y + 1 - topo(1:(tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + !topo(1:(tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) deallocate(tile) @@ -835,7 +821,7 @@ end function get_dom_landmask - subroutine get_box_points(lat, lon, nx, ny, path, sub_path) + subroutine get_box_points(lat, lon, nx, ny, path, sub_path, tilesHead) implicit none character(len=*), intent(in) :: path @@ -843,10 +829,21 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path) real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny + TYPE(mpas_gwd_tile_type), POINTER, intent(in) :: tilesHead + TYPE(mpas_gwd_tile_type), POINTER :: thisTile - integer :: i, j, ii, jj, ic, jc, l + integer :: i, j, ii, jj, ic, jc, l, tile_offset real (kind=RKIND) :: sg_lat + + if (associated(box)) deallocate(box) + allocate(box(nx,ny)) + + if (associated(box_landuse)) deallocate(box_landuse) + allocate(box_landuse(nx,ny)) + + if (associated(dxm)) deallocate(dxm) + allocate(dxm(nx,ny)) ! ! ! Find coordinates in global topography array of the box center @@ -861,6 +858,7 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path) ! Extract sub-array (box) from global array; must properly account for ! the periodicity in the longitude coordinate, as well as the poles ! + box_mean = 0.0 l = 0 do j=1,ny do i=1,nx @@ -884,15 +882,18 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path) end do l = l + 1 - thisTile = get_tile_from_box_point(tilesHead, ii, jj, nx, ny) + + if (ii - topo_shift > 0 ) then + ii = ii - topo_shift + else + ii = ii + topo_shift + end if + + thisTile => get_tile_from_box_point(tilesHead, ii, jj, nx, ny, path, sub_path) ! ibox = (j-1)*nx + i - ! if (ix - topo_shift > 0 ) then - ! ix = ix - topo_shift - ! else - ! ix = ix + topo_shift - ! end if + ! tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 ! tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 @@ -900,19 +901,24 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path) ! itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) - - ! tile_base = (itile-1)*tile_x*tile_y + 1 - tile_offset = (ii - thisTile%tile_start_x) + (jx - thisTile%tile_start_y) * tile_x - box(i,j) = thisTile%topo(tile_offset) - box_landuse(i,j) = landuse(tile_offset) - sg_lat = (start_lat + (real(jx-1,RKIND) + 0.5) / pts_per_degree) / rad2deg ! Add 0.5 for cell center + ! tile_base = (itile-1)*tile_x*tile_y + 1 + tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 + if (tile_offset < 0) then + call mpas_log_write('In box_new for (i,j):($i,$i), tile_offset:$i tile_start_end $i, $i', intArgs=(/ii,jj,tile_offset, thisTile%tile_start_x, thisTile%tile_start_y/)) + call mpas_log_write('In box_new for (i,j):($i,$i), lat:$r, lon:$r', intArgs=(/ii,jj/), realArgs=(/lat*rad2deg, lon*rad2deg /)) + end if + box(i,j) = thisTile%array(tile_offset) + !box_landuse(i,j) = landuse(tile_offset) + 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 + box_mean = box_mean / real(nx*ny, RKIND) + end subroutine get_box_points @@ -989,7 +995,7 @@ subroutine get_box(nx, ny) tile_offset = (ix - local_tile_x(itile)) + (jx - local_tile_y(itile)) * tile_x box(i,j) = topo(tile_base + tile_offset) - box_landuse(i,j) = landuse(tile_base + tile_offset) + !box_landuse(i,j) = landuse(tile_base + tile_offset) sg_lat = (start_lat + (real(jx-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) From b2ae7647c597837b6b45c9ff32a2392e68643bd0 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Wed, 4 Sep 2024 21:27:36 +0000 Subject: [PATCH 11/27] topo and land working --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 125 ++++++++++++++----- 1 file changed, 96 insertions(+), 29 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 579fd3290c..0221a42dec 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -21,7 +21,8 @@ module mpas_init_atm_gwd type :: mpas_gwd_tile_type - real (kind=R4KIND), dimension(:), pointer :: array + real (kind=R4KIND), dimension(:), pointer :: topo_array + real (kind=R4KIND), dimension(:), pointer :: land_array integer :: tile_start_x, tile_start_y ! linked list next pointer type (mpas_gwd_tile_type), pointer :: next => null() @@ -293,7 +294,7 @@ function compute_gwd_fields(domain) result(iErr) ! subroutines to compute each sub-grid orography statistic ! var2d(iCell) = get_var(nx, ny) - !con(iCell) = get_con(nx, ny) + con(iCell) = get_con(nx, ny) oa1(iCell) = get_oa1(nx, ny) oa2(iCell) = get_oa2(nx, ny) oa3(iCell) = get_oa3(nx, ny) @@ -309,7 +310,7 @@ function compute_gwd_fields(domain) result(iErr) ol3(iCell) = get_ol3(nx, ny) ol4(iCell) = get_ol4(nx, ny) - !hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell + hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell ! elvmax_1 = get_elvmax(nx,ny) ! htheta_1 = get_htheta(nx,ny) @@ -492,14 +493,14 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path character(len=*), intent(in) :: path character(len=*), intent(in) :: sub_path - integer :: ix, jx, ii, l,i, iCell, tile_start_x, tile_start_y + integer :: ix, jx, ii, l,i, iCell, tile_start_x, tile_start_y, tile_start_x_topo - ! if (box_x > topo_shift) then - ! tile_start_x = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 - ! else - ! tile_start_x = floor( real(box_x + topo_shift - 1) / real(tile_x)) * tile_x + 1 - ! end if + if (box_x > topo_shift) then + tile_start_x_topo = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 + else + tile_start_x_topo = floor( real(box_x + topo_shift - 1) / real(tile_x)) * tile_x + 1 + end if tile_start_x = floor( real(box_x - 1) / real(tile_x)) * tile_x + 1 tile_start_y = floor( real(box_y - 1) / real(tile_y)) * tile_y + 1 @@ -518,35 +519,37 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path enddo ! associated(thisTile) if (.not. associated(thisTile)) then - thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, path, sub_path) + thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) end if ! end function get_tile_from_box_point - function add_tile(tilesHead, tile_start_x, tile_start_y, path, sub_path) result(newTile) + function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) result(newTile) implicit none integer :: iErr - integer, intent(in) :: tile_start_x, tile_start_y + integer, intent(in) :: tile_start_x, tile_start_y, tile_start_x_topo type(mpas_gwd_tile_type), pointer :: tilesHead type(mpas_gwd_tile_type), pointer :: newTile character(len=*), intent(in) :: path character(len=*), intent(in) :: sub_path allocate(newTile) - allocate(newTile%array(topo_x*topo_y)) + allocate(newTile%topo_array(topo_x*topo_y)) + allocate(newTile%land_array(topo_x*topo_y)) newTile%tile_start_x = tile_start_x newTile%tile_start_y = tile_start_y newTile%next => tilesHead - iErr = read_30s_topo_tile(path, sub_path, newTile%array, newTile%tile_start_x, newTile%tile_start_y) + iErr = read_30s_topo_tile(path, sub_path, newTile%topo_array, tile_start_x_topo, newTile%tile_start_y) + iErr = read_30s_landuse_tile(path, sub_path, newTile%land_array, newTile%tile_start_x, newTile%tile_start_y) tilesHead => newTile cum_local_tiles = cum_local_tiles + 1 - call mpas_log_write('In add_tile for tile_start:($i,$i), cum_local_tiles:$i ', intArgs=(/tile_start_x, tile_start_y,cum_local_tiles /)) + !call mpas_log_write('In add_tile for tile_start:($i,$i), cum_local_tiles:$i ', intArgs=(/tile_start_x, tile_start_y,cum_local_tiles /)) end function add_tile @@ -588,6 +591,14 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re topo_shift = 0 + ! if (tile_start_x - topo_shift > 0 ) then + ! ix = tile_start_x - topo_shift + ! else + ! ix = tile_start_x + topo_shift + ! end if + iy = tile_start_y + ix = tile_start_x + !allocate(topo(tile_x*tile_y)) ! @@ -598,8 +609,10 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re topo_shift = topo_x / 2 end if - 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), 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) call mpas_f_to_c_string(filename, c_filename) call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & wordsize, istatus) @@ -610,7 +623,6 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re return end if - ix_shift = (il-1)*tile_x*tile_y + 1 !topo(1:(tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) @@ -621,6 +633,61 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re end function read_30s_topo_tile + function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start_y) result(iErr) + + implicit none + + character(len=*), intent(in) :: path + character(len=*), intent(in) :: sub_path + REAL(kind=R4KIND), DIMENSION(:), POINTER, intent(in) :: landuse + integer, intent(in) :: tile_start_x + integer, intent(in) :: tile_start_y + integer :: iErr + + integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second landuse + integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second landuse + + integer (c_int) :: istatus + integer :: ix, iy, ix_shift, il + integer (c_int) :: isigned, endian, wordsize, nx, ny, nz + real (c_float) :: scalefactor + real (c_float), dimension(:,:,:), pointer, contiguous :: tile + type (c_ptr) :: tile_ptr + character(len=StrKIND) :: filename + character(kind=c_char), dimension(StrKIND+1) :: c_filename + character(len=StrKIND):: message + + allocate(tile(tile_x,tile_y,1)) + tile_ptr = c_loc(tile) + + isigned = 1 + endian = 0 + wordsize = 1 + scalefactor = 1.0 + nx = tile_x + ny = tile_y + 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) + call mpas_f_to_c_string(filename, c_filename) + call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & + wordsize, istatus) + tile(:,:,:) = tile(:,:,:) * scalefactor + if (istatus /= 0) then + call mpas_log_write('Error reading landuse tile '//trim(filename)) + iErr = 1 + return + end if + + landuse = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) + + deallocate(tile) + + iErr = 0 + + end function read_30s_landuse_tile + !*********************************************************************** ! @@ -883,11 +950,11 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path, tilesHead) l = l + 1 - if (ii - topo_shift > 0 ) then - ii = ii - topo_shift - else - ii = ii + topo_shift - end if + ! if (ii - topo_shift > 0 ) then + ! ii = ii - topo_shift + ! else + ! ii = ii + topo_shift + ! end if thisTile => get_tile_from_box_point(tilesHead, ii, jj, nx, ny, path, sub_path) @@ -904,12 +971,12 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path, tilesHead) ! tile_base = (itile-1)*tile_x*tile_y + 1 tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 - if (tile_offset < 0) then - call mpas_log_write('In box_new for (i,j):($i,$i), tile_offset:$i tile_start_end $i, $i', intArgs=(/ii,jj,tile_offset, thisTile%tile_start_x, thisTile%tile_start_y/)) - call mpas_log_write('In box_new for (i,j):($i,$i), lat:$r, lon:$r', intArgs=(/ii,jj/), realArgs=(/lat*rad2deg, lon*rad2deg /)) - end if - box(i,j) = thisTile%array(tile_offset) - !box_landuse(i,j) = landuse(tile_offset) + !if (tile_offset < 0) then + ! call mpas_log_write('In box_new for (i,j):($i,$i), tile_offset:$i tile_start_end $i, $i', intArgs=(/ii,jj,tile_offset, thisTile%tile_start_x, thisTile%tile_start_y/)) + ! call mpas_log_write('In box_new for (i,j):($i,$i), lat:$r, lon:$r', intArgs=(/ii,jj/), realArgs=(/lat*rad2deg, lon*rad2deg /)) + !end if + box(i,j) = thisTile%topo_array(tile_offset) + box_landuse(i,j) = thisTile%land_array(tile_offset) 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) From db6a95000fff5c34307c507266dc43cc188bd3c7 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Wed, 4 Sep 2024 22:46:23 +0000 Subject: [PATCH 12/27] some cleanup --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 394 +++---------------- 1 file changed, 53 insertions(+), 341 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0221a42dec..a2f45f4e6f 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -24,7 +24,7 @@ module mpas_init_atm_gwd real (kind=R4KIND), dimension(:), pointer :: topo_array real (kind=R4KIND), dimension(:), pointer :: land_array integer :: tile_start_x, tile_start_y - ! linked list next pointer + ! linked list next pointer type (mpas_gwd_tile_type), pointer :: next => null() end type mpas_gwd_tile_type @@ -68,12 +68,10 @@ end subroutine read_geogrid ! Nominal delta-x (in meters) for sub-grid topography cells real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) - real (kind=R4KIND), dimension(:), pointer :: topo ! Global 30-arc-second topography real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell real (kind=RKIND) :: box_mean ! Mean value of topography in box !integer :: nx, ny ! Dimensions of box covering grid cell - integer (kind=I1KIND), dimension(:), pointer :: landuse ! Global 30-arc-second landuse integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell ! NB: At present, only the USGS GLCC land cover dataset is supported, so we can assume 16 == water @@ -83,8 +81,6 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height - integer, dimension(:), pointer :: local_tile_x, local_tile_x_land ! - integer, dimension(:), pointer :: local_tile_y ! integer, dimension(:), pointer :: local_box_x ! integer, dimension(:), pointer :: local_box_y ! integer:: max_local_tiles, cum_local_tiles @@ -245,9 +241,6 @@ function compute_gwd_fields(domain) result(iErr) tilesHead => null() max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) - allocate(local_tile_x(max_local_tiles)) - allocate(local_tile_x_land(max_local_tiles)) - allocate(local_tile_y(max_local_tiles)) topo_shift = topo_x / 2 ! @@ -286,7 +279,7 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box_points(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path, tilesHead) + call get_box(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path, tilesHead) !call get_box(nx, ny, geog_data_path, geog_sub_path) ! @@ -338,9 +331,6 @@ function compute_gwd_fields(domain) result(iErr) end if end do - - !deallocate(topo) - !deallocate(landuse) deallocate(hlanduse) iErr = 0 @@ -376,52 +366,6 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) end subroutine get_box_size_from_lat_lon - - function unique_tile(ix, jx, len) result (is_unique) - - implicit none - - integer, intent(in) :: ix - integer, intent(in) :: jx - integer, intent(in) :: len - - integer :: i - logical :: is_unique - - do i = 1, len - if(local_tile_x(i)==ix .and. local_tile_y(i)==jx) then - is_unique = .false. - return - end if - end do - - is_unique = .true. - - end function unique_tile - - function locate_tile_in_set(tile_index_i, tile_index_j, is_topo) result (tile_index) - - implicit none - - integer, intent(in) :: tile_index_i - integer, intent(in) :: tile_index_j - logical, intent(in) :: is_topo - - integer :: itile, tile_index - - - do itile = 1, cum_local_tiles - if(local_tile_x(itile)==tile_index_i .and. local_tile_y(itile)==tile_index_j) then - tile_index = itile - return - end if - end do - - tile_index = -1 - - end function locate_tile_in_set - - function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_index) implicit none @@ -449,39 +393,6 @@ function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_i end function locate_tile_in_list - function determine_complete_tile_set(nx, ny) result(iErr) - - implicit none - integer :: iErr - integer, intent(in) :: nx, ny - - integer :: ix, jx, ii, l,i, iCell, itile, jtile - - ii = cum_local_tiles - do l=1,nx*ny - - if (local_box_x(l) > topo_shift) then - itile = floor( real(local_box_x(l) - topo_shift - 1) / real(tile_x)) * tile_x + 1 - else - itile = floor( real(local_box_x(l) + topo_shift - 1) / real(tile_x)) * tile_x + 1 - end if - jtile = floor( real(local_box_y(l) - 1) / real(tile_y)) * tile_y + 1 - - if(unique_tile(itile, jtile, ii)) then - ii = ii + 1 - local_tile_x(ii) = itile - local_tile_x_land(ii) = floor( real(local_box_x(l) - 1) / real(tile_x)) * tile_x + 1 - local_tile_y(ii) = jtile - !call mpas_log_write('---- Adding tile to local tile list ($i, land_x:$i, $i)', intArgs=(/local_tile_x(ii),local_tile_x_land(ii), local_tile_y(ii)/)) - !call mpas_log_write('-------- local box_x : $i local_box_x(l) - topo_shift: $i', intArgs=(/local_box_x(l),local_box_x(l) - topo_shift/)) - end if - end do - - cum_local_tiles = ii - - ! - end function determine_complete_tile_set - function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path) result(thisTile) @@ -553,6 +464,18 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path end function add_tile + !*********************************************************************** + ! + ! function read_global_30s_topo + ! + !> \brief Reads global 30-arc-second topography into 'topo' module variable + !> \author Michael Duda + !> \date 28 August 2017 + !> \details + !> This subroutine reads the global 30-arc-second topography from the subdirectory + !> identified by the 'sub_path' argument within the 'path' provided as the first argument. + ! + !----------------------------------------------------------------------- function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) result(iErr) implicit none @@ -623,7 +546,6 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re return end if - !topo(1:(tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) deallocate(tile) @@ -632,7 +554,18 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re end function read_30s_topo_tile - + !*********************************************************************** + ! + ! function read_global_30s_landuse + ! + !> \brief Reads global 30-arc-second landuse into 'landuse' module variable + !> \author Michael Duda + !> \date 14 March 2017 + !> \details + !> This subroutine reads the global 30-arc-second USGS landuse from + !> the subdirectory 'landuse_30s' of the path provided as an argument. + ! + !----------------------------------------------------------------------- function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start_y) result(iErr) implicit none @@ -689,165 +622,11 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start end function read_30s_landuse_tile - !*********************************************************************** - ! - ! function read_global_30s_topo - ! - !> \brief Reads global 30-arc-second topography into 'topo' module variable - !> \author Michael Duda - !> \date 28 August 2017 - !> \details - !> This subroutine reads the global 30-arc-second topography from the subdirectory - !> identified by the 'sub_path' argument within the 'path' provided as the first argument. - ! - !----------------------------------------------------------------------- - function read_global_30s_topo(path, sub_path) result(iErr) - - implicit none - - character(len=*), intent(in) :: path - character(len=*), intent(in) :: sub_path - - integer :: iErr - - 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 :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile - - integer (c_int) :: istatus - integer :: ix, iy, ix_shift, il - integer (c_int) :: isigned, endian, wordsize, nx, ny, nz - real (c_float) :: scalefactor - real (c_float), dimension(:,:,:), pointer, contiguous :: tile - type (c_ptr) :: tile_ptr - character(len=StrKIND) :: filename - character(len=StrKIND):: message - character(kind=c_char), dimension(StrKIND+1) :: c_filename - - allocate(tile(tile_x+2*tile_bdr,tile_y+2*tile_bdr,1)) - tile_ptr = c_loc(tile) - - isigned = 1 - endian = 0 - wordsize = 2 - scalefactor = 1.0 - nx = tile_x + 2*tile_bdr - ny = tile_y + 2*tile_bdr - nz = 1 - - topo_shift = 0 - - allocate(topo(tile_x*tile_y*cum_local_tiles)) - - ! - ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location - ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 - ! - if (trim(sub_path) == 'topo_gmted2010_30s/') then - topo_shift = topo_x / 2 - end if - - do il=1,cum_local_tiles - ix = local_tile_x(il) - iy = local_tile_y(il) - - 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) - call mpas_f_to_c_string(filename, c_filename) - call read_geogrid(c_filename, tile_ptr, nx, ny, nz, isigned, endian, & - wordsize, istatus) - tile(:,:,:) = tile(:,:,:) * scalefactor - if (istatus /= 0) then - call mpas_log_write('Error reading topography tile '//trim(filename), messageType=MPAS_LOG_ERR) - iErr = 1 - return - end if - - ix_shift = (il-1)*tile_x*tile_y + 1 - topo(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) - - end do - - deallocate(tile) - - iErr = 0 - - end function read_global_30s_topo - - - !*********************************************************************** - ! - ! function read_global_30s_landuse - ! - !> \brief Reads global 30-arc-second landuse into 'landuse' module variable - !> \author Michael Duda - !> \date 14 March 2017 - !> \details - !> This subroutine reads the global 30-arc-second USGS landuse from - !> the subdirectory 'landuse_30s' of the path provided as an argument. - ! - !----------------------------------------------------------------------- - function read_global_30s_landuse(path) result(iErr) - - implicit none - - character(len=*), intent(in) :: path - - integer :: iErr - - integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second landuse - integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second landuse - - integer (c_int) :: istatus - integer :: ix, iy, ix_shift, il - integer (c_int) :: isigned, endian, wordsize, nx, ny, nz - real (c_float) :: scalefactor - real (c_float), dimension(:,:,:), pointer, contiguous :: tile - type (c_ptr) :: tile_ptr - character(len=StrKIND) :: filename - character(kind=c_char), dimension(StrKIND+1) :: c_filename - character(len=StrKIND):: message - - allocate(tile(tile_x,tile_y,1)) - tile_ptr = c_loc(tile) - - isigned = 1 - endian = 0 - wordsize = 1 - scalefactor = 1.0 - nx = tile_x - ny = tile_y - nz = 1 - - allocate(landuse(tile_x*tile_y*cum_local_tiles)) - - do il=1,cum_local_tiles - ix = local_tile_x_land(il) - iy = local_tile_y(il) - - write(filename,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(path)//'/landuse_30s/', ix, '-', (ix+tile_x-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) - tile(:,:,:) = tile(:,:,:) * scalefactor - if (istatus /= 0) then - call mpas_log_write('Error reading landuse tile '//trim(filename)) - iErr = 1 - return - end if - - ix_shift = (il-1)*tile_x*tile_y + 1 - landuse(ix_shift:(ix_shift+tile_x*tile_y-1)) = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) - - end do - - deallocate(tile) - - iErr = 0 + - end function read_global_30s_landuse + + !*********************************************************************** ! @@ -887,8 +666,27 @@ integer (kind=I1KIND) function get_dom_landmask(nx, ny) end function get_dom_landmask - - subroutine get_box_points(lat, lon, nx, ny, path, sub_path, tilesHead) + !*********************************************************************** + ! + ! subroutine get_box + ! + !> \brief Cuts out a rectangular box of data centered at a given (lat,lon) + !> \author Michael Duda + !> \date 29 May 2015 + !> \details + !> This subroutine extracts a rectangular sub-array of the 30-arc-second + !> global topography dataset, stored in the module variable 'topo'; the + !> sub-array will be centered at the (lat,lon) specified in the argument + !> list, and will have a width and height large enough to span 'dx' meters. + !> The extracted sub-array is stored in the module variable 'box', and the + !> dimensions of this sub-array are stored in the module variables 'nx' and + !> 'ny'. + !> Since the mean value of the terrain in a grid cell is needed by many of + !> the GWDO statistics computations, this mean value is also computed by + !> this subroutine and stored in the module variable 'box_mean'. + ! + !----------------------------------------------------------------------- + subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) implicit none character(len=*), intent(in) :: path @@ -986,98 +784,12 @@ subroutine get_box_points(lat, lon, nx, ny, path, sub_path, tilesHead) box_mean = box_mean / real(nx*ny, RKIND) - end subroutine get_box_points - - - - !*********************************************************************** - ! - ! subroutine get_box - ! - !> \brief Cuts out a rectangular box of data centered at a given (lat,lon) - !> \author Michael Duda - !> \date 29 May 2015 - !> \details - !> This subroutine extracts a rectangular sub-array of the 30-arc-second - !> global topography dataset, stored in the module variable 'topo'; the - !> sub-array will be centered at the (lat,lon) specified in the argument - !> list, and will have a width and height large enough to span 'dx' meters. - !> The extracted sub-array is stored in the module variable 'box', and the - !> dimensions of this sub-array are stored in the module variables 'nx' and - !> 'ny'. - !> Since the mean value of the terrain in a grid cell is needed by many of - !> the GWDO statistics computations, this mean value is also computed by - !> this subroutine and stored in the module variable 'box_mean'. - ! - !----------------------------------------------------------------------- - subroutine get_box(nx, ny) - - implicit none - - integer, intent(in) :: nx, ny - integer :: itile, ix_shift, ix, ix_land, ix_topo, jx, ibox, tile_index_i, tile_index_j - integer :: tile_base, tile_offset, i, j - real (kind=RKIND) :: sg_lat - character(len=StrKIND):: message - - - if (associated(box)) deallocate(box) - allocate(box(nx,ny)) - - if (associated(box_landuse)) deallocate(box_landuse) - allocate(box_landuse(nx,ny)) - - if (associated(dxm)) deallocate(dxm) - allocate(dxm(nx,ny)) - - ! - ! Extract sub-array (box) from global array; must properly account for - ! the periodicity in the longitude coordinate, as well as the poles - ! - box_mean = 0.0 - ibox = 1 - do j=1,ny - do i=1,nx - - ibox = (j-1)*nx + i - ! pixels on the global topo grid - ix = local_box_x(ibox) !+ topo_shift - ix_topo = local_box_x(ibox) - jx = local_box_y(ibox) - - if (ix - topo_shift > 0 ) then - ix = ix - topo_shift - else - ix = ix + topo_shift - end if - - tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 - tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 - - itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) - - !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) - - tile_base = (itile-1)*tile_x*tile_y + 1 - tile_offset = (ix - local_tile_x(itile)) + (jx - local_tile_y(itile)) * tile_x - - box(i,j) = topo(tile_base + tile_offset) - !box_landuse(i,j) = landuse(tile_base + tile_offset) - sg_lat = (start_lat + (real(jx-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) - - ibox = ibox + 1 - end do - end do + end subroutine get_box - ! - ! Compute mean topography in the extracted box - ! - box_mean = box_mean / real(nx*ny, RKIND) - end subroutine get_box + + !*********************************************************************** From c48ddb76ec00e8498ea26f4fa118c9ab1e1c95ff Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 5 Sep 2024 18:28:07 +0000 Subject: [PATCH 13/27] cleanup and adding comments --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 175 ++++++++----------- 1 file changed, 72 insertions(+), 103 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index a2f45f4e6f..af9aadf4de 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -19,11 +19,15 @@ module mpas_init_atm_gwd private + ! 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 real (kind=R4KIND), dimension(:), pointer :: land_array - integer :: tile_start_x, tile_start_y + ! 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 ! linked list next pointer type (mpas_gwd_tile_type), pointer :: next => null() @@ -63,7 +67,8 @@ end subroutine read_geogrid real (kind=RKIND) :: start_lat real (kind=RKIND) :: start_lon - integer :: topo_shift ! special handling + ! To introduce an offset in the x-coordinate in the case of GMTED2010 topo data, as the dataset starts at 0.0 longitude + integer :: topo_shift = 0 ! Nominal delta-x (in meters) for sub-grid topography cells real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) @@ -71,7 +76,6 @@ end subroutine read_geogrid real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell real (kind=RKIND) :: box_mean ! Mean value of topography in box - !integer :: nx, ny ! Dimensions of box covering grid cell integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell ! NB: At present, only the USGS GLCC land cover dataset is supported, so we can assume 16 == water @@ -81,9 +85,11 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height - integer, dimension(:), pointer :: local_box_x ! - integer, dimension(:), pointer :: local_box_y ! - integer:: max_local_tiles, cum_local_tiles + ! For each point in the box, contains x and y coordinates of the global pixel + integer, dimension(:), pointer :: local_box_x + integer, dimension(:), pointer :: local_box_y + integer :: max_local_tiles ! Max possible tiles required by current rank + integer :: cum_local_tiles ! Total number of tiles required by current rank contains @@ -167,12 +173,14 @@ 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. start_lat = -90.0_RKIND start_lon = -180.0_RKIND + ! so we introduce an offset in the x-coordinate of topo_x/2 + topo_shift = topo_x / 2 case('default') call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) call mpas_log_write('Invalid topography dataset '''//trim(config_topo_data) & @@ -210,11 +218,9 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'gamma', hgamma) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) - nx = topo_x - ny = topo_y - - allocate(local_box_x(nx*ny)) - allocate(local_box_y(nx*ny)) + + allocate(local_box_x(topo_x*topo_y)) + allocate(local_box_y(topo_x*topo_y)) allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells @@ -242,7 +248,7 @@ function compute_gwd_fields(domain) result(iErr) max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) - topo_shift = topo_x / 2 + ! ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. @@ -337,6 +343,18 @@ function compute_gwd_fields(domain) result(iErr) end function compute_gwd_fields + + !*********************************************************************** + ! + ! subroutine get_box_size_from_lat_lon + ! + !> \brief Routine to obtain box size given dx, lat, lon (in radians) + !> \author Abishek Gopal + !> \date 05 Sep 2024 + !> \details + !> Routine to obtain box size given dx, lat, lon (in radians) + ! + !----------------------------------------------------------------------- subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) implicit none @@ -366,39 +384,26 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) end subroutine get_box_size_from_lat_lon - function locate_tile_in_list(tile_index_i, tile_index_j, is_topo) result (tile_index) - - implicit none - type(mpas_gwd_tile_type), pointer :: tilesHead - type(mpas_gwd_tile_type), pointer :: thisTile - integer, intent(in) :: tile_index_i - integer, intent(in) :: tile_index_j - logical, intent(in) :: is_topo - - integer :: itile, tile_index - ! loop over tiles - thisTile => tilesHead - - do while (associated(thisTile)) - - if(thisTile%tile_start_x==tile_index_i .and. thisTile%tile_start_y==tile_index_j) then - tile_index = itile - return - end if - - thisTile => thisTile % next - - enddo ! associated(thisTile) - - end function locate_tile_in_list - - - function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path) result(thisTile) + !*********************************************************************** + ! + ! function get_tile_from_box_point + ! + !> \brief Routine to obtain a tile given a box pixel + !> \author Abishek Gopal + !> \date 05 Sep 2024 + !> \details + !> Routine to obtain a tile of type mpas_gwd_tile_type, given the linked + ! list of tiles tilesHead, box coordinates box_x, box_y, and path to + ! static dataset. The function first searches the linked list to locate + ! the tile, and if search fails, adds the target tile to the linked list + ! after reading in the data for the tile from disk + !----------------------------------------------------------------------- + function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result(thisTile) implicit none integer :: iErr - integer, intent(in) :: box_x, box_y, nx, ny + 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 @@ -437,6 +442,18 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, nx, ny, path, sub_path end function get_tile_from_box_point + + !*********************************************************************** + ! + ! function get_tile_from_box_point + ! + !> \brief Routine to obtain box size given dx, lat, lon (in radians) + !> \author Abishek Gopal + !> \date 05 Sep 2024 + !> \details + !> Routine to obtain box size given dx, lat, lon (in radians) + ! + !----------------------------------------------------------------------- function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) result(newTile) implicit none @@ -466,13 +483,13 @@ end function add_tile !*********************************************************************** ! - ! function read_global_30s_topo + ! function read_30s_topo_tile ! - !> \brief Reads global 30-arc-second topography into 'topo' module variable + !> \brief Reads a single tile of the global 30-arc-second topography into memory !> \author Michael Duda !> \date 28 August 2017 !> \details - !> This subroutine reads the global 30-arc-second topography from the subdirectory + !> This subroutine reads a single tile of the global 30-arc-second topography from the subdirectory !> identified by the 'sub_path' argument within the 'path' provided as the first argument. ! !----------------------------------------------------------------------- @@ -512,26 +529,9 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re ny = tile_y + 2*tile_bdr nz = 1 - topo_shift = 0 - - ! if (tile_start_x - topo_shift > 0 ) then - ! ix = tile_start_x - topo_shift - ! else - ! ix = tile_start_x + topo_shift - ! end if iy = tile_start_y ix = tile_start_x - !allocate(topo(tile_x*tile_y)) - - ! - ! For GMTED2010 data, the dataset starts at 0.0 longitude, but we need to shift the starting location - ! in the topo array to -180.0, so we introduce an offset in the x-coordinate of topo_x/2 - ! - if (trim(sub_path) == 'topo_gmted2010_30s/') then - topo_shift = topo_x / 2 - end if - !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), '.', & @@ -556,14 +556,14 @@ end function read_30s_topo_tile !*********************************************************************** ! - ! function read_global_30s_landuse + ! function read_30s_landuse_tile ! - !> \brief Reads global 30-arc-second landuse into 'landuse' module variable + !> \brief Reads a single tile of the global 30-arc-second landuse into memory !> \author Michael Duda !> \date 14 March 2017 !> \details - !> This subroutine reads the global 30-arc-second USGS landuse from - !> the subdirectory 'landuse_30s' of the path provided as an argument. + !> This subroutine reads the a single tile of global 30-arc-second USGS landuse + !> from the subdirectory 'landuse_30s' of the path provided as an argument. ! !----------------------------------------------------------------------- function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start_y) result(iErr) @@ -622,12 +622,6 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start end function read_30s_landuse_tile - - - - - - !*********************************************************************** ! ! function get_dom_landmask @@ -671,8 +665,8 @@ end function get_dom_landmask ! subroutine get_box ! !> \brief Cuts out a rectangular box of data centered at a given (lat,lon) - !> \author Michael Duda - !> \date 29 May 2015 + !> \author Michael Duda, Abishek Gopal + !> \date Sep 2024 !> \details !> This subroutine extracts a rectangular sub-array of the 30-arc-second !> global topography dataset, stored in the module variable 'topo'; the @@ -747,32 +741,12 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) end do l = l + 1 + ! Obtain tile for given box pixel from the linked list of tiles (tilesHead), + ! which would involve reading in the data from disk if said tile is not already in memory + thisTile => get_tile_from_box_point(tilesHead, ii, jj, path, sub_path) - ! if (ii - topo_shift > 0 ) then - ! ii = ii - topo_shift - ! else - ! ii = ii + topo_shift - ! end if - - thisTile => get_tile_from_box_point(tilesHead, ii, jj, nx, ny, path, sub_path) - - ! ibox = (j-1)*nx + i - - - - ! tile_index_i = floor( real(ix - 1) / real(tile_x)) * tile_x + 1 - ! tile_index_j = floor( real(jx - 1) / real(tile_y)) * tile_y + 1 - - ! itile = locate_tile_in_set(tile_index_i, tile_index_j, .true.) - - !call mpas_log_write('In box_new for (i,j):($i,$i), (ix,jx):($i,$i), tile_index_i,tile_index_j:($i,$i), itile:$i ', intArgs=(/i,j,ix,jx,tile_index_i,tile_index_j,itile/)) - - ! tile_base = (itile-1)*tile_x*tile_y + 1 tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 - !if (tile_offset < 0) then - ! call mpas_log_write('In box_new for (i,j):($i,$i), tile_offset:$i tile_start_end $i, $i', intArgs=(/ii,jj,tile_offset, thisTile%tile_start_x, thisTile%tile_start_y/)) - ! call mpas_log_write('In box_new for (i,j):($i,$i), lat:$r, lon:$r', intArgs=(/ii,jj/), realArgs=(/lat*rad2deg, lon*rad2deg /)) - !end if + box(i,j) = thisTile%topo_array(tile_offset) box_landuse(i,j) = thisTile%land_array(tile_offset) sg_lat = (start_lat + (real(jj-1,RKIND) + 0.5) / pts_per_degree) / rad2deg ! Add 0.5 for cell center @@ -787,11 +761,6 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) end subroutine get_box - - - - - !*********************************************************************** ! ! function get_var From 979145e09b31a849d2e7ff932ff0c79e155cd9ff Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 5 Sep 2024 21:30:18 +0000 Subject: [PATCH 14/27] some more cleanup + function to dealloc tiles --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 93 ++++++++++++++------ 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index af9aadf4de..ffc3d17a98 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -140,7 +140,7 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - TYPE(mpas_gwd_tile_type), POINTER :: tilesHead + TYPE(mpas_gwd_tile_type), POINTER :: tilesHead ! Pointer to linked list of tiles ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell @@ -173,7 +173,6 @@ 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. @@ -218,11 +217,12 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'gamma', hgamma) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) - + allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells + + ! Allocate fields to hold all points in a box. Using a very conservative + ! box size estimate. allocate(local_box_x(topo_x*topo_y)) allocate(local_box_y(topo_x*topo_y)) - - allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells ! @@ -254,8 +254,6 @@ function compute_gwd_fields(domain) result(iErr) ! grid cell in the mesh. ! cum_local_tiles=0 - - call mpas_log_write(' End of determine_complete_tile_set. We need ($i / $i) tiles ', intArgs=(/cum_local_tiles, max_local_tiles/)) do iCell=1,nCells @@ -273,20 +271,19 @@ function compute_gwd_fields(domain) result(iErr) end if dc = dc * config_gwd_cell_scaling - call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) - ! ! Cut out a rectangular piece of the global 30-arc-second topography - ! data that is centered at the lat/lon of the current cell being + ! data that is centered at the lat/lon (in radians) of the current cell being ! processed and that is just large enough to cover the cell. The ! rectangular array of topography data is stored in the module - ! variable 'box', and the dimensions of this array are given by the - ! module variables 'nx' and 'ny'. The get_box() routine also + ! variable 'box', and the dimensions of this array are obtained from + ! the routine get_box_size_from_lat_lon and store in the + ! local variables 'nx' and 'ny'. The get_box() routine also ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! + call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) call get_box(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path, tilesHead) - !call get_box(nx, ny, geog_data_path, geog_sub_path) ! ! With a box of 30-arc-second data for the current grid cell, call @@ -339,6 +336,8 @@ function compute_gwd_fields(domain) result(iErr) deallocate(hlanduse) + iErr = remove_tiles(tilesHead) + iErr = 0 end function compute_gwd_fields @@ -348,11 +347,12 @@ end function compute_gwd_fields ! ! subroutine get_box_size_from_lat_lon ! - !> \brief Routine to obtain box size given dx, lat, lon (in radians) + !> \brief Routine to obtain box size given the mean diameter (meters), lat, lon (radians) !> \author Abishek Gopal !> \date 05 Sep 2024 !> \details - !> Routine to obtain box size given dx, lat, lon (in radians) + !> Routine to obtain box size (nx, ny) given the mean diameter of the grid cell (meters), + ! and the latitude and longitude coordinates (radians) ! !----------------------------------------------------------------------- subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) @@ -402,16 +402,17 @@ end subroutine get_box_size_from_lat_lon function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result(thisTile) implicit none - integer :: iErr 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 - integer :: ix, jx, ii, l,i, iCell, tile_start_x, tile_start_y, tile_start_x_topo + 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 + ! in certain datasets. We use tile_start_x, tile_start_y to search for tiles and open landmask tiles, + ! whereas tile_start_x_topo is only required to open the correct topo tiles from disk if (box_x > topo_shift) then tile_start_x_topo = floor( real(box_x - topo_shift - 1) / real(tile_x)) * tile_x + 1 else @@ -421,9 +422,9 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result tile_start_y = floor( real(box_y - 1) / real(tile_y)) * tile_y + 1 - ! loop over tiles + thisTile => tilesHead - + ! loop over tiles do while (associated(thisTile)) if(thisTile%tile_start_x==tile_start_x .and. thisTile%tile_start_y==tile_start_y) then @@ -434,6 +435,7 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result enddo ! associated(thisTile) + ! Couldn't find such a tile, so we add the tile to the front of the linked list if (.not. associated(thisTile)) then thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) end if @@ -445,14 +447,14 @@ end function get_tile_from_box_point !*********************************************************************** ! - ! function get_tile_from_box_point + ! function add_tile ! - !> \brief Routine to obtain box size given dx, lat, lon (in radians) + !> \brief Routine to read in a new topo and landmask tile !> \author Abishek Gopal !> \date 05 Sep 2024 !> \details - !> Routine to obtain box size given dx, lat, lon (in radians) - ! + !> Routine to read in a new topo and landmask tile, given the tile + ! coordinates !----------------------------------------------------------------------- function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) result(newTile) @@ -472,15 +474,56 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path newTile%next => tilesHead iErr = read_30s_topo_tile(path, sub_path, newTile%topo_array, tile_start_x_topo, newTile%tile_start_y) + if (iErr /= 0) then + call mpas_log_write('Error reading global 30-arc-sec topography for GWD statistics', messageType=MPAS_LOG_ERR) + return + end if + iErr = read_30s_landuse_tile(path, sub_path, newTile%land_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 + end if + tilesHead => newTile cum_local_tiles = cum_local_tiles + 1 - !call mpas_log_write('In add_tile for tile_start:($i,$i), cum_local_tiles:$i ', intArgs=(/tile_start_x, tile_start_y,cum_local_tiles /)) end function add_tile + + !*********************************************************************** + ! + ! function remove_tiles + ! + !> \brief Routine to deallocate all tiles in the list + !> \author Abishek Gopal + !> \date 05 Sep 2024 + !> \details + !> Routine to deallocate all tiles in the list + ! + !----------------------------------------------------------------------- + function remove_tiles(tilesHead) result(iErr) + + implicit none + integer :: iErr + type(mpas_gwd_tile_type), pointer :: tilesHead + type(mpas_gwd_tile_type), pointer :: thisTile + + ! loop over tiles + do while (associated(tilesHead)) + thisTile => tilesHead + tilesHead => thisTile % next + deallocate(thisTile) + deallocate(thisTile%topo_array) + deallocate(thisTile%land_array) + enddo ! associated(thisTile) + + iErr = 0 + + end function remove_tiles + !*********************************************************************** ! ! function read_30s_topo_tile From 7d2ccd8a72aa8105fea30445b63e42a04224b2b7 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 5 Sep 2024 21:47:09 +0000 Subject: [PATCH 15/27] fix --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index ffc3d17a98..0b7cd5c355 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -515,9 +515,9 @@ function remove_tiles(tilesHead) result(iErr) do while (associated(tilesHead)) thisTile => tilesHead tilesHead => thisTile % next - deallocate(thisTile) deallocate(thisTile%topo_array) - deallocate(thisTile%land_array) + deallocate(thisTile%land_array) + deallocate(thisTile) enddo ! associated(thisTile) iErr = 0 From 9577a5cce40800d778695b22e8d5e19e9270e8ea Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 6 Sep 2024 05:35:25 +0000 Subject: [PATCH 16/27] with lat,lon in deg to compare results --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0b7cd5c355..789699d5c0 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -282,8 +282,8 @@ function compute_gwd_fields(domain) result(iErr) ! computes the mean elevation in the array and stores that value in ! the module variable 'box_mean'. ! - call get_box_size_from_lat_lon(latCell(iCell), lonCell(iCell), dc, nx, ny) - call get_box(latCell(iCell),lonCell(iCell), nx, ny, geog_data_path, geog_sub_path, tilesHead) + 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) ! ! With a box of 30-arc-second data for the current grid cell, call @@ -369,8 +369,8 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) ! ! Get number of points to extract in the zonal direction ! - if (cos(lat) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat))) + if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) else nx = topo_x / 2 end if @@ -750,8 +750,8 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) ! ! Find coordinates in global topography array of the box center ! - ic = nint((lon*rad2deg - start_lon) * pts_per_degree) + 1 - jc = nint((lat*rad2deg - start_lat) * pts_per_degree) + 1 + ic = nint((lon - start_lon) * pts_per_degree) + 1 + jc = nint((lat - start_lat) * pts_per_degree) + 1 if (ic <= 0) ic = ic + topo_x if (ic > topo_x) ic = ic - topo_x From e3eb6e6e68bc434bb62f8b43278610f9f7ff23a2 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 12 Sep 2024 23:23:16 +0000 Subject: [PATCH 17/27] Addressing review comments. - box, box_landuse, dxm are now local variables - landuse changed to be of type I1KIND - cleanup --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 207 +++++++++++-------- 1 file changed, 117 insertions(+), 90 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 789699d5c0..9732d559db 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -19,11 +19,13 @@ module mpas_init_atm_gwd private + integer, parameter :: I1KIND = selected_int_kind(2) + ! 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 - real (kind=R4KIND), dimension(:), pointer :: land_array + integer (kind=I1KIND), dimension(:), pointer :: land_array ! 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. @@ -49,8 +51,6 @@ subroutine read_geogrid(fname, rarray, nx, ny, nz, isigned, endian, & end subroutine read_geogrid end interface - integer, parameter :: I1KIND = selected_int_kind(2) - 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 @@ -73,11 +73,6 @@ end subroutine read_geogrid ! Nominal delta-x (in meters) for sub-grid topography cells real (kind=RKIND), parameter :: sg_delta = 2.0 * Pi * Re / (360.0_RKIND * real(pts_per_degree,RKIND)) - real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell - real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell - real (kind=RKIND) :: box_mean ! Mean value of topography in box - integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell - ! NB: At present, only the USGS GLCC land cover dataset is supported, so we can assume 16 == water ! See the read_global_30s_landuse function integer (kind=I1KIND), parameter :: WATER = 16 @@ -88,9 +83,7 @@ end subroutine read_geogrid ! For each point in the box, contains x and y coordinates of the global pixel integer, dimension(:), pointer :: local_box_x integer, dimension(:), pointer :: local_box_y - integer :: max_local_tiles ! Max possible tiles required by current rank - integer :: cum_local_tiles ! Total number of tiles required by current rank - + contains @@ -140,15 +133,18 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - TYPE(mpas_gwd_tile_type), POINTER :: tilesHead ! Pointer to linked list of tiles + TYPE(mpas_gwd_tile_type), pointer :: tilesHead ! Pointer to linked list of tiles ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse real (kind=RKIND) :: sum_var - integer :: nx, ny - character(len=StrKIND):: message + real (kind=RKIND), dimension(:,:), pointer :: box => null() ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer :: dxm => null() ! Size (meters) in zonal direction of a grid cell + real (kind=RKIND) :: box_mean ! Mean value of topography in box + integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse => null() ! Subset of landuse covering a grid cell + integer :: nx, ny call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh) call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state) @@ -212,10 +208,10 @@ function compute_gwd_fields(domain) result(iErr) call mpas_pool_get_array(mesh, 'oa2', oa2) call mpas_pool_get_array(mesh, 'oa3', oa3) call mpas_pool_get_array(mesh, 'oa4', oa4) -! call mpas_pool_get_array(mesh, 'elvmax', elvmax) -! call mpas_pool_get_array(mesh, 'theta', htheta) -! call mpas_pool_get_array(mesh, 'gamma', hgamma) -! call mpas_pool_get_array(mesh, 'sigma', hsigma) + ! call mpas_pool_get_array(mesh, 'elvmax', elvmax) + ! call mpas_pool_get_array(mesh, 'theta', htheta) + ! call mpas_pool_get_array(mesh, 'gamma', hgamma) + ! call mpas_pool_get_array(mesh, 'sigma', hsigma) allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells @@ -246,14 +242,10 @@ function compute_gwd_fields(domain) result(iErr) tilesHead => null() - max_local_tiles = (topo_x/tile_x) * (topo_y/tile_y) - - ! ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. ! - cum_local_tiles=0 do iCell=1,nCells @@ -283,35 +275,37 @@ function compute_gwd_fields(domain) result(iErr) ! the module variable 'box_mean'. ! 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) + + call get_box(latCell(iCell)*rad2deg,lonCell(iCell)*rad2deg, nx, ny, & + & 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 ! subroutines to compute each sub-grid orography statistic ! - var2d(iCell) = get_var(nx, ny) - con(iCell) = get_con(nx, ny) - oa1(iCell) = get_oa1(nx, ny) - oa2(iCell) = get_oa2(nx, ny) - oa3(iCell) = get_oa3(nx, ny) - oa4(iCell) = get_oa4(nx, ny) + var2d(iCell) = get_var(box, box_mean, nx, ny) + con(iCell) = get_con(box, box_landuse, box_mean, nx, ny) + oa1(iCell) = get_oa1(box, box_mean, nx, ny) + oa2(iCell) = get_oa2(box, box_mean, nx, ny) + oa3(iCell) = get_oa3(box, box_mean, nx, ny) + oa4(iCell) = get_oa4(box, box_mean, nx, ny) ! Critical height, to be used in OL computation ! See Appendix of Kim, Y-J, 1996: Representation of Sub-Grid Scale Orographic Effects ! in a General Circulation Model. J. Climate, 9, 2698-2717. hc = 1116.2_RKIND - 0.878_RKIND * var2d(iCell) - ol1(iCell) = get_ol1(nx, ny) - ol2(iCell) = get_ol2(nx, ny) - ol3(iCell) = get_ol3(nx, ny) - ol4(iCell) = get_ol4(nx, ny) + ol1(iCell) = get_ol1(box, nx, ny) + ol2(iCell) = get_ol2(box, nx, ny) + ol3(iCell) = get_ol3(box, nx, ny) + ol4(iCell) = get_ol4(box, nx, ny) - hlanduse(iCell) = get_dom_landmask(nx, ny) ! get dominant land mask in cell + hlanduse(iCell) = get_dom_landmask(box_landuse, nx, ny) ! get dominant land mask in cell - ! elvmax_1 = get_elvmax(nx,ny) - ! htheta_1 = get_htheta(nx,ny) - ! hgamma_1 = get_hgamma(nx,ny) - ! hsigma_1 = get_hsigma(nx,ny) + ! elvmax(iCell) = get_elvmax(box, nx, ny) + ! htheta(iCell) = get_htheta(box, dxm, nx, ny) + ! hgamma(iCell) = get_hgamma(box, dxm, nx, ny) + ! hsigma(iCell) = get_hsigma(box, dxm, nx, ny) end do @@ -338,8 +332,6 @@ function compute_gwd_fields(domain) result(iErr) iErr = remove_tiles(tilesHead) - iErr = 0 - end function compute_gwd_fields @@ -362,7 +354,6 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) real (kind=RKIND), intent(in) :: lat real (kind=RKIND), intent(in) :: lon real (kind=RKIND), intent(in) :: dx - integer, intent(out) :: nx integer, intent(out) :: ny @@ -402,6 +393,7 @@ end subroutine get_box_size_from_lat_lon function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result(thisTile) implicit none + integer, intent(in) :: box_x, box_y type(mpas_gwd_tile_type), pointer, intent(in) :: tilesHead type(mpas_gwd_tile_type), pointer :: thisTile @@ -420,8 +412,6 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result end if tile_start_x = floor( real(box_x - 1) / real(tile_x)) * tile_x + 1 tile_start_y = floor( real(box_y - 1) / real(tile_y)) * tile_y + 1 - - thisTile => tilesHead ! loop over tiles @@ -433,14 +423,13 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result thisTile => thisTile % next - enddo ! associated(thisTile) + end do ! associated(thisTile) ! Couldn't find such a tile, so we add the tile to the front of the linked list if (.not. associated(thisTile)) then thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) end if - ! end function get_tile_from_box_point @@ -459,12 +448,13 @@ end function get_tile_from_box_point function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) result(newTile) implicit none + integer :: iErr integer, intent(in) :: tile_start_x, tile_start_y, tile_start_x_topo type(mpas_gwd_tile_type), pointer :: tilesHead - type(mpas_gwd_tile_type), pointer :: newTile character(len=*), intent(in) :: path character(len=*), intent(in) :: sub_path + type(mpas_gwd_tile_type), pointer :: newTile allocate(newTile) allocate(newTile%topo_array(topo_x*topo_y)) @@ -487,9 +477,6 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path tilesHead => newTile - cum_local_tiles = cum_local_tiles + 1 - - end function add_tile @@ -507,8 +494,10 @@ end function add_tile function remove_tiles(tilesHead) result(iErr) implicit none + integer :: iErr type(mpas_gwd_tile_type), pointer :: tilesHead + type(mpas_gwd_tile_type), pointer :: thisTile ! loop over tiles @@ -542,15 +531,14 @@ 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, intent(in) :: topo + real(kind=R4KIND), dimension(:), pointer :: topo integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y - integer :: iErr + integer :: iErr 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 :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile - integer (c_int) :: istatus integer :: ix, iy, ix_shift, il integer (c_int) :: isigned, endian, wordsize, nx, ny, nz @@ -558,7 +546,6 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re real (c_float), dimension(:,:,:), pointer, contiguous :: tile type (c_ptr) :: tile_ptr character(len=StrKIND) :: filename - character(len=StrKIND):: message character(kind=c_char), dimension(StrKIND+1) :: c_filename allocate(tile(tile_x+2*tile_bdr,tile_y+2*tile_bdr,1)) @@ -615,14 +602,11 @@ 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 - REAL(kind=R4KIND), DIMENSION(:), POINTER, intent(in) :: landuse + integer (kind=I1KIND), dimension(:), pointer :: landuse integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y + integer :: iErr - - integer, parameter :: tile_x = 1200 ! x-dimension of each tile of global 30-arc-second landuse - integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second landuse - integer (c_int) :: istatus integer :: ix, iy, ix_shift, il integer (c_int) :: isigned, endian, wordsize, nx, ny, nz @@ -631,7 +615,6 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start type (c_ptr) :: tile_ptr character(len=StrKIND) :: filename character(kind=c_char), dimension(StrKIND+1) :: c_filename - character(len=StrKIND):: message allocate(tile(tile_x,tile_y,1)) tile_ptr = c_loc(tile) @@ -656,7 +639,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start return end if - landuse = reshape(int(tile(1:tile_x,1:tile_y,1), kind=I1KIND),(/tile_x*tile_y/)) + landuse = int(reshape(tile,(/tile_x*tile_y/)), kind=I1KIND) deallocate(tile) @@ -675,9 +658,11 @@ end function read_30s_landuse_tile !> \details 1 = land, 0 = water ! !----------------------------------------------------------------------- - integer (kind=I1KIND) function get_dom_landmask(nx, ny) + integer (kind=I1KIND) function get_dom_landmask(box_landuse, nx, ny) implicit none + + integer (kind=I1KIND), dimension(:,:), pointer, intent(in) :: box_landuse ! Subset of landuse covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -723,21 +708,24 @@ end function get_dom_landmask !> this subroutine and stored in the module variable 'box_mean'. ! !----------------------------------------------------------------------- - subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) + subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse, dxm, box_mean) implicit none - character(len=*), intent(in) :: path - character(len=*), intent(in) :: sub_path - real (kind=RKIND), intent(in) :: lat, lon + real (kind=RKIND), intent(in) :: lat, lon integer, intent(in) :: nx, ny - TYPE(mpas_gwd_tile_type), POINTER, intent(in) :: tilesHead - TYPE(mpas_gwd_tile_type), POINTER :: thisTile - - integer :: i, j, ii, jj, ic, jc, l, tile_offset + character(len=*), intent(in) :: path + character(len=*), intent(in) :: sub_path + TYPE(mpas_gwd_tile_type), pointer, intent(in) :: 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 + real (kind=RKIND), intent(inout) :: box_mean ! Mean value of topography in box + + TYPE(mpas_gwd_tile_type), pointer :: thisTile + integer :: i, j, ii, jj, ic, jc, tile_offset real (kind=RKIND) :: sg_lat - if (associated(box)) deallocate(box) allocate(box(nx,ny)) @@ -761,7 +749,6 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) ! the periodicity in the longitude coordinate, as well as the poles ! box_mean = 0.0 - l = 0 do j=1,ny do i=1,nx @@ -782,7 +769,6 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) do while (ii > topo_x) ii = ii - topo_x end do - l = l + 1 ! Obtain tile for given box pixel from the linked list of tiles (tilesHead), ! which would involve reading in the data from disk if said tile is not already in memory @@ -795,10 +781,13 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead) 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 + ! box_mean = box_mean / real(nx*ny, RKIND) end subroutine get_box @@ -814,12 +803,15 @@ end subroutine get_box !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_var(nx, ny) + real (kind=RKIND) function get_var(box, box_mean, nx, ny) implicit none - integer :: i, j + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny + + integer :: i, j real (kind=RKIND) :: s2 s2 = 0.0 @@ -845,12 +837,16 @@ end function get_var !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_con(nx, ny) + real (kind=RKIND) function get_con(box, box_landuse, box_mean, nx, ny) implicit none - integer :: i, j + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + integer (kind=I1KIND), dimension(:,:), pointer, intent(in) :: box_landuse ! Subset of landuse covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny + + integer :: i, j real (kind=RKIND) :: s2, s4, var, xland, mean_land, mean_water, oro s2 = 0.0 @@ -923,12 +919,15 @@ end function get_con !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa1(nx, ny) + real (kind=RKIND) function get_oa1(box, box_mean, nx, ny) implicit none - integer :: i, j + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny + + integer :: i, j integer :: nu, nd nu = 0 @@ -963,9 +962,12 @@ end function get_oa1 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa2(nx, ny) + real (kind=RKIND) function get_oa2(box, box_mean, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny integer :: i, j @@ -1005,9 +1007,12 @@ end function get_oa2 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa3(nx, ny) + real (kind=RKIND) function get_oa3(box, box_mean, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny integer :: i, j @@ -1048,9 +1053,12 @@ end function get_oa3 !> the comment from N. Wood in the footnote of Kim and Doyle (QRJMS, 2005). ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_oa4(nx, ny) + real (kind=RKIND) function get_oa4(box, box_mean, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny integer :: i, j @@ -1089,9 +1097,11 @@ end function get_oa4 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol1(nx, ny) + real (kind=RKIND) function get_ol1(box, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1123,9 +1133,11 @@ end function get_ol1 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol2(nx, ny) + real (kind=RKIND) function get_ol2(box, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1157,9 +1169,11 @@ end function get_ol2 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol3(nx, ny) + real (kind=RKIND) function get_ol3(box, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1197,9 +1211,11 @@ end function get_ol3 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_ol4(nx, ny) + real (kind=RKIND) function get_ol4(box, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1237,9 +1253,11 @@ end function get_ol4 !> \details ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_elvmax(nx, ny) + real (kind=RKIND) function get_elvmax(box, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1267,9 +1285,12 @@ end function get_elvmax !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_htheta(nx, ny) + real (kind=RKIND) function get_htheta(box, dxm, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: dxm ! Size (meters) in zonal direction of a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1316,9 +1337,12 @@ end function get_htheta !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_hgamma(nx, ny) + real (kind=RKIND) function get_hgamma(box, dxm, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: dxm ! Size (meters) in zonal direction of a grid cell integer, intent(in) :: nx, ny integer :: i, j @@ -1370,9 +1394,12 @@ end function get_hgamma !> \details Computation following Lott and Miller (QJRMS 1997) ! !----------------------------------------------------------------------- - real (kind=RKIND) function get_hsigma(nx, ny) + real (kind=RKIND) function get_hsigma(box, dxm, nx, ny) implicit none + + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: dxm ! Size (meters) in zonal direction of a grid cell integer, intent(in) :: nx, ny integer :: i, j From 1bc486091add7e0e9c5a2aeb38a872e527d04978 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 13 Sep 2024 10:51:04 -0600 Subject: [PATCH 18/27] cleaning up unused variables --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 9732d559db..bfff8351e4 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -246,7 +246,6 @@ function compute_gwd_fields(domain) result(iErr) ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. ! - do iCell=1,nCells ! @@ -306,7 +305,6 @@ function compute_gwd_fields(domain) result(iErr) ! htheta(iCell) = get_htheta(box, dxm, nx, ny) ! hgamma(iCell) = get_hgamma(box, dxm, nx, ny) ! hsigma(iCell) = get_hsigma(box, dxm, nx, ny) - end do @@ -540,7 +538,7 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re integer, parameter :: tile_y = 1200 ! y-dimension of each tile of global 30-arc-second topography integer, parameter :: tile_bdr = 3 ! number of layers of border/halo points surrounding each tile integer (c_int) :: istatus - integer :: ix, iy, ix_shift, il + integer :: ix, iy integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile @@ -608,7 +606,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start integer :: iErr integer (c_int) :: istatus - integer :: ix, iy, ix_shift, il + integer :: ix, iy integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile From 9fc6ab6d4a0a2d8a908ef89004a82f855d1a6598 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 13 Sep 2024 20:06:58 +0000 Subject: [PATCH 19/27] - Change land and topo arrays to be 2D and fixing allocate size --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 39 ++++++++------------ 1 file changed, 15 insertions(+), 24 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index bfff8351e4..fff3ef63fa 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -24,8 +24,8 @@ 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 + integer (kind=I1KIND), dimension(:,:), pointer :: land_array ! 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. @@ -79,10 +79,6 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height - - ! For each point in the box, contains x and y coordinates of the global pixel - integer, dimension(:), pointer :: local_box_x - integer, dimension(:), pointer :: local_box_y contains @@ -214,12 +210,7 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells - - ! Allocate fields to hold all points in a box. Using a very conservative - ! box size estimate. - allocate(local_box_x(topo_x*topo_y)) - allocate(local_box_y(topo_x*topo_y)) - + ! ! It is possible that this code is called before the mesh fields have been scaled @@ -455,8 +446,8 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path type(mpas_gwd_tile_type), pointer :: newTile allocate(newTile) - allocate(newTile%topo_array(topo_x*topo_y)) - allocate(newTile%land_array(topo_x*topo_y)) + allocate(newTile%topo_array(tile_x,tile_y)) + allocate(newTile%land_array(tile_x,tile_y)) newTile%tile_start_x = tile_start_x newTile%tile_start_y = tile_start_y newTile%next => tilesHead @@ -529,7 +520,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 :: topo integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -574,7 +565,7 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re return end if - topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + topo = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) deallocate(tile) @@ -600,7 +591,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 :: landuse integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -637,7 +628,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start return end if - landuse = int(reshape(tile,(/tile_x*tile_y/)), kind=I1KIND) + landuse = int(tile(:,:,1), kind=I1KIND) deallocate(tile) @@ -721,7 +712,7 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse real (kind=RKIND), intent(inout) :: box_mean ! Mean value of topography in box TYPE(mpas_gwd_tile_type), pointer :: thisTile - integer :: i, j, ii, jj, ic, jc, tile_offset + integer :: i, j, ii, jj, ic, jc, ix, jx, tile_offset real (kind=RKIND) :: sg_lat if (associated(box)) deallocate(box) @@ -771,11 +762,11 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse ! Obtain tile for given box pixel from the linked list of tiles (tilesHead), ! which would involve reading in the data from disk if said tile is not already in memory thisTile => get_tile_from_box_point(tilesHead, ii, jj, path, sub_path) - - tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 - - box(i,j) = thisTile%topo_array(tile_offset) - box_landuse(i,j) = thisTile%land_array(tile_offset) + + 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) 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) From 47c990003f3b3988f08ef3233d7bfc8ceff1ecaf Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 13 Sep 2024 20:50:36 +0000 Subject: [PATCH 20/27] tmp --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 30 +++++++++++--------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index fff3ef63fa..c01ef2b591 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -24,8 +24,8 @@ 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 + integer (kind=I1KIND), dimension(:), pointer :: land_array ! 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. @@ -446,8 +446,8 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path type(mpas_gwd_tile_type), pointer :: newTile allocate(newTile) - allocate(newTile%topo_array(tile_x,tile_y)) - allocate(newTile%land_array(tile_x,tile_y)) + allocate(newTile%topo_array(tile_x*tile_y)) + allocate(newTile%land_array(tile_x*tile_y)) newTile%tile_start_x = tile_start_x newTile%tile_start_y = tile_start_y newTile%next => tilesHead @@ -520,7 +520,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 :: topo integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -565,8 +565,8 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re return end if - topo = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) - + !topo = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) + topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) deallocate(tile) iErr = 0 @@ -591,7 +591,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 :: landuse integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -628,7 +628,8 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start return end if - landuse = int(tile(:,:,1), kind=I1KIND) + !landuse = int(tile(:,:,1), kind=I1KIND) + landuse = int(reshape(tile,(/tile_x*tile_y/)), kind=I1KIND) deallocate(tile) @@ -763,10 +764,13 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse ! which would involve reading in the data from disk if said tile is not already in memory thisTile => get_tile_from_box_point(tilesHead, ii, jj, path, sub_path) - 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) + !ix = (ii - thisTile%tile_start_x) + 1 + !jx = (jj - thisTile%tile_start_y) + 1 + tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 + box(i,j) = thisTile%topo_array(tile_offset) + box_landuse(i,j) = thisTile%land_array(tile_offset) + !box(i,j) = thisTile%topo_array(ix, jx) + !box_landuse(i,j) = thisTile%land_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) From 010d89220b0eb68371be8ff47482e6a62d765eed Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Fri, 13 Sep 2024 21:04:35 +0000 Subject: [PATCH 21/27] reverting tmp changes --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 30 +++++++++----------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index c01ef2b591..fff3ef63fa 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -24,8 +24,8 @@ 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 + integer (kind=I1KIND), dimension(:,:), pointer :: land_array ! 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. @@ -446,8 +446,8 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path type(mpas_gwd_tile_type), pointer :: newTile allocate(newTile) - allocate(newTile%topo_array(tile_x*tile_y)) - allocate(newTile%land_array(tile_x*tile_y)) + allocate(newTile%topo_array(tile_x,tile_y)) + allocate(newTile%land_array(tile_x,tile_y)) newTile%tile_start_x = tile_start_x newTile%tile_start_y = tile_start_y newTile%next => tilesHead @@ -520,7 +520,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 :: topo integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -565,8 +565,8 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re return end if - !topo = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) - topo = reshape(tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1),(/tile_x*tile_y/)) + topo = tile((tile_bdr+1):(tile_x+tile_bdr),(tile_bdr+1):(tile_y+tile_bdr),1) + deallocate(tile) iErr = 0 @@ -591,7 +591,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 :: landuse integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y @@ -628,8 +628,7 @@ function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start return end if - !landuse = int(tile(:,:,1), kind=I1KIND) - landuse = int(reshape(tile,(/tile_x*tile_y/)), kind=I1KIND) + landuse = int(tile(:,:,1), kind=I1KIND) deallocate(tile) @@ -764,13 +763,10 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse ! which would involve reading in the data from disk if said tile is not already in memory thisTile => get_tile_from_box_point(tilesHead, ii, jj, path, sub_path) - !ix = (ii - thisTile%tile_start_x) + 1 - !jx = (jj - thisTile%tile_start_y) + 1 - tile_offset = (ii - thisTile%tile_start_x) + (jj - thisTile%tile_start_y) * tile_x + 1 - box(i,j) = thisTile%topo_array(tile_offset) - box_landuse(i,j) = thisTile%land_array(tile_offset) - !box(i,j) = thisTile%topo_array(ix, jx) - !box_landuse(i,j) = thisTile%land_array(ix, jx) + 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) 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) From c333251b68b3e403de73709241b0efa1e1031d98 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Thu, 19 Sep 2024 21:37:06 +0000 Subject: [PATCH 22/27] Addressing review comments --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 62 +++++++++++--------- 1 file changed, 34 insertions(+), 28 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index fff3ef63fa..c8413dda63 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -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 @@ -79,7 +79,7 @@ end subroutine read_geogrid integer (kind=I1KIND), dimension(:), pointer :: hlanduse ! Dominant land mask (0 or 1) real (kind=RKIND) :: hc ! critical height - + contains @@ -129,19 +129,24 @@ function compute_gwd_fields(domain) result(iErr) character(len=StrKIND) :: geog_sub_path character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash - TYPE(mpas_gwd_tile_type), pointer :: tilesHead ! Pointer to linked list of tiles + type(mpas_gwd_tile_type), pointer :: tilesHead ! Pointer to linked list of tiles ! Variables for smoothing variance integer, dimension(:,:), pointer:: cellsOnCell integer (kind=I1KIND) :: sum_landuse real (kind=RKIND) :: sum_var - real (kind=RKIND), dimension(:,:), pointer :: box => null() ! Subset of topography covering a grid cell - real (kind=RKIND), dimension(:,:), pointer :: dxm => null() ! Size (meters) in zonal direction of a grid cell - real (kind=RKIND) :: box_mean ! Mean value of topography in box - integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse => null() ! Subset of landuse covering a grid cell + real (kind=RKIND), dimension(:,:), pointer :: box ! Subset of topography covering a grid cell + real (kind=RKIND), dimension(:,:), pointer :: dxm ! Size (meters) in zonal direction of a grid cell + real (kind=RKIND) :: box_mean ! Mean value of topography in box + integer (kind=I1KIND), dimension(:,:), pointer :: box_landuse ! Subset of landuse covering a grid cell integer :: nx, ny + box => null() + dxm => null() + box_landuse => null() + tilesHead => null() + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh) call mpas_pool_get_subpool(domain % blocklist % structs, 'state', state) @@ -231,7 +236,6 @@ function compute_gwd_fields(domain) result(iErr) call mpas_log_write('in the computation of GWD static fields.') end if - tilesHead => null() ! ! Main loop to compute each of the GWDO fields for every horizontal @@ -259,10 +263,12 @@ function compute_gwd_fields(domain) result(iErr) ! processed and that is just large enough to cover the cell. The ! rectangular array of topography data is stored in the module ! variable 'box', and the dimensions of this array are obtained from - ! the routine get_box_size_from_lat_lon and store in the + ! the routine get_box_size_from_lat_lon() and stored in the ! local variables 'nx' and 'ny'. The get_box() routine also ! computes the mean elevation in the array and stores that value in - ! the module variable 'box_mean'. + ! the module variable 'box_mean'. 'tilesHead' points to the head of the linked + ! list of tiles, which is used by get_box() and its internal subroutines to search + ! for tile data and add new tiles to the head of this list as necessary. ! call get_box_size_from_lat_lon(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc, nx, ny) @@ -319,7 +325,7 @@ function compute_gwd_fields(domain) result(iErr) deallocate(hlanduse) - iErr = remove_tiles(tilesHead) + iErr = free_tile_list(tilesHead) end function compute_gwd_fields @@ -406,7 +412,7 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result ! loop over tiles do while (associated(thisTile)) - if(thisTile%tile_start_x==tile_start_x .and. thisTile%tile_start_y==tile_start_y) then + if (thisTile%tile_start_x==tile_start_x .and. thisTile%tile_start_y==tile_start_y) then exit end if @@ -414,7 +420,7 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result end do ! associated(thisTile) - ! Couldn't find such a tile, so we add the tile to the front of the linked list + ! Could not find such a tile, so we add the tile to the front of the linked list if (.not. associated(thisTile)) then thisTile => add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path, sub_path) end if @@ -471,7 +477,7 @@ end function add_tile !*********************************************************************** ! - ! function remove_tiles + ! function free_tile_list ! !> \brief Routine to deallocate all tiles in the list !> \author Abishek Gopal @@ -480,7 +486,7 @@ end function add_tile !> Routine to deallocate all tiles in the list ! !----------------------------------------------------------------------- - function remove_tiles(tilesHead) result(iErr) + function free_tile_list(tilesHead) result(iErr) implicit none @@ -500,7 +506,7 @@ function remove_tiles(tilesHead) result(iErr) iErr = 0 - end function remove_tiles + end function free_tile_list !*********************************************************************** ! @@ -705,13 +711,13 @@ 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 - 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 - real (kind=RKIND), intent(inout) :: box_mean ! Mean value of topography in box + type(mpas_gwd_tile_type), pointer, intent(in) :: 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 + real (kind=RKIND), intent(inout) :: box_mean ! Mean value of topography in box - TYPE(mpas_gwd_tile_type), pointer :: thisTile + type(mpas_gwd_tile_type), pointer :: thisTile integer :: i, j, ii, jj, ic, jc, ix, jx, tile_offset real (kind=RKIND) :: sg_lat @@ -830,8 +836,8 @@ real (kind=RKIND) function get_con(box, box_landuse, box_mean, nx, ny) implicit none - real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell - integer (kind=I1KIND), dimension(:,:), pointer, intent(in) :: box_landuse ! Subset of landuse covering a grid cell + real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell + integer (kind=I1KIND), dimension(:,:), pointer, intent(in) :: box_landuse ! Subset of landuse covering a grid cell real (kind=RKIND), intent(in) :: box_mean integer, intent(in) :: nx, ny From 53d9332ed0a5d0bb3d50b294ba37724b682847dd Mon Sep 17 00:00:00 2001 From: abishekg7 <56273301+abishekg7@users.noreply.github.com> Date: Thu, 19 Sep 2024 18:02:45 -0600 Subject: [PATCH 23/27] Removing whitespaces --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index c8413dda63..93a35f0bbd 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -215,7 +215,6 @@ function compute_gwd_fields(domain) result(iErr) ! call mpas_pool_get_array(mesh, 'sigma', hsigma) allocate(hlanduse(nCells+1)) ! +1, since we access hlanduse(cellsOnCell(i,iCell)) later on for iCell=1,nCells - ! ! It is possible that this code is called before the mesh fields have been scaled @@ -236,7 +235,6 @@ function compute_gwd_fields(domain) result(iErr) call mpas_log_write('in the computation of GWD static fields.') end if - ! ! Main loop to compute each of the GWDO fields for every horizontal ! grid cell in the mesh. From 0c4e0cf168cc7321cbf4ef41c4b8e18a3a659e4b Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 22 Oct 2024 16:45:11 -0600 Subject: [PATCH 24/27] Addressing review comments --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 66 +++++++++++--------- 1 file changed, 35 insertions(+), 31 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 93a35f0bbd..d3bb7569cf 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -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() @@ -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 @@ -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. @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -488,9 +492,10 @@ 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 @@ -498,9 +503,9 @@ function free_tile_list(tilesHead) result(iErr) 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 @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 @@ -770,7 +773,7 @@ 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) @@ -778,6 +781,7 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse end do end do + ! ! Compute mean topography in the extracted box ! From d0df0ebb5cdae0b72189114c2447dad9db384474 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 22 Oct 2024 17:00:53 -0600 Subject: [PATCH 25/27] Removing intent for return pointers. Removing whitespaces --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index d3bb7569cf..0ea9a12f7a 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -55,10 +55,10 @@ end subroutine read_geogrid 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 @@ -392,7 +392,7 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result character(len=*), intent(in) :: path character(len=*), intent(in) :: sub_path - type(mpas_gwd_tile_type), pointer, intent(out) :: thisTile + type(mpas_gwd_tile_type), pointer :: thisTile integer :: tile_start_x, tile_start_y, tile_start_x_topo @@ -449,7 +449,7 @@ function add_tile(tilesHead, tile_start_x, tile_start_y, tile_start_x_topo, path character(len=*), intent(in) :: path character(len=*), intent(in) :: sub_path - type(mpas_gwd_tile_type), pointer, intent(out) :: newTile + type(mpas_gwd_tile_type), pointer :: newTile integer :: iErr From 25ecf82d80ffbe048ad06d788db7b0e766d0a293 Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 5 Nov 2024 23:09:14 +0000 Subject: [PATCH 26/27] addressing review comments --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 26 +++++++++----------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index 0ea9a12f7a..ed4e4246d7 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -260,16 +260,16 @@ function compute_gwd_fields(domain) result(iErr) ! Cut out a rectangular piece of the global 30-arc-second topography ! data that is centered at the lat/lon (in radians) of the current cell being ! processed and that is just large enough to cover the cell. The - ! rectangular array of topography data is stored in the module + ! rectangular array of topography data is stored in the local ! variable 'box', and the dimensions of this array are obtained from - ! the routine get_box_size_from_lat_lon() and stored in the + ! the routine get_box_size_from_lat() and stored in the ! local variables 'nx' and 'ny'. The get_box() routine also ! computes the mean elevation in the array and stores that value in - ! the module variable 'box_mean'. 'tilesHead' points to the head of the linked + ! the local variable 'box_mean'. 'tilesHead' points to the head of the linked ! list of tiles, which is used by get_box() and its internal subroutines to search ! for tile data and add new tiles to the head of this list as necessary. ! - call get_box_size_from_lat_lon(latCell(iCell)*rad2deg, lonCell(iCell)*rad2deg, dc, nx, ny) + call get_box_size_from_lat(latCell(iCell), 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) @@ -338,15 +338,14 @@ end function compute_gwd_fields !> \date 05 Sep 2024 !> \details !> Routine to obtain box size (nx, ny) given the mean diameter of the grid cell (meters), - ! and the latitude and longitude coordinates (radians) + ! and the latitude (radians) ! !----------------------------------------------------------------------- - subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) + subroutine get_box_size_from_lat(lat, dx, nx, ny) implicit none real (kind=RKIND), intent(in) :: lat - real (kind=RKIND), intent(in) :: lon real (kind=RKIND), intent(in) :: dx integer, intent(out) :: nx integer, intent(out) :: ny @@ -354,8 +353,8 @@ subroutine get_box_size_from_lat_lon(lat, lon, dx, nx, ny) ! ! Get number of points to extract in the zonal direction ! - if (cos(lat/rad2deg) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then - nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat/rad2deg))) + if (cos(lat) > (2.0 * pts_per_degree * dx * 180.0) / (real(topo_x,RKIND) * Pi * Re)) then + nx = ceiling((180.0 * dx * pts_per_degree) / (Pi * Re * cos(lat))) else nx = topo_x / 2 end if @@ -365,8 +364,7 @@ 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 - + end subroutine get_box_size_from_lat !*********************************************************************** @@ -427,7 +425,6 @@ function get_tile_from_box_point(tilesHead, box_x, box_y, path, sub_path) result end function get_tile_from_box_point - !*********************************************************************** ! ! function add_tile @@ -511,6 +508,7 @@ function free_tile_list(tilesHead) result(iErr) end function free_tile_list + !*********************************************************************** ! ! function read_30s_topo_tile @@ -580,6 +578,7 @@ function read_30s_topo_tile(path, sub_path, topo, tile_start_x, tile_start_y) re end function read_30s_topo_tile + !*********************************************************************** ! ! function read_30s_landuse_tile @@ -592,12 +591,11 @@ end function read_30s_topo_tile !> from the subdirectory 'landuse_30s' of the path provided as an argument. ! !----------------------------------------------------------------------- - function read_30s_landuse_tile(path, sub_path, landuse, tile_start_x, tile_start_y) result(iErr) + function read_30s_landuse_tile(path, landuse, tile_start_x, tile_start_y) result(iErr) implicit none character(len=*), intent(in) :: path - character(len=*), intent(in) :: sub_path integer (kind=I1KIND), dimension(:,:), pointer, intent(inout) :: landuse integer, intent(in) :: tile_start_x integer, intent(in) :: tile_start_y From 650a60c94ce883434aaa177ce261ba8b286c673f Mon Sep 17 00:00:00 2001 From: Abishek Gopal Date: Tue, 5 Nov 2024 23:18:59 +0000 Subject: [PATCH 27/27] few more fixes --- src/core_init_atmosphere/mpas_init_atm_gwd.F | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_gwd.F b/src/core_init_atmosphere/mpas_init_atm_gwd.F index ed4e4246d7..8af9619a76 100644 --- a/src/core_init_atmosphere/mpas_init_atm_gwd.F +++ b/src/core_init_atmosphere/mpas_init_atm_gwd.F @@ -463,7 +463,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%landuse_array, newTile%tile_start_x, newTile%tile_start_y) + iErr = read_30s_landuse_tile(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 @@ -602,7 +602,6 @@ function read_30s_landuse_tile(path, landuse, tile_start_x, tile_start_y) result integer :: iErr integer (c_int) :: istatus - integer :: ix, iy integer (c_int) :: isigned, endian, wordsize, nx, ny, nz real (c_float) :: scalefactor real (c_float), dimension(:,:,:), pointer, contiguous :: tile @@ -717,7 +716,7 @@ subroutine get_box(lat, lon, nx, ny, path, sub_path, tilesHead, box, box_landuse real (kind=RKIND), intent(inout) :: box_mean ! Mean value of topography in box type(mpas_gwd_tile_type), pointer :: thisTile - integer :: i, j, ii, jj, ic, jc, ix, jx, tile_offset + integer :: i, j, ii, jj, ic, jc, ix, jx real (kind=RKIND) :: sg_lat if (associated(box)) deallocate(box)