      program proc

      open(unit=1,file='cal_stage.dat',form='formatted',status='old')

  1   read(1,*,end=2)itile,iseas,xlon,xlat,stg0,stg5,stg4,stg3,stgtot,
     .               s0,s5,s4,s3
      read(1,*)stot,nobs,avgdp,tarea

      call mercgom2(1,xlon,xlat,xpos,ypos)

      write(10+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               stg0,stg0,stg0,.1
      write(20+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               stg5,stg5,stg5,.1
      write(30+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               stg4,stg4,stg4,.1
      write(40+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               stg3,stg3,stg3,.1
      write(50+iseas,105)xpos,ypos,10.,float(30+(iseas-1)*60),
     .               stgtot,stgtot,stgtot,.1

      goto 1

  2   continue

 105  format(12e13.5)

c
c output format x,y,z,year day,t,s,st,.1
c

      end

C=======================================================================
C========= ROUTINES OF MERCATOR PROJECTION =============================
C=======================================================================
C SUBROUTINE MERCGOM1
C Provided to NML by D. Greenberg
C Latest Revision by Christopher E. Naimie
C Dartmouth College - 25 NOV 1991
C
C This set of routines was sent to Dartmouth by D. Greenberg.
C   In its present form it performs a mercator projection
C   which is valid for the gom1,gb1, and bf1 meshes to within
C   0.0025 deg (5th sig fig.) for longitude and latitude measurements
C   and to within 200 meters(4th sig fig.) for x and y measurements.
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
C Units: XLON is the longitude in degrees, positive east of Grenwich.
C        XLAT is the latitude in degrees, positive north of the Equator.
C        the x and y coordinates (XPOS,YPOS) are in meters.
C-----------------------------------------------------------------------

      SUBROUTINE MERCGOM1(NOPT,XLON,XLAT,XPOS,YPOS)
C
C      scale is specified at the equator
C      R = degrees/radian   S = cm per degree of latitude
C
      DATA A,B/-1.57079632,0.0067226/
      DATA R,S/57.2957795,1.111338e7/
C Set constants
      SCALE=1353100.3
      xoffset=17.1108
      yoffset=28.5281
C
      if(kkk.NE.64)then
         SF=SCALE/S
         RLAT=42.35
         RLON=-71.0333
      kkk=64
      endif
C
      IF(NOPT.EQ.1)then
         YPOS=((VB(XLAT/R)-VB(RLAT/R))*R/SF+yoffset)*10000.0
         XPOS=((XLON-RLON)/SF+xoffset)*10000.0
         RETURN
C
      else IF(NOPT.EQ.2)then
         XLON=((XPOS/10000.0-xoffset)*SF)+RLON
         ABA=((YPOS/10000.0-yoffset)*SF/R)+VB(RLAT/R)
         ABB=A+2.0*ATAN(EXP(ABA))
         ABC=A+2.0*ATAN(EXP(ABA+B*SIN(ABB)))
         XLAT=R*ABC
         RETURN
      else
         PRINT*,'INVALID MERCATER OPTION'
      endif
      RETURN
      END
C=======================================================================
      FUNCTION VB(X)
C *
C * INCREASING LATITUDE IN A MERCATOR PROJECTION
C * X IN RADIANS, VB IN RADIANS
C
      VB=ALOG(TAN(0.78539816+.5*X))-0.0067226*SIN(X)
      return
      END
C=======================================================================
C*****67****************************************************************
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***********************************************************************
