!
!           PROGRAMA DE LEITURA de arquivos tipo ETAGRDhh.tmhh
!
      program leit
!      use variables

      integer, parameter         :: kmax=50

      integer               :: dt,kgtype,imdlty,igout,jgout,im,jm,NRSTRT,NPINCR
      real                  :: dx,lonw,alatvt,lats,dy
      character(len=20)     :: field
      character(len=6)      :: nufile,outype,proj,readco,readll,itag,dataset,datset
      character(len=12)     :: compl
      logical               :: ltsoil, lqsoil,north
      logical               :: infloop

!
! Dimension 2D variables
!
      real,allocatable,dimension(:,:)      :: lon       !LONGITUDE
      real,allocatable,dimension(:,:)      :: lat       !LATITUDE
      real,allocatable,dimension(:,:)      :: mask      !LAND/SEA MASK
      real,allocatable,dimension(:,:)      :: cssf      !AVE SFC SENHEAT FX
      real,allocatable,dimension(:,:)      :: cmsf      !AVE SFC MOMENTUM FX
      real,allocatable,dimension(:,:)      :: ps        !SURFACE PRESSURE
      real,allocatable,dimension(:,:)      :: ttprec    !ACM TOTAL PRECIP
      real,allocatable,dimension(:,:)      :: zorl      !ROUGHNESS LENGTH
      real,allocatable,dimension(:,:)      :: GRIB
      real,allocatable,dimension(:,:)      :: tsh       !SHELTER TEMPERATURE
      real,allocatable,dimension(:,:)      :: zs        !SURFACE HEIGHT
      real,allocatable,dimension(:,:)      :: tdsh      !SHELTER DEWPOINT
      real,allocatable,dimension(:,:)      :: uanem     !U WIND AT ANEMOM HT
      real,allocatable,dimension(:,:)      :: vanem     !V WIND AT ANEMOM HT
						
!
!
! Dimension 3D variables
!

      real,allocatable,dimension(:,:,:) :: phi   !HEIGHT ON ETA SFCS
      real,allocatable,dimension(:,:,:) :: peta  !PRESS ON ETA SFCS
      real,allocatable,dimension(:,:,:) :: q      !SPEC HUM ON ETA SFCS
      real,allocatable,dimension(:,:,:) :: tpot  !POT TEMP ON ETA SFCS
      real,allocatable,dimension(:,:,:) :: tke   !TRBLNT KE ON ETA SFC
      real,allocatable,dimension(:,:,:) :: u     !U WIND ON ETA SFCS
      real,allocatable,dimension(:,:,:) :: v     !V WIND ON ETA SFCS

      infloop=.true.
      read(5,*)ITAG,NRSTRT,NPINCR,dataset
!GSM	  dataset=TRIM(dataset)
      write(6,*)ITAG,NRSTRT,NPINCR,dataset

!     READ OUTPUT GRID SPECIFICATIONS.
!
      OPEN(18,FILE='cntrl.parm_'//dataset,status='old')
      rewind(18)
      read(18,1000) kgtype
      read(18,1000) imdlty
      read(18,1030) datset
      read(18,1030) outype
      read(18,1030) nufile
      read(18,1030) proj
      read(18,1010) north
      read(18,1000) im
      read(18,1000) jm
      read(18,1020) dx
      read(18,1020) lonw
      read(18,1020) alatvt
      read(18,1020) lats
      read(18,1020) dy
      read(18,1030) readll
      read(18,1030) readco
 1000 format(T28,I5)
 1010 format(T28,L1)
 1020 format(T28,F11.6)
 1030 format(T28,A6)
      print*,"im: ",im," jm: ",jm,"kmax: ",kmax


! Allocate variables
      allocate(GRIB(im,jm))
      allocate(lon(im,jm))
      allocate(lat(im,jm))
      allocate(tsh(im,jm))
      allocate(tdsh(im,jm))
      allocate(zs(im,jm))
      allocate(uanem(im,jm))
      allocate(vanem(im,jm))
      allocate(mask(im,jm))
      allocate(ps(im,jm))
      allocate(ttprec(im,jm))
      allocate(cssf(im,jm))
      allocate(cmsf(im,jm))
      allocate(zorl(im,jm))


!       Variaveis 3D
      allocate(peta(im,jm,kmax))
      allocate(phi(im,jm,kmax))
      allocate(tpot(im,jm,kmax))
      allocate(q(im,jm,kmax))
      allocate(u(im,jm,kmax))
      allocate(v(im,jm,kmax))
      allocate(tke(im,jm,kmax))


      OPEN(unit=51,file=datset//itag//'.T00S',status='old',form='unformatted')

! get I and J indexes for the smaller domain
!      IW= ABS((-83) - (-83))/ 0.40  + 1
!      IE= ABS((-83) - (-25.8))/ 0.40 + 1
!      JS= ABS((-50.2) - (-50.2))/ 0.40 + 1
!      JN= ABS((-50.2) - (12.2))/ 0.40 + 1
      IW= 1
      IE= im
      JS= 1
      JN= jm

!
! get the nx and ny  number of points in x and y direction
      nx = ie - iw + 1
      ny = jn - js + 1
!
! get the exactly lat/lon of the smaller domain
!
      clonll= -lonw + ((iw-1) * dx)
      clonur= -lonw + ((ie-1) * dx)
      clatll= lats + ((js-1) * dy)
      clatur= lats + ((jn-1) * dy)

      print*, iw, ie, js, jn
      print*, clonll, clonur, clatll, clatur

      READ(51)   IHRST,IMM,IDD,IYY,IHH
      print*, ' reading header'
      READ(51)   KGTYP,PROJ,NORTH,IMOUT,JMOUT,POLEI,POLEJ,ALATVT,ALONVT,XMESHL
      kqs=0
      kts=0
      kt=0
      kq=0
      kr=0
      kp=0
      kw=0
      ku=0
      kv=0
      kc=0
      ltsoil=.false.
      lqsoil=.false.

      do while (infloop)
        READ(51,END=999) FIELD,SFC
        print*, field, sfc
        READ(51) GRIB
        if (field.eq.    'LONGITUDE           ') then; lon=grib 
	elseif (field.eq.     'LATITUDE            ') then; lat=grib
	elseif (field.eq.     'SURFACE HEIGHT      ') then; zs=grib
	elseif (field.eq.     'SHELTER TEMPERATURE ') then; tsh=grib
	elseif (field.eq.     'SHELTER DEWPOINT    ') then; tdsh=grib
	elseif (field.eq.     'U WIND AT ANEMOM HT ') then; uanem=grib
 elseif (field.eq.     'V WIND AT ANEMOM HT ') then; vanem=grib
	elseif (field.eq.     'LAND/SEA MASK       ') then; mask=grib
	elseif (field.eq.     'SURFACE PRESSURE    ') then; ps=grib
        elseif (field.eq.'ACM TOTAL PRECIP    ') then; ttprec=grib/1000
        elseif (field.eq.'AVE SFC SENHEAT FX  ') then; cssf=grib
        elseif (field.eq.'AVE SFC MOMENTUM FX ') then; cmsf=grib
        elseif (field.eq.'ROUGHNESS LENGTH    ') then; zorl=grib
                         
!
!  retrieving 3d fields
!                        
        elseif (field.eq.'PRESS ON ETA SFCS   ') then; kw=kw+1; peta(1:im,1:jm,kw)=grib
        elseif (field.eq.'HEIGHT ON ETA SFCS  ') then; kp=kp+1; phi(1:im,1:jm,kp)=grib
        elseif (field.eq.'POT TEMP ON ETA SFCS') then; kt=kt+1; tpot(1:im,1:jm,kt)=grib
        elseif (field.eq.'SPEC HUM ON ETA SFCS') then; kq=kq+1; q(1:im,1:jm,kq)=grib
        elseif (field.eq.'U WIND ON ETA SFCS  ') then; ku=ku+1; u(1:im,1:jm,ku)=grib
        elseif (field.eq.'V WIND ON ETA SFCS  ') then; kv=kv+1; v(1:im,1:jm,kv)=grib
        elseif (field.eq.'TRBLNT KE ON ETA SFC') then; kc=kc+1; tke(1:im,1:jm,kc)=grib
        endif
      enddo
 999  continue


!
!
!  save fields for Grads plotting. First 2D, then 3D fields
!
      OPEN(unit=17,file='latlon_'//itag,STATUS='UNKNOWN', &
     &        form='unformatted',access='direct' ,           &
     &        recl=4*nx*ny)

       irec=1
       write(17,rec=irec)((lon(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((lat(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((zs(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((tsh(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((tdsh(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((uanem(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((vanem(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((mask(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((ps(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       where(ttprec<=0.) ttprec=0.
       write(17,rec=irec)((ttprec(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((cssf(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((cmsf(i,j),i=iw,ie),j=js,jn)
       irec=irec+1
       write(17,rec=irec)((zorl(i,j),i=iw,ie),j=js,jn)


      print*, ' 2D irec=',irec-1
      irec=irec+1
      do k=kmax,1,-1
       write(17,rec=irec)((peta(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((phi(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((tpot(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((q(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((u(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((v(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      do k=kmax,1,-1
       write(17,rec=irec)((tke(i,j,k),i=iw,ie),j=js,jn)
       irec=irec+1
      enddo

      stop
      end program leit

