c Fortran program example to compute mean and std. deviation
c for precipitation data in homework #3

c Binary file version


c Define array limits
c For the HW #3 data these correspond to the data dimensions
c in space and time

      parameter (xpoints=81)
      parameter (ypoints=51)
      parameter (months=12)
      parameter (years=55)


c Counting dimesions for space and time

      integer i,j,p,q

c Other variable declarations

      real missingval

c Precipitation declared as four dimensional array
c All other arrays necessary for computation 3-D

      
      real precip(xpoints,ypoints,months,years)
      real missingarray(xpoints,ypoints,months)
      real sumprecip(xpoints,ypoints,months)
      real sumprecipsq(xpoints,ypoints,months)
      real meanprecip(xpoints,ypoints,months)
      real stdevprecip(xpoints,ypoints,months)
 

c Input and output files
c  In Fortran, files need reference numbers
      

      open (10,file='precipmonthly_usmx_le.bin',
     + form='unformatted', access='direct', 
     + recl=4*xpoints*ypoints*months*years)

      open (11,file='Mon_means_f.dat', 
     + form='unformatted', access='direct',
     + recl=4*xpoints*ypoints*months)

      open (12,file='Mon_std_f.dat',
     + form='unformatted', access='direct',
     + recl=4*xpoints*ypoints*months)

c Declare missing value

      missingval=-999.9

c Initialize missing array

      do p=1,months
        do j=1,ypoints
          do i=1,xpoints
            missingarray(i,j,p)=0
          end do
        end do
      end do

c Read in data

      read(10,rec=1) ((((precip(i,j,p,q),i=1,xpoints),
     + j=1,ypoints),p=1,months),q=1,years)
      
c Keep track of missing values
 
      do q=1,years
        do p=1,months
          do j=1,ypoints
            do i=1,xpoints
              if(precip(i,j,p,q).eq.missingval) then
                missingarray(i,j,p)=1
              end if
            end do
          end do
        end do
      end do


c Initialize summation arrays

      do p=1,months
        do j=1,ypoints
          do i=1,xpoints
            sumprecip(i,j,p)=0
            sumprecipsq(i,j,p)=0
          end do
        end do
      end do

c Compute mean and std. deviation, ignoring missing values

      do p=1,months
        do q=1,years
          do j=1,ypoints
            do i=1,xpoints
              if(missingarray(i,j,p).eq.0) then
                sumprecip(i,j,p)=sumprecip(i,j,p)+
     +             precip(i,j,p,q)
                sumprecipsq(i,j,p)=sumprecipsq(i,j,p)+
     +             precip(i,j,p,q)*precip(i,j,p,q)
              end if
            end do
          end do
        end do
      end do


      do p=1,months
        do j=1,ypoints
          do i=1,xpoints
            if(missingarray(i,j,p).eq.0) then
              meanprecip(i,j,p)=sumprecip(i,j,p)/years
              stdevprecip(i,j,p)=sqrt((years*sumprecipsq(i,j,p)-
     +         sumprecip(i,j,p)*sumprecip(i,j,p))/(years*
     +         (years-1)))
            end if
            if(missingarray(i,j,p).eq.1) then
              meanprecip(i,j,p)=missingval
              stdevprecip(i,j,p)=missingval
            end if
          end do
        end do
      end do 


c Write out data to binary file


      write(11,rec=1) (((meanprecip(i,j,p),i=1,xpoints)
     + ,j=1,ypoints),p=1,months)
      write(12,rec=1) (((stdevprecip(i,j,k),i=1,xpoints)
     + ,j=1,ypoints),p=1,months)

      close(10)
      close(11)
      close(12)      


      end