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

       v(i,j,L) = vphy(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)

       vbar(i) = vbar(i) + v(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))  )

        delvv = v(i,Lm+1-L)-v(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))  )

        delvv = v(i,Lm+1-L)-v(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,

        rfv(i,j,L) = - rf(i,j,L) * v(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(vz,igrd,v,nchp,ijall,istrip,nlay,nn)

       v(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) = ( v(i,l)-tmpstrip(i) )

      _RL v(imt,0:Nrp1)

         v(i,Nrp1) = v(i,Nr)

            KRi_range(k) = p5 + SIGN(p5,v(i,k))

     &                     ( p5 + SIGN(p5,(Riinfty-v(i,k))) )

         zwork  = KRi_range(1) * v(i,1)

         v(i,1) = p2 * v(i,1) +

     &            KRi_range(1) * KRi_range(2) * v(i,2)

         v(i,1) = v(i,1) / zflag

            zflag = v(i,k)

            v(i,k) = p2 * v(i,k) +

     &           KRi_range(k) * KRi_range(kp1) * v(i,kp1) +

            v(i,k) = v(i,k) / zflag

     U     v,

      _RL u, v, density

              v =(vvel(i,j,kn0,bi,bj)+vvel(i,j+1,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))

      fcv   = 0.0

          fcloc     = fcv

            rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)v*maskC(i,j,k-1,bi,bj)

     &              rA(i,j,bi,bj)v

      fcref = fcv

            fcpertplus = fcv

               fcpertminus = fcv

          CALL WRITE_FLD_XYZ_RL( 'U.',suff,uVelv,myIter,myThid)

          CALL WRITE_FLD_XYZ_RL( 'V.',suff,vVelv,myIter,myThid)

          CALL WRITE_FLD_XYZ_RL( 'T.',suff,thetav,myIter,myThid)

          CALL WRITE_FLD_XYZ_RL( 'S.',suff,saltv,myIter,myThid)

          CALL WRITE_FLD_XY_RL('Eta.',suff,etaNv,myIter,myThid)

          CALL WRITE_FLD_XYZ_RL( 'W.',suff,wVelv,myIter,myThid)

            CALL WRITE_FLD_XYZ_RL('PH.',suff,totPhiHydv,myIter,myThid)

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

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

     &                 ncv, v, ldv, iparam, ipntr, workd, workl,

         call DSEUPD ( rvec, 'All', select, s, v, ldv, sigma, 

     &        bmat, n, which, nev, tol, resid, ncv, v, ldv, 

                  tmpvec1(l) = v(l,j)