do L=1,lm-1

      pkmin = pkz(1,1,L)

      pkmax = pkz(1,1,L+1)

         if( pkz(i,1,L)  .lt.pkmin ) pkmin = pkz(i,1,L)

         if( pkz(i,1,L+1).gt.pkmax ) pkmax = pkz(i,1,L+1)

         if( pk.le.pkz(i,1,L+1) .and. pk.ge.pkz(i,1,L) ) then

         temp      = ( pkz(i,1,L)-pk ) / ( pkz(i,1,L)-pkz(i,1,L+1) )

                 if( qinp(i,1,L)  .ne.undef  .and. 

     .               qinp(i,1,L+1).ne.undef ) then

         qprs(i,1) = qinp(i,1,L+1)*temp + qinp(i,1,L)*(1.-temp)

            else if( qinp(i,1,L+1).ne.undef  .and. temp.ge.0.5 ) then

         qprs(i,1) = qinp(i,1,L+1)

            else if( qinp(i,1,L)  .ne.undef  .and. temp.le.0.5 ) then

         qprs(i,1) = qinp(i,1,L)

      DO 120 L=2,IPPH

         LC = IPP2-L

            C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)

               C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)

      DO 127 L=2,IPPH

         LC = IPP2-L

            CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)

               CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)

      DO 120 L=2,IPPH

         LC = IPP2-L

            C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)

               C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)

      DO 127 L=2,IPPH

         LC = IPP2-L

            CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)

               CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)

      do L=1,lm

        tmpdiag(i,j) = rfu(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'RFU     ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = rfv(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'RFV     ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = rft(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'RFT     ',L,1,3,bi,bj,myid)

      integer  i,j,L

      do L=1,lm

        dz(i,j,L) = cpog * (pekap(i,j,L+1)-pekap(i,j,L)) * t(i,j,L) *

     .                (1.+virtcon*s(i,j,L))

      do L=lm-1,1,-1

        z(i,j,L) = z(i,j,L+1) + 0.5 * (dz(i,j,L)+dz(i,j,L+1))

      do L=1,lm

        rf(i,j,L) = dampcoef*(1+tanh((z(i,j,L)-50000.)/5000.))/86400.

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

        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

     .                        /pkap(i,j,L)

      do L = 1,Lm

      tv(i,j,L) = tz(i,j,L)*pkz(i,j,L)*(1.+.609*qz(i,j,L))

      do L = 1,Lm

       dragu(i,j,L) = 0.

       dragv(i,j,L) = 0.

       dragt(i,j,L) = 0.

      do L = 1,Lm

         dragu(i,j,L) = sign( min(0.006 _d 0,abs(dragu(i,j,L))), dragu(i

     $        ,j,L) ) 

         dragv(i,j,L) = sign( min(0.006 _d 0,abs(dragv(i,j,L))), dragv(i

     $        ,j,L) ) 

      dragt(i,j,L) = -( uz(i,j,L)*dragu(i,j,L)+vz(i,j,L)*dragv(i,j,L) )

       dudt(i,j,L) = dudt(i,j,L) + dragu(i,j,L)

       dvdt(i,j,L) = dvdt(i,j,L) + dragv(i,j,L)

       dtdt(i,j,L) = dtdt(i,j,L) + dragt(i,j,L)*pz(i,j)/pkz(i,j,L)

      do L = 1,Lm

        tmpdiag(i,j) = dragu(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'GWDU    ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = dragv(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'GWDV    ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = dragt(i,j,L)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'GWDT    ',L,1,3,bi,bj,myid)

      do L = 1,Lm

      sumu(i,j) = sumu(i,j) + dragu(i,j,L)*dpres(i,j,L)/pz(i,j)

      do L = 1,Lm

      sumu(i,j) = sumu(i,j) + dragv(i,j,L)*dpres(i,j,L)/pz(i,j)

      integer  i,L

       ro(i,L) = pl(i,Lm+1-L)/(rgas*t(i,Lm+1-L))

       roiave  = ( 1./ro(i,L-1) + 1./ro(i,L) )*0.5

       dz(i,L) = agrav*roiave*( pl(i,Lm+2-L)-pl(i,Lm+1-L) )

      do  L = 1,nbase(i)-1

       robar(i) = robar(i) + ro(i,L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))

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

       vbar(i) = vbar(i) + v(i,Lm+1-L) * (ple(i,Lm+2-L)-ple(i,Lm+1-L))

        vai1 = (t(i,Lm+1-L)-t(i,Lm+2-L))/dz(i,L)+gocp

        vai2    = 2.0*grav/( t(i,Lm+1-L)+t(i,Lm+2-L) )

        bv(i,L) = sqrt(vsqua)

       nbar(i) = nbar(i) + bv(i,L)*(pl(i,Lm+2-L)-pl(i,Lm+1-L))

         tensio(i,L) = tensio(i,1)

         xtens(i,L) = xtens(i,1)

         ytens(i,L) = ytens(i,1)

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

        roave = 0.5*(ro(i,L-1)+ro(i,L)) * 100.0     

         tensio(i,L) = tensio(i,L-1)

        fro2 = bv(i,L)/(akwnmb*roave*velco*velco*velco)*tensio(i,L-1)

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

        delvv = v(i,Lm+1-L)-v(i,Lm+2-L)

         delz = dz(i,L)

         vsqua = bv(i,L)*bv(i,L)

         tensio(i,L) = tensio(i,L-1)

         tensio(i,L) = crif2/fro2*tensio(i,L-1)

         tensio(i,L) = tensio(i,L-1)

        xtens(i,L) = tensio(i,L)*cos(ang(i))

        ytens(i,L) = tensio(i,L)*sin(ang(i))

        tensio(i,L) = 0.0

         tensio(i,L) = tensio(i,L-1)

        roave = 0.5*(ro(i,L-1)+ro(i,L)) * 100.0

        vai1  = (t(i,Lm+1-L)-t(i,Lm+2-L))/dz(i,L)+gocp

         tensio(i,L) = 0.0

        vai2  = 2.0*grav/(t(i,Lm+1-L)+t(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))  )

         tensio(i,L) = 0.0

        fro2 = vaisd/(akwnmb*roave*velco*velco*velco)*tensio(i,L-1)

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

        delvv = v(i,Lm+1-L)-v(i,Lm+2-L)

         delz   = dz(i,L)

         tensio(i,L) = 0.0

         tensio(i,L) = crif2/fro2*tensio(i,L-1)

         tensio(i,L) = tensio(i,L-1)

        xtens(i,L) = tensio(i,L)*cos(ang(i))

        ytens(i,L) = tensio(i,L)*sin(ang(i))

       coef = -grav*ps(i)/dpres(i,Lm+1-L)

       dudt(i,Lm+1-L) = coef*(xtens(i,L+1)-xtens(i,L))

       dvdt(i,Lm+1-L) = coef*(ytens(i,L+1)-ytens(i,L))

         dudt(i,Lm+1-L) = dudt(i,Lm)

         dvdt(i,Lm+1-L) = dvdt(i,Lm)

       dudt(i,L) = - dudt(i,L)/ps(i)*0.01

       dvdt(i,L) = - dvdt(i,L)/ps(i)*0.01

      integer n,i,j,L

      do L=1,Nrphys

        tmpdiag(i,j) = (moistu (i,j,L,bi,bj)+turbu(i,j,L,bi,bj) )*86400

       call DIAGNOSTICS_FILL(tmpdiag,'DIABU   ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = (moistv (i,j,L,bi,bj)+turbv(i,j,L,bi,bj) )*86400

       call DIAGNOSTICS_FILL(tmpdiag,'DIABV   ',L,1,3,bi,bj,myid)

     .   ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +

     .      lwdt(i,j,L,bi,bj) +

     .      dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +

     .      swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )

     .      * pk(i,j,L,bi,bj)*pinv(i,j)*86400

       call DIAGNOSTICS_FILL(tmpdiag,'DIABT   ',L,1,3,bi,bj,myid)

     . ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) * 

       call DIAGNOSTICS_FILL(tmpdiag,'DIABQ   ',L,1,3,bi,bj,myid)

     .            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )

       call DIAGNOSTICS_FILL(tmpdiag,'RADLW   ',L,1,3,bi,bj,myid)

     .            dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )

       call DIAGNOSTICS_FILL(tmpdiag,'LWCLR   ',L,1,3,bi,bj,myid)

       call DIAGNOSTICS_FILL(tmpdiag,'RADSW   ',L,1,3,bi,bj,myid)

       call DIAGNOSTICS_FILL(tmpdiag,'SWCLR   ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = uphy(i,j,L)

       call DIAGNOSTICS_FILL(tmpdiag,'UWND    ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = vphy(i,j,L)

       call DIAGNOSTICS_FILL(tmpdiag,'VWND    ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = thphy(i,j,L)*pk(i,j,L,bi,bj)

       call DIAGNOSTICS_FILL(tmpdiag,'TMPU    ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = qq(i,j,L,bi,bj)

       call DIAGNOSTICS_FILL(tmpdiag,'TKE     ',L,1,3,bi,bj,myid)

        tmpdiag(i,j) = sphy(i,j,L) * 1000.

       call DIAGNOSTICS_FILL(tmpdiag,'SPHU    ',L,1,3,bi,bj,myid)

       do L=1,Nrphys

     .             moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

       do L=1,Nrphys

     .             turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

       do L=1,Nrphys

        qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +

     .  dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )

     .             *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

       do L=1,Nrphys

     .             swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

      integer  i,j,L

       do L=1,Nrphys

     .             sphy(i,j,L)*dp(i,j,L,bi,bj)

           cldtot(i,j,L)=min(1.0 _d 0,max(cldlw(i,j,L),fccave(i,j,L)/

           cldmxo(i,j,L) =  min( 1.0 _d 0,      clwmo(i,j,L) )

           lwlz(i,j,L) =  lwlz(i,j,L) + qliqave(i,j,L)/imstturb

        do L =1,lm

         cldtot(i,j,L) =  min( 1.0 _d 0,cldlw(i,j,L) )

         cldmxo(i,j,L) =  min( 1.0 _d 0,clwmo(i,j,L) )

       DO L = 1,lm

        ADELPL(I,L) = convrt / ( ple(I,L+1)-ple(I,L) )

       DO L = 1,lm

        CLRO(I,L) = min( 1.0 _d 0,clro(i,L) )

        CLMO(I,L) = min( 1.0 _d 0,clmo(i,L) )

      DO L = 1,lm

      TZL(I,L) = TZL(I,L)*pk(I,L)

      do L = 1,lm

      tauc(i,L,n) = tauc(i,L,n)*0.75

      do L = 1,lm

      taual(i,L,n) = 0.

      ssaal(i,L,n) = 0.

      asyal(i,L,n) = 0.

      do L = 1,lm

       cwc(i,L,n) = 0.

      reff(i,L,n) = 0.

       do L = 1,lm

         dtrad(i,L) = (   flx(i,L)-   flx(i,L+1))*adelpl(i,L)

         tmpstrip(i,L) = flx(i,L)

         dtdtg(i,L) = ( dfdts(i,L)- dfdts(i,L+1))*adelpl(i,L)

        dtradc(i,L) = (flxclr(i,L)-flxclr(i,L+1))*adelpl(i,L)

      DO L = 1,lm

      DTRAD (I,L) = ple(i,lm+1) * DTRAD (I,L)/pk(I,L)

      DTRADC(I,L) = ple(i,lm+1) * DTRADC(I,L)/pk(I,L)

       dtdtg(I,L) = ple(i,lm+1) * dtdtg (I,L)/pk(I,L)

      do L = 1,lm

         tmpdiag(i,j) = tz(i,j,L)*pkz(i,j,L)