program mpi_timavg_3d 
      implicit none
      include 'mpif.h'
!
!     Take time average for grib I format file
!
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real, allocatable :: data(:)
      real, allocatable :: fld(:)
      real, allocatable :: field(:)
      real  fac,xmin,xmax

      integer kpds(200),kgds(200),ipds(200),igds(200)
      integer jpds(200),jgds(200),idate(4)
      integer ibds(9)
      integer, allocatable :: ivr(:)
      integer i,j,k,jj,kk,n,ntimes,nhours,lupgb,lupgi,iret,lskip
      integer ijdim,nvar,kpds16,ndata,ntasks,nrank,lupgbo
      integer iy,im,id,ih,icnt,nt,ift,m1,m2,item

      character*300, cfile(2)
      character*300, allocatable :: cfiles(:,:)
      character*15  grib_out
      character     chead*81,chd*8,ctype*6,crun*6,cprog*9
      character pds(200)
      logical fcst_avrg
      integer iys, ims, ids, ihs, iye, ime, ide, ihe, iyb,dbgme
      integer ierr 

      logical*1, allocatable :: lbmi(:),lbmo(:)

      namelist/namin/ ntimes, nhours, lupgi, fcst_avrg                  &
     &,               iys, ims, ids, ihs, iye, ime, ide, ihe            &
     &,               dbgme
      namelist/files/ cfile

      call start()
!
      call mpi_init(iret)
      if (iret .ne. 0) then
        print*,' iret = ',iret,' from mpi_init' ;stop 1
      endif
      call mpi_comm_rank(MPI_COMM_WORLD,nrank,iret)
      if (iret .ne. 0) then
        print*,'MPI_COMM_WORLD, iret = ',MPI_COMM_WORLD,iret,' from mpi_comm_rank';stop 1
      endif
      call mpi_comm_size(MPI_COMM_WORLD,ntasks,iret)
      if (iret .ne. 0) then
        print*,' MPI_COMM_WORLD,iret = ',MPI_COMM_WORLD,iret,' from mpi_comm_size' ;stop 1
      endif
      print*,' task ',nrank,' of ',ntasks

      open(5,file="input1",status="old",form="formatted",access="sequential")!

!     dbgme = 15
      dbgme = 90
      fcst_avrg = .false.
      lupgb=11    
      lupgi=12 ! default to be replaced on input

      read(5,namin)
      if (nrank == 0) write(*,namin)
!     write(0,namin)

      allocate ( cfiles(2,ntimes),stat=iret)
      if (iret.ne.0) then
        print*, 'allocate ( cfiles(2,ntimes),stat=',iret,')'
        call mpi_abort(MPI_COMM_WORLD,1,ierr)
      endif

      do n = 1,ntimes
        read(5,files)
        cfiles(1,n) = cfile(1)
        if (lupgi .eq. 0 ) then
          if (nrank == 0) print*,n,lupgb,cfiles(1,n)
        else
          cfiles(2,n) = cfile(2)
          if (nrank == 0) print*,n,lupgb,cfiles(1,n),cfiles(2,n)
        endif
!                                                 open all the files
        call baopenr(lupgb,cfiles(1,n),iret)
        if (iret.ne.0) then
          write(6,*) 'open failed for n,lupgb,iret =', n,lupgb,iret,cfiles(1,n)
          call mpi_abort(MPI_COMM_WORLD,111,ierr)
        endif
        if (lupgi .ne. 0) then
          call baopenr(lupgi,cfiles(2,n),iret)
          if (iret.ne.0) then
            write(6,*) 'open failed for n,lupgi,iret =', n,lupgi,iret,cfiles(2,n)
            call mpi_abort(MPI_COMM_WORLD,111,ierr)
          endif
        endif
        lupgb = lupgb + 2
        if( lupgi .ne. 0 ) lupgi = lupgi + 2
      enddo

      lupgbo = 20 + ntimes + ntimes + nrank
      write(grib_out,'(''gribavg.'',i3.3)') nrank
      call baopenw(lupgbo,grib_out,iret)
      if (nrank == 0) print *, ' opening ',grib_out,' for write, iret,  &
     &                lupgbo =', iret,lupgbo

      ibds=0
      lskip=0

!     file names input from datacard files

!caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
! count the records, nvar
!caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa

      iret  = 0
      k     = 0
      nvar  = 0
      lupgb = 11
      if(lupgi.ne.0) lupgi=12
      n = 1
      do while (iret == 0 ) 
        j    = -1 - k
        ipds = -1
        igds = -1
        call getgbh(lupgb,lupgi,j,ipds,igds,j,ijdim,k,jpds,jgds,iret)
        if (iret .eq. 0) nvar = nvar + 1

!       print *,' nvar=',nvar,' nrank=',nrank,' ijdim=',ijdim
!
!  store the jpds5-7 from the k record for the first file 
!
        if (nvar == 1) then
!         if (nrank == 0) print*,'ijdim = ',ijdim
!         if (nrank == dbgme) print*,'ijdim = ',ijdim
          allocate ( field(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
          allocate ( fld(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
          allocate ( data(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
          allocate ( ivr(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
          allocate ( lbmi(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
          allocate ( lbmo(ijdim),stat=iret)
          if (iret.ne.0) call mpi_abort(MPI_COMM_WORLD,222,ierr)
        endif
      enddo
      if (nrank == 0) print*, nvar, ' records in ',cfiles(1,1)


!cc   take average after this loop, here ntimes is total files

      m1 = nint(float(nvar* nrank   ) / float(ntasks)) + 1
      m2 = nint(float(nvar*(nrank+1)) / float(ntasks))

      print*,' m1=',m1,' m2=',m2,' nrank=',nrank,' ntasks=',ntasks

      do kk = m1,m2                     ! first do loop for kk record
        k    = kk
        fld  = 0.0
        ivr  = 0
        nt   = 0

        lupgb = 11
        if(lupgi.ne.0) lupgi = 12
        lbmo = .true.
        fld  = 0.
        ivr  = 0.

        j = -k
        do n = 1, ntimes                   ! second do loop for n file

          ipds = -1
          igds = -1
          call getgb(lupgb,lupgi,ijdim,j,ipds,igds,ndata,lskip,jpds,jgds,lbmi,data,iret)
          if (iret .ne. 0 ) then
  	    print*,'getgb failed for lupgb,lupgi,k,n ' ,lupgb,lupgi,k,n,' iret=',iret
            call mpi_abort(MPI_COMM_WORLD,777,ierr)
	    cycle
          endif
          field = data 

          if (n.eq.1) then
!  
!  store the date-time info for file 1 and give to the average
!
            iy = jpds(8)
            im = jpds(9)
            id = jpds(10)
            ih = jpds(11)
            if (iy == 0) iy = 100
            iyb = (jpds(21)-1)*100 + iy
      print *,' iyb=',iyb,' jpds=',jpds(21),' iy=',iy
!           IF(jPDS(16).GE.2 .AND. jPDS(16).LE.5) THEN
!             IFT = jPDS(15)
!           ELSE
!             IFT = jPDS(14)
!           ENDIF
            print*,' Avg for k = ',k,jpds(5),jpds(6),jpds(7),iy,im,id,ih
            kpds16 = jpds(16)
          endif

          if (jpds(4) == 192) lbmo = lbmi
	
          nt = nt + 1
		
!cc  consideration of bitmap 

!         if (nrank == dbgme) print *,' jpds(4)=',jpds(4)
		
          do i = 1,ijdim
            if (jpds(4)  == 192) then
              if (i == 1 .and. nrank == 0)                              &
     &          print*, ' bit map considered at record ',k
              if (lbmo(i)) then
                ivr(i)  = ivr(i)  + 1
                fld(i)  = fld(i)  + field(i)
              endif
            else
              ivr(i)  = ivr(i)  + 1
              fld(i)  = fld(i)  + field(i)
            endif 
          enddo

!         if(jpds(5).eq.59) then
!	    print *, ' precip rate at 140,73 =', field(140*73)
!	    print *, ' accumulated precip rate=', fld(140*73)
!         endif
	

          lupgb=lupgb+2
          if( lupgi .ne. 0 ) lupgi=lupgi+2
        enddo                             ! end second do loop for n
!caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
!
!cbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb

!     start to write out average value

        ndata = ijdim
        data  = 0.
        lbmo  = .false.
        do i=1,ijdim  
          if (ivr(i) > 0) then
            fac     = 1. / float(ivr(i))
            data(i) = fld(i)*fac
            lbmo(i) = .true.
          endif
        enddo

!     if (nrank == dbgme) then
!       print *,'After averaging data=',data(36481),' fld=',fld(36481),'ivr=',ivr(36481)
!       do j=1,190
!         do i=1,384
!           write (333,*)' data=',data((j-1)*384+i),' i=',i,'j=',j, lbmo((j-1)*384+i)
!         enddo
!       enddo
!     endif
!       if(kpds(5).eq.59) then
!         print *, ' the averaged precip rate =', data(ijdim/2)
!       endif

        kgds     = jgds
        kpds     = jpds
        kpds(8)  = iy
        kpds(9)  = im
        kpds(10) = id
        kpds(11) = ih
!
        if (fcst_avrg) then
          print *,' iy=',iy,' iys=',iys , ' iyb=',iyb
          KPDS(13) = 3                ! Forecast time unit - month
          KPDS(14) = (IYS-IYB)*12 + IMS - IM ! Start month
          KPDS(15) =  KPDS(14) + 1           ! End   month
          KPDS(16) = 3                ! Time range inidicator - average
          KPDS(17) = nt               ! Number of time included in avrg
          KPDS(20) = 0                ! Number of missing time
          KPDS(21) = ((IYB-1)/100) + 1
        else

!         kpds(16) = kpds16
!         IF(KPDS(16).GE.2 .AND. KPDS(16).LE.5) THEN
!           kpds(15) = ift
!         ELSE
!           kpds(14) = ift
!         ENDIF

!     print *,' nt=',nt

          kpds(17) = nt
!
          if(kpds(16) == 10) then
            if(kpds(14) == 0) then        ! Average of analysis
              kpds(16) = 123              ! Octet no. 21
            else                          ! Average of forecasts
              kpds(16) = 113              ! Octet no. 21
            endif
            item = mod(nhours,24)
            if (item == 0 ) then       !  for time mean for a given cycle
              kpds(11) = ih
              kpds(15) = 24
            elseif (nhours < 24) then  !  for time mean of daily mean
              kpds(11) = 0
              kpds(15) = nhours
            else
              print *,' Bad choice of nhours =', nhours
            endif
          elseif (kpds(16) == 3 ) then    ! Average
            if (nhours == 24) then
              kpds(16) = 130              ! Octet no. 21
            elseif (nhours == 6) then
              kpds(16) = 138              ! Octet no. 21
            elseif (nhours == 12) then
              kpds(16) = 140              ! Octet no. 21
            else
              print *,' No valid kpds16 for this nhours=',nhours        &
     &,                ' using kpds(16)=113'
              kpds(16) = 113
            endif
          elseif (kpds(16) == 4 ) then    ! Accumulation
            if (nhours == 24) then
              kpds(16) = 128              ! Octet no. 21
            elseif (nhours == 6) then
              kpds(16) = 137              ! Octet no. 21
            elseif (nhours == 12) then
              kpds(16) = 139              ! Octet no. 21
            else
              print *,' No valid kpds16 for this nhours=',nhours        &
     &,                ' using kpds(16)=113'
              kpds(16) = 113
            endif
          else
            print *,' The input data does not have supported kpds16=',kpds16
          endif
!
        endif
!
!
!... add precision to monthly average...
        kpds(22) = kpds(22) + 1
!
!... add precision to the snow....
!          if((kpds(5) == 65).and.(kpds(6) == 1).and.(kpds(7) == 0))    &
!    &         kpds(22) = kpds(22) + 1

!     if (nrank == dbgme) then
!       write(334,*) data
!       write(334) data
!       do j=1,190
!         do i=1,384
!           write (335,*)' data=',data((j-1)*384+i),' i=',i,'j=',j,lbmo((j-1)*384+i)
!         enddo
!       enddo
!      endif

!       print *,' calling putgb for lupgbo=',lupgbo,' nrank=',nrank

        call putgb(lupgbo,ndata,kpds,kgds,lbmo,data,iret)

        if (iret .ne. 0) then
          print*,'putgb failed for lupgbo,k,jpds(5),jpds(6),jpds(7),iret ',lupgbo,k,jpds(5),jpds(6),jpds(7),iret
          call mpi_abort(MPI_COMM_WORLD,333,ierr)
        endif

!cbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb
      enddo                             ! end first do loop for kk
      deallocate ( field,stat=iret)
      deallocate ( fld,stat=iret)
      deallocate ( data,stat=iret)
      deallocate ( ivr,stat=iret)
      deallocate ( lbmi,stat=iret)
      deallocate ( lbmo,stat=iret)
      lupgb = 11
      if(lupgi.ne.0) lupgi = 12
      do n = 1, ntimes                   ! second do loop for n file
        call baclose(lupgb,iret)
        if (iret.ne.0) then
          write(6,*) 'baclose failure ',cfiles(1,1),' for 1,lupgb,iret =', 1,lupgb,iret
          call mpi_abort(MPI_COMM_WORLD,333,ierr)
        endif
        if ( lupgi .ne. 0 )  then
          call baclose(lupgi,iret)
          if (iret.ne.0) then
            write(6,*) 'baclose failure ',cfiles(2,1),' for 2,lupgb,iret =', 2,lupgb,iret
            call mpi_abort(MPI_COMM_WORLD,333,ierr)
          endif
          lupgi=lupgi+2
        endif
        lupgb = lupgb + 2
      enddo
      call baclose(lupgbo,iret)
      if (iret.ne.0) then
        write(6,*) 'baclose failure lupgbo,iret = ', lupgbo,iret
        call mpi_abort(MPI_COMM_WORLD,333,ierr)
      endif
      call mpi_finalize(iret)
!
      if (nrank == 0) call summary()
!

      stop
      end