&             exp( (24.2-0.085*t)*(pressc-1.0)/(83.143*tk) )

     &             exp( (16.4-0.040*t)*(pressc-1.0)/(83.143*tk) )

           dv = -29.48 + 0.1622*t + 2.608d-3*t*t

           xvalue = ( (48.8d0 - 0.53d0*t)*zdum

     &                 + (-0.00588d0 + 0.0001845d0*t)*zdum*zdum)

     &            / (188.93d0*(t + 273.15d0))

     I                       t,s,diclocal,pt,sit,ta,

     I                       t,s,diclocal,pt,sit,ta,

        _RL  t, s, pt, sit, ta

     I                       t,s,diclocal,pt,sit,ta,

        _RL  t, s, pt, sit, ta

        _RL  t, s, pt, sit, ta

        _RL  t

           t = ttemp(i,j,1,bi,bj)

           tk = 273.15 + t

        _RL  t

           t = ttemp(i,j,Klevel,bi,bj)

           tk = 273.15 + t

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

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

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

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

       ro(i,L) = pl(i,Lm+1-L)/(rgas*t(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) )

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

      if(t.gt.0.) then

      if( t.le.0. .and. t.gt.-40.0 ) then

        w = (40.0 + t)/40.0

      if( t.le.-40.0 .and. t.ge.-70.0 ) then

      if(t.lt.-70.0) then

       qx = f0+t*(f1+t*(f2+t*(f3+t*(f4+t*(f5+t*f6)))))

      dqx = g0+t*(g1+t*(g2+t*(g3+t*(g4+t*(g5+t*g6)))))

      _RL    getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq

      t = th(i,L)*plk(i,L)

      call QSAT ( t,p,qs,dqsdt,.true. )

      krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )

      kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )

     .    / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )

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

       t(i,j,L) = thphy(i,j,L,bi,bj)/p0kappa

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

     .   p,pl,ple,dpres,pkht,pkl,t,chemq,chemo3,co2,

     .   p,pl,ple,dpres,pkht,pkl,t,chemq,chemo3,co2,

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

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

      _RL t, t2, t3, t4, s, s3o2

            t  = tFld(i,j,k,bi,bj) 

            t2 = t*t

            t3 = t2*t

            t4 = t3*t

     &           + eosJMDCFw(2)*t

     &           + eosJMDCFw(6)*t4*t

     &           + eosJMDCSw(2)*t

     &           + eosJMDCSw(7)*t

      _RL t, t2, t3, t4, s, s3o2, p, p2

            t  = tFld(i,j,k,bi,bj) 

            t2 = t*t

            t3 = t2*t

            t4 = t3*t

     &           + eosJMDCKFw(2)*t

     &           + eosJMDCKSw(2)*t

     &           + eosJMDCKSw(6)*t

     &           + eosJMDCKP(2)*t

     &           + eosJMDCKP(6)*t

     &           + eosJMDCKP(10)*t

     &           + eosJMDCKP(13)*t

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

       t(i,j,L) = thphy(i,j,L,bi,bj)/p0kappa

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

      real*8 x1, x2, xs, t           

      t = sqrt(-2.0*log(xs)/xs)

      port_rand_norm = t*x1

      _RL s, t, u, v, density

      z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t)

      subroutine CUBIC( t, f, fp, ta, fa, fpa, tlower, tupper )

            if (fp.lt.0.) t = tupper

            if (fp.ge.0.) t = tlower

            if (fp.lt.0.) t = tupper

            if (fp.ge.0.) t = tlower

      if (t-ta .lt. 0.0) then

      sign = (t-ta)/abs(t-ta)

         t    = t + fp*(ta-t)/sngl(b+discri)

         if (abs((t-ta)*anum).lt.(tupper-tlower)*abs(den)) then

            t = t + anum*(ta-t)/den

            t = tupper

      t = max( t, tlower )

      t = min( t, tupper )

      double precision    t, f, fp, ta, fa, fpa, tlower, tupper