       parameter (ngx=501, ngy=191)
       character*12 cdate
       dimension cd(ngx,ngy),sig_wht(ngx,ngy),
     + pper(ngx,ngy),z0(ngx,ngy),beta(ngx,ngy)
       dimension sig_wht0(ngx,ngy),rms(ngx,ngy),
     +  diff(ngx,ngy),pper0(ngx,ngy),beta_ini(ngx,ngy),
     +  bias(ngx,ngy),beta0(ngx,ngy)
       dimension swh_m0(ngx,ngy),swh_m1(ngx,ngy)


        rms = 0.
        bias = 0.
c SPECIAL VALUE IS 9999.

c cdate = yyyymmddhhmm

       open(60,file='WAVES.INPUT_1',form='unformatted')
       open(61,file='WAVES.INPUT_0',form='unformatted')
       open(80,file='beta_diff18')


       iinc = 0
       swh_m0 = 0.
       swh_m1 = 0.

c      do itimes = 1,165
       do itimes = 1,48 
       read(60) cdate
       print *,cdate
       read(60) ((cd(i,j),i=1,ngx),j=1,ngy)
       read(60) ((sig_wht(i,j),i=1,ngx),j=1,ngy)
       read(60) ((pper(i,j),i=1,ngx),j=1,ngy)
       read(60) ((z0(i,j),i=1,ngx),j=1,ngy)
       read(60) ((beta(i,j),i=1,ngx),j=1,ngy)

       read(61) cdate
       read(61) ((cd(i,j),i=1,ngx),j=1,ngy)
       read(61) ((sig_wht0(i,j),i=1,ngx),j=1,ngy)
       read(61) ((pper0(i,j),i=1,ngx),j=1,ngy)
       read(61) ((z0(i,j),i=1,ngx),j=1,ngy)
       read(61) ((beta0(i,j),i=1,ngx),j=1,ngy)

        if (itimes.ge.22 .and. itimes.le.45) then
        iinc = iinc + 1
         do i = 1,ngx
         do j = 1,ngy
          if (sig_wht0(i,j).ne.9999.) then
         swh_m0(i,j) = swh_m0(i,j) + sig_wht0(i,j)
         swh_m1(i,j) = swh_m1(i,j) + sig_wht(i,j)
          end if 
         end do
         end do
        end if



c       if (itimes.eq.1) then
c        beta_ini = beta
c       else
c       beta = beta - beta_ini
        do i = 1,ngx
        do j = 1,ngy
        if (beta(i,j).ne.9999.) then
c       beta(i,j) = beta(i,j) - 0.018
        end if
        end do
        end do
          write(80,300) ((beta(i,j),i=1,ngx),j=1,ngy)
c       end if

           do i = 1,ngx
           do j = 1,ngy
            diff(i,j) = sig_wht(i,j)-sig_wht0(i,j)
           end do
           end do
            amax = -999.
           do i = 1,ngx
           do j = 1,ngy
           val=abs(diff(i,j))
            if (val.gt.amax) amax=val
           end do
           end do
            print *,cdate,' ','max val:',amax
 

c         if (cdate.eq.'110322000000') then
c          print *,'done'
c          do i = 1,ngx
c          do j = 1,ngy
c           diff(i,j) = 9999.
c           if (sig_wht(i,j).ne.9999.) then
c           diff(i,j) = sig_wht(i,j)-sig_wht0(i,j)
c           end if
c          end do
c          end do
c          open(30,file='wave_res')
c         write(30,300) ((diff(i,j),i=1,ngx),j=1,ngy)
c         stop
c         end if
          
           do i = 1,ngx
           do j = 1,ngy
           rms(i,j)=rms(i,j)+(sig_wht(i,j)-sig_wht0(i,j))**2.
           bias(i,j)=bias(i,j)+(sig_wht(i,j)-sig_wht0(i,j))
c          rms(i,j)=rms(i,j)+(beta(i,j)-beta0(i,j))**2.
c          bias(i,j)=bias(i,j)+beta0(i,j)
           end do
           end do
          

       end do

          do i = 1,ngx
          do j = 1,ngy
          if (sig_wht(i,j).eq.9999.) then
          rms(i,j)=9999.
          bias(i,j)=9999.
          else
c         rms(i,j)=sqrt(rms(i,j)/165.)
c         bias(i,j)=bias(i,j)/165. 
          rms(i,j)=sqrt(rms(i,j)/48.)
          bias(i,j)=bias(i,j)/48. 
          end if
          end do
          end do

          close(60)
          close(61)

           open(30,file='wave_res')
          write(30,300) ((bias(i,j),i=1,ngx),j=1,ngy)
          write(30,300) ((rms(i,j),i=1,ngx),j=1,ngy)
          close(30)
 300  FORMAT(8e10.3)

          print *, iinc
         do i = 1,ngx
         do j = 1,ngy
        if (sig_wht0(i,j).ne.9999.) then
         swh_m0(i,j) = swh_m0(i,j)/float(iinc)
         swh_m1(i,j) = swh_m1(i,j)/float(iinc)
        else
         swh_m0(i,j)=9999.
         swh_m1(i,j)=9999.
        end if
         end do
         end do

           open(30,file='SWH_MOD')
          write(30,300) ((swh_m0(i,j),i=1,ngx),j=1,ngy)
          write(30,300) ((swh_m1(i,j),i=1,ngx),j=1,ngy)
          close(30)




       stop
       end

