      program read_gfs

      implicit none
      character(len=255)              :: gribfile,outdir
!_____________________________________________________________________________________________________________
      integer                         :: nx,ny,nz,np
      integer                         :: dx,dy 
      integer                         :: proj,LatS,LonW,LatN,LonE
      integer                         :: i,j,k,l,n,n2

      REAL, PARAMETER                 :: Undef=9.9990003E+20 
      real,    dimension(:,:,:) ,allocatable    :: ht                        !Isobaric heights (m)
      real,    dimension(:,:,:) ,allocatable    :: tp                        !Isobaric temps (K)
      real,    dimension(:,:,:) ,allocatable    :: th                        !Isobaric theta (K)
      real,    dimension(:,:,:) ,allocatable    :: uw                        !Isobaric u-wind (m/s)
      real,    dimension(:,:,:) ,allocatable    :: vw                        !Isobaric v-wind (m/s)
      real,    dimension(:,:,:) ,allocatable    :: rh                        !Isobaric rh,mr (%,kg/kg)
      real,    dimension(:)     ,allocatable    :: pr,pri                    !Isobaric pressures (mb)
!_____________________________________________________________________________________________________________
      real,    dimension(:,:,:) ,allocatable    :: slp
                                                                             !slp(i,j,1)=LCMSK;slp(i,j,2)=PSLM
                                                                             !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC
                                                                             !slp(i,j,5 & 6) are 0-10 cm STC and SMC
                                                                             !slp(i,j,7 & 8) are 10-200 cm STC and SMC
!_____________________________________________________________________________________________________________
      character(len=15)               :: ModelDriver
!_____________________________________________________________________________________________________________

      open(1,file='InputModelInf.txt',form='formatted',status='old')
      rewind 1
      read(1,'(a15)')ModelDriver
      read(1,*)nx,ny,nz
      read(1,*)proj
      read(1,*)LatS
      read(1,*)LonW
      read(1,*)LatN
      read(1,*)LonE
      read(1,*)dx
      read(1,*)dy

      allocate(ht(nx,ny,nz))
      allocate(tp(nx,ny,nz))                       
      allocate(th(nx,ny,nz))                     
      allocate(uw(nx,ny,nz))                 
      allocate(vw(nx,ny,nz))                
      allocate(rh(nx,ny,nz)) 
      allocate(slp(nx,ny,12)) 
      allocate(pr(nz))
      allocate(pri(nz)) 

      read(5,'(a)') gribfile

      do k=1,nz
        read(1,*)pr(k)
        print*,pr(k)
      enddo
      close(1)

! *** Fill pressure levels.

        write(6,*) 'using hardwired stuff! ', nx,ny,nz
        write(6,*) (pr(k),k=1,nz)

! *** Read in degrib data.
!

        n=index(outdir,' ')-1
        n2=index(ModelDriver,' ')-1
        write(6,*)n2
        write(6,*)ModelDriver(1:n2)
        print*, "MODELO UTILIZADO=",ModelDriver(1:n2)
        
      n=index(gribfile,' ')-1
      write(*,*) 'OPEN FILE ',gribfile(1:n)

      open(12,file=gribfile(1:n),form='unformatted',status='old')
						
!!!  Dragan 12/07/2018    open(15,form='unformatted', access='direct',recl=1440*721*4)
 
      open(15,form='unformatted', access='direct',recl=nx*ny*4)
      
      do l=1,nz

        read(12) ht(:,:,l)
        print*,' ht l ',l,ht(481,241,l),ht(241,481,l)

      enddo

      do l=1,nz

        read(12) uw(:,:,l)
        print*,'uw l ',l,uw(481,241,l),ht(241,481,l)

      enddo

      do l=1,nz

        read(12) vw(:,:,l)
        print*,'vw l ',l, vw(481,241,l),ht(241,481,l)

      enddo

      do l=1,nz

        read(12) tp(:,:,l)
        print*,'tp l ',l, tp(481,241,l),ht(241,481,l)

      enddo

      do l=1,nz

        read(12) rh(:,:,l)
        print*,'rh l ',l, rh(481,241,l),ht(241,481,l)

      enddo

!   Missing RH at 20hPa 
!   
!      rh(:,:,nz)=-99999.
!      rh(:,:,nz-1)=-99999.
       
      do l=1,12

         if (l.eq.2.or.l.eq.3) then

             read(12) slp(:,:,l)
             slp(:,:,l)=slp(:,:,l)/100

          else

            read(12) slp(:,:,l)

          endif

       print*,'slp l ',l, slp(481,241,l),ht(241,481,l)
read_gfs.f90
      enddo
      close(12)

!#############################################################################################################

       do i=1,nx
           do j=1,ny
              do l=6,12
                 if(l.eq.6.or.l.eq.8.or.l.eq.10.or.l.eq.12) then

                      if (slp(i,j,l).eq.Undef) slp(i,j,l)=0
                      if (slp(i,j,l).gt.0.5) slp(i,j,l)=0.3
 
                  endif
              enddo
           enddo
         enddo

!#############################################################################################################

!------Over the Ocean, the value of Temp/Soil data are equal 9.9990003E+20 leaving shadows at the edges,
!------so, to solve the problem we fill with neighbors value


       do l=5,11
          if(l.eq.5.or.l.eq.7.or.l.eq.9.or.l.eq.11) then

              do np=1,200
                 do I=1,nx
                   do J=1,ny
                      if ( slp(I,J,l).eq.Undef.and.slp(I-1,J,l).ne.Undef) then
                           slp(I,J,l)=slp(I-1,J,l)
                      endif
                   enddo
                 enddo

                 do I=nx-1, 1, -1
                   do J=1,ny
                      if ( slp(I,J,l).eq.Undef.and.slp(I+1,J,l).ne.Undef) then
                           slp(I,J,l)=slp(I+1,J,l)
                      endif
                   enddo
                 enddo
              enddo


              do np=1,200
                 do I=1,nx
                   do J=1,ny
                      if ( slp(I,J,l).eq.Undef.and.slp(I,J-1,l).ne.Undef) then
                           slp(I,J,l)=slp(I,J-1,l)
                      endif
                   enddo
                 enddo


                 do I=1,nx
                   do J=ny-1, 1, -1
                      if ( slp(I,J,l).eq.Undef.and.slp(I,J+1,l).ne.Undef) then
                           slp(I,J,l)=slp(I,J+1,l)
                      endif
                   enddo
                 enddo
              enddo
        endif
       enddo

!!##################################################  WRITE ##################################################

      do l=1,nz
         write(15) uw(:,:,l)
        print*,' uww l ',l,uw(481,241,l),uw(241,481,l)
      enddo


      do l=1,nz
         write(15) vw(:,:,l)
        print*,' vww l ',l,vw(481,241,l),vw(241,481,l)
      enddo


      do l=1,nz
         write(15) ht(:,:,l)
        print*,' htw l ',l,ht(481,241,l),ht(241,481,l)
      enddo


      do l=1,nz
         write(15) tp(:,:,l)
        print*,' tpw l ',l,tp(481,241,l),tp(241,481,l)
      enddo


      do l=1,nz
         write(15) rh(:,:,l)
        print*,' rhw l ',l,rh(481,241,l),rh(241,481,l)
      enddo


      do l=1,12
         write(15) slp(:,:,l)
        print*,' slpw l ',l,slp(481,241,l),slp(241,481,l)
      enddo


      end program read_gfs

