C =============================================================================
C
      SUBROUTINE ADJUST_LL(IDAY,IYEAR,UT,GLAT,GLONG,SZA,STL)
C
C ROUTINE TO:
C - ESTIMATE THE SOLAR ZENITH ANGLE (SZA)
C - IF SZA IS LARGER THAN SPECIFIED TWILIGHT VALUE (SZAL), THEN:
C     COMPUTE THE LATITUDE & LONGITUDE FOR THE POINT WHERE SZA=SZAL
C     ALONG THE GREAT CIRCLE CONTAINING THE SUBSOLAR POINT
C     (SHADOW HEIGHT DEFINED RELATIVE TO 90 km LEVEL)
C - COMPUTE THE APPARENT SOLAR LOCAL TIME (STL) FOR MSIS CALLS
C updated: 25 November 97 (jbishop)
C *** not exact, but serviceable ***
C
      PI    = 2.0 * asin(1.0)
      R2D   = 180.0 / PI                    ! radian-to-degree conversion factor
      SZAL  = 100.50                        ! SZA where shadow hgt = 200km 
      cSZAL = cos(SZAL / R2D)
      sSZAL = sin(SZAL / R2D)
C
C place GLAT & GLONG in proper ranges: (90,-90) & (0,360)
      CALL RANGE(GLAT,-90.0,90.0,IU_LAT,ID_LAT,RANGE_LAT)
      CALL RANGE(GLONG,0.0,360.0,IU_LON,ID_LON,RANGE_LON)
      IF (GLAT.LE.-89.90)  GLAT  = -89.90
      IF (GLAT.GE. 89.90)  GLAT  =  89.90
      IF (GLONG.LE.  0.10) GLONG =   0.10
      IF (GLONG.GE.359.90) GLONG = 359.90
C
C crude conversion from Julian calendar to Pagan calender
      YEAR = 365.0
      IF (MOD(IYEAR,4).EQ.0) YEAR = 366.0
      DAY_PC = FLOAT(IDAY + 11) - 0.50
      CALL RANGE(DAY_PC,0.0,YEAR,IU_DAY,ID_DAY,RANGE_DAY)
C
C ecliptic plane coordinates
      EthetaS = 90.0 / R2D                  ! solar colatitude
      EthetaG = 23.441 / R2D                ! colatitude of rotation pole
C . . . "longitudes" in ecliptic plane relative to solar direction
      EphiS   = 0.0                         ! reference "longitude" (solar)
C . . . daily "longitude" of rotation pole (treated as an interior angle)
      EphiG   = abs(180.0 - 360.0*(DAY_PC/YEAR)) / R2D
C
C colatitude of the subsolar point in geocentric coordinate system
      cGthetaS = sin(EthetaG) * cos(EphiG)
      GthetaS  = acos(cGthetaS)
      sGthetaS = sin(GthetaS)
C
C observer location in geocentric coordinate system
      GthetaO = (90.0 - GLAT) / R2D         ! colatitude of observer
      GphiO   = GLONG / R2D                 ! east longitude of observer
      alpha   = (GLONG + UT*15.0 - 180.0) / R2D
C . . . alpha is the observer "longitude" wrt the subsolar meridian
C . . . (treated as an interior angle)
C alpha phase information
      CALL RANGE(alpha,-pi,pi,IU_a,ID_a,RANGE_a)
      IM_a = 1
      if (alpha.lt.0.0) then
        IM_a  = -1
        alpha = abs(alpha)
      end if
C
C solar zenith angle
      cSZA = cos(GthetaO)*cGthetaS + sin(GthetaO)*sGthetaS*cos(alpha)
      SZA  = acos(cSZA)
      sSZA = sin(SZA)
      SZA  = SZA * R2D
C
C compare with specified twilight value 
C find effective illumination lat & long if needed
      IF (SZA .GT. SZAL) THEN
c       write(06,01)SZA,SZAL
c01     format(1x,'ADJUSTING SZA:  ',f8.3,' --> ',f8.3)
C beta: interior angle at subsolar point between great circles containing
C       rotation pole and observer location
        cbeta = (cos(GthetaO) - cGthetaS*cSZA)/(sGthetaS*sSZA)
C colatitude of shadow illumination point
        cGthetaN = cSZAL*cGthetaS + sSZAL*sGthetaS*cbeta
        GthetaN  = acos(cGthetaN)
        sGthetaN = sin(GthetaN)
        GLATN    = 90.0 - GthetaN*R2D
c       write(06,02)GLAT,GLATN
c02     format(1x,'ADJUSTING GLAT: ',f8.3,' --> ',f8.3)
        GLAT = GLATN
C interior ("longitude") angle
        alN = acos((cSZAL - cGthetaN*cGthetaS)/(sGthetaN*sGthetaS))
        IF (alN.GE.alpha) THEN
          write(06,99)alpha,alN
 99       format(1x,'APPARENT ERROR IN LONGITUDE ADJUSTMENT:',2f8.3)
        END IF
C longitude of shadow illumination point
        alN    = IM_a*alN + (ID_a - IU_a)*RANGE_a
        GLONGN = alN*R2D + 180.0 - UT*15.0
        CALL RANGE(GLONGN,0.0,360.0,IU_LON,ID_LON,RANGE_LON)
c       write(06,03)GLONG,GLONGN
c03     format(1x,'ADJUSTING GLONG:',f8.3,' --> ',f8.3)
        GLONG  = GLONGN
      ELSE
c       write(06,04)SZA
c04     format(1x,'NOT ADJUSTING SZA:',f8.3)
      ENDIF
C
      STL = UT + GLONG/15.0
      CALL RANGE(STL,0.0,24.0,IU_T,ID_T,RANGE_T)
C
      RETURN
      END
C
C -----------------------------------------------------------------------------
C
      SUBROUTINE RANGE(T,T0,T1,IU,ID,DT)  ! put periodic 't' into range t0...t1
C 
      DT = T1 - T0
      IU = 0
      ID = 0
      DO WHILE (T.LT.T0 .OR. T.GT.T1)
        IF (T .LT. T0) THEN
          T  = T + DT
          IU = IU + 1
        ELSEIF (T .GT. T1) THEN
          T  = T - DT
          ID = ID + 1
        ENDIF
      ENDDO
      RETURN
      END
C
C =============================================================================