**********************************************************************
*                                                                    *
*        The functions used above.                                   *
*                                                                    *
*       Version 0.5 - Last Modified 5/26/04                          *
*                                                                    *
**********************************************************************
**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*                                                                    *
*                                                                    *
**********************************************************************

        REAL function c1(i,j)

        if (af(i,j).ge.0.0) then
          c1=min(r1(i,j+1),r0(i,j))
         else
          c1=min(r1(i,j),r0(i,j+1))
          end if
        return
        end

**********************************************************************

        REAL function c2(i,j)

        if (ag(i,j).ge.0.0) then
          c2=min(r1(i+1,j),r0(i,j))
         else
          c2=min(r1(i,j),r0(i+1,j))
          end if
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*                                                                    *
*                                                                    *
**********************************************************************

        REAL function r1(i,j)

        INCLUDE "fpvar.inc"

        p(i,j)=max(0.0,ag(i-1,j))-min(0.0,ag(i,j))+
     $  max(0.0,af(i,j-1))-min(0.0,af(i,j))
        q(i,j)=(wmax(i,j)-dntd(i,j))*dxz
        if (p(i,j).gt.0.0) then
          r1=min(1.0,q(i,j)/p(i,j))
         else
          r1=0.0
          end if
        return
        end

**********************************************************************

        REAL function r0(i,j)

        INCLUDE "fpvar.inc"

        p(i,j)=max(0.0,ag(i,j))-min(0.0,ag(i-1,j))+
     $  max(0.0,af(i,j))-min(0.0,af(i,j-1))
        q(i,j)=(dntd(i,j)-wmin(i,j))*dxz
        if (p(i,j).gt.0.0) then
          r0=min(1.0,q(i,j)/p(i,j))
         else
          r0=0.0
          end if
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*       Set of Velocity*Density Calculation routines:                *
*                                                                    *
*       F is for the x-direction for a given cell,                   *
*       G is for the z-direction                                     *
*                                                                    *
**********************************************************************

        REAL function f(i,j)

        INCLUDE "fpvar.inc"

        vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*sini/((1.0+rvi(i)**2)*b)
     $  +sini*cosi*gg/((1.0+rvi(i)*rvi(i))*vin(i))
     $  -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
        f=vix1*n(i,j)
        return
        end

**********************************************************************

        REAL function g(i,j)

        INCLUDE "fpvar.inc"

        viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b))
     $  +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b))
     $  -eoy*cosi/((1.0+rvi(i)**2)*b)
     $  -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i))
     $  +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b)
     $  -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b)
        g=viz1*n(i,j)
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*                                                                    *
*                                                                    *
**********************************************************************
 
       REAL function FL(I,J)

        INCLUDE "fpvar.inc"

        fl=(0.5*dt*(f(i,j+1)+f(i,j))-0.25*dx*(n(i,j+1)-n(i,j)))*dz
        return
        end

**********************************************************************

        REAL function GL(I,J)

        INCLUDE "fpvar.inc"

        gl=(0.5*dt*(G(i+1,j)+G(i,j))-0.25*dz*(N(i+1,j)-n(i,j)))*dx
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*                                                                    *
*                                                                    *
**********************************************************************

        REAL FUNCTION fh(I,J)

        INCLUDE "fpvar.inc"

        fh=(7.0*(F(i,j+1)+F(i,j))-(F(i,j+2)+F(i,j-1)))*DT*dz/12.0
        RETURN
        END

**********************************************************************

        REAL function gh(i,j)

        INCLUDE "fpvar.inc"

        gh=(7.0*(g(i+1,j)+g(i,j))-(g(i+2,j)+g(i-1,j)))*dt*dx/12.0
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*                                                                    *
*                                                                    *
**********************************************************************

        REAL function ag(i,j)

        INCLUDE "fpvar.inc"

        ag=gh(i,j)-gl(i,j)
        t1=ag*(dntd(i+1,j)-dntd(i,j))
        t2=ag*(dntd(i+2,j)-dntd(i+1,j))
        t3=ag*(dntd(i,j)-dntd(i-1,j))
        if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) ag=0.0
        return
        end

**********************************************************************

        REAL function af(i,j)

        INCLUDE "fpvar.inc"

        AF=FH(I,J)-FL(I,J)
        t1=af*(dntd(i,j+1)-dntd(i,j))
        t2=af*(dntd(i,j+2)-dntd(i,j+1))
        t3=af*(dntd(i,j)-dntd(i,j-1))
        if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) af=0.0
        return
        end

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*       Set of max/min functions designed to retrieve max/min value  *
*       of n or dntd  for (i,j) and neighbouring cells.              *
*                                                                    *
*       Used in functions r0 and r1                                  *
*                                                                    *
**********************************************************************

        REAL FUNCTION WMAX(I,J)

        INCLUDE "fpvar.inc"

        WMAX=MAX(N(I-1,J),N(I,J),N(I,J-1),N(I,J+1),N(I+1,J),
     $   DNTD(I-1,J),DNTD(I,J),DNTD(I,J-1),DNTD(I,J+1),DNTD(I+1,J))
        RETURN
        END

**********************************************************************

        REAL FUNCTION WMIN(I,J)

        INCLUDE "fpvar.inc"

        WMIN=MIN(N(I-1,J),N(I,J),N(I,J-1),N(I,J+1),N(I+1,J),
     $   DNTD(I-1,J),DNTD(I,J),DNTD(I,J-1),DNTD(I,J+1),DNTD(I+1,J))
        RETURN
        END

**********************************************************************
*////////////////////////////////////////////////////////////////////*
**********************************************************************
*                                                                    *
*       Set of functions to change between degrees and radians.      *
*                                                                    *
**********************************************************************

        REAL FUNCTION TORAD(X)

        REAL X
        PARAMETER (PI=3.141592654)
        TORAD=X*PI/180.
        RETURN
        END

**********************************************************************

        REAL FUNCTION TODEG(X)

        REAL X
        PARAMETER (PI=3.141592654)
        TODEG=X*180./PI
        RETURN
        END