       program create_firstguess
!***********************************************************************
!   prgmmr: pondeca           org: np20                date: 2004-10-13
!
! abstract:
! generate wrf format first guess for surface analysis. input is 
! binary file contructed from first guess GRIB1 file using 'wgrib'
!
! program history log:
!   2004-10-13  pondeca
!   2007-06-11  pondeca - rewrite to handle not just the ndfd conus
!                         grid but also alaska, hawaii, guam and 
!                         puerto rico ndfd grids
!   2010-04-08  pondeca - set XICE to zero. it was set to landmask
!   2011-02-03  zhu     - add changes for pblh, gust, vis
!***********************************************************************

       implicit none

       character(60) cgrid
       include './param.incl'

       integer(4),parameter:: nt=11! # of horiz. slabs from wgrib
       integer(4),parameter:: nsig=1       
       real(4),parameter:: gravity=9.81
       real(4),parameter:: qmin=1.e-06
       real(4),parameter:: nhwrfmax=10

       integer(4) nx,ny
       integer(4) i,j,n
       integer(4) iyear,imonth,iday,ihour,iminute,isecond
       integer(4) ihwrfyear,ihwrfmonth,ihwrfday,ihwrfhour, & 
                  ihwrfminute,ihwrfsecond
       integer(4) ihdrbuf(512)
       integer(4),dimension(5):: iadate, iadate2

       real(4) ds
       real(4) rmax,rmin
       real(4) wgrid0

       integer(4),allocatable,dimension(:,:)::ifield,ivgtyp,isltyp

       real(4),allocatable,dimension(:,:)::glat,glon,field,sst,tslb, &
                                      tsk,vegfra,smois,qgrid,tdgrid, & 
                                      psfcgrid,ugrid,vgrid,landmask, & 
                                      dx,dy,mapfac,phbgrid,tgrid,gust, & 
                                      dist,vis,pblh,hgtgrid

       real(4),allocatable,dimension(:,:)::bias_qgrid,bias_psfcgrid, & 
                               bias_ugrid,bias_vgrid,bias_tgrid,bias_tdgrid, &
                               bias_gust,bias_vis,bias_pblh

       character(4) var

       real(4) pbiascor,tbiascor,qbiascor,ubiascor,vbiascor,tdbiascor,& 
               gustbiascor,visbiascor,pblhbiascor 
       logical lpbiascor,ltbiascor,lqbiascor,lubiascor,lvbiascor,ltdbiascor, &
               lgustbiascor,lvisbiascor,lpblhbiascor
       logical lbiascor
       logical mkrjlists
       logical fexist
       logical hwrfblend

       namelist/gridname/cgrid,lbiascor,lpbiascor,pbiascor,ltbiascor,tbiascor, &
                         lqbiascor,qbiascor,lubiascor,ubiascor,lvbiascor,vbiascor, & 
                         ltdbiascor,tdbiascor,lgustbiascor,gustbiascor,lvisbiascor, &
                         visbiascor,lpblhbiascor,pblhbiascor,mkrjlists,hwrfblend

       namelist/timeinfo/iyear,imonth,iday,ihour, & 
                         iminute,isecond

       namelist/hwrfinfo/ihwrfyear,ihwrfmonth,ihwrfday,ihwrfhour, & 
                         ihwrfminute,ihwrfsecond

       data cgrid/'conus'/
       data lbiascor/.false./
       data pbiascor/0.0/
       data tbiascor/0.0/
       data qbiascor/0.0/
       data ubiascor/0.0/
       data vbiascor/0.0/
       data tdbiascor/0.0/
       data gustbiascor/0.0/
       data visbiascor/0.0/
       data pblhbiascor/0.0/
       data mkrjlists/.false./
       data iyear/2005/
       data imonth/6/
       data iday/17/
       data ihour/12/
       data iminute/0/
       data isecond/0/
       data hwrfblend/.false./
       data ihwrfyear/-9999/
       data ihwrfmonth/-99/
       data ihwrfday/-99/
       data ihwrfhour/-99/
       data ihwrfminute/0/
       data ihwrfsecond/0/
!
!====================================================================
!==> find out what grid this is
!====================================================================
       open (55,file='gridname_input',form='formatted')
       read(55,gridname)
       close(55)
       print*,'in create_firstguess: cgrid is =',trim(cgrid)
       print*,'in create_firstguess: lbiascor=',lbiascor
       print*,'in create_firstguess: lpbiascor,pbiascor=',  lpbiascor,pbiascor
       print*,'in create_firstguess: ltbiascor,tbiascor=',  ltbiascor,tbiascor
       print*,'in create_firstguess: lqbiascor,qbiascor=',  lqbiascor,qbiascor
       print*,'in create_firstguess: lubiascor,ubiascor=',  lubiascor,ubiascor
       print*,'in create_firstguess: lvbiascor,vbiascor=',  lvbiascor,vbiascor
       print*,'in create_firstguess: ltdbiascor,tdbiascor=',ltdbiascor,tdbiascor
       print*,'in create_firstguess: lgustbiascor,gustbiascor=',lgustbiascor,gustbiascor
       print*,'in create_firstguess: lvisbiascor,visbiascor=',lvisbiascor,visbiascor
       print*,'in create_firstguess: lpblhbiascor,pblhbiascor=',lpblhbiascor,pblhbiascor
       print*,'in create_firstguess: mkrjlists=',mkrjlists
       print*,'in create_firstguess: hwrfblend=',hwrfblend

       if (hwrfblend) then
           inquire(file='hwrfinfo_input',exist=fexist)
           if(fexist) then
             open (55,file='hwrfinfo_input',form='formatted')
             read (55,hwrfinfo)
             close(55)
             print*,'in create_firstguess: ihwrfyear=',ihwrfyear
             print*,'in create_firstguess: ihwrfmonth=',ihwrfmonth
             print*,'in create_firstguess: ihwrfday=',ihwrfday
             print*,'in create_firstguess: ihwrfhour=',ihwrfhour
          endif
       endif
!====================================================================
!==> read in analysis time info from namelist
!====================================================================
       read(9,timeinfo)

       print*,'in create_firstguess: iyear,imonth,iday,ihour, & 
                 & iminute,isecond=',iyear,imonth,iday,ihour, & 
                   iminute,isecond 

       print*,'in create_firstguess: recomputeq =',recomputeq

       call domain_dims(cgrid,nx,ny,ds)
       print*,'in create_firstguess: nx,ny,ds=',nx,ny,ds
!
!====================================================================
!==> allocate fields 
!====================================================================
       allocate(ifield(nx,ny))
       allocate(ivgtyp(nx,ny))
       allocate(isltyp(nx,ny))
       allocate(glat(nx,ny))
       allocate(glon(nx,ny))
       allocate(field(nx,ny))
       allocate(sst(nx,ny))
       allocate(tslb(nx,ny))
       allocate(tsk(nx,ny))
       allocate(vegfra(nx,ny))
       allocate(smois(nx,ny))
       allocate(qgrid(nx,ny))
       allocate(tdgrid(nx,ny))
       allocate(psfcgrid(nx,ny))
       allocate(ugrid(nx,ny))
       allocate(vgrid(nx,ny))
       allocate(tgrid(nx,ny))
       allocate(landmask(nx,ny))
       allocate(dx(nx,ny))
       allocate(dy(nx,ny))
       allocate(mapfac(nx,ny))
       allocate(phbgrid(nx,ny))
       allocate(gust(nx,ny))
       allocate(dist(nx,ny))
       allocate(vis(nx,ny))
       allocate(pblh(nx,ny))
       allocate(hgtgrid(nx,ny))
!
!====================================================================
!=> generate file for surface analysis in " wrf first-guess format"
!====================================================================
       rewind(30)
       read(30) mapfac
       read(30) glat
       read(30) glon
       read(30) landmask
       write(*,*) maxval(mapfac),minval(mapfac)
       write(*,*) maxval(glat),minval(glat)
       write(*,*) maxval(glon),minval(glon)
       write(*,*) maxval(landmask),minval(landmask)

       dx=ds*mapfac
       dy=ds*mapfac

       ihdrbuf=0

       write(88) ihdrbuf
       write(88) iyear,imonth,iday,ihour,iminute,isecond,nx,ny,nsig
       write(88) dx,dy 
       write(88) glat  
       write(88) glon 

       rewind(20)
       do 300 n=1,nt
          print*,'in create_firstguess: read slab# ',n
          read (20) field
          print*,'rmin,rmax=',minval(field),maxval(field)

          if (n.eq.1) psfcgrid=field

          if (n.eq.2) then
            hgtgrid=field
            phbgrid=field*gravity                        ! PHB (zsfc*g)
          endif

          if (n.eq.3) then
            tgrid=field  ! T(k)                          ! TEMP (sensible)
            sst=field
            tslb=field
            tsk=field
            do j=1,ny
            do i=1,nx
               if (landmask(i,j).gt.0.5) sst(i,j)=0.
            enddo
            enddo
          endif

          if (n.eq.4) tdgrid=field
          if (n.eq.5) ugrid=field
          if (n.eq.6) vgrid=field

          if (n.eq.7) then 
             qgrid=field
             if (recomputeq) call td_to_q(tdgrid,psfcgrid,qgrid,nx,ny)
          endif

          if (n.eq.8) gust=field
          if (n.eq.9) dist=field
          if (n.eq.10) then
             print*,'vis,min,max,before capping: ', minval(vis),maxval(vis)
             do j=1,ny
             do i=1,nx
                vis(i,j)=min(20000.,field(i,j))
             enddo
             enddo
             print*,'vis,min,max,after capping: ', minval(vis),maxval(vis)
          endif
          if (n.eq.11) then 
             pblh=field
             call pblh_80km(pblh,nx,ny)
          end if
300     continue
        close(20)

        do j=1,ny
        do i=1,nx
           wgrid0=ugrid(i,j)*ugrid(i,j)+vgrid(i,j)*vgrid(i,j)
           if (wgrid0 > 0.) wgrid0=sqrt(wgrid0)
           if (gust(i,j) < wgrid0) gust(i,j)=wgrid0
        enddo
        enddo

        open (70,file='slabs2_nobiasc.dat',form='unformatted')
        write(70) psfcgrid
        write(70) hgtgrid
        write(70) tgrid
        write(70) tdgrid
        write(70) ugrid
        write(70) vgrid
        write(70) qgrid
        write(70) gust
        write(70) dist
        write(70) vis
        write(70) pblh
        close(70)

        if (lbiascor) then
           inquire(file='rtma_biascor_in.dat',exist=fexist)
           if(fexist) then
               allocate(bias_qgrid(nx,ny))
               allocate(bias_psfcgrid(nx,ny))
               allocate(bias_ugrid(nx,ny))
               allocate(bias_vgrid(nx,ny))
               allocate(bias_tgrid(nx,ny))
               allocate(bias_tdgrid(nx,ny))
               allocate(bias_gust(nx,ny))
               allocate(bias_vis(nx,ny))
               allocate(bias_pblh(nx,ny))

               bias_psfcgrid=0.
               bias_tgrid=0.
               bias_qgrid=0.
               bias_ugrid=0.
               bias_vgrid=0.
               bias_tdgrid=0.
               bias_gust=0.
               bias_vis=0.
               bias_pblh=0.

               open (60,file='rtma_biascor_in.dat',form='unformatted')
               read(60) bias_psfcgrid
               read(60) bias_tgrid
               read(60) bias_qgrid
               read(60) bias_ugrid
               read(60) bias_vgrid
               read(60) bias_tdgrid

               if (lpbiascor)   psfcgrid = psfcgrid + bias_psfcgrid
               if (ltbiascor)   tgrid = tgrid + bias_tgrid
               if (lqbiascor)   qgrid = max(qmin, (qgrid + bias_qgrid))    
               if (lubiascor)   ugrid = ugrid + bias_ugrid
               if (lvbiascor)   vgrid = vgrid + bias_vgrid
!              if (ltdbiascor) tdgrid = tdgrid + bias_tdgrid
!              if (lgustbiascor) gust = gust + bias_gust
!              if (lvisbiascor) vis = vis + bias_vis
!              if (lpblhbiascor) pblh = pblh + bias_pblh

!              tdgrid = min(tgrid,tdgrid) 

!              if (recomputeq) call td_to_q(tdgrid,psfcgrid,qgrid,nx,ny)

               close(60)
               deallocate(bias_qgrid)
               deallocate(bias_psfcgrid)
               deallocate(bias_ugrid)
               deallocate(bias_vgrid)
               deallocate(bias_tgrid)
               deallocate(bias_tdgrid)
               deallocate(bias_gust)
               deallocate(bias_vis)
               deallocate(bias_pblh)
             else
             print*,'in create_firstguess: lbias=.true., but found no input bias file.'
             print*,'assume zero input bias'
           endif
        endif

        if (hwrfblend) then 
           open (70,file='slabs2_pre_hwrf_blending.dat',form='unformatted')
           write(70) psfcgrid
           write(70) hgtgrid
           write(70) tgrid
           write(70) tdgrid
           write(70) ugrid
           write(70) vgrid
           write(70) qgrid
           write(70) gust
           write(70) dist
           write(70) vis
           write(70) pblh
           close(70)

           iadate(1)=iyear   ;  iadate2(1)=ihwrfyear
           iadate(2)=imonth  ;  iadate2(2)=ihwrfmonth
           iadate(3)=iday    ;  iadate2(3)=ihwrfday
           iadate(4)=ihour   ;  iadate2(4)=ihwrfhour
           iadate(5)=iminute ;  iadate2(5)=ihwrfminute

           
           do n=1,nhwrfmax
              call blendwind(cgrid,glat,glon,iadate,iadate2,ugrid,vgrid,gust,nx,ny,n)
           enddo
        endif
 
        field=0.
        ifield=0

        write(88) psfcgrid       !  psfc0
        write(88) phbgrid        !  PHB (zsfc*g)
        write(88) tgrid          !  T(k)  ! TEMP (sensible)
        write(88) qgrid          !  Q(k)
        write(88) ugrid          !  U(K)
        write(88) vgrid          !  V(K)
        write(88) landmask       !  LANDMASK  (0=water and >0.5 for land) 
        write(88) field          !  XICE
        write(88) sst            !  SST
        write(88) ifield         !  IVGTYP
        write(88) ifield         !  ISLTYP
        write(88) field          !  VEGFRA
        write(88) field          !  SNOW
        write(88) ugrid          !  U10
        write(88) vgrid          !  V10
        write(88) field          !  SMOIS
        write(88) tslb           !  TSLB
        write(88) tsk            !  TSK
        write(88) gust           !  GUST
        write(88) vis            !  VIS
        write(88) pblh           !  PBLH

        open(70,file='rtma_terrain.dat',form='unformatted')
           write(70) phbgrid/9.8
        close(70)
        open(70,file='rtma_slmask.dat',form='unformatted')
           write(70) LANDMASK
        close(70)

        open (70,file='slabs2.dat',form='unformatted')
        write(70) psfcgrid
        write(70) hgtgrid
        write(70) tgrid
        write(70) tdgrid
        write(70) ugrid
        write(70) vgrid
        write(70) qgrid
        write(70) gust
        write(70) dist
        write(70) vis
        write(70) pblh
        close(70)


        deallocate(ifield)
        deallocate(ivgtyp)
        deallocate(isltyp)
        deallocate(glat)
        deallocate(glon)
        deallocate(field)
        deallocate(sst)
        deallocate(tslb)
        deallocate(tsk)
        deallocate(vegfra)
        deallocate(smois)
        deallocate(qgrid)
        deallocate(tdgrid)
        deallocate(psfcgrid)
        deallocate(ugrid)
        deallocate(vgrid)
        deallocate(landmask)
        deallocate(dx)
        deallocate(dy)
        deallocate(mapfac)
        deallocate(phbgrid)
        deallocate(tgrid)
        deallocate(gust)
        deallocate(dist)
        deallocate(vis)
        deallocate(pblh)
        deallocate(hgtgrid)

        if (mkrjlists) then
          fexist=.false.
          inquire(file='bigrjlist.txt',exist=fexist)
          if (fexist)      call separate_rjlists
          if (.not.fexist) call join_rjlists
          var='t'    ;   call process_allrjlists(var)
          var='q'    ;   call process_allrjlists(var)
          var='p'    ;   call process_allrjlists(var)
          var='w'    ;   call process_allrjlists(var)
          var='mass' ;   call process_allrjlists(var)
        endif
      end 
!=======================================================================
!=======================================================================
!***********************************************************************
          subroutine td_to_q(td,p,q,nx,ny)
!
!    prgmmr: pondeca           org: np20                date: 2004-10-13
!
!  abstract: given dewpt in K and pressure in Pa, compute
!            specific humidity in kg/kg. this is done by inverting:
!            
!                qv=q(i,j)/(1.-q(i,j))
!                e=p(i,j)/100.*qv/(eps+qv)
!                eln=alog(e)
!                td(i,j) = (243.5*eln-440.8)/(19.48-eln)+273.15
!
! program history log:
!   2004-10-13  pondeca
!***********************************************************************
          implicit none
 
          real(4),parameter::alpha=243.5
          real(4),parameter::beta=440.8
          real(4),parameter::gamma=19.48
          real(4),parameter::eta=273.15
          real(4),parameter::eps=0.62197  !=Rd/Rv

          integer(4) nx,ny,i,j 
          real(4),dimension(nx,ny),intent(in)::td,p
          real(4),dimension(nx,ny),intent(inout)::q
          real(4) loge,e,qv

          do j=1,ny
          do i=1,nx
             loge=(gamma*(td(i,j)-eta)+beta)/(alpha+td(i,j)-eta) 
             e=exp(loge)
             qv=100.*e*eps/(p(i,j)-100.*e)
             q(i,j)=qv/(1.+qv)
          enddo
          enddo

          return
          end
!=======================================================================
!=======================================================================
         subroutine pblh_80km(grd,nx,ny)
         implicit none
         integer(4) nx,ny,i,j,n
         real(4),dimension(nx,ny),intent(inout):: grd
         real(4),dimension(0:nx+1,0:ny+1):: grd2
         real(4) corn,cent,side,temp,c2,c3,c4,rterm

!        Set weights for nine point smoother
         corn=0.3
         cent=1.0
         side=0.5
         c4=4.0
         c3=3.0
         c2=2.0
         rterm = 1.0/(cent + c4*side + c4*corn)

!        30 passes of 9pt smoother
         do n=1,30

            do i=1,nx
               do j=1,ny
                  grd2(i,j)=grd(i,j)
               end do
            end do
            grd2(:,0)=grd(:,1)
            grd2(:,ny+1)=grd(:,ny)
            grd2(0,:)=grd(nx,:)
            grd2(nx+1,:)=grd(1,:)

            do i=1,nx
               do j=1,ny
                  temp = cent*grd2(i,j) + side*(grd2(i+1,j) + &
                         grd2(i-1,j) + grd2(i,j+1) + grd2(i,j-1)) + &
                         corn*(grd2(i+1,j+1) + grd2(i+1,j-1) + grd2(i-1,j-1) + &
                         grd2(i-1,j+1))
                  grd(i,j) = temp*rterm         
               end do
            end do

         end do

         return
         end subroutine pblh_80km
