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

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

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

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

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

      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 t, t2, t3, t4, s, s3o2

            t  = tFld(i,j)

            t2 = t*t

            t3 = t2*t

            t4 = t3*t

     &               t* ( eosF2 +

     &               t* ( eosF3 +

     &               t* ( eosF4 +

     &               t* ( eosF5 +

     &               t*   eosF6  ))))

     &           + eosJMDCFw(2)*t

     &           + eosJMDCFw(6)*t4*t

     &               t * ( eosS2 +

     &               t * ( eosS3 +

     &               t * ( eosS4 +

     &               t *   eosS5  ))))

     &               t * ( eosS7 +

     &               t *   eosS8  ))

     &           + eosJMDCSw(2)*t

     &           + eosJMDCSw(7)*t

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

            t  = tFld(i,j)

            t2 = t*t

            t3 = t2*t

            t4 = t3*t

     &                t * ( eosF2 +

     &                t * ( eosF3 +

     &                t * ( eosF4 +

     &                t *   eosF5  )))

     &           + eosJMDCKFw(2)*t

     &                t * ( eosS2 +

     &                t * ( eosS3 +

     &                t *   eosS4  )))

     &                t * ( eosS6 +

     &                t *   eosS7  ))

     &                t * ( eosP2 +

     &                t * ( eosP3 +

     &                t *   eosP4  )))

     &                t * ( eosP6 +

     &                t *   eosP7  ))

     &                t * ( eosP10 +

     &                t *   eosP11  ))

     &                t * ( eosP13 +

     &                t *   eosP14 ))

     &           + eosJMDCKSw(2)*t

     &           + eosJMDCKSw(6)*t

     &           + eosJMDCKP(2)*t

     &           + eosJMDCKP(6)*t

     &           + eosJMDCKP(10)*t

     &           + eosJMDCKP(13)*t

      _RL rs(4*imax+1,MAX_NO_THREADS),t,ro

      t = 1.0 _d 0/ ro

         vv(j,1,bi,bj) = vv(j,1,bi,bj)*t

       call SEAICE_SCALPROD(n, im+1, j, i1, vv, vv, t, myThid)

       hh(j,i,myThid) = t

          vv(k,i1,bi,bj) = vv(k,i1,bi,bj) - t*vv(k,j,bi,bj)

      call SEAICE_SCALPROD(n, im+1, i1, i1, vv, vv, t, myThid)

      t = sqrt(t)

      hh(i1,i,myThid) = t

      if (t .ne. 0.0 _d 0) then

       t = 1.0 _d 0 / t

          vv(k,i1,bi,bj) = vv(k,i1,bi,bj)*t

        t = hh(k1,i,myThid)

        hh(k1,i,myThid) = c(k1,myThid)*t+s(k1,myThid)*hh(k,i,myThid)

        hh(k,i,myThid) = -s(k1,myThid)*t+c(k1,myThid)*hh(k,i,myThid)

       t=rs(k,myThid)

        t = t-hh(k,j,myThid)*rs(j,myThid)

       rs(k,myThid) = t/hh(k,k,myThid)

       t = rs(j,myThid)

          sol(k,bi,bj) = sol(k,bi,bj) + t*w(k,j,bi,bj)

       t = rs(j,myThid)

       if (j .eq. 1)  t = t - 1.0 _d 0

          vv(k,1,bi,bj) = vv(k,1,bi,bj) + t*vv(k,j,bi,bj)

      subroutine SEAICE_SCALPROD(n,im,i1,i2,dx,dy,t,myThid)

      _RL t

      t     = 0. _d 0

      CALL GLOBAL_SUM_TILE_RL( dtemp,t,myThid )

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

      real*8 x1, x2, xs, t

      t = SQRT(-2.0*LOG(xs)/xs)

      port_rand_norm = t*x1

           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

        _RL  t

        _RL  t, s, pt, sit, ta

           t = ttemp(i,j)

           tk = 273.15 _d 0 + t

        _RL  t

           t = ttemp(i,j)

           tk = 273.15 + t

     &             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

        t=Tair(i,j,bi,bj)

        tta=t+Celsius2K

            deltap = t  - to + gamma_blk*(zref-ht)

      _RL t

      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

           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

        _RL  t

           t = ttemp(i,j)

           tk = 273.15 _d 0 + t

        _RL  t, s, pt, sit, ta

        _RL  t

           t = ttemp(i,j)

           tk = 273.15 + t

     &             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