_RL u(sNx,sNy,Nrphys),v(sNx,sNy,Nrphys),t(sNx,sNy,Nrphys)

       u(i,j,L) = uphy(i,j,L,bi,bj)

     & p(1,1,bi,bj),u,v,t,q,pl(1,1,1,bi,bj),ple(1,1,1,bi,bj),

      call FIZHI_STEP_DIAG(myid,p,u,v,t,q,qq,pkl,dpres,

     &  ptop,p,u,v,t,q,ntracer,ptracer,pl,ple,dpres,pkht,pkl,

      call GWDRAG (myid,p,pl,ple,dpres,pkl,u,v,t,q,phis_var,

      call RAYLEIGH(myid,pl,pkl,pkht,surfz,u,v,t,q,im,jm,lm,

     &    p,pl,ple,dpres,pkht,pkl,u,v,t,q,bi,bj,ntracer,ptracer,qq,

      _RL u(im,jm,lm),v(im,jm,lm),t(im,jm,lm)

     & p,u,v,t,q,pl,ple,dpres,pkht,pkl,surfz,fracland,landtype,radswt,

      SUBROUTINE GWDD ( ps,u,v,t,dudt,dvdt,xdrag,ydrag,

      _RL u(irun,Lm), v(irun,Lm), t(irun,Lm)

       ubar(i) = ubar(i) + u(i,Lm+1-L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))

        velco = 0.5*( (u(i,Lm+1-L)*ubar(i) + v(i,Lm+1-L)*vbar(i))

     .            + (u(i,Lm+2-L)*ubar(i) + v(i,Lm+2-L)*vbar(i))  )

        deluu = u(i,Lm+1-L)-u(i,Lm+2-L)

        velco = 0.5*( (u(i,Lm+1-L)*ubar(i) + v(i,Lm+1-L)*vbar(i))

     .            + (u(i,Lm+2-L)*ubar(i) + v(i,Lm+2-L)*vbar(i))  )

        deluu = u(i,Lm+1-L)-u(i,Lm+2-L)

      _RL u(im,jm,lm),v(im,jm,lm),t(im,jm,lm),s(im,jm,lm)

      subroutine RAYLEIGH(myid,pres,pkap,pekap,zsurf,u,v,t,s,im,jm,lm,

        rfu(i,j,L) = - rf(i,j,L) * u(i,j,L)

        rft(i,j,L) = -(u(i,j,L)*rfu(i,j,L) + v(i,j,L)*rfv(i,j,L) )*cpinv

      _RL u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)

       call STRIP2TILE(uz,igrd,u,nchp,ijall,istrip,nlay,nn)

       u(i,nlay+1) = 0.

      call TRBFLX(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,

       wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))

      tends(i) = ( u(i,l)-tmpstrip(i) )

     &       u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) + 

     &       u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)

            uy = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) + 

     &       u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) + 

     &       u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) + 

     &       u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)

     I    u, 

      _RL u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

            uq = u(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) + 

     &       u(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) + 

     &       u(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) + 

     &       u(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)

            ux = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) + 

     &       u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) + 

      _RL u(sNx,sNy,Nrphys),v(sNx,sNy,Nrphys),t(sNx,sNy,Nrphys)

     &                 v(ldv,maxncv), u(ldu, maxnev),

               call DCOPY( m, tmpvec2, 1, u(l,j), 1 )

               temp = one / dnrm2( m, u(l,j), 1 )

               print *, 'ph-print ', j, dnrm2( m, u(l,j), 1 )

               call DSCAL(m , temp, u(l,j), 1 )

               call DAXPY(m, -s(j,1), u(1,j), 1, tmpvec2, 1)

      u = windSq(i,j)

      Du= SQRT(u + Wg**2 )  !include gustiness in wind spd. difference

      u = SQRT(u)

       Du=SQRT(u**2 + Wg**2)        !include gustiness in wind spd.

      _RL u,q,Tas,tta,zi,es,qs,tsw

      _RL recipNym1, local, u, ssqa

             u = uWind(i,j,bi,bj)**2 + vWind(i,j,bi,bj)**2

             u = SQRT(u)

             wavesp(i,j,bi,bj) = 0.729 _d 0 * u

             wavesh(i,j,bi,bj) = 0.018 _d 0 * u*u*(1. + .015 _d 0 *u)

              u = uWind(i,j,bi,bj)**2 + vWind(i,j,bi,bj)**2

              u = SQRT(u)

              wavesp(i,j,bi,bj)=0.729 _d 0 * u

              wavesh(i,j,bi,bj)=0.018 _d 0 * u*u*(1. + .015 _d 0 *u)