c Computes modip latitude	for a given geographical grid
c It requires the file "dip.asc"
c
c Inputs:
c  la(*) = geographical longitude [0° to 360°] (r*8)
c  fi(*) = geographical latitude in degrees (r*8)  
c		 
c Output: 
c   fm(*) = modip in degrees (r*8) 
c
c Created by Claudio Brunini, 5/07.
c
      subroutine modip(la,fi,fm)
      implicit real*8 (a-h,o-z)
      real*8 la
	dimension z0(73,73)
	save z0
	pi=4.d0*datan(1.d0)
      d2r=pi/180.d0
c
c Reads tabulated	dip
c		
	if(z0(1,1).ge.0.d0)then
        open(10,file='dip.asc')
	  read(10,*)
	  read(10,*)
	  read(10,*)
	  y=-90.d0
	  i=1
	  do while (y.le.90.d0)
	    x=-180.d0
	    j=1
	    do while(x.le.180.d0)
	      read(10,'(26x,f8.2)') z0(i,j)
	      x=x+5.d0
	      j=j+1
	    enddo
	    y=y+2.5d0
	    i=i+1
	  enddo
   	  close(10)
	endif
c
c Calculates modip
c
      x=la
	if(x.gt.180.d0) x=x-360.d0
	y=fi
      kx=idint(x/5.d0+37.d0)
      ky=idint(y/2.5d0+37.d0)
      z11=z0(ky,kx )
      z12=z0(ky,kx+1)
      z21=z0(ky+1,kx)
      z22=z0(ky+1,kx+1)
      dx=x-dfloat(kx-37)*5.d0
      z1=dx*(z12-z11)/5.d0+z11
      z2=dx*(z22-z21)/5.d0+z21
      dy=y-dfloat(ky-37)*2.5d0 
      z=dy*(z2-z1)/2.5d0+z1
      fm=datan2(z*d2r,dsqrt(dcos(y*d2r)))/d2r

      return
      end