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

      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 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_R4( x )

      _RS x

         if ( x .gt. smoothAbsFuncRange ) then

            smoothAbs_R4 = x

         else if ( x .lt. -smoothAbsFuncRange ) then

            smoothAbs_R4 = -x

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

      _RL function smoothAbs_R8( x )

      _RL x

         if ( x .ge. smoothAbsFuncRange ) then

            smoothAbs_R8 = x

         else if ( x .le. -smoothAbsFuncRange ) then

            smoothAbs_R8 = -x

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

                 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 )

                 x      = sqrt(xsq)

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

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

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

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