c---------------------------------------------------------------------
      program merg
c     programmer:      Bob Joyce
c     date created:    9/21/2000 NOTE! for sgi109
c     purpose:         daily and syn period geo sat overlap calb study
c     Date: 10/10/2001 Purpose is to determine daily geo-stationary
c                      inter-satellite calibration rates.  Runs once
c                      a day and updates the 21-day-back-filter file
c                      which is an input file for the 4-km IR merge
c                      file which does inter-satellite calibration
c                      operationally
c---------------------------------------------------------------------

      parameter (I5MAX=720, J5MAX=240, IP5MAX=6, IP5LEV=2,
     2 ITMPBINMX=34, ZACTLMX=26.5, ZALMBMN=40.0, ISYNMX=8,
     2 ZALMBMX=84.0, XMISSNG=-9999., ISTDYMX=5,
     3 TMPCTLMN=160., IPMX=3, MINOBS=10, 
     4 IONE=1, TMPINCR=5.0, ZAINCR=2.0, HALFSYN=1.5,
     5 SYNPER=3.0, ISDAY=28, ISYEAR=1999, 
     & ISTMN=2, ISTDY1=26)

c      Control TEMP binnned from 160 to 330 degrees every 5

      real*4 xsatlon(5), stat(ITMPBINMX,ISYNMX,ISTDYMX,IPMX),
     1 arms(ISTDYMX), abs1(ISTDYMX), stdycnt(ISTDYMX),
     & statdy(ITMPBINMX,ISTDYMX,IPMX)

      character*110  mergein, mergeinc, errorout, errordyout
      character*3    month(12), mnth

      integer*2 iprod(I5MAX,J5MAX,IP5LEV), yrdays(13), lyrdays(13)

      
      logical there

      character yyyymmddhh*10

c---------------------------------------------------------------------
c updated H-8 value (141.0) and Met-8 value (41.5) --S.Wu (2018.04.13)
c---------------------------------------------------------------------
c GOES-West (GOES-17) :: 137.0W (=223.0) --S.Wu (2019.06.10)
c---------------------------------------------------------------------
      data xsatlon / 141.0, 223.0, 285.0, 0.0, 41.5 /

c     data xsatlon / 141.0, 225.0, 285.0, 0.0, 41.5 /
c     data xsatlon / 140.0, 225.0, 285.0, 0.0, 57.0 /
c     data xsatlon / 155.0, 225.0, 285.0, 0.0, 63.0 /

      data month/"Jan","Feb","Mar","Apr","May","Jun",
     1           "Jul","Aug","Sep","Oct","Nov","Dec"/
      data  yrdays/0,31,59,90,120,151,181,212,243,273,304,334,365/
      data lyrdays/0,31,60,91,121,152,182,213,244,274,305,335,366/


c  Get high-level path name from appropriate ".alias" file (dictated by the Unix script that runs this pgm)

      iwrec = 0
      icrec = 0

c    Make call to determine Julian Day of dataset first day

      call date2doy(ISYEAR,ISTMN,ISTDY1,istjday)

c============================================= swu
c      open(29, file='../input/calb_datefile')
c      read (29,'(4x,a3,i3,i3,11x,i4)') mnth, idy, ihr, iyr
c      do 107 imth=1,12
c 107   if (mnth.eq.month(imth)) goto 108
c 108  continue
c
      if( iargc().lt.1 ) then
        write(*, *)
        write(*, *) '###################################'
        write(*, *) '#                                 #'
        write(*, *) '#  parameter <yyyymmddhh> needed  #'
        write(*, *) '#                                 #'
        write(*, *) '###################################'
        write(*, *)
        stop
      endif
      call getarg(1, yyyymmddhh)
      read(yyyymmddhh, '(i4.4,3i2.2)') iyr, imth, idy, ihr
c
c== swu

      call date2doy(iyr, imth, idy, ijday)

      do iyr1 = ISYEAR, iyr
 
       if (mod(iyr1,4).eq.0.and.mod(iyr1,100).ne.0
     &              .or.mod(iyr1,1000).eq.0) then
        lyear = IONE
       else
        lyear = 0
       endif
 
       if (iyr1.ne.iyr) then
        if (lyear.eq.0) then
         icrec = icrec + yrdays(13)
        else
         icrec = icrec + lyrdays(13)
        endif
       else
        if (lyear.eq.0) then
         icrec = icrec + ijday - istjday + 1
        else
         icrec = icrec + ijday - istjday + 1
        endif
       endif
 
       write (*,*) ' iyr1, icrec ', iyr1, icrec

      enddo

      write (*,*) ' icrec ', icrec
c     if (icrec.eq.icrec) stop 979

      iyrmndy = iyr*10000+imth*100+idy

c     Open output statistics file
 
      errorout='../library/' //
     1    'geo-sat_daily-syn-ir-calb-diff_02261999-current'
      errordyout='../library/' //
     1    'geo-sat_daily-ir-calb-diff_02261999-current'
 
      open (unit=60, file = errorout, access = 'direct',
     1        form='unformatted', recl=34*8*5*3*4)

      open (unit=62, file = errordyout, access = 'direct',
     1        form='unformatted', recl=34*5*3*4)
 
 
c    Read and write current daily + synoptic calibration file
     
      do irrec = IONE, icrec - 1
c      write (*,*) ' irrec ', irrec
       read (60,rec=irrec) stat
       write (60,rec=irrec) stat
       read (62,rec=irrec) statdy
       write (62,rec=irrec) statdy
      enddo
 
c    Process current day
 
      write (*,*) ' imth, idy ', imth, idy

c     Initialize arrays

      do 11 istdy = 1, ISTDYMX
       stdycnt(istdy) = 0.
       arms(istdy) = 0.
       abs1(istdy)  = 0.
       do 11 itbn = 1, ITMPBINMX
        do 11 iprd = 1, IPMX
         statdy(itbn,istdy,iprd) = 0.
         do 11 isyn = 1, ISYNMX
 11   stat(itbn,isyn,istdy,iprd) = 0.

c     Do hourly loop

c     do 20 ihr = 0, 17, 1
      do 20 ihr = 0, 23, 1
 
c      write (6,'(/,'' Proc '',3i9/)') iyrmndy, ihr, ihlfhr
  
       if (icrec.ge.784) then
c--------------->123456789|123456789|123456789|123456789|123456789|123<--
       mergeinc='../output/uncalb/half_deg/globl_YYYY0M0D0H_0.5-deg_i2'
        else
       mergeinc='../output/merged/half_deg/globl_YYYY0M0D0H_0.5-deg_i2'
       endif

       write(mergeinc(33:42), '(i8,i2.2)') iyrmndy, ihr

       write (*,*) ' mergeinc ', mergeinc

c     Open input 0.5 degree product data set

       do 25 ihlfhr = 0, 30, 30

      if (ihlfhr.eq.0) then

       xhlfhr = 0.0

c      write (*,*) mergeinc

       inquire (file=mergeinc,exist=there)

       if (.not.there) goto 20
 
       open (unit=10, file = mergeinc, access = 'direct',
     1              form='unformatted', recl=720*240*2*2)

       read (10, rec=1) iprod

      else

       xhlfhr = 0.5
       read (10, rec=2) iprod

      endif

c     Determine which synoptic frame hour belongs

      xhour = ihr + xhlfhr
      isyn  = (xhour + HALFSYN)/SYNPER + 1
      if (isyn.gt.ISYNMX) isyn  = 1

c     Determine calibration differences in temp

      do 500 j5 = 1, J5MAX

       x5lat =  (1-j5)*0.5 + 59.75
       do 500 i5 = 1, I5MAX
        x5lon = (i5-1)*0.5 + 0.25
        tmpe = iprod(i5,j5,2)/10.
        tmpo = iprod(i5,j5,1)/10.

c        write(*,*) tmpe, tmpo
 
        if(tmpe.GT.0. .and. tmpo.GT.0.)then

c     Determine Odd Satellite ID by longitude location
c     Modify to account for GOES-9 at 155E replacing 
c     GMS-5 at 140E

         if (i5.le.140) then
          isato = 5
c        elseif (i5.ge.141.and.i5.le.420) then
         elseif (i5.ge.141.and.i5.le.440) then
          isato = 1
         elseif (i5.ge.431.and.i5.le.710) then
          isato = 3
         endif

c     Determine Even Satellite ID by longitude location

         if (i5.le.140) then
          isate = 4
         elseif (i5.ge.141.and.i5.le.266) then
          isate = 5
         elseif (i5.ge.311.and.i5.le.585) then
          isate = 2
         elseif (i5.ge.586) then
          isate = 4
         endif
 
         ctmp = (tmpe+tmpo)/2.
         itbn = (ctmp-TMPCTLMN)/5. + 1.
         satpsne = xsatlon(isate)
         satpsno = xsatlon(isato)
 
c     Determine Study site (1-5)

         if (isato.eq.1.and.isate.eq.2) then
          istdy=1
          diff = tmpo-tmpe
         elseif (isato.eq.3.and.isate.eq.2) then
          istdy=2
          diff = tmpe-tmpo
         elseif (isato.eq.3.and.isate.eq.4) then
          istdy=3
          diff = tmpo-tmpe
         elseif (isato.eq.5.and.isate.eq.4) then
          istdy=4
          diff = tmpe-tmpo
         elseif (isato.eq.1.and.isate.eq.5) then
          istdy=5
          diff = tmpe-tmpo
         endif

         call satang(x5lat,x5lon,satpsne,szae,angndr) 
         call satang(x5lat,x5lon,satpsno,szao,angndr) 

         szadiff = szae-szao

c       write(*,*) itbn, isyn, istdy, szadiff

         if (szadiff.le.5.and.szadiff.ge.-5.and.x5lat.le.40.and.
     &                                            x5lat.ge.-40) then

          stat(itbn,isyn,istdy,3)=stat(itbn,isyn,istdy,3) + 1
 
          statdy(itbn,istdy,3)=statdy(itbn,istdy,3) + 1
          stdycnt(istdy)    = stdycnt(istdy) + 1

c       write(*,*) '----------1'

          diffsqrd = diff*diff
          stat(itbn,isyn,istdy,1) = stat(itbn,isyn,istdy,1) + diff
c       write(*,*) '----------2'

          statdy(itbn,istdy,1) = statdy(itbn,istdy,1) + diff
          abs1(istdy) = abs1(istdy) + diff
          stat(itbn,isyn,istdy,2) = stat(itbn,isyn,istdy,2) + diffsqrd
          statdy(itbn,istdy,2) = statdy(itbn,istdy,2) + diffsqrd
          arms(istdy) = arms(istdy) + diffsqrd
         endif
c       write(*,*) 'HERE==========================='

        endif

 500  continue

 25   continue

      close (10)
 


 20   continue

      do 700 istdy = IONE, ISTDYMX
       do 707 itbn = IONE, ITMPBINMX
        if (statdy(itbn,istdy,3).ge.MINOBS) then
         statdy(itbn,istdy,1)=statdy(itbn,istdy,1)/
     1                           statdy(itbn,istdy,3)
         statdy(itbn,istdy,2)=sqrt(statdy(itbn,istdy,2)/
     1                           statdy(itbn,istdy,3))
        else
         statdy(itbn,istdy,1)=XMISSNG
         statdy(itbn,istdy,2)=XMISSNG
         statdy(itbn,istdy,3)=XMISSNG
        endif
        do 707 isyn = IONE, ISYNMX
         if (stat(itbn,isyn,istdy,3).ge.MINOBS) then
          stat(itbn,isyn,istdy,1)=stat(itbn,isyn,istdy,1)/
     1                            stat(itbn,isyn,istdy,3)
          stat(itbn,isyn,istdy,2)=sqrt(stat(itbn,isyn,istdy,2)/
     1                            stat(itbn,isyn,istdy,3))
         else
          stat(itbn,isyn,istdy,1)=XMISSNG
          stat(itbn,isyn,istdy,2)=XMISSNG
          stat(itbn,isyn,istdy,3)=XMISSNG
         endif
 707  continue

      arms(istdy) = sqrt(arms(istdy)/stdycnt(istdy))
      abs1(istdy)  = abs1(istdy)/stdycnt(istdy)

 700  continue

      write (60,rec=icrec) stat
      write (62,rec=icrec) statdy

 15   continue

      call backfilter(icrec)

      call meanfilter(icrec)
c      stop 

      end

 
      subroutine backfilter(icrec)
 
c     programmer:      Bob Joyce
c     date created:    December 4, 2000
c     purpose:         day a 21 day backwards running mean


      real XMISSNG

      integer DAYS

      parameter (ITMPBINMX=34, IPMX=3, DAYMN=7, XMISSNG=-9999.0, 
     & ISTDYMX=5, IOUTMX=6, DAYS=21, IONE=1, TMPINCR=5.0)

c    Control TEMP binnned from 160 to 330 degrees every 5

      real*4 statdy(ITMPBINMX,ISTDYMX,IPMX), 
     &  statoutcnt(ITMPBINMX,IOUTMX,IPMX), studycnt(ITMPBINMX,IPMX),
     &  statout(ITMPBINMX,IOUTMX,IPMX)
     
      character*150 error21back

      irec = 0

      write (*,*) ' in backfilter '

c    Open backfilter output statistics file
c    Currently both files go back to Feb 26, 1999

      error21back='../library/' //
     1    'geo-sat_calb-diff_02261999-current-21day-back-filter'

      open (unit=70, file = error21back, access = 'direct',
     1            form='unformatted',recl=34*6*3*4)

c    Do Daily loop

      do iyrmndy = IONE, icrec

c    Initialize 21-day backfilter averaging arrays

       do itbn = IONE, ITMPBINMX
        do iprd = IONE, IPMX
         studycnt(itbn,iprd) = 0.
         do istdy = IONE, IOUTMX
          statout(itbn,istdy,iprd) = 0.
          statoutcnt(itbn,istdy,iprd) = 0.
         enddo
        enddo
       enddo

       iareccnt = 0

c     Read in 21 days of calibration data for each backfilter day

       do iarec = iyrmndy - DAYS, iyrmndy - IONE

        if (iarec.ge.IONE.and.iarec.le.icrec) then
 
         read (62, rec = iarec) statdy

         iareccnt = iareccnt + IONE

c        write (*,*) ' iyrmndy icrec iarec ', iyrmndy,icrec,iarec

         do itbn = IONE, ITMPBINMX
          do iprd = IONE, IPMX
           do istdy = IONE, ISTDYMX
            if (statdy(itbn,istdy,iprd).gt.XMISSNG) then
             statout(itbn,istdy,iprd) =
     &       statout(itbn,istdy,iprd)+statdy(itbn,istdy,iprd)
             statoutcnt(itbn,istdy,iprd) =
     &       statoutcnt(itbn,istdy,iprd) + 1.
             studycnt(itbn,iprd) = studycnt(itbn,iprd) + 1.
            endif
           enddo
          enddo
         enddo
        endif
       enddo

       do itbn = IONE, ITMPBINMX
        do iprd = IONE, IPMX
         istdycnt = 0
         do istdy = IONE, ISTDYMX
          if (statoutcnt(itbn,istdy,iprd).ge.DAYMN) then
           statout(itbn,istdy,iprd) =
     &     statout(itbn,istdy,iprd)/statoutcnt(itbn,istdy,iprd)
           istdycnt = istdycnt + IONE
           statout(itbn,6,iprd) =
     &      statout(itbn,6,iprd) + statout(itbn,istdy,iprd)
          else
           statout(itbn,istdy,iprd) = XMISSNG
          endif
         enddo
         if (istdycnt.ne.ISTDYMX) statout(itbn,6,iprd) = XMISSNG
        enddo
       enddo

c    Write out whole record

       write (70, rec = iyrmndy) statout

      enddo

      return

      end
 

      subroutine meanfilter(icrec)

c     programmer:      Bob Joyce
c     date created:    December 4, 2000
c     purpose:         create a 31 day running mean calibration filter

      real XMISSNG

      parameter (ITMPBINMX=34, IPMX=3, DAYMN=7, XMISSNG=-9999.0, 
     & ISTDYMX=5, IOUTMX=6, MDAYS=15, IONE=1, TMPINCR=5.0)

c    Control TEMP binnned from 160 to 330 degrees every 5

      real*4 statdy(ITMPBINMX,ISTDYMX,IPMX), 
     &  statoutcnt(ITMPBINMX,IOUTMX,IPMX), studycnt(ITMPBINMX,IPMX),
     &  statout(ITMPBINMX,IOUTMX,IPMX)
     
      character*150 error31mean

      irec = 0

      write (*,*) ' in meanfilter '

c    Open backfilter output statistics file
c    Currently both files go back to Jan 28, 2000

      error31mean='../library/' //
     1    'geo-sat_calb-diff_02261999-current-31day-mean-filter'

      open (unit=80, file = error31mean, access = 'direct',
     1                   form='unformatted', recl=34*6*3*4)

c    Do Daily loop

      do iyrmndy = IONE, icrec

c    Initialize 31-day mean_filter averaging arrays

       do itbn = IONE, ITMPBINMX
        do iprd = IONE, IPMX
         studycnt(itbn,iprd) = 0.
         do istdy = IONE, IOUTMX
          statout(itbn,istdy,iprd) = 0.
          statoutcnt(itbn,istdy,iprd) = 0.
         enddo
        enddo
       enddo

       iareccnt = 0

c     Read in 21 days of calibration data for each backfilter day

       do iarec = iyrmndy - MDAYS, iyrmndy + MDAYS

        if (iarec.ge.IONE.and.iarec.le.icrec) then

         read (62, rec = iarec) statdy

         iareccnt = iareccnt + IONE

c        write (*,*) ' iyrmndy icrec iarec ', iyrmndy,icrec,iarec

         do itbn = IONE, ITMPBINMX
          do iprd = IONE, IPMX
           do istdy = IONE, ISTDYMX
            if (statdy(itbn,istdy,iprd).gt.XMISSNG) then
             statout(itbn,istdy,iprd) =
     &       statout(itbn,istdy,iprd)+statdy(itbn,istdy,iprd)
             statoutcnt(itbn,istdy,iprd) =
     &       statoutcnt(itbn,istdy,iprd) + 1.
             studycnt(itbn,iprd) = studycnt(itbn,iprd) + 1.
            endif
           enddo
          enddo
         enddo
        endif
       enddo

       do itbn = IONE, ITMPBINMX
        do iprd = IONE, IPMX
         istdycnt = 0
         do istdy = IONE, ISTDYMX
          if (statoutcnt(itbn,istdy,iprd).ge.DAYMN) then
           statout(itbn,istdy,iprd) =
     &     statout(itbn,istdy,iprd)/statoutcnt(itbn,istdy,iprd)
           istdycnt = istdycnt + IONE
           statout(itbn,6,iprd) =
     &      statout(itbn,6,iprd) + statout(itbn,istdy,iprd)
          else
           statout(itbn,istdy,iprd) = XMISSNG
          endif
         enddo
         if (istdycnt.ne.ISTDYMX) statout(itbn,6,iprd) = XMISSNG
        enddo
       enddo

c    Write out whole record

       write (80, rec = iyrmndy) statout

      enddo

      return

      end

 

      SUBROUTINE satang(XLAT,XLON,SUBLON,ANGZEN,ANGNDR)

      DATA TWOPI/6.28318/, R/6371./, H/35680./

      DEGRAD=360./TWOPI

      ZLON=XLON
      IF(ZLON.LT.0.) ZLON=360.+ZLON

      ZSUBLN=SUBLON
      IF(ZSUBLN.LT.0.) ZSUBLN=360.+ZSUBLN

      DIFLON=ABS(ZLON-ZSUBLN)/DEGRAD
      DIFLAT=ABS(XLAT)/DEGRAD

      COSPSI=COS(DIFLON)*COS(DIFLAT)
      PSI=ACOS(COSPSI)

      SSQR=R*R + (R+H)*(R+H) - 2*R*(R+H)*COSPSI
      S=SQRT(SSQR)

      SINNDR=(R/S)*SIN(PSI)
      SINZEN=((R+H)/R)*SINNDR

      ANGZEN=ASIN(SINZEN)*DEGRAD
      ANGNDR=ASIN(SINNDR)*DEGRAD

      RETURN
      END
 

      SUBROUTINE date2doy(NYY,NMM,NDD,NCEN)
C
c         computes day-of-year from year(NYY), month(NMM), and day(NDD)
c
      INTEGER*2  NTAB(12)
C
      DATA  NTAB/0,31,59,90,120,151,181,212,243,273,304,334/
C
      LP = 0
      IF (MOD(NYY,4).EQ.0 )  LP = 1
      NLP = NYY / 4
      IF (LP.EQ.1 .AND. NLP.GT.0)  NLP = NLP - 1

      IF (LP.EQ.1 .AND. NMM.LT.3)  LP  = 0
      JDAY = NTAB(NMM) + NDD + LP
      NCEN =  JDAY
      RETURN
      END
