PROGRAM dir1deg,1
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM DIR1DEG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!
! PURPOSE:
!
! A stand alone program which rewrites 1 degree terrain elevation data
! from elev.dat recieved from NCAR DATA SERVICES into a direct access
! file using either 2 of 8 byte integer format depending upon the type
! of machine used. The purpose of writing a direct access data format
! is to expedite the run time of the ARPSTERN11.F terrain processing
! program for initializing the ARPS terrain field.
! In addition, a dir1deg.hdr file (header file) will be created and
! used by arpstern11.f.
!
! INPUT:
! elev.dat (NCAR DATA FILE 5 MINUTE AND 1 DEGREE
! RESOLUTION)
!
! arpstern.input ARPSTERN11.F input file
!
!
!
!
! OUTPUT:
! dir1deg.dat direct access 1 degree unformatted
! data file with xxxx seperate records.
!
! dir1deg.hdr header file which has record numbers
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber
! 1/12/94
!
! Modification history:
!
! Dan Weber 7/06/94.
! Documentation added.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!
! INPUT:
!
INTEGER :: ndx ! Number of data points in the longitudinal
! and latitudinal direction
INTEGER :: jlim ! Limit of the number of reads in the
! elev.dat data file
INTEGER :: klim ! Number of reads at a given
! longitude (determined by elev.dat = 6)
INTEGER :: ilim ! Number of 80 character lines per record
! in elev.dat (=20)
INTEGER :: ilon ! Maximum longitude
PARAMETER (ilim=20,ilon=360,jlim=100000,klim=6,ndx=2)
INTEGER :: rlat(ilim) ! Latitude read in from the ascii data file
INTEGER :: rlon(ilim) ! Longitude read in from the ascii data file
INTEGER :: flat(jlim) ! Latitude corresponding a record number
INTEGER :: flon(jlim) ! Longitude corresponding a record number
INTEGER :: qsize(ilim) ! Quad size of the block (1= 1 deg by 1 deg and
! 30' by 30', 5 = 5' by 5')
INTEGER :: nlat(ilon) ! Search Latitude
INTEGER :: nlon(ilon) ! Search Longitude
!
! LOCAL VARIABLES:
!
INTEGER :: irec,i,ii,j,k,kk,l,num,iloc
CHARACTER (LEN=80 ) :: tdatadir ! Directory in which the data file
! dma_elev.dat is stored.
CHARACTER (LEN=80 ) :: terndir ! Directory in which the original ASCII terrain
! data files are stored.
INTEGER :: lterndir ! Length of the non-blank part of string terndir.
INTEGER :: analtype ! Dummy variable used for namelist data
INTEGER :: mapproj ! Dummy variable used for namelist data
INTEGER :: itertype ! Dummy variable used for namelist data
INTEGER :: rmsopt ! Dummy variable used for namelist data
INTEGER :: comtype ! Computer type for program to run on...
! for IBM, comtype = 1
! CRAY, = 4
CHARACTER (LEN=1) :: dum1 ! Dummy variable used to read ascii header
CHARACTER (LEN=2) :: dum2 ! Dummy variable used to read ascii header
CHARACTER (LEN=8) :: dum3 ! Dummy variable used to read ascii header
CHARACTER (LEN=9) :: dum4 ! Dummy variable used to read ascii header
INTEGER :: units(ilim) ! Units of measure for the terrain data
! units = 1 meters
! = 2 feet
! = 3 fathoms
INTEGER*2 TYPE(klim,ilim) ! Type of area the block defines See
! appendix fro definition of type.
INTEGER :: ielev(klim,ilim) ! Intermediate terrain block data put in
! terms of increasing record number order.
INTEGER*2 itp(4,jlim) ! Intermediate terrain values in record
! number order
!
! OUTPUT:
!
INTEGER*2 iw(ndx,ndx) ! Terrain data written in unformatted
! form for entire 1 degree by 1 degree
! block to dir1deg.f file.
!
!-----------------------------------------------------------------------
!
! Defining namelist:
!
!-----------------------------------------------------------------------
NAMELIST /terraind/analtype,mapproj,itertype,rmsopt,comtype, &
tdatadir,terndir
!-----------------------------------------------------------------------
!
! Begin executable portion of the code.
!
!-----------------------------------------------------------------------
!
! Read the arpstern.input file for the terraind namelist information.
! (note: only need comtype for this program to run properly.)
READ(5,terraind)
lterndir = LEN(terndir)
CALL strlnth
( terndir, lterndir )
IF( lterndir == 0 ) THEN
terndir = '.'
lterndir=1
END IF
!
! open elev.dat, dir1deg.dat and dir1deg.hdr files...
!
OPEN(10,FILE=terndir(1:lterndir)//'/elev.dat',STATUS='old')
OPEN(11,FILE='dir1deg.dat',ACCESS='direct',FORM='unformatted', &
STATUS='unknown',RECL=8*comtype)
OPEN(12,FILE='dir1deg.hdr',FORM='formatted',STATUS='unknown')
!
! Initializing counters....
!
irec=0
num=1
nlat(1)=89 ! latitude of the southwestern part of the northern
! most block of data.
nlon(1)=0
!
! Starting the main loop which reads the elev.dat data set
! and check to see if the data is 5 minute or 1 degree (or
! 30 minute) resolution.
!
DO j=1,jlim
!
! Reading 1 record at a time. (note: each record consists of
! 20 lines with 80 charactes...)
!
READ(10,'(20(A2,I5,A1,I5,I1,6(i1,a1,i6),A9,I1,a8))',END=99) &
(dum2,rlat(i),dum1,rlon(i),units(i),TYPE(1,i),dum1,ielev(1,i) &
,TYPE(2,i),dum1,ielev(2,i),TYPE(3,i),dum1,ielev(3,i),TYPE(4,i) &
,dum1,ielev(4,i),TYPE(5,i),dum1,ielev(5,i),TYPE(6,i),dum1 &
,ielev(6,i),dum4,qsize(i),dum3,i=1,ilim)
! print *,'ok after read',j
!
! Test for 1 degree block type
!
IF(qsize(1) == 1)THEN ! qsize=1 for 1 degree, =5 for 5 minute data
!
! write 1 degree header to the dir1deg.hdr file and the 1 degree data
! to the dir1deg.dat file. Note: the latitude is for the nortwest
! corner of the block, this is inconsistant with the dir1km.dat storage
! format. Thus, 1 will be subtracted from the latitude to provide
! consistancy. The first record will be 89 south and 0 east. the file
! will store 89 south and 1 east and so on until 89 south and 359 east
! is written. Then the next latitude to the north will be written
! and so on, all the way up to the north pole.
!
DO i=1,ilim
IF(nlat(num) == (rlat(i)/100-1))THEN !testing if rlat is same as
nlon(num)=1+nlon(num) ! previous one
ELSE
num=num+1
nlat(num)=rlat(i)/100-1
nlon(num)=1
END IF
irec=irec+1
flat(irec)=rlat(i)/100-1
flon(irec)=rlon(i)/100
IF(TYPE(1,i) == 9)THEN !quad is all ocean and the elev is set to zero
DO k=2,5
ielev(k,i)=200 ! and all components of ielev
END DO
ELSE IF(TYPE(1,i) == 0.OR.TYPE(1,i) == 4)THEN !quad contains land/ice
IF(ielev(2,i) == 0.AND.ielev(3,i) == 0.AND.ielev(4,i) == 0 &
.AND.ielev(5,i) == 0.AND.ielev(1,i) > 0)THEN ! packing data
DO k=2,5
ielev(k,i)=(ielev(1,i)+4000)/20
END DO
ELSE ! setting elevation to sea level (0.0 meters)
DO k=2,5
IF(ielev(k,i) < 0)THEN
ielev(k,i)=200
ELSE
ielev(k,i)=(ielev(k,i)+4000)/20 ! putting in packed form
END IF
END DO
END IF
ELSE IF(TYPE(1,i) == 8)THEN ! quad contains some negative land
DO k=2,5
ielev(k,i)=(ielev(k,i)+4000)/20 ! putting in packed form
END DO
ELSE IF(TYPE(1,i) == 5.OR.TYPE(1,i) == 3.OR.TYPE(1,i) == 2)THEN
! land/ocean mix and
! ice/ocean mix.
IF(ielev(2,i) == 0.AND.ielev(3,i) == 0.AND.ielev(4,i) == 0.AND. &
ielev(5,i) == 0.AND.ielev(1,i) > 0)THEN
DO k=2,5
ielev(k,i)=(ielev(1,i)+4000)/20 ! packing data...
END DO
ELSE
DO k=2,5
IF(ielev(k,i) < 0)ielev(k,i)=0
ielev(k,i)=(ielev(k,i)+4000)/20 ! putting in packed form
END DO
END IF
ELSE ! type(k,i) out of accepted values
PRINT *,'type of elevation unknown'
STOP
END IF
DO k=2,5
itp(k-1,irec)=ielev(k,i)
END DO
END DO
END IF
END DO
99 PRINT *,j,jlim,irec ! print the number of records
iloc=0
kk=0
DO i=num,1,-1
iloc=iloc+nlon(i)
PRINT *,'latitude, number of blocks found of the latitude circle' &
,nlat(i),nlon(i),iloc
DO j=1,nlon(i)
kk=kk+1
ii=irec-iloc+j
!
! putting lat in the same format as the 30 second file.
!
iw(1,1)=itp(3,ii)
iw(2,1)=itp(4,ii)
iw(1,2)=itp(1,ii)
iw(2,2)=itp(2,ii)
WRITE(12,'(i6,i6,i6)') flat(ii),flon(ii),kk ! writing data to
! the dir1deg.hdr header file
WRITE(11,REC=kk) iw ! writing terrain data to the
! dir1deg.dat data file
DO l=1,2
DO k=1,2
IF(iw(k,l) < 200)PRINT *,flat(ii),flon(ii),kk ! printing negative
! value locations....
END DO
END DO
END DO
END DO
STOP
END PROGRAM dir1deg
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE STRLNTH ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE strlnth( string, length ) 176
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Return the length of the non-blank part of a character string.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/20/91
!
! MODIFICATION HISTORY:
!
! 5/05/92 (M. Xue)
! Added full documentation.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! string A character string
! length The declared length of the character string 'string'.
!
! OUTPUT:
!
! length The length of the non-blank part of the string.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=* ) :: string ! A character string for the name of
! this run.
INTEGER :: length ! The length of the non-blank part
! of a string.
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i = length,1,-1
IF(string(i:i) /= ' ') EXIT
END DO
! 200 CONTINUE
length = MAX(1,i)
RETURN
END SUBROUTINE strlnth