#include "dims.h"
      subroutine sunloc
      use cons_module,only: jmax,dlamda,dphi,pi
      use init_module,only: secs,iday
      implicit none
#include "params.h"
      real :: rlatm,rlonm,dipmag,decmag,sndec,csdec,sn2dec,sncsdc,
     |  rlatmp,rlonmp,dlons
      COMMON/TRGM/RLATM(ZIMXP,ZJMX),RLONM(ZIMXP,ZJMX),
     1DIPMAG(ZIMXP,ZJMX),DECMAG(ZIMXP,ZJMX),SNDEC(ZIMXP,ZJMX),
     2CSDEC(ZIMXP,ZJMX),SN2DEC(ZIMXP,ZJMX),SNCSDC(ZIMXP,ZJMX)
     3,RLATMP(ZJMX),RLONMP(ZJMX),DLONS(ZJMX)
C     ****
C     ****     CALCULATE SUN'S LONGITUDE, DLONS, IN DIPOLE COORDINATES.
C     ****     UT IS UNIVERSAL TIME IN SECONDS.
C     ****     IDAY IS DAY OF YEAR
C     ****
!
! Local:
      integer :: j,isun,jsun
      real :: W(4),glats,glons,pisun,pjsun,sndlons,csdlons
C     ****
C     ****     CALCULATE SUN'S GEOGRAPHIC COORDINATES.
C     ****
      GLATS=ASIN(.398749*SIN(2.*pi*(IDAY-80)/365.))
      GLONS=pi*(1.-2.*SECS/86400.)
C     ****
C     ****     CALCULATE SUN'S APEX LONGITUDE USING BILINEAR
C     ****       INTERPOLATION
C     ****
C     ****     FILL ARRAY DLONS(ZJMX) WITH THIS VALUE
C     ****
      PISUN = (GLONS+pi)/dlamda+1.
      PJSUN = (GLATS+.5*(pi-dphi))/dphi+1.
      ISUN = PISUN
      JSUN = PJSUN
      PISUN = PISUN-FLOAT(ISUN)
      PJSUN = PJSUN-FLOAT(JSUN)
      SNDLONS = (1.-PISUN)*(1.-PJSUN)*SIN(RLONM(ISUN+2,JSUN))+
     1           PISUN*(1.-PJSUN)*SIN(RLONM(ISUN+3,JSUN))+
     2           PISUN*PJSUN*SIN(RLONM(ISUN+3,JSUN+1))+
     3           (1.-PISUN)*PJSUN*SIN(RLONM(ISUN+2,JSUN+1))
      CSDLONS = (1.-PISUN)*(1.-PJSUN)*COS(RLONM(ISUN+2,JSUN))+
     1           PISUN*(1.-PJSUN)*COS(RLONM(ISUN+3,JSUN))+
     2           PISUN*PJSUN*COS(RLONM(ISUN+3,JSUN+1))+
     3           (1.-PISUN)*PJSUN*COS(RLONM(ISUN+2,JSUN+1))
      DLONS(1) = ATAN2(SNDLONS,CSDLONS)
      DO 1 J = 2,JMAX
        DLONS(J) = DLONS(1)
    1 CONTINUE
      RETURN
      END
C