!-------------------------------------------------------------------------
!     NASA GSFC Land Information Systems LIS 2.3                         !
!-------------------------------------------------------------------------
!BOP
! !ROUTINE: retgdas.F90
!
! !DESCRIPTION:
!  Defines forcing parameters, retrieves the fields using calls to
!  getgb, and interpolates the fields to LDAS specifications
!
! !REVISION HISTORY:
!  14 Dec 2000: Urszula Jambor; Rewrote geteta.f to use GDAS forcing in GLDAS
!  15 Mar 2001: Jon Gottschalck; Added additional parameters and octets in 
!               which to search in GRIB files
!  01 Jun 2001: Urszula Jambor; Added option to get forcing from different 
!               files (F00 instantaneous and F06 six hour means)
!  29 Jan 2003: Urszula Jambor; Rewrote code, uses GETGB call to replace
!               ungribgdas.  Interpolation now occurs in interp_gdas.  
!               Using GETGB avoids problems with the Oct2002 GDAS 
!               grid update.
!  12 Nov 2003: Matt Rodell; Check to make sure input file exists before
!		opening and thereby creating a new, empty file.
!  14 Nov 2003: Matt Rodell; Ensure lugb varies in call to baopen
!
! !INTERFACE:
subroutine retgdas( order, ld, gindex, name, ferror,try )
! !USES:  
  use lis_module      ! LDAS non-model-specific 1-D variables
  use lisdrv_module, only : lis
  use time_manager
  use baseforcing_module, only: glbdata1, glbdata2
  use gdasdomain_module, only : gdasdrv
  implicit none
! !ARGUMENTS:
  type (lisdec) ld
  integer :: gindex(ld%d%lnc, ld%d%lnr)
  integer :: order    ! 1 indicates lesser interp. bdry, 2 indicates greater
  character(len=80) :: name, nameF06
  integer :: F00flag  ! if 1, need for data from 2 files (name, nameF06)
  integer :: ferror   ! set to zero if there's an error
  integer :: try
!EOP
!==== Local Variables=======================
  
  integer agrmetFlag, cmapFlag

  character(len=80) :: fname
  integer :: iv, c, r, ii, mx, my,i
  integer :: errorflag
  integer :: endloop, nforce
  integer :: j, lugb,iret,gbret,jret,jpds(200),jgds(200)
  integer :: lubi,kf,k,kpds(200),kgds(200)
  integer :: ngdas
  integer, dimension(gdasdrv%nmif) :: pds5, pds6, pds7, pds16

  logical*1, allocatable :: lb(:)
  logical :: file_exists
  real, allocatable :: f(:)
  real, allocatable :: varfield(:,:)
  real :: ism
  integer :: count

! IPOLATES VARIABLES, JESSE 20050225

  integer                :: ip, ipopt(20)
  integer                :: kgdsi(25), kgdso(25)
  integer                :: mi, mo
  integer,   parameter   :: km = 1
  integer,   parameter   :: ibi = 1
  integer                :: no
  real,      allocatable :: rlat(:)
  real,      allocatable :: rlon(:)
  integer                :: ibo
  logical*1, allocatable :: lo(:)
  real,      allocatable :: g(:)

!=== End Variable Definition =======================
!BOC
  ngdas = (gdasdrv%ncold*gdasdrv%nrold)
!--------------------------------------------------------------------------
! Set the GRIB parameter specifiers
!  FORCING() ARRAY: \\
!  1. T 2m    Temperature interpolated to 2 metres [$K$] \\
!  2. q 2m    Instantaneous specific humidity interpolated to 2 metres[$kg/kg$] \\
!  3. radswg  Downward shortwave flux at the ground [$W/m^2$] \\
!  4. lwgdwn  Downward longwave radiation at the ground [$W/m^2$] \\
!  5. u 10m   Instantaneous zonal wind interpolated to 10 metres [$m/s$] \\
!  6. v 10m   Instantaneous meridional wind interpolated to 10 metres[$m/s$] \\
!  7. ps      Instantaneous Surface Pressure [$Pa$] \\
!  8. preacc  Total precipitation [$mm/s$] \\
!  9. gfrac   Greenness fraction (0-1)
! 10. albedo  Surface albedo (0-1)
! 11. zlvl    Height of atmospheric forcing [m]
! 12. z0      Surface roughness [m]
! 13. ch      Surface exchange coefficient [m/s]
! 14. t1      Skin Temperature [K]
! 15. sneqv   SWE [m]
! 16. snowh   Snow Depth [m]
!--------------------------------------------------------------------------
!J  if (get_nstep() .eq. 0) then
!J     pds5 = (/011,051,204,205,033,034,001,059,214,084,144,144, 011,011, 065/) !parameter
!J     pds7 = (/002,002,000,000,010,010,000,000,000,000,010,2760,010,2760,000/) !htlev2
!J     nforce = gdasdrv%nmif
!J  else
     pds5 = (/ 011,051,204,205,033,034,001,059,087,084,007,083,208,011,065,066 /) !parameter
     pds6 = (/ 109,109,001,001,109,109,001,001,001,001,109,001,001,001,001,001 /) !level
     pds7 = (/ 001,001,000,000,001,001,000,000,000,000,001,000,000,000,000,000 /) !height
     pds16= (/ 010,010,010,010,010,010,010,003,010,003,010,010,010,010,010,010 /) !3-ave; 10-fcst
     nforce = gdasdrv%nmif
!J  endif

  ferror = 1  
!--------------------------------------------------------------------------
! if there's a problem then ferror is set to zero
!--------------------------------------------------------------------------
  iv = 0
  errorflag = 0 
  endloop = 0
  allocate(lb(gdasdrv%ncold*gdasdrv%nrold))
  allocate( f(gdasdrv%ncold*gdasdrv%nrold))

  mi = ngdas
  mo = ld%d%lnc*ld%d%lnr
  allocate(rlat(mo))
  allocate(rlon(mo))
  allocate(  lo(mo))
  allocate(   g(mo))

  allocate(varfield(ld%d%lnc,ld%d%lnr))

  do
     if ( endloop == 1 ) exit
     iv = iv+1
     fname = name
     inquire (file=fname, exist=file_exists)
     if (file_exists) then      
!--------------------------------------------------------------------------
! Set up to open file and retrieve specified field 
!--------------------------------------------------------------------------
        lugb = iv +try
        j = 0
        jpds = -1
        jpds(5) = pds5(iv)
        jpds(6) = pds6(iv)
        jpds(7) = pds7(iv)
        jpds(16)= pds16(iv)

        call baopenr(lugb,fname,iret)
  print*,"J---retgdas()---baopenr()---iret= ", iret

        if(iret==0) then 
           call getgb(lugb,lubi,ngdas,j,jpds,jgds,kf,k,kpds,kgds,lb,f,gbret)
        else 
           gbret = 99
        endif
  print*,"J---retgdas()---getgb()----gbret= ", gbret
        call baclose(lugb,jret)
     else
        ferror = 0
        deallocate(f)
        deallocate(lb)
        deallocate(rlat)
        deallocate(rlon)
        deallocate(lo)
        deallocate(g)
        deallocate(varfield)
        return
     endif
!--------------------------------------------------------------------------
! If field successfully retrieved, interplate to LIS domain
!--------------------------------------------------------------------------
! JESSE 20050225 USE IPOLATES
!--------------------------------------------------------------------------

 if (gbret==0) then

  ip = 0
  ipopt = 0

  write(*,'(3I7)') iv, mi, mo

  write(*,'(11I7)') kpds(1:16)
  kgdsi(1:25) = kgds(1:25)
  write(*,'(11I7)') kgdsi(1:11)

  kgdso = 0
  do i = 1, 10
     if( i.EQ.4 .OR. i.EQ.7 .OR. i.EQ.8 .OR. i.EQ.9 ) then
       kgdso(i) = int(lis%d%gridDesc(i)*1000)
     else
       kgdso(i) = int(lis%d%gridDesc(i))
     endif
  enddo
  kgdso(11) = 0
  kgdso(20) = 255
  write(*,'(11I7)') kgdso(1:11)

!------------------------------------------------------
! IPOLATE FORCING TO RUN DOMAIN
!------------------------------------------------------
  iret = 0
!  call ipolates (ip,ipopt,kgdsi,kgdso,mi,mo, &
!       km,ibi,lb,f,no,rlat,rlon,ibo,lo,g,iret)
  g = f

  if(iret .NE. 0) then
        print*, "IPOLATES ERROR!! PROGRAM STOP!!"
        call exit(iret)
  end if
!------------------------------------------------------
! END JESSE 20050225 IPOLATES
!------------------------------------------------------

     count = 0
     do r = ld%d%lnr, 1, -1
        do c = 1, ld%d%lnc
           varfield(c,r) = g(c+count)
        end do
        count = count + ld%d%lnc
     end do
 else
  errorflag = 1
 endif !endif (gbret==0)

     if ( errorflag == 1 ) then 
        endloop = 1
        ferror = 0
     else
        count = 0
        do r = 1, ld%d%lnr
           do c = 1, ld%d%lnc
             if (gindex(c,r).ne. -1) then 
              if ( order == 1 ) then 
                 glbdata1(iv,gindex(c,r)) = varfield(c,r)
              else
                 glbdata2(iv,gindex(c,r)) = varfield(c,r)
              end if
             endif
          end do
          count = count +ld%d%lnc
       end do
    end if

     if ( errorflag == 1 ) then
        print *, 'Could not find correct forcing parameter in file ',name
     end if 
     if ( iv == nforce ) endloop = 1
  end do
  deallocate(lb)
  deallocate(f)
  deallocate(rlat)
  deallocate(rlon)
  deallocate(lo)
  deallocate(g)
  deallocate(varfield)

  return
!EOC
end subroutine retgdas


!BOP
! !ROUTINE: interp_gdas
!
! !DESCRIPTION:
!   This subroutine interpolates a given GDAS field 
!   to the LIS domain.  Special treatment for some
!   initialization fields.
!   Code based on old ungribgdas.f
!
! !INTERFACE:
subroutine interp_gdas(kpds,kgds,ngdas,f,lb,lis_gds,nc,nr, &
     varfield)
! !USES:
!  use def_ipMod, only : w110,w120,w210,w220,n110,n120,n210,n220,& 
!       rlat0,rlon0, &
!       w113,w123,w213,w223,n113,n123,n213,n223,rlat3,rlon3            
  use bilinear_interpMod
  use conserv_interpMod
  use gdasdomain_module, only : mi
!EOP
!w113,w123,w213,w223,n113,n123,n213,n223,rlat3, & 
!       rlon3
  implicit none

!=== Begin variable declarations
  integer :: nc, nr, ngdas, nglis
  integer :: kpds(200),kgds(200), lis_gds(200)
  integer :: ip, ipopt(20),ibi,km,iret
  integer :: no, ibo
  integer :: count,i,j,v

  real :: f(ngdas)
  real :: ism, udef
  real, dimension(nc,nr) :: varfield, geogtemp
  real, dimension(nc*nr) :: lis1d

  logical*1 :: geogmask(nc,nr)
  logical*1 :: lb(ngdas)
  logical*1 :: lo(nc*nr)

!=== End variable declarations
!BOC
!-----------------------------------------------------------------------
! Setting interpolation options (ip=0,bilinear)
! (km=1, one parameter, ibi=1,use undefined bitmap
! (needed for soil moisture and temperature only)
! Use budget bilinear (ip=3) for precip forcing fields
!-----------------------------------------------------------------------
  nglis = nc*nr
  if (kpds(5)==59 .or. kpds(5)==214) then
     ip=3
     ipopt(1)=-1
     ipopt(2)=-1
     km=1
     ibi=1          
  else
     ip=0
     do i=1,20
        ipopt(i)=0
     enddo
     km=1
     ibi=1
  endif
!-----------------------------------------------------------------------
! Initialize output bitmap. Important for soil moisture and temp.
!-----------------------------------------------------------------------
  lo = .true.
!-----------------------------------------------------------------------  
! Interpolate to LIS grid
!-----------------------------------------------------------------------  
  if (kpds(5)==59 .or. kpds(5)==214) then     
     call polates3(lis_gds,ibi,lb,f,ibo,lo,lis1d,mi, & 
          rlat3,rlon3,w113,w123,w213,w223,n113,n123,n213,n223,iret)
  else 
     call polates0 (lis_gds,ibi,lb,f,ibo,lo,lis1d,mi,&
          rlat0, rlon0,w110,w120,w210,w220,n110,n120,n210,n220,iret)
  endif
!-----------------------------------------------------------------------    
! Create 2D array for main program. Also define a "soil" mask
! due to different geography between GDAS & LDAS. For LDAS land 
! points not included in GDAS geography dataset only.
!-----------------------------------------------------------------------    
  count = 0
  do j = 1, nr
     do i = 1, nc
        varfield(i,j) = lis1d(i+count)
        geogmask(i,j) = lo(i+count)
     enddo
     count = count + nc
  enddo
!-----------------------------------------------------------------------    
! Save air tempertaure interpolated field for later use in
! initialization of soil temp where geography differs 
! between GDAS and LDAS
!-----------------------------------------------------------------------    
  if (kpds(5) .eq. 11 .and. kpds(6) .eq. 105) then
     do i = 1, nc
        do j = 1, nr
           geogtemp(i,j) = varfield(i,j)
        enddo
     enddo
  endif
!EOC
end subroutine interp_gdas