gnorm-gh.f

c     calculate normal gravity 1980
      implicit real*8 (w)
      parameter (wrad = 3.14159265359d0 / 180.d0)
      parameter (wc0 = 978032.67715d0, wc2 = 5163.07d0, wc4 = 22.76d0)
c
      dimension nn(30000),kyr(30000),kdy(30000),mmin(30000)
      dimension kdp(30000),kmgt(30000),kmga(30000),ggr(30000)
      dimension ggto(30000),gga(30000),kflg(30000)
      dimension hdg(30000),spd(30000)
      dimension wwlat(30000),wwlon(30000)
      dimension wolat(30000),wolon(30000)
c
      open(unit=10,file='H:\Myhome\gh90a-wgs.txt',
     +   status='old')
      open(unit=20,file='H:\Myhome\gh90a.datas1',
     +   status='old')
      open(unit=30,file='H:\Myhome\gh90a.bs1')
c      
      i=0
    1 continue
          read(20,*,end=99) no,iyr,jdy,min,alt,aln,
     +   spd0,hdg0,idp,mgt,mga,gr,gto,ga,iflg
         if(no.eq.0) then
           go to 1
         else
           i=i+1
            nn(i)=no
            kyr(i)=iyr
            kdy(i)=jdy
            mmin(i)=min
            kdp(i)=idp
            kmgt(i)=mgt
            kmga(i)=mga
            ggr(i)=gr
            ggto(i)=gto
            gga(i)=ga
            kflg(i)=iflg
           spd(i)=spd0
           hdg(i)=hdg0
           wolat(i)=dble(alt)
           wolon(i)=dble(aln)
         end if
         go to 1
   99 close(20)
      jnum=i
      write(6,*) 'read end-20'
      k=0
   11 continue
      read(10,*,end=9) wgslat,wgslon
      k=k+1
      wwlat(k)=dble(wgslat)
      wwlon(k)=dble(wgslon)
      go to 11
    9 close(10)
      knum=k
      write(6,*) jnum,knum
      if(jnum.ne.knum) then
         write(6,*) 'caution!!'
         go to 89
      end if
c
      do 46 ll=1,knum
       head=hdg(ll)
       speed=spd(ll)
c      wgs=>1, Tokyo=>2
            wlatw=wwlat(ll)
            jj=1
            call etvoes(jj,wlatw,speed,head,etvn,wgnrmn)
c
            wlato=wolat(ll)
            jj=2
            call etvoes(jj,wlato,speed,head,etvo,wgnrmo)
c
      if (ggr(ll).ge.90000) then
          ggo=999999.9
          ganew=9999.9
       else
          ggo=ggto(ll)-etvo+etvn
          ganew=ggo-sngl(wgnrmn)+0.87
      end if
c
      alt=sngl(wwlat(ll))
      aln=sngl(wwlon(ll))
      write(30,207) nn(ll),kyr(ll),kdy(ll),mmin(ll),alt,aln,
     +   spd(ll),hdg(ll),kdp(ll),kmgt(ll),kmga(ll),ggr(ll),ggo,ganew,
     +   kflg(ll)
   46 continue
   89 continue
      close(30)
  207 format(i7,i3,i4,i5,f10.5,f10.5,f6.2,f6.1,i6,i6,i6,f8.1,
     +   f10.1,f8.1,i2)
c
      stop
      end
c
      subroutine etvoes(jj,wlat,speed,head,etvn,wgnorm)
      implicit real*8 (w)
      parameter (wrad = 3.14159265359d0 / 180.d0)
      parameter (wc0 = 978032.67715d0, wc2 = 5163.07d0, wc4 = 22.76d0)
c
C     wlat is degrees
      wslat = dsin(wlat*wrad)
      wslat2 = wslat*wslat
      wslat4 = wslat2*wslat2
c
      wgnorm = wc0 + wc2*wslat2 + wc4*wslat4
c
      wclat = dcos(wlat*wrad)
      wshd=dsin(head*wrad)
      vv=speed
      if(jj.eq.1) then
       etvn=(7.3*wclat*wshd+ 0.00392*vv)*vv
      else
       etvn=(7.5*wclat*wshd+0.0042*vv)*vv
      end if
      return
      end