program create_monthly_zonal_csrbs

c     Code to create monthly zonal means of CSRBs from monthly 1 deg.
c     by 1 deg. CSRBS.

c *** CAUTION: The 1-degree radiance and BT biases ('monthly_rad_diff',
c     'monthly_bt_diff' below) have been smoothed by 
c     averaging 3x3 groups of values.  However, 'ndays' refers to the
c     original number of days for each grid box, and 'monthly_sdev' is
c     also the original standard deviation of gridded values.

      implicit none

c----------------------------------------------------------------------

      integer nlon_bins
      parameter (nlon_bins = 360)

      integer nlat_bins
      parameter (nlat_bins = 180)

      integer nbands
      parameter (nbands = 8)

      integer nscns
      parameter (nscns = 3)

      integer nbins
      parameter(nbins = nlat_bins * nbands * nscns)

c----------------------------------------------------------------------

c     Input file name that contains 1-degree gridded monthly mean CSRBs.
      character*120 Input_File, parm, Output_File

c     Input monthly mean data arrays.
      real*4 ndays(nlon_bins, nlat_bins, nbands, nscns)
      real*4 monthly_rad_diff(nlon_bins, nlat_bins, nbands, nscns)
      real*4 monthly_bt_diff(nlon_bins, nlat_bins, nbands, nscns)
      real*4 monthly_sdev(nlon_bins, nlat_bins, nbands, nscns)

c----------------------------------------------------------------------

c     Output data arrays.
      real*4 zonal_csrb_rads(nlat_bins, nbands, nscns)
      real*4 smoothed_zonal_csrb_rads(nlat_bins, nbands, nscns)
      real*4 zonal_csrb_bts(nlat_bins, nbands, nscns)
      real*4 zonal_csrb_sdev(nlat_bins, nbands, nscns)
     
c----------------------------------------------------------------------

      integer i, j, k, m, jj, jk
      real*4 n, sum_rads, sum_bts, sum_sdev, rlat, tnum, tsum
      
c----------------------------------------------------------------------

c     Get input and output file names.

      call getarg(1,parm)
      read(parm,'(a120)') Input_File
      write(*,'(2x,'' Input File is: '',A70)') Input_File

      call getarg(2,parm)
      read(parm,'(a120)') Output_File
      write(*,'(2x,'' Output File is: '',A70)') Output_File

c----------------------------------------------------------------------

c     Open, read, and close input file.

      open(30, file=Input_File, form='unformatted', access='sequential')
      read(30) ndays, monthly_rad_diff, monthly_bt_diff, monthly_sdev
      close(30)

c----------------------------------------------------------------------

c     Loop through the 1-degree gridded CSRBs.

      do m = 1, nscns
        do k = 1, nbands
          do j = 1, nlat_bins

            n = 0.0
            sum_rads = 0.0
            sum_bts = 0.0
            sum_sdev = 0.0

            do i = 1, nlon_bins

              if(monthly_rad_diff(i,j,k,m) .gt. -900.0) then
                n = n + 1.0
                sum_rads = sum_rads + monthly_rad_diff(i,j,k,m)
                sum_bts = sum_bts + monthly_bt_diff(i,j,k,m)
                sum_sdev = sum_sdev + monthly_sdev(i,j,k,m)
              end if
c             if(j .eq. 6) then
c               write(*,'(i5,8f12.3)') i, ndays(i,j,k,m), monthly_rad_diff(i,j,k,m),
c    *             monthly_bt_diff(i,j,k,m), monthly_sdev(i,j,k,m), n, sum_rads,
c    *             sum_bts, sum_sdev
c             end if

            enddo

            if(n .gt. 0.0) then
              zonal_csrb_rads(j,k,m) = sum_rads / n
              zonal_csrb_bts(j,k,m) = sum_bts / n
              zonal_csrb_sdev(j,k,m) = sum_sdev / n
            else
              zonal_csrb_rads(j,k,m) = -999.0
              zonal_csrb_bts(j,k,m) = -999.0
              zonal_csrb_sdev(j,k,m) = -999.0
            end if

c           write(*,'(3i5,7f12.3)') m, k, j, n, sum_rads, sum_bts, sum_sdev,
c    *         zonal_csrb_rads(j,k,m), zonal_csrb_bts(j,k,m), zonal_csrb_sdev(j,k,m)

          enddo

c----------------------------------------------------------------------

c         Get 5-zone runnimg mean. Values at end points (zones 1, 2, 
c         179, 180) are copied from zones 3 and 178.

          do jk = 3, nlat_bins-2

            tsum = 0.0
            tnum = 0.0

c           Form running means.
            do jj = jk-2, jk+2

              if(zonal_csrb_rads(jj,k,m) .gt. -900.0) then
                tsum = tsum + zonal_csrb_rads(jj,k,m)
                tnum = tnum + 1.0
              end if

            enddo

            if(tnum .gt. 0.0) then
              smoothed_zonal_csrb_rads(jk,k,m) = tsum / tnum
            else
              smoothed_zonal_csrb_rads(jk,k,m) = -999.0
            end if

c           Fill in bands 1, 2, 179, and 180.
            if(jk .eq. 3) then
              smoothed_zonal_csrb_rads(1,k,m) = smoothed_zonal_csrb_rads(3,k,m)
              smoothed_zonal_csrb_rads(2,k,m) = smoothed_zonal_csrb_rads(3,k,m)
            else if(jk .eq. 178) then
              smoothed_zonal_csrb_rads(179,k,m) = smoothed_zonal_csrb_rads(178,k,m)
              smoothed_zonal_csrb_rads(180,k,m) = smoothed_zonal_csrb_rads(178,k,m)
            end if

c           write(*,*)'Original biases: '
c           write(*,'(3i5,8f12.3)') m, k, jk, zonal_csrb_rads(jk,k,m)
c           write(*,*)'Smoothed biases: '
c           write(*,'(i5,10f12.3)') jk, tnum, tsum, smoothed_zonal_csrb_rads(jk,k,m)

          enddo

c----------------------------------------------------------------------

        enddo
      enddo

c     m =1 (ocean), m=2 (AN land), m=3 (DN land)
      do m = 1, 1
        do j = 1, nlat_bins

          rlat = (90.0 - (j*1.0)) + 0.5
c         write(*,'(10f12.3)') j*1.0,rlat,(zonal_csrb_bts(j,k,m),k=1,8)
          write(*,'(10f12.3)') j*1.0,rlat,(smoothed_zonal_csrb_rads(j,k,m),k=1,8)
c         write(*,'(10f12.3)') j*1.0,rlat,(zonal_csrb_sdev(j,k,m),k=1,8)

        enddo
      enddo

c----------------------------------------------------------------------

c     Output CSRBs to binary file.

      open(20, file=Output_File, access='sequential', form='unformatted')
      write(20) smoothed_zonal_csrb_rads
      close(20)

c----------------------------------------------------------------------

     

      end