From e3e8f3278dfc315b00a7e79f42b03d570fa03784 Mon Sep 17 00:00:00 2001 From: David George Date: Sun, 27 May 2018 13:11:16 -0700 Subject: [PATCH 1/4] beginning to create integrate_module which contains raster_set data type. --- src/2d/shallow/integrate_module.f90 | 244 ++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 src/2d/shallow/integrate_module.f90 diff --git a/src/2d/shallow/integrate_module.f90 b/src/2d/shallow/integrate_module.f90 new file mode 100644 index 000000000..002b268fe --- /dev/null +++ b/src/2d/shallow/integrate_module.f90 @@ -0,0 +1,244 @@ +! ============================================================================ +! Program: integrate_module +! File: integrate_module.f90 +! Created: 2018-05-27 +! +! ============================================================================ +! +! +! Distributed under the terms of the Berkeley Software Distribution (BSD) +! license +! http://www.opensource.org/licenses/ +! this module contains integration routines that can fill grids by +! integrating multiple raster data sets (ie. sets of DEM file formats) +! It works for qinit, topo, and aux data types. +! ============================================================================ +! Module contains a few utilities for integrating. +! Defines a raster set data type. A single raster_set is a set of DEMs +! that define a particular field (eg. topo, component of q or aux) +! ============================================================================ +module integrate_module + + + implicit none + + + type raster_set + ! data for a set of rasters/DEMs for initializing grids (topo, qinit, aux) + ! Work array for all data for a particular field for all t + real(kind=8), allocatable :: rasterwork(:) + + !raster data + character(len=150), allocatable :: rasterfname(:) + integer :: mrasterfiles,mrastersize + real(kind=8), allocatable :: xlowraster(:), ylowraster(:), tlowraster(:) + real(kind=8), allocatable :: xhiraster(:), yhiraster(:), thiraster(:) + real(kind=8), allocatable :: dxraster(:), dyraster(:) + real(kind=8), allocatable :: rastertime(:) + integer, allocatable :: mxraster(:), myraster(:) + + integer, allocatable :: i0raster(:), mraster(:), mrasterorder(:) + integer, allocatable :: minlevelraster(:), maxlevelraster(:), irasterfiletype(:) + integer, allocatable :: rasterID(:),raster0save(:),irasterfield(:) + logical :: raster_finalized +contains + + +recursive subroutine topoarea(x1,x2,y1,y2,m,area) + + ! Compute the area of overlap of topo with the rectangle (x1,x2) x (y1,y2) + ! using topo arrays indexed mtopoorder(mtopofiles) through mtopoorder(m) + ! (coarse to fine). + + ! The main call to this subroutine has corners of a physical domain for + ! the rectangle and m = 1 in order to compute the area of overlap of + ! domain by all topo arrays. Used to check inputs and insure topo + ! covers domain. + + ! The recursive strategy is to first compute the area using only topo + ! arrays mtopoorder(mtopofiles) to mtopoorder(m+1), + ! and then apply corrections due to adding topo array mtopoorder(m). + + ! Corrections are needed if the new topo array intersects the grid cell. + ! Let the intersection be (x1m,x2m) x (y1m,y2m). + ! Two corrections are needed, first to subtract out the area over + ! the rectangle (x1m,x2m) x (y1m,y2m) computed using + ! topo arrays mtopoorder(mtopofiles) to mtopoorder(m+1), + ! and then adding in the area over this same region using + ! topo array mtopoorder(m). + + ! Based on the recursive routine rectintegral that integrates + ! topo over grid cells using a similar strategy. + + implicit none + + ! arguments + real (kind=8), intent(in) :: x1,x2,y1,y2 + integer, intent(in) :: m + real (kind=8), intent(out) :: area + + ! local + real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m, & + y1m,y2m, area1,area2,area_m + integer :: mfid, indicator, i0 + real(kind=8), external :: topointegral + + + mfid = mtopoorder(m) + i0=i0topo(mfid) + + if (m == mtopofiles) then + ! innermost step of recursion reaches this point. + ! only using coarsest topo grid -- compute directly... + call intersection(indicator,area,xmlo,xmhi, & + ymlo,ymhi, x1,x2,y1,y2, & + xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + + else + ! recursive call to compute area using one fewer topo grids: + call topoarea(x1,x2,y1,y2,m+1,area1) + + ! region of intersection of cell with new topo grid: + call intersection(indicator,area_m,x1m,x2m, & + y1m,y2m, x1,x2,y1,y2, & + xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + + + if (area_m > 0) then + + ! correction to subtract out from previous set of topo grids: + call topoarea(x1m,x2m,y1m,y2m,m+1,area2) + + ! adjust integral due to corrections for new topo grid: + area = area1 - area2 + area_m + else + area = area1 + endif + endif + +end subroutine topoarea + + + +recursive subroutine rectintegral(x1,x2,y1,y2,m,integral) + + ! Compute the integral of topo over the rectangle (x1,x2) x (y1,y2) + ! using topo arrays indexed mtopoorder(mtopofiles) through mtopoorder(m) + ! (coarse to fine). + + ! The main call to this subroutine has corners of a grid cell for the + ! rectangle and m = 1 in order to compute the integral over the cell + ! using all topo arrays. + + ! The recursive strategy is to first compute the integral using only topo + ! arrays mtopoorder(mtopofiles) to mtopoorder(m+1), + ! and then apply corrections due to adding topo array mtopoorder(m). + + ! Corrections are needed if the new topo array intersects the grid cell. + ! Let the intersection be (x1m,x2m) x (y1m,y2m). + ! Two corrections are needed, first to subtract out the integral over + ! the rectangle (x1m,x2m) x (y1m,y2m) computed using + ! topo arrays mtopoorder(mtopofiles) to mtopoorder(m+1), + ! and then adding in the integral over this same region using + ! topo array mtopoorder(m). + + ! Note that the function topointegral returns the integral over the + ! rectangle based on a single topo array, and that routine calls + ! bilinearintegral. + + + implicit none + + ! arguments + real (kind=8), intent(in) :: x1,x2,y1,y2 + integer, intent(in) :: m + real (kind=8), intent(out) :: integral + + ! local + real(kind=8) :: xmlo,xmhi,ymlo,ymhi,area,x1m,x2m, & + y1m,y2m, int1,int2,int3 + integer :: mfid, indicator, mp1fid, i0 + real(kind=8), external :: topointegral + + + mfid = mtopoorder(m) + i0=i0topo(mfid) + + if (m == mtopofiles) then + ! innermost step of recursion reaches this point. + ! only using coarsest topo grid -- compute directly... + call intersection(indicator,area,xmlo,xmhi, & + ymlo,ymhi, x1,x2,y1,y2, & + xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + + if (indicator.eq.1) then + ! cell overlaps the file + ! integrate surface over intersection of grid and cell + integral = topointegral( xmlo,xmhi,ymlo, & + ymhi,xlowtopo(mfid),ylowtopo(mfid),dxtopo(mfid), & + dytopo(mfid),mxtopo(mfid),mytopo(mfid),topowork(i0),1) + else + integral = 0.d0 + endif + + else + ! recursive call to compute area using one fewer topo grids: + call rectintegral(x1,x2,y1,y2,m+1,int1) + + ! region of intersection of cell with new topo grid: + call intersection(indicator,area,x1m,x2m, & + y1m,y2m, x1,x2,y1,y2, & + xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + + + if (area > 0) then + + ! correction to subtract out from previous set of topo grids: + call rectintegral(x1m,x2m,y1m,y2m,m+1,int2) + + ! correction to add in for new topo grid: + int3 = topointegral(x1m,x2m, y1m,y2m, & + xlowtopo(mfid),ylowtopo(mfid),dxtopo(mfid), & + dytopo(mfid),mxtopo(mfid),mytopo(mfid),topowork(i0),1) + + ! adjust integral due to corrections for new topo grid: + integral = int1 - int2 + int3 + else + integral = int1 + endif + endif + +end subroutine rectintegral + + +subroutine intersection(indicator,area,xintlo,xinthi, & + yintlo,yinthi,x1lo,x1hi,y1lo,y1hi,x2lo,x2hi,y2lo,y2hi) + + ! find the intersection of two rectangles, return the intersection + ! and it's area, and indicator =1 + ! if there is no intersection, indicator =0 + + implicit none + + integer, intent(out) :: indicator + + real(kind=8), intent(in) :: x1lo,x1hi,y1lo,y1hi,x2lo,x2hi,y2lo,y2hi + real(kind=8), intent(out) :: area,xintlo,xinthi,yintlo,yinthi + + xintlo=dmax1(x1lo,x2lo) + xinthi=dmin1(x1hi,x2hi) + yintlo=dmax1(y1lo,y2lo) + yinthi=dmin1(y1hi,y2hi) + + + if (xinthi.gt.xintlo.and.yinthi.gt.yintlo) then + area = (xinthi-xintlo)*(yinthi-yintlo) + indicator = 1 + else + area = 0.d0 + indicator = 0 + endif + +end subroutine intersection + +end module topo_module From d79469c30abd8d22004133e4abbd327147f4aaac Mon Sep 17 00:00:00 2001 From: David George Date: Sun, 27 May 2018 13:13:10 -0700 Subject: [PATCH 2/4] edits to new integrate_module which contains raster_set data type. --- src/2d/shallow/integrate_module.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/2d/shallow/integrate_module.f90 b/src/2d/shallow/integrate_module.f90 index 002b268fe..e002bbc14 100644 --- a/src/2d/shallow/integrate_module.f90 +++ b/src/2d/shallow/integrate_module.f90 @@ -14,8 +14,8 @@ ! It works for qinit, topo, and aux data types. ! ============================================================================ ! Module contains a few utilities for integrating. -! Defines a raster set data type. A single raster_set is a set of DEMs -! that define a particular field (eg. topo, component of q or aux) +! Defines a raster set data type. A single raster_set is a set of +! rasters/DEMs that define a particular field (eg. topo, component of q or aux) ! ============================================================================ module integrate_module From 31367c23960920ca66c71604c7672793df41b809 Mon Sep 17 00:00:00 2001 From: David George Date: Sun, 27 May 2018 21:58:06 -0700 Subject: [PATCH 3/4] integrate_module still wip. --- src/2d/shallow/integrate_module.f90 | 310 +++++++++++++++++++++++----- 1 file changed, 255 insertions(+), 55 deletions(-) diff --git a/src/2d/shallow/integrate_module.f90 b/src/2d/shallow/integrate_module.f90 index e002bbc14..dbd1b0096 100644 --- a/src/2d/shallow/integrate_module.f90 +++ b/src/2d/shallow/integrate_module.f90 @@ -14,8 +14,9 @@ ! It works for qinit, topo, and aux data types. ! ============================================================================ ! Module contains a few utilities for integrating. -! Defines a raster set data type. A single raster_set is a set of +! Defines raster data types. A single raster_collection is a set of ! rasters/DEMs that define a particular field (eg. topo, component of q or aux) +! a single raster_type is raster data input (eg. a topo file) ! ============================================================================ module integrate_module @@ -23,30 +24,239 @@ module integrate_module implicit none - type raster_set - ! data for a set of rasters/DEMs for initializing grids (topo, qinit, aux) - ! Work array for all data for a particular field for all t - real(kind=8), allocatable :: rasterwork(:) - - !raster data - character(len=150), allocatable :: rasterfname(:) - integer :: mrasterfiles,mrastersize - real(kind=8), allocatable :: xlowraster(:), ylowraster(:), tlowraster(:) - real(kind=8), allocatable :: xhiraster(:), yhiraster(:), thiraster(:) - real(kind=8), allocatable :: dxraster(:), dyraster(:) - real(kind=8), allocatable :: rastertime(:) - integer, allocatable :: mxraster(:), myraster(:) - - integer, allocatable :: i0raster(:), mraster(:), mrasterorder(:) - integer, allocatable :: minlevelraster(:), maxlevelraster(:), irasterfiletype(:) - integer, allocatable :: rasterID(:),raster0save(:),irasterfield(:) - logical :: raster_finalized + type raster_data_type + ! data for a set of rasters/DEMs for initializing grids (topo, qinit, aux) + ! Work array for all data for a particular field for all t + real(kind=8), allocatable :: data(:,:) + + !raster data + character(len=150), :: fname + + real(kind=8) :: xlow, ylow, tlow, xhi, yhi, thi + real(kind=8) :: dx, dy + real(kind=8) :: rastertime + + integer :: mx, my, morder + integer :: minlevel, maxlevel + integer :: ifiletype, ID, save0data, iq, iaux, ifieldID + logical :: finalized ! for this raster grid + end type raster_data_type + + type raster_collection_type + ! data for a set of rasters/DEMs for initializing grids (topo, qinit, aux) + ! All arrays for all data for a particular field for all t + type(raster_data_type), allocatable :: rasters(:) + integer, allocatable :: raster_IDs(:) + integer, allocatable :: raster_resolution_order(:) + + integer, :: num_raster_sets + integer, :: iq, iaux, ifieldID + logical :: finalized ! for all raster grids in collection + end type raster_collection_type + contains + subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell, + & yjp,raster_collection) + + implicit none + + ! input + real(kind=8), intent(in) :: xim,xcell,xip,yjm,ycell,yjp + type(raster_collection_type), intent(in) :: raster_collection + ! output + real(kind=8), intent(out) :: cellint + + !local + integer :: mrasters,rasternum + type(raster_data_type) :: raster_data + + !############################################################################## + ! cellrasterintegrate integrates a unique surface, + ! defined from data from a raster collection, + ! (using the finest data available in any region) + ! over a cell + ! + ! The rectangle has coords: + ! xim <= x <= xip, yjm <= y <= yjp, + ! with center (x,y) = (xcell, ycell) + + !The intersection (with one particular raster data) has coords: + !xintlo <= x <= xinthi, yintlo <= y <= yinthi + + !initialize the integral of the surface + cellint=0.d0 + + !determine the type of integration + im = 1 + + ! first see if the grid cell is entirely in a fine topofile + ! if so we are done...if not use recursive strategy for grids + !as coarse or coarser than first finest intersection + + mrasters = raster_collection%num_raster_sets + do m = 1,mrasters + !look at raster files, from fine to coarse + rasternum = raster_collection%raster_resolution_order(m)) + raster_data = raster_collection(rasternum) + !check for intersection of grid cell and this raster set + cellarea = (xip-xim)*(yjp-yjm) + call intersection(indicator,area,xmlo,xmhi,ymlo,ymhi,xim,xip,yjm,yjp,raster_data) + if (indicator.eq.1) then !cell overlaps raster set + if (area.eq.cellarea) then !cell is entirely in raster set + ! (should we check if they agree to some tolerance??) + !integrate surface and get out of here + cellint = cellint + rasterintegral(xmlo,xmhi,ymlo,ymhi,raster_data,im) + return + else + go to 222 + endif + endif + enddo + + + 222 continue + + ! this grid cell intersects only raster set m and perhaps coarser: + call rectintegral(xim,xip,yjm,yjp,m,cellint,raster_collection) + + return + end + + end subroutine cellgridintegrate + + !------------------------------------------------------------------------------- + + function rasterintegral(xim,xip,yjm,yjp,intmethod,raster_data) + !=========================================================================== + + !###################################################################### + ! rasterintegral integrates a surface over a cell + ! that is the intersection with a single raster data set + ! the surface integrated is defined by a piecewise bilinear through the + ! nodes of the raster data set + + ! The rectangular intersection has coords: + ! xim <= x <= xip, yjm <= y <= yjp + ! + ! The Cartesian grid has coords: + ! xxlow <= x <= xxhi, yylow <= y <= yyhi, with grid cell size dxx by dyy + ! and mxx by myy cells. + ! + ! written by David L. George + ! Seattle, WA 7/16/08 + !########################################################################### + + implicit none + + ! input + real(kind=8), intent(in) :: xim,xip,yjm,yjp + integer, intent(in) :: intmethod + type(raster_data_type) :: raster_data -recursive subroutine topoarea(x1,x2,y1,y2,m,area) + + ! output + real(kind=8), intent(out) :: cellint + + ! local + real(kind=8) :: theintegral,xxhi,xxlow,yyhi,yylow,dxx,dyy,dx,dy + real(kind=8) :: z11,z22,z12,z21 + integer, :: mxx,myy + integer, :: jj,jjstart,jjend + integer, :: ii,iistar,iiend + integer, :: djjstart,djjend,diistart,diiend + integer, :: jjz1,jjz2 + + xxlow = raster_data%xlow + yylow = raster_data%ylow + dxx = raster_data%dx + dyy = raster_data%dy + mxx = raster_data%mx + myy = raster_data%my + + !initialize: + theintegral = 0.d0 + + xxhi=xxlow+(mxx-1)*dxx + yyhi=yylow+(myy-1)*dyy + + !TEST FOR SMALL ROUNDING ERROR========== + if ((xim-xxlow).lt.0.d0.or.(xip-xxhi).gt.0.d0) then + xim=max(xxlow,xim) + xip=min(xxhi,xip) + endif + + if ((yjm-yylow).lt.0.d0.or.(yjp-yyhi).gt.0.d0) then + yjm=max(yylow,yjm) + yjp=min(yyhi,yjp) + endif + + dx=xip-xim + dy=yjp-yjm + + !INTEGRATE PIECEWISE BILINEAR OVER RECTANGULAR REGION==== + if (intmethod.eq.1) then !use bilinear method + + !find indices that include the rectangular region + djjstart=(yjm-yylow)/dyy + jjstart=idint(djjstart)+1 + + diistart=(xim-xxlow)/dxx + iistart=idint(diistart)+1 + + diiend=(xip-xxlow)/dxx + iiend=ceiling(diiend) + 1 + + djjend=(yjp-yylow)/dyy + jjend=ceiling(djjend)+1 + + iistart=max(iistart,1) + jjstart=max(jjstart,1) + iiend=min(mxx,iiend) + jjend=min(myy,jjend) + + + do jj=jjstart,jjend-1 + y1=yylow + (jj-1.d0)*dyy + y2=yylow + (jj)*dyy + !the array zz is indexed from north to south: + ! jjz is the actual index of interest in the array zz + jjz1= myy-jj+1 + jjz2= jjz1-1 + + do ii=iistart,iiend-1 + x1=xxlow + (ii-1.d0)*dxx + x2=xxlow + (ii)*dxx + + z11 = raster_data%data(ii,jjz1) + z12 = raster_data%data(ii,jjz2) + z21 = raster_data%data(ii+1,jjz1) + z22 = raster_data%data(ii+1,jjz2) + + if (coordinate_system.eq.1) then !cartesian rectangle + theintegral = theintegral + & + bilinearintegral(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) + elseif (coordinate_system.eq.2) then !integrate on surface of sphere + theintegral = theintegral + & + bilinearintegral_s(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) + else + write(*,*) 'TOPOINTEGRAL: coordinate_system error' + endif + enddo + enddo - ! Compute the area of overlap of topo with the rectangle (x1,x2) x (y1,y2) + else + write(*,*) 'RASTERINTEGRAL: only intmethod = 1,2 is supported' + endif + + topointegral= theintegral + return + end function rasterintegral + + +recursive subroutine topoarea(x1,x2,y1,y2,m,area,raster_collection) + + ! Compute the area of overlap of top with the rectangle (x1,x2) x (y1,y2) ! using topo arrays indexed mtopoorder(mtopofiles) through mtopoorder(m) ! (coarse to fine). @@ -77,37 +287,31 @@ recursive subroutine topoarea(x1,x2,y1,y2,m,area) integer, intent(in) :: m real (kind=8), intent(out) :: area - ! local - real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m, & - y1m,y2m, area1,area2,area_m - integer :: mfid, indicator, i0 - real(kind=8), external :: topointegral - + type(raster_collection_type), intent(in) :: raster_collection - mfid = mtopoorder(m) - i0=i0topo(mfid) + ! local + real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m,y1m,y2m,area1,area2,area_m + integer :: indicator - if (m == mtopofiles) then - ! innermost step of recursion reaches this point. - ! only using coarsest topo grid -- compute directly... - call intersection(indicator,area,xmlo,xmhi, & - ymlo,ymhi, x1,x2,y1,y2, & - xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + mrasters = raster_collection%num_raster_sets + rasternum = raster_collection%raster_resolution_order(m) + raster_data = raster_collection(rasternum) + if (m == mrasters) then + ! innermost step of recursion reaches this point. + ! only using coarsest topo grid -- compute directly... + call intersection(indicator,area,xmlo,xmhi,ymlo,ymhi,x1,x2,y1,y2,raster_data) else ! recursive call to compute area using one fewer topo grids: - call topoarea(x1,x2,y1,y2,m+1,area1) + call topoarea(x1,x2,y1,y2,m+1,area1,raster_collection) ! region of intersection of cell with new topo grid: - call intersection(indicator,area_m,x1m,x2m, & - y1m,y2m, x1,x2,y1,y2, & - xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) - - + call intersection(indicator,area_m,x1m,x2m,y1m,y2m,x1,x2,y1,y2,raster_data) + if (area_m > 0) then ! correction to subtract out from previous set of topo grids: - call topoarea(x1m,x2m,y1m,y2m,m+1,area2) + call topoarea(x1m,x2m,y1m,y2m,m+1,area2,raster_collection) ! adjust integral due to corrections for new topo grid: area = area1 - area2 + area_m @@ -120,7 +324,7 @@ end subroutine topoarea -recursive subroutine rectintegral(x1,x2,y1,y2,m,integral) +recursive subroutine rectintegral(x1,x2,y1,y2,m,integral,raster_collection) ! Compute the integral of topo over the rectangle (x1,x2) x (y1,y2) ! using topo arrays indexed mtopoorder(mtopofiles) through mtopoorder(m) @@ -155,21 +359,17 @@ recursive subroutine rectintegral(x1,x2,y1,y2,m,integral) real (kind=8), intent(out) :: integral ! local - real(kind=8) :: xmlo,xmhi,ymlo,ymhi,area,x1m,x2m, & - y1m,y2m, int1,int2,int3 - integer :: mfid, indicator, mp1fid, i0 - real(kind=8), external :: topointegral + real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m,y1m,y2m,int1,int2,int3,area + integer :: indicator + mrasters = raster_collection%num_raster_sets + rasternum = raster_collection%raster_resolution_order(m) + raster_data = raster_collection(rasternum) - mfid = mtopoorder(m) - i0=i0topo(mfid) - - if (m == mtopofiles) then + if (m == mrasters) then ! innermost step of recursion reaches this point. ! only using coarsest topo grid -- compute directly... - call intersection(indicator,area,xmlo,xmhi, & - ymlo,ymhi, x1,x2,y1,y2, & - xlowtopo(mfid),xhitopo(mfid),ylowtopo(mfid),yhitopo(mfid)) + call intersection(indicator,area,xmlo,xmhi,ymlo,ymhi,x1,x2,y1,y2,raster_data) if (indicator.eq.1) then ! cell overlaps the file @@ -241,4 +441,4 @@ subroutine intersection(indicator,area,xintlo,xinthi, & end subroutine intersection -end module topo_module +end module integrate_module From a22fdb93f788c112123cbbe16a940f4f22ad73eb Mon Sep 17 00:00:00 2001 From: David George Date: Sun, 27 May 2018 22:36:18 -0700 Subject: [PATCH 4/4] integrate_module still WIP. --- src/2d/shallow/integrate_module.f90 | 299 ++++++++++++++-------------- 1 file changed, 146 insertions(+), 153 deletions(-) diff --git a/src/2d/shallow/integrate_module.f90 b/src/2d/shallow/integrate_module.f90 index dbd1b0096..cc3962a68 100644 --- a/src/2d/shallow/integrate_module.f90 +++ b/src/2d/shallow/integrate_module.f90 @@ -20,13 +20,10 @@ ! ============================================================================ module integrate_module - implicit none - type raster_data_type - ! data for a set of rasters/DEMs for initializing grids (topo, qinit, aux) - ! Work array for all data for a particular field for all t + ! data for a single raster/DEM for initializing grids (topo, qinit, aux) real(kind=8), allocatable :: data(:,:) !raster data @@ -49,18 +46,18 @@ module integrate_module integer, allocatable :: raster_IDs(:) integer, allocatable :: raster_resolution_order(:) - integer, :: num_raster_sets + integer, :: num_rasters integer, :: iq, iaux, ifieldID logical :: finalized ! for all raster grids in collection end type raster_collection_type contains - subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell, - & yjp,raster_collection) - implicit none + !------------------------------------------------------------------------------------ + subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell,yjp,raster_collection) + implicit none ! input real(kind=8), intent(in) :: xim,xcell,xip,yjm,ycell,yjp type(raster_collection_type), intent(in) :: raster_collection @@ -94,7 +91,7 @@ subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell, ! if so we are done...if not use recursive strategy for grids !as coarse or coarser than first finest intersection - mrasters = raster_collection%num_raster_sets + mrasters = raster_collection%num_rasters do m = 1,mrasters !look at raster files, from fine to coarse rasternum = raster_collection%raster_resolution_order(m)) @@ -109,15 +106,13 @@ subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell, cellint = cellint + rasterintegral(xmlo,xmhi,ymlo,ymhi,raster_data,im) return else - go to 222 + exit endif endif enddo - - 222 continue - ! this grid cell intersects only raster set m and perhaps coarser: + ! consider only rasters > m in resolution order call rectintegral(xim,xip,yjm,yjp,m,cellint,raster_collection) return @@ -125,151 +120,22 @@ subroutine cellrasterintegrate(cellint,xim,xcell,xip,yjm,ycell, end subroutine cellgridintegrate - !------------------------------------------------------------------------------- - - function rasterintegral(xim,xip,yjm,yjp,intmethod,raster_data) - !=========================================================================== - - !###################################################################### - ! rasterintegral integrates a surface over a cell - ! that is the intersection with a single raster data set - ! the surface integrated is defined by a piecewise bilinear through the - ! nodes of the raster data set - - ! The rectangular intersection has coords: - ! xim <= x <= xip, yjm <= y <= yjp - ! - ! The Cartesian grid has coords: - ! xxlow <= x <= xxhi, yylow <= y <= yyhi, with grid cell size dxx by dyy - ! and mxx by myy cells. - ! - ! written by David L. George - ! Seattle, WA 7/16/08 - !########################################################################### - - implicit none - - ! input - real(kind=8), intent(in) :: xim,xip,yjm,yjp - integer, intent(in) :: intmethod - type(raster_data_type) :: raster_data - - - ! output - real(kind=8), intent(out) :: cellint - - ! local - real(kind=8) :: theintegral,xxhi,xxlow,yyhi,yylow,dxx,dyy,dx,dy - real(kind=8) :: z11,z22,z12,z21 - integer, :: mxx,myy - integer, :: jj,jjstart,jjend - integer, :: ii,iistar,iiend - integer, :: djjstart,djjend,diistart,diiend - integer, :: jjz1,jjz2 - - xxlow = raster_data%xlow - yylow = raster_data%ylow - dxx = raster_data%dx - dyy = raster_data%dy - mxx = raster_data%mx - myy = raster_data%my - - !initialize: - theintegral = 0.d0 - - xxhi=xxlow+(mxx-1)*dxx - yyhi=yylow+(myy-1)*dyy - - !TEST FOR SMALL ROUNDING ERROR========== - if ((xim-xxlow).lt.0.d0.or.(xip-xxhi).gt.0.d0) then - xim=max(xxlow,xim) - xip=min(xxhi,xip) - endif - - if ((yjm-yylow).lt.0.d0.or.(yjp-yyhi).gt.0.d0) then - yjm=max(yylow,yjm) - yjp=min(yyhi,yjp) - endif - - dx=xip-xim - dy=yjp-yjm - - !INTEGRATE PIECEWISE BILINEAR OVER RECTANGULAR REGION==== - if (intmethod.eq.1) then !use bilinear method - - !find indices that include the rectangular region - djjstart=(yjm-yylow)/dyy - jjstart=idint(djjstart)+1 - - diistart=(xim-xxlow)/dxx - iistart=idint(diistart)+1 - - diiend=(xip-xxlow)/dxx - iiend=ceiling(diiend) + 1 - - djjend=(yjp-yylow)/dyy - jjend=ceiling(djjend)+1 - - iistart=max(iistart,1) - jjstart=max(jjstart,1) - iiend=min(mxx,iiend) - jjend=min(myy,jjend) - - - do jj=jjstart,jjend-1 - y1=yylow + (jj-1.d0)*dyy - y2=yylow + (jj)*dyy - !the array zz is indexed from north to south: - ! jjz is the actual index of interest in the array zz - jjz1= myy-jj+1 - jjz2= jjz1-1 - - do ii=iistart,iiend-1 - x1=xxlow + (ii-1.d0)*dxx - x2=xxlow + (ii)*dxx - - z11 = raster_data%data(ii,jjz1) - z12 = raster_data%data(ii,jjz2) - z21 = raster_data%data(ii+1,jjz1) - z22 = raster_data%data(ii+1,jjz2) - - if (coordinate_system.eq.1) then !cartesian rectangle - theintegral = theintegral + & - bilinearintegral(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) - elseif (coordinate_system.eq.2) then !integrate on surface of sphere - theintegral = theintegral + & - bilinearintegral_s(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) - else - write(*,*) 'TOPOINTEGRAL: coordinate_system error' - endif - enddo - enddo - - else - write(*,*) 'RASTERINTEGRAL: only intmethod = 1,2 is supported' - endif - - topointegral= theintegral - return - end function rasterintegral - - -recursive subroutine topoarea(x1,x2,y1,y2,m,area,raster_collection) + !---------------------------------------------------------------------------- + recursive subroutine topoarea(x1,x2,y1,y2,m,area,raster_collection) - ! Compute the area of overlap of top with the rectangle (x1,x2) x (y1,y2) - ! using topo arrays indexed mtopoorder(mtopofiles) through mtopoorder(m) - ! (coarse to fine). + ! Compute the area of overlap of rasters with the rectangle (x1,x2) x (y1,y2) + ! using rasters indexed numrasters down to m. (coarse to fine). ! The main call to this subroutine has corners of a physical domain for ! the rectangle and m = 1 in order to compute the area of overlap of ! domain by all topo arrays. Used to check inputs and insure topo ! covers domain. - ! The recursive strategy is to first compute the area using only topo - ! arrays mtopoorder(mtopofiles) to mtopoorder(m+1), - ! and then apply corrections due to adding topo array mtopoorder(m). + ! The recursive strategy is to first compute the area using only rasters + ! numrasters to m + 1, + ! and then apply corrections due to adding raster m. - ! Corrections are needed if the new topo array intersects the grid cell. + ! Corrections are needed if the new raster array intersects the grid cell. ! Let the intersection be (x1m,x2m) x (y1m,y2m). ! Two corrections are needed, first to subtract out the area over ! the rectangle (x1m,x2m) x (y1m,y2m) computed using @@ -293,7 +159,7 @@ recursive subroutine topoarea(x1,x2,y1,y2,m,area,raster_collection) real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m,y1m,y2m,area1,area2,area_m integer :: indicator - mrasters = raster_collection%num_raster_sets + mrasters = raster_collection%num_rasters rasternum = raster_collection%raster_resolution_order(m) raster_data = raster_collection(rasternum) @@ -322,7 +188,6 @@ recursive subroutine topoarea(x1,x2,y1,y2,m,area,raster_collection) end subroutine topoarea - recursive subroutine rectintegral(x1,x2,y1,y2,m,integral,raster_collection) @@ -362,7 +227,7 @@ recursive subroutine rectintegral(x1,x2,y1,y2,m,integral,raster_collection) real(kind=8) :: xmlo,xmhi,ymlo,ymhi,x1m,x2m,y1m,y2m,int1,int2,int3,area integer :: indicator - mrasters = raster_collection%num_raster_sets + mrasters = raster_collection%num_rasters rasternum = raster_collection%raster_resolution_order(m) raster_data = raster_collection(rasternum) @@ -410,6 +275,134 @@ recursive subroutine rectintegral(x1,x2,y1,y2,m,integral,raster_collection) end subroutine rectintegral + !------------------------------------------------------------------------------- + + function rasterintegral(xim,xip,yjm,yjp,intmethod,raster_data) + !=========================================================================== + + !###################################################################### + ! rasterintegral integrates a surface over a cell + ! that is the intersection with a single raster data set + ! the surface integrated is defined by a piecewise bilinear through the + ! nodes of the raster data set + + ! The rectangular intersection has coords: + ! xim <= x <= xip, yjm <= y <= yjp + ! + ! The Cartesian grid has coords: + ! xxlow <= x <= xxhi, yylow <= y <= yyhi, with grid cell size dxx by dyy + ! and mxx by myy cells. + ! + ! written by David L. George + ! Seattle, WA 7/16/08 + !########################################################################### + + implicit none + + ! input + real(kind=8), intent(in) :: xim,xip,yjm,yjp + integer, intent(in) :: intmethod + type(raster_data_type) :: raster_data + + + ! output + real(kind=8), intent(out) :: cellint + + ! local + real(kind=8) :: theintegral,xxhi,xxlow,yyhi,yylow,dxx,dyy,dx,dy + real(kind=8) :: z11,z22,z12,z21 + integer, :: mxx,myy + integer, :: jj,jjstart,jjend + integer, :: ii,iistar,iiend + integer, :: djjstart,djjend,diistart,diiend + integer, :: jjz1,jjz2 + + xxlow = raster_data%xlow + yylow = raster_data%ylow + dxx = raster_data%dx + dyy = raster_data%dy + mxx = raster_data%mx + myy = raster_data%my + + !initialize: + theintegral = 0.d0 + + xxhi=xxlow+(mxx-1)*dxx + yyhi=yylow+(myy-1)*dyy + + !TEST FOR SMALL ROUNDING ERROR========== + if ((xim-xxlow).lt.0.d0.or.(xip-xxhi).gt.0.d0) then + xim=max(xxlow,xim) + xip=min(xxhi,xip) + endif + + if ((yjm-yylow).lt.0.d0.or.(yjp-yyhi).gt.0.d0) then + yjm=max(yylow,yjm) + yjp=min(yyhi,yjp) + endif + + dx=xip-xim + dy=yjp-yjm + + !INTEGRATE PIECEWISE BILINEAR OVER RECTANGULAR REGION==== + if (intmethod.eq.1) then !use bilinear method + + !find indices that include the rectangular region + djjstart=(yjm-yylow)/dyy + jjstart=idint(djjstart)+1 + + diistart=(xim-xxlow)/dxx + iistart=idint(diistart)+1 + + diiend=(xip-xxlow)/dxx + iiend=ceiling(diiend) + 1 + + djjend=(yjp-yylow)/dyy + jjend=ceiling(djjend)+1 + + iistart=max(iistart,1) + jjstart=max(jjstart,1) + iiend=min(mxx,iiend) + jjend=min(myy,jjend) + + + do jj=jjstart,jjend-1 + y1=yylow + (jj-1.d0)*dyy + y2=yylow + (jj)*dyy + !the array zz is indexed from north to south: + ! jjz is the actual index of interest in the array zz + jjz1= myy-jj+1 + jjz2= jjz1-1 + + do ii=iistart,iiend-1 + x1=xxlow + (ii-1.d0)*dxx + x2=xxlow + (ii)*dxx + + z11 = raster_data%data(ii,jjz1) + z12 = raster_data%data(ii,jjz2) + z21 = raster_data%data(ii+1,jjz1) + z22 = raster_data%data(ii+1,jjz2) + + if (coordinate_system.eq.1) then !cartesian rectangle + theintegral = theintegral + & + bilinearintegral(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) + elseif (coordinate_system.eq.2) then !integrate on surface of sphere + theintegral = theintegral + & + bilinearintegral_s(xim,xip,yjm,yjp,x1,x2,y1,y2,dxx,dyy,z11,z12,z21,z22) + else + write(*,*) 'TOPOINTEGRAL: coordinate_system error' + endif + enddo + enddo + + else + write(*,*) 'RASTERINTEGRAL: only intmethod = 1,2 is supported' + endif + + topointegral= theintegral + return + end function rasterintegral + subroutine intersection(indicator,area,xintlo,xinthi, & yintlo,yinthi,x1lo,x1hi,y1lo,y1hi,x2lo,x2hi,y2lo,y2hi)