c  gfortran dumpMapLRO_fixErrors.f -O2
c
c 
c Modify to smooth over errors caused by rawMosaic.f
c Sept 10 2018

      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                startFill, endFill


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

      CHARACTER*6           MAP0
      CHARACTER*72          LMRKFILE
      CHARACTER*72          outName
      character*5          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        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. -qsz) then
         WRITE(6,*) 'I,J,baseEmm', I, J, baseEmm (I,J)
        ENDIF
       ENDDO
      READ (10,"(f12.5)"), baseEmm (Iend+1,J)
      ENDDO


      checkI = -195
      checkJ = -112
      WRITE(6,*) 'offset should be', baseEmm (checkI, checkJ)

c ########################################
c Determine and fix errors in emission
c ########################################
c     Wide field of view
      Istart = -qsz
      Iend = qsz
      Jstart = -qsz
      Jend = qsz

      fix = 0                                                           0 is no fix, 1 is fix

c J is line, I is pixel
       J=Jstart
      DO J=Jstart,Jend 
       DO I=Istart,Iend

c Check to see if we should modify the value
        IF (fix .eq. 0) then
         if (baseEmm (I,J) .eq. baseEmm (checkI, checkJ)) then
          if ((baseEmm (I-1,J) .ne. baseEmm (checkI, checkJ))
     +        .AND. (baseEmm (I-1,J) .ne. 0.0*PI/180)) then
c        identified start of error, record last "good" value
           counter=1
           startFill = baseEmm (I-1,J)
           IstartFill = I
           DO WHILE ((baseEmm(I+counter,J).eq.baseEmm(checkI,checkJ))
     +               .AND.(counter.lt.61))
            counter=counter+1
           ENDDO                                                         DO WHILE
           IF (counter .le. 59) then
            endFill = baseEmm (I+counter,J)
            fix=1
            WRITE (6,*) 's, e, c, I, J', startFill, endFill, counter, I,
     +   J
           ENDIF                                                         if endFill
          endif                                                          if baseEmm (I-1,J)
         endif                                                           if baseEmm (I,J)
        ENDIF                                                           end search for fix
        IF (fix .eq. 1) then
         diff=I-IstartFill
         baseEmm (I,J) = startFill +((endFill-startFill)/(counter+1))
     +                   *(ABS(diff)+1)
         WRITE (6,*) ABS(diff)
        ENDIF
        IF ((IstartFill + counter -1) .eq. I) then
         fix = 0                                                         Get out of fixing loop
        ENDIF

       ENDDO                                                            Istart to Iend
      ENDDO




c ########################################
c Print Grid
c ########################################

      outName=MAP0//'-fixed-e.txt'
      CALL write_array (baseEmm, 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,"(f12.5)", 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


 
