c  gfortran dumpMapLRO_RG2.f -O2
c
c  Version 1.0 - 17 February 2012 - Eric Palmer (from Gaskell)
c		This outputs a bigmap as a simple csv textfile, grid format
c	Version 1.1 - 31 May 2012
c		Outputs the flux on the surface based on solar lat/lon
c  Version 1.2 - 11 July 2013
c     Outputs flux, equlibTemp, i, e, a
c  Version 1.3 -  4 Sep 2013
c  Version 1.4 - 9 Jan 2014 - just added cubes and i,e into degrees
c  Version 1.5 - 14 Sept 2018 - writes out the flip values as a grid

      IMPLICIT NONE

      INTEGER               NTMP
      PARAMETER            (NTMP=2000)

      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      UX(3)
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)

      REAL*4                HT0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                AL0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                slope(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                inc(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                emm(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                baseEmm(-NTMP:NTMP,-NTMP:NTMP)
      REAL                  tempA (-NTMP:NTMP,-NTMP:NTMP)
      REAL                  fluxA (-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                thermSpec(30000)
      REAL*4                dHx, dHy
      REAL*4                mVal, albScale
      REAL*4                minV, delta, holdMin
      REAL*4                a, b, c, cosA, cosC
      REAL*4                bigC, bigB, doubleB, bigD
      REAL*4                aPrime, bPrime, cPrime, dPrime
      REAL*8                SCx, SCz, Nx, Ny, Nz
      REAL*8                temp, SCm, Nm
      REAL*4                angle, avgTemp, avgAlb, avgInc
      REAL*4                azmuth, elevation
      REAL*4                sunLat, sunLon, objLat, objLon
      REAL*4                spaceLat, spaceLon
      REAL*4                Pi
      REAL*4                LonInput
      REAL*4                solarDist
      REAL*4                g_version
      REAL*4                targetLat, targetLon
      REAL*4                subsolarLat, subsolarLon
      REAL*4                  phi1, phi2, tanAlpha
      REAL*4               MevflipX(-NTMP:NTMP,-NTMP:NTMP)

      INTEGER               QSZ
      INTEGER               QSZLESS
      INTEGER               I
      INTEGER               J
      INTEGER               K
      INTEGER               Istart, Iend, Jstart, Jend
      INTEGER               dX, dY, d
      INTEGER               cnt
      integer               fluxF, tempF, iF, eF, aF, slopeF
      integer               whichImage
      integer               argc
      integer               evinitX, evflipX, imgBreak, flip, nslope

      CHARACTER*6           MAP0
      CHARACTER*72          LMRKFILE
      CHARACTER*72          outName
      character*6          whichImgStr
      character*80          argv
      CHARACTER*72          EmmName
      g_version = 1.4

      argc = IArgC()
      if (argc .NE. 3) then
         WRITE(6,*) '# Input input map name'
         READ(5,FMT='(A6)') MAP0
         WRITE(6,*) '# Which image (3611 324 396 484 047): '
         READ(5,*) whichImgStr
c         WRITE(6,*) 'whichImgStr', whichImgStr
c        READ(5,*) whichImage
         WRITE(6,*) '# Longitude (1000 is default): '
         READ(5,*) LonInput
       else
         call getArg (1, MAP0)
         call getArg (2, whichImgStr)
c         READ(argv, *) whichImage
         call getArg (3, argv)
         READ(argv, *) LonInput
      endif


      LMRKFILE='MAPFILES/'//MAP0//'.MAP'
      CALL READ_MAP(LMRKFILE,NTMP,QSZ,SCALE,V,UX,UY,UZ,HT0,AL0)
      Pi = 3.1415926

      WRITE(6,*) '# QSize ', QSZ
      if (qsz .gt. NTMP) then
        write (*,*) "Qsize is too big for this program: ", qsz, NTMP
        stop
      endif


c ########################################
c Populate baseEmm with fake values
c ########################################
c     Creat array for testing purposes. Replace with actual baseEmission.txt
c     before generating final values.

      Istart = -qsz
      Iend = qsz
      Jstart = -qsz
      Jend = qsz

      DO J=Jstart,Jend
      DO I=Istart,Istart/3

       baseEmm (I, J) = 0.0 * Pi / 180;

      ENDDO
      ENDDO

      DO J=Jstart,Jend
      DO I=Istart/3,Iend/3

       baseEmm (I, J) = 0.0 * Pi / 180;

      ENDDO
      ENDDO
      DO J=Jstart,Jend
      DO I=Iend/3,Iend

       baseEmm (I, J) = 0.0 * Pi / 180;

      ENDDO
      ENDDO


c      WRITE(6,*) 'Test', baseEmm (2, 1)

c ########################################
c Reading raw Emission file
c ########################################


      EmmName='../img-' // whichImgStr // '/rawEmm.txt'
      WRITE (6,*) whichImgStr 
      OPEN (UNIT=10, FILE=EmmName, STATUS='OLD')

      Istart = -qsz
      Iend = qsz
      Jstart = -qsz
      Jend = qsz

      DO J=Jstart,Jend
       DO I=Istart,Iend-1
        READ (10,"(f12.5)",ADVANCE='NO'), baseEmm (I,J)
        if (J .eq. 43) then
         WRITE(6,*) 'I,J,baseEmm', I, J, baseEmm (I,J)
        ENDIF
        baseEmm (I,J) = baseEmm (I,J) * Pi / 180 
c        baseEmm (I,J) = 0

       ENDDO
      READ (10,"(f12.5)"), baseEmm (Iend+1,J)
      baseEmm (Iend+1,J) = baseEmm (Iend+1,J) *Pi/180 
c      baseEmm (Iend+1,J) = 0
      ENDDO



      WRITE(6,*) 'From file', baseEmm (0, 2)




c ########################################
c Set location --- this is all Tsiolkovsky
c ########################################

      targetLat = 90 - (-19.5753)                                        3611LE
      targetLon = 128.441


c ####################
c ######3 Update location -- Karpinskiy 
c ####################
      evflipX = 1                                                       Always start at 1
      flip = 1                                                          flip=1: check for inflection; flip=0: don't check for inflection
      nslope = -1                                                       nslope; determine if slope is negative or just 0
      imgBreak = -2000                                                  Line number where image is no longer at edge of bigmap -2,000 means not used

      targetLat = 90 - (7.013)                                          3611LE
      targetLon = 300.802                                               Converted from W long to E long 

      spaceLon = 0
      spaceLat = 0

      solarDist = 0

c      if (whichImgStr .EQ. "") then                                     RG2 
c         write (*,*) "# Which is "
c         subsolarLat = 90 - ()
c         subsolarLon = 
c         solarDist = 1
c         evinitX = 
c        S/C to at start (), inflection
c         flip=
c      endif

      if (whichImgStr .EQ. "00406L") then                                     RG2 
         write (*,*) "# Which is 00406L"
         subsolarLat = 90 - (0.51977745054861)
         subsolarLon = 9.8264237635933 
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "01166L") then                                     RG2 
         write (*,*) "# Which is 01166L"
         subsolarLat = 90 - (0.81026907908644)
         subsolarLon = 317.41949921189
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "01166R") then                                     RG2 
         write (*,*) "# Which is 01166R"
         subsolarLat = 90 - (0.81026907908644)
         subsolarLon = 317.41949921189
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "01476L") then                                     RG2 
         write (*,*) "# Which is 01476L"
         subsolarLat = 90 - (-0.74605290624147)
         subsolarLon = 315.05684062282
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), inflection
         flip=1
      endif

      if (whichImgStr .EQ. "02532R") then                                     RG2 
         write (*,*) "# Which is 02532R"
         subsolarLat = 90 - (-1.5069388875502)
         subsolarLon = 288.76116897645
         solarDist = 1
         evinitX = -1
c        S/C to at west start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "03387L") then                                     RG2 
         write (*,*) "# Which is 03387L"
         subsolarLat = 90 - (1.2240730010547)
         subsolarLon = 234.78419506045
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "03387R") then                                     RG2 
         write (*,*) "# Which is 03387R"
         subsolarLat = 90 - (1.2240730010547)
         subsolarLon = 234.78419506045
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "05753L") then                                     RG2 
         write (*,*) "# Which is 05753L"
         subsolarLat = 90 - (1.5244751624177)
         subsolarLon = 258.82905828898
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "23648L") then                                     RG2 
         write (*,*) "# Which is 23648L"
         subsolarLat = 90 - (-1.5113063243896)
         subsolarLon = 306.4086552046
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "23648R") then                                     RG2 
         write (*,*) "# Which is 23648R"
         subsolarLat = 90 - (-1.5113063243896)
         subsolarLon = 306.4086552046
         solarDist = 1
         evinitX = -1
c        S/C to east at start (+), flip at start (-), inflection
         flip=1
      endif

      if (whichImgStr .EQ. "27006L") then                                     RG2 
         write (*,*) "# Which is 27006L"
         subsolarLat = 90 - (-1.4104384869374)
         subsolarLon = 287.94646092923
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "34115R") then                                     RG2 
         write (*,*) "# Which is 34115R"
         subsolarLat = 90 - (-1.4098346098625)
         subsolarLon = 286.94239639932
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "36848L") then                                     RG2 
         write (*,*) "# Which is 36848L"
         subsolarLat = 90 - (0.12983134700247)
         subsolarLon = 227.80073723438
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif



      if (whichImgStr .EQ. "39671R") then                                     RG2 
         write (*,*) "# Which is 39671R"
         subsolarLat = 90 - (-0.24792159944025)
         subsolarLon = 342.57024302524
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "44013L") then                                     RG2 
         write (*,*) "# Which is 44013L"
         subsolarLat = 90 - (0.12752807864847)
         subsolarLon = 226.78554877913
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "46783L") then                                     RG2 
         write (*,*) "# Which is 46783L"
         subsolarLat = 90 - (-0.24556512225161)
         subsolarLon = 341.5697079068
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "50376L") then                                     RG2 
         write (*,*) "# Which is 50376L"
         subsolarLat = 90 - (0.10138843204491)
         subsolarLon = 345.90118459577
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "51166L") then                                     RG2 
         write (*,*) "# Which is 51166L"
         subsolarLat = 90 - (0.12522606941235)
         subsolarLon = 225.77199079771
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "52583R") then                                     RG2 
         write (*,*) "# Which is 52583R"
         subsolarLat = 90 - (-1.1396139218671)
         subsolarLon = 316.92778937333
         solarDist = 1
         evinitX = -1
c        S/C to east at start (+), flip at start (-) inflection
         flip=1
      endif

      if (whichImgStr .EQ. "53942R") then                                     RG2 
         write (*,*) "# Which is 53942R"
         subsolarLat = 90 - (1.22708138979)
         subsolarLon = 290.41023591436
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), inflection
c       Only a sliver in the upper right is east, rest of image is west
         imgBreak = -431
         flip=1
      endif

      if (whichImgStr .EQ. "64615L") then                                     RG2 
         write (*,*) "# Which is 64615L"
         subsolarLat = 90 - (0.095700623472075)
         subsolarLon = 343.89581210945
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "68263L") then                                     RG2 
         write (*,*) "# Which is 68263L"
         subsolarLat = 90 - (0.2865466791321)
         subsolarLon = 12.166609304412
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "70190L") then                                     RG2 
         write (*,*) "# Which is 70190L"
         subsolarLat = 90 - (1.5287095592712)
         subsolarLon = 263.83430013832
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "73381R") then                                     RG2 
         write (*,*) "# Which is 73381R"
         subsolarLat = 90 - (-1.5484455216581)
         subsolarLon = 260.70223131944
         solarDist = 1
         evinitX = -1
c        S/C to east at start (+), flip at start (-), inflection
         flip=1
      endif

      if (whichImgStr .EQ. "86730R") then                                     RG2 
         write (*,*) "# Which is 86730R"
         subsolarLat = 90 - (-0.94338895933418)
         subsolarLon = 314.60070253532
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "87807R") then                                     RG2 
         write (*,*) "# Which is 87807R"
         subsolarLat = 90 - (0.49732805487309)
         subsolarLon = 349.5995322891
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "91527L") then                                     RG2 
         write (*,*) "# Which is 91527L"
         subsolarLat = 90 - (1.5261585445841)
         subsolarLon = 260.83128568039
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "95388L") then                                     RG2 
         write (*,*) "# Which is 95388L"
         subsolarLat = 90 - (-1.5076894988216)
         subsolarLon = 289.76967067447
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif

      if (whichImgStr .EQ. "98641L") then                                     RG2 
         write (*,*) "# Which is 98641L"
         subsolarLat = 90 - (1.5253150196408)
         subsolarLon = 259.83026896324
         solarDist = 1
         evinitX = -1
c        S/C to west at start (), no inflection
         flip=0
      endif




      if (whichImgStr .EQ. "5634L") then                                Kar 
         write (*,*) "# Which is 5634L"
         subsolarLat = 90 - (-0.3278127349197)
         subsolarLon = 201.16252431651
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif



      if (whichImgStr .EQ. "7899R") then                                Kar 
         write (*,*) "# Which is 7899R"
         subsolarLat = 90 - (-0.77673380353124)
         subsolarLon = 88.649639918528
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif



      if (whichImgStr .EQ. "3150R") then                                Kar 
         write (*,*) "# Which is 3150R"
         subsolarLat = 90 - (1.5538621855453)
         subsolarLon = 140.60305506547
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "6002L") then                                Kar 
         write (*,*) "# Which is 6002L"
         subsolarLat = 90 - (1.5536242633519)
         subsolarLon = 141.61267107067
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "5426R") then                                Kar 
         write (*,*) "# Which is 5426R"
         subsolarLat = 90 - (0.71898867807073)
         subsolarLon = 84.619710079862
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "8277L") then                                Kar 
         write (*,*) "# Which is 8277L"
         subsolarLat = 90 - (0.71724412060989)
         subsolarLon = 85.62769603618
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "9193L") then                                Kar 
         write (*,*) "# Which is 9193L"
         subsolarLat = 90 - (0.97665202421357)
         subsolarLon = 159.93557951658
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "1974R") then                                Kar 
         write (*,*) "# Which is 1974R"
         subsolarLat = 90 - (0.97844172889283)
         subsolarLon = 160.95533463746
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "2701L") then                                Kar 
         write (*,*) "# Which is 2701L"
         subsolarLat = 90 - (1.4159704015655)
         subsolarLon = 133.04009362394
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "3453L") then                                Kar 
         write (*,*) "# Which is 3453L"
         subsolarLat = 90 - (-1.4378707588627)
         subsolarLon = 171.53195424355
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), inflection
         flip=0
      endif


      if (whichImgStr .EQ. "6668R") then                                Kar 
         write (*,*) "# Which is 6668R"
         subsolarLat = 90 - (-1.4385758664828)
         subsolarLon = 172.49201162202
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "4495R") then                                Kar 
         write (*,*) "# Which is 4495R"
         subsolarLat = 90 - (-1.5179158057432)
         subsolarLon = 144.27548457952
         solarDist = 1
         if (MAP0 .EQ. "KAR512") then
         evinitX = 1
c        S/C to east at start (+), flip at image start
         evinitX = -1
c        inflection point occurs in lower left, but not top left
         imgBreak = -80
         endif
         if (MAP0 .EQ. "KAR522") then
          evinitX = 1
c         S/C to east at start (+), no flip at start
c        image no longer at edge of bigmap
          imgBreak = 7
         else
          evinitX = 1
c         S/C to east at start (+), flip at start
          evinitX = -1
         endif

      endif


      if (whichImgStr .EQ. "4495L") then                                Kar 
         write (*,*) "# Which is 4495L"
         subsolarLat = 90 - (-1.5179158057432)
         subsolarLon = 144.27548457952 
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "2919R") then                                Kar 
         write (*,*) "# Which is 2919R"
         subsolarLat = 90 - (-0.2495808884134)
         subsolarLon = 217.09964982763
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "5991R") then                                Kar 
         write (*,*) "# Which is 5991R"
         subsolarLat = 90 - (0.70502327783508)
         subsolarLon = 237.11276963531
         solarDist = 1
         if (MAP0 .EQ. "KAR512") then
         evinitX = 1
c        S/C to east at start (+), no flip at start 
c        image no longer at edge of bigmap
         imgBreak=410
         endif
         if (MAP0 .EQ. "KAR522") then
          evinitX = 1
c         S/C to east at start (+), no flip at start
c        image no longer at edge of bigmap
          imgBreak = -476
         else
          evinitX = 1
c         S/C to east at start (+), flip at start
          evinitX = -1
         endif
      endif



      if (whichImgStr .EQ. "9559L") then                                Kar 
         write (*,*) "# Which is 9559L"
         subsolarLat = 90 - (-0.030238435728367)
         subsolarLon = 217.80651629117
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "4525L") then                                Kar 
         write (*,*) "# Which is 4525L"
         subsolarLat = 90 - (0.26330895354472)
         subsolarLon = 101.70917444942
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif



      if (whichImgStr .EQ. "9879L") then                                Kar 
         write (*,*) "# Which is 9879L"
         subsolarLat = 90 - (-0.79504629017298)
         subsolarLon = 236.57773368435
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif


      if (whichImgStr .EQ. "9667L") then                                Kar 
         write (*,*) "# Which is 9667L"
         subsolarLat = 90 - (1.4029867070325)
         subsolarLon = 108.4010127248
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), no inflection
         flip=0
      endif




      if (whichImgStr .EQ. "7794L") then                                Kar 
         write (*,*) "# Which is 7794L"
         subsolarLat = 90 - (1.53)
         subsolarLon = 136.51
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+)
      endif

      if (whichImgStr .EQ. "4024R") then                                Kar 
         write (*,*) "# Which is 4024R"
         subsolarLat = 90 - (-1.3857058015866)
         subsolarLon = 169.74778851451
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), no inflection
         flip=0
      endif



      if (whichImgStr .EQ. "3608L") then                                Kar
         write (*,*) "# Which is 3608L"
         subsolarLat = 90 - (-0.77263867127241)
         subsolarLon = 90.671414899152
         solarDist = 1
         if (MAP0 .EQ. "KAR512") then
          evinitX = 1
c         S/C to east at start (+), no inflection
          evinitX = 1
         endif
         if (MAP0 .EQ. "KAR522") then
          evinitX = 1
c         S/C to east at start (+), no inflection
c         disable search for inflection
          flip=0
          evinitX = 1
         endif
         if (MAP0 .EQ. "KAR532") then
          evinitX = 1
c         S/C to east at start (+), flip when image starts (-) 
          evinitX = -1
         endif
         if (MAP0 .EQ. "KAR542") then
          evinitX = 1
c         S/C to east at start (+), flip when images starts (-)
          evinitX = -1
         else
          write (*,*) "WRONG MAP NAME!"
         endif
      endif

      if (whichImgStr .EQ. "1721L") then                                Kar
         write (*,*) "# Which is 1721L"
         subsolarLat = 90 - (-1.5285378266391)
         subsolarLon = 181.04286212524
         solarDist = 1
         if (MAP0 .EQ. "KAR532") then
          evinitX = -1
c         S/C to west at start (-), flip when image starts (-) 
          evinitX = 1
         endif
         if (MAP0 .EQ. "KAR542") then
          evinitX = -1
c         S/C to west at start (-), flip when images starts (-)
          evinitX = 1
         else
          write (*,*) "WRONG MAP NAME!"
         endif
      endif





      if (whichImgStr .EQ. "0544L") then                                Ing
         write (*,*) "# Which is 0544L"
         subsolarLat = 90 - (-1.3740831506497)
         subsolarLon = 160.33160152693
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), flip when image starts (-)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "0544R") then                                Ing
         write (*,*) "# Which is 0544R"
         subsolarLat = 90 - (-1.3740831506497)
         subsolarLon = 160.33160152693
         solarDist = 1
         evinitX = 1
c        S/C to east at start (+), flip when image starts (-)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "1128L") then                                Ing
         write (*,*) "# Which is 1128L"
         subsolarLat = 90 - (-0.74336506853272)
         subsolarLon = 245.41937722448
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), flip when images starts (-)
         evinitX = 1
      endif

      if (whichImgStr .EQ. "1128R") then                                Ing
         write (*,*) "# Which is 1128R"
         subsolarLat = 90 - (-0.74336506853272)
         subsolarLon = 245.41937722448
         solarDist = 1
         evinitX = -1
c        S/C to west at start (-), flip when image starts (-)
         evinitX = 1
      endif

      if (whichImgStr .EQ. "1607R") then                                Ing
         write (*,*) "# Which is 1607R"
         subsolarLat = 90 - (-0.60294927647943)
         subsolarLon = 83.931636669465
         solarDist = 1
         evinitX = 1
c        Inflection point is at ~330 pixels above middle, isn't one above that 
c        Need to disable the search for the inflection point for this image 
         flip=0
c        S/C east/west switches after first ~5 pixels, 99% of image is to West (-)
         evinitX = -1
      endif


      if (whichImgStr .EQ. "2292L") then                                Ing
         write (*,*) "# Which is 2292L"
         subsolarLat = 90 - (-0.7893601030968)
         subsolarLon = 233.41712965056
         solarDist = 1
         evinitX = 1
c        S/C to the west (-), flip when image starts (-)
         evinitX = 1
      endif

      if (whichImgStr .EQ. "2292R") then                                Ing
         write (*,*) "# Which is 2292R"
         subsolarLat = 90 - (-0.7893601030968)
         subsolarLon = 233.41712965056
         solarDist = 1
         evinitX = 1
c        S/C to the east (+), flip when image starts (-)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "3890L") then                                Ing
         write (*,*) "# Which is 3890L"
         subsolarLat = 90 - (0.55419872375988)
         subsolarLon = 243.4260606642
         solarDist = 1
         evinitX = -1
c       S/C to the west (-), flip when image starts (-) 
         evinitX = 1
      endif

      if (whichImgStr .EQ. "3890R") then                                Ing
         write (*,*) "# Which is 3890R"
         subsolarLat = 90 - (0.55419872375988)
         subsolarLon = 243.4260606642
         solarDist = 1
         evinitX = 1
c        Inflection about 200 pixles from bottom, only a few pixels
         flip=0
c        S/C in west for 99% of images (-)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "4018R") then                                Ing
         write (*,*) "# Which is 4018R"
         subsolarLat = 90 - (-0.74137451075919)
         subsolarLon = 246.42604538998
         solarDist = 1
         evinitX = 1
c        S/C in east (+), flip when image starts (-) 
         evinitX = -1
      endif

      if (whichImgStr .EQ. "4503L") then                                Ing
         write (*,*) "# Which is 4503L"
         subsolarLat = 90 - (-0.74242832660511)
         subsolarLon = 184.97054700561
         solarDist = 1
         evinitX = -1
c        S/C is in the west (-), no inflection (+)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "5227L") then                                Ing
         write (*,*) "# Which is 5227L"
         subsolarLat = 90 - (-1.3003342783437)
         subsolarLon = 106.70855819665
         solarDist = 1
         evinitX = 1
c        S/C is in east (+), flip when image starts (-)
         evinitX = -1
      endif

      if (whichImgStr .EQ. "5229R") then                                Ing
         write (*,*) "# Which is 5229R"
         subsolarLat = 90 - (1.5542431343956)
         subsolarLon = 138.90056433247
         solarDist = 1
         evinitX = -1
c        S/C is in west (-), no inflection point
         evinitX = -1
      endif

      if (whichImgStr .EQ. "7476L") then                                Ing
         write (*,*) "# Which is 7476L"
         subsolarLat = 90 - (0.72192225278666)
         subsolarLon = 82.925974813148
         solarDist = 1
         evinitX = 1
c        S/C is in the east (+), no inflection point        
         evinitX = 1
      endif

      if (whichImgStr .EQ. "7476R") then                                Ing
         write (*,*) "# Which is 7476R"
         subsolarLat = 90 - (0.72192225278666)
         subsolarLon = 82.925974813148
         solarDist = 1
         evinitX = 1
c        S/C is in the east (+), flip when image starts (-)
c        Flips in many places, need to investigate more
         evinitX = -1
      endif

      if (whichImgStr .EQ. "7807R") then                                Ing
         write (*,*) "# Which is 7807R"
         subsolarLat = 90 - (-0.0050630550821127)
         subsolarLon = 213.17267063247
         solarDist = 1
         evinitX = 1
c        S/C in east (+), flip when image starts (-)
         evinitX = -1
      endif



















      if (whichImage .EQ. 3611) then
         write (*,*) "# Which is 3611"
         subsolarLat = 90 - (-1.452126)
         subsolarLon = 120.436497
         solarDist = 1
      endif

      if (whichImage .EQ. 324) then
         write (*,*) "# Which is 324"
         subsolarLat = 90 - (-.2249707)
         subsolarLon = 67.522966
         solarDist = 1
         spaceLat = 0
         spaceLon = -8.28346
      endif

      if (whichImage .EQ. 396) then
         write (*,*) "# Which is 396"
         subsolarLat = 90 - (1.4524)
         subsolarLon = 119.4785
         solarDist = 1
         spaceLat = 0
         spaceLon = 0
      endif

      if (whichImage .EQ. 484) then
         write (*,*) "# Which is 484"
         subsolarLat = 90 - (-0.227120507)
         subsolarLon = 66.511393
         solarDist = 1
         spaceLon = 11.18313
      endif

      if (whichImage .EQ. 047) then
         write (*,*) "# Which is 047"
         subsolarLat = 90 - (-1.538407)
         subsolarLon = 147.21658
         spaceLon = 1.1106531
      endif
      if (whichImage .EQ. 048) then
         write (*,*) "# Which is 048"
         subsolarLat = 90 - (-1.2326614326948)
         subsolarLon = 79.22908512374
         spaceLon = targetLon - 127.7208
      endif
      if (whichImage .EQ. 261) then
         write (*,*) "# Which is 261"
         subsolarLat = 90 - (-1.2306523046095)
         subsolarLon = 80.188202798133
         spaceLon = targetLon - 128.75202930418
      endif
      if (whichImage .EQ. 506) then
         write (*,*) "# Which is 506"
         subsolarLat = 90 - (-1.147997924108)
         subsolarLon = 166.44806133895
         spaceLon = targetLon - 128.08030736391
      endif
      if (whichImage .EQ. 567) then
         write (*,*) "# Which is 567"
         subsolarLat = 90 - (-0.91302955683526)
         subsolarLon = 93.377705578176
         spaceLon = targetLon - 128.54458371738
      endif
      if (whichImage .EQ. 691) then
         write (*,*) "# Which is 691"
         subsolarLat = 90 - (0.16956523258729)
         subsolarLon = 190.04470902811
         spaceLon = targetLon - 128.21639323159
      endif
      if (whichImage .EQ. 723) then
         write (*,*) "# Which is 723"
         subsolarLat = 90 - (-1.1490872543041)
         subsolarLon = 167.40619347527
         spaceLon = targetLon - 128.75202930418
      endif
      if (whichImage .EQ. 5062) then
         write (*,*) "# Which is 5062"
         subsolarLat = 90 - (-1.147997924108)
         subsolarLon = 166.44806133895
         spaceLon = targetLon - 128.08030736391
      endif

      sunLat = subsolarLat - targetLat
c     1000 means use exact observation conditions
      if (LonInput .EQ. 1000) then
         sunLon = subsolarLon - targetLon
      else
         sunLon = LonInput
      endif

      write (*,*) "################## Parameters ###############"
      write (*,*) "# Target: ", targetLat, targetLon
      write (*,*) "# SubSolar: ", subsolarLat, subsolarLon
      write (*,*) "# Delta Solar: ", sunLat, sunLon
      write (*,*) "# Spacecraft: ", spaceLat, spaceLon
      write (*,*) "# solarDist: ", solarDist

      targetLat = targetLat / 180.0 * Pi
      targetLon = targetLon / 180.0 * Pi
      subsolarLat = subsolarLat / 180.0 * Pi
      subsolarLon = subsolarLon / 180.0 * Pi
      sunLat = sunLat / 180.0 * Pi
      sunLon = sunLon / 180.0 * Pi
      spaceLon = spaceLon / 180.0 * Pi
      spaceLat = spaceLat / 180.0 * Pi


c ########################################
c Angles
c ########################################

c     How low in the horizon is the Sun
c		Calculated by taking the atan of sunLat and sunLon

c     Elevation: how high the Sun is above the horzion
c     Get C (the angular distance of how far "down" the sun is

c           Great circle formula - cosine rule for angles
      cosA = sin(targetLat) * sin(subsolarLat)
     +            + cos (targetLat) * cos (subsolarLat) * cos (sunLon)

      A = acos (cosA)
      write (*,*) "# Zenith Angle: ", A * 180/Pi
      elevation = -1 * (Pi/2.0 -  A )
      write (*,*) "# look angle: ", elevation * 180 / Pi

c     Three step caluclation of zenith angle - allows for non-equator point
c        lon is all based on the target being zero
c        latitude must be fully handled
      a = sunLon
      b = subsolarLat
      bigC = Pi/2

      c = acos ( cos(a) * cos (b) + sin(a)*sin(b) * cos (bigC))
      write (*,*) "##########################"
      write (*,*) "# step1, c: ", c*180/Pi

      bigB = asin ( sin (bigC) * sin (b) / sin (c))
      if (targetLat < 0) then
         doubleB = Pi/2 + bigB
      else
         doubleB = Pi/2 - bigB
      endif

      write (*,*) "##########################"
      write (*,*) "# step2, bigB: ", bigB*180/Pi
      write (*,*) "# step2, doubleB: ", doubleB*180/Pi

c     We assign c to be "d" - they are the same
      dPrime = c
      aPrime = targetLat
      bPrime = acos (cos(dPrime) * cos(aPrime)
     +                + sin(dPrime) * sin(aPrime) * cos(doubleB))
      write (*,*) "##########################"
      write (*,*) "# step3, bPrime: ", bPrime*180/Pi
      write (*,*) "##########################"

c     This is the azmuth towards the Sun. The direction we look
c        for cliffs throwing shadows.  Angles is typical grid
c        This angle is in geographic coordinates (N is 0)
c        It gets the ground position of the subsolar location
c        azmuth is used in shadow
      B = sunLon
      C = sunLat
      bigC = (cos (C) - cos(A)*cos(B)) / sin(A) / sin(B)
      bigC = acos (bigC);
      azmuth = bigC
      if (subsolarLat .lt. targetLat) then
         azmuth = -1 * azmuth
      endif
      if (azmuth .LT. 0) azmuth = azmuth + 2 * Pi
      write (*,*) "# old azmuth: ", azmuth * 180.0/Pi


c     Use great-circle navigation and solve for sigma_12
c        It returns "bearing", which is 0 is north
c        To convert it to plotting azmuith = 90 - heading
      phi1 = targetLat
      phi2 = subsolarLat
      azmuth = atan2 ( sin(sunLon) * cos (phi2),
     +      cos (phi1)*sin(phi2) - sin(phi1) * cos(phi2) * cos (sunLon))
      azmuth = Pi/2 - azmuth
      if (azmuth .LT. 0) azmuth = azmuth + 2 * Pi
      write (*,*) "# plot azmiuth: ", azmuth*180/Pi

c ########################################
c Slope
c ########################################
c		 Loop to sample a large area around what I want to see

c     Wide field of view
      Istart = -qsz
      Iend = qsz
      Jstart = -qsz
      Jend = qsz

      WRITE(6,*) '# Version: ', g_version
      WRITE(6,*) '# Calculating slope'
      DO J=Jstart,Jend
c New line, so reset evflipX, and nslope
      evflipX = 1
      nslope =-1
c     If part of image is at edge of bigmap and part is not...
      if (imgBreak .gt. -2000) then
       if (J .ge. imgBreak) then
        evflipX = -evflipX
       end if
      end if

      DO I=Istart,Iend

c     This assumes 5 meter grid spacing and a skip over the center point
      dHx = (HT0 (I+1,J) - HT0 (I-1,J) ) /2
      dHy = (HT0 (I, J+1) - HT0 (I, J-1) ) /2

c      Slope in percent
c      write (*,*) "#dHx ", dHx;
c      write (*,*) "#dHy ", dHy;

      dHx = atan (dHx)
      dHy = atan (dHy)

c      Display the slope in degrees
c      write (*,*) "#dHx degrees: ", dHx * 180/Pi
c      write (*,*) "#dHy degrees: ", dHy * 180/Pi

c      Vector (does not need the end of loop LF)
c      write (*,*) I-Istart+1, J-Jstart+1, dHx, dHy

c     Incidence angle, cos(i)
      objLat = dHy;
      objLon = dHx;
      A = sunLon - dHx;                                   
      A = sunLon + dHx;                                                 vesta
      B = Pi/2.0 - dHy
      C = Pi/2.0 - sunLat
      cosA = cos(B) * cos(C) + sin (B) * sin(C) * cos (A)
      if (cosA .lt. 0) then
         cosA = 0
      endif


      slope(I, J) = acos (cosA);
      inc(I, J) = acos (cosA);

c     Emission angle
c old      A = spaceLon - objLon;
c old      B = Pi/2.0 - objLat
c old      C = Pi/2.0 - spaceLat
c old      cosA = cos(B) * cos(C) + sin (B) * sin(C) * cos (A)
c old      emm (I, J) = acos (cosA);

c     New Emission angle calculation (JRW)
c     Use base emission angle from ISIS, then use vectors to modify the base
c     angle. First change emission angle into spacecraft vector that has x
c     and z. (y is North/South and for LRO is always 0). Need to determine if
c     x is pos or neg since emission angle alone doesn't determine East or
c     West. Determine height difference in y and x, the cross product of y
c     and x is the normal at that pixel. Use dot-product of the normal and
c     spacecraft vector to get the angle betweeen them; this final value is
c     the modified Emission angle. 



c     Calculate spacecraft vector

c Somewhere need to do a check to see if left of images is pos or neg x (i.e. S/C west or east of topography
c Need to add variables for these two values so it will fit on a single line
      if (flip .eq. 1)  then
       if (J .eq. 43) WRITE(6,*) 'I', I, 'nslope', nslope
       if (baseEmm(I,J)-baseEmm(I-1,J) .le. 0) then
c        WRITE(6,*) baseEmm(I+1,J), baseEmm(I,J), 'diff', baseEmm(I+1,J) - baseEmm(I,J)
        if (baseEmm(I,J)-baseEmm(I-1,J) .lt. 0) then
         if (J .eq. 43) WRITE(6,*) 'I', I
         nslope = -1
        end if
        if (baseEmm(I+1,J)-baseEmm(I,J) .gt. 0) then
         if (nslope .eq. -1) then
          evflipX = -evflipX
          nslope = 1
          if (J .eq. 43) then
           WRITE(6,*) 'X of s/c vector has flipped. I,J=', I, J
           WRITE(6,*) 'I-1', baseEmm(I-1,J), 'I', baseEmm(I,J)
           WRITE(6,*) 'I+1', baseEmm(I+1,J)
          end if
         end if
        end if
       end if
      end if


c     Record evflipX for transparency
      MevflipX(I,J) = evflipX

      SCx = evinitX * evflipX * tan (baseEmm(I,J)) * cos (baseEmm(I,J))
      SCz = cos (baseEmm(I,J))
      SCm = SQRT(SCx*SCx + SCz*SCz)

c     Calculate normal at pixel
c     Bring back dHx and dHy

c     Change HT0 values so we can see the math easily
c
c      HT0(I+1,J)=2
c      HT0(I-1,J)=2
c      HT0(I,J+1)=2
c      HT0(I,J-1)=0
c
c     END Changing of HT0 values

      dHx = (HT0 (I+1,J) - HT0 (I-1,J) )
      dHy = (HT0 (I, J+1) - HT0 (I, J-1) )


c     Cross product results from y slope vector and x slope vector
      Nx = -dHx * 2.0
      Ny = -dHy * 2.0
      Nz = 4.0
      Nm = SQRT(Nx*Nx + Ny*Ny + Nz*Nz)

c     Testing only below take out before final analysis!
      if (I .eq. -96) then
       if (J .eq. 335) then
        emm(I,J) = ACOS((Nx*SCx + Nz*SCz)/(Nm*SCm))*180.0/Pi
        WRITE(6,*) 'Normal Vector', Nx, Ny, Nz, Nm
        WRITE(6,*) 'Normal Vector Sq', Nx*Nx, Ny*Ny, Nz*Nz
        WRITE(6,*) 'HTO J direction', HT0(I, J+1), HT0(I, J-1)
        WRITE(6,*) 'HTO I direction', HT0(I+1,J), HT0(I-1,J)

        WRITE(6,*) 'emm I direction', emm(I-1,J), emm(I,J), emm(I+1,J)
        WRITE(6,*) 'emm J direction', emm(I,J-1), emm(I,J), emm(I,J+1)
        WRITE(6,*) 'emm I direction', emm(I-2,J), emm(I,J), emm(I+2,J)
        WRITE(6,*) 'emm J direction', emm(I,J-2), emm(I,J), emm(I,J+2)
        WRITE(6,*) 'emm I direction', emm(I-3,J), emm(I,J), emm(I+3,J)
        WRITE(6,*) 'emm J direction', emm(I,J-3), emm(I,J), emm(I,J+1)


        WRITE(6,*) 'S/C vector', SCx, SCz, SCm



        WRITE(6,*) 'emm at I J', emm(I,J)
        WRITE(6,*) 'emm at I J', emm(I,J)
        WRITE(6,*) (Nx*SCx + Nz*SCz)/(Nm*SCm)
        WRITE(6,*) (Nx*SCx)/(Nm*SCm), (Nz*SCz)/(Nm*SCm)
      endif
      endif

c      Nx = 0
c      Ny = Nz  
c      Nm = SQRT(Nx*Nx + Ny*Ny + Nz*Nz)
c     End Testing

c     normal vector dotted with spacecraft vector (SCy is zero)
      emm(I,J) = ACOS((Nx * SCx + Nz * SCz) / (Nm*SCm))

c     Solar vector and object vector
c      write (*,*) "Sun (C): ", sunLat*180/Pi, sunLon*180/Pi, 
c     +      "b", C*180/Pi
c      write (*,*) "Object (B): ", objLat*180/Pi, objLon*180/Pi, 
c     +      "c", B*180/Pi


c      write (*,*) "cosA", cosA
c      write (*,*) "cosA angle ", acos(cosA)*180/Pi

      ENDDO
      ENDDO

c ########################################
c       Calculate flux
c ########################################

      WRITE(6,*) '# Calculating flux and temperature'
      DO J=Jstart,Jend
      DO I=Istart,Iend
      if (slope(I,J) .lt. 89.0/180 * Pi) then
         fluxA (I,J) = cos(slope(I,J)) * 1370 / solarDist / solarDist *
     +    (1-AL0(I,J)) 
      if (fluxA (I,J) .lt. 0) fluxA (I,J) = 0
      endif

      ENDDO
      ENDDO


c ########################################
c Average values
c ########################################
      avgInc = 0
      avgAlb = 0
      cnt = 0
      DO J=Jstart,Jend
        DO I=Istart,Iend
          cnt = cnt+1
          avgInc = avgInc + slope(I,J);
          avgAlb = avgAlb + AL0(I,J) 

        ENDDO
c        write (*,*)
      ENDDO

      avgAlb = avgAlb/cnt
      avgInc = avgInc/cnt
      write (*,"(A,f6.2)") " # Mean Albedo: ", avgAlb
      write (*,*) "# Mean Incidence angle: ", avgInc*180/Pi
      avgTemp = sqrt (sqrt (cos(avgInc) * 1370 /solarDist/solarDist
     +          * (1-avgAlb) / 4.56e-8))

c ########################################
c Print Grid
c ########################################
      DO J=Jstart,Jend
        DO I=Istart,Iend
          inc(I,J) = inc(I,J) * 180 / 3.1415926
          emm(I,J) = emm(I,J) * 180 / 3.1415926

        ENDDO
      ENDDO
      outName=MAP0//'-flux.txt'
c      CALL write_array (fluxA, NTMP,QSZ, outName)
      outName=MAP0//'-i.txt'
      CALL write_array (inc, NTMP,QSZ, outName)
      outName=MAP0//'-e.txt'
      CALL write_array (emm, NTMP,QSZ, outName)
      outName=MAP0//'-a.txt'
      CALL write_array (AL0, NTMP,QSZ, outName)
      outName=MAP0//'-evflipX.txt'
      CALL write_array (MevflipX, NTMP,QSZ, outName)

      STOP
      END

c   ................................................
      SUBROUTINE planck(currTemp, spec)
c   ................................................
      IMPLICIT NONE
      
      REAL               currTemp
      REAL*4             spec(30000)
      INTEGER            wave
      REAL*4             h, c, boltz, flux, emm, st_boltz, tWave
      REAL*4             top, bottom, eTerm

      h = 6.626069e-34
      c = 299792458
      boltz = 1.380648e-23
      flux = 1370
      emm = 0.94
      st_boltz = 5.67e-8

      DO wave=1,30000
        tWave = wave/1.0e9
        top = 2.0 * h * c * c
        eTerm = h * c / tWave / boltz / currTemp
        bottom = tWave**5 * (exp (eTerm) - 1)
       
        flux = top/bottom
        spec(wave) = spec(wave) + flux

c        write (*,*) wave/1000.0, flux
      enddo

  
      return
      END

c   ................................................
      SUBROUTINE READ_MAP(LMRKFILE,NTMP,QSZ,SCALE,V,UX,UY,UZ,HT,ALB)
c   ................................................

      IMPLICIT NONE
      
      INTEGER               NTMP
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      INTEGER               K
      INTEGER               K0
      INTEGER               IX(24)
      INTEGER               JX(24)
      INTEGER               NREC

      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      UX(3)
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)

      REAL*4                HT(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                ALB(-NTMP:NTMP,-NTMP:NTMP)

      CHARACTER*72          BLINE
      CHARACTER*72          LMRKFILE

      CHARACTER*1           CH1
      CHARACTER*2           CH2, CH2F
      CHARACTER*4           CH4, CH4F
      INTEGER*2             IX2
      REAL*4                RL4
      EQUIVALENCE          (IX2,CH2)
      EQUIVALENCE          (RL4,CH4)

      CHARACTER*2           C2
      INTEGER*2             I2
      EQUIVALENCE          (I2,C2)
      LOGICAL               LFLAG

      c2='69'
      LFLAG=.TRUE.
      if(i2.eq.13881) LFLAG=.FALSE.

      OPEN(UNIT=20, FILE=LMRKFILE, ACCESS='DIRECT',
     .     RECL=72, status='OLD')

c       Read is the scale fo the bigmap
        READ(20,REC=1) BLINE
        CH4f=BLINE(7:10)
        call flip(4,lflag,ch4f,ch4)
        SCALE=RL4

c       Read is the size of the bigmap
        QSZ=ICHAR(BLINE(11:11))
     .     +ICHAR(BLINE(12:12))*256

c       Loop over K four times (for values for V, and U[XYZ])
c       Load up V, UX, UY and UZ
        DO K=1,3
          CH4f=BLINE(12+4*K:15+4*K)
          call flip(4,lflag,ch4f,ch4)
          V(K)=RL4
          CH4f=BLINE(24+4*K:27+4*K)
          call flip(4,lflag,ch4f,ch4)
          UX(K)=RL4
          CH4f=BLINE(36+4*K:39+4*K)
          call flip(4,lflag,ch4f,ch4)
          UY(K)=RL4
          CH4f=BLINE(48+4*K:51+4*K)
          call flip(4,lflag,ch4f,ch4)
          UZ(K)=RL4
        ENDDO


c       Load the H Scale
        CH4f=BLINE(64:67)
        call flip(4,lflag,ch4f,ch4)
        HSCALE=RL4

c       Clear the memory
        DO I=-QSZ,QSZ
          DO J=-QSZ,QSZ
            ALB(I,J)=0
            HT(I,J)=0
          ENDDO
        ENDDO
        
c       Read stuff in, 24 bytes at a time
        NREC=1        
        K=0
        DO J=-QSZ,QSZ
          DO I=-QSZ,QSZ
            K=K+1
            IX(K)=I
            JX(K)=J
            IF(K.EQ.24) THEN
              NREC=NREC+1
              READ(20,REC=NREC) BLINE
              DO K=1,24
                CH1=BLINE(3*K:3*K)
                IF(ICHAR(CH1).NE.0) THEN
                  CH2f=BLINE(3*K-2:3*K-1)
                  call flip(2,lflag,ch2f,ch2)
                  HT(IX(K),JX(K))=HSCALE*IX2
                  ALB(IX(K),JX(K))=.01*ICHAR(CH1)
                ENDIF
              ENDDO
              K=0
            ENDIF
          enddo
        enddo

c       Read the remainder of a 24 byte line
        IF(K.NE.0) THEN
          K0=K
          NREC=NREC+1
          READ(20,REC=NREC) BLINE
          DO K=1,K0
            CH1=BLINE(3*K:3*K)
            CH2f=BLINE(3*K-2:3*K-1)
            call flip(2,lflag,ch2f,ch2)
            HT(IX(K),JX(K))=HSCALE*IX2
            ALB(IX(K),JX(K))=.01*ICHAR(CH1)
          ENDDO
        ENDIF

      CLOSE(UNIT=20)
      
      RETURN
      END     


c   ..................................................
      subroutine flip(n,lflag,ch1,ch2)
c   ..................................................

      integer*4        n, i
      character*(*)    ch1, ch2
      logical          lflag

      if(lflag) then
        do i=1,n
          ch2(i:i)=ch1(n-i+1:n-i+1)
        enddo
      else
        ch2=ch1
      endif
      
      return
      end


c   ................................................
      SUBROUTINE write_array(temp0, NTMP, QSZ, fName)

c   ................................................

      IMPLICIT NONE
      
      INTEGER               NTMP

      INTEGER               NREC
      INTEGER               LENSIZE
      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      REAL                temp0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                sunLon


      CHARACTER*72          fName

      WRITE(6,*) '# Output: ', fName
      OPEN(UNIT=10, FILE=fName)
      DO J=-qsz,qsz
        DO I=-qsz,qsz
          write (10,"(f15.4)", advance='no') temp0(I,J)
        enddo
        write (10,*)
      enddo
      CLOSE(UNIT=10)

      RETURN
      END     

c   ................................................
      SUBROUTINE read_array(temp0, NTMP, QSZ, fName)

c   ................................................

      IMPLICIT NONE

      INTEGER               NTMP

      INTEGER               NREC
      INTEGER               LENSIZE
      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      REAL                temp0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                sunLon


      CHARACTER*72          fName

      OPEN(UNIT=10, FILE=fName)
      DO J=-qsz,qsz
        DO I=-qsz,qsz
          write (10,"(f15.4)", advance='no') temp0(I,J)
        enddo
        write (10,*)
      enddo
      CLOSE(UNIT=10)

      RETURN
      END


 
