      parameter (nl=128,nt=96)
c
      real yrPEE1(NL+1,NT/2+1), yrPEE2(NL+1,NT/2+1)
      real sumPEE(0:NL+2,1:NT/2+2)
      real smth1(2:nt/2+1)
      real smth2(1:nl+1)

	open (1,file='./gdat/out.pres',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
	open (2,file='./gdat/out.prea',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
c
	open (11,file='./gdat/out_log.pres',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
	open (12,file='./gdat/out_log.prea',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
c
      read (1,rec=1) yrPEE1
	do i=1,nl+1
	 yrPEE1(i,1)=-999.
	do j=2,nt/2+1
	  yrPEE1(i,j)=alog10(yrPEE1(i,j))
	enddo
	enddo
	
      read (2,rec=1) yrPEE2
	do i=1,nl+1
	 yrPEE2(i,1)=-999.
	do j=2,nt/2+1
	  yrPEE2(i,j)=alog10(yrPEE2(i,j))
	enddo
	enddo	
c
      write (11,rec=1) yrPEE1
      write (12,rec=1) yrPEE2
c
	open (13,file='./gdat/out_log.prcp.sum',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
c
	do i=1,nl+1
	  sumPEE(i,1)=-999.
	do j=2,nt/2+1
	  if(yrPEE1(i,j).ne.-999.) then
	    sumPEE(i,j)=0.5*yrPEE1(i,j)+0.5*yrPEE2(i,j)
	  else
	    sumPEE(i,j)=-999.
	  endif
	enddo
	enddo	
      write (13,rec=1) ((sumPEE(i,j),i=1,nl+1),j=1,nt/2+1)
c
c! set boundary values
	do i=0,nl+2
	 sumPEE(i,1) = sumPEE(i,2)
	 sumPEE(i,nt/2+2) = sumPEE (i,nt/2+1)
	enddo
c
	do j=1,nt/2+2
	 sumPEE(0,j) = sumPEE(1,j)
	 sumPEE(nl+2,j) = sumPEE (nl+1,j)
	enddo
c
	open (14,file='./gdat/out_log.prcp.sum_121',access='direct',
     *   form='unformatted',recl=(NL+1)*(NT/2+1),status='unknown')
c
c! 1-2-1 for frequency
	itr=0
1111    continue
	itr=itr+1
	print*,'itr=',itr

	do 1110 i=1,nl+1

	do j=2,nt/2+1
	  smth1(j)=0.
	enddo

	do j=2,nt/2+1
	 smth1(j)=0.25*sumPEE(i,j-1)+0.5*sumPEE(i,j)+0.25*sumPEE(i,j+1)
	enddo
c
	do j=2,nt/2+1
	 sumPEE(i,j)=smth1(j)
	enddo

1110	continue
c
	if (itr.eq.10) goto 1112
	goto 1111
c
	 
1112    continue
c
c! 1-2-1 for wavenumber
c
	itr=0
1114    continue
	itr=itr+1
	print*,'itr=',itr

	do 1113 j=1,nt/2+1

        do i=1,nl+1
          smth2(i)=0.
	enddo

	do i=1,nl+1
	 smth2(i)=0.25*sumPEE(i-1,j)+0.5*sumPEE(i,j)+0.25*sumPEE(i+1,j)
	enddo
c
	do i=1,nl+1
	 sumPEE(i,j)=smth2(i)
	enddo

1113	continue
c
	if (itr.eq.40) goto 1115
	goto 1114
c
	 
1115    continue
c
998     continue
c
! flagged for undefined
	do i=1,nl+1
	do j=1,nt/2+1
	  if(yrPEE1(i,j).eq.-999.) sumPEE(i,j)=-999.
	enddo
	enddo
c
        write (14,rec=1) ((sumPEE(i,j),i=1,nl+1),j=1,nt/2+1)
999	stop
	end	
