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

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

     . xlons,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) )

      _RL u, v, density

              u =(uvel(i,j,kn0,bi,bj)+uvel(i+1,j,kn0,bi,bj))*0.5 _d 0

              xoamc = xoamc + ( v*SIN(lon)-u*SIN(lat)*COS(lon))

              yoamc = yoamc + (-v*COS(lon)-u*SIN(lat)*SIN(lon))

              zoamc = zoamc +   u*COS(lat)

      _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)