topo_mod.f90.html | ![]() |
Source file: topo_mod.f90 | |
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/geoclaw/2d/lib | |
Converted: Mon Jan 21 2013 at 20:15:59 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
! ============================================================================ ! Program: /Users/mandli/src/claw44/branches/ktm-geoclaw-mod/2dxy/lib ! File: topo_mod ! Created: 2010-04-22 ! Author: Kyle Mandli ! ============================================================================ ! Copyright (C) 2010-04-22 Kyle Mandli! ! Distributed under the terms of the Berkeley Software Distribution (BSD) ! license ! http://www.opensource.org/licenses/ ! ============================================================================ ! Module for topography data ! ============================================================================ module topo_module implicit none ! Work array double precision, allocatable :: topowork(:) ! Topography file data character*150, allocatable :: topofname(:) integer :: mtopofiles,mtoposize double precision, allocatable :: xlowtopo(:), ylowtopo(:), tlowtopo(:) double precision, allocatable :: xhitopo(:), yhitopo(:), thitopo(:) double precision, allocatable :: dxtopo(:), dytopo(:) integer, allocatable :: mxtopo(:), mytopo(:) integer, allocatable :: i0topo(:), mtopo(:), mtopoorder(:) integer, allocatable :: minleveltopo(:), maxleveltopo(:), itopotype(:) ! Moving topography support integer :: imovetopo contains ! ======================================================================== ! Read topography files as specified in settopo.data ! ! Each topography file has a type stored in topotype(i). ! topotype = 1: standard GIS format: 3 columns: lon,lat,height(m) ! topotype = 2: Header as in DEM file, height(m) one value per line ! topotype = 3: Header as in DEM file, height(m) one row per line ! For other formats modify readtopo routine. ! ! advancing northwest to northeast then from north to south. Values should ! be uniformly spaced. ! ! Finest value of topography in a given region will be used for ! computation ! ======================================================================== subroutine set_topo(fname) use geoclaw_module implicit none ! Input arguments character*25, intent(in), optional :: fname ! Locals integer, parameter :: iunit = 7 integer :: i,j,itopo,finer_than,rank double precision :: area_i,area_j,x_junk,y_junk character*25 :: file_name logical :: found_file ! Open and begin parameter file output write(GEO_PARM_UNIT,*) ' ' write(GEO_PARM_UNIT,*) '--------------------------------------------' write(GEO_PARM_UNIT,*) 'SETTOPO:' write(GEO_PARM_UNIT,*) '---------' if (present(fname)) then file_name = fname else file_name = 'settopo.data' endif inquire(file=file_name,exist=found_file) if (.not. found_file) then print *, 'You must provide a file ', file_name stop endif call opendatafile(iunit, file_name) read(iunit,*) mtopofiles if (mtopofiles == 0) then write(GEO_PARM_UNIT,*) ' mtopofiles = 0' write(GEO_PARM_UNIT,*) ' No topo files specified, ' write(GEO_PARM_UNIT,*) ' will set B(x,y) = 0 in setaux' return endif write(GEO_PARM_UNIT,*) ' mtopofiles = ',mtopofiles ! Read and allocate data parameters for each file allocate(mxtopo(mtopofiles),mytopo(mtopofiles)) allocate(xlowtopo(mtopofiles),ylowtopo(mtopofiles)) allocate(tlowtopo(mtopofiles),xhitopo(mtopofiles),yhitopo(mtopofiles)) allocate(thitopo(mtopofiles),dxtopo(mtopofiles),dytopo(mtopofiles)) allocate(topofname(mtopofiles),itopotype(mtopofiles)) allocate(minleveltopo(mtopofiles),maxleveltopo(mtopofiles)) allocate(i0topo(mtopofiles),mtopo(mtopofiles),mtopoorder(mtopofiles)) do i=1,mtopofiles read(iunit,*) topofname(i) read(iunit,*) itopotype(i),minleveltopo(i), maxleveltopo(i), & tlowtopo(i),thitopo(i) write(GEO_PARM_UNIT,*) ' ' write(GEO_PARM_UNIT,*) ' ',topofname(i) write(GEO_PARM_UNIT,*) ' itopotype = ', itopotype(i) write(GEO_PARM_UNIT,*) ' minlevel, maxlevel = ', & minleveltopo(i), maxleveltopo(i) write(GEO_PARM_UNIT,*) ' tlow, thi = ', tlowtopo(i),thitopo(i) if (abs(itopotype(i)) == 1) then print *, 'WARNING: topotype 1 has been deprecated' print *, 'converting to topotype > 1 is encouraged' print *, 'python tools for converting files are provided' endif call read_topo_header(topofname(i),itopotype(i),mxtopo(i), & mytopo(i),xlowtopo(i),ylowtopo(i),xhitopo(i),yhitopo(i), & dxtopo(i),dytopo(i)) mtopo(i) = mxtopo(i)*mytopo(i) enddo ! Indexing into work array i0topo(1)=1 if (mtopofiles > 1) then do i=2,mtopofiles i0topo(i)=i0topo(i-1) + mtopo(i-1) enddo endif ! Read and allocate topography for each file mtoposize = sum(mtopo) allocate(topowork(mtoposize)) do i=1,mtopofiles call read_topo(mxtopo(i),mytopo(i),itopotype(i),topofname(i), & topowork(i0topo(i):i0topo(i)+mtopo(i)-1)) enddo ! topography order...This determines which order to process topography ! ! The finest topography will be given priority in any region ! mtopoorder(rank) = i means that i'th topography file has rank rank, ! where the file with rank=1 is the finest and considered first. do i=1,mtopofiles finer_than = 0 do j=1,mtopofiles if (j /= i) then area_i=dxtopo(i)*dytopo(i) area_j=dxtopo(j)*dytopo(j) if (area_i < area_j) finer_than = finer_than + 1 ! if two files have the same resolution, order is arbitrarily chosen if ((area_i == area_j).and.(j < i)) then finer_than = finer_than + 1 endif endif enddo ! # ifinerthan tells how many other files i is finer than rank = mtopofiles - finer_than mtopoorder(rank) = i enddo write(GEO_PARM_UNIT,*) ' ' write(GEO_PARM_UNIT,*) ' Ranking of topography files', & ' finest to coarsest: ', & (mtopoorder(rank),rank=1,mtopofiles) write(GEO_PARM_UNIT,*) ' ' end subroutine set_topo ! ======================================================================== ! read_topo(mx,my,dx,dy,xlow,xhi,ylow,yhi,itopo,fname,topo_type) ! ! Read topo file. ! New feature: topo_type < 0 means z values need to be negated. ! ======================================================================== subroutine read_topo(mx,my,topo_type,fname,topo) use geoclaw_module implicit none ! Arguments integer, intent(in) :: mx,my,topo_type character*150, intent(in) :: fname double precision, intent(inout) :: topo(1:mx*my) ! Locals integer, parameter :: iunit = 19, miss_unit = 17 double precision, parameter :: topo_missing = -150.d0 logical, parameter :: maketype2 = .false. integer :: i,j,num_points,missing,status,topo_start double precision :: no_data_value,x,y,z, topoi print *, ' ' print *, 'Reading topography file ', fname open(unit=iunit, file=fname, status='unknown',form='formatted') select case(abs(topo_type)) ! ASCII file with x,y,z values on each line. ! (progressing from upper left corner across rows, then down) ! Assumes a uniform rectangular grid of data values. case(1) i = 0 status = 0 do while (status == 0) i = i + 1 read(iunit,fmt=*,iostat=status) x,y,topoi if ((i .gt. mx*my) .and. (status == 0)) then write(6,*) '*** Error: i > mx*my = ',mx*my write(6,*) '*** i, mx, my: ',i,mx,my write(6,*) '*** status = ',status stop endif if (status == 0) then topo(i) = topoi endif enddo ! ================================================================ ! ASCII file with header followed by z data ! (progressing from upper left corner across rows, then down) ! one value per line if topo_type=2 or ! mx values per line if topo_type=3 ! ================================================================ case(2:3) ! Read header do i=1,5 read(iunit,*) enddo read(iunit,*) no_data_value ! Read in data missing = 0 select case(abs(topo_type)) case(2) do i=1,mx*my read(iunit,*) topo(i) if (topo(i) == no_data_value) then missing = missing + 1 topo(i) = topo_missing endif enddo case(3) do j=1,my read(iunit,*) (topo((j-1)*mx + i),i=1,mx) do i=1,mx if (topo((j-1)*mx + i) == no_data_value) then missing = missing + 1 topo((j-1)*mx + i) = topo_missing endif enddo enddo end select ! Write a warning if we found and missing values if (missing > 0) then print *, ' WARNING...some missing data values this file' print *, ' ',missing,' missing data values' print *, ' (see fort.missing)' print *, ' These values have arbitrarily been set to ',& topo_missing endif end select close(unit=iunit) ! Handle negative topo types if (topo_type < 0) then forall(i=1:mx*my) topo(i) = -topo(i) end forall endif ! ==================================================================== ! when topo_type=1 and data has x,y,z columns, ! set maketype2 to true to create a file new.topo_type2 with only z ! values, so next time it will take less time to read in. ! only works if dx = dy. ! if ((topo_type == 1).and.maketype2) then ! open(unit=29,file='new.tt2',status='unknown',form='formatted') ! write(29,*) mx, ' mx' ! write(29,*) my, ' my' ! write(29,*) xll, ' xllcorner' ! write(29,*) yll, ' yllcorner' ! write(29,*) dx, ' cellsize' ! write(29,*) -9999, ' nodataval' ! do i=topo_start,itopo ! write(29,'(d18.8)') topowork(i) ! enddo ! close(unit=29) ! print *, 'New topo file new.tt2 created' ! if (dx.ne.dy) then ! print *, ' ** Warning, dx must equal dy for this' ! print *, ' dx = ',dx,' dy = ',dy ! endif ! endif ! ! ==================================================================== end subroutine read_topo ! ======================================================================== ! subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy) ! ======================================================================== ! Read topo file header to determine space needed in allocatable array ! ! :Input: ! - fname - (char) Name of file ! - topo_type - (int) Type of topography file (-3 < topo_type < 3) ! ! :Output: ! - mx,my - (int) Number of grid points ! - xll,yll,xhi,yhi - (float) Lower and upper coordinates for grid ! - dx,dy - (float) Spatial resolution of grid ! ======================================================================== subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy) use geoclaw_module implicit none ! Input and Output character*150, intent(in) :: fname integer, intent(in) :: topo_type integer, intent(out) :: mx,my double precision, intent(out) :: xll,yll,xhi,yhi,dx,dy ! Local integer, parameter :: iunit = 19 integer :: topo_size, status double precision :: x,y,z,nodata_value logical :: found_file inquire(file=fname,exist=found_file) if (.not. found_file) then print *, 'Missing topography file:' print *, ' ', fname stop endif open(unit=iunit, file=fname, status='unknown',form='formatted') select case(abs(topo_type)) ! ASCII file with 3 columns ! determine data size case(1) ! Initial size variables topo_size = 0 mx = 0 ! Read in first values, determines xlow and yhi read(iunit,*) xll,yhi topo_size = topo_size + 1 mx = mx + 1 ! Go through first row figuring out mx, continue to count y = yhi do while (yhi == y) read(iunit,*) x,y,z topo_size = topo_size + 1 mx = mx + 1 enddo mx = mx - 1 ! Continue to count the rest of the lines status = 0 do while (status == 0) read(iunit,fmt=*,iostat=status) x,y,z topo_size = topo_size + 1 enddo if (status > 0) then print *, "IO error occured in ",fname,", aborting!" stop endif ! Calculate remaining values my = topo_size / mx xhi = x yll = y dx = (xhi-xll) / (mx-1) dy = (yhi-yll) / (my-1) ! ASCII file with header followed by z data case(2:3) read(iunit,*) mx read(iunit,*) my read(iunit,*) xll read(iunit,*) yll read(iunit,*) dx read(iunit,*) nodata_value dy = dx xhi = xll + (mx-1)*dx yhi = yll + (my-1)*dy case default print *, 'ERROR: Unrecognized topo_type' print *, ' topo_type = ',topo_type print *, ' for topography file:' print *, ' ', fname stop end select close(iunit) write(GEO_PARM_UNIT,*) ' mx = ',mx,' x = (',xll,',',xhi,')' write(GEO_PARM_UNIT,*) ' my = ',my,' y = (',yll,',',yhi,')' write(GEO_PARM_UNIT,*) ' dx, dy (meters/degrees) = ', dx,dy end subroutine read_topo_header end module topo_module