      program extract_ann_signal
       
       implicit real*8 (a-h,o-z)
       integer*4 ind1, ind2
       character*56 fname,oname
       real*4 dsin(5000,5000),dcos(5000,5000),damp(5000,5000)
       real*4 dpha(5000,5000)
       topi = 2.*3.1415926536d0

c      Main
       write(*,*) ' Enter decimal year (eg. 0.87/1 or 0.1/1)'
       read(*,*) dec
       time_rad =  topi*dec
       write(*,*) ' Enter latitude (-90 to 90)'
       read(*,*) glat      
       write(*,*) ' Enter longitude (0 to 360) '
       read(*,*) glon     
   
       open (20,file = 'MONARCH-A_ANN_cos.grd',form = 'formatted')
       open (21,file = 'MONARCH-A_ANN_sin.grd',form = 'formatted')
       open (22,file = 'MONARCH-A_ANN_amp.grd',form = 'formatted')
       open (23,file = 'MONARCH-A_ANN_pha.grd',form = 'formatted')
       read(20,*) alat1,alat2,alon1,alon2,dlat,dlon
       read(21,*) alat1,alat2,alon1,alon2,dlat,dlon
       read(22,*) alat1,alat2,alon1,alon2,dlat,dlon
       read(23,*) alat1,alat2,alon1,alon2,dlat,dlon
       nlon = (alon2-alon1)/dlon + 1
       nlat = (alat2-alat1)/dlat + 1
c      read arrays
       do i = 1,nlat	 
 	 read(20,*) (dcos(i,j),j=1,nlon)
	 read(21,*) (dsin(i,j),j=1,nlon)	
	 read(22,*) (damp(i,j),j=1,nlon)
	 read(23,*) (dpha(i,j),j=1,nlon)
       end do
       
c      bilinear interpolation in points
       x1 = int(glon)
       x2 = int(glon)+1
       if (x1 .eq. 360) x1=0
       if (x2 .eq. 361) x2=1
       y1 = int(glat) 
       y2 = int(glat) + 1
       ind1 = 91
       ind2 = 90
       if(glat < 0 ) then
       y1 = int(glat) - 1
       y2 = int(glat)          
       end if
       yfac1 = glat-y1
       yfac2 = y2-glat
  
c      sin
       q_x1y1=  dsin(ind1-y1,x1+1)*(x2-glon)*yfac2
       q_x2y1=  dsin(ind1-y1,x1+2)*(glon-x1)*yfac2
       q_x1y2=  dsin(ind2-y1,x1+1)*(x2-glon)*yfac1      
       q_x2y2=  dsin(ind2-y1,x1+2)*(glon-x1)*yfac1
       f_xy_sin = q_x1y1+q_x1y2+q_x2y1+q_x2y2
c      cos       
       q_x1y1=  dcos(ind1-y1,x1+1)*(x2-glon)*yfac2
       q_x2y1=  dcos(ind1-y1,x1+2)*(glon-x1)*yfac2
       q_x1y2=  dcos(ind2-y1,x1+1)*(x2-glon)*yfac1      
       q_x2y2=  dcos(ind2-y1,x1+2)*(glon-x1)*yfac1
       f_xy_cos = q_x1y1+q_x1y2+q_x2y1+q_x2y2
c      amp
       q_x1y1=  damp(ind1-y1,x1+1)*(x2-glon)*yfac2
       q_x2y1=  damp(ind1-y1,x1+2)*(glon-x1)*yfac2
       q_x1y2=  damp(ind2-y1,x1+1)*(x2-glon)*yfac1     
       q_x2y2=  damp(ind2-y1,x1+2)*(glon-x1)*yfac1
       f_xy_amp = q_x1y1+q_x1y2+q_x2y1+q_x2y2
      
c      annual signal and phase      
       f_xy_pha = 360*ATAN2(f_xy_sin,f_xy_cos)/topi
       if (f_xy_pha .ge. 180) f_xy_pha = f_xy_pha -360 
       dhann = f_xy_amp*cos(time_rad-f_xy_pha*topi/360)
       
c      write to screen    
       write(*,*) 'Sin= ', f_xy_sin,' Cos= ', f_xy_cos	
       write(*,*) 'Lat= ', glat,' Lon= ', glon, ' dhann(t)= ', dhann
       write(*,*) 'Amp(lat,lon)= ', f_xy_amp 
       write(*,*) 'Phase(lat,lon)= ', f_xy_pha
       
       end
