_RL x                     ! stability function  [-]

             x      = SQRT(xsq)

     &                ( LOG( (exf_one + exf_two*x + xsq)*

     &                  -exf_two*ATAN(x) + exf_half*pi )

       _RL FUNCTION LAGRAN(i,x,a,sp)

        _RS x

          numer = numer*(x    - a(k))

      _RL  sdf(m,n),sclr(m,n),csm(m,n),x

           x=so2(i,j,k)*csm(i,j)

           df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*sqrt(x)))

         x = (1.-rsirdf(i,j))*fdifir(i,j) + (1.-rsirbm(i,j))*fdirir(i,j)

         x = (-sdf(i,j)-df(i,j,np+1)*sclr(i,j))/x

         fdirir(i,j)=fdirir(i,j)*(1.+x)

         fdifir(i,j)=fdifir(i,j)*(1.+x)

      _RL x,c,xx

      x=MAX(3.d0,xx)

      c= 1.d0-(2.d0/pi)*ASIN(1.d0/x)

      IDEMIX_gofx2 = 2.d0/pi/c*0.9d0*x**(-2.d0/3.d0)*(1.-EXP(-x/4.3d0))

      _RL FUNCTION IDEMIX_hofx1(x)

      _RL x

     &              ASIN(1.d0/MAX(1.01d0,x)))*(x-1.d0)/(x+1.d0)

     I                          xx,n,m,sNx,sNy,x,

      _RL     x(snx+1,sny+1)

          IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.

     &         ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN

        IF ( x(i,j).LT.xx(2) ) THEN

        ELSE IF ( x(i,j).GE.xx(n) ) THEN

        IF (x(i,j) .GE. xx(n)) THEN

        ELSE IF (x(i,j) .LT. xx(2)) THEN

        ELSE IF ( (x(i,j) .GE. xx(k(i,j)))

     &    .AND.   (x(i,j) .LT. xx(k(i,j)+1)) ) THEN

        ELSE IF (x(i,j) .GE. xx(k(i,j))) THEN

         DO WHILE (x(i,j) .GE. xx(k(i,j)+1))

        ELSE IF (x(i,j) .LT. xx(k(i,j)+1)) THEN

         DO WHILE (x(i,j) .LT. xx(k(i,j)))

      x = -x*V350P/(1028.1063+x)

      x = x*(1.0 - P/B) + (V350P+x)*P*SR/(B*(B+SR))

      SBO_RHO = 1028.106331 + P/B/SR - x / (SR*(SR+x))

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

      T = T + 0.29289322*(x-Q)

      Q = 0.58578644*x + 0.121320344*Q

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

      T = T + 1.707106781*(x-Q)

      Q = 3.414213562*x - 4.121320344*Q

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

      T = T + (x-2.0*Q)/6.0

      x = (4.8314E-4 * S +

      _RL x2, x, y, cp0

      x    = 0. _d 0

      if (x2.gt.0. _d 0) x = sqrt(x2)

     &     x*(937.2099110620707 _d 0 + y*(588.1802812170108 _d 0 +

     &     x*(-1687.914374187449 _d 0 + x*(246.9598888781377 _d 0 + 

     &     x*(123.59576582457964 _d 0 - 48.5891069025409 _d 0*x)) + 

      _RL sfac, x2, x, y, g03, g08

      x    = 0. _d 0

      if (x2.gt.0. _d 0) x = sqrt(x2)

      g08 = x2*(1760.062705994408  + x*(-86.1329351956084 +

     &     x*( -137.1145018408982  + y*(296.20061691375236 +

          x = xC(i,j,bi,bj)

           IF ((x .gt. min_x_noflow_SOUTH) .AND.

     &         (x .le. max_x_noflow_SOUTH)) THEN

           IF ((x .gt. min_x_noStress_SOUTH) .AND.

     &         (x .le. max_x_noStress_SOUTH)) THEN

           IF ((x .gt. min_x_FluxBdry_SOUTH) .AND.

     &         (x .le. max_x_FluxBdry_SOUTH)) THEN

      _RL x, y

           IF ((x .gt. min_x_Dirich_SOUTH) .AND.

     &         (x .le. max_x_Dirich_SOUTH)) THEN

           IF ((x .gt. min_x_CFBC_SOUTH) .AND.

     &         (x .le. max_x_CFBC_SOUTH)) THEN

           IF ((x .gt. min_x_noflow_NORTH) .AND.

     &         (x .le. max_x_noflow_NORTH)) THEN

           IF ((x .gt. min_x_noStress_NORTH) .AND.

     &         (x .le. max_x_noStress_NORTH)) THEN

           IF ((x .gt. min_x_FluxBdry_NORTH) .AND.

     &         (x .le. max_x_FluxBdry_NORTH)) THEN

           IF ((x .gt. min_x_Dirich_NORTH) .AND.

     &         (x .le. max_x_Dirich_NORTH)) THEN

           IF ((x .gt. min_x_CFBC_NORTH) .AND.

     &         (x .le. max_x_CFBC_NORTH)) THEN

      _RL slope_pos, c1, x, y, lenx, leny

           x = xC(i,j,bi,bj)

     &        (1+sin(2*streamice_kx_b_init*PI*x/lenx)*

           x = xC(i,j,bi,bj)

     &        sin(2*streamice_kx_b_init*PI*x/lenx)))

      _RL x                  ! stability function  [-]

           x      = SQRT(xsq)

     &             *( LOG( (exf_one + exf_two*x + xsq)

     &                -exf_two*ATAN(x) + exf_half*pi )

           x=(xC(I,1,bi,bj)-Xcenter)/Lrho +Width/2. _d 0

           IF ((x.GE.0.).AND.(x.LE.Width)) THEN

            hinit=Dinf*(exp(-x)) - Dmax

     &     -  sqrt(gp_inflow*Dinf)*exp(-x)

      _RL Width,x,Xcenter,Fz,zt,Rit

      _RS FUNCTION SMOOTHABS_RS( x )

      _RS x

         IF ( x .GT. smoothAbsFuncRange ) THEN

            SMOOTHABS_RS = x

         ELSEIF ( x .LT. -smoothAbsFuncRange ) THEN

            SMOOTHABS_RS = -x

            SMOOTHABS_RS = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf

      _RL FUNCTION SMOOTHABS_RL( x )

      _RL x

         IF ( x .GE. smoothAbsFuncRange ) THEN

            SMOOTHABS_RL = x

         ELSEIF ( x .LE. -smoothAbsFuncRange ) THEN

            SMOOTHABS_RL = -x

            SMOOTHABS_RL = log(.5*(exp(x*sf)+exp(-x*sf)))*rsf

          x = xguess

            x2=x*x

            x3=x2*x

            a = x3 + k1plocal*x2 + k12p*x + k123p

            da = 3.0*x2 + 2.0*k1plocal*x + k12p

            b = x2 + k1local*x + k12

            db = 2.0*x + k1local

            fn = k1local*x*diclocal/b +

     &        btlocal/(1.0 + x/kblocal) +

     &        kwlocal/x +

     &        pt*k12p*x/a +

     &        sit/(1.0 + x/ksilocal) -

     &        x/c -

     &        stlocal/(1.0 + kslocal/x/c) -

     &        ftlocal/(1.0 + kflocal/x) -

            df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -

     &        btlocal/kblocal/(1.0+x/kblocal)**2. -

     &        (pt*k12p*(a - x*da))/a2 -

     &        sit/ksilocal/(1.0+x/ksilocal)**2. +

     &        stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +

     &        ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -

     &        pt*x2*(3.0*a-x*da)/a2

            x = x + deltax

              if(hstep .eq. 1)x = xhi

              if(hstep .eq. 2)x = xlo

              if(hstep .eq. 3)x = xmid

              x2=x*x

              x3=x2*x

              a = x3 + k1plocal*x2 + k12p*x + k123p

              da = 3.0*x2 + 2.0*k1plocal*x + k12p

              b = x2 + k1local*x + k12

              db = 2.0*x + k1local

              fn = k1local*x*diclocal/b +

     &          btlocal/(1.0 + x/kblocal) +

     &          kwlocal/x +

     &          pt*k12p*x/a +

     &          sit/(1.0 + x/ksilocal) -

     &          x/c -

     &          stlocal/(1.0 + kslocal/x/c) -

     &          ftlocal/(1.0 + kflocal/x) -

          x = xmid

        htotal = x

        _RL  x

                 x      = sqrt(xsq)

     &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +

     &                     2. _d 0*atan(x) + pi*.5 _d 0)

      _RL x                  ! stability function  [-]

      _RL x                  ! stability function  [-]

                 x      = SQRT(xsq)

     &                    ( LOG( (one + two*x + xsq)*(one+xsq)*.125 )

     &                     -two*ATAN(x) + half*pi )

      _RL x,y,psik,psic,f,c

       x = (1. - 15.*zL)**.25                   !Kansas unstable

       psik=2.*LOG((1.+x)/2.)+LOG((1.+x*x)/2.)-2.*ATAN(x)+2.*ATAN(oneRL)

      _RL x,y,psik,psic,f,c

       x = (1. - 15.*zL)**.5                    !Kansas unstable

       psik = 2.*LOG((1.+x)/2.)

                 x      = sqrt(xsq)

     &                    (2. _d 0*log(5. _d -1*(1. _d 0+x)) +

     &                     2. _d 0*atan(x) + pi*.5 _d 0)

                 x      = sqrt(xsq)

      _RL huol,stable,xsq,x,psimh,psixh

      subroutine DOSTORE( n, x, store, j )

          tmpx(i) = x(i)

          x(i) = tmpx(i)

      double precision    x(n)

        _RL  x

          x = xguess

            x2=x*x

            x3=x2*x

            a = x3 + k1plocal*x2 + k12p*x + k123p

            da = 3.0*x2 + 2.0*k1plocal*x + k12p

            b = x2 + k1local*x + k12

            db = 2.0*x + k1local

            fn = k1local*x*diclocal/b +

     &        btlocal/(1.0 + x/kblocal) +

     &        kwlocal/x +

     &        pt*k12p*x/a +

     &        sit/(1.0 + x/ksilocal) -

     &        x/c -

     &        stlocal/(1.0 + kslocal/x/c) -

     &        ftlocal/(1.0 + kflocal/x) -

            df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -

     &        btlocal/kblocal/(1.0+x/kblocal)**2. -

     &        (pt*k12p*(a - x*da))/a2 -

     &        sit/ksilocal/(1.0+x/ksilocal)**2. +

     &        stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +

     &        ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -

     &        pt*x2*(3.0*a-x*da)/a2

            x = x + deltax

              if(hstep .eq. 1)x = xhi

              if(hstep .eq. 2)x = xlo

              if(hstep .eq. 3)x = xmid

              x2=x*x

              x3=x2*x

              a = x3 + k1plocal*x2 + k12p*x + k123p

              da = 3.0*x2 + 2.0*k1plocal*x + k12p

              b = x2 + k1local*x + k12

              db = 2.0*x + k1local

              fn = k1local*x*diclocal/b +

     &          btlocal/(1.0 + x/kblocal) +

     &          kwlocal/x +

     &          pt*k12p*x/a +

     &          sit/(1.0 + x/ksilocal) -

     &          x/c -

     &          stlocal/(1.0 + kslocal/x/c) -

     &          ftlocal/(1.0 + kflocal/x) -

          x = xmid

        htotal = x

      _RS y(1-OLy:sNy+OLy), x(1-OLx:sNx+OLx)

         x(i) = 0.0

         x(i) = i

         xa(i) = x(i) - MOD(x(i),10.0)

         temp2(i,j) = 0.1*( (x(i)-xa(i))*temp(INT(xa2(i)),j)+

     &        (10.0-x(i)+xa(i))*temp(INT(xa(i)),j) )