PROGRAM arpstrn,18
!
!
!-----------------------------------------------------------------------
!
! Program for creating a terrain file for ARPS from
! USGS 3 arc second data base (covers continetal US and Alaska)
! or 5min global data base from NCAR.
!
! Author: Ming Xue, CAPS. University of Oklahma.
!
! First written: 11/25/1995.
! Last Modified: 12/05/1995.
!
! 11/1/1998 (Ming Xue).
! Significantly rewritten to also handle 5min global data.
!
! The code for accessing USGS FTP server and reading the DEM data
! patches were written by Burt Freeman and P. Patnaik of S-Cubed.
!
! We assume you have ZXPLOT Version 3.0 or later installed.
! Otherwise, link nozxplot.f.
!
! 2001/06/18 (Gene Bassett)
! Cleaned up namelist input of grid and map projections parameters.
!
!-----------------------------------------------------------------------
!
!
! INSTRUCTIONS:
!
! A number functions are controled by Makfile.
!
! To compile, do 'make arpstrn' in ARPS root directory.
!
! You can also do the compilation and linking directly in this
! directory, making use of pre-built arps library:
!
! zxpostf90 -I../../include -o arpstrn arpstrn.f -L../../lib -larps
!
! Here we assumed libarps.a has been produced by 'makearps arps'.
! We also assume you have ZXPLOT Version 3.0 or later installed.
! If not, do
!
! f90 -I../../include -o arpstrn arpstrn.f nozxplot.f -L../../lib -larps
!
! then no plot of the analyzed terrain field will be produced.
!
! To execute, do 'arpstrn < arpstrn.input >! arpstrn.output'
!
! When you use the 3 arc second data, a temporary directory
! whose name is specified for dir_trndata in the input file
! will be used to store the original data. The required
! data will be downloaded from anonymous ftp server
! edcftp.cr.usgs.gov unless you already have them in that directory.
!
! To allow automatic login, you need to add
!
! machine edcftp.cr.usgs.gov login anonymous password your_address
!
! in your $HOME/.netrc file. Make .netrc readable to the owner only
! i.e., do 'chmod 600 .netrc'.
!
! To produce a tar file of the entire package, do 'make arpstrn.tar'.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INTEGER :: nx, ny
REAL, allocatable :: x(:) ! The x-coord. of the u grid point (m)
REAL, allocatable :: y(:) ! The y-coord. of the v grid point (m)
REAL, allocatable :: xh(:,:) ! X-coord. of the terrain (scalar) point.
REAL, allocatable :: yh(:,:) ! Y-coord. of the terrain (scalar) point.
REAL, allocatable :: xh1d(:)
REAL, allocatable :: yh1d(:)
REAL, allocatable :: lath(:,:) ! Latitude of each terrain point
REAL, allocatable :: lonh(:,:) ! Longitude of each terrain point
REAL, allocatable :: hterain(:,:) ! The height of the terrain.
REAL, allocatable :: tem (:,:)
REAL, allocatable :: tem2 (:,:)
REAL, allocatable :: xw (:)
REAL, allocatable :: yw (:)
INTEGER :: nx_trn,ny_trn
REAL, allocatable :: h_trn(:,:)
NAMELIST /grid_dims/ nx,ny
NAMELIST /jobname/ runname
NAMELIST /grid/ dx,dy,ctrlat,ctrlon
NAMELIST /projection/ mapproj,trulat1,trulat2,trulon,sclfct
INTEGER :: trnanxopt ! Terrain analysis scheme option
INTEGER :: trndataopt ! Terrain data set option
INTEGER :: nsmth ! Number of times 9-point smoother is applied
! to the terrain field
CHARACTER (LEN=120) :: usgs_dem_index_fn
CHARACTER (LEN=80) :: fn_trndata
CHARACTER (LEN=120) :: dir_trndata
REAL :: hctrinc
INTEGER :: lat_sample,lon_sample ! Data sampling frequency in seconds
REAL :: latd_sample,lond_sample ! Data sampling frequency in degrees
NAMELIST /dem_trn/ trnanxopt,trndataopt,nsmth, &
lat_sample,lon_sample,dir_trndata,fn_trndata, &
usgs_dem_index_fn,hctrinc
INTEGER :: ldir_trndata,lfn
REAL :: latgrid,longrid
INTEGER :: maxmap
PARAMETER(maxmap=10)
CHARACTER (LEN=132) :: mapfile(maxmap)
INTEGER :: lmapfile,nmapfile
NAMELIST /map_plot/latgrid,longrid,nmapfile,mapfile
COMMON /mappar2/ latgrid,longrid
NAMELIST /trn_output/ terndmp
INTEGER :: trndatares ! Terrain data resolution in seconds
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
REAL :: trulat (2) ! Latitude at which the map projection
! is true (degees N)
REAL :: swx,swy,ctrx,ctry ! Temporary variables.
REAL :: latmax,latmin,lonmax,lonmin ! max/min lat/lon of the projected grid
INTEGER :: swlat_trn,swlon_trn,nelat_trn,nelon_trn
! SW and NE coner lat/lon of the
! the assembled terrain patch.
INTEGER :: i,j,ih,jh
REAL :: lath1,lonh1,w1,w2,w3,w4,a,b,c,d,h1,h2,h3
REAL :: badvalue
REAL :: hmax,hmin,swlat_offset,swlon_offset
INTEGER :: inith_set, badvfound
LOGICAL :: fexist
INTEGER :: imin,imax,jmin,jmax,kmin,kmax, istatus
!
!-----------------------------------------------------------------------
!
! TEMPORARY ARRAYS:
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Set up the default values for all the variables to be read in
! using the namelist method. In case the user does not specify a
! particular value, this value will be used.
!
!-----------------------------------------------------------------------
!
nx = 67
ny = 67
runname = 'may20'
dx = 1000.0
dy = 1000.0
ctrlat = 35.0
ctrlon = -100.0
mapproj = 0
trulat1 = 30.0
trulat2 = 60.0
trulon = -100.0
sclfct = 1.0
trnanxopt = 1
trndataopt = 4
fn_trndata = 'tbase_global_5min.data'
dir_trndata= '/work/official/arpstern.data/'
usgs_dem_index_fn = '../src/arpstrn/usgs_dem.index'
nsmth = 1
lat_sample = 30
lon_sample = 30
hctrinc = 250.0
latgrid = 10.
longrid = 10.
nmapfile = 1
mapfile(1) = '/home/mxue/zxplot3/data/us_state.mapdata'
mapfile(2) = '/home/mxue/zxplot3/data/world_us_country.mapdata'
terndmp = 1
!
!-----------------------------------------------------------------------
!
! Read and Print of Namelist variables and parameters.
!
!-----------------------------------------------------------------------
!
READ(5,grid_dims,ERR=980,END=990)
allocate( allocate(
allocate(xh(nx,ny))
allocate(yh(nx,ny))
allocate(xh1d(0:nx))
allocate(yh1d(0:ny))
allocate(lath(0:nx,0:ny))
allocate(lonh(0:nx,0:ny))
allocate(hterain(nx,ny))
allocate(tem (nx,ny))
allocate(tem2 (nx,ny))
allocate(xw (8*nx))
allocate(yw (8*nx))
!
! Default value of -9999.0 for the flag is used.
!
READ (5,jobname,ERR=980,END=990)
WRITE(6,'(/a,a)') 'The name of this run is: ', runname
CALL gtlfnkey
( runname, lfnkey )
READ(5,grid,ERR=980,END=990)
WRITE(6,'(/a,f10.3,a)')'Input dx was ',dx,' meters'
WRITE(6,'(/a,f10.3,a)')'Input dy was ',dy,' meters'
READ(5,projection,ERR=980,END=990)
WRITE(6,'(/a,i4)') 'Input mapproj was ',mapproj
WRITE(6,'(/a,f10.3,a)') &
'Input trulat1 was ',trulat1,' degree North'
WRITE(6,'(/a,f10.3,a)') &
'Input trulat2 was ',trulat2,' degree North'
WRITE(6,'(/a,f10.3)') &
'The latitude of the center of the model domain was ',ctrlat
WRITE(6,'(/a,f10.3)') &
'The longitude of the center of the model domain was ',ctrlon
WRITE(6,'(/a,f10.3,a)') &
'Input trulon was ',trulon,' degree East'
WRITE(6,'(/a,e15.5)') &
'Input sclfct was ',sclfct
IF ( mapproj == 0 ) THEN
trulat1 = ctrlat
trulat2 = ctrlat
trulon = ctrlon
END IF
trnanxopt = 1
trndataopt= 2
nsmth = 0
fn_trndata = 'tbase_global_5min.data'
usgs_dem_index_fn = 'usgs_dem.index'
READ(5,dem_trn,ERR=980,END=990)
IF( trndataopt == 1) THEN
trndatares = 3600
lon_sample = max(3600,lon_sample)
lat_sample = max(3600,lat_sample)
ELSE IF( trndataopt == 2) THEN
trndatares = 300
lon_sample = max(300,lon_sample)
lat_sample = max(300,lat_sample)
ELSE IF( trndataopt == 3) THEN
trndatares = 30
lon_sample = max(30,lon_sample)
lat_sample = max(30,lat_sample)
ELSE IF( trndataopt == 4) THEN
trndatares = 3
lon_sample = max(3,lon_sample)
lat_sample = max(3,lat_sample)
END IF
WRITE(6,'(/a,i5,a)') &
'The sample internal in longitude dir. was set to ',lon_sample, &
' seconds.'
WRITE(6,'(/a,i5,a)') &
'The sample internal in latitude dir. was set to ',lat_sample, &
' seconds.'
IF(MOD(lat_sample,trndatares) /= 0 .OR. MOD(lon_sample,trndatares) /= 0) THEN
lon_sample = MAX(1,lon_sample/trndatares) * trndatares
lat_sample = MAX(1,lat_sample/trndatares) * trndatares
WRITE(6,'(/a,i4,a,/a,i6,a,i6/))') &
'lat_sample and lon_sample have to be divisible by ', &
trndatares,'.', &
'They were reset to lon_sample=',lon_sample, &
' lat_sample=',lat_sample
END IF
IF(MOD(3600,lat_sample) /= 0 .OR. MOD(3600,lon_sample) /= 0) THEN
WRITE(6,'(/a,/a/))') &
'3600 has to be divisible by both lat_sample and lon_sample.', &
'Please reset them in the input file. Job stopped.'
STOP
END IF
latgrid = 0.0
longrid = 0.0
READ(5,map_plot)
READ(5,trn_output)
WRITE(6,'(/a,i4)') &
'Terrain dumping option set to ',terndmp
WRITE(6,'(/a/)')'Finished reading namelist input file.'
ldir_trndata = 120
CALL strlnth( dir_trndata, ldir_trndata)
IF( dir_trndata(ldir_trndata:ldir_trndata) /= '/') THEN
ldir_trndata=ldir_trndata+1
dir_trndata(ldir_trndata:ldir_trndata)='/'
END IF
!
!-----------------------------------------------------------------------
!
! setting ctrlon in terms of degrees east
!
!-----------------------------------------------------------------------
!
! IF(ctrlon.lt.0)THEN
! ctrlon=360+ctrlon
! END IF
!
!-----------------------------------------------------------------------
!
! test to see if reset value is within longitudinal limits..
!
!-----------------------------------------------------------------------
!
! IF(ctrlon.lt.0.or.ctrlon.gt.360)THEN
! print *,'ctrlon is incorrectly set in terrain.input, please
! : reset, use degrees east'
! STOP
! END IF
!
!-----------------------------------------------------------------------
!
! test to see if the ctrlat is within limits...
!
!-----------------------------------------------------------------------
!
IF(ctrlat < 0.OR.ctrlat > 90)THEN
PRINT *,'ctrlat in terrain.input is too large or too small. ', &
'must be set between 0. AND +90. degree north.'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! setting trulon in terms of degrees east
!
!-----------------------------------------------------------------------
!
! IF(trulon.lt.0)THEN
! trulon=360+trulon
! END IF
!
!-----------------------------------------------------------------------
!
! test to see if reset value is within longitudinal limits..
!
!-----------------------------------------------------------------------
!
! IF(trulon.lt.0.or.trulon.gt.360)THEN
! print *,'trulon is incorrectly set in terrain.input, please
! & reset, use degrees east'
! STOP
! END IF
!
!-----------------------------------------------------------------------
!
! test to see if trulat is within limits...
!
!-----------------------------------------------------------------------
!
IF(trulat1 < 0.OR.trulat1 > 90)THEN
PRINT *,'Input parameter trulat1 is too large or too small.', &
'must be set between 0. AND +90. degree north.'
STOP
END IF
IF(trulat2 < 0.OR.trulat2 > 90)THEN
PRINT *,'Input parameter trulat2 is too large or too small.', &
'must be set between 0. AND +90. degree north.'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Set up map projection
!
!-----------------------------------------------------------------------
!
trulat(1) = trulat1
trulat(2) = trulat2
CALL setmapr( mapproj, 1.0 , trulat , trulon)
! set up parameters for map projection
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
swx = ctrx - (FLOAT(nx-3)/2.) * dx
swy = ctry - (FLOAT(ny-3)/2.) * dy
CALL setorig( 1, swx, swy)
!
!-----------------------------------------------------------------------
!
! Set the horizontal coordinates of the model grid
!
!-----------------------------------------------------------------------
!
DO i=1,nx
x(i) = dx * (i-2)
END DO
DO j=1,ny
y(j) = dy * (j-2)
END DO
DO j=1,ny
DO i=1,nx
xh(i,j)=(i-1.5)*dx
yh(i,j)=(j-1.5)*dy
END DO
END DO
DO j=0,ny
DO i=0,nx
xh1d(i)=(i-1.5)*dx
yh1d(j)=(j-1.5)*dy
END DO
END DO
PRINT*,' dx,dy=',dx,dy
!
!-----------------------------------------------------------------------
!
! Final out the lat/lon of each scalar (terrain data) point.
!
!-----------------------------------------------------------------------
!
! Include an extra zone so that data grid covering
! (lonmax-lonmin)*(latmax-latmin) is large enough area for 2nd order
! interpolation
!
CALL xytoll(nx+1,ny+1,xh1d,yh1d,lath,lonh)
!
!-----------------------------------------------------------------------
!
! Final out the max/min lat/lon of each scalar (terrain data) point.
!
!-----------------------------------------------------------------------
CALL a3dmax0(lath,0,nx,0,nx,0,ny,0,ny,1,1,1,1,latmax,latmin)
CALL a3dmax0(lonh,0,nx,0,nx,0,ny,0,ny,1,1,1,1,lonmax,lonmin)
IF( trndataopt == 3 ) THEN
swlat_trn = INT((1000.0+latmin-15.0/3600.0)*0.1)*10-1000
swlon_trn = INT((1000.0+lonmin-15.0/3600.0)*0.1)*10-1000
nelat_trn = INT((1000.0+latmax+15.0/3600.0)*0.1)*10-1000
nelon_trn = INT((1000.0+lonmax+15.0/3600.0)*0.1)*10-1000
! IF( nelat_trn /= latmax ) nelat_trn = nelat_trn + 10
! IF( nelon_trn /= lonmax ) nelon_trn = nelon_trn + 10
IF( nelat_trn < latmax+15.0/3600.0 ) nelat_trn = nelat_trn + 10
IF( nelon_trn < lonmax+15.0/3600.0 ) nelon_trn = nelon_trn + 10
swlat_offset = 15.0/3600.0
swlon_offset = 15.0/3600.0
ELSE
swlat_trn = INT(1000.0+latmin)-1000
swlon_trn = INT(1000.0+lonmin)-1000
nelat_trn = INT(1000.0+latmax)-1000
nelon_trn = INT(1000.0+lonmax)-1000
IF( nelat_trn /= latmax ) nelat_trn = nelat_trn + 1
IF( nelon_trn /= lonmax ) nelon_trn = nelon_trn + 1
swlat_offset = 0.0
swlon_offset = 0.0
ENDIF
swlat_trn=MAX( -90, swlat_trn)
nelat_trn=MIN( 90, nelat_trn)
swlon_trn=MAX(-180, swlon_trn)
nelon_trn=MIN( 180, nelon_trn)
PRINT*,'lonmin, lonmax=',lonmin, lonmax
PRINT*,'latmin, latmax=',latmin, latmax
PRINT*,'sw lat/lon, ne lat/lon=', &
swlat_trn,swlon_trn,nelat_trn,nelon_trn
IF( trndataopt == 3 ) THEN
nx_trn = (3600*(nelon_trn-swlon_trn)-30.0)/lon_sample+1
ny_trn = (3600*(nelat_trn-swlat_trn)-30.0)/lat_sample+1
ELSE
nx_trn = (3600*(nelon_trn-swlon_trn))/lon_sample+1
ny_trn = (3600*(nelat_trn-swlat_trn))/lat_sample+1
ENDIF
Write(6,'(/2(a,i7),a/)') &
'Allocating array h_trn of size ',nx_trn,'x',ny_trn,'.'
allocate(h_trn(nx_trn,ny_trn),STAT=istatus)
IF( istatus /= 0 ) THEN
Write(6,'(/a/)')'Allocation of array h_trn failed. Program stopped.'
STOP
END IF
badvalue = -9999.0
DO i=1,nx_trn
DO j=1,ny_trn
h_trn(i,j)=badvalue
END DO
END DO
!-----------------------------------------------------------------------
!
! Fill in array h_trn, by reading in 3 arc second DEM/terrain data.
!
!-----------------------------------------------------------------------
IF( trndataopt == 1) THEN
! CALL FILLHTRN1DEG(h_trn,nx_trn,ny_trn, &
! swlon_trn,nelon_trn,lon_sample,swlat_trn,nelat_trn,lat_sample,&
! dir_trndata(1:ldir_trndata),badvalue)
ELSE IF( trndataopt == 2) THEN
CALL fillhtrn5min(h_trn,nx_trn,ny_trn, &
swlon_trn,nelon_trn,lon_sample,swlat_trn,nelat_trn,lat_sample, &
dir_trndata(1:ldir_trndata),fn_trndata,badvalue)
ELSE IF( trndataopt == 3) THEN
print*,'nx_trn,ny_trn,nx_trn',nx_trn,ny_trn
CALL FILLHTRN30SEC(h_trn,nx_trn,ny_trn, &
swlon_trn,nelon_trn,lon_sample,swlat_trn,nelat_trn,lat_sample,&
dir_trndata(1:ldir_trndata),badvalue)
ELSE IF( trndataopt == 4) THEN
lfn = len_trim(usgs_dem_index_fn)
CALL fillhtrn3sec(h_trn,nx_trn,ny_trn, &
swlon_trn,nelon_trn,lon_sample,swlat_trn,nelat_trn,lat_sample, &
dir_trndata(1:ldir_trndata),usgs_dem_index_fn(1:lfn),badvalue)
ELSE
WRITE(6,'(1x,a/1x,a,i3,a)') 'trndataopt must be from 1 to 4.', &
'The input value was ', trndataopt,', Program stopped.'
STOP 991
END IF
CALL a3dmax(h_trn,1,nx_trn,1,nx_trn,1,ny_trn,1,ny_trn, &
1,1,1,1, &
hmax,hmin,imax,jmax,kmax, imin,jmin,kmin)
WRITE(6,'(2(1x,a,g25.12,2(a,i5)))') &
'hmin =',hmin,' at i=',imin,', j=',jmin, &
'hmax =',hmax,' at i=',imax,', j=',jmax
WRITE(6,'(/a)') 'Finished reading terrain data.'
! print*,' nx_trn, ny_trn =', nx_trn, ny_trn
inith_set = 0
badvfound = 0
DO i=1,nx_trn
DO j=1,ny_trn
IF( h_trn(i,j) /= badvalue)THEN
IF(inith_set == 0 ) THEN
hmax = h_trn(i,j)
hmin = h_trn(i,j)
inith_set = 1
ELSE
hmax = MAX(hmax, h_trn(i,j))
hmin = MIN(hmin, h_trn(i,j))
END IF
ELSE
badvfound = 1
END IF
END DO
END DO
IF( badvfound == 1) THEN
WRITE(6,'(/a,/a)') &
'Missing value(s) found in sampled terrain array.', &
'The missing data need to be filled. Job stopped.'
STOP
END IF
IF(inith_set == 0 ) THEN ! If still not set
hmax = 0.0
hmin = 0.0
END IF
WRITE(6,'(/a,2f10.3)') &
'Min and max in the sampled terrain data array=',hmin, hmax
!
!-----------------------------------------------------------------------
!
! Initialize ARPS terrain array hterain from h_trn using bi-linear
! interpolation.
!
!-----------------------------------------------------------------------
!
PRINT*,'swlon_trn, swlat_trn=',swlon_trn, swlat_trn
PRINT*,'lon_sample,lat_sample=',lon_sample,lat_sample
lond_sample = lon_sample/3600.0
latd_sample = lat_sample/3600.0
IF( trnanxopt == 1) THEN ! 1st order interpolation
DO i=1,nx-1
DO j=1,ny-1
ih=MIN(nx_trn-1, &
MAX(1,INT((lonh(i,j)-(swlon_trn+swlon_offset))*3600.0/lon_sample)+1))
jh=MIN(ny_trn-1, &
MAX(1,INT((lath(i,j)-(swlat_trn+swlat_offset))*3600.0/lat_sample)+1))
! print*,'i,j,ih,jh=', i,j,ih,jh,h_trn(ih,jh)
IF(h_trn(ih ,jh ) == badvalue .OR. h_trn(ih+1,jh) == badvalue .OR. &
h_trn(ih+1,jh+1) == badvalue .OR. &
h_trn(ih ,jh+1) == badvalue) THEN
hterain(i,j)=badvalue
ELSE
lonh1=lonh(i,j)-(swlon_trn+swlon_offset)-(ih-1)*lond_sample
lath1=lath(i,j)-(swlat_trn+swlat_offset)-(jh-1)*latd_sample
w1=(latd_sample-lath1)*(lond_sample-lonh1)
w2=(latd_sample-lath1)*lonh1
w3=lath1*lonh1
w4=lath1*(lond_sample-lonh1)
hterain(i,j)=(w1*h_trn(ih ,jh )+w2*h_trn(ih+1,jh) &
+w3*h_trn(ih+1,jh+1)+w4*h_trn(ih,jh+1)) &
/(latd_sample*lond_sample)
END IF
END DO
END DO
ELSE ! 2nd order interpolation
DO i=1,nx-1
DO j=1,ny-1
ih=MIN(nx_trn-1, &
MAX(2,INT((lonh(i,j)-(swlon_trn+swlon_offset))*3600.0/lon_sample)+1))
jh=MIN(ny_trn-1, &
MAX(2,INT((lath(i,j)-(swlat_trn+swlat_offset))*3600.0/lat_sample)+1))
IF(h_trn(ih ,jh ) == badvalue .OR. h_trn(ih+1,jh) == badvalue .OR. &
h_trn(ih+1,jh+1) == badvalue .OR. &
h_trn(ih ,jh+1) == badvalue) THEN
hterain(i,j)=badvalue
ELSE
IF(h_trn(ih-1,jh+1) == badvalue .OR. &
h_trn(ih-1,jh) == badvalue .OR. &
h_trn(ih-1,jh-1) == badvalue .OR. &
h_trn(ih+1,jh-1) == badvalue .OR. &
h_trn(ih ,jh-1) == badvalue) THEN ! Use 1st order
hterain(i,j)=badvalue
lonh1=lonh(i,j)-(swlon_trn+swlon_offset)-(ih-1)*lond_sample
lath1=lath(i,j)-(swlat_trn+swlat_offset)-(jh-1)*latd_sample
w1=(latd_sample-lath1)*(lond_sample-lonh1)
w2=(latd_sample-lath1)*lonh1
w3=lath1*lonh1
w4=lath1*(lond_sample-lonh1)
hterain(i,j)=(w1*h_trn(ih ,jh )+w2*h_trn(ih+1,jh) &
+w3*h_trn(ih+1,jh+1)+w4*h_trn(ih,jh+1)) &
/(latd_sample*lond_sample)
ELSE ! now real 2nd order interpolation
d=lonh(i,j)-(swlon_trn+swlon_offset)-(ih-1)*lond_sample
a=-lond_sample
b=0.0
c=+lond_sample
w1= (d-b)*(d-c)/((a-b)*(a-c))
w2= (d-a)*(d-c)/((b-a)*(b-c))
w3= (d-a)*(d-b)/((c-a)*(c-b))
h1=w1*h_trn(ih-1,jh-1)+w2*h_trn(ih,jh-1) &
+w3*h_trn(ih+1,jh-1)
h2=w1*h_trn(ih-1,jh )+w2*h_trn(ih,jh ) &
+w3*h_trn(ih+1,jh )
h3=w1*h_trn(ih-1,jh+1)+w2*h_trn(ih,jh+1) &
+w3*h_trn(ih+1,jh+1)
d=lath(i,j)-(swlat_trn+swlat_offset)-(jh-1)*latd_sample
a=-latd_sample
b=0.0
c=+latd_sample
w1= (d-b)*(d-c)/((a-b)*(a-c))
w2= (d-a)*(d-c)/((b-a)*(b-c))
w3= (d-a)*(d-b)/((c-a)*(c-b))
hterain(i,j)=w1*h1+w2*h2+w3*h3
END IF
END IF
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in Great Lakes. The original data set corresponds to the
! bottom of the lakes.
!
!-----------------------------------------------------------------------
!
DO j=1,ny
DO i=1,nx
tem(i,j) = lath(i,j)
tem2(i,j) = lonh(i,j)
END DO
END DO
CALL fix_lake_eliv(hterain,tem,tem2,nx,ny)
CALL a3dmax(hterain,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
hmax,hmin,imax,jmax,kmax, imin,jmin,kmin)
WRITE(6,'(2(1x,a,g25.12,2(a,i5)))') &
'htrnmin =',hmin,' at i=',imin,', j=',jmin, &
'htrnmax =',hmax,' at i=',imax,', j=',jmax
!-----------------------------------------------------------------------
!
! Apply smoothing on terrain height
!
!-----------------------------------------------------------------------
DO i=1,nsmth
CALL smth(hterain,nx,1,nx-1,ny,1,ny-1,2,tem)
! call SMTH(hterain,nx,1,nx-1,ny,1,ny-1,1,tem)
END DO
!
!-----------------------------------------------------------------------
!
! Write out the terrain file on ARPS grid.
!
!-----------------------------------------------------------------------
!
hmax=hterain(1,1)
hmin=hterain(1,1)
DO i=1,nx-1
DO j=1,ny-1
hterain(i,j) = MAX(0.0,hterain(i,j))
hmax=MAX(hmax,hterain(i,j))
hmin=MIN(hmin,hterain(i,j))
END DO
END DO
WRITE(6,'(/a,2f10.3)') &
'Max/min in the analyzed terrain array was ',hmax,hmin
IF (terndmp /= 0) THEN
CALL writtrn(nx,ny,runname(1:lfnkey)//'.trndata',dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,hterain)
END IF
DO j=1,ny
DO i=1,nx
xh(i,j)=xh(i,j)*0.001
yh(i,j)=yh(i,j)*0.001
! hterain(i,j)=hterain(i,j)*0.001
END DO
END DO
!wdt skip graphics
STOP
!wdt update
!
!-----------------------------------------------------------------------
!
! Generate graphical output.
!
!-----------------------------------------------------------------------
!
CALL xdevic
CALL xctrbadv( 1 ) ! Turn on missing value checking for contouring
CALL xvtrbadv( 1 ) ! Turn on missing value checking for vector plotting
CALL xsetclrs(6)
CALL xaxfmt('(I6)')
CALL xbadval ( badvalue ) ! Set the missing value flag to -9999.0
xl = (nx-3)*dx * 0.001
yl = (ny-3)*dy * 0.001
xorig = 0.0
yorig = 0.0
CALL xstpjgrd(mapproj,trulat1,trulat2,trulon, &
ctrlat,ctrlon,xl,yl,xorig,yorig)
CALL plot_trn(hterain,xh,yh,nx,1,nx-1,ny,1,ny-1, &
-dx*0.001,(nx-2)*dx*0.001,dx*0.001, &
-dy*0.001,(ny-2)*dy*0.001,dy*0.001, hctrinc, tem,xw,yw)
DO i=1,nmapfile
lmapfile=LEN(mapfile(i))
CALL strlnth(mapfile(i), lmapfile)
WRITE(6,'(1x,a,a)') 'Input was ',mapfile(i)(1:lmapfile)
INQUIRE(FILE=mapfile(i)(1:lmapfile), EXIST = fexist )
IF( .NOT.fexist) THEN
WRITE(6,'(a)') 'Map file '//mapfile(i)(1:lmapfile)// &
' not found. Program stopped in ARPSPLT.'
STOP 001
END IF
CALL xdrawmap(10,mapfile(i)(1:lmapfile),latgrid,longrid)
END DO
! ctrlat = 40.28
! ctrlon = -111.56
! CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
! call xcolor(1)
! call xthick(3)
! call xpenup(ctrx*0.001,ctry*0.001-0.03*(ny-3)*dy*0.001)
! call xpendn(ctrx*0.001,ctry*0.001+0.03*(ny-3)*dy*0.001)
! call xpenup(ctrx*0.001-0.03*(nx-3)*dx*0.001,ctry*0.001)
! call xpendn(ctrx*0.001+0.03*(nx-3)*dx*0.001,ctry*0.001)
! ctrlat = 40.5
! ctrlon = -111.5
! CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
! call xcolor(1)
! call xthick(3)
! call xpenup(ctrx*0.001,ctry*0.001-0.03*(ny-3)*dy*0.001)
! call xpendn(ctrx*0.001,ctry*0.001+0.03*(ny-3)*dy*0.001)
! call xpenup(ctrx*0.001-0.03*(nx-3)*dx*0.001,ctry*0.001)
! call xpendn(ctrx*0.001+0.03*(nx-3)*dx*0.001,ctry*0.001)
!
CALL xgrend
STOP
980 CONTINUE
WRITE(6,'(/1x,a,a)') &
'Error occured when reading namelist input file. Program stopped.'
STOP
990 CONTINUE
WRITE(6,'(/1x,a,a)') &
'End of file reached when reading namelist input file. ', &
'Program stopped.'
STOP
END PROGRAM arpstrn
SUBROUTINE fillhtrn3sec(h_trn,nx_trn,ny_trn, & 1,3
swlon_trn,nelon_trn,lon_sample, &
swlat_trn,nelat_trn,lat_sample,dir_trndata, &
usgs_dem_index_fn,badvalue)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/25/95
!
! MODIFICATION HISTORY:
!
! 2000/02/03 Gene Bassett
! Added subroutine to fix the elevations in the Great Lakes.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: nlat,nlon ! Number of data points in the 3"
! one degree DEM data block.
INTEGER, ALLOCATABLE :: h_block(:,:) ! Array to store 1x1 degree DEM patch
INTEGER :: nx_trn,ny_trn
REAL :: h_trn(nx_trn,ny_trn)
CHARACTER (LEN=80) :: trn_name
INTEGER :: lat_sample,lon_sample ! Data sampling frequency in data index
INTEGER :: swlat_trn,swlon_trn,nelat_trn,nelon_trn
! SW and NE coner lat/lon of the
! the assembled terrain patch.
INTEGER :: istatus,istatus1,istatus2,ii,jj,i,j,ilat,ilon
REAL :: badvalue, tmax0,tmin0
CHARACTER (LEN=3) :: dtflag
LOGICAL :: existfile
CHARACTER (LEN=*) :: dir_trndata,usgs_dem_index_fn
INTEGER :: lat_skip, lon_skip,ldir_trndata
!
!-----------------------------------------------------------------------
!
! Make sure directory dir_trndata already exists.
!
!-----------------------------------------------------------------------
ldir_trndata = len_trim(dir_trndata)
INQUIRE(FILE=dir_trndata(1:ldir_trndata),EXIST=existfile)
IF( .NOT. existfile ) THEN
WRITE(6,'(/1x,3a,3(/1x,a))') &
'Directory ',dir_trndata(1:ldir_trndata),' does not exist.', &
'Create the directory or create a link to the data directory.', &
'This directory has to be writable, unless all required DEM data', &
'exists locally.'
END IF
nlat = 1201
nlon = 1201
ALLOCATE (h_block(nlat,nlon))
DO ilat=swlat_trn,nelat_trn-1
DO ilon=swlon_trn,nelon_trn-1
!-----------------------------------------------------------------------
!
! Find the name of DEM file whose SW corner lat/lon = ilat/ilon.
!
!-----------------------------------------------------------------------
CALL lookup_trn(ilat,ilon,usgs_dem_index_fn, &
trn_name,dtflag,istatus)
IF( istatus /= 0) GO TO 310
!-----------------------------------------------------------------------
!
! Acquire file trn_name from anonymous FTP server
! edcftp.cr.usgs.gov if it does not already exist locally in dir_trndata
!
!-----------------------------------------------------------------------
CALL get_trn(dir_trndata(1:ldir_trndata), &
trn_name,dtflag,istatus1)
IF(istatus1 /= 0) GO TO 310
!
!-----------------------------------------------------------------------
!
! Read in 1x1 degree DEM data block into array h_block.
!
!-----------------------------------------------------------------------
!
CALL read_trn(dir_trndata(1:ldir_trndata)//trn_name, &
h_block,nlat,nlon,istatus2)
IF(istatus2 /= 0) GO TO 310
tmax0=h_block(1,1)
tmin0=h_block(1,1)
DO i=1,1201
DO j=1,1201
tmax0=MAX(tmax0,FLOAT(h_block(i,j)))
tmin0=MIN(tmin0,FLOAT(h_block(i,j)))
END DO
END DO
PRINT*,'max/min in data file ',dir_trndata(1:ldir_trndata), &
'=',tmax0,tmin0
!
!-----------------------------------------------------------------------
!
! Transfer selected data points from h_block into h_trn.
!
!-----------------------------------------------------------------------
!
lat_skip = max(1,lat_sample/3)
lon_skip = max(1,lon_sample/3)
ii = (nlon-1)/lon_skip*(ilon-swlon_trn)
jj = (nlat-1)/lat_skip*(ilat-swlat_trn)
DO j=1,(nlat-1)/lat_skip+1
DO i=1,(nlon-1)/lon_skip+1
h_trn(ii+i,jj+j)= &
FLOAT(h_block((i-1)*lon_skip+1,(j-1)*lat_skip+1))
END DO
END DO
CYCLE
310 CONTINUE
WRITE(6,'(/1x,a,2i5,a,/1x,a)') &
'DEM file with S-W lat/lon=',ilat,ilon,' was not found or not', &
'read correctly. The DEM array is filled with missing values.'
ii = (nlon-1)/lon_skip*(ilon-swlon_trn)
jj = (nlat-1)/lat_skip*(ilat-swlat_trn)
DO j=1,(nlat-1)/lat_skip+1
DO i=1,(nlon-1)/lon_skip+1
h_trn(ii+i,jj+j)= badvalue
END DO
END DO
END DO
END DO
DEALLOCATE(h_block)
RETURN
END SUBROUTINE fillhtrn3sec
SUBROUTINE fillhtrn5min(h_trn,nx_trn,ny_trn, & 1,2
swlon_trn,nelon_trn,lon_sample, &
swlat_trn,nelat_trn,lat_sample,dir_trndata,fn_trndata, &
badvalue)
!-----------------------------------------------------------------------
!
! read the TerrainBase data
! This original program modified by LDI
! This promgram used to make the 5min topo input data using NCAR 5min topo
! output is the 5 min topo for ARPS terrain program
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: itopo(4320) ! Array for topo along one latitude circle
INTEGER :: imin,imax,jmin,jmax,i,j
INTEGER :: nx_trn,ny_trn
REAL :: h_trn(nx_trn,ny_trn)
INTEGER :: swlat_trn,swlon_trn,nelat_trn,nelon_trn
! SW and NE coner lat/lon of the assembled terrain patch.
INTEGER :: lon_sample, lat_sample
CHARACTER (LEN=80) :: fn_trndata
CHARACTER (LEN=120) :: dir_trndata
INTEGER :: ldir_trndata,lfn_trndata
REAL :: badvalue
REAL :: hmax, hmin
INTEGER :: inith_set , iskip, jskip, ii
LOGICAL existfile
imin = swlon_trn*12+1
imax = nelon_trn*12+1
jmax = (90-swlat_trn)*12+1
jmin = (90-nelat_trn)*12+1
PRINT*,'imin,imax,jmin,jmax=',imin,imax,jmin,jmax
ldir_trndata = 120
CALL strlnth( dir_trndata, ldir_trndata)
lfn_trndata = 80
CALL strlnth( fn_trndata, lfn_trndata)
INQUIRE(FILE=dir_trndata(1:ldir_trndata)//fn_trndata(1:lfn_trndata), &
EXIST=existfile)
IF ( .not. existfile) THEN
WRITE(6,'(a,a,/a)') &
'File ',dir_trndata(1:ldir_trndata)//fn_trndata(1:lfn_trndata), &
' not found. Program stopped in subroutine FILLHTRN5MIN.'
STOP
END IF
OPEN(10,FILE=dir_trndata(1:ldir_trndata)// &
fn_trndata(1:lfn_trndata),STATUS='old',FORM='formatted')
iskip = lon_sample/300
jskip = lat_sample/300
DO j=1,jmax ! jmax =< 2160
READ(10,100)(itopo(i),i=1,4320)
IF(j >= jmin.AND. &
MOD((jmax-jmin+1)-(j-jmin+1),jskip) == 0)THEN ! Fill in array h_trn
DO i=imin, imax, iskip
ii = i
IF( i < 0) ii=4320+i
h_trn((i-imin)/iskip+1, &
((jmax-jmin+1)-(j-jmin+1))/jskip+1)=itopo(ii)
END DO
END IF
END DO
100 FORMAT(20I6)
! CALL XDEVIC
IF( .false.) THEN ! Plotting is true
nx_trn = (12*(nelon_trn-swlon_trn))/iskip+1
ny_trn = (12*(nelat_trn-swlat_trn))/jskip+1
! print*,'nx_trn,ny_trn,imax-imin+1,jmax-jmin+1=',
! : nx_trn,ny_trn,imax-imin+1,jmax-jmin+1
inith_set = 0
DO i=1,nx_trn
DO j=1,ny_trn
IF( h_trn(i,j) /= badvalue)THEN
IF(inith_set == 0 ) THEN
hmax = h_trn(i,j)
hmin = h_trn(i,j)
inith_set = 1
ELSE
hmax = MAX(hmax, h_trn(i,j))
hmin = MIN(hmin, h_trn(i,j))
END IF
END IF
END DO
END DO
IF(inith_set == 0 ) THEN ! If still not set
hmax = 0.0
hmin = 0.0
END IF
PRINT*,' hmin, hmax =', hmin, hmax
! hmax = h_trn(1,1)
! hmin = h_trn(1,1)
! DO 1203 i=1,nx_trn
! DO 1203 j=1,ny_trn
! hmax = max(hmax, h_trn(i,j))
! hmin = min(hmin, h_trn(i,j))
! 1203 CONTINUE
! print*,' hmin, hmax =', hmin, hmax
END IF
CLOSE(2)
RETURN
END SUBROUTINE fillhtrn5min
SUBROUTINE lookup_trn(swlat,swlon,usgs_dem_index_fn, & 1
trn_name,dtflag,ierr)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Find a 1x1 degree DEM file whose SW corner is at (swlat,swlon).
! dtflag= 'US' or 'AK', indicating if the data file is in Alaska.
! Note: The longitude W is given as positive numbers in the Index.
! This should be fine as long as this program is used for
! accessing US 3" data only. The code should be modified slightly
! to handle global data.
!
! INPUT : swlat, swlon,usgs_dem_index_fn
! OUTPUT: trn_name, dtflag, ierr
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: swlat,swlon
CHARACTER (LEN=80) :: trn_name,name1
CHARACTER (LEN=3) :: dtflag
INTEGER :: ierr
INTEGER :: selat,selon,lat,lon,ll
REAL :: alat,alon
CHARACTER (LEN=*) :: usgs_dem_index_fn
!
OPEN (50,FILE=trim(usgs_dem_index_fn),STATUS='old',IOSTAT=ierr)
IF( ierr /= 0) THEN
WRITE(6,'(/1x,a,a,a,/1x,a)') 'Opening file ', &
trim(usgs_dem_index_fn),' failed. Program stopped.'
STOP
END IF
READ (50,*)
!
selat = swlat
selon = -(swlon + 1)
WRITE(6,'(/1x,2(a,i4))') &
'Searching for a file with the SE corner latitude=',selat, &
' longitude=', selon
100 READ (50,'(a3,a30,2f10.1)',END=300) dtflag,name1, alat, alon
lat=nint(alat)
lon=nint(alon)
IF (lat == selat.AND.lon == selon) THEN
trn_name=name1
ll=INDEX(trn_name,'.gz')
WRITE(6,'(1x,a,a,a)') 'Found the entry ',trn_name(1:ll+2), &
' with the matching lat/lon.'
CLOSE (50)
RETURN
ELSE
GO TO 100
END IF
300 CONTINUE
CLOSE (50)
PRINT*,'No entry in data base for this latitude and longitude.'
ierr = 1
RETURN
END SUBROUTINE lookup_trn
SUBROUTINE get_trn(dir_trndata,trn_name,dtflag,istatus) 1,3
IMPLICIT NONE
CHARACTER (LEN=120) :: dir_trndata
CHARACTER (LEN=80) :: trn_name
CHARACTER (LEN=3) :: dtflag
INTEGER :: istatus,ldir_trndata
CHARACTER (LEN=120) :: ch1,ch2,ch3
CHARACTER (LEN=80) :: dr
LOGICAL :: existfile,existfile1
INTEGER :: l1
l1=INDEX(trn_name,'.gz')
IF( dtflag == 'AK') THEN
dr='Alaska/'//CHAR(ICHAR(trn_name(1:1))-32)
ELSE
dr=trn_name(1:1)
dr=CHAR(ICHAR(dr(1:1))-32)
END IF
ch1='cd '//dr
ch2='dir '//trn_name(1:l1+2)
ch3='get '//trn_name(1:l1+2)
!
! check to see if dem file pre-exists
!
istatus = 1
ldir_trndata=120
CALL strlnth(dir_trndata,ldir_trndata)
INQUIRE(FILE=dir_trndata(1:ldir_trndata)//trn_name(1:l1+2), &
EXIST=existfile)
INQUIRE(FILE=dir_trndata(1:ldir_trndata)//trn_name(1:l1-1), &
EXIST=existfile1)
IF (existfile) THEN
PRINT *,trn_name(1:l1+2),' exists locally'
ELSE IF (existfile1 ) THEN
PRINT *,trn_name(1:l1-1),' exists locally'
ELSE
!
! write the ftp input file
!
OPEN (29,FILE='ftp.in')
WRITE (29,1) 'binary'
WRITE (29,1) 'cd pub/data/DEM/250'
WRITE (29,1) ch1
WRITE (29,1) ch2
WRITE (29,1) ch3
WRITE (29,1) 'quit'
CLOSE (29)
1 FORMAT(a)
!
! invoke ftp command
!
ch1='ftp edcftp.cr.usgs.gov< ftp.in'
PRINT *,'Connecting to ftp site to get the file...'
CALL unixcmd(ch1)
CALL unixcmd( &
'mv '//trn_name(1:l1+2) //' '//dir_trndata(1:ldir_trndata))
PRINT *,'File ',trn_name(1:l1+2), &
' transferred to the local machine now.'
END IF
istatus=0
!
RETURN
END SUBROUTINE get_trn
SUBROUTINE read_trn(trn_fullname,h_block,nlat,nlon,istatus) 1,7
IMPLICIT NONE
CHARACTER (LEN=128) :: trn_fullname
INTEGER :: nlat,nlon
! integer*2 h_block(nlat,nlon)
INTEGER :: h_block(nlat,nlon)
INTEGER :: istatus
CHARACTER (LEN=80) :: ch1,ch2
CHARACTER (LEN=120) :: ch4
DIMENSION dum15(15),dum6(6)
CHARACTER (LEN=144) :: head
INTEGER :: i,j,ii,ijump,m,n,l1,ltd,lnd,idum
REAL :: xjumpinv,tmax0,tmin0,tmax,tmin,dum,dum6,dum15,amx,amn
REAL :: selnsec,seltsec
LOGICAL :: existfile
istatus=1
l1=INDEX(trn_fullname,'.gz')
INQUIRE(FILE='tmp1.file', EXIST = existfile)
IF(existfile)THEN
ch2='/bin/rm tmp1.file'
CALL unixcmd(ch2)
ENDIF
INQUIRE(FILE=trn_fullname(1:l1-1),EXIST=existfile)
IF(existfile) THEN ! unzipped version exists
ch1='cp '//trn_fullname(1:l1-1)//' tmp1.file'
CALL unixcmd(ch1)
ELSE ! need to unzip
ch1='gunzip -c '//trn_fullname(1:l1+2)//'> tmp1.file'
PRINT *,'Unzipping the terrain file..'
CALL unixcmd(ch1)
PRINT *,'Terrain file unzipping competed..'
END IF
INQUIRE(FILE='tmp2.file', EXIST = existfile)
IF(existfile)THEN
ch2='/bin/rm tmp2.file'
CALL unixcmd(ch2)
ENDIF
ch4='dd if=tmp1.file of=tmp2.file ibs=4096 cbs=1024 conv=unblock'
PRINT *,'Converting the terrain file..'
CALL unixcmd(ch4)
PRINT *,'Conversion completed..'
INQUIRE(FILE='tmp1.file', EXIST = existfile)
IF(existfile)THEN
ch2='/bin/rm tmp1.file'
CALL unixcmd(ch2)
ENDIF
!
! Read in terrain data
!
OPEN (29,FILE='tmp2.file')
READ (29,11) head,idum,idum,idum,idum,(dum15(i),i=1,15), &
idum,idum,idum,(dum6(i),i=1,6),selnsec,seltsec,tmin,tmax,dum, &
idum,dum,dum,dum,idum,m
11 FORMAT(a144,4I6,15E24.0,3I6,11E24.0,i6,3E12.0,2I6)
PRINT*,'Min and max in the input DEM file header=', &
tmin,tmax
ltd=seltsec/3600.
lnd=selnsec/3600.
WRITE (6,'(/1x,a,2i5)') &
'South-east corner latitude and longitude from the file ',ltd,lnd
WRITE (6,'(/1x,a)') &
'Reading terrain data and store into the terrain array.'
IF( m == 1201) THEN
ijump=1
ELSE IF( m == 601) THEN
ijump=2
ELSE IF( m == 401) THEN
ijump=3
ELSE
WRITE(6,'(/1x,a,/1x,i6,a)') &
'The DEM file must have 1201, 601 or 401 data in the lon. dir.' &
,m,' data points were found. Program stopped.'
STOP
END IF
DO i=1,1201,ijump
READ (29,2) idum,idum,n,idum,dum,dum,dum,amx,amn, &
(h_block(i,j),j=1,n)
END DO
2 FORMAT (4I6,5E24.0,146I6,/(170I6))
IF( n /= 1201) THEN
WRITE(6,'(/1x,a,/1x,i6,a)') &
'The DEM file must have 1201 data in the lat. dir.' &
,n,' data points were found. Program stopped.'
STOP
END IF
IF( .true. ) THEN
tmax0=h_block(1,1)
tmin0=h_block(1,1)
DO i=1,1201,ijump
DO j=1,n
tmax0=MAX(tmax0,FLOAT(h_block(i,j)))
tmin0=MIN(tmin0,FLOAT(h_block(i,j)))
END DO
END DO
PRINT*,'READTRN: max/min in the data read in=',tmax0,tmin0
END IF
!
! Linearly interpolate terrain data to all 3 second grid points
! when the input data resolution is coarser than 3 second.
!
xjumpinv=1.0/ijump
DO ii=1,ijump-1
DO i=1,1201-ijump,ijump
DO j=1,n
h_block(i+ii,j)=h_block(i,j)+(h_block(i+ijump,j) &
-h_block(i,j))*ii*xjumpinv
END DO
END DO
END DO
IF( .false.) THEN
tmax0=h_block(1,1)
tmin0=h_block(1,1)
DO i=1,1201
DO j=1,n
tmax0=MAX(tmax0,FLOAT(h_block(i,j)))
tmin0=MIN(tmin0,FLOAT(h_block(i,j)))
END DO
END DO
PRINT*,'max/min in the data read in=',tmax0,tmin0
END IF
PRINT*,'Terrain data stored into the terrain array.'
!
CLOSE (29)
CALL unixcmd('/bin/rm tmp2.file')
istatus = 0
!
RETURN
END SUBROUTINE read_trn
SUBROUTINE fillhtrn30sec(h_trn,nx_trn,ny_trn, & 1,2
swlon_trn,nelon_trn,lon_sample, &
swlat_trn,nelat_trn,lat_sample,dir_trndata,badvalue)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 1/31/2002 (based on fillhtrn3sec)
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES:
!
!-----------------------------------------------------------------------
!
INTEGER :: nlat,nlon ! Number of data points in the 3"
! one degree DEM data block.
INTEGER, ALLOCATABLE :: h_block(:,:) ! Array to store 1x1 degree DEM patch
INTEGER :: nx_trn,ny_trn
REAL :: h_trn(nx_trn,ny_trn)
CHARACTER (LEN=80) :: trn_name
INTEGER :: lat_sample,lon_sample ! Data sampling frequency in data index
INTEGER :: swlat_trn,swlon_trn,nelat_trn,nelon_trn
! SW and NE coner lat/lon of the
! the assembled terrain patch.
INTEGER :: istatus,istatus1,istatus2,ii,jj,i,j,ilat,ilon
REAL :: badvalue, tmax0,tmin0
CHARACTER (LEN=3) :: dtflag
LOGICAL :: existfile
CHARACTER (LEN=*) :: dir_trndata
INTEGER :: lat_skip, lon_skip,ldir_trndata
!
!-----------------------------------------------------------------------
!
! Make sure directory dir_trndata already exists.
!
!-----------------------------------------------------------------------
print*,'Inside fillhtrn30sec '
ldir_trndata = len_trim(dir_trndata)
! INQUIRE(FILE=dir_trndata(1:ldir_trndata),EXIST=existfile)
! IF( .NOT. existfile ) THEN
! WRITE(6,'(/1x,3a,3(/1x,a))') &
! 'Directory ',dir_trndata(1:ldir_trndata),' does not exist.', &
! 'Create the directory or create a link to the data directory.', &
! 'This directory has to be writable, unless all required DEM data', &
! 'exists locally.'
! END IF
nlat = 1200
nlon = 1200
ALLOCATE(h_block(nlat,nlon))
print*,'done allocating h_block'
print*,'swlat_trn,nelat_trn,swlon_trn,nelon_trn=',swlat_trn,nelat_trn,swlon_trn,nelon_trn
DO ilat=swlat_trn,nelat_trn-10,10
DO ilon=swlon_trn,nelon_trn-10,10
print*,' '
print*,'SW corner lat, lon of 10x10 degree patch to be read=',ilat, ilon
!-----------------------------------------------------------------------
! Find the name of the 10x10 degree file whose SW corner lat/lon = ilat/ilon.
!-----------------------------------------------------------------------
IF(ilon.lt.0.0.and.ilat.ge.0.0)write(trn_name,'(a,i3.3,a,i3.3,a)') &
'W',ABS(ilon),'N',ABS(ilat),'.30s_ascdat'
IF(ilon.lt.0.0.and.ilat.lt.0.0)write(trn_name,'(a,i3.3,a,i3.3,a)') &
'W',ABS(ilon),'S',ABS(ilat),'.30s_ascdat'
IF(ilon.ge.0.0.and.ilat.ge.0.0)write(trn_name,'(a,i3.3,a,i3.3,a)') &
'E',ABS(ilon),'N',ABS(ilat),'.30s_ascdat'
IF(ilon.ge.0.0.and.ilat.lt.0.0)write(trn_name,'(a,i3.3,a,i3.3,a)') &
'E',ABS(ilon),'S',ABS(ilat),'.30s_ascdat'
!-----------------------------------------------------------------------
!
! Acquire file trn_name from CAPS anonymous FTP server
! ftp://ftp.caps.ou.edu/pub/ARPS/ARPS.data/arpstopo30.data if it does
! not already exist locally in dir_trndata
!
!-----------------------------------------------------------------------
CALL get_30s_trn(trim(dir_trndata),trim(trn_name),istatus1)
IF(istatus1 /= 0) GO TO 310
!
!-----------------------------------------------------------------------
!
! Read in 10x10 degree DEM data block into array h_block.
!
!-----------------------------------------------------------------------
!
print*,'To read data from ', dir_trndata(1:ldir_trndata)//trim(trn_name)
CALL read_30s_trn(dir_trndata(1:ldir_trndata)//trim(trn_name), &
h_block,nlat,nlon,istatus2)
print*,'Done reading data in ', dir_trndata(1:ldir_trndata)//trim(trn_name)
DO i=1,nlon
DO j=1,nlat
IF(h_block(i,j).eq.-9999)h_block(i,j)=0.0 ! in ocean by definition
! of the data set
ENDDO
ENDDO
IF(istatus2 /= 0) GO TO 310
!
!-----------------------------------------------------------------------
!
! Transfer selected data points from h_block into h_trn.
!
!-----------------------------------------------------------------------
!
lat_skip = max(1,lat_sample/30)
lon_skip = max(1,lon_sample/30)
ii = nlon/lon_skip*((ilon-swlon_trn)/10)
jj = nlat/lat_skip*((ilat-swlat_trn)/10)
DO j=1,nlat/lat_skip
DO i=1,nlon/lon_skip
h_trn(ii+i,jj+j)= &
FLOAT(h_block((i-1)*lon_skip+1,(j-1)*lat_skip+1))
END DO
END DO
print*,'data values transferred from h_block to h_trn'
CYCLE
310 CONTINUE
WRITE(6,'(/1x,a,2i5,a,/1x,a)') &
'30second data file with S-W lat/lon=',ilat,ilon,' was not found or not', &
'read correctly. The data array is filled with missing values.'
ii = nlon/lon_skip*((ilon-swlon_trn)/10)
jj = nlat/lat_skip*((ilat-swlat_trn)/10)
DO j=1,nlat/lat_skip
DO i=1,nlon/lon_skip
h_trn(ii+i,jj+j)= badvalue
END DO
END DO
END DO
END DO
DEALLOCATE( h_block )
print*,'end of fillhtrn30sec '
RETURN
END SUBROUTINE fillhtrn30sec
SUBROUTINE get_30s_trn(dir_trndata,trn_name,istatus) 1,3
IMPLICIT NONE
CHARACTER (LEN=*) :: dir_trndata
CHARACTER (LEN=*) :: trn_name
INTEGER :: istatus,ldir_trndata
CHARACTER (LEN=120) :: ch1
LOGICAL :: existfile,existfile1
INTEGER :: nunit
ch1='get '//trim(trn_name)//'.gz'
!
! check to see if dem file pre-exists
!
istatus = 1
INQUIRE(FILE=trim(dir_trndata)//trim(trn_name),EXIST=existfile)
INQUIRE(FILE=trim(dir_trndata)//trim(trn_name)//'.gz',EXIST=existfile1)
IF (existfile) THEN
PRINT *,trim(trn_name),' exists locally'
ELSE IF (existfile1 ) THEN
PRINT *,trim(trn_name)//'.gz',' exists locally'
ELSE
!
PRINT *,trim(trn_name),' or its compressed version not found in '
PRINT *,'in directory ',trim(dir_trndata), &
'. Will try to ftp it from ftp.caps.ou.edu.'
Print *,'It requires that you have .netrc setup in your home directory'
Print *,'for automatic anonymous ftp access to ftp.caps.ou.edu. '
! write the ftp input file
!
call getunit( nunit)
OPEN (nunit,FILE='ftp.in')
WRITE (nunit,1) 'binary'
WRITE (nunit,1) 'cd pub/ARPS/ARPS.data/arpstopo30.data'
WRITE (nunit,1) ch1
WRITE (nunit,1) 'quit'
CLOSE (nunit)
call retunit(nunit)
1 FORMAT(a)
!
! invoke ftp command
!
ch1='ftp ftp.caps.ou.edu < ftp.in'
PRINT *,'Connecting to ftp site to get the file...'
CALL unixcmd(ch1)
CALL unixcmd('mv '//trim(trn_name)//'.gz '//trim(dir_trndata))
PRINT *,'File ',trim(trn_name)//'.gz ', &
' transferred to the local machine now.'
END IF
istatus=0
!
RETURN
END SUBROUTINE get_30s_trn
SUBROUTINE read_30s_trn(trn_fullname,h_block,nlat,nlon,istatus) 1,4
IMPLICIT NONE
CHARACTER (LEN=*) :: trn_fullname
INTEGER :: nlat,nlon
INTEGER :: h_block(nlat,nlon)
INTEGER :: istatus
INTEGER :: i,j
LOGICAL :: existfile
character(len=132) :: ch1
character(len=100) :: filename
Doubleprecision :: ll_lon,ll_lat,lon_inc,lat_inc
INTEGER :: imissing
IF( nlat .ne. 1200 .or. nlon.ne. 1200 ) then
print*,'Array size for array h_block, nlat and nlon must be '
print*,'equal to 1200 inside read_30s_trn. '
print*,'Input nlat and nlon =', nlat, nlon
istatus = 1
RETURN
END IF
istatus=1
print*,'inside read_30s_trn, trn_fullname=', trim(trn_fullname)
INQUIRE(FILE=trim(trn_fullname),EXIST=existfile)
IF(existfile) THEN ! unzipped version exists
ch1='cp '//trim(trn_fullname)//' tmp1.file'
CALL unixcmd(ch1)
ELSE ! need to unzip
INQUIRE(FILE=trim(trn_fullname)//'.gz',EXIST=existfile)
IF(existfile) then
ch1='gunzip -c '//trim(trn_fullname)//'>tmp1.file'
PRINT *,'Unzipping the terrain file..'
CALL unixcmd(ch1)
PRINT *,'Terrain file unzipping competed..'
ELSE
print*,'File ',trim(trn_fullname),' or its gzipped version not found,'
print*,'Terrain data read failed.'
istatus = 1
STOP
RETURN
ENDIF
END IF
!
! Read in terrain data
!
call read_10x10_topo30s('tmp1.file',h_block,ll_lon, ll_lat, &
lon_inc, lat_inc,imissing)
write(6,'(4(f20.14),3i7)') ll_lon, ll_lat, lon_inc, lat_inc, imissing
!
CALL unixcmd('/bin/rm tmp1.file')
istatus = 0
!
RETURN
END SUBROUTINE read_30s_trn
SUBROUTINE read_10x10_topo30s(filename,it1,ll_lon,ll_lat, & 1
lon_inc,lat_inc,imissing)
!
! Read in 10x10 degrees 30 second encoded ASCII terrain file filename
! and store the values in it1
!
! Author: Ming Xue (1/31/2002)
!
IMPLICIT none
character (len=*) :: filename
INTEGER :: it1(1200,1200)
Doubleprecision :: ll_lon,ll_lat,lon_inc,lat_inc
INTEGER :: imissing
character(len=1), ALLOCATABLE :: ichar1(:,:)
character(len=1), ALLOCATABLE :: ichar2(:,:)
INTEGER :: itopomax, itopomin
INTEGER :: itopomax1,itopomin1
INTEGER :: i,j,ii
real :: topomean1
ALLOCATE(ichar1(1200,1200),ichar2(1200,1200))
open(21,file=trim(filename),status='old',form='formatted')
read(21,'(4(f20.14),3i7)')ll_lon,ll_lat,lon_inc,lat_inc, &
itopomin,itopomax,imissing
do j=1,1200
do i=1,1200,20
read(21,'(40a1)') (ichar1(i-1+ii,j),ichar2(i-1+ii,j),ii=1,20)
enddo
enddo
close(21)
do j=1,1200
do i=1,1200
it1(i,j) = (ichar(ichar1(i,j))-32)*110+ ichar(ichar2(i,j))-32 - 999
if( it1(i,j) .eq. -999) it1(i,j) = -9999
enddo
enddo
itopomin1 = 100000
itopomax1 = -100000
topomean1 = 0.0
do j=1,1200
do i=1,1200
if( it1(i,j).ne.-9999 .and. it1(i,j).gt.itopomax1) itopomax1=it1(i,j)
if( it1(i,j).ne.-9999 .and. it1(i,j).lt.itopomin1) itopomin1=it1(i,j)
if( it1(i,j).ne.-9999) topomean1 = topomean1 + it1(i,j)
enddo
enddo
topomean1 = topomean1/(1200*1200)
if( itopomin1 .eq. 100000 ) itopomin1 = -9999
if( itopomax1 .eq.-100000 ) itopomax1 = -9999
write(6,'(a,2i7,f20.14)') &
'Min., Max and mean of the 10x10 degree block are ', &
itopomin1, itopomax1, topomean1
IF( itopomin1 .ne. itopomin1 .or.itopomax1.ne.itopomax ) then
print*,'Warning: min or max in data read in do not agree with '
print*,'those in the header. Min., Max in thedata header were ', &
itopomin, itopomax
ENDIF
DEALLOCATE(ichar1,ichar2)
END SUBROUTINE read_10x10_topo30s
SUBROUTINE plot_trn(z,x,y,m,m1,m2,n,n1,n2,x1,x2,xstep, & 1,2
y1,y2,ystep,zinc,tem,xw,yw)
IMPLICIT NONE
INTEGER :: m,m1,m2,n,n1,n2
REAL :: x1,x2,xstep,y1,y2,ystep
REAL :: z(m,n),x(m,n),y(m,n),cl(100)
REAL :: tem(m,n),xw(8*m),yw(8*m)
REAL :: zinc,zmin, zmax, amax, amin
CHARACTER (LEN=80) :: ch
INTEGER :: lch, mode, ncl
REAL :: pl, pr, pb, pt, px, py, pxc, pyc, xs, ys
!
! Setup device and define plotting space.
!
pl = 0.15
pr = 0.85
pb = 0.15
pt = 0.85
px = pr - pl
py = pt - pb
pxc = (pr+pl)/2
pyc = (pb+pt)/2
xs = x2-x1
ys = y2-y1
IF( py/px >= ys/xs ) THEN
py = ys/xs*px
CALL xpspac(pl, pr, pyc-py/2, pyc+py/2 )
ELSE
px = xs/ys*py
CALL xpspac(pxc-px/2, pxc+px/2, pb, pt)
END IF
CALL xmap(x1,x2,y1,y2)
!
! Call XCONTA and XCLIMT to plot a contour map with default options.
!
cl(1)=0.0
cl(2)=zinc
CALL xnctrs(10,100 )
mode=1
PRINT*,'m,m1,m2, n,n1,n2=' , m,m1,m2, n,n1,n2
CALL a3dmax0(z,1,m,m1,m2,1,n,n1,n2,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'Min. terrain= ', amin,', Max. terrain=',amax
CALL xctrclr(46,68)
CALL xctrlim(0.0, amax)
CALL xctrlim(0.0, 0.0)
! call xcontcopt(2)
CALL xcolfil( tem,xw,yw,m,m2-m1+1,n2-n1+1, cl, ncl,mode)
CALL xchsiz( 0.025 * (y2-y1))
CALL xcpalet(2,-9999.9, -9999.0)
CALL xclfmt('(I4)')
CALL xctrclr(1,1)
! CALL XCONTA(Z(m1,n1),X(m1,n1),Y(m1,n1),tem,M,m2-m1+1,n2-n1+1,
! : CL, NCL,MODE)
!
zmax=cl(ncl)
zmin=cl(1)
zinc=cl( MIN(2,ncl) )-cl(1)
CALL xchsiz(0.03*(y2-y1))
WRITE(ch,'(''HMIN='',F6.1,'' HMAX='',F6.1,'' INC='',F5.1)') &
amin, amax, zinc
lch = 80
CALL strlnth(ch,lch)
CALL xcolor(1)
CALL xcharc((x1+x2)*0.5,y2+0.01*(y2-y1),ch(1:lch))
!
! Draw a border and plot ticks and labels on the border.
!
CALL xaxsca(x1,x2,xstep,y1,y2,ystep)
RETURN
END SUBROUTINE plot_trn
SUBROUTINE smth(a,m,m1,m2,n,n1,n2, ismooth,tem) 1
! Two-Dimension smooth use nine-points
IMPLICIT NONE
INTEGER :: m,n,m1,m2,n1,n2
INTEGER :: ismooth
INTEGER :: i,j,im,ip,jm,jp
REAL :: s
REAL :: a(m,n)
REAL :: tem(m,n)
IF(ismooth == 0) RETURN
DO j=1,n
DO i=1,m
tem(i,j) = a(i,j)
END DO
END DO
IF(ismooth == 1) THEN ! 5 point smoothing
PRINT*,'Performing 5-point smoothing on input array'
s=0.5
DO j=n1,n2
DO i=m1,m2
im=MAX(m1,i-1)
ip=MIN(m2,i+1)
jm=MAX(n1,j-1)
jp=MIN(n2,j+1)
tem(i,j)=a(i,j)*(1.-s) + &
( a(ip,j)+a(im,j)+a(i,jp)+a(i,jm) )*0.25*s
END DO
END DO
ELSE IF(ismooth == 2) THEN ! 9 point smoothing
PRINT*,'Performing 9-point smoothing on input array'
s=0.5
DO j=n1,n2
DO i=m1,m2
im=MAX(m1,i-1)
ip=MIN(m2,i+1)
jm=MAX(n1,j-1)
jp=MIN(n2,j+1)
tem(i,j)=a(i,j)*(1.-s)**2 + &
( a(ip,j)+a(im,j)+a(i,jp)+a(i,jm) )*0.5*s*(1.-s) &
+ ( a(im,jm)+a(im,jp)+a(ip,jm)+a(ip,jp) ) &
*0.25*s**2
END DO
END DO
ELSE IF(ismooth == 3) THEN ! horizontal smoothing only
PRINT*,'Performing horizontal smoothing on input array'
s=0.5
DO j=n1,n2
DO i=m1,m2
im=MAX(m1,i-1)
ip=MIN(m2,i+1)
tem(i,j)=a(i,j)*(1.-s)+(a(ip,j)+a(im,j))*0.5*s
END DO
END DO
ELSE IF(ismooth == 4) THEN ! vertical smoothing only
PRINT*,'Performing vertical smoothing on input array'
s=0.5
DO j=n1,n2
DO i=m1,m2
jm=MAX(n1,j-1)
jp=MIN(n2,j+1)
tem(i,j)=a(i,j)*(1.-s)+(a(i,jp)+a(i,jm))*0.5*s
END DO
END DO
END IF
DO j=n1,n2
DO i=m1,m2
a(i,j) = tem(i,j)
END DO
END DO
RETURN
END SUBROUTINE smth