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

      s = ( (k-krd)/(kmm-krd) )

      if (s.GT.0. ) f = max( 0.0 _d 0, min( 1.0 _d 0, exp(-s)) )

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

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

     I            s,

      H(1,1,1,1) = s(1,1,1,1) + 

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

            s  = sFld(i,j)

           IF ( s .GT. 0. _d 0 ) THEN

            s3o2 = s*SQRT(s)

            s    = 0. _d 0

            rsalt =  s * ( eosS1 +

     &          + s*s  *   eosS9

     &         s*(

     &           + eosJMDCSw(9)*s*s

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

            s  = sFld(i,j)

          IF ( s .GT. 0. _d 0 ) THEN

            s3o2 = s*SQRT(s)

            s    = 0. _d 0

            bMsalt =  s * ( eosS1 +

     &             + p*s* ( eosP5 +

     &            + p2*s* ( eosP12 +

     &         s*( eosJMDCKSw(1)

     &     + p*s*( eosJMDCKP(5)

     &    + p2*s*( eosJMDCKP(12)

      _RL c(4*imax,MAX_NO_THREADS),s(4*imax,MAX_NO_THREADS)

      COMMON /SEAICE_FMRES_LOC_RL/ hh, c, s, rs, eps1

        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)

      s(i,myThid)   = hh(i1,i,myThid)/gam

      rs(i1,myThid) = -s(i,myThid)*rs(i,myThid)

     &     +           s(i,myThid)*hh(i1,i,myThid)

       rs(jj-1,myThid) = -s(jj-1,myThid)*rs(jj,myThid)

     &                 s(maxncv,2), resid(maxn), ax(maxm)

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

               print *, 'ph-ev ', j, s(j,1), s(j,2)

               s(j,1) = sqrt(s(j,1))

               call DAXPY(m, -s(j,1), u(1,j), 1, tmpvec2, 1)

               s(j,2) = dnrm2(m, tmpvec2, 1)

            call DMOUT(6, nconv, 2, s, maxncv, -6,

     &          log(1.0 - 0.001005*s))

     &          log(1.0 - 0.001005*s) +

         tmpa3 = -(0.07711*s) + (0.0041249*s15)

         tmpa3 = -(0.10018*s) + (0.0059415*s15)

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

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

        _RL  t, s, pt, sit, ta

        _RL  s

        _RL  t, s, pt, sit, ta

           s = stemp(i,j)

           is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)

           s2=s*s

           sqrts=sqrt(s)

           s15=s**1.5 _d 0

           scl=s/1.80655 _d 0

     &          s * (.025695 _d 0 - .025225 _d 0*tk100 +

     &        s * (0.023517 _d 0 - 0.023656 _d 0*tk100 +

     &          0.0118 _d 0 * s + 0.000116 _d 0*s2))

     &          0.0184 _d 0*s + 0.000118 _d 0*s2))

     &          77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk +

     &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +

     &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *

     &          (-0.65643 _d 0*invtk - 0.01844 _d 0)*s)

     &          (0.37335 _d 0*invtk - 0.05778 _d 0) * s)

     &          sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s)

     &          log(1.0 _d 0-0.001005 _d 0*s))

     &          * sqrts - 0.01615 _d 0 * s)

     &          log(1.0 _d 0 - 0.001005 _d 0*s))

     &   1.525 _d 0*sqrtis + log(1.0 _d 0 - 0.001005 _d 0*s) +

        _RL  s

           s = stemp(i,j)

           is=19.924*s/(1000.-1.005*s)

           s2=s*s

           sqrts=sqrt(s)

           s15=s**1.5

           scl=s/1.80655

     &          s * (.025695 - .025225*tk100 +

     &        s * (0.023517 - 0.023656*tk100 +

     &          0.0118 * s + 0.000116*s2))

     &          0.0184*s + 0.000118*s2))

           akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s +

     &          (148.0248 + 137.1942*sqrts + 1.62142*s) +

     &          (-24.4344 - 25.085*sqrts - 0.2474*s) *

     &          (-0.65643*invtk - 0.01844)*s)

     &          (0.37335*invtk - 0.05778) * s)

     &          sqrts + (-44.99486*invtk - 0.09984) * s)

     &          log(1.0-0.001005*s))

     &          sqrts - 0.01615 * s)

         tmpa3 = -(0.07711*s) + (0.0041249*s15)

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

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

        _RL  t, s, pt, sit, ta

        _RL  s

           s = stemp(i,j)

           is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)

           s2=s*s

           sqrts=sqrt(s)

           s15=s**1.5 _d 0

           scl=s/1.80655 _d 0

        _RL  t, s, pt, sit, ta

     &          s * (.025695 _d 0 - .025225 _d 0*tk100 +

     &        s * (0.023517 _d 0 - 0.023656 _d 0*tk100 +

     &          0.0118 _d 0 * s + 0.000116 _d 0*s2))

     &          0.0184 _d 0*s + 0.000118 _d 0*s2))

     &          77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk +

     &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +

     &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *

     &          (-0.65643 _d 0*invtk - 0.01844 _d 0)*s)

     &          (0.37335 _d 0*invtk - 0.05778 _d 0) * s)

     &          sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s)

     &          log(1.0 _d 0-0.001005 _d 0*s))

     &          * sqrts - 0.01615 _d 0 * s)

     &          log(1.0 _d 0 - 0.001005 _d 0*s))

     &   1.525 _d 0*sqrtis + log(1.0 _d 0 - 0.001005 _d 0*s) +

        _RL  s

           s = stemp(i,j)

           is=19.924*s/(1000.-1.005*s)

           s2=s*s

           sqrts=sqrt(s)

           s15=s**1.5

           scl=s/1.80655

     &          s * (.025695 - .025225*tk100 +

     &        s * (0.023517 - 0.023656*tk100 +

     &          0.0118 * s + 0.000116*s2))

     &          0.0184*s + 0.000118*s2))

           akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s +

     &          (148.0248 + 137.1942*sqrts + 1.62142*s) +

     &          (-24.4344 - 25.085*sqrts - 0.2474*s) *

     &          (-0.65643*invtk - 0.01844)*s)

     &          (0.37335*invtk - 0.05778) * s)

     &          sqrts + (-44.99486*invtk - 0.09984) * s)

     &          log(1.0-0.001005*s))

     &          sqrts - 0.01615 * s)

     &          log(1.0 - 0.001005*s))

     &          log(1.0 - 0.001005*s) +