text:
c-----------------------------------------------------------------------
      subroutine uzawa_gmres(res,h1,h2,h2inv,intype,iter)

c     Solve the pressure equation by right-preconditioned
c     GMRES iteration.
c     intype =  0  (steady)
c     intype =  1  (explicit)
c     intype = -1  (implicit)

      include 'SIZE'
      include 'TOTAL'
      include 'GMRES'
      common  /ctolpr/ divex
      common  /cprint/ ifprint
      logical          ifprint
      real             res  (lx2*ly2*lz2*lelv)
      real             h1   (lx1,ly1,lz1,lelv)
      real             h2   (lx1,ly1,lz1,lelv)
      real             h2inv(lx1,ly1,lz1,lelv)

      common /scrmg/    wp (lx2,ly2,lz2,lelv)

      common /ctmp0/   wk1(lgmres),wk2(lgmres)
      common /cgmres1/ y(lgmres)

      real alpha, l, temp
      integer j,m
c
      logical iflag
      save    iflag
      data    iflag /.false./
      real    norm_fac
      save    norm_fac
c
      real*8 etime1,dnekclock
c
      if(.not.iflag) then
         iflag=.true.
         call uzawa_gmres_split0(ml,mu,bm2,bm2inv,nx2*ny2*nz2*nelv)
         norm_fac = 1./sqrt(volvm2)
      endif
c
      etime1 = dnekclock()
      etime_p = 0.
      divex = 0.
      iter  = 0
      m = lgmres
c
      call chktcg2(tolps,res,iconv)
      if (param(21).gt.0.and.tolps.gt.abs(param(21)))
     $   tolps = abs(param(21))
c     if (param(21).lt.0) tolps = abs(param(21))
      if (istep.eq.0) tolps = 1.e-4
      tolpss = tolps
c
      ntot2  = nx2*ny2*nz2*nelv
c
      iconv = 0
      call rzero(x,ntot2)

      do while(iconv.eq.0.and.iter.lt.100)

         if(iter.eq.0) then
                                                  !      -1
            call col3(r,ml,res,ntot2)             ! r = L  res
c           call copy(r,res,ntot2)
         else
            !update residual
            call copy(r,res,ntot2)                ! r = res
            call cdabdtp(w,x,h1,h2,h2inv,intype)  ! w = A x
            call add2s2(r,w,-1.,ntot2)            ! r = r - w
                                                  !      -1
            call col2(r,ml,ntot2)                 ! r = L   r
         endif
                                                  !            ______
         gamma(1) = sqrt(glsc2(r,r,ntot2))        ! gamma  = \/ (r,r)
                                                  !      1
         if(iter.eq.0) then
            div0 = gamma(1)*norm_fac
            if (param(21).lt.0) tolpss=abs(param(21))*div0
         endif

         !check for lucky convergence
         rnorm = 0.
         if(gamma(1) .eq. 0.) goto 9000
         temp = 1./gamma(1)
         call cmult2(v(1,1),r,temp,ntot2)         ! v  = r / gamma
                                                  !  1            1
         do j=1,m
            iter = iter+1
                                                  !       -1
            call col3(w,mu,v(1,j),ntot2)          ! w  = U   v
                                                  !           j

            etime2 = dnekclock()
            if(param(43).eq.1) then
               call uzprec(z(1,j),w,h1,h2,intype,wp)
            else                                  !       -1
               call hsmg_solve(z(1,j),w)          ! z  = M   w
c              call copy(z(1,j),w,ntot2)          ! z  = M   w
            endif
            etime_p = etime_p + dnekclock()-etime2

            call cdabdtp(w,z(1,j),                ! w = A z
     $                   h1,h2,h2inv,intype)      !        j

                                                  !      -1
            call col2(w,ml,ntot2)                 ! w = L   w

c           !modified Gram-Schmidt
c           do i=1,j
c              h(i,j)=glsc2(w,v(1,i),ntot2)       ! h    = (w,v )
c                                                 !  i,j       i
c              call add2s2(w,v(1,i),-h(i,j),ntot2)! w = w - h    v
c           enddo                                 !          i,j  i

c           2-PASS GS, 1st pass:

            do i=1,j
               h(i,j)=vlsc2(w,v(1,i),ntot2)       ! h    = (w,v )
            enddo                                 !  i,j       i

            call gop(h(1,j),wk1,'+  ',j)          ! sum over P procs

            do i=1,j
               call add2s2(w,v(1,i),-h(i,j),ntot2)! w = w - h    v
            enddo                                 !          i,j  i

c           2-PASS GS, 2nd pass:
c
c           do i=1,j
c              wk1(i)=vlsc2(w,v(1,i),ntot2)       ! h    = (w,v )
c           enddo                                 !  i,j       i
c                                                 !
c           call gop(wk1,wk2,'+  ',j)             ! sum over P procs
c
c           do i=1,j
c              call add2s2(w,v(1,i),-wk1(i),ntot2)! w = w - h    v
c              h(i,j) = h(i,j) + wk1(i)           !          i,j  i
c           enddo

            !apply Givens rotations to new column
            do i=1,j-1
               temp = h(i,j)
               h(i  ,j)=  c(i)*temp + s(i)*h(i+1,j)
               h(i+1,j)= -s(i)*temp + c(i)*h(i+1,j)
            enddo
                                                  !            ______
            alpha = sqrt(glsc2(w,w,ntot2))        ! alpha =  \/ (w,w)
            rnorm = 0.
            if(alpha.eq.0.) goto 900  !converged
            l = sqrt(h(j,j)*h(j,j)+alpha*alpha)
            temp = 1./l
            c(j) = h(j,j) * temp
            s(j) = alpha  * temp
            h(j,j) = l
            gamma(j+1) = -s(j) * gamma(j)
            gamma(j)   =  c(j) * gamma(j)

c            call outmat(h,m,j,' h    ',j)

            rnorm = abs(gamma(j+1))*norm_fac
            ratio = rnorm/div0
            if (ifprint.and.nid.eq.0)
     $         write (6,66) iter,tolpss,rnorm,div0,ratio,istep
   66       format(i5,1p4e12.5,i8,' Divergence')

#ifndef TST_WSCAL
            if (rnorm .lt. tolpss) goto 900  !converged
#else
            if (iter.gt.param(151)-1) goto 900
#endif
            if (j.eq.m) goto 1000 !not converged, restart

            temp = 1./alpha
            call cmult2(v(1,j+1),w,temp,ntot2)   ! v    = w / alpha
                                                 !  j+1
         enddo
  900    iconv = 1
 1000    continue
         !back substitution
         !     -1
         !c = H   gamma
         do k=j,1,-1
            temp = gamma(k)
            do i=j,k+1,-1
               temp = temp - h(k,i)*c(i)
            enddo
            c(k) = temp/h(k,k)
         enddo
         !sum up Arnoldi vectors
         do i=1,j
            call add2s2(x,z(1,i),c(i),ntot2)     ! x = x + c  z
                                                 !          i  i
         enddo
c        if(iconv.eq.1) call dbg_write(x,nx2,ny2,nz2,nelv,'esol',3)
      enddo
 9000 continue
c
      divex = rnorm
c     iter = iter - 1
c
c     DIAGNOSTICS
c      call copy   (w,x,ntot2)
       call ortho  (w) ! Orthogonalize wrt null space, if present
c      call copy(r,res,ntot2) !r = res
c      call cdabdtp(r,w,h1,h2,h2inv,intype)  ! r = A w
c      do i=1,ntot2
c         r(i) = res(i) - r(i)               ! r = res - r
c      enddo
c      call uzawa_gmres_temp(r,bm2inv,ntot2)
c                                               !            ______
c      gamma(1) = sqrt(glsc2(r,r,ntot2)/volvm2) ! gamma  = \/ (r,r)
c                                               !      1
c      print *, 'GMRES end resid:',gamma(1)
c     END DIAGNOSTICS
      call copy(res,x,ntot2)

      call ortho (res)  ! Orthogonalize wrt null space, if present

      etime1 = dnekclock()-etime1
      if (nid.eq.0) write(6,9999) istep,iter,divex,tolpss,div0,etime_p,
     &                            etime1
c     call flush_hack
 9999 format(i10,' U-PRES gmres:',I7,1p5e12.4)

      return
      end

c-----------------------------------------------------------------------

      subroutine uzawa_gmres_split0(l,u,b,binv,n)
      integer n
      real l(n),u(n),b(n),binv(n)
      integer i
      do i=1,n
         l(i)=sqrt(binv(i))
         u(i)=sqrt(b(i))
         if(abs(u(i)*l(i)-1.0).gt.1e-13) print *, i, u(i)*l(i)
      enddo
      return
      end

c-----------------------------------------------------------------------
      subroutine uzawa_gmres_split(l,u,b,binv,n)
      integer n
      real l(n),u(n),b(n),binv(n)
      integer i
      do i=1,n
c        l(i)=sqrt(binv(i))
c        u(i)=sqrt(b(i))

c        u(i)=sqrt(b(i))
c        l(i)=1./u(i)

c        l(i)=sqrt(binv(i))
         l(i)=1.
         u(i)=1./l(i)

c        if(abs(u(i)*l(i)-1.0).gt.1e-13)write(6,*) i,u(i)*l(i),' gmr_sp'
      enddo
      return
      end

c-----------------------------------------------------------------------
      subroutine uzawa_gmres_temp(a,b,n)
      integer n
      real a(n),b(n)
      integer i
      do i=1,n
         a(i)=sqrt(b(i))*a(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine ax(w,x,h1,h2,n)
      include 'SIZE'
      include 'TOTAL'

c
c     w = A*x for pressure iteration
c

      integer n
      real w(n),x(n),h1(n),h2(n)

      imsh = 1
      isd  = 1
      call axhelm (w,x,h1,h2,imsh,isd)
      call dssum  (w,nx1,ny1,nz1)
      call col2   (w,pmask,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine hmh_gmres(res,h1,h2,wt,iter)

c     Solve the Helmholtz equation by right-preconditioned
c     GMRES iteration.

      include 'SIZE'
      include 'TOTAL'
      include 'FDMH1'
      include 'GMRES'
      common  /ctolpr/ divex
      common  /cprint/ ifprint
      logical          ifprint
      real             res  (lx1*ly1*lz1*lelv)
      real             h1   (lx1,ly1,lz1,lelv)
      real             h2   (lx1,ly1,lz1,lelv)
      real             wt   (lx1,ly1,lz1,lelv)

      common /scrcg/ d(lx1*ly1*lz1*lelv),wk(lx1*ly1*lz1*lelv)

      common /cgmres1/ y(lgmres)
      common /ctmp0/   wk1(lgmres),wk2(lgmres)
      real alpha, l, temp
      integer outer

      logical iflag,if_hyb
      save    iflag,if_hyb
c     data    iflag,if_hyb  /.false. , .true. /
      data    iflag,if_hyb  /.false. , .false. /
      real    norm_fac
      save    norm_fac

      real*8 etime1,dnekclock

      n = nx1*ny1*nz1*nelv

      etime1 = dnekclock()
      etime_p = 0.
      divex = 0.
      iter  = 0
      m     = lgmres

      if(.not.iflag) then
         iflag=.true.
         call uzawa_gmres_split(ml,mu,bm1,binvm1,nx1*ny1*nz1*nelv)
         norm_fac = 1./sqrt(volvm1)
      endif

      if (param(100).ne.2) call set_fdm_prec_h1b(d,h1,h2,nelv)

      call chktcg1(tolps,res,h1,h2,pmask,vmult,1,1)
      if (param(21).gt.0.and.tolps.gt.abs(param(21)))
     $   tolps = abs(param(21))
      if (istep.eq.0) tolps = 1.e-4
      tolpss = tolps
c
      iconv = 0
      call rzero(x,n)

      outer = 0
      do while (iconv.eq.0.and.iter.lt.500)
         outer = outer+1

         if(iter.eq.0) then               !      -1
            call col3(r,ml,res,n)         ! r = L  res
c           call copy(r,res,n)
         else
            !update residual
            call copy  (r,res,n)                  ! r = res
            call ax    (w,x,h1,h2,n)              ! w = A x
            call add2s2(r,w,-1.,n)                ! r = r - w
                                                  !      -1
            call col2(r,ml,n)                     ! r = L   r
         endif
                                                  !            ______
         gamma(1) = sqrt(glsc3(r,r,wt,n))         ! gamma  = \/ (r,r)
                                                  !      1
         if(iter.eq.0) then
            div0 = gamma(1)*norm_fac
            if (param(21).lt.0) tolpss=abs(param(21))*div0
         endif

         !check for lucky convergence
         rnorm = 0.
         if(gamma(1) .eq. 0.) goto 9000
         temp = 1./gamma(1)
         call cmult2(v(1,1),r,temp,n)             ! v  = r / gamma
                                                  !  1            1
         do j=1,m
            iter = iter+1
                                                  !       -1
            call col3(w,mu,v(1,j),n)              ! w  = U   v
                                                  !           j

c . . . . . Overlapping Schwarz + coarse-grid . . . . . . .

            etime2 = dnekclock()

c           if (outer.gt.2) if_hyb = .true.       ! Slow outer convergence
            if (ifmgrid) then
               call h1mg_solve(z(1,j),w,if_hyb)   ! z  = M   w
            else                                  !  j
               kfldfdm = ndim+1
               if (param(100).eq.2) then
                   call h1_overlap_2 (z(1,j),w,pmask)
               else
                   call fdm_h1
     $               (z(1,j),w,d,pmask,vmult,nelv,ktype(1,1,kfldfdm),wk)
               endif
               call crs_solve_h1 (wk,w)           ! z  = M   w
               call add2         (z(1,j),wk,n)    !  j
            endif

            call ortho        (z(1,j)) ! Orthogonalize wrt null space, if present
            etime_p = etime_p + dnekclock()-etime2
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

            call ax  (w,z(1,j),h1,h2,n)           ! w = A z
                                                  !        j

                                                  !      -1
            call col2(w,ml,n)                     ! w = L   w

c           !modified Gram-Schmidt

c           do i=1,j
c              h(i,j)=glsc3(w,v(1,i),wt,n)        ! h    = (w,v )
c                                                 !  i,j       i

c              call add2s2(w,v(1,i),-h(i,j),n)    ! w = w - h    v
c           enddo                                 !          i,j  i

c           2-PASS GS, 1st pass:

            do i=1,j
               h(i,j)=vlsc3(w,v(1,i),wt,n)        ! h    = (w,v )
            enddo                                 !  i,j       i

            call gop(h(1,j),wk1,'+  ',j)          ! sum over P procs

            do i=1,j
               call add2s2(w,v(1,i),-h(i,j),n)    ! w = w - h    v
            enddo                                 !          i,j  i

c           2-PASS GS, 2nd pass:
c
c           do i=1,j
c              wk1(i)=vlsc3(w,v(1,i),wt,n)        ! h    = (w,v )
c           enddo                                 !  i,j       i
c                                                 !
c           call gop(wk1,wk2,'+  ',j)             ! sum over P procs
c
c           do i=1,j
c              call add2s2(w,v(1,i),-wk1(i),n)    ! w = w - h    v
c              h(i,j) = h(i,j) + wk1(i)           !          i,j  i
c           enddo

            !apply Givens rotations to new column
            do i=1,j-1
               temp = h(i,j)
               h(i  ,j)=  c(i)*temp + s(i)*h(i+1,j)
               h(i+1,j)= -s(i)*temp + c(i)*h(i+1,j)
            enddo
                                                 !            ______
            alpha = sqrt(glsc3(w,w,wt,n))        ! alpha =  \/ (w,w)
            rnorm = 0.
            if(alpha.eq.0.) goto 900  !converged
            l = sqrt(h(j,j)*h(j,j)+alpha*alpha)
            temp = 1./l
            c(j) = h(j,j) * temp
            s(j) = alpha  * temp
            h(j,j) = l
            gamma(j+1) = -s(j) * gamma(j)
            gamma(j)   =  c(j) * gamma(j)

            rnorm = abs(gamma(j+1))*norm_fac
            ratio = rnorm/div0
            if (ifprint.and.nid.eq.0)
     $         write (6,66) iter,tolpss,rnorm,div0,ratio,istep
   66       format(i5,1p4e12.5,i8,' Divergence')

#ifndef TST_WSCAL
            if (rnorm .lt. tolpss) goto 900  !converged
#else
            if (iter.gt.param(151)-1) goto 900
#endif
            if (j.eq.m) goto 1000 !not converged, restart

            temp = 1./alpha
            call cmult2(v(1,j+1),w,temp,n)   ! v    = w / alpha
                                             !  j+1
         enddo
  900    iconv = 1
 1000    continue
         !back substitution
         !     -1
         !c = H   gamma
         do k=j,1,-1
            temp = gamma(k)
            do i=j,k+1,-1
               temp = temp - h(k,i)*c(i)
            enddo
            c(k) = temp/h(k,k)
         enddo
         !sum up Arnoldi vectors
         do i=1,j
            call add2s2(x,z(1,i),c(i),n)     ! x = x + c  z
         enddo                               !          i  i
c        if(iconv.eq.1) call dbg_write(x,nx1,ny1,nz1,nelv,'esol',3)
      enddo
 9000 continue

      divex = rnorm
      call copy(res,x,n)

      call ortho   (res) ! Orthogonalize wrt null space, if present

      etime1 = dnekclock()-etime1
      if (nid.eq.0) write(6,9999) istep,iter,divex,tolpss,div0,etime_p,
     &                            etime1,if_hyb
c     call flush_hack
 9999 format(i9,' PRES gmres:',i5,1p5e12.4,1x,l4)

      if (outer.le.2) if_hyb = .false.

      return
      end
c-----------------------------------------------------------------------
      subroutine set_overlap2
c
c     Sets up the gather scatter and the SEM operators
c
      include 'SIZE'
      include 'TOTAL'

      common /c_is1/ glo_num(lxs*lys*lzs*lelv)
      integer*8 glo_num
      common /ivrtx/ vertex ((2**ldim)*lelt)
      common /handle/ gsh_dd
      integer vertex,gsh_dd

      mz = ndim-2
      nx = nx1+2
      ny = ny1+2
      nz = nz1+2*mz
      call setupds_no_crn(gsh_dd,nx,ny,nz,nelv,nelgv,vertex,glo_num)
      call swap_lengths ! Set up Matrices for FDM
      call gen_fast_g

      return
      end
c-----------------------------------------------------------------------
      subroutine h1_overlap_2(u,v,mask)
c
c     Local overlapping Schwarz solves with overlap of 2
c
      include 'SIZE'
      include 'TOTAL'

      common /cwork1/ v1(lxs,lys,lzs,lelt)
      common /handle/ gsh_dd
      integer gsh_dd

      real u(lx1,lx1,lz1,1),v(1),mask(1)
      integer e

      n = nx1*ny1*nz1*nelfld(ifield)

      call dd_swap_vals(v1,v,gsh_dd)

      iz1 = 0
      if (if3d) iz1=1
      do ie=1,nelfld(ifield)
         do iz=1,nz1
         do iy=1,ny1
         do ix=1,nx1
            u(ix,iy,iz,ie) = v1(ix+1,iy+1,iz+iz1,ie)
         enddo
         enddo
         enddo
      enddo

      call dssum (u,nx1,ny1,nz1)
      call col2  (u,mask,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine dd_swap_vals(v1,v0,gsh_dd)

      include 'SIZE'
      include 'TOTAL'

      common /work1/ w1(lxs,lys,lzs)           ! work arrarys for locals
     $              ,w2(lxs,lys,lzs)
      integer gsh_dd

      real v1(lxs,lys,lzs,lelt)
      real v0(lx1,ly1,lz1,lelt)

      integer e

      mz = ndim-2

      nx = nx1+2
      ny = ny1+2
      nz = nz1+2*mz

      n  = nx1*ny1*nz1*nelv
      m  = (nx1+2)*(ny1+2)*(nz1+2*mz)*nelv

      do e=1,nelv
         call rzero_g        (v1,e,nx,ny,nz)
         call fill_interior_g(v1,v0,e,nx1,nz1,mz,nelv)    ! v0      --> v1(int)
         call dface_ext_g    (v1,2,e,nx,nz)               ! v1(int) --> v1(face)
      enddo
c
c                 ~ ~T
c     This is the Q Q  part
c
      call gs_op(gsh_dd,v1,1,1,0)  ! 1 ==> +         ! swap v1 & add vals

      do e =1,nelv
         call dface_add1si_g  (v1,-1.,2,e,nx,nz)
         call fastdm1_g       (v1(1,1,1,e),e,w1,w2)
         call s_face_to_int2_g(v1,-1.,2,e,nx,nz)
      enddo
c
c     Exchange/add elemental solutions
c
      call gs_op              (gsh_dd,v1,1,1,0)
      do e =1,nelv
         call s_face_to_int2_g(v1,1.,2,e,nx,nz)
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine gen_fast_g
c
c     Generate fast diagonalization matrices for each element
c
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'SOLN'
      include 'WZ'

      parameter (lxss=lxs*lxs)
      common /fastg/  sr(lxss,2,lelv),ss(lxss,2,lelv),st(lxss,2,lelv)
     $             ,  df(lxs*lys*lzs,lelv)

      common /ctmpf/  lr(2*lx1+4),ls(2*lx1+4),lt(2*lx1+4)
     $              , llr(lelt),lls(lelt),llt(lelt)
     $              , lmr(lelt),lms(lelt),lmt(lelt)
     $              , lrr(lelt),lrs(lelt),lrt(lelt)
      real lr ,ls ,lt
      real llr,lls,llt
      real lmr,lms,lmt
      real lrr,lrs,lrt

      integer lbr,rbr,lbs,rbs,lbt,rbt

      call load_semhat_weighted   !   Fills the SEMHAT arrays

      ierr = 0
      do ie=1,nelv

         call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,3,ierr)
c
c        Set up matrices for each element.
c
         call set_up_fast_1D_sem_g( sr(1,1,ie),lr,nr ,lbr,rbr
     $                      ,llr(ie),lmr(ie),lrr(ie),ie)
         call set_up_fast_1D_sem_g( ss(1,1,ie),ls,ns ,lbs,rbs
     $                      ,lls(ie),lms(ie),lrs(ie),ie)
         if (if3d) then
            call set_up_fast_1D_sem_g( st(1,1,ie),lt,nt ,lbt,rbt
     $                      ,llt(ie),lmt(ie),lrt(ie),ie)
         endif
c
c        Set up diagonal inverse
c
         if (if3d) then
            eps = 1.e-5 * (vlmax(lr(2),nr-2)
     $                  +  vlmax(ls(2),ns-2) + vlmax(lt(2),nt-2))
            l   = 1
            do k=1,nt
            do j=1,ns
            do i=1,nr
               diag = lr(i) + ls(j) + lt(k)
               if (diag.gt.eps) then
                  df(l,ie) = 1.0/diag
               else
c                 write(6,3) ie,'Reset Eig in gen fast:',i,j,k,l
c    $                         ,eps,diag,lr(i),ls(j),lt(k)
c   3             format(i6,1x,a21,4i5,1p5e12.4)
                  df(l,ie) = 0.0
               endif
               l = l+1
            enddo
            enddo
            enddo
         else
            eps = 1.e-5*(vlmax(lr(2),nr-2) + vlmax(ls(2),ns-2))
            l   = 1
            do j=1,ns
            do i=1,nr
               diag = lr(i) + ls(j)

               if (diag.gt.eps) then
                  df(l,ie) = 1.0/diag

               else
c                 write(6,2) ie,'Reset Eig in gen fast:',i,j,l
c    $                         ,eps,diag,lr(i),ls(j)
c   2             format(i6,1x,a21,3i5,1p4e12.4)
                  df(l,ie) = 0.0
               endif
               l = l+1
            enddo
            enddo
         endif
c
c        Next element ....
c
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine set_up_fast_1D_sem_g(s,lam,n,lbc,rbc,ll,lm,lr,ie)
      include 'SIZE'
      include 'SEMHAT'

      parameter (lr3=2*lxs*lxs)
      common /fast1dsem/ g(lr3),w(lr3)

      real g,w
      real s(1),lam(1),ll,lm,lr
      integer lbc,rbc

      integer bb0,bb1,eb0,eb1,n,n1
      logical l,r

      n=nx1-1
      !bcs on E are from normal vel component
      if(lbc.eq.2 .or. lbc.eq.3) then !wall,sym  - Neumann
         eb0=0
         bb0=0
      else !outflow,element                      - Dirichlet
         eb0=1
         bb0=1
      endif
      if(rbc.eq.2 .or. rbc.eq.3) then !wall,sym  - Neumann
         eb1=n
         bb1=n
      else !outflow,element                      - Dirichlet
         eb1=n-1
         bb1=n-1
      endif
      l = (lbc.eq.0)
      r = (rbc.eq.0)

c     calculate A tilde operator
      call set_up_fast_1D_sem_op_a(s,eb0,eb1,l,r,ll,lm,lr,ah)
c     calculate B tilde operator
      call set_up_fast_1D_sem_op_b(g,bb0,bb1,l,r,ll,lm,lr,bh)

      n=n+3
c     call outmat   (s,n,n,'  A   ',ie)
c     call outmat   (g,n,n,'  B   ',ie)
      call generalev(s,g,lam,n,w)
      if(.not.l) call row_zero(s,n,n,1)
      if(.not.r) call row_zero(s,n,n,n)
      call transpose(s(n*n+1),n,s,n) ! compute the transpose of s
c     call outmat   (s,n,n,'  S   ',ie)
c     call outmat   (s(n*n+1),n,n,'  St  ',1)
c     call exitt

      return
      end
c-----------------------------------------------------------------------
      subroutine set_up_fast_1D_sem_op_a(g,b0,b1,l
     $                                   ,r,ll,lm,lr,ah)
c            -1 T
c     G = J B  J
c
c     gives the inexact restriction of this matrix to
c     an element plus two node on either side
c
c     g - the output matrix
c     b0, b1 - the range for Bhat indices for the element
c              (enforces boundary conditions)
c     l, r - whether there is a left or right neighbor
c     ll,lm,lr - lengths of left, middle, and right elements
c     ah - hat matrix for A or B
c
c     result is inexact because:
c        neighbor's boundary condition at far end unknown
c        length of neighbor's neighbor unknown
c        (these contribs should be small for large N and
c         elements of nearly equal size)
c
      include 'SIZE'
      real g(0:lx1+1,0:lx1+1)
      real ah(0:lx1-1,0:lx1-1)
      real ll,lm,lr
      integer b0,b1
      logical l,r

      real bl,bm,br
      integer n

      n =nx1-1

c
c     compute the weight of A hat
c
      if (l) bl = 2. /ll
             bm = 2. /lm
      if (r) br = 2. /lr

      call rzero(g,(n+3)*(n+3))
      do j=1,n+1
         do i=1,n+1
            g(i,j) = ah(i-1,j-1)*bm
         enddo
      enddo

      if (l) then
         g(0,0) = g(0,0) + bl*ah(n-1,n-1)
         g(0,1) = g(0,1) + bl*ah(n-1,n  )
         g(1,0) = g(1,0) + bl*ah(n  ,n-1)
         g(1,1) = g(1,1) + bl*ah(n  ,n  )
      elseif (b0.eq.0) then          !Neumann BC
         g(0,0) = 1.
      else                           !Dirichlet BC
         g(0,0)     = 1.
         g(1,1)     = 1.
         do i=2,n+1
            g(i,1)  = 0.
            g(1,i)  = 0.
         enddo
      endif

      if (r) then
         g(n+1,n+1) = g(n+1,n+1) + br*ah(0,0)
         g(n+1,n+2) = g(n+1,n+2) + br*ah(0,1)
         g(n+2,n+1) = g(n+2,n+1) + br*ah(0,1)
         g(n+2,n+2) = g(n+2,n+2) + br*ah(1,1)
      elseif (b1.eq.n) then          !Neumann BC
         g(n+2,n+2)=1.
      else                           !Dirichlet BC
         g(n+2,n+2)   = 1.
         g(n+1,n+1)   = 1.
         do i=1,n
            g(i,n+1)  = 0.
            g(n+1,i)  = 0.
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine set_up_fast_1D_sem_op_b(g,b0,b1,l
     $                                   ,r,ll,lm,lr,bh)
c            -1 T
c     G = J B  J
c
c     gives the inexact restriction of this matrix to
c     an element plus two node on either side
c
c     g - the output matrix
c     b0, b1 - the range for Bhat indices for the element
c              (enforces boundary conditions)
c     l, r - whether there is a left or right neighbor
c     ll,lm,lr - lengths of left, middle, and right elements
c     bh - hat matrix for  B
c
c     result is inexact because:
c        neighbor's boundary condition at far end unknown
c        length of neighbor's neighbor unknown
c        (these contribs should be small for large N and
c         elements of nearly equal size)
c
      include 'SIZE'
      real g(0:lx1+1,0:lx1+1)
      real bh(0:lx1-1)
      real ll,lm,lr
      integer b0,b1
      logical l,r

      real bl,bm,br
      integer n

      n =nx1-1
c
c     compute the weight of B hat
c
      if (l) bl = ll / 2.
             bm = lm / 2.
      if (r) br = lr / 2.

      call rzero(g,(n+3)*(n+3))
      do i=1,n+1
            g(i,i) = bh(i-1)*bm
      enddo

      if (l) then
         g(0,0) = g(0,0) + bl*bh(n-1)
         g(1,1) = g(1,1) + bl*bh(n  )
      elseif (b0.eq.0) then          !Neumann BC
         g(0,0) = 1.
      else                           !Dirichlet BC
         g(0,0)     = 1.
         g(1,1)     = 1.
         do i=2,n+1
            g(i,1)  = 0.
            g(1,i)  = 0.
         enddo
      endif

      if (r) then
         g(n+1,n+1) = g(n+1,n+1) + br*bh(0)
         g(n+2,n+2) = g(n+2,n+2) + br*bh(1)
      elseif (b1.eq.n) then          !Neumann BC
         g(n+2,n+2)=1.
      else                           !Dirichlet BC
         g(n+2,n+2)   = 1.
         g(n+1,n+1)   = 1.
         do i=1,n
            g(i,n+1)  = 0.
            g(n+1,i)  = 0.
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine fill_interior_g(v1,v,e,nx,nz,iz1,nel)

      real v1(nx+2,nx+2,nz+2*iz1,nel) ! iz1=ndim-2
      real v (nx  ,nx  ,nz      ,nel)
      integer e

      ny = nx

      do iz=1,nz
      do iy=1,ny
      do ix=1,nx
         v1(ix+1,iy+1,iz+iz1,e) = v(ix,iy,iz,e)
      enddo
      enddo
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine dface_ext_g(x,t,e,nx,nz)
c     Extend interior to face of element

      include 'SIZE'
      include 'INPUT'
      real x(nx,nx,nz,1)
      integer e,t

      ny = nx

      if (if3d) then
        do iz=2,nz-1
        do ix=2,nx-1
            x(ix,1 ,iz,e) = x(ix, 1+t,iz,e)
            x(ix,ny,iz,e) = x(ix,ny-t,iz,e)
        enddo
        enddo

        do iz=2,nz-1
        do iy=2,ny-1
            x(1 ,iy,iz,e) = x( 1+t,iy,iz,e)
            x(nx,iy,iz,e) = x(nx-t,iy,iz,e)
        enddo
        enddo

        do iy=2,ny-1
        do ix=2,nx-1
            x(ix,iy,1 ,e) = x(ix,iy, 1+t,e)
            x(ix,iy,nz,e) = x(ix,iy,nz-t,e)
        enddo
        enddo
      else
         do ix=2,nx-1
            x(ix,1 ,1,e) = x(ix, 1+t,1,e)
            x(ix,ny,1,e) = x(ix,ny-t,1,e)
         enddo
         do iy=2,ny-1
            x(1 ,iy,1,e) = x( 1+t,iy,1,e)
            x(nx,iy,1,e) = x(nx-t,iy,1,e)
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine dface_add1si_g(x,c,t,e,nx,nz)
c     Scale interior and add to face of element

      include 'SIZE'
      include 'INPUT'
      real x(nx,nx,nz,1)

      integer e,t

      ny = nx

      if (if3d) then

         do iz=2,nz-1
         do ix=2,nx-1
            x(ix,1 ,iz,e) = x(ix,1 ,iz,e) + c*x(ix, 1+t,iz,e)
            x(ix,ny,iz,e) = x(ix,ny,iz,e) + c*x(ix,ny-t,iz,e)
         enddo
         enddo

         do iz=2,nz-1
         do iy=2,ny-1
            x(1 ,iy,iz,e) = x(1 ,iy,iz,e) + c*x( 1+t,iy,iz,e)
            x(nx,iy,iz,e) = x(nx,iy,iz,e) + c*x(nx-t,iy,iz,e)
         enddo
         enddo

         do iy=2,ny-1
         do ix=2,nx-1
            x(ix,iy,1 ,e) = x(ix,iy,1 ,e) + c*x(ix,iy, 1+t,e)
            x(ix,iy,nz,e) = x(ix,iy,nz,e) + c*x(ix,iy,nz-t,e)
         enddo
         enddo

      else  ! 2D

         do ix=2,nx-1
            x(ix,1 ,1,e) = x(ix,1 ,1,e) + c*x(ix, 1+t,1,e)
            x(ix,ny,1,e) = x(ix,ny,1,e) + c*x(ix,ny-t,1,e)
         enddo
         do iy=2,ny-1
            x(1  ,iy,1,e) = x(1 ,iy,1,e) + c*x( 1+t,iy,1,e)
            x(nx ,iy,1,e) = x(nx,iy,1,e) + c*x(nx-t,iy,1,e)
         enddo

      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine fastdm1_g(R,ie,w1,w2)
c
c     Fast diagonalization solver for FEM on mesh 1
c
      include 'SIZE'

      parameter (lxss=lxs*lxs)
      common /fastg/  sr(lxss,2,lelv),ss(lxss,2,lelv),st(lxss,2,lelv)
     $             ,  df(lxs*lys*lzs,lelv)

      parameter (lxyz = lxs*lys*lzs)

      real r(1),w1(1),w2(1)

      nx = nx1+2
c
c      T
c     S  r
      call tensr3 (w1,nx,r ,nx,sr(1,2,ie),ss(1,1,ie),st(1,1,ie),w2)
c
c
c      -1 T
c     D  S  r
c
      call col2   (w1,df(1,ie),lxyz)
c
c
c        -1 T
c     S D  S  r
c
      call tensr3 (r ,nx,w1,nx,sr(1,1,ie),ss(1,2,ie),st(1,2,ie),w2)

      return
      end
c-----------------------------------------------------------------------
      subroutine s_face_to_int2_g(x,c,t,e,nx,nz)
c
c     Scale face and add to interior of element
c
      include 'SIZE'
      include 'INPUT'
      real x(nx,nx,nz,1)
      integer t,e

      ny=nx

      if (if3d) then

         do iz=2,nz-1
         do ix=2,nx-1
            x(ix, 1+t,iz,e) = c*x(ix,1 ,iz,e) + x(ix, 1+t,iz,e)
            x(ix,ny-t,iz,e) = c*x(ix,ny,iz,e) + x(ix,ny-t,iz,e)
         enddo
         enddo

         do iz=2,nz-1
         do iy=2,ny-1
            x( 1+t,iy,iz,e) = c*x(1 ,iy,iz,e) + x( 1+t,iy,iz,e)
            x(nx-t,iy,iz,e) = c*x(nx,iy,iz,e) + x(nx-t,iy,iz,e)
         enddo
         enddo

         do iy=2,ny-1
         do ix=2,nx-1
            x(ix,iy, 1+t,e) = c*x(ix,iy,1 ,e) + x(ix,iy, 1+t,e)
            x(ix,iy,nz-t,e) = c*x(ix,iy,nz,e) + x(ix,iy,nz-t,e)
         enddo
         enddo

      else
c     2D
         do ix=2,nx-1
            x(ix, 1+t,1,e) = c*x(ix,1 ,1,e) + x(ix, 1+t,1,e)
            x(ix,ny-t,1,e) = c*x(ix,ny,1,e) + x(ix,ny-t,1,e)
         enddo
         do iy=2,ny-1
            x( 1+t,iy,1,e) = c*x(1 ,iy,1,e) + x( 1+t,iy,1,e)
            x(nx-t,iy,1,e) = c*x(nx,iy,1,e) + x(nx-t,iy,1,e)
         enddo
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine outfldr_g(x,txt10,nx,nz,ichk)
      INCLUDE 'SIZE'
      INCLUDE 'TSTEP'
      real x(nx,nx,nz,lelt)
      character*10 txt10

      integer idum,e
      save idum
      data idum /3/

      if (idum.lt.0) return

      ny = nx

      mtot = nx*ny*nz*nelv
      if (nx.gt.8.or.nelv.gt.16) return
      xmin = glmin(x,mtot)
      xmax = glmax(x,mtot)

      nell = nelt
      rnel = nell
      snel = sqrt(rnel)+.1
      ne   = snel
      ne1  = nell-ne+1
      k = 1
      do ie=1,1
         ne = 0
         write(6,116) txt10,k,ie,xmin,xmax,istep,time,nx
         do l=12,0,-4
            write(6,117)
            do j=ny,1,-1
              if (nx.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.8) write(6,118) ((x(i,j,k,e+l),i=1,nx),e=1,4)
            enddo
         enddo
      enddo

  102 FORMAT(4(2f9.5,2x))
  103 FORMAT(4(3f9.5,2x))
  104 FORMAT(4(4f7.3,2x))
  105 FORMAT(5f9.5,10x,5f9.5)
  106 FORMAT(6f7.1,5x,6f7.1,5x,6f7.1,5x,6f7.1)
  107 FORMAT(7f8.4,5x,7f8.4)
  108 FORMAT(8f8.4,4x,8f8.4)
  118 FORMAT(8f12.9)

  116 FORMAT(  /,5X,'     ^              ',/,
     $    5X,'   Y |              ',/,
     $    5X,'     |              ',A10,/,
     $    5X,'     +---->         ','Plane = ',I2,'/',I2,2x,2e12.4,/,
     $    5X,'       X            ','Step  =',I9,f15.5,i2)
  117 FORMAT(' ')

      if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
      return
      end
c-----------------------------------------------------------------------
      subroutine outfldi_g(x,txt10,nx,nz,ichk)
      INCLUDE 'SIZE'
      INCLUDE 'TSTEP'
      integer x(nx,nx,nz,lelt)
      character*10 txt10

      integer idum,e,xmin,xmax
      save idum
      data idum /3/

      if (idum.lt.0) return

      ny = nx

      mtot = nx*ny*nz*nelv
      if (nx.gt.8.or.nelv.gt.16) return
      xmin = iglmin(x,mtot)
      xmax = iglmax(x,mtot)

      nell = nelt
      rnel = nell
      snel = sqrt(rnel)+.1
      ne   = snel
      ne1  = nell-ne+1
      k = 1
      do ie=1,1
         ne = 0
         write(6,116) txt10,k,ie,xmin,xmax,istep,time,nx
         do l=12,0,-4
            write(6,117)
            do j=ny,1,-1
              if (nx.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,nx),e=1,4)
              if (nx.eq.8) write(6,118) ((x(i,j,k,e+l),i=1,nx),e=1,4)
            enddo
         enddo
      enddo

  102 FORMAT(4(2i9,2x))
  103 FORMAT(4(3i9,2x))
  104 FORMAT(4(4i7,2x))
  105 FORMAT(5i9,10x,5i9)
  106 FORMAT(6i7,5x,6i7,5x,6i7,5x,6i7)
  107 FORMAT(7i8,5x,7i8)
  108 FORMAT(8i8,4x,8i8)
  118 FORMAT(8i12)

  116 FORMAT(  /,5X,'     ^              ',/,
     $    5X,'   Y |              ',/,
     $    5X,'     |              ',A10,/,
     $    5X,'     +---->         ','Plane = ',I2,'/',I2,2x,2i12,/,
     $    5X,'       X            ','Step  =',I9,f15.5,i2)
  117 FORMAT(' ')

      if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
      return
      end
c-----------------------------------------------------------------------
      subroutine setupds_no_crn(gs_h,nx,ny,nz,nel,melg,vertex,glo_num)
      include 'SIZE'
      include 'INPUT'
      include 'PARALLEL'
      include 'NONCON'
      integer   gs_h,vertex(1),e
      integer*8 ngv,glo_num(nx,ny,nz,nel)

      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

c     set up the global numbering
      call set_vert(glo_num,ngv,nx,nel,vertex,.false.)

c     zero out corners
      mz1 = max(1,nz-1)
      do e=1,nel
      do k=1,nz,mz1
      do j=1,ny,ny-1
      do i=1,nx,nx-1
         glo_num(i,j,k,e) = 0
      enddo
      enddo
      enddo
      enddo

      ntot = nx*ny*nz*nel

      t0 = dnekclock()
      call gs_setup(gs_h,glo_num,ntot,nekcomm,mp) ! initialize gs code
      t1 = dnekclock()

      et = t1-t0

c     call gs_chkr(glo_num)

      if (nid.eq.0) then
         write(6,1) et,nx,nel,ntot,ngv,gs_h
    1    format('   gs_init time',1pe11.4,' seconds ',i3,4i10)
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine rzero_g(a,e,nx,ny,nz)

      real a(nx,ny,nz,1)
      integer e

      do  i=1,nx
      do  j=1,ny
      do  k=1,nz
          a(i,j,k,e) = 0.0
      enddo
      enddo
      enddo

      return
      END
c-----------------------------------------------------------------------


text:
      SUBROUTINE SFCBKGD(bgmodel,t,q,height,tsfc,qsfc_i,
     &                   tdsfc_i,tdsfc_o,ter,p,IMX,JMX,KX,psfc,
     &                   nx_pr,ny_pr,istatus)

c inputs are from laps analyzed 2- 3d fields
C     INPUT       t        ANALYZED TEMPERATURE     3D
C                 q        Analyzed MIXING RATIO    3D
C                 height   ANALYZED HEIGHT          3D
C                 qsfc_i   Q Surface                2D
C                 ter      TERRAIN                  2D
C                 IMX      DIMENSION E-W
C                 JMX      DIMENSION N-S
C                 KX       NUMBER OF VERTICAL LEVELS
C                 p        LAPS Pressure levels     1D/3D (Pa)
c
c Compute surface variables on hi-res terrain using first 3D model level
c above the LAPS terrain. This routine is also called for reduced pressure
c where 'ter' is passed in as a constant array.

c
c J. Smart    09-22-98:	Original Version: This is used to compute
c                       sfc p for lgb when using NOGAPS1.0 deg since
c                       this field currently does not come with the
c                       model grids at AFWA.
c    "        02-01-99: Recompute tdsfc with new psfc and tsfc for
c                       consistency.
c    "        11-18-99: Put psfc,tsfc, and tdsfc comps into subroutine
c
      IMPLICIT NONE

      INTEGER    I
      INTEGER    IMX
      INTEGER    J
      INTEGER    JMX
      INTEGER    K
      INTEGER    KX
      INTEGER    nx_pr,ny_pr,ip,jp
      INTEGER    it, idebug, istatus, ishow_timer
      INTEGER    istat_qsfc_i, istat_tdsfc_i
      INTEGER    bgmodel

      LOGICAL    lfndz

      REAL       HEIGHT      ( IMX , JMX, KX )
      REAL       P           ( nx_pr,ny_pr,KX ) ! Pa
      REAL       PSFC        ( IMX , JMX )
      REAL       TSFC        ( IMX , JMX )
      REAL       Q           ( IMX , JMX , KX )
      REAL       rh3d        ( IMX , JMX , KX )
      REAL       rh2d        ( IMX , JMX )
      REAL       qsfc_i      ( imx,  jmx )
      REAL       tdsfc_i     ( imx,  jmx )
      REAL       tdsfc_o     ( imx,  jmx )
      REAL       T           ( IMX , JMX , KX )
      REAL       esat
      REAL       TER         ( IMX , JMX )
      REAL       make_td
      REAL       make_rh
      REAL       make_ssh
      REAL       t_ref,badflag
      REAL       tdsfc_o_min, tdsfc_o_max
      REAL       qsfc_i_min, qsfc_i_max
      REAL       p_mb
      REAL       r_missing_data
c
c if bgmodel = 6 or 8 then tdsfc_i is used
c if bgmodel = 4      then qsfc_i is rh (WFO - RUC)
c if bgmodel = 3      then qsfc_i is rh and q is rh
c if bgmodel = 9      then no surface fields input. Compute all from 3d
c                     fields. q3d used. (NOS - ETA)
c otherwise qsfc_i is used directly with a test for g/kg or dimensionless
c
      call get_r_missing_data(r_missing_data,istatus)

      write(6,*)
      write(6,*)' Subroutine sfcbkgd, bgmodel = ',bgmodel

      qsfc_i_min = minval(qsfc_i)
      qsfc_i_max = maxval(qsfc_i)
      write(6,*)' qsfc_i range = ',qsfc_i_min,qsfc_i_max
      if(qsfc_i_min .eq. r_missing_data .AND.
     1   qsfc_i_max .eq. r_missing_data       )then
         write(6,*)' NOTE: qsfc_i has missing data values'
         istat_qsfc_i = 0
      else
         istat_qsfc_i = 1
      endif

      write(6,*)' tdsfc_i range = ',minval(tdsfc_i),maxval(tdsfc_i)
      if(minval(tdsfc_i) .eq. r_missing_data .AND.
     &   maxval(tdsfc_i) .eq. r_missing_data       )then
         write(6,*)' NOTE: tdsfc_i has missing data values'
         istat_tdsfc_i = 0
      else
         istat_tdsfc_i = 1
      endif

      if(istat_qsfc_i .eq. 0 .and. istat_tdsfc_i .eq. 0)then
         write(6,*)
     1           ' WARNING: both qsfc_i and tdsfc_i have missing values'
      endif

      write(6,*)' ter range = ',minval(ter),maxval(ter)

      t_ref=-132.0
      if(bgmodel.eq.3.or.bgmodel.eq.9)then
         write(6,*)' bgmodel is ',bgmodel,' convert 3D q to rh'
         do k=1,kx
            do j=1,jmx
            do i=1,imx
               ip = min(i,nx_pr)
               jp = min(j,ny_pr)
               p_mb = p(ip,jp,k) / 100.
               rh3d(i,j,k)=make_rh(p_mb,t(i,j,k)-273.15,q(i,j,k)*1000.
     +,t_ref)*100.
            enddo
            enddo
         enddo

         badflag=0.
         write(6,*)
     1  ' Interp 3D T and RH to hi-res terrain, set tdsfc_o array to RH'
         write(6,*)' height bottom level range = '
     1             ,minval(height(:,:,1)),maxval(height(:,:,1))
         call interp_to_sfc(ter,rh3d,height,imx,jmx,kx,
     &                      badflag,tdsfc_o) ! Here tdsfc_o temporarily is RH
         call interp_to_sfc(ter,t,height,imx,jmx,kx,badflag,
     &                      tsfc)
      endif

      istatus = ishow_timer()

      write(6,*)' Compute sfc fields using 3D model data'

      idebug = 1

      do j=1,jmx
      do i=1,imx
         k=0
         lfndz=.false.
         do while(.not.lfndz)
            k=k+1
            if(height(i,j,k).gt.ter(i,j))then
               lfndz=.true.

               if(k.gt.0)then

                  call compute_sfc_bgfields(bgmodel,imx,jmx,kx,i,j,k
     &                ,ter(i,j),height,t,p,q,t_ref,psfc(i,j),tsfc(i,j)
     &                ,qsfc_i(i,j),qsfc_i_min,qsfc_i_max,istat_qsfc_i
     &                ,tdsfc_i(i,j),istat_tdsfc_i
     &                ,tdsfc_o(i,j),idebug,ip,jp,nx_pr,ny_pr)

                  idebug = 0

               endif

            endif
            if(k.eq.kx)lfndz=.true.
         enddo
      enddo
      enddo

      istatus = 1

      tdsfc_o_min = minval(tdsfc_o)
      tdsfc_o_max = maxval(tdsfc_o)
      if(tdsfc_o_max .gt. 1000.)then
          write(6,*)' ERROR: tdsfc_o is out of bounds'
          istatus = 0
      endif

      write(6,*)' tdsfc_o range = ',tdsfc_o_min,tdsfc_o_max
      write(6,*)' psfc range = ',minval(psfc),maxval(psfc)

      if(minval(psfc) .lt. 50000.)then
          write(6,*)' ERROR: psfc is out of bounds'
          istatus = 0
      endif

      write(6,*)' returning from sfcbkgd...'
      write(6,*)

      return
      end
c
c----------------------------------------------------------------------------
c
      subroutine compute_sfc_bgfields(bgm,nx,ny,nz,i,j,k,ter,height
     &,t,p,q,t_ref,psfc_pa,tsfc,qsfc,qsfc_i_min,qsfc_i_max,istat_qsfc_i
     &,tdsfc_i,istat_tdsfc_i,tdsfc
     &,idebug,ip,jp,nx_pr,ny_pr)
c
c J. Smart 11/18/99 put existing code in subroutine for other process use in laps
c
      implicit none

      integer nx,ny,nz
      integer i,j,k            !I, i,j,k coordinate of point for calculation
      integer bgm              !I, model type {if = 0, then tdsfc input = qsfc}
      integer nx_pr,ny_pr,ip,jp
      integer istatus, idebug, init_td
      integer istat_qsfc_i,istat_tdsfc_i
      data init_td/0/
      save init_td

      real   p(nx_pr,ny_pr,nz) !I, pressure of levels (Pa)
      real   ter               !I, Terrain height at i,j
      real   t_ref             !I, reference temp for library moisture conv routines

      real   height(nx,ny,nz)  !I, heights of pressure levels 3d
      real   t(nx,ny,nz)       !I, temperatures 3d
      real   q(nx,ny,nz)       !I, specific humidity 3d
      real   psfc_pa           !O, output surface pressure, pa
      real   tsfc              !I/O input sfc T, output recomputed T
      real   qsfc              !I   input sfc Q (g/kg)
      real   tdsfc_i           !I   input sfc Td
      real   tdsfc             !O   hi-res Td (K)

      real   qsfc_l            !L   surface spec hum, Input as q or computed internally (dimensionless)
      real   psfc_mb           !L

      real   tbar,tsfc_c
      real   td1,td2,tdsfc_c
      real   p_mb,p_mb_p1,p_mb_m1
      real   G,R
      real   ssh2,make_ssh,make_td
      real   dz,dzp,dtdz
      real   tvsfc,tvk,tbarv
      real   r_missing_data
      real   qsfc_i_min, qsfc_i_max

      parameter (G         = 9.8,
     &           R         = 287.04)
c
c if first guess values are missing data then return missing data
c
      call get_r_missing_data(r_missing_data,istatus)

      ip = min(i,nx_pr)
      jp = min(j,ny_pr)
      p_mb = p(ip,jp,k) / 100.

      if(tsfc.lt.500.0.and.t(i,j,k).lt.500.0) then
c
c first guess psfc_mb without moisture consideration
c
          tbar=(tsfc+t(i,j,k))*0.5
          dz=height(i,j,k)-ter
          psfc_mb=p_mb*exp(G/(R*tbar)*dz)

!         Calculate qsfc_l according to model (bgm=0 for reduced P)
!         This is used below for virtual temperature and sfc P reduction
          if(istat_qsfc_i .eq. 1)then
           if(bgm.eq.0.or.bgm.eq.6.or.
     &        bgm.eq.8.or.bgm.eq.12)then                   ! qsfc is Td
             tsfc_c  = tsfc-273.15
             tdsfc_c = qsfc-273.15
             if(tdsfc_c .lt. -200.)then ! a la ssh2 error check
                 if(init_td .le. 100)write(6,*)
     1           ' WARNING: setting qsfc_l to zero, tdsfc_c = ',tdsfc_c
                 init_td = init_td + 1
!                tdsfc_c = tsfc_c - 30.
                 qsfc_l = 0.0
             else
                 if(idebug .eq. 1)write(6,*)' qsfc is Td: calling ssh2'
                 qsfc_l=ssh2(psfc_mb,tsfc_c,tdsfc_c,t_ref)*.001 ! kg/kg
             endif

           elseif(bgm.eq.3.or.bgm.eq.4.or.bgm.eq.9)then    ! qsfc is RH
!    1                                .or.bgm.eq.13)then   ! qsfc is RH
             if(idebug .eq. 1)write(6,*)' qsfc is RH: calling make_ssh'
             qsfc_l=make_ssh(psfc_mb,tsfc-273.15
     &                      ,qsfc/100.,t_ref)*.001 !kg/kg

           else                                            ! qsfc is qsfc
             if(qsfc_i_max .gt. .050)then                 ! likely g/kg
                 qsfc_l=qsfc*.001                         ! make dimensionless
                 if(idebug .eq. 1)then
                     write(6,*)
     1                      ' compute_sfc_bgfields: qsfc_l = qsfc*.001 '
     1                        ,qsfc_l,qsfc
                     write(6,*)' thus qsfc_i is assumed g/kg'
                 endif
             else
                 qsfc_l=qsfc                              ! keep dimensionless
                 if(idebug .eq. 1)then
                     write(6,*)' compute_sfc_bgfields: qsfc_l = qsfc '
     1                        ,qsfc_l,qsfc
                     write(6,*)' thus qsfc_i is assumed dimensionless'
                 endif
             endif
           endif

           if(qsfc_l .lt. 0.0 .or. qsfc_l .gt. .050)then
             write(6,*)' WARNING in compute_sfc_bgfields, qsfc_l = '
     1                 ,qsfc_l
             psfc_mb =r_missing_data
             tsfc =r_missing_data
             tdsfc=r_missing_data
           endif

          elseif(istat_tdsfc_i .eq. 1)then ! q is missing, so use td more directly
           tsfc_c  = tsfc-273.15
           tdsfc_c = tdsfc_i-273.15
           qsfc_l=ssh2(psfc_mb,tsfc_c,tdsfc_c,t_ref)*.001 ! kg/kg

          else ! both q and td are missing
           qsfc_l = 0.

          endif ! valid range for qsfc_i

c pressure
          tvsfc=tsfc*(1.+0.608*qsfc_l)
          tvk=t(i,j,k)*(1.+0.608*q(i,j,k))
          tbarv=(tvsfc+tvk)*.5
          psfc_mb=(p_mb*exp(G/(R*tbarv)*dz))

!         if(j .eq. ny/2)then
!             write(6,101)i,tvsfc,tvk,tbarv,height(i,j,k),dz,p_mb,psfc_pa
!    1                   ,(psfc_pa/100.-p_mb) / dz
!101          format(' i,tvsfc,tvk,tbarv,ht,dz,pk,psfc',i4,6f8.2,f9.2
!    1                                                 ,e15.6)
!         endif

          if(k.gt.1)then ! calculate tsfc,qsfc

             dzp=height(i,j,k)-height(i,j,k-1)
c temp
             dtdz=(t(i,j,k-1)-t(i,j,k))/dzp
             tsfc=t(i,j,k)+dtdz*dz

c dew point temp
             td2=make_td(p_mb,t(i,j,k)-273.15,q(i,j,k)*1000.,t_ref)
             p_mb_m1 = p(ip,jp,k-1) / 100.
             td1=make_td(p_mb_m1,t(i,j,k-1)-273.15,q(i,j,k-1)*1000.
     .,t_ref)
             tdsfc=td2+((td1-td2)/dzp)*dz+273.15

             if(idebug .eq. 1)then
                 write(6,*)' k/tdsfc (C) = ',k,tdsfc
             endif
c
          else ! k=1: calculate tsfc,qsfc and psfc_mb

             dzp=height(i,j,k+1)-height(i,j,k)
             p_mb_p1 = p(ip,jp,k+1) / 100.
             td2=make_td(p_mb,t(i,j,k)-273.15
     &,q(i,j,k)*1000.,t_ref)+273.15
             td1=make_td(p_mb_p1,t(i,j,k+1)-273.15
     &,q(i,j,k+1)*1000.,t_ref)+273.15
             dtdz=(t(i,j,k)-t(i,j,k+1))/dzp
             tsfc=t(i,j,k)+dtdz*dz
             tbar=(tsfc+t(i,j,k))*0.5
             psfc_mb=p_mb*exp(G/(R*tbar)*dz)
             tdsfc=td2+((td1-td2)/dzp)*dz

             if(idebug .eq. 1)then
                 write(6,*)' k/tdsfc = ',k,tdsfc
             endif

          endif

          if(psfc_mb .ne. r_missing_data)then
             psfc_pa = psfc_mb * 100.
          else
             psfc_pa = r_missing_data
          endif

       else
          psfc_pa =r_missing_data
          tsfc =r_missing_data
          tdsfc=r_missing_data
       endif

       return

       end

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
C    @(#)chbsis.f	20.3 2/13/96
      subroutine chbsis

      include 'ipfinc/parametr.inc'

      include 'ipfinc/smallp.inc'

      double precision boundj, dritem, rl, rr, si, xofneg, xxx

      itr = itr + 1
      mosneg = 0
      howneg = 0.0
      xofneg = 0.0
      dritem = 0.0
      tol1 = tol(1)
      if (inrev.eq.1) go to 90
      if (r.eq.0.0) go to 40
      do 30 k = 1,size
      xr(k) = xr(k) - r * gr(k) * xkpos
      xxx = xr(k)
      if (abs(xxx).le.tol1) xr(k) = 0.0
      j = xbasis(k)
      if (j.le.n) go to 10
      xxx = xr(k)
      i = j - n
      si = s(i)
      if (si .eq. 0.0 .and. xxx .gt. 0.0 .or. si .eq. -1.0) xxx = -xxx
      go to 20

   10 boundj = bound(j)
      if (abs(boundj-xxx).le.tol1) xr(k) = boundj
      xxx = xr(k)
      if (xxx .le. boundj. or. boundj .eq. -1.0) go to 20
      xxx = boundj - xxx

   20 if (k.eq.neginv) xofneg = xxx
      if (xxx.ge.howneg.or.k.eq.neginv) go to 30
      mosneg = k
      dritem = 1.0
      if (xr(k).ge.0.0) dritem =-1.0
      howneg = xxx

   30 continue
      if (newy.ne.-1) go to 40
      it = inbase(newx)
      inbase(newx) = -1
      if (it.eq.-1) inbase(newx) = 0
      ixout = newx
      obj = obj - r * yaminc
      go to 120

   40 ixout = xbasis(newy)
      if (ixout.gt.n) go to 50
      inbase(ixout) = 0
      if (gr(newy) * xkpos.lt.0.0.and.newy.ne.neginv) inbase(ixout) = -1
      if (newy.eq.neginv.and.xr(newy).gt.0.0) inbase(ixout) = -1
   50 if (newx.gt.n) go to 60
      ihold = inbase(newx)
      inbase(newx) = newy
   60 xbasis(newy) = newx
      if (newx.gt.n) numslk = numslk + 1
      if (ixout.gt.n) numslk = numslk - 1
      xr(newy) = r
      if (newx.le.n) go to 70
      i = newx - n
      if (s(i).eq.-1.0) xr(newy) = -r
      go to 80

   70 if (ihold.eq.-1) xr(newy) = bound(newx) - r
   80 obj = obj - r * yaminc
   90 rr = 1.0/ gr(newy)
      do 110 l = 1,size
         if (abs(inv(newy,l)).lt.small) go to 110
         rl = inv(newy,l) * rr
         do 100 k = 1,size
            inv(k,l) = inv(k,l) - rl * gr(k)
  100    continue
         inv(newy,l) = rl
         if (inrev.ne.1) yr(l) = yr(l) - rl * yaminc * xkpos
  110 continue

  120 if (r.eq.0.0 .or. xofneg.lt.0.0 .and. newy.ne.neginv) go to 130
      neginv = mosneg
      driver = dritem

  130 return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
cc Copyright (C) 2009: Vladimir Rokhlin
cc
cc This software is being released under a modified FreeBSD license
cc (see COPYING in home directory).
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    $Date: 2010-01-05 12:57:52 -0500 (Tue, 05 Jan 2010) $
c    $Revision: 782 $
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       Printing routines
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C
C
        SUBROUTINE PRINI(IP1,IQ1)
        CHARACTER *1 MES(1), AA(1)
         save
        REAL *4 A(1)
        REAL *8 A2(1)
        REAL *8 A4(1)
ccc        INTEGER *4 IA(1)
        INTEGER IA(1)
        INTEGER *4 IA1(1)
        INTEGER *2 IA2(1)
        data IP/0/,IQ/0/
        IP=IP1
        IQ=IQ1
        RETURN

C
C
C
C
C
        ENTRY PRIN(MES,A,N)
        CALL  MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1200)(A(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1200)(A(J),J=1,N)
 1200 FORMAT(6(2X,E11.5))
         RETURN
C
C
C
C
        ENTRY PRIN2(MES,A2,N)
        CALL MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1400)(A2(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1400)(A2(J),J=1,N)
 1400 FORMAT(6(2X,E11.5))
        RETURN
C
C
C
C
        ENTRY PRINQ(MES,A4,N)
        CALL MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1500)(A4(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1500)(A4(J),J=1,N)
 1500 FORMAT(6(2X,e11.5))
        RETURN
C
C
C
C
        ENTRY PRINF(MES,IA,N)
        CALL MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1600)(IA(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1600)(IA(J),J=1,N)
 1600 FORMAT(10(1X,I7))
        RETURN
C
C
C
C
        ENTRY PRINF1(MES,IA1,N)
        CALL MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1600)(IA1(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1600)(IA1(J),J=1,N)
        RETURN
C
C
C
C
        ENTRY PRINF2(MES,IA2,N)
        CALL MESSPR(MES,IP,IQ)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,1600)(IA2(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,1600)(IA2(J),J=1,N)
        RETURN
C
C
C
C
        ENTRY PRINA(MES,AA,N)
        CALL MESSPR(MES,IP,IQ)
 2000 FORMAT(1X,80A1)
        IF(IP.NE.0 .AND. N.NE.0) WRITE(IP,2000)(AA(J),J=1,N)
        IF(IQ.NE.0 .AND. N.NE.0) WRITE(IQ,2000)(AA(J),J=1,N)
        RETURN
        END
c
c
c
c
c
        SUBROUTINE MESSPR(MES,IP,IQ)
        CHARACTER *1 MES(1),AST
        DATA AST/'*'/
C
C         DETERMINE THE LENGTH OF THE MESSAGE
C
        I=0
        DO 1400 I=1,10000
        IF(MES(I).EQ.AST) GOTO 1600
        I1=I
 1400 CONTINUE
 1600 CONTINUE
         IF ( (I1.NE.0) .AND. (IP.NE.0) )
     1     WRITE(IP,1800) (MES(I),I=1,I1)
         IF ( (I1.NE.0) .AND. (IQ.NE.0) )
     1     WRITE(IQ,1800) (MES(I),I=1,I1)
 1800 FORMAT(1X,80A1)
         RETURN
         END
c
c
c
c
c
        subroutine msgmerge(a,b,c)
        character *1 a(1),b(1),c(1),ast
        data ast/'*'/
c
        do 1200 i=1,1000
c
        if(a(i) .eq. ast) goto 1400
        c(i)=a(i)
        iadd=i
 1200 continue
c
 1400 continue
c
        do 1800 i=1,1000
c
        c(iadd+i)=b(i)
        if(b(i) .eq. ast) return
 1800 continue
        return
        end
c

max_stars_repo_name: zgimbutas/triasymq

max_issues_repo_name: zgimbutas/triasymq

max_forks_repo_name: zgimbutas/triasymq