PROGRAM dir5min,1

!##################################################################
!##################################################################
!######                                                      ######
!######                 PROGRAM DIR5MIN                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!
!  PURPOSE:
!
!  A stand alone program which rewrites 5 minute 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 dir5min.hdr file (header file) will be created and
!  used by arpstern11.f.
!
!  INPUT:
!            elev.dat          (NCAR DATA FILE 5 MINUTE RESOLUTION)
!
!            arpstern.input    ARPSTERN11.F input file
!
!
!
!
!  OUTPUT:
!            dir5min.dat       direct access 5 minute unformatted
!                              data file with 4118 seperate records.
!
!            dir5min.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 5 minute 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=12)

  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 to the record number
  INTEGER :: flon(jlim)       !  Longitude corresponding to the 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)       !  Latitude
  INTEGER :: nlon(ilon)       !  Longitude

  INTEGER :: TYPE(klim,ilim)  ! Type of area the block defines.  See
                              !  appendix for definition of type.

  INTEGER :: relev(klim,ilim) ! Elevation value read in from the elev.dat
                              !  data file

  INTEGER :: ielev(jlim,ndx,ndx) ! Intermediate terrain block data put
                                 !  in terms of increasing record number order

!
!  LOCAL VARIABLES:
!

  INTEGER :: irec,i,ii,j,k,kk,l,num,iloc
  INTEGER :: item1,item2,item3
  INTEGER :: comtype          !  Computer type for program to run on...
                              !  for IBM, comtype = 1
                              !     CRAY,         = 4
  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
  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
!
!  OUTPUT:
!
  INTEGER*2 iw(ndx,ndx)    !  Terrain data written in unformatted
                           !  form for entire 1 degree by 1 degree
                           !  block to dir5min.f file.

!-----------------------------------------------------------------------
!
!  Defining namelist:
!
!-----------------------------------------------------------------------

  NAMELIST /terraind/analtype,mapproj,itertype,rmsopt,                  &
      comtype,tdatadir,terndir

!-----------------------------------------------------------------------
!
!  Begin executable portion of the code.
!
!-----------------------------------------------------------------------

  READ(5,terraind)
!
!  open elev.dat, dir5min.dat and dir5min.hdr files...
!

  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='dir5min.dat',ACCESS='direct',FORM='unformatted',        &
      STATUS='unknown',RECL=288*comtype)
  OPEN(12,FILE='dir5min.hdr',FORM='formatted',STATUS='unknown')

!
!  Initializing counters.....
!

  irec=0            ! number of 5 minute data records read
  num=1
  nlat(num)=74      ! no 5 minute data above this latitude (deg north)
  nlon(num)=0       ! starting longitude for each latitude
  item3=0

!
!    Starting the main loop.  Note that the entire data file is read
!    and the data stored and then the data is written to file outside
!    this loop....
!

  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,relev(1,i)   &
        ,TYPE(2,i),dum1,relev(2,i),TYPE(3,i),dum1,relev(3,i),TYPE(4,i)  &
        ,dum1,relev(4,i),TYPE(5,i),dum1,relev(5,i),TYPE(6,i),dum1       &
        ,relev(6,i),dum4,qsize(i),dum3,i=1,ilim)

!
!  test to see if the block type is 5 minute
!

    IF(qsize(1) /= 1)THEN !qsize =1 for 1 degree, =5 for 5 minute data

!
!  Test to find the character of the data (ocean, land...)
!  and set ocean values to zero elevation.
!

      DO k=1,klim
        DO i=1,ilim
          IF(TYPE(k,i) == 9)THEN ! quad is ocean and the elev is set to zero.
            relev(k,i)=200       ! 200(packed) = 0.0 meters unpacked

          ELSE IF(TYPE(k,i) == 0.OR.TYPE(k,i) == 4)THEN !quad contains land/ice
            relev(k,i)=(relev(k,i)+4000)/20          ! putting in packed form

          ELSE IF(TYPE(k,i) == 8)THEN       ! quad contains some negative land
            relev(k,i)=(relev(k,i)+4000)/20          ! putting in packed form

          ELSE IF(TYPE(k,i) == 5.OR.TYPE(k,i) == 3.OR.                  &
                    TYPE(k,i) == 2)THEN                  !land/ocean mix
            IF(relev(k,i) < 0)THEN
              relev(k,i)=200                         ! ice/ocean mix
            ELSE
              relev(k,i)=(relev(k,i)+4000)/20        ! putting in packed form
            END IF
          ELSE                               ! type(k,i) does not match
                                             ! defined list

            PRINT *,'type of elevation unknown,irec=',irec
            STOP
          END IF
        END DO
      END DO

!
!  Sort the data into southwest to northeast order from northwest to
!  southeast order.
!
!  Note the data is read from the northwest corner of the 1 degree
!  blocks.  This data will be sorted and stored with the (1,1) element
!  as the southwest corner of the block similar to the 30 second data
!  set format.
!

      DO i=1,ilim
        IF(item3 /= rlat(i) ) THEN

        IF( rlat(i) == 99999 ) CYCLE
 
          item1=MOD(rlat(i),100)
          item2=MOD(rlon(i),100)
          IF(item1 == 0.AND.item2 == 0)THEN    ! starting a new block...
            kk=ndx
            irec=irec+1
            flat(irec)=(rlat(i)-100)/100      ! adjusting the latitude...
            flon(irec)=rlon(i)/100            ! setting the longitude...
            DO k=1,ndx/2
              ielev(irec,k,kk)=relev(k,i)
            END DO

            IF(nlat(num) == flat(irec))THEN   ! counting the number of blocks
              nlon(num)=nlon(num)+1           ! for a latitude circle...
            ELSE
              num=num+1
              nlat(num)=flat(irec)
              nlon(num)=1
            END IF
          ELSE IF(item1 == 0.AND.item2 /= 0)THEN ! setting an existing block..
            DO k=ndx/2+1,ndx
              kk=ndx
              ielev(irec,k,kk)=relev(k-6,i)
            END DO
          ELSE IF(item1 /= 0.AND.item2 == 0)THEN ! setting an existing block..
            kk=kk-1
            DO k=1,ndx/2
              ielev(irec,k,kk)=relev(k,i)
            END DO
          ELSE IF(item1 /= 0.AND.item2 /= 0)THEN ! setting an existing block..
            kk=kk-1
            DO k=ndx/2+1,ndx
              ielev(irec,k,kk)=relev(k-6,i)
            END DO
          ELSE
            PRINT *,'lat,lon does not follow format',rlat(i),rlon(i)
          END IF
          item3=rlat(i)
        ELSE
          PRINT *,'duplicate found at',rlat(i),rlon(i)
        END IF
      END DO ! i loop

    END IF

  END DO  ! j loop 

!
!  End of read/sort loop (main loop)
!

!
!  write 5 minute header to the dir5min.hdr file and the 5 minute data
!  to the dir5min.dat file.  Note: the latitude is for the northwest
!  corner of the block, this is inconsistant with the dir30sec.dat storage
!  format.  The 5 minute data will be written following the 30 second
!  format with the southwestern most data block written first and then
!  the next block to the east and so on.  The next latitiude to the north
!  will then be written.
!

  99   PRINT *,'print main loop index, and irec',jlim,irec
  iloc=0
  kk=0

  DO i=num,1,-1
    iloc=iloc+nlon(i)
    PRINT *,'latitude, number of blocks on the latitude circle, ',      &
            'running sum of block written', nlat(i),nlon(i),iloc
    DO j=1,nlon(i)
      kk=kk+1
      ii=irec-iloc+j
      DO l=1,ndx
        DO k=1,ndx
          iw(k,l)=ielev(ii,k,l)
        END DO
      END DO

!
!   Writing the block data to the unformatted file and the header
!   number corresponding to the block location to the header file.
!

      WRITE(12,'(i6,i6,i6)') flat(ii),flon(ii),kk
      WRITE(11,REC=kk) iw

!
!  Testing for negative values...
!

      DO l=1,ndx
        DO k=1,ndx
          IF(iw(k,l) < 200)PRINT *,'negative elev.',flat(ii),flon(ii),kk
        END DO
      END DO

    END DO
  END DO

!
!  Closing all files and exiting the program....
!

  CLOSE(10)
  CLOSE(11)
  CLOSE(12)

  STOP
END PROGRAM dir5min
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                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