!
      subroutine othend
!
! 6/30/98 btf:
!
! For each geographic grid point, determine the geographic coordinates
! of the point at the opposite end of the magnetic field line. This
! routine is called once per model run from subroutine start.
! Results are placed in /othrend_com/ in othend.h.
!
! This subroutine is based on a program written by Cicely Ridley
! (see ~ecridley/othend or ~foster/dyn/othend.f). This routine reads
! magnetic coordinates and magnetic field data from the mss file 
! /FOSTER/magdat/mg8. This mss file is equivalent to the now deleted 
! file /ECRIDLEY/ECR90/ECRMG8, and was written by ~foster/dyn/dyn16.f, 
! based on ~ecridley/tgcm/dyn16.Z.
!
      include "params.h"
      PARAMETER (IMAXG=ZIMX, JMAXG=ZJMX)
      PARAMETER (IMAXGP=IMAXG+1, JMAXGP=JMAXG+1)
      include "othend.h"   ! output rlamda, rtheta
      include "fieldz.h"   ! /FIELD/ read from mss_magfile
      include "cterp.h"    ! /TERP/  read from mss_magfile
      include "trig.h"     ! /TRIG/  read from mss_magfile
      include "strt.h"     ! for tmpdir
!
      character(len=80) :: mss_magfile = '/FOSTER/magdat/mg8'
      character(len=8 ) :: diskfile = 'mg8.dat '
      integer :: lu,ios,len,i,j
      real :: err,rtd
      real :: XXG(IMAXGP,0:JMAXGP,2),YYG(IMAXGP,0:JMAXGP,2),
     +  ZZG(IMAXGP,0:JMAXGP,2),XXM(IMAXMP,JMAXM),YYM(IMAXMP,JMAXM),
     +  ZZM(IMAXMP,JMAXM)
      real,external :: unit
!
! Get mag data file from mss:
!
      call getms(mss_magfile,diskfile,tmpdir,' ')
      lu = nextlu()
      open(lu,file=diskfile,status='OLD',err=900,iostat=ios,
     +    form='unformatted')
!
! Read mag data into /FIELD/ (fieldz.h), /TERP/ (cterp.h), and
! /TRIG/ (trig.h):
!
      buffer in(lu,1) (alatm(1,0),rjacd(imaxgp,jmaxgp)) ! /FIELD/
      err = unit(lu)
      if (err /= -1.) then
        len = length(lu) 
        write(6,"(/'>>> othend: Error reading /FIELD/: err=',f9.2,
     +    ' len=',i8)") err,len
        stop 'othend' 
      endif
!
      buffer in(lu,1) (ig(1,1),djm(imaxgp,jmaxgp)) ! /TERP/
      err = unit(lu)
      if (err /= -1.) then
        len = length(lu) 
        write(6,"(/'>>> othend: Error reading /TERP/: err=',f9.2,
     +    ' len=',i8)") err,len
        stop 'othend' 
      endif
!
      buffer in(lu,1) (cslatm(1,1),snlong(imaxgp)) ! /TRIG/
      err = unit(lu)
      if (err /= -1.) then
        len = length(lu) 
        write(6,"(/'>>> othend: Error reading /TRIG/: err=',f9.2,
     +    ' len=',i8)") err,len
        stop 'othend' 
      endif
      close(lu)
C     ****
C     ****     SET UP ARRAYS OF GEOGRAPHIC COORDINATES, XXG, YYG, ZZG
C     ****
      DO 2 J = 1,JMAXG
	DO 2 I = 1,IMAXGP
	  XXG(I,J,1) = CSLATG(J)*SNLONG(I)
	  YYG(I,J,1) = CSLATG(J)*CSLONG(I)
	  ZZG(I,J,1) = SNLATG(J)
    2 CONTINUE
C     ****
C     ****     OPLAR POINTS
C     ****
      DO 3 I = 1,IMAXGP
	XXG(I,0,1) = 0.
	XXG(I,JMAXGP,1) = 0.
	YYG(I,0,1) = 0.
	YYG(I,JMAXGP,1) = 0.
	ZZG(I,0,1) = -1.
	ZZG(I,JMAXGP,1) = 1.
    3 CONTINUE
C     ****
C     ****     TRANSFORM XXG, YYG, ZZG TO GEOMAGNETIC COORDINATES
C     ****
      DO 4 J = 1,JMAXM
	CALL GRDINT(XXM(1,J),XXG(1,0,1),IG,JG,WT,IMAXGP,IMAXMP,IMAXG,
     1    JMAXG+2,IMAXM,JMAXM,J)
	CALL GRDINT(YYM(1,J),YYG(1,0,1),IG,JG,WT,IMAXGP,IMAXMP,IMAXG,
     1    JMAXG+2,IMAXM,JMAXM,J)
	CALL GRDINT(ZZM(1,J),ZZG(1,0,1),IG,JG,WT,IMAXGP,IMAXMP,IMAXG,
     1    JMAXG+2,IMAXM,JMAXM,J)
    4 CONTINUE
C     ****
C     ****     PERIODIC POINTS
C     ****
      DO 5 J = 1,JMAXM
	XXM(IMAXMP,J) = XXM(1,J)
	YYM(IMAXMP,J) = YYM(1,J)
	ZZM(IMAXMP,J) = ZZM(1,J)
    5 CONTINUE
C     ****
C     ****     NOW REVERSE TRANSFORMATION
C     ****
      CALL GRDTRP(XXM(1,1),XXG(1,0,2),IM(1,0),JM(1,0),DIM(1,0),DJM(1,0),
     1  IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM)
      CALL GRDTRP(YYM(1,1),YYG(1,0,2),IM(1,0),JM(1,0),DIM(1,0),DJM(1,0),
     1  IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM)
      CALL GRDTRP(ZZM(1,1),ZZG(1,0,2),IM(1,0),JM(1,0),DIM(1,0),DJM(1,0),
     1  IMAXGP,IMAXMP,IMAXG,JMAXG+2,IMAXM,JMAXM)
C     ****
C     ****     PERIODIC POINTS
C     ****
      DO 7 J = 0,JMAXGP
	XXG(IMAXGP,J,2) = XXG(1,J,2)
	YYG(IMAXGP,J,2) = YYG(1,J,2)
	ZZG(IMAXGP,J,2) = ZZG(1,J,2)
    7 CONTINUE
C     ****
C     ****     CONVERT (XXG,YYG,ZZG,2) TO (LAMDA,THETA) (degrees)
C     ****
      RTD = 180./(4.*ATAN(1.))
      DO 6 J = 1,JMAXG
	DO 6 I = 1,IMAXGP
	  RLAMDA(I,J) = ATAN2(XXG(I,J,2),YYG(I,J,2))*rtd
	  RTHETA(I,J) = ATAN2(ZZG(I,J,2),SQRT(XXG(I,J,2)**2+
     1      YYG(I,J,2)**2))*rtd
    6 CONTINUE
!
! Write to stdout:
!     do j=1,jmaxg
!       write(6,"(/'othend: j=',i2,' (jmaxg=',i2,') rlamda(*,j)=',/
!    +    (6e12.4))") j,jmaxg,rlamda(:,j)
!       write(6,"( 'othend: j=',i2,' (jmaxg=',i2,') rtheta(*,j)=',/
!    +    (6e12.4))") j,jmaxg,rtheta(:,j)
!     enddo
!
! Contour results:
!     call mkcon(rlamda,imaxgp,jmaxg,0.,0,'longitude','latitude',
!    +  'RLAMDA from othend')
!     call mkcon(rtheta,imaxgp,jmaxg,0.,0,'longitude','latitude',
!    +  'RTHETA from othend')
      return
!
! Problem opening disk file:
!
 900  continue
      write(6,"('>>> othend: error opening disk file ',a,
     +  ' ios=',i4)") diskfile(1:len_trim(diskfile)),ios
      write(6,"('    mss path = ',a)") 
     +  mss_magfile(1:len_trim(mss_magfile))
      stop 'othend'
      end