implicit real*8 (a-h,o-z)
      pi=4.d0*atan(1.d0)
      print*,'Input latitude and longitude and height in ETRS89/WGS84'
      read(5,*)elat,elong,ez
      elat=elat*pi/180
      elong=elong*pi/180
      call etrslg(elat,elong,x,y)
      call etrsos(x,y,ez)
      stop
      end
c
c
      subroutine etrslg(elat,elong,x,y)
      implicit real*8 (a-h,o-z)
      real*8 nu,nu1,n,m1,m2,m3,m4,m,i,ii,iii,iiia,iv,v,vi,north,east,
     1       north0,east0
      data a,b,f0,flat0,flong0,north0,east0
     1     /6378137.d0,6356752.3141d0,0.9996012717d0,
     1     49.d0,-2.d0,100000.d0,400000.d0/
      pi=4.d0*atan(1.d0)
      elat0=flat0*pi/180
      elong0=flong0*pi/180
      e2=1.d0-(b/a)**2
c     print*,'e2 = ',e2
      nu1=(1.d0-e2*sin(elat)**2)
      nu=a*f0/sqrt(nu1)
c     print*,'nu = ',nu
      rho=a*f0*(1.d0-e2)/nu1/sqrt(nu1)
c     print*,'rho = ',rho
      eta2=nu/rho-1.d0
c     print*,'eta2 = ',eta2
      n=(a-b)/(a+b)
      m1=(1.d0+n+1.25d0*n**2+1.25d0*n**3)*(elat-elat0)
      m2=(3*n+3*n**2+2.625*n**3)*sin(elat-elat0)*cos(elat+elat0)
      m3=1.875*n**2*(1.d0+n)*sin(2*(elat-elat0))*cos(2*(elat+elat0))
      m4=35.d0/24.d0*n**3*sin(3*(elat-elat0))*cos(3*(elat+elat0))
      m=b*f0*(m1-m2+m3-m4)
c     print*,'m = ',m
      i=m-north0
      ii=nu*sin(elat)*cos(elat)/2
      iii=nu*sin(elat)*cos(elat)**3*(5.d0-tan(elat)**2+9*eta2)/24
      iiia=nu*sin(elat)*cos(elat)**5*(61.d0-58.d0*tan(elat)**2+
     1                                tan(elat)**4)/720
      iv=nu*cos(elat)
      v=nu*cos(elat)**3*(nu/rho-tan(elat)**2)/6
      vi=nu*cos(elat)**5*(5.d0-18.d0*tan(elat)**2+tan(elat)**4+
     1                    14.d0*eta2-58.d0*tan(elat)**2*eta2)/120
      north=i+ii*(elong-elong0)**2+iii*(elong-elong0)**4+
     1      iiia*(elong-elong0)**6
      east=east0+iv*(elong-elong0)+v*(elong-elong0)**3+
     1     vi*(elong-elong0)**5
c     print*,'i = ',i
c     print*,'ii = ',ii
c     print*,'iii = ',iii
c     print*,'iiia = ',iiia
c     print*,'iv = ',iv
c     print*,'v = ',v
c     print*,'vi = ',vi
c     print*,'north = ',north
c     print*,'east = ',east
      x=east
      y=north
      return
      end
c
c
      subroutine etrsos(x,y,ez)
      implicit real*8 (a-h,o-z)
      external grdbox
      integer eindex,nindex,rec1,rec2,rec3,rec4,ieast,inorth
      real*8 east,north
      character box*2,grdbox*2
      eindex=int(x/1000)
      nindex=int(y/1000)
      x0=eindex*1000
      y0=nindex*1000
      rec1=eindex+nindex*701+1
      rec2=eindex+1+nindex*701+1
      rec3=eindex+1+(nindex+1)*701+1
      rec4=eindex+(nindex+1)*701+1
      open(unit=10,file='OSTN02_OSGM02_GB.txt',status='old')
      do i=1,rec1-1
        read(10,*)index
      enddo
      read(10,*)index,ie,in,se0,sn0,sg0
      read(10,*)index,ie,in,se1,sn1,sg1
      do i=1,699
        read(10,*)index
      enddo
      read(10,*)index,ie,in,se3,sn3,sg3
      read(10,*)index,ie,in,se2,sn2,sg2
      close(10)
c     print*,se0,sn0,sg0
c     print*,se1,sn1,sg1
c     print*,se2,sn2,sg2
c     print*,se3,sn3,sg3
      dx=x-x0
      dy=y-y0
c     print*,'dx = ',dx
c     print*,'dy = ',dy
      t=dx/1000
      u=dy/1000
c     print*,'t = ',t
c     print*,'u = ',u
      se=(1-t)*(1-u)*se0+t*(1-u)*se1+t*u*se2+(1-t)*u*se3
      sn=(1-t)*(1-u)*sn0+t*(1-u)*sn1+t*u*sn2+(1-t)*u*sn3
      sg=(1-t)*(1-u)*sg0+t*(1-u)*sg1+t*u*sg2+(1-t)*u*sg3
c     print*,'se = ',se
c     print*,'sn = ',sn
c     print*,'sg = ',sg
      east=x+se
      north=y+sn
      z=ez-sg
      print*,'National Grid/OSGB36 are:'
      print*,'East = ',east
      print*,'North = ',north
      print*,'Height = ' ,z
      ieast=int(east)
      inorth=int(north)
      box=grdbox(ieast,inorth)
      ieast=ieast-int(ieast/100000)*100000
      inorth=inorth-int(inorth/100000)*100000
      print*,box,ieast,inorth
      return
      end
c
c
      function grdbox(east,north)
      character grdbox*2,grid(7,13)*2
      integer east,north,ebox,norbox
      data grid/'SV','SW','SX','SY','SZ','TV','TW',
     1          'SQ','SR','SS','ST','SU','TQ','TR',
     2          'SL','SM','SN','SO','SP','TL','TM',
     3          'SF','SG','SH','SJ','SK','TF','TG',
     4          'SA','SB','SC','SD','SE','TA','TB',
     5          'NV','NW','NX','NY','NZ','OV','OW',
     6          'NQ','NR','NS','NT','NU','OQ','OR',
     7          'NL','NM','NN','NO','NP','OL','OM',
     8          'NF','NG','NH','NJ','NK','OF','OG',
     9          'NA','NB','NC','ND','NE','OA','OB',
     a          'HV','HW','HX','HY','HZ','JV','JW',
     b          'HQ','HR','HS','HT','HU','JQ','JR',
     c          'HL','HM','HN','HO','HP','JL','JM'/
      ebox=int(east/100000)+1
      norbox=int(north/100000)+1
      grdbox=grid(ebox,norbox)
      return
      end