C***********************************************************************
C***********************************************************************
      SUBROUTINE POIVAL3D(nndim,nnvdim,nedim,in,jjj,nelem,basis,data,
     &atpoint)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the real nodal based 
C            quantity stored in the array data, on element nelem at 
C            vertical node jjj.
C
C restrictions: Applicable only 3-D data which has been computed using 
C            linear triangular elements in the horizontal and linear 
C            elements in the vertical (where the nodes line up in the 
C            vertical). 
C
C inputs:  in(3,nedim) is the element incidence list
C            NOTE:  The incidence list must be entered in a
C                     counterclockwise order (LABLIB standard)
C          nedim is the dimensioning parameter for the incidence list
C          nndim is the horizontal dimensioning parameter for data
C          nnvdim is the vertical dimensioning parameter for data
C          data is the quantity to be interpolated
C          nelem is the element number on which the desired point lies 
C          basis(3) are the values of the basis functions at the point
C            where the nodal based quantity is desired
C
C outputs: atpoint is the real value of data at the desired point
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      integer in(3,nedim)     
      real data(nndim,nnvdim),basis(3)
C
C START OF EXECUTABLE ROUTINE
      atpoint=data(in(1,nelem),jjj)*basis(1)
     &+data(in(2,nelem),jjj)*basis(2)+data(in(3,nelem),jjj)*basis(3)
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE BASIS2D(nndim,nedim,in,x,y,nn,ne,xval,yval,
     &nclosest,nelem,basis)
C-----------------------------------------------------------------------
C purpose: This subroutine determines the element on which the point
C            (xval,yval) resides and determines the value of the basis
C            functions at the point.  In the first pass, only elements 
C            containing the closest node to (xval,yval) are tested to 
C            see if (xval,yval) lies on them.  If this fails, all
C            elements are tested to see if (xval,yval) lies on them.
C
C
C restrictions: Applicable only for  2-D triangular linear elements
C
C calling sequence: CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C                   CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
C                   CALL PHI2D(xlocal,ylocal,xval,yval,basis)
C
C inputs:  nn is the number of nodes in the mesh
C          x(nn),y(nn) 1-D arrays containing the are the x and y
C            coordinates of each nodes
C          ne is the number of elements in the mesh
C          in(3,ne) is the LABLIB Standard incidence list
C          xval,yval is the location of the point 
C          nedim is the size of the arrays in the calling routine
C
C outputs: basis is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C          nelem is the element on which (xval,yval),resides
C          nclosest is the closest node to (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      integer in(3,nedim)
      real x(nndim),y(nndim),basis(3)
C
C LOCAL VARIABLES
      dimension xlocal(3),ylocal(3)
C
C START OF EXECUTABLE CODE
C
C  Find the closest node to (xval,yval)
      CALL CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C
C  Find if an element containing the closest node also contains the
C    point (xval,yval)
C
      do 10 i=1,ne
         n1=in(1,i)
         n2=in(2,i)
         n3=in(3,i)
         if(nclosest.EQ.n1.OR.nclosest.EQ.n2.OR.nclosest.EQ.n3)then
            xlocal(1)=x(n1)
            ylocal(1)=y(n1)
            xlocal(2)=x(n2)
            ylocal(2)=y(n2)
            xlocal(3)=x(n3)
            ylocal(3)=y(n3)
            CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
            if(nflag.EQ.1) then
               nelem=i
               CALL PHI2D(xlocal,ylocal,xval,yval,basis)
               return
            endif
         endif
 10   continue
C
C  Find if any element contains the point (xval,yval)
C
      do 20 i=1,ne
         xlocal(1)=x(in(1,i))
         ylocal(1)=y(in(1,i))
         xlocal(2)=x(in(2,i))
         ylocal(2)=y(in(2,i))
         xlocal(3)=x(in(3,i))
         ylocal(3)=y(in(3,i))
         CALL POINTONEL(xlocal,ylocal,xval,yval,nflag)
         if(nflag.EQ.1) then
            nelem=i
            CALL PHI2D(xlocal,ylocal,xval,yval,basis)
            return
         endif
 20   continue
      nelem=0
C
C END OF ROUTINE
C
      return
      end
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE CLOSESTNODE(nndim,x,y,nn,xval,yval,nclosest)
C-----------------------------------------------------------------------
C purpose: This subroutine finds the node closest to a point specified
C            by the coordinates (xval,yval)
C  
C restrictions:  all nodes must lie in a flat plane
C
C inputs:  x(nn) and y(nn) are 1-D arrays containing the x and y
C            coordinates of each node
C          nn is the number of nodes
C          xval,yval are the point coordinates
C
C outputs: nclosest is the node who's (x,y) coordinates are closest to
C            (xval,yval)
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension x(nndim),y(nndim)
C
C START OF EXECUTABLE CODE
C
C  Find node closest to (xval,yval)
      dclosest=SQRT((x(1)-xval)**2.0+(y(1)-yval)**2.0)
      nclosest=1
      do 25 j=2,nn
         dval=SQRT((x(j)-xval)**2.0+(y(j)-yval)**2.0)
         if(dval.LT.dclosest)then
            dclosest=dval
            nclosest=j
         endif
 25   continue
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE POINTONEL(xlocal,ylocal,xval,yval,nflag)
C-----------------------------------------------------------------------
C purpose: This Subroutine determines if a point (xval,yval) is on a 
C            local 2-D triangular linear element
C    
C restrictions: Applicable only for  2-D triangular linear element
C
C inputs:  xlocal(3) and y(3) are 1-D arrays containing the x and y 
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xval,yval is the location of the point 
C
C outputs: nflag is set equal to 1 if the point (xval,yval) is on the 
C            local element.
C          nflag is set equal to 0 if the point (xval,yval) is not on 
C            the local element.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xlocal(3),ylocal(3)
C
C LOCAL VARIABLES
      dimension vx(3),vy(3)

      real epsilon 
      parameter (epsilon = 1e-6)
C
C  Calculate the x and y components of vectors pointing from (xval,yval)
C    to each node on the element
      do 10 i=1,3
         delx=xlocal(i)-xval
         dely=ylocal(i)-yval
         d=SQRT(delx**2.0+dely**2.0)
	 if((abs(delx).lt.epsilon).and.(abs(dely).lt.epsilon)) then 
	    theta = 0.0 
	 else
            theta=ATAN2(dely,delx)
	 endif
         vx(i)=d*COS(theta)
         vy(i)=d*SIN(theta)
 10   continue
C
C  Determine if the point is on the element by calculating the
C    crossproducts of neighboring vectors in a direction
C    which yeild all positive numbers if x,y is on the element.
C    (ie: v1Xv2>0, v2Xv3>0, and v3Xv1>0 => x,y is on element)
      do 20 i=1,3
         j=i+1
         if(i.EQ.3)j=1
         crossprod=vx(i)*vy(j)-vy(i)*vx(j)
         if(crossprod.LT.0.0)then
            nflag=0
            go to 21
         endif
 20   continue
      nflag=1
 21   continue
C
C END OF ROUTINE
      return
      end
C****67***************************************************************72
C***********************************************************************
      SUBROUTINE PHI2D(xlocal,ylocal,xval,yval,phi)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the basis functions 
C            alive on a 2-D triangular linear element at a 
C            point (xval,yval) on a local element.
C
C restrictions: Applicable only for  2-D triangular linear element
C
C inputs:  xlocal(3) and ylocal(3) are 1-D arrays containing the x and y
C            coordinates of each local node on this element.
C            NOTE:  The node locations must be entered in a
C                     counterclockwise order (LABLIB standard)
C          xval,yval is the location of the point 
C
C outputs: phi is a 1-D array containing the value of the 3 basis
C            functions at the point xval,yval.
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           22 APRIL 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension xlocal(3),ylocal(3),phi(3)
C
C START OF EXECUTABLE CODE
C
C Calculate the area of the element
      area=0.5*(xlocal(1)*(ylocal(2)-ylocal(3))
     &+xlocal(2)*(ylocal(3)-ylocal(1))+xlocal(3)*(ylocal(1)-ylocal(2)))
C
C Calculate the value of the basis functions at xval,yval
      do 20 i=1,3
         j=i+1
         k=i+2
         if(i.EQ.2)k=1
         if(i.EQ.3)then
            j=1
            k=2
         endif
         a=(xlocal(j)*ylocal(k)-xlocal(k)*ylocal(j))/(2.0*area)
         b=(ylocal(j)-ylocal(k))/(2.0*area)
         c=-(xlocal(j)-xlocal(k))/(2.0*area)
         phi(i)=a+b*xval+c*yval
 20   continue
C
C END OF ROUTINE
      return
      end
C***********************************************************************
C***********************************************************************
C SUBROUTINE mercgom2
C Christopher E. Naimie
C Dartmouth College
C Latest version - 3 MAY 1993
C
C This subroutine performs a mercator projection, using 
C   Boston, MA., USA as the (x,y) origin.
C 
C Units: dlong is the longitude in degrees, positive east of Grenwich.
C        dlat is the latitude in degrees, positive north of the Equator.
C        x and y are in MKS units.
C
C The constants used are:
C      longitude at Boston = -71.03 deg E
C      latitude at Boston  =  42.35 deg N
C      radius of the earth =  6.3675E+6 meters 
C
C The projection type is set by nopt:
C   nopt=1  =>  convert provided (dlong,dlat) to (x,y)
C   nopt=2  =>  convert provided (x,y) to (dlong,dlat)
C-----------------------------------------------------------------------
      SUBROUTINE mercgom2(nopt,dlong,dlat,x,y)
C Set constants and (xo,yo) coordinates of Boston
      r=6.3675E+6                            
      reflong=rad(-71.03)                    
      reflat=rad(42.35)                      
      xo=r*COS(reflat)*reflong
      yo=r*COS(reflat)*ALOG((1.0+SIN(reflat))/COS(reflat))
C If nopt=1:  Convert (long,lat) to (x,y) and subtract (xo,yo)
      if(nopt.EQ.1)then
         rlong=rad(dlong)              
         rlat=rad(dlat)
         x=r*COS(reflat)*rlong-xo
         y=r*COS(reflat)*ALOG((1.0+SIN(rlat))/COS(rlat))-yo             
         return
C If nopt=2:  Convert (x+xo,y+xo) to (long,lat)
      else if(nopt.EQ.2)then
         rlong=(x+xo)/(r*COS(reflat))
         call newtlat(rlat,y+yo,reflat,r)
         dlong=deg(rlong)
         dlat=deg(rlat)
         return
      else
         write(*,*)'SUBROUTINE mercbos requires nopt be set to 1 or 2.'
         stop
      endif
      end      
C*****67****************************************************************
C***********************************************************************
      SUBROUTINE NEWTLAT(rlat,y,reflat,r)
C-----------------------------------------------------------------------
C Christopher E. Naimie
C Dartmouth College - 25 NOV 1991
C
C This subroutine uses the Newton-Raphson Method to calculate the 
C   latitude (in radians), given the following information:
C      r=radius of the earth (m)
C      y=y coordinate of the point who's latitude is desired (m/sec)
C      reflat=reference latitude where y is zero (radians)
C The latitude at y is returned as rlat (in radians).
C-----------------------------------------------------------------------
C Set constants
      tol=0.0001
      n=1000
C  Enter iteration loop
      po=reflat
      do 10 i=1,n
         p1=po-((1.0+SIN(po))/COS(po)-EXP(y/(r*COS(reflat))))
     &         /((1.0+SIN(po))/((COS(po))**2.0))
         if(ABS(p1-po).lt.tol)then
           rlat=p1
           return
         endif
         po=p1
 10   continue
C
C  If loop has completed, the method failed.
C
      write(*,*)' A root was not found after ',n,' iterations.'
      return
      end
C****67****************************************************************
C****67****************************************************************
C Convert degrees to radians
      real function rad(degrees)
      pi=3.14159265359
      rad=degrees*pi/180.0
      return
      end
C*****67****************************************************************
C***********************************************************************
C Convert radians to degrees
      real function deg(radians)
      pi=3.14159265359
      deg=radians*180.0/pi
      return
      end
C*****67****************************************************************
C***********************************************************************
