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

c ASCII 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.asci',
     + form='formatted')

      open (11,file='Mon_means_f.asci', 
     + form='formatted')

      open (12,file='Mon_std_f.asci',
     + form='formatted')

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

      do q=1,years
        do p=1,months
          do j=1,ypoints
            do i=1,xpoints
              read(10,*) precip(i,j,p,q)
            end do
          end do
        end do
      end do

      
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 asci file

      do p=1,months
        do j=1,ypoints
          do i=1,xpoints
            write(11,*) meanprecip(i,j,p)
            write(12,*) stdevprecip(i,j,k)
          end do
        end do
      end do


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


      end