real*8 function lagran(i,x,a,sp)

        real*8 a(4)

          denom = denom*(a(i) - a(k))

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

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

            a2=a*a

     &        pt*k12p*x/a +

     &        2.0*pt*k123p/a +

     &        pt*x3/a -

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

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

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

              a2=a*a

     &          pt*k12p*x/a +

     &          2.0*pt*k123p/a +

     &          pt*x3/a -

        _RL  a

      _RL a,g,agrav,akwnmb

      a      = 1.0

        gstar = g*frsf*frsf/(frsf*frsf+a*a)

      subroutine STRIPIT(a,b,irun,ia,ib,l,k)

      _RL a(ia,l), b(ib,l)

        b(i,j) = a(i+offset-1,j)

        b(i,j) = a(i+offset-1,j)

        b(Lena+i,j) = a(Lena+offset-1,j)

      subroutine STRIPITINT(a,b,irun,ia,ib,l,k)

      integer ia,ib,irun,l,k,a(ia,l),b(ib,l)

        b(i,j) = a(i+offset-1,j)

        b(i,j) = a(i+offset-1,j)

        b(Lena+i,j) = a(Lena+offset-1,j)

      subroutine PASTIT(b,a,ib,ia,irun,L,k)

      _RL a(ia,l), b(ib,l)

      a(i+offset-1,j) = b(i,j)

      subroutine PSTBITINT(b,a,ib,ia,irun,l,k)

      _RL a(ia,l)

      a(i+offset-1,j) = a(i+offset-1,j) + float(b(i,j))

      subroutine PSTBMPIT(b,a,ib,ia,irun,l,k)

      _RL a(ia,l), b(ib,l)

      a(i+offset-1,j) = a(i+offset-1,j) + b(i,j)

      subroutine STRIP2TILE(a,indx,b,irun,ia,ib,levs,npeice)

      _RL a(ia,levs), b(ib,levs)

        b(i,k) = a(indx(i+offset-1),k)

        b(i,k) = a(indx(i+offset-1),k)

        b(Lena+i,k) = a(indx(Lena+offset-1),k)

      subroutine PASTE2GRD_OLD(b,indx,chfr,ib,numpts,a,ia,levs,npeice)

      _RL a(ia,levs), b(ib,levs), chfr(ib)

       a(indx(i+offset-1),L) = a(indx(i+offset-1),L) + b(i,L)*chfr(i)

      subroutine PASTE2GRD (b,indx,chfr,ib,numpts,a,ia,levs,npeice,

      _RL a(ia,levs), b(ib,levs), chfr(ib)

        if( a(indx(i+offset-1),L).eq.undef .or. b(i,L).eq.undef ) then

         a(indx(i+offset-1),L) = undef

         a(indx(i+offset-1),L)=a(indx(i+offset-1),L) + b(i,L)*chfr(i)

        a(indx(i+offset-1),L)=a(indx(i+offset-1),L) + b(i,L)*chfr(i)

      _RL a(lma)

      data a  / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,

         acr(i) = a(L+1)*temp + a(L)*(1-temp)

      if( pl(i).lt.p(1)   ) acr(i) = a(1)

      if( pl(i).gt.p(lma) ) acr(i) = a(lma)

      a  = getcon('EARTH RADIUS')

      dxu(i,j) = a*cos(phi(j))*dlam(i)

      dyv(i,j) = a*dphi(j)

      dyv(i,1)    = a*(dphi(1)   +0.5*dphi(2)   )

      dyv(i,jm-1) = a*(dphi(jm-1)+0.5*dphi(jm-2))

      _RL     pi,a,getcon,undef

      _RS function smoothMin_R4( a, b )

      _RS a, b

      smoothMin_R4 = .5*( a+b - smoothAbs_R4(a-b) )

      _RL function smoothMin_R8( a, b )

      _RL a, b

      smoothMin_R8 = .5*( a+b - smoothAbs_R8(a-b) )

          a(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k,bi,bj)*recip_drF(k)

          IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.

          b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)

          IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.) 

     &        bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )

     &        (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))

      _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)

         a(i,j,1) = 0. _d 0 

         a(i,j,1) = 0. _d 0 

         a(i,j,k) = -deltaTggl90

          IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.

          b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)

     I     a, b, c,

      _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)

     I                     a, b, c,

     &        (gXNm1(i,j,k) - a(i,j,k)*gYnm1(i,j,k-1))

      _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)

          IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.) 

     &        bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )

  These lines are here to deliberately cause a compile-time error.