program PWREFIND
C     Calculating averages 

	integer RD, RDS,rns,runs,ILAT,kgit
     	PARAMETER (runs=10)
      PARAMETER (ILAT=36)
      parameter (kgit=56)
	parameter (nvar=2)

      real, dimension(ILAT,kgit,nvar)::aus,rea,var,num
      real, dimension(ILAT,kgit,nvar)::auh,aul,rea2,Dis
      real, dimension(runs,ILAT,kgit,nvar)::aus1,aus2
	integer spantime,starttime

	integer h,p
	CHARACTER*1 jfilein(runs),jfilein2(runs)

	spantime = 373	   ! num steps
	starttime = 361

c=========================input file & computing area====================================

      DATA jfilein/'0','1','2','3','4','5'
     &,'6','7','8','9'/
      DATA jfilein2/'0','1','2','3','4','5'
     &,'6','7','8','9'/ 




c================== DATA OUTPUT =================================

      open (71,status='unknown',form='unformatted',file=
     +'resid_wqbo.dx', access='direct',
     + recl=4*ILAT*kgit*nvar)
      open (72,status='unknown',form='unformatted',file=
     +'resid_eqbo.dx', access='direct',
     + recl=4*ILAT*kgit*nvar)
      open (73,status='unknown',form='unformatted',file=
     +'resid_disp.dx', access='direct',
     + recl=4*ILAT*kgit*nvar)

	do rns = 1,runs
	  RD = rns+20
	  RDS = rns+120

c================== DATA INPUT =================================

      OPEN (RD,FILE= "../../Muam_Mar_wQBO/"//jfilein(rns)//"/resid.dx",
     *form='unformatted',
     *access='direct', status='old',recl=4*ILAT*kgit*nvar) 
      OPEN(RDS,FILE= "../../Muam_Mar_eQBO/"//jfilein2(rns)//"/resid.dx",
     *form='unformatted',
     *access='direct', status='old',recl=4*ILAT*kgit*nvar) 
	  print*,"open: ../../" //jfilein(rns)// "/resid.dx"

	enddo




   	DO 84 l=1,spantime	!timeloop for intervals

	print*, 'step : ', l, 'curr : ', starttime
	

	do rns = 1,runs	

	  RD = rns+20
	  RDS = rns+120

	read(RD,rec=starttime) rea
	read(RDS,rec=starttime) rea2


			do p=1,nvar
                do j=1,ILAT
                   do h=1,kgit

			aus1(rns,j,h,p)=rea(j,h,p)
	        aus2(rns,j,h,p)=rea2(j,h,p)

				end do
			  end do
			end do


	enddo

c================================================================

	do p=1,nvar
        do j=1,ILAT
          do h=1,kgit
	         do rns = 1,runs

c-------aver---------	
	auh(j,h,p)=auh(j,h,p)+aus1(rns,j,h,p)
      aul(j,h,p)=aul(j,h,p)+aus2(rns,j,h,p) 

c-------disp---------	
	aus(j,h,p)=aus(j,h,p)+(aus1(rns,j,h,p)-aus2(rns,j,h,p))
      var(j,h,p)=var(j,h,p)+(aus1(rns,j,h,p)-aus2(rns,j,h,p))**2

	num(j,h,p) = num(j,h,p) + 1

	         end do

          end do
        end do
	end do

!=====================zanulyaem
	do p=1,nvar
        do j=1,ILAT
          do h=1,kgit
	         do rns = 1,runs 
			             aus1(rns,j,h,p)=0
			             aus2(rns,j,h,p)=0
	         enddo
			rea(j,h,p)=0
			rea2(j,h,p)=0				
          end do
        end do
	end do
!===============================


	starttime = starttime + 1

  84	ENDDO  ! time loop

 	print*,'num vars = ', num(1,1,1)

	do p=1,nvar
        do j=1,ILAT
          do h=1,kgit

        if(num(j,h,p).gt.0) then

	aus(j,h,p)=aus(j,h,p)/num(j,h,p)
	var(j,h,p)=var(j,h,p)/num(j,h,p)

	auh(j,h,p)=auh(j,h,p)/num(j,h,p)
	aul(j,h,p)=aul(j,h,p)/num(j,h,p)

        endif 

	Dis(j,h,p)=sqrt(abs(var(j,h,p)-aus(j,h,p)*aus(j,h,p)))

          end do
        end do
	end do


	write(71,rec=1) auh
	write(72,rec=1) aul
	write(73,rec=1) Dis



c============reset all, close all etc ==============================
 
 	do rns = 1,runs
  	RD = rns+20
        close(RD)
  	RDS = rns+120
        close(RDS)
	end do
        close(71)
        close(72)
        close(73)




      STOP
      END
*************************************************************************