subroutine box_ave(fld_in,area_in,undef,gxout,gyout, $ nii,nji,nio,njo,fld_out,iunit_diag,vote,istat, $ rmin_vote_max,rmin_vote_min) parameter (na=300) dimension fld_in(nii,nji),area_in(nii,nji),fld_out(nio,njo), $ gxout(nio+1),gyout(njo+1) dimension area_box(na),fld_box(na),dxdy_box(na),ifld_rank(na), $ fld_cand(na),area_cand(na),dxdy_cand(na),rmin_dxdy_vote(3) logical cyclicx,diag_out,vote diag_out=.false. cccc diag_out=.true. istat=1 C C tolerance for checking whether grid is undefined (zero sfc area) C eps=1.0e-12 C C minimum fractional area of the output grid box C in order to have a winner in the voting C rmin_dxdy_vote(1)=rmin_vote_max rmin_dxdy_vote(2)=(rmin_vote_max+rmin_vote_min)*0.5 rmin_dxdy_vote(3)=rmin_vote_min ibb=1 iee=nio jbb=1 jee=njo do j=jbb,jee do i=ibb,iee ib=int(gxout(i))+1 ie=int(gxout(i+1))+1 jb=int(gyout(j))+1 je=int(gyout(j+1))+1 C C check for exceeding n pole C if(je.gt.nji) je=nji C C cyclic continuity in x C cyclicx=.false. if(ie.lt.ib) then ie=ie+nii cyclicx=.true. endif if(diag_out) then write(iunit_diag,*) ' ' write(iunit_diag,*) $ 'i,j,ib,ie,jb,je = ',i,j,ib,ie,jb,je endif C C initialize the counter for intersecting grid boxes C icnt=0 C C CASE 1: only one input grid box in output grid box C if(ib.eq.ie.and.jb.eq.je) then icnt=1 dxdy_box(icnt)=1.0 area_box(icnt)=area_in(ib,jb) if(area_box(icnt).eq.0.0) dxdy_box(icnt)=0.0 fld_box(icnt)=fld_in(ib,jb) if(diag_out) then write(iunit_diag,*) $ 'ib=ie and jb=je ',area_box(icnt),dxdy_box(icnt), $ fld_box(icnt) endif else if(ib.eq.ie) then C C CASE 2: intersecting boxes in y only C ii=ib dx=gxout(i+1)-gxout(i) dxout=1.0 do jj=jb,je icnt=icnt+1 if(icnt.gt.na) then istat=0 go to 999 endif if(jj.eq.jb) then dy=float(jj)-gyout(j) else if(jj.eq.je) then dy=gyout(j+1)-float(jj-1) else dy=1.0 endif if(jj.eq.jb.or.jj.eq.je) then dyout=dy/(gyout(j+1)-gyout(j)) else dyout=1.0 endif dxdy_box(icnt)=dxout*dyout area_box(icnt)=dx*dy*area_in(ii,jj) if(area_in(ii,jj).eq.0.0) dxdy_box(icnt)=0.0 fld_box(icnt)=fld_in(ii,jj) if(diag_out) then write(iunit_diag,*) $ 'ib=ie ',ii,jj,dx,dy,area_in(ii,jj) write(iunit_diag,*) $ 'area_box ... ',area_box(icnt),fld_box(icnt) write(iunit_diag,*) $ 'dx dyout ... ',dxout,dyout,dxdy_box(icnt) endif end do else if(jb.eq.je) then C C CASE 3: intersecting boxes in x only C jj=jb dy=gyout(j+1)-gyout(j) dyout=1.0 do ii=ib,ie icnt=icnt+1 if(icnt.gt.na) then istat=0 go to 999 endif ii0=ii if(cyclicx.and.ii0.gt.nii) ii0=ii-nii if(ii.eq.ib) then dx=float(ii)-gxout(i) else if(ii.eq.ie) then x0=float(ii-1) if(cyclicx) x0=float(ii0)-1.0 dx=gxout(i+1)-x0 else dx=1.0 endif if(ii.eq.ib.or.ii.eq.ie) then dxout=dx/(gxout(i+1)-gxout(i)) else dxout=1.0 endif dxdy_box(icnt)=dxout*dyout area_box(icnt)=dx*dy*area_in(ii0,jj) if(area_in(ii0,jj).eq.0.0) dxdy_box(icnt)=0.0 fld_box(icnt)=fld_in(ii0,jj) if(diag_out) then write(iunit_diag,*) $ 'jb=je ',ii,ii0,jj,dx,dy,area_in(ii0,jj) write(iunit_diag,*) $ 'area_box ... ',area_box(icnt),fld_box(icnt) write(iunit_diag,*) $ 'dx dyout ... ',dxout,dyout,dxdy_box(icnt) endif end do else C C CASE 4: intersecting boxes in both directions C do jj=jb,je if(jj.eq.jb) then dy=float(jj)-gyout(j) else if(jj.eq.je) then dy=gyout(j+1)-float(jj-1) else dy=1.0 endif if(jj.eq.jb.or.jj.eq.je) then dyout=dy/(gyout(j+1)-gyout(j)) else dyout=1.0 endif do ii=ib,ie icnt=icnt+1 if(icnt.gt.na) then istat=0 go to 999 endif ii0=ii if(cyclicx.and.ii0.gt.nii) ii0=ii-nii if(ii.eq.ib) then dx=float(ii)-gxout(i) else if(ii.eq.ie) then x0=float(ii-1) if(cyclicx) x0=float(ii0)-1.0 dx=gxout(i+1)-x0 else dx=1.0 endif if(ii.eq.ib.or.ii.eq.ie) then dxout=dx/(gxout(i+1)-gxout(i)) else dxout=1.0 endif dxdy_box(icnt)=dxout*dyout area_box(icnt)=dx*dy*area_in(ii0,jj) if(area_in(ii0,jj).eq.0.0) dxdy_box(icnt)=0.0 fld_box(icnt)=fld_in(ii0,jj) if(diag_out) then write(iunit_diag,*) 'ib.ne.ib.and.jb.ne.je ', $ ii,ii0,jj,dx,dy,area_in(ii0,jj) write(iunit_diag,*) $ 'area_box ... ',area_box(icnt),fld_box(icnt) write(iunit_diag,*) $ 'dx dyout ... ',dxout,dyout,dxdy_box(icnt) endif end do enddo endif C C integrate or vote for the average value C if(vote) then C C voting routine; first get total area C tot_area=0.0 do ii=1,icnt tot_area=tot_area+area_box(ii) end do if(diag_out) then write(iunit_diag,*) 'i,j,icnt,tot_area = ', $ i,j,icnt,tot_area endif if(tot_area.le.eps) then fld_out(i,j)=undef go to 100 endif if(icnt.eq.1) then C C USSR election -- only one "candidate" to vote for C C check if the the total area is greater C than the minimum required to hold the election (e.g., 0.5) C if(diag_out) then write(iunit_diag,*) $ 'USSR: ',dxdy_box(1),rmin_dxdy_vote(1),fld_box(1) endif if(dxdy_box(1).lt.rmin_dxdy_vote(1)) then fld_out(i,j)=undef else fld_out(i,j)=fld_box(1) endif go to 100 else if(icnt.eq.2) then C C USA election -- two-party, two-candidate race; area wins C if(diag_out) then write(iunit_diag,*)'USA: ',dxdy_box(1),dxdy_box(2), $ rmin_vote endif if(dxdy_box(1).eq.0.0.or.dxdy_box(2).eq.0.0) then rmin_vote=rmin_dxdy_vote(1) else if(fld_box(1).eq.fld_box(2)) then rmin_vote=rmin_dxdy_vote(1) else rmin_vote=rmin_dxdy_vote(2) endif if((dxdy_box(1).ge.dxdy_box(2)).and. $ (dxdy_box(1).gt.rmin_vote)) then fld_out(i,j)=fld_box(1) else if(dxdy_box(2).gt.rmin_vote) then fld_out(i,j)=fld_box(2) C C case where both candidates are the same in the two-persion race C else if(fld_box(1).eq.fld_box(2).and. $ ((dxdy_box(1)+dxdy_box(2)).gt.rmin_vote)) then fld_out(i,j)=fld_box(2) else fld_out(i,j)=undef endif else C C a wide open election - three or more candidates C C sort the data by surface area using C the numercial recipes routine -- indexx C call indexx(icnt,fld_box,ifld_rank) if(diag_out) then write(iunit_diag,*) $ 'fld_box = ',(fld_box(ipp),ipp=1,icnt) write(iunit_diag,*) $ 'area_box = ',(area_box(ipp),ipp=1,icnt) write(iunit_diag,*) $ 'ifld_rank = ',(ifld_rank(ipp),ipp=1,icnt) write(iunit_diag,*) endif C C the indexes are in reverse order, with the C biggest fld element in the last element of the array C first check if the biggest is in the majority C C set up the candidates C ncand=1 it1=ifld_rank(1) area_cand(ncand)=area_box(it1) fld_cand(ncand)=fld_box(it1) dxdy_cand(ncand)=dxdy_box(it1) do ii=2,icnt i1=ii-1 i2=ii it1=ifld_rank(i1) it2=ifld_rank(i2) if(fld_box(it1).eq.fld_box(it2)) then area_cand(ncand)=area_cand(ncand)+area_box(it2) dxdy_cand(ncand)=dxdy_cand(ncand)+dxdy_box(it2) else ncand=ncand+1 area_cand(ncand)=area_box(it2) dxdy_cand(ncand)=dxdy_box(it2) fld_cand(ncand)=fld_box(it2) endif end do if(diag_out) then write(iunit_diag,*) 'ncand = ',ncand write(iunit_diag,*) $ 'fld_cand = ',(fld_cand(ipp),ipp=1,ncand) write(iunit_diag,*) $ 'area_cand = ',(area_cand(ipp),ipp=1,ncand) write(iunit_diag,*) $ 'dxdy_cand = ',(dxdy_cand(ipp),ipp=1,ncand) endif C C if one candidate, all done C if(ncand.eq.1) then if(dxdy_cand(1).gt.rmin_dxdy_vote(3)) then fld_out(i,j)=fld_cand(1) else fld_out(i,j)=undef endif go to 100 else C C the candidate with the most area is the winner C area_max=0.0 do ii=1,ncand if(area_cand(ii).gt.area_max) then iamax=ii area_max=area_cand(ii) endif end do if(ncand.le.2) then rmin_vote=rmin_dxdy_vote(ncand) else rmin_vote=rmin_dxdy_vote(3) endif if(dxdy_cand(iamax).gt.rmin_vote) then fld_out(i,j)=fld_cand(iamax) else fld_out(i,j)=undef endif if(diag_out) then write(iunit_diag,*) ' ' write(iunit_diag,*) 'the winner is...', $ iamax,area_max,dxdy_cand(iamax), $ rmin_vote,fld_out(i,j) endif go to 100 endif endif else C C area integrate C tot_fld=0.0 tot_area=0.0 C do ii=1,icnt tot_fld=tot_fld+fld_box(ii)*area_box(ii) tot_area=tot_area+area_box(ii) end do if(tot_area.gt.eps) then fld_out(i,j)=tot_fld/tot_area else fld_out(i,j)=undef endif endif 100 continue if(diag_out) then write(iunit_diag,*) $ 'i,j fld_out(i,j) = ',i,j,fld_out(i,j) endif end do end do 999 continue return end