       subroutine cnv_to_grib2_mpi(cgrid,is,ie,js,je,nx,ny, & 
                  lbiascor,igust,ivis,ipblh,idist,mype,npe)

       use mpi
       implicit none

!Declare passed variables
       character(60),intent(in):: cgrid
       integer(4),intent(in):: is,ie,js,je
       integer(4),intent(in):: nx,ny
       integer(4),intent(in):: mype,npe
       integer(4),intent(in):: igust,ivis,ipblh,idist
       logical,intent(in):: lbiascor

!Declare local parameters
       integer(4),parameter:: nvarmax=30  !maximum number of rtma 
                                          !flds that can be processed
       real(4),parameter::qmin=1.e-06

!Declare local variables
       integer(4),parameter:: nvarbckg0=9
       integer(4),parameter:: nvaranl0=9
       integer(4),parameter:: nvarerr0=9

       integer(4) nvarbckg,nvaranl,nvarerr
       integer(4) nvar2d

       character(60),allocatable::varname(:)
       character(60),allocatable::typeproc(:)

       character(60) varnamebckg(nvarmax)
       character(60) varnameanl(nvarmax)
       character(60) varnameerr(nvarmax)

       character(60) varnamebckg0(nvarbckg0)
       character(60) varnameanl0(nvaranl0)
       character(60) varnameerr0(nvarerr0)

       character(60) filename, filename2, filename3

       integer(4) ioan1,n,n0,ncase,i,j,nn,k,j0

       integer(4) itime0(6),itime(6),nlon,nlat,nsig
       integer(4) itimem1(6),itimep1(6)
 
       integer(4) iyear01,imonth01,iday01,ihour01, &
                 iminute01,isecond01,ifcsthr

       integer(4) iyearm1,imonthm1,idaym1,ihourm1, &
                 iminutem1,isecondm1,ifcsthrm1

       integer(4) iyearp1,imonthp1,idayp1,ihourp1, &
                 iminutep1,isecondp1,ifcsthrp1

       integer(4) istat

       integer(4), allocatable,dimension(:,:)::   ifield
       real(4),    allocatable,dimension(:,:)::   field,field2
       real(4),    allocatable,dimension(:,:)::   glon
       real(4),    allocatable,dimension(:,:,:):: gndfd,gndfd7,gndfd8,gndfd9
       real(4),    allocatable,dimension(:,:)::   psfc1,psfc2,fis0, &
                   fis1,fis2,t1,t2,q1,q2,u1,u2,v1,v2,td1,td2,td0, & 
                   gust1,gust2,dist1,dist2,vis1,vis2,pblh1,pblh2, & 
                   wspd,wdir,ue,ve

       logical guessfcst
       logical guessfcstm1
       logical guessfcstp1

       logical fexist
       logical usefgat

       integer(4) isign !significance of reference time
                  !!!(0 for analysis and anlerr, and 1 for forecast)

       integer(4) ityped !type of processed data
                  !!!(0 for analysis and anlerr,  and 1 for forecast)
 
       real(4) w0
       real(4) amin1,amin2
       real(4) amax1,amax2

       integer(4) mx,my
       integer(4) iopt,iopthgt,npass,npasshgt
       real(4) smcf1,smcf2
       character(60) cgridhalfres
       logical lambconform
       logical polarstereo
       logical lhalfres
       logical one2one

       integer(4) ix,jy
       character(60) cgridnew
       integer(4),parameter::nx_cohres=2145
       integer(4),parameter::ny_cohres=1377
       integer(4),parameter::nx_nwrfc=709
       integer(4),parameter::ny_nwrfc=795

       integer(4) ierror
       integer(4) is2,ie2,js2,je2
       integer(4) is3,ie3,js3,je3

       namelist/ptimeinfo/iyear01,imonth01,iday01,ihour01, &
        iminute01,isecond01,ifcsthr,guessfcst, & ! timeinfo for reference time
        usefgat                                  ! ifcsthr is the forecast hour 
                                                 ! guessfcst=.true. if guess is
                                                 ! a forecast field

       namelist/ptimeinfom1/iyearm1,imonthm1,idaym1,ihourm1, &
        iminutem1,isecondm1,ifcsthrm1,guessfcstm1  

       namelist/ptimeinfop1/iyearp1,imonthp1,idayp1,ihourp1, &
        iminutep1,isecondp1,ifcsthrp1,guessfcstp1  

       namelist/hgridoptions/lhalfres,iopt,iopthgt,smcf1,smcf2,one2one,npass,npasshgt
!
!===================================================================================================
!==> use data statements to set default values for various variables
!===================================================================================================
!
       data lhalfres/.true./
       data iopt/1/
       data iopthgt/0/
       data smcf1/0.50/
       data smcf2/0.52/
       data one2one/.false./
       data npass/1/
       data npasshgt/0/

       data iyear01/2005/
       data imonth01/6/
       data iday01/17/
       data ihour01/12/
       data iminute01/0/
       data isecond01/0/
       data ifcsthr/1/
       data guessfcst/.true./
       data usefgat/.false./

       data iyearm1/2005/
       data imonthm1/6/
       data idaym1/17/
       data ihourm1/11/
       data iminutem1/0/
       data isecondm1/0/
       data ifcsthrm1/0/
       data guessfcstm1/.false./

       data iyearp1/2005/
       data imonthp1/6/
       data idayp1/17/
       data ihourp1/11/
       data iminutep1/0/
       data isecondp1/0/
       data ifcsthrp1/2/
       data guessfcstp1/.true./

       data varnamebckg0/ &
       'hgt', &
       'psfc', &
       't2m', &
       'td2m', &
       'u10m', &
       'v10m', &
       'q2m', &
       'wdir10m', &
       'wspd10m'/

       data varnameanl0/ &
       'hgt', &
       'psfc', &
       't2m', &
       'td2m', &
       'u10m', &
       'v10m', &
       'q2m', &
       'wdir10m', &
       'wspd10m'/

       data varnameerr0/ &
       'hgt', &
       'anlerr-psfc', &
       'anlerr-t2m', &
       'anlerr-td2m', &
       'anlerr-u10m', &
       'anlerr-v10m', &
       'anlerr-q2m', &
       'anlerr-wdir10m', &
       'anlerr-wspd10m'/
!
!===================================================================================================
!==> load fixed field names
!===================================================================================================
!
       varnamebckg(1:nvarbckg0)  =  varnamebckg0(1:nvarbckg0)  ; varnamebckg(nvarbckg0+1:nvarmax)='   '
       varnameanl(1:nvaranl0  )  =  varnameanl0(1:nvaranl0  )  ; varnameanl(nvaranl0+1  :nvarmax)='   '
       varnameerr(1:nvarerr0  )  =  varnameerr0(1:nvarerr0  )  ; varnameerr(nvarerr0+1  :nvarmax)='   '
!
!===================================================================================================
!==> more fields may be available. find out and add to varnamebckg, varnameanl & varnameerr if applicable
!===================================================================================================
!
       if (mype==0) print*,'in cnv_to_grib2: igust,idist,ivis,ipblh=',igust,idist,ivis,ipblh

       n=0 
       if (igust > 0) then
          n=n+1
          varnamebckg(nvarbckg0+n)='gust'
          varnameanl(nvaranl0+n)='gust'
          varnameerr(nvarerr0+n)='anlerr-gust'
       endif
       if (idist > 0) then
          n=n+1
          varnamebckg(nvarbckg0+n)='dist'
          varnameanl(nvaranl0+n)='dist'
          varnameerr(nvarerr0+n)='anlerr-dist'
       endif
       if (ivis > 0) then
          n=n+1
          varnamebckg(nvarbckg0+n)='vis'
          varnameanl(nvaranl0+n)='vis'
          varnameerr(nvarerr0+n)='anlerr-vis'
       endif
       if (ipblh > 0) then
          n=n+1
          varnamebckg(nvarbckg0+n)='pblh'
          varnameanl(nvaranl0+n)='pblh'
          varnameerr(nvarerr0+n)='anlerr-pblh'
       endif

       nvarbckg=nvarbckg0+n
       nvaranl=nvaranl0+n
       nvarerr=nvarerr0+n

       if (mype==0) print*,'in cnv_to_grib2: after accounting for new variables: nvarbckg,nvaranl,nvarerr=',& 
                                                                                  nvarbckg,nvaranl,nvarerr
!
!===================================================================================================
!==>  print out some useful information
!===================================================================================================
!
       read(9,ptimeinfo)   !fetch usefgat from here
       close(9)

       if (mype==0) then
          print*,'in cnv_to_grib2: cgrid=',trim(cgrid)
          print*,'in cnv_to_grib2: lbiascor=',lbiascor
          print*,'in cnv_to_grib2: usefgat=',usefgat
          print*,'in cnv_to_grib2: nx,ny=',nx,ny
       endif
!
!===================================================================================================
!==> read in namelist of smoothing options for half-resolution grids
!===================================================================================================
!
       inquire(file='hgridoptions_input',exist=fexist)
       if(fexist) then
         open (55,file='hgridoptions_input',form='formatted')
         read (55,hgridoptions)
         close(55)
       endif

       if (mype==0) then
          print*,'in cnv_to_grib2: lhalfres=',lhalfres
          print*,'in cnv_to_grib2: iopt=',iopt
          print*,'in cnv_to_grib2: iopthgt=',iopthgt
          print*,'in cnv_to_grib2: smcf1,smcf2=',smcf1,smcf2
          print*,'in cnv_to_grib2: one2one=',one2one
          print*,'in cnv_to_grib2: npass=',npass
          print*,'in cnv_to_grib2: npasshgt=',npasshgt
       endif
!
!===================================================================================================
!==> allocate local fields
!===================================================================================================
!
       allocate(ifield(nx,ny))
       allocate(field(nx,ny))
       allocate(field2(nx,ny))

       allocate( glon  (is:ie,js:je) )
       allocate( psfc1 (is:ie,js:je) )
       allocate( psfc2 (is:ie,js:je) )
       allocate( fis0  (is:ie,js:je) )
       allocate( fis1  (is:ie,js:je) )
       allocate( fis2  (is:ie,js:je) )
       allocate( t1    (is:ie,js:je) )
       allocate( t2    (is:ie,js:je) )
       allocate( q1    (is:ie,js:je) )
       allocate( q2    (is:ie,js:je) )
       allocate( u1    (is:ie,js:je) )
       allocate( u2    (is:ie,js:je) )
       allocate( ue    (is:ie,js:je) )
       allocate( v1    (is:ie,js:je) ) 
       allocate( v2    (is:ie,js:je) )
       allocate( ve    (is:ie,js:je) )
       allocate( td1   (is:ie,js:je) )
       allocate( td2   (is:ie,js:je) )
       allocate( td0   (is:ie,js:je) )
       allocate( wspd  (is:ie,js:je) )
       allocate( wdir  (is:ie,js:je) )

       allocate( gust1 (is:ie,js:je) )
       allocate( gust2 (is:ie,js:je) )
       allocate( dist1 (is:ie,js:je) )
       allocate( dist2 (is:ie,js:je) )
       allocate( vis1  (is:ie,js:je) )
       allocate( vis2  (is:ie,js:je) )
       allocate( pblh1 (is:ie,js:je) )
       allocate( pblh2 (is:ie,js:je) )
!
!===================================================================================================
!==> read in the guess and analysis fields
!===================================================================================================
!
!==>guess
       ioan1=10
       open (ioan1,file='sigges',form='unformatted')

       read(ioan1)
       read(ioan1) itime0,nlon,nlat,nsig
       read(ioan1) field,field2  !dx,by
       read(ioan1) field  !!glat
       read(ioan1) field  !!glon

       glon (is:ie,js:je) = field(is:ie,js:je) 

       if (mype==0) then
          print*,'in cnv_to_grib2, first/itime0=',itime0
          print*,'in cnv_to_grib2, first/nlon,nlat,nsig=', &
                                         nlon,nlat,nsig
       endif

       read(ioan1) field ; psfc1 (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; fis1  (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; t1    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; q1    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; u1    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; v1    (is:ie,js:je) = field(is:ie,js:je)

       do n=1,12                     !must jump a total of 12 records to get to gust1
         if (n==4 .or. n==5) then    !(see subroutine convert_binary_2d)
            read(ioan1) ifield
           else
            read(ioan1) field
         endif
       enddo

       read(ioan1) field ; gust1 (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; vis1  (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; pblh1 (is:ie,js:je) = field(is:ie,js:je)

       dist1=0.              !use this for now
       close(ioan1)

!==>analysis
       open (ioan1,file='siganl',form='unformatted')

       read(ioan1)
       read(ioan1) itime0,nlon,nlat,nsig
       read(ioan1) field,field2  !dx,dy
       read(ioan1) field  !!glat
       read(ioan1) field  !!glon

       if (mype==0) then
          print*,'in cnv_to_grib2, second/itime0=',itime0
          print*,'in cnv_to_grib2, second/nlon,nlat,nsig=', &
                                          nlon,nlat,nsig
       endif

       read(ioan1) field ; psfc2 (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; fis2  (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; t2    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; q2    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; u2    (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; v2    (is:ie,js:je) = field(is:ie,js:je)

       do n=1,12
          read(ioan1) field       !must jump 14 records to get to gust2
       enddo                      !(see subroutine convert_binary_2d)

       read(ioan1) field ; gust2 (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; vis2  (is:ie,js:je) = field(is:ie,js:je)
       read(ioan1) field ; pblh2 (is:ie,js:je) = field(is:ie,js:je)

       dist2=0.              !use this for now

       do j=js,je
       do i=is,ie
          w0=u2(i,j)*u2(i,j)+v2(i,j)*v2(i,j)
          if (w0 > 0.) w0=sqrt(w0)
          if (gust2(i,j) < w0) gust2(i,j)=w0
       enddo
       enddo

       do j=js,je
       do i=is,ie
          q2(i,j)=max(qmin,q2(i,j))
       enddo
       enddo
       close(ioan1)
!
!===================================================================================================
!==> recompute the analysis dew point to be the guess dew point td0 plus the dew point increment,
!    whereby the increment is the difference between the calculated analysis dew point td2 (via
!    subroutine get_dewpt) and the calculated guess dew point td1. A one-pass smoothing is applied 
!    to that difference field
!===================================================================================================
!
       open (11,file='slabs2.dat',form='unformatted')
       do n=1,4
         read(11) field
       enddo
       close(11)

       td0(is:ie,js:je)=field(is:ie,js:je)

       !-----------------------------------------------------
       call get_dewpt(psfc1,q1,t1,td1,is,ie,js,je)

       amin1=minval(td1(is:ie,js:je))
       amax1=maxval(td1(is:ie,js:je))

       call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror)
       call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror)

       if (mype==0) print*,'in cnv_to_grib2: td1 guess/min,max=', & 
                                                    amin2,amax2
       !-----------------------------------------------------
       call get_dewpt(psfc2,q2,t2,td2,is,ie,js,je)

       !dew point increment
       do j=js,je
       do i=is,ie
          td2(i,j) = max ( -20., td2(i,j)-td1(i,j) )
       enddo
       enddo

       amin1=minval(td2(is:ie,js:je))
       amax1=maxval(td2(is:ie,js:je))

       call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror)
       call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror)

       if (mype==0) print*,'in cnv_to_grib2: td2 inc/min,max=', & 
                                                   amin2,amax2

       !-----------------------------------------------------
       field(:,:)=0.

       do j=js,je
       do i=is,ie
          field(i,j)=td2(i,j)
       enddo
       enddo

       call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror)

       call smther_one(field2,1,nx,1,ny,2)

       do j=js,je
       do i=is,ie
          td2(i,j)=field2(i,j)
       enddo
       enddo
       deallocate(field2)

       do j=js,je
       do i=is,ie
          td2(i,j) = min ( t2(i,j), td2(i,j)+td0(i,j) )   !td1
       enddo
       enddo

       amin1=minval(td2(is:ie,js:je))
       amax1=maxval(td2(is:ie,js:je))

       call mpi_allreduce(amin1,amin2,1,mpi_real4,mpi_min,mpi_comm_world,ierror)
       call mpi_allreduce(amax1,amax2,1,mpi_real4,mpi_max,mpi_comm_world,ierror)

       if (mype==0) print*,'in cnv_to_grib2: td2 anl/min,max=', & 
                                                   amin2,amax2
!
!===================================================================================================
!==> populate gndfd with background fields
!   (the wind components are relative to the grid).
!===================================================================================================
  do 100 ncase=1,3 ! bckg , analysis , anlerr

       if (ncase.eq.1) then ; n0=nvarbckg0 ;  nvar2d=nvarbckg ; endif
       if (ncase.eq.2) then ; n0=nvaranl0  ;  nvar2d=nvaranl  ; endif
       if (ncase.eq.3) then ; n0=nvarerr0  ;  nvar2d=nvarerr  ; endif

       if (mype==0) print*,'in cnv_to_grib2: ncase,n0,nvar2d=',ncase,n0,nvar2d

       allocate(gndfd(is:ie,js:je,nvar2d))
       allocate(typeproc(nvar2d))
       allocate(varname(nvar2d))

       ioan1=10
       open (ioan1,file='slabs2.dat',form='unformatted')

       do n=1,2!
          read(ioan1) field
       enddo
       fis0(is:ie,js:je)=field(is:ie,js:je)
       close(ioan1)

       if (ncase.eq.1) then
               call wspdwdir(cgrid,glon,u1,v1,ue,ve,wspd,wdir,is,ie,js,je,mype,npe)

               do j=js,je
               do i=is,ie
                  gndfd(i,j,1)=fis0(i,j)
                  gndfd(i,j,2)=psfc1(i,j)
                  gndfd(i,j,3)=t1(i,j)
                  gndfd(i,j,4)=td0(i,j)        !td1(i,j)
                  gndfd(i,j,5)=u1(i,j)
                  gndfd(i,j,6)=v1(i,j)
                  gndfd(i,j,7)=q1(i,j)
                  gndfd(i,j,8)=wdir(i,j)
                  gndfd(i,j,9)=wspd(i,j)
               enddo
               enddo

              n=n0
              if (igust > 0) then
                 n=n+1
                 gndfd(is:ie,js:je,n)=gust1(is:ie,js:je)
              endif
              if (idist > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=dist1(is:ie,js:je)
              endif
              if (ivis > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=vis1(is:ie,js:je)
              endif
              if (ipblh > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=pblh1(is:ie,js:je)
              endif

              if (lbiascor) then
                 allocate(gndfd9(is:ie,js:je,nvar2d))

                 open (70,file='slabs2_nobiasc.dat',form='unformatted')

                 read(70) field ; gndfd9(is:ie,js:je,2) = field(is:ie,js:je)    !psfc
                 read(70) field ; gndfd9(is:ie,js:je,1) = field(is:ie,js:je)    !fis
                 read(70) field ; gndfd9(is:ie,js:je,3) = field(is:ie,js:je)    !t
                 read(70) field ; gndfd9(is:ie,js:je,4) = field(is:ie,js:je)    !td
                 read(70) field ; gndfd9(is:ie,js:je,5) = field(is:ie,js:je)    !u
                 read(70) field ; gndfd9(is:ie,js:je,6) = field(is:ie,js:je)    !v
                 read(70) field ; gndfd9(is:ie,js:je,7) = field(is:ie,js:je)    !q

                 call wspdwdir(cgrid,glon,gndfd9(is:ie,js:je,5),gndfd9(is:ie,js:je,6), & 
                               ue,ve,wspd,wdir,is,ie,js,je,mype,npe)

                 gndfd9(is:ie,js:je,8) = wdir(is:ie,js:je)    !wdir
                 gndfd9(is:ie,js:je,9) = wspd(is:ie,js:je)    !wspd

                 n=n0
                 if (igust > 0) then
                    n=n+1 
                    read(70) field ; gndfd9(is:ie,js:je,n) = field(is:ie,js:je)   !gust
                 endif
                 if (idist > 0) then
                    n=n+1
                    read(70) field ; gndfd9(is:ie,js:je,n) = field(is:ie,js:je)   !dist
                 endif
                 if (ivis > 0) then
                    n=n+1
                    read(70) field ; gndfd9(is:ie,js:je,n) = field(is:ie,js:je)   !vis
                 endif
                 if (ipblh > 0) then
                    n=n+1
                    read(70) field ; gndfd9(is:ie,js:je,n) = field(is:ie,js:je)   !pblh
                 endif
                 close(70)
              endif

              if (usefgat) then
                 allocate(gndfd7(is:ie,js:je,nvar2d))
                 allocate(gndfd8(is:ie,js:je,nvar2d))

                 open (70,file='slabsm1.dat',form='unformatted')
                 open (71,file='slabsp1.dat',form='unformatted')

                 read(70) field ; gndfd7(is:ie,js:je,2) = field(is:ie,js:je)    !psfc
                 read(70) field ; gndfd7(is:ie,js:je,1) = field(is:ie,js:je)    !fis
                 read(70) field ; gndfd7(is:ie,js:je,3) = field(is:ie,js:je)    !t
                 read(70) field ; gndfd7(is:ie,js:je,4) = field(is:ie,js:je)    !td
                 read(70) field ; gndfd7(is:ie,js:je,5) = field(is:ie,js:je)    !u
                 read(70) field ; gndfd7(is:ie,js:je,6) = field(is:ie,js:je)    !v
                 read(70) field ; gndfd7(is:ie,js:je,7) = field(is:ie,js:je)    !q

                 call wspdwdir(cgrid,glon,gndfd7(is:ie,js:je,5),gndfd7(is:ie,js:je,6), & 
                               ue,ve,wspd,wdir,is,ie,js,je,mype,npe)

                 gndfd7(is:ie,js:je,8) = wdir(is:ie,js:je)    !wdir
                 gndfd7(is:ie,js:je,9) = wspd(is:ie,js:je)    !wspd

                 read(71) field ; gndfd8(is:ie,js:je,2) = field(is:ie,js:je)    !psfc
                 read(71) field ; gndfd8(is:ie,js:je,1) = field(is:ie,js:je)    !fis
                 read(71) field ; gndfd8(is:ie,js:je,3) = field(is:ie,js:je)    !t
                 read(71) field ; gndfd8(is:ie,js:je,4) = field(is:ie,js:je)    !td
                 read(71) field ; gndfd8(is:ie,js:je,5) = field(is:ie,js:je)    !u
                 read(71) field ; gndfd8(is:ie,js:je,6) = field(is:ie,js:je)    !v
                 read(71) field ; gndfd8(is:ie,js:je,7) = field(is:ie,js:je)    !q

                 call wspdwdir(cgrid,glon,gndfd8(is:ie,js:je,5),gndfd8(is:ie,js:je,6), & 
                               ue,ve,wspd,wdir,is,ie,js,je,mype,npe)

                 gndfd8(is:ie,js:je,8) = wdir(is:ie,js:je)    !wdir
                 gndfd8(is:ie,js:je,9) = wspd(is:ie,js:je)    !wspd

                 n=n0
                 if (igust > 0) then
                    n=n+1 
                    read(70) field ; gndfd7(is:ie,js:je,n) = field(is:ie,js:je)   !gust
                    read(71) field ; gndfd8(is:ie,js:je,n) = field(is:ie,js:je)   !gust
                 endif
                 if (idist > 0) then
                    n=n+1
                    read(70) field ; gndfd7(is:ie,js:je,n) = field(is:ie,js:je)   !dist
                    read(71) field ; gndfd8(is:ie,js:je,n) = field(is:ie,js:je)   !dist
                 endif
                 if (ivis > 0) then
                    n=n+1
                    read(70) field ; gndfd7(is:ie,js:je,n) = field(is:ie,js:je)   !vis
                    read(71) field ; gndfd8(is:ie,js:je,n) = field(is:ie,js:je)   !vis
                 endif
                 if (ipblh > 0) then
                    n=n+1
                    read(70) field ; gndfd7(is:ie,js:je,n) = field(is:ie,js:je)   !pblh
                    read(71) field ; gndfd8(is:ie,js:je,n) = field(is:ie,js:je)   !pblh
                 endif
                 close(70)
                 close(71)
              endif

        else if (ncase.eq.2) then

              call wspdwdir(cgrid,glon,u2,v2,ue,ve,wspd,wdir,is,ie,js,je,mype,npe)

              do j=js,je
              do i=is,ie
                 gndfd(i,j,1)=fis0(i,j)
                 gndfd(i,j,2)=psfc2(i,j)
                 gndfd(i,j,3)=t2(i,j)
                 gndfd(i,j,4)=td2(i,j)
                 gndfd(i,j,5)=u2(i,j)
                 gndfd(i,j,6)=v2(i,j)
                 gndfd(i,j,7)=q2(i,j)
                 gndfd(i,j,8)=wdir(i,j)
                 gndfd(i,j,9)=wspd(i,j)
              enddo
              enddo

                        !********************************************************************
                        !ADD THIS TO REPRODUCE WSPD AND WDIR FROM OPS RTMA BIT BY BIT:
                        !RATHER THAN ROTATING THE FULL FIELDS (U2,V2) AND CALLING SPEEDDIR,
                        !THE OPS RTMA ROTATES THE GUESS AND THE INCREMENT SEPARATELY (A LINEAR 
                        !OPERATION!), ADDS THE ROTATED FIELDS TOGETHER, AND THEN CALLS SPEEDDIR.
                        !THE DIFFERENCES IN THE WDIR AND WSPD COMPUTED FROM THE TWO METHODS 
                        !ARE SIMPLY DUE TO THE USE OF FINITE PRECISION NUMERICS, 
                        !AND THEY ARE INSIGNIFICANT. THIS BLOCK WAS ONLY INCLUDED TO MAKE SURE
                        !THAT NO BUG CREPT INTO THE CODE WHILE ADDING MPI TO IT. REMOVE THE
                        !BLOCK EVENTUALLY.

                         u2(is:ie,js:je)=u2(is:ie,js:je)-u1(is:ie,js:je)
                         v2(is:ie,js:je)=v2(is:ie,js:je)-v1(is:ie,js:je)

                         call wspdwdir(cgrid,glon,u1,v1,ue,ve,wspd,wdir,is,ie,js,je,mype,npe)
                         u1(is:ie,js:je)=ue(is:ie,js:je)
                         v1(is:ie,js:je)=ve(is:ie,js:je)

                         call wspdwdir(cgrid,glon,u2,v2,ue,ve,wspd,wdir,is,ie,js,je,mype,npe)
 
                         ue(is:ie,js:je)=u1(is:ie,js:je)+ue(is:ie,js:je)
                         ve(is:ie,js:je)=v1(is:ie,js:je)+ve(is:ie,js:je)

                         do j=js,je
                         do i=is,ie
                            call speeddir(ue(i,j),ve(i,j),wdir(i,j),wspd(i,j))
                         enddo
                         enddo

                         do j=js,je
                         do i=is,ie
                            gndfd(i,j,8)=wdir(i,j)
                            gndfd(i,j,9)=wspd(i,j)
                         enddo
                         enddo
                        !********************************************************************

              n=n0
              if (igust > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=gust2(is:ie,js:je)
              endif
              if (idist > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=dist2(is:ie,js:je)
              endif
              if (ivis > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=vis2(is:ie,js:je)
              endif
              if (ipblh > 0) then
                  n=n+1
                  gndfd(is:ie,js:je,n)=pblh2(is:ie,js:je)
              endif

        else if (ncase.eq.3) then

              inquire(file='bckg_z.dat',exist=fexist)

              if (fexist) then
                 allocate(field2(ny,nx)) !transposed field

                 open(94,file='bckg_z.dat',form='unformatted')
                 read(94) field2
                 field=transpose(field2)
                 gndfd(is:ie,js:je,1)=field(is:ie,js:je)
                 close(94)

                 deallocate(field2)
               else
                 gndfd(is:ie,js:je,1)=fis0(is:ie,js:je)
              endif

              open (53,file='errfield.dat',form='unformatted')

              do n=2,9     !perr,terr,tderr,uerr,verr,qerr,wdirerr2,wspderr
                read(53) field
                gndfd(is:ie,js:je,n)=field(is:ie,js:je)
              enddo

              read(53) field           !skip uerr2
              read(53) field           !skip verr2
              read(53) field           !skip wdirerr
         
              do n=n0+1,nvar2d           !gusterr,disterr,viserr,pblherr
                read(53,iostat=istat) field
                if (istat/=0) exit
                gndfd(is:ie,js:je,n)=field(is:ie,js:je)
              enddo
              close(53)
       endif
!===================================================================================================
!==> prepare to write grib2 file
!===================================================================================================
       if (ncase.eq.1) then 
          read(9,ptimeinfo)
          if (mype==0) & 
                print*,'in cnv_to_grib2,iyear01,imonth01,iday01,ihour01, &
     &          iminute01,isecond01,ifcsthr,guessfcst=',iyear01,imonth01, &
                iday01,ihour01,iminute01,isecond01,ifcsthr,guessfcst
          close(9)
          if (guessfcst) then
              itime(1)=iyear01
              itime(2)=imonth01
              itime(3)=iday01
              itime(4)=ihour01
              itime(5)=iminute01
              itime(6)=isecond01
              isign=1
              ityped=1
              typeproc(1:nvar2d)='forecast'
            else
              itime(1:6)=itime0(1:6)
              isign=0
              ityped=0
              typeproc(1:nvar2d)='analysis'
          endif
          varname(1:nvar2d)=varnamebckg(1:nvar2d)
       endif
       
       if (ncase.eq.2) then 
          itime(1:6)=itime0(1:6)
          isign=0
          ityped=0
          typeproc(1:nvar2d)='analysis'
          varname(1:nvar2d)=varnameanl(1:nvar2d)
       endif

       if (ncase.eq.3) then 
          itime(1:6)=itime0(1:6)
          isign=0
          ityped=0
          typeproc(1:nvar2d)='analysis error'
          varname(1:nvar2d)=varnameerr(1:nvar2d)
       endif

       if (ncase.eq.1) filename='bckg.grib2'
       if (ncase.eq.2) filename='anl.grib2'
       if (ncase.eq.3) filename='anlerr.grib2'

!===================================================================================================
!==> write grib2 file
!===================================================================================================
       ! first  create binary output for further processing if desired 
       if (mype==0) then
          open (41,file=trim(filename)//'.bin',form='unformatted')
          write(41) nvar2d
          write(41) varname,typeproc
          write(41) isign,ityped,ifcsthr
          write(41) itime
       endif

       allocate(field2(nx,ny))
       do n=1,nvar2d
          field(:,:)=0.
          do j=js,je
          do i=is,ie
             field(i,j)=gndfd(i,j,n)
          enddo
          enddo

          call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror)
          if (mype==0) write(41) field2
       enddo
       close(41)
       deallocate(field2)

       call write_ndfd_grib2(gndfd,is,ie,js,je,nx,ny,nvar2d,& 
            varname,typeproc,isign, & 
            itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
            ityped,ifcsthr,cgrid,filename,mype,npe)

       if (ncase.eq.1) then 
         if (lbiascor) then
            filename='bckg.grib2_nobiasc'
             call write_ndfd_grib2(gndfd9,is,ie,js,je,nx,ny,nvar2d,& 
                  varname,typeproc,isign, & 
                  itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
                  ityped,ifcsthr,cgrid,filename,mype,npe)
         endif

         if (usefgat) then
            read(19,ptimeinfom1)
            close(19)

            if (mype==0) & 
            print*,'in cnv_to_grib2,iyearm1,imonthm1,idaym1,ihourm1, &
     &      iminutem1,isecondm1,ifcsthrm1,guessfcstm1=',iyearm1,imonthm1, &
            idaym1,ihourm1,iminutem1,isecondm1,ifcsthrm1,guessfcstm1

            itime(1)=iyearm1
            itime(2)=imonthm1
            itime(3)=idaym1
            itime(4)=ihourm1
            itime(5)=iminutem1
            itime(6)=isecondm1
            if (guessfcstm1) then 
               typeproc(1:nvar2d)='forecast'
               isign=1
               ityped=1
              else
               typeproc(1:nvar2d)='analysis'
               isign=0
               ityped=0
            endif
            filename='bckg.grib2_m1'
            call write_ndfd_grib2(gndfd7,is,ie,js,je,nx,ny,nvar2d,& 
                 varname,typeproc,isign, & 
                 itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
                 ityped,ifcsthrm1,cgrid,filename,mype,npe)

            read(29,ptimeinfop1)
            close(29)

            if (mype==0) & 
            print*,'in cnv_to_grib2,iyearp1,imonthp1,idayp1,ihourp1, &
     &      iminutep1,isecondp1,ifcsthrp1,guessfcstp1=',iyearp1,imonthp1, &
            idayp1,ihourp1,iminutep1,isecondp1,ifcsthrp1,guessfcstp1

            itime(1)=iyearp1
            itime(2)=imonthp1
            itime(3)=idayp1
            itime(4)=ihourp1
            itime(5)=iminutep1
            itime(6)=isecondp1
            typeproc(1:nvar2d)='forecast'
            isign=1
            ityped=1
            filename='bckg.grib2_p1'
            call write_ndfd_grib2(gndfd8,is,ie,js,je,nx,ny,nvar2d,& 
                 varname,typeproc,isign, & 
                 itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
                 ityped,ifcsthrp1,cgrid,filename,mype,npe)
         endif
       endif

       deallocate(gndfd)
       deallocate(typeproc)
       deallocate(varname)
100 continue
!===================================================================================================
!==> deallocate fields
!===================================================================================================
       deallocate(field)
       deallocate(ifield)
       deallocate(glon)
       deallocate(psfc1)
       deallocate(psfc2)
       deallocate(fis0)
       deallocate(fis1)
       deallocate(fis2)
       deallocate(t1)
       deallocate(t2)
       deallocate(q1)
       deallocate(q2)
       deallocate(u1)
       deallocate(ue)
       deallocate(u2)
       deallocate(v1) 
       deallocate(v2)
       deallocate(ve)
       deallocate(td1)
       deallocate(td2)
       deallocate(td0)
       deallocate(gust1)
       deallocate(gust2)
       deallocate(dist1)
       deallocate(dist2)
       deallocate(vis1)
       deallocate(vis2)
       deallocate(pblh1)
       deallocate(pblh2)
       deallocate(wspd)
       deallocate(wdir)

       if (lbiascor) deallocate(gndfd9)

       if (usefgat) then
          deallocate(gndfd7)
          deallocate(gndfd8)
       endif
       call mpi_barrier(mpi_comm_world,ierror)
!
!===================================================================================================
!==> if cgrid is 'cohresext', map output to the true conus2p5 ndfd grid
!    and the nwrfc grid
!===================================================================================================

       if (trim(cgrid)=='cohresext') then

          do ncase=1,3 ! bckg , analysis , anlerr
   
             if (ncase.eq.1) filename='bckg.grib2'
             if (ncase.eq.2) filename='anl.grib2'
             if (ncase.eq.3) filename='anlerr.grib2'

             open (41,file=trim(filename)//'.bin',form='unformatted')

             do nn=1,2          !map to cohres and nwrfc
                if (nn==1) then
                   ix=nx_cohres
                   jy=ny_cohres
                   cgridnew='cohres'
                   filename2=trim(filename)//'_map2cohres'
                   cgridhalfres='conus'
                 endif

                if (nn==2) then
                   ix=nx_nwrfc
                   jy=ny_nwrfc
                   cgridnew='nwrfc'
                   filename2=trim(filename)//'_map2nwrfc'
                   cgridhalfres='nwrfc5p0'
                endif

                call horiz_domain_partition(ix,jy,mype,npe,is2,ie2,js2,je2)

                print*,'in cnv_to_grib2:filename2,is2,ie2,js2,je2=', & 
                                  trim(filename2),is2,ie2,js2,je2

                allocate(gndfd8(is2:ie2,js2:je2,nvar2d))
                allocate(field(nx,ny))

                rewind(41)
                read(41) nvar2d

                allocate(typeproc(nvar2d))
                allocate(varname(nvar2d))

                read(41) varname,typeproc
                read(41) isign,ityped,ifcsthr
                read(41) itime

                do k=1,nvar2d
                   read(41) field
                   if (nn==1) then 
                       do j=js2,je2
                       do i=is2,ie2
                          gndfd8(i,j,k)=field(i,j)
                       enddo
                       enddo
                   endif

                   if (nn==2) then 
                       j0=ny-jy
                       do j=js2,je2
                       do i=is2,ie2
                          gndfd8(i,j,k)=field(i,j0+j) !note that js2=1 and je2=ny_nwrfc
                       enddo
                       enddo
                   endif
                enddo

                if (nn==2) then 
                   if (mype==0) print*,'in cnv_to_grib2: for nwrfc grid, j0=',j0
                endif

                call write_ndfd_grib2(gndfd8,is2,ie2,js2,je2,ix,jy,nvar2d,& 
                     varname,typeproc,isign, & 
                     itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
                     ityped,ifcsthr,cgridnew,filename2,mype,npe)

                if (nn==1 .and. lhalfres) then !note: no need for half grid for nwrfc
                   mx=(ix-1+0.0001)/2+1 
                   my=(jy-1+00.001)/2+1 

                   if (mype==0) print*,'in cnv_to_grib2: & 
                      &dimensions for half-resolution grid ', & 
                      trim(cgridhalfres),' are mx,my=',mx,my

                   call horiz_domain_partition(mx,my,mype,npe,is3,ie3,js3,je3)

                   print*,'in cnv_to_grib2: cgridhalfres,is3,ie3,js3,je3=', & 
                                        trim(cgridhalfres),is3,ie3,js3,je3

                   allocate(gndfd(is3:ie3,js3:je3,nvar2d))

                   call resolution_halfing(gndfd8,ix,jy,nvar2d,is2,ie2,js2,je2, & 
                              gndfd,mx,my,is3,ie3,js3,je3,varname,iopt,iopthgt, & 
                              smcf1,smcf2,one2one,npass,npasshgt,mype,npe)

                   filename3=trim(filename2)//'_halfres'

                   call write_ndfd_grib2(gndfd,is3,ie3,js3,je3,mx,my,nvar2d,&
                        varname,typeproc,isign, &
                        itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), &
                        ityped,ifcsthr,cgridhalfres,filename3,mype,npe)

                    deallocate(gndfd)
                endif

                 deallocate(typeproc)
                 deallocate(varname)
                 deallocate(gndfd8)
                 deallocate(field)
             enddo
             close(41)
          enddo
       endif

!===================================================================================================
!==> generate half resolution analysis when cgrid is cohres or akhres.
!    note that for cgrid=cohresext this is taken care of in the previous block
!===================================================================================================

       if (lhalfres) then
          if (trim(cgrid)=='cohres' .or. trim(cgrid)=='akhres') then

             if (trim(cgrid)=='cohres') cgridhalfres='conus'
             if (trim(cgrid)=='akhres') cgridhalfres='alaska'

             do ncase=1,3 ! bckg or analysis , anlerr
   
                if (ncase.eq.1) filename='bckg.grib2'
                if (ncase.eq.2) filename='anl.grib2'
                if (ncase.eq.3) filename='anlerr.grib2'

                open (41,file=trim(filename)//'.bin',form='unformatted')

                read(41) nvar2d

                allocate(typeproc(nvar2d))
                allocate(varname(nvar2d))

                read(41) varname,typeproc
                read(41) isign,ityped,ifcsthr
                read(41) itime

                allocate(gndfd9(is:ie,js:je,nvar2d))
                allocate(field(nx,ny))

                do k=1,nvar2d
                   read(41) field
                   do j=js,je
                   do i=is,ie
                      gndfd9(i,j,k)=field(i,j)
                   enddo
                   enddo
                enddo
                close(41)
       
                mx=(nx-1+0.0001)/2+1 
                my=(ny-1+00.001)/2+1 

                if (mype==0) print*,'in cnv_to_grib2: & 
                   &dimensions for half-resolution grid ', & 
                   trim(cgridhalfres),' are mx,my=',mx,my

                call horiz_domain_partition(mx,my,mype,npe,is3,ie3,js3,je3)

                print*,'in cnv_to_grib2: cgridhalfres,is3,ie3,js3,je3=', & 
                                     trim(cgridhalfres),is3,ie3,js3,je3

                allocate(gndfd(is3:ie3,js3:je3,nvar2d))

                call resolution_halfing(gndfd9,nx,ny,nvar2d,is,ie,js,je, & 
                     gndfd,mx,my,is3,ie3,js3,je3,varname,iopt,iopthgt, & 
                     smcf1,smcf2,one2one,npass,npasshgt,mype,npe)

                filename=trim(filename)//'_halfres'

                call write_ndfd_grib2(gndfd,is3,ie3,js3,je3,mx,my,nvar2d,& 
                     varname,typeproc,isign, & 
                     itime(1),itime(2),itime(3),itime(4),itime(5),itime(6), & 
                     ityped,ifcsthr,cgridhalfres,filename,mype,npe)

                deallocate(field)
                deallocate(gndfd9)
                deallocate(typeproc)
                deallocate(varname)
                deallocate(gndfd)
             enddo
          endif
       endif

       return
       end
!===================================================================================================
!===================================================================================================
!***************************************************************************************************
      subroutine resolution_halfing(gin,nx,ny,nk,is,ie,js,je, & 
                                   gout,mx,my,m1,m2,n1,n2,varname,iopt,iopthgt, & 
                                   smcf1,smcf2,one2one,npass,npasshgt,mype,npe)
!***************************************************************************************************
!abstract: half the resolution of the input hres grid
! iopt=0 ==> select every other grid point of hres grid
! iopt=1 ==> 0.5x+0.125*(values at four nearest grid points of hres grid)
! iopt=2 ==> 0.5x+0.125*(values at four vertice points of 
!                        enclosing square of hres grid)
! note: x is the hres field value at the targetted grid point
! note: for wdir10m, code always used iopt=0
! 
! iopthgt: for terrain height, override iopt with iopthgt
!
! smcf1,smcf2 : coefficients used with smoother-desmoother
! one2one     : if .true. perform 1-2-1 smoothing rather than smoothing-desmoothing
! npass       : number of soothing passes
! npasshgt    : for terrain height, override npass with npasshgt
!***************************************************************************************************
      use mpi

      implicit none

!Declare passed variables
      integer(4),intent(in)::nx,ny,nk,is,ie,js,je
      integer(4),intent(in)::mx,my,m1,m2,n1,n2
      integer(4),intent(in)::iopt,iopthgt,npass,npasshgt
      integer(4),intent(in)::mype,npe
      real(4),intent(in)::gin(is:ie,js:je,nk)
      real(4),intent(in)::smcf1,smcf2
      character(60),intent(in)::varname(nk)
      logical,intent(in)::one2one
      real(4),intent(out)::gout(m1:m2,n1:n2,nk)

!Declare local variables
      integer(4) iopt0,npass0
      integer(4) i,j,i2,j2,k
      integer(4) im,ip,jm,jp
      integer(4) ierror

      real(4),allocatable,dimension(:,:):: hfine,hfine2
      real(4),allocatable,dimension(:,:):: hcoarse,hcoarse2
!---------------------------------------------------------------------------------------------------
      allocate(hfine(nx,ny))
      allocate(hfine2(nx,ny))
      allocate(hcoarse(mx,my))
      allocate(hcoarse2(mx,my))

      do k=1,nk
         if (mype==0) print*,'in resolution_halfing: k, varname=',k,varname(k)

         hfine2(:,:)=0.

         do j=js,je
         do i=is,ie
            hfine2(i,j)=gin(i,j,k)
         enddo
         enddo

         call mpi_allreduce(hfine2,hfine,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror)

         hcoarse(:,:)=0.

         if (trim(varname(k))=='hgt') then 
             iopt0=iopthgt

           else if(trim(varname(k))=='wdir10m') then 
             iopt0=0  !won't try to average wdir

           else
             iopt0=iopt
         endif

         if (iopt0.ne.0 .and. iopt0.ne.1  .and. iopt0.ne.2) then 
            if (mype==0) then
               print*,'in resolution_halfing: warning ... invalid iopt0=',iopt0,' value'
               print*,'in resolution_halfing: reset to iopt0=1'
            endif
            iopt0=1 
         endif

         do j2=n1,n2
            j=2*j2-1
            jm=max(j-1,1)
            jp=min(j+1,ny)

            do i2=m1,m2
               i=2*i2-1
               im=max(i-1,1)
               ip=min(i+1,nx)

               if (iopt0==0) then 
                  hcoarse(i2,j2)=hfine(i,j) 
                elseif (iopt0==1) then
                  hcoarse(i2,j2)=0.5*hfine(i,j) + & 
                                0.125*(hfine(i,jm)+hfine(i,jp)+hfine(im,j)+hfine(ip,j))
                elseif (iopt0==2) then
                  hcoarse(i2,j2)=0.5*hfine(i,j) + & 
                                0.125*(hfine(im,jm)+hfine(ip,jm)+hfine(im,jp)+hfine(ip,jp))
               endif
            enddo
         enddo

         if (trim(varname(k))=='hgt') then
             npass0=npasshgt
           else if(trim(varname(k))=='wdir10m') then 
             npass0=0
           else
             npass0=npass
         endif

         call mpi_allreduce(hcoarse,hcoarse2,mx*my,mpi_real4,mpi_sum,mpi_comm_world,ierror)
         call smoother_desmoother(hcoarse2,1,mx,1,my,smcf1,smcf2,one2one,npass0)

         do j=n1,n2
         do i=m1,m2
            gout(i,j,k)=hcoarse2(i,j)
         enddo
         enddo
      enddo

      deallocate(hfine)
      deallocate(hfine2)
      deallocate(hcoarse)
      deallocate(hcoarse2)

      return
      end
!***************************************************************
      subroutine get_earth_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, &
                                 olon,dg2rad)
!
!     Abstract: Rotates wind relative to grid on Lambert Conformal
!               or polar streo projection to get winds with respect
!               to the true north.
!     Adapted from Stan Benjamin's code snippet

      implicit none
      real(4) un,vn,ut,vt,ROTCON_P,LON_XX_P,olon, &
              dg2rad,angle2,sinx2,cosx2

!**  ROTCON_P          R  WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO
!**                         AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL
!**  LAT_TAN_P         R  LATITUDE AT LAMBERT CONFORMAL PROJECTION
!**                         IS TRUE (DEG)
!**  LON_XX_P          R  MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG)

           angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad
           sinx2 = sin(angle2)
           cosx2 = cos(angle2)
           un = cosx2*ut+sinx2*vt
           vn =-sinx2*ut+cosx2*vt

           return
           end
!***************************************************************
      subroutine get_grid_winds(un,vn,ut,vt,ROTCON_P,LON_XX_P, &
                                 olon,dg2rad)
!
!     Abstract: Rotates wind relative to true north to 
!               get winds with respect to grid on Lambert Conformal 
!               or polar streo projection.
!     Adapted from Stan Benjamin's code snippet which performs the
!     inverse operation

      implicit none
      real(4) un,vn,ut,vt,ROTCON_P,LON_XX_P,olon, &
              dg2rad,angle2,sinx2,cosx2

!**  ROTCON_P          R  WIND ROTATION CONSTANT, = 1 FOR POLAR STEREO
!**                         AND SIN(LAT_TAN_P) FOR LAMBERT CONFORMAL
!**  LAT_TAN_P         R  LATITUDE AT LAMBERT CONFORMAL PROJECTION
!**                         IS TRUE (DEG)
!**  LON_XX_P          R  MERIDIAN ALIGNED WITH CARTESIAN X-AXIS(DEG)

           angle2 = rotcon_p*(olon-lon_xx_p)*dg2rad
           sinx2 = sin(angle2)
           cosx2 = cos(angle2)
           ut = cosx2*un-sinx2*vn
           vt = sinx2*un+cosx2*vn

           return
           end
!*******************************************************************************
!*******************************************************************************
       subroutine wspdwdir(cgrid,glon,u,v,ue,ve,w,wd,i1,i2,j1,j2,mype,npe)

       implicit none

! Declare passed variables
       character(60),intent(in):: cgrid
       integer(4),intent(in)::i1,i2,j1,j2
       integer(4),intent(in)::mype,npe

       real(4),dimension(i1:i2,j1:j2),intent(in) :: glon,u,v
       real(4),dimension(i1:i2,j1:j2),intent(out):: ue,ve,w,wd

! Declare local variables
       real(8) alat18,elon18,dx8,elonv8,alatan8
       real(4) dg2rad,elonv4,alatan4,xn
       real(4) ylon
    
       integer(4) i,j

       logical lambconform
       logical polarstereo


       call proj_info(cgrid,dx8,alat18,elon18,elonv8,alatan8)

       if (mype==0) then
          print*,'in wspdwdir: cgrid=',trim(cgrid)
          print*,'in wspdwdir: alat18=',alat18
          print*,'in wspdwdir: elon18=',elon18
          print*,'in wspdwdir: dx8=',dx8
          print*,'in wspdwdir: elonv8=',elonv8
          print*,'in wspdwdir: alatan8=',alatan8
       endif

       lambconform=trim(cgrid)=='conus'.or.trim(cgrid)=='cohres'.or. & 
                   trim(cgrid)=='hrrr'.or.trim(cgrid)=='cohresext'

       polarstereo=trim(cgrid)=='alaska'.or.trim(cgrid)=='akhres'.or. &
                   trim(cgrid)=='juneau'

       dg2rad=atan(1.)/45.
       elonv4=elonv8
       alatan4=alatan8

       xn=1.
       if (lambconform) xn=sin(alatan4*dg2rad)

       if (mype==0) print*,'in wspdwdir:  lambconform,polarstereo,xn=', &
                                          lambconform,polarstereo,xn


       if (lambconform .or. polarstereo) then
           do j=j1,j2
           do i=i1,i2
              ylon=glon(i,j)/dg2rad                                          !CHECK SIGN OF GLON FOR CONUS + ALASKA 
              call get_earth_winds(ue(i,j),ve(i,j),u(i,j),v(i,j),xn,elonv4-360.,ylon,dg2rad)
           enddo
           enddo
         else
           ue(i1:i2,j1:j2)=u(i1:i2,j1:j2)
           ve(i1:i2,j1:j2)=v(i1:i2,j1:j2)
       endif

       do j=j1,j2
       do i=i1,i2
          call speeddir(ue(i,j),ve(i,j),wd(i,j),w(i,j))
       enddo
       enddo

       return
       end
!*******************************************************************************
!*******************************************************************************
       subroutine speeddir(ue,ve,wdir,wspd)
! 
       implicit none

!Declare passed variables
       real(4),intent(in):: ue,ve
       real(4),intent(out):: wdir,wspd

!Declare local variables
       real(4) dg2rad,aearth

       dg2rad=atan(1.)/45.

       wspd=ue*ue+ve*ve
       if (wspd.ne.0.0) wspd=sqrt(wspd)

       if (wspd .eq. 0.0) then
           wdir=0.!   calm !99999.0
         else
           if (ve.eq.0.0) then
              if (ue.gt.0.0) wdir = 270.
              if (ue.lt.0.0) wdir =  90.
           else
              aearth = atan(ue/ve)/dg2rad
              if (ue.le.0.0 .and. ve.le.0.0 ) wdir = aearth
              if (ue.le.0.0 .and. ve.ge.0.0 ) wdir = aearth + 180.0
              if (ue.ge.0.0 .and. ve.ge.0.0 ) wdir = aearth + 180.0
              if (ue.ge.0.0 .and. ve.le.0.0 ) wdir = aearth + 360.0
           endif
       endif

       return
       end
!*******************************************************************************
!*******************************************************************************
       subroutine get_dewpt(p,q,t,td,is,ie,js,je)
!
!==> input variables
!    q-specific humidity ; p-pressure in Pa 
!    t-temperature in K
                                                                                                                  
!==> output variable
!    td-dewpoint in K
                                                                                                                  
       implicit none

       integer(4),intent(in):: is,ie,js,je
       real(4),dimension(is:ie,js:je),intent(in)::q,p,t
       real(4),dimension(is:ie,js:je),intent(out)::td

       real(4),parameter::eps=0.62197   !=Rd/Rv

       integer(4) i,j
       real(4) e, qv, eln
                                                                                                                  
       do j=js,je
       do i=is,ie
          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
          td(i,j) = min(t(i,j),td(i,j))
       enddo
       enddo
        
       return
       end
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       subroutine write_ndfd_grib2(g,i1,i2,j1,j2,nx,ny,nflds,varname,typeproc, & 
            isign,y0,m0,d0,h0,m1,s0,ityped,ifcsthr,clocation,filename,mype,npe)

!******************************************************************
!  prgmmr: pondeca           org: np20         date: 2006-03-03   *
!                                                                 *
!  abstract:                                                      *
!  use steve gilbert's gribcreate, addgrid, addfield, and gribend *
!  subroutines to write NDFD files in Grib2 format                * 
!                                                                 *
!  program history log:                                           *
!    2006-03-03  pondeca                                          *
!                                                                 *
!  input argument list:                                           *
!                                                                 *              
! 1. g(i1:i2,j1:j2,nflds): floating point array that contains all *
!    the nflds fields on the mpi subdomains of dimensions         *
!    (i1:i2,j1:j2). The global dimensions are (nx,ny)             *
!                                                                 *
! 2. varname(nflds): character array that conatains the names of  *
!    all the nflds fields from g(nx,ny,nflds)                     *
!                                                                 *
! 3. typeproc(nflds): character array that contains process types *
!    (analysis, forecast, or analysis error)                      *
!                                                                 *
! 4. y0,m0,d0,h0,m1,s0: reference time (year, month, day, hour,   *
!    minutes and seconds)                                         *
!                                                                 *
! 5. isign: significance of reference time. Values of interest to *
!    NDFD are 0 for analysis and 1 for start of forecast.         *
!                                                                 *
! 6. ityped: type of processed data. Values of interest to        *
!    NDFD are 0 for analysis products and 1 for forecast products *
!                                                                 *
! 7. ifcsthr: forecast hour from reference time                   *
!                                                                 *
! 8. clocation: string variable. denotes grid name. supported     *
!    grids are "conus", "alaska" , "hawaii", "guam",  "prico".    *
!    "cohres", "akhres" , "hrrr" , "juneau" , "cohresext" ,       *
!    "conusext" , and "nwrfc"                                     *
!                                                                 *
!  output argument list:                                          *
! 1. filename: Name of output Grib2 file                          *
!                                                                 *
! attributes:                                                     *
!   language: f90                                                 *
!   machine:  ibm RS/6000 SP                                      *
!******************************************************************
      use mpi

      implicit none

      integer(4),intent(in):: mype,npe
      integer(4),intent(in):: i1,i2,j1,j2
      integer(4),intent(in):: nx,ny,nflds,isign, & 
                              y0,m0,d0,h0,m1,s0,ityped,ifcsthr
 
      integer(4), parameter::max_bytes=200*130000
      integer(4), parameter::idefnum=1
      integer(4), parameter::ipdstmplen=15
      integer(4), parameter::idrstmplen=5

      integer(4) listsec0(2)
      integer(4) listsec1(13)
      integer(4) ierr,n,i,j,ij
      integer(4) lengrib,dscal
      integer(4) lunout

      integer(4) igds(5),igdstmplen
      integer(4) ideflist(idefnum)     
      integer(4) ipdstmpl(ipdstmplen)
      integer(4) idrstmpl(idrstmplen)
      integer(4) idrsnum,ibmap,numcoord,ipdsnum
      integer(4),allocatable,dimension(:):: igdstmpl

      character*60 filename
      character*60 varname(nflds)
      character*60 typeproc(nflds)
      character*1 cgrib(max_bytes)
      character*60 clocation

      logical*1 bmap(nx,ny)
 
      real(4) g(i1:i2,j1:j2,nflds)
      real(4) coordlist(1)

      real(4),allocatable,dimension(:):: fld
      real(4),allocatable,dimension(:,:):: field,field2

      integer(4) ierror

      allocate(fld(nx*ny))
      allocate(field(nx,ny))
      allocate(field2(nx,ny))

      if (mype==0) print*,'---- starting write_ndfd_grib2 -------'

!==>determine ndfd grid location and set igdstmplen dimension
      if (mype==0) print*,'grid name is=',trim(clocation)

      if (trim(clocation)=="puerto_rico") clocation="prico"

      if (trim(clocation) /= "conus"  .and. & 
          trim(clocation) /= "alaska" .and. &
          trim(clocation) /= "hawaii" .and. &
          trim(clocation) /= "guam"  .and. &
          trim(clocation) /= "prico" .and. &
          trim(clocation) /= "cohres" .and. &
          trim(clocation) /= "akhres" .and. &
          trim(clocation) /= "hrrr" .and. &
          trim(clocation) /= "dtc" .and. &
          trim(clocation) /= "juneau" .and. &
          trim(clocation) /= "cohresext" .and. &
          trim(clocation) /= "conusext" .and. &
          trim(clocation) /= "nwrfc" ) then
          if (mype==0) print*,'grid not supported. aborting ...'
          call abort
      endif

       if (trim(clocation)=="conus")  igdstmplen=22
       if (trim(clocation)=="alaska") igdstmplen=18
       if (trim(clocation)=="hawaii") igdstmplen=19
       if (trim(clocation)=="guam")   igdstmplen=19
       if (trim(clocation)=="prico" ) igdstmplen=19
       if (trim(clocation)=="cohres") igdstmplen=22
       if (trim(clocation)=="akhres") igdstmplen=18
       if (trim(clocation)=="hrrr" )  igdstmplen=22
       if (trim(clocation)=="rap" )  igdstmplen=22
       if (trim(clocation)=="juneau") igdstmplen=18
       if (trim(clocation)=="cohresext") igdstmplen=22
       if (trim(clocation)=="conusext")  igdstmplen=22
       if (trim(clocation)=="nwrfc")  igdstmplen=22

       allocate(igdstmpl(igdstmplen))

!==>open file that will store grib2 messages

     if (mype==0) then
        call getlun90(lunout,1)

        print*,'about to open',lunout,trim(filename)

        call baopenw(lunout,trim(filename),ierr)

        print*,'write_ndfd_grib2:  opened ',lunout, &
               'for grib2 data  ',trim(filename), &
               'return code is ',ierr
      endif


      do n=1,nflds

!==>initialize new GRIB2 message and pack
! GRIB2 sections 0 (Indicator Section) and 1 (Identification 
! Section)

       listsec0(1)=0 ! Discipline-GRIB Master Table Number (see Code Table 0.0)
       listsec0(2)=2 ! GRIB Edition Number (currently 2)
       
       listsec1(1)=7       ! Id of orginating centre (Common Code Table C-1)
       listsec1(2)=4 !"EMC"! Id of orginating sub-centre (local table)/Table C of ON388
       listsec1(3)=1       ! GRIB Master Tables Version Number (Code Table 1.0)
       listsec1(4)=1 !per Brent! GRIB Local Tables Version Number (Code Table 1.1)
       listsec1(5)=isign   ! Significance of Reference Time (Code Table 1.2)
       listsec1(6)=y0      ! Reference Time - Year (4 digits)
       listsec1(7)=m0      ! Reference Time - Month
       listsec1(8)=d0      ! Reference Time - Day
       listsec1(9)=h0      ! Reference Time - Hour
       listsec1(10)=m1     ! Reference Time - Minute
       listsec1(11)=s0     ! Reference Time - Second
       listsec1(12)=0      ! Production status of data (Code Table 1.3)
       listsec1(13)=ityped ! Type of processed data (Code Table 1.4)
        
       call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr)
       if (mype==0) print*,'gribcreate status=',ierr

!==> Pack up Grid Definition Section (Section 3) add to GRIB2 message.

       igds(1)=0      !Source of grid definition (see Code Table 3.0)
       igds(2)=nx*ny  !Number of grid points in the defined grid.
       igds(3)=0      !Number of octets needed for each additional grid points definition
       igds(4)=0      !Interpretation of list for optional points definition (Code Table 3.11)

       if (trim(clocation)=="conus")  igds(5)=30 !Grid Definition Template Number (Code Table 3.1) 
       if (trim(clocation)=="alaska") igds(5)=20 
       if (trim(clocation)=="hawaii") igds(5)=10
       if (trim(clocation)=="guam"  ) igds(5)=10
       if (trim(clocation)=="prico")  igds(5)=10
       if (trim(clocation)=="cohres") igds(5)=30
       if (trim(clocation)=="akhres") igds(5)=20 
       if (trim(clocation)=="hrrr")   igds(5)=30
       if (trim(clocation)=="juneau") igds(5)=20 
       if (trim(clocation)=="cohresext") igds(5)=30
       if (trim(clocation)=="conusext")  igds(5)=30
       if (trim(clocation)=="nwrfc")  igds(5)=30

       if (trim(clocation)=="conus")  call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) 
       if (trim(clocation)=="alaska") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation) 
       if (trim(clocation)=="hawaii") call apply_template_310_ndfd(igdstmpl,igdstmplen,"h") 
       if (trim(clocation)=="guam")   call apply_template_310_ndfd(igdstmpl,igdstmplen,"g") 
       if (trim(clocation)=="prico")  call apply_template_310_ndfd(igdstmpl,igdstmplen,"p") 
       if (trim(clocation)=="cohres") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) 
       if (trim(clocation)=="akhres") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation)
       if (trim(clocation)=="hrrr")   call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation)
       if (trim(clocation)=="juneau") call apply_template_320_ndfd(igdstmpl,igdstmplen,clocation)
       if (trim(clocation)=="cohresext") call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) 
       if (trim(clocation)=="conusext")  call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation) 
       if (trim(clocation)=="nwrfc")  call apply_template_330_ndfd(igdstmpl,igdstmplen,clocation)

       ideflist=0     !Used if igds(3) .ne. 0. Dummy array otherwise

       call addgrid(cgrib,max_bytes,igds,igdstmpl,igdstmplen,ideflist,idefnum,ierr)
       if (mype==0) print*,'addgrid status=',ierr

       
!==> pack up sections 4 through 7 for a given field and add them to a GRIB2 message.  
! They are Product Definition Section, Data Representation Section, Bit-Map Section 
! and Data Section, respectively.


       ipdsnum=0    !Product Definition Template Number ( see Code Table 4.0)

       if (mype==0) print*,'write_ndfd_grib2:, n,varname,typeproc=', & 
                                               n,varname(n),typeproc(n)

       call apply_template_40_ndfd(ipdstmpl,ipdstmplen,ifcsthr,& 
                                varname(n),typeproc(n))

       numcoord=0
       coordlist=0. !needed for hybrid vertical coordinate
       idrsnum=0    !Data Representation Template Number ( see Code Table 5.0 )

       call apply_template_50_ndfd(idrstmpl,idrstmplen,varname(n))
       dscal=idrstmpl(3)   !NO NEED I BELIEVE /MPondeca, May 15, 2007

       field=0.
       do j=j1,j2
       do i=i1,i2
          field(i,j)=g(i,j,n)
       enddo
       enddo

       call mpi_allreduce(field,field2,nx*ny,mpi_real4,mpi_sum,mpi_comm_world,ierror)

       ij=0 
       do j=1,ny
       do i=1,nx
          ij=ij+1
          fld(ij)=field2(i,j)
       enddo
       enddo

       ibmap=255     ! Bitmap indicator ( see Code Table 6.0 )

       call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmplen, &
                          coordlist,numcoord,idrsnum,idrstmpl, &
                          idrstmplen,fld,nx*ny,ibmap,bmap,ierr)
       if (mype==0) print*,'addfield status=',ierr

!     enddo                                                         !CORRECT WAY. UNCOMMENT! 

!==> finalize  GRIB message after all grids
! and fields have been added.  It adds the End Section ( "7777" )

      call gribend(cgrib,max_bytes,lengrib,ierr)
      if (mype==0) print*,'gribend status=',ierr
      if (mype==0) print*,'length of the final GRIB2 message in octets =',lengrib

      if (mype==0) call wryte(lunout, lengrib, cgrib)

      enddo                                                          !INCORRET WAY. DELETE THIS LINE

      if (mype==0) print*,'---- exiting write_ndfd_grib2 -------'

      deallocate(igdstmpl)
      deallocate(fld)
      deallocate(field)
      deallocate(field2)

      call mpi_barrier(mpi_comm_world,ierror)

      return
      end
!===========================================================================
!===========================================================================
       subroutine apply_template_310_ndfd(ifield3,len3,c1) 

       implicit none

       integer(4) nx,ny,lat1,lon1,lat2,lon2,lad,ds1
       integer(4) len3
       integer(4) ifield3(len3)
       character(1) c1

       if (c1=="h") then  !hawaii
           nx=321
           ny=225
!          lat1=18066780   !lat of 1st grd pt in micro-deg
!          lon1=198374755  !east-long of 1st grd pt in micro-deg 
!          lat2=23082000   !lat of last grd pt in micro-deg
!          lon2=206031000  !east-long of last grd pt in micro-deg 

           lat1=18072699   !lat of 1st grd pt in micro-deg
           lon1=198474999  !east-long of 1st grd pt in micro-deg
           lat2=23087799   !lat of last grd pt in micro-deg
           lon2=206130999  !east-long of last grd pt in micro-deg
           lad=20000000    !lat at which projection intersects earth
           ds1=2500000     !grid spacing in x and y
         elseif (c1=="g") then  !guam
           nx=193
           ny=193
           lat1=12349884
           lon1=143686538
           lat2=16794399
           lon2=148280000
           lad=20000000
           ds1=2500000
         elseif (c1=="p") then  !puerto rico
           nx=177
           ny=129
           lat1=16828685
           lon1=291804687
           lat2=19747399
           lon2=296027600
           lad=20000000
           ds1=2500000
       endif

                                                                                                        
       ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
       ifield3(2) = 0
       ifield3(3) = 0
       ifield3(4) = 0
       ifield3(5) = 0
       ifield3(6) = 0
       ifield3(7) = 0
       ifield3(8) = nx
       ifield3(9) = ny
       ifield3(10) = lat1
       ifield3(11) = lon1
       ifield3(12) = 56! 8! 56  !CHECK/MPondeca !file from NDFD site using "0"
       ifield3(13) = lad
       ifield3(14) = lat2
       ifield3(15) = lon2
       ifield3(16) = 64   !CHECK/MPondeca! Web notes have value of 80
       ifield3(17) = 0
       ifield3(18) = ds1
       ifield3(19) = ds1
 
       return
       end
!===========================================================================
!===========================================================================
       subroutine apply_template_320_ndfd(ifield3,len3,cgrid) 

       implicit none

       integer(4) len3
       integer(4) ifield3(len3)
       character*60  cgrid

       integer(4) nx,ny   !number of grid points in x and y-direction
       integer(4) ds1     !grid spacing in x and y
       integer(4) lat1    !lat of 1st grd pt in micro-deg
       integer(4) lon1    !east-long of 1st grd pt in micro-deg
       integer(4) latan1  !true lat in micro-deg
       integer(4) lonv    !y-axis || to long-circle at this long

       if (trim(cgrid)=='alaska') then
          nx=825
          ny=553
          lat1=40530101
          lon1=181429000
          latan1=60000000
          lonv=210000000
          ds1=5953125

        elseif (trim(cgrid)=='akhres') then
          nx=1649
          ny=1105
          lat1=40530101
          lon1=181429000
          latan1=60000000
          lonv=210000000
          ds1=2976563

        elseif (trim(cgrid)=='juneau') then
          nx=655
          ny=855
          lat1=51500000
          lon1=217500000
          latan1=60000000
          lonv=225000000
          ds1=1448281

       endif

       ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
       ifield3(2) = 0
       ifield3(3) = 0
       ifield3(4) = 0
       ifield3(5) = 0
       ifield3(6) = 0
       ifield3(7) = 0
       ifield3(8) = nx
       ifield3(9) = ny
       ifield3(10) = lat1
       ifield3(11) = lon1
       ifield3(12) = 8! 56 !CHECK/MPondeca !NDFD file from Web uses 0
       ifield3(13) = latan1
       ifield3(14) = lonv
       ifield3(15) = ds1
       ifield3(16) = ds1
       ifield3(17) = 0
       ifield3(18) = 64  !CHECK/MPondeca !NDFD file from Web uses 80
 
       return
       end
!===========================================================================
!===========================================================================
!===========================================================================
       subroutine apply_template_330_ndfd(ifield3,len3,cgrid) 

       implicit none

       integer(4) len3
       integer(4) ifield3(len3)
       character*60  cgrid

       integer(4) nx,ny   !number of grid points in x and y-direction
       integer(4) ds1     !grid spacing in mm
       integer(4) lat1    !lat of 1st grd pt in micro-deg
       integer(4) lon1    !east-long of 1st grd pt in micro-deg
       integer(4) latan1  !true lat in micro-deg
       integer(4) lonv    !y-axis || to long-circle at this long

       if (trim(cgrid)=='conus') then
          nx=1073
          ny=689 
          ds1=5079406
          lat1=20191999
          lon1=238445999
          latan1=25000000
          lonv=265000000

        else if (trim(cgrid)=='cohres') then
          nx=2145
          ny=1377
          ds1=2539703
          lat1=20191999
          lon1=238445999
          latan1=25000000
          lonv=265000000

        else if (trim(cgrid)=='hrrr') then  ! hrrr
          nx=1799
          ny=1059
          lat1=21138000
          lon1=237280000
          latan1=38500000
          lonv=262500000
          ds1=3000000

        else if (trim(cgrid)=='cohresext') then
          nx=2145
          ny=1597
          ds1=2539703
          lat1=20191999
          lon1=238445999
          latan1=25000000
          lonv=265000000

        else if (trim(cgrid)=='conusext') then
          nx=1073
          ny=799
          ds1=5079406
          lat1=20191999
          lon1=238445999
          latan1=25000000
          lonv=265000000

       else if (trim(cgrid)=='nwrfc') then
          nx=709
          ny=795
          ds1=2539703
          lat1=37979684
          lon1=234042704
          latan1=25000000
          lonv=265000000

       end if

       ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
       ifield3(2) = 0
       ifield3(3) = 0
       ifield3(4) = 0
       ifield3(5) = 0
       ifield3(6) = 0
       ifield3(7) = 0
       ifield3(8) = nx
       ifield3(9) = ny
       ifield3(10) = lat1
       ifield3(11) = lon1
       ifield3(12) = 8! 56 
       ifield3(13) = latan1! Or is it 35N ?
       ifield3(14) = lonv
       ifield3(15) = ds1
       ifield3(16) = ds1
       ifield3(17) = 0
       ifield3(18) = 64
       ifield3(19) = latan1
       ifield3(20) = latan1
       ifield3(21) = -90000000
       ifield3(22) = 0
 
       return
       end
!===========================================================================
       subroutine apply_template_40_ndfd(ifield4,len4,ifcsthr,var,typep) 

       implicit none

       integer(4) len4,ifcsthr
       integer(4) ifield4(len4)
   
       character*60 var
       character*60 typep
       
!==> ifield4(1):parameter category (see Code Table 4.1)
!==> ifield4(2):parameter number (see Code Table 4.2)

       if (trim(var) .eq. 'psfc' .or. & 
           trim(var) .eq. 'anlerr-psfc') then
         ifield4(1) = 3
         ifield4(2) = 0
       elseif (trim(var) .eq. 't2m' .or. & 
               trim(var) .eq. 'anlerr-t2m' ) then
         ifield4(1) = 0
         ifield4(2) = 0
       elseif (trim(var) .eq. 'td2m' .or. & 
               trim(var) .eq. 'anlerr-td2m' ) then
         ifield4(1) = 0
         ifield4(2) = 6
       elseif (trim(var) .eq. 'u10m' .or. & 
               trim(var) .eq. 'anlerr-u10m') then
         ifield4(1) = 2
         ifield4(2) = 2
       elseif (trim(var) .eq. 'v10m' .or. & 
               trim(var) .eq. 'anlerr-v10m') then
         ifield4(1) = 2
         ifield4(2) = 3
       elseif (trim(var) .eq. 'wdir10m' .or. & 
               trim(var) .eq. 'anlerr-wdir10m') then
         ifield4(1) = 2
         ifield4(2) = 0
       elseif (trim(var) .eq. 'wspd10m' .or. & 
               trim(var) .eq. 'anlerr-wspd10m') then
         ifield4(1) = 2
         ifield4(2) = 1
       elseif (trim(var) .eq. 'q2m' .or. & 
               trim(var) .eq. 'anlerr-q2m') then
         ifield4(1) = 1
         ifield4(2) = 0
       elseif (trim(var) .eq. 'hgt')then 
         ifield4(1) = 3
         ifield4(2) = 5
       elseif (trim(var) .eq. 'gust' .or. &
               trim(var) .eq. 'anlerr-gust')then
         ifield4(1) = 2
         ifield4(2) = 22
       elseif (trim(var) .eq. 'vis' .or. &
               trim(var) .eq. 'anlerr-vis')then
         ifield4(1) = 19
         ifield4(2) = 0
       elseif (trim(var) .eq. 'pblh' .or. &
               trim(var) .eq. 'anlerr-pblh')then
         ifield4(1) = 19
         ifield4(2) = 12
        else
         print*,'define param category and param number for field ',var
         print*,'sorry, ... aborting in apply_template_40_ndfd' 
         call abort
       endif 

!==> ifield4(3):type of generating process (see Code Table 4.3)
       if (trim(typep) .eq. 'analysis') then
         ifield4(3) = 0
        elseif (trim(typep) .eq. 'forecast') then
         ifield4(3) = 2
        elseif (trim(typep) .eq. 'analysis error') then
         ifield4(3) = 7
        else
         print*,'for RTMA,typep must be either analysis, forecast, &
     &           or analysis error' 
         print*,'sorry, ... aborting in apply_template_40_ndfd' 
         call abort
       endif

!==> ifield4(4):background generating process identifier 
!                         (defined by originating Center)
       ifield4(4) = 0 !hasn't been defined yet 

!==> ifield4(5):analysis or forecast generating process identifier 
!                         (defined by originating Center)
       ifield4(5) = 109

!==> ifield4(6):hours of observational data cutoff after reference time 
!==> ifield4(7):minutes of observational data cutoff after reference time 

       ifield4(6) = 0   !per steve
       ifield4(7) = 0   

!==> ifield4(8):indicator of unit of time range (see Code Table 4.4) 
       ifield4(8) = 1   

!==> ifield4(9):forecast time in units defined by ifield4(8) 
       ifield4(9) = 0
       if (trim(typep) .eq. 'forecast') ifield4(9) = ifcsthr

!==> ifield4(10):type of first fixed surface (see Code Table 4.5)
!==> ifield4(11):scale factor of first fixed surface
!==> ifield4(12):scaled value of first fixed surface

       ifield4(11) = 0 !Because not saving any precision

       if (trim(var) .eq. 'psfc' .or. trim(var) .eq. 'anlerr-psfc' .or. & 
           trim(var) .eq. 'hgt' .or.  & 
           trim(var) .eq. 'vis' .or. trim(var) .eq. 'anlerr-vis' ) then
         ifield4(10) = 1
         ifield4(12) = 0
        else
         ifield4(10) = 103

         if (trim(var) .eq. 't2m' .or. & 
             trim(var) .eq. 'anlerr-t2m') ifield4(12) = 2

         if (trim(var) .eq. 'td2m' .or. & 
             trim(var) .eq. 'anlerr-td2m') ifield4(12) = 2

         if (trim(var) .eq. 'u10m' .or. & 
             trim(var) .eq. 'anlerr-u10m') ifield4(12) = 10

         if (trim(var) .eq. 'v10m' .or. & 
             trim(var) .eq. 'anlerr-v10m') ifield4(12) = 10

         if (trim(var) .eq. 'wdir10m' .or. & 
             trim(var) .eq. 'anlerr-wdir10m') ifield4(12) = 10

         if (trim(var) .eq. 'wspd10m' .or. & 
             trim(var) .eq. 'anlerr-wspd10m') ifield4(12) = 10

         if (trim(var) .eq. 'q2m' .or. & 
             trim(var) .eq. 'anlerr-q2m') ifield4(12) = 2

         if (trim(var) .eq. 'gust' .or. &
             trim(var) .eq. 'anlerr-gust') ifield4(12) = 10

         if (trim(var) .eq. 'pblh' .or. &
             trim(var) .eq. 'anlerr-pblh') ifield4(12) = 220
       endif

!==> ifield4(13):type of second fixed surface(See Code Table 4.5)
!==> ifield4(14):scale factor of second fixed surface
!==> ifield4(15):scaled value of second fixed surface

       ifield4(13) = 255
       ifield4(14) = 0
       ifield4(15) = 0

       return
       end
!===========================================================================
       subroutine apply_template_50_ndfd(ifield5,len5,var) 

!********************************************************************************
! ifield5(1): reference value(R) (IEEE 32-bit floating-point value)             *
! ifield5(2): binary scale factor (E)                                           *
! ifield5(3): decimal scale factor (D)                                          *
! ifield5(4): number of bits used for each packed value for simple packing      *
!             or for each group reference value for complex packing or          *
!             spatial differencing                                              *
! ifield5(5): type of original field values (See Code Table 5.1)                *
!********************************************************************************
       implicit none

       integer(4) len5
       integer(4) ifield5(len5)
       character*60 var

   
       ifield5(1)=0 !Any value. Will be overwritten
       ifield5(2)=0 

       if (trim(var) .eq. 'psfc') then
         ifield5(3) = 0! -1
        elseif (trim(var) .eq. 't2m') then
         ifield5(3) = 2!  1
        elseif (trim(var) .eq. 'td2m') then
         ifield5(3) = 2!  2
        elseif (trim(var) .eq. 'u10m') then
         ifield5(3) = 2!  1
        elseif (trim(var) .eq. 'v10m') then
         ifield5(3) = 2!  1
        elseif (trim(var) .eq. 'wdir10m') then
         ifield5(3) = 1!  0
        elseif (trim(var) .eq. 'wspd10m') then
         ifield5(3) = 2!  1
        elseif (trim(var) .eq. 'q2m') then
         ifield5(3) = 6
        elseif (trim(var) .eq. 'gust') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'vis') then
         ifield5(3) = 0
        elseif (trim(var) .eq. 'pblh') then
         ifield5(3) = 0
        elseif (trim(var) .eq. 'anlerr-psfc') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-t2m') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-td2m') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-u10m') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-v10m') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-q2m') then
         ifield5(3) = 6
        elseif (trim(var) .eq. 'anlerr-wdir10m') then
         ifield5(3) = 1!  0
        elseif (trim(var) .eq. 'anlerr-wspd10m') then
         ifield5(3) = 2!  1
        elseif (trim(var) .eq. 'hgt') then
         ifield5(3) = 1
        elseif (trim(var) .eq. 'anlerr-gust') then
         ifield5(3) = 2
        elseif (trim(var) .eq. 'anlerr-vis') then
         ifield5(3) = 0
        elseif (trim(var) .eq. 'anlerr-pblh') then
         ifield5(3) = 0
        else
         print*,'please define decimal scale factor ',var
         print*,'sorry, ... aborting in apply_template_50_ndfd'
         call abort
        endif

       ifield5(4) = 0 !Must reset to 0
       ifield5(5) = 0
      
       return
       end
!===========================================================================
!===========================================================================
!$$$  SUBPROGRAM DOCUMENTATION BLOCK 
!                .      .    .                                       . 
! SUBPROGRAM:    GETLUN      GET UNIQUE LOGICAL UNIT NUMBERS
!   PRGMMR: SMITH, TRACY     ORG: FSL/PROFS  DATE: 90-06-15 
! 
! ABSTRACT: THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE
!   OR RETURNS THEM TO THE POOL FOR CLFILE.
! 
! PROGRAM HISTORY LOG: 
! FORTRAN 90 VERSION IS GETLUN90:  PONDECA,      DATE: 2006-03-08
! 
! USAGE:    CALL GETLUN(LUN,OPTN) 
!   INPUT ARGUMENT LIST: 
!     LUN      - INTEGER  LOGICAL UNIT NUMBER
!     OPTN     - INTEGER  CNCT=1, DSCT=2.
!                IF CONNECTING A FILE(CNCT) SET THE NUMBER TO
!                NEGATIVE SO IT WON'T BE USED UNTIL AFTER
!                DSCT SETS IT POSITIVE.
! 
!   OUTPUT ARGUMENT LIST:   
!     LUN      - INTEGER  LOGICAL UNIT NUMBER
! 
! REMARKS: 
! 
! ATTRIBUTES: 
!   LANGUAGE: FORTRAN-90
!   MACHINE:  NAS-9000 
!$$$ 

       SUBROUTINE GETLUN90(LUN,OPTN)
!* THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE
!* OR RETURNS THEM TO THE POOL FOR CLFILE
       IMPLICIT NONE
       INTEGER LUN,NUM(80),OPTN,CNCT,DSCT,I
       PARAMETER (CNCT=1,DSCT=2)
       SAVE NUM
! 
       DATA NUM/99,98,97,96,95,94,93,92,91,90, &
                  89,88,87,86,85,84,83,82,81,80, &
                  79,78,77,76,75,74,73,72,71,70, &
                  69,68,67,66,65,64,63,62,61,60, &
                  59,58,57,56,55,54,53,52,51,50, &
                  49,48,47,46,45,44,43,42,41,40, &
                  39,38,37,36,35,34,33,32,31,30, &
                  29,28,27,26,25,24,23,22,21,20/
!* START
       IF(OPTN.EQ.CNCT) THEN
        DO 10 I=1,80
         IF(NUM(I).GT.0) THEN
           LUN=NUM(I)
           NUM(I)=-NUM(I)
           GOTO 20
         ENDIF
10      CONTINUE
         PRINT*, 'NEED MORE THAN 80 UNIT NUMBERS'
20      CONTINUE

       ELSE IF(OPTN.EQ.DSCT) THEN
!* MAKE THE NUMBER AVAILABLE BY SETTING POSITIVE
         DO 30 I=1,80
           IF(LUN.EQ.-NUM(I)) NUM(I)=ABS(NUM(I))
30       CONTINUE
       END IF
       RETURN
       END
!===========================================================================
!=======================================================================
!=======================================================================
subroutine smther_one_same(g1,is,ie,js,je,ns)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   smther_one
! prgmmr: pondeca          org: np22                date: 2006-08-01
!
! abstract: apply 1-2-1 smoother in each direction of data slab
!
! program history log:
!   2006-08-01  pondeca
!
!   input argument list:
!    g1        - 2d array of field to be smoothed
!    is,ie     - first and last i values of data slab
!    js,je     - first and last j values of data slab
!    ns        - number of passes
!
!   output argument list:
!    g1        - smoothed 2d field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!

  implicit none

  integer(4)  is, ie, js, je
  integer(4)  i,j,l,ip,im,jp,jm
  integer(4), intent(in) :: ns

  real(4), dimension(is:ie, js:je), intent(inout) :: g1
                            !  on input: original data slab
                            !  on ouput: filtered data slab
 

  real(4), allocatable:: g2(:,:)

   allocate(g2(is:ie,js:je))
   do l=1,ns

     do j=js,je
      do i=is,ie
       ip=min(i+1,ie) ; im=max(is,i-1)
         g2(i,j)=.25*(g1(ip,j)+g1(im,j))+.5*g1(i,j)
      end do
     end do

     do i=is,ie
      do j=js,je
       jp=min(j+1,je) ; jm=max(js,j-1)
       g1(i,j)=.25*(g2(i,jp)+g2(i,jm))+.5*g2(i,j)
      end do
     end do

   end do
   deallocate(g2)

   return
end subroutine smther_one_same
!=======================================================================
!=======================================================================
subroutine smoother_desmoother(g1,is,ie,js,je,smcf1,smcf2,one2one,ns)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   smoother_desmoother
! prgmmr: pondeca          org: np22                date: 2006-08-01
!
! abstract: apply smoother-desmoother to 2-d data slab
!
! program history log:
!   2011-09-15  pondeca
!
!   input argument list:
!    g1        - 2d array of field to be smoothed
!    is,ie     - first and last i values of data slab
!    js,je     - first and last j values of data slab
!    ns        - number of passes
!
!   output argument list:
!    g1        - smoothed 2d field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  implicit none

  integer, parameter :: r_single = selected_real_kind(6)  ! single precision
  integer, parameter :: i_kind = selected_int_kind(8)     ! long  integer

  integer(i_kind)                        , intent(in   ) :: is, ie, js, je
  integer(i_kind)                        , intent(in   ) :: ns

  real(r_single), dimension(is:ie, js:je), intent(inout) :: g1
                                   !  on input: original data slab
                                   !  on ouput: filtered data slab

  real(r_single), intent(in) :: smcf1 ! smoothing coefficient
  real(r_single), intent(in) :: smcf2 ! desmoothing coefficient
  logical,intent(in) :: one2one       ! if .true., perform 1-2-1 smoothing instead

  integer(i_kind)  i,j,l,ip,im,jp,jm,istep
  real(r_single), allocatable:: g2(:,:)
  real(r_single) smcf0

  real(r_single),parameter::half=0.5_r_single
  real(r_single),parameter::quarter=0.25_r_single

  allocate(g2(is:ie,js:je))
  do l=1,ns

     if (.not.one2one) then 

        do istep=1,2
           if (istep==1) smcf0=+abs(smcf1)
           if (istep==2) smcf0=-abs(smcf2)
           do j=js,je
              do i=is,ie
                 ip=min(i+1,ie) ; im=max(is,i-1)
                 g2(i,j)=g1(i,j)+smcf0*(0.5*(g1(ip,j)+g1(im,j))-g1(i,j)) 
              end do
           end do

           do i=is,ie
              do j=js,je
                 jp=min(j+1,je) ; jm=max(js,j-1)
                 g1(i,j)=g2(i,j)+smcf0*(0.5*(g2(i,jp)+g2(i,jm))-g2(i,j))
              end do
           end do
        end do

       elseif (one2one) then 

        do j=js,je
           do i=is,ie
              ip=min(i+1,ie) ; im=max(is,i-1)
              g2(i,j)=quarter*(g1(ip,j)+g1(im,j))+half*g1(i,j)
           end do
        end do

        do i=is,ie
           do j=js,je
              jp=min(j+1,je) ; jm=max(js,j-1)
              g1(i,j)=quarter*(g2(i,jp)+g2(i,jm))+half*g2(i,j)
           end do
        end do
     endif

  end do
  deallocate(g2)

  return
end subroutine smoother_desmoother
!=======================================================================
!=======================================================================
