subroutine shap_lev(zz,m,n,key,nord) c c======================================================================= c === c PE model version of the Shapiro Filter. === c === c Arguments: === c === c M the size of first dimension of the array to filter === c N the size of second dimension of the array to filter === c NORD the order of the filter to apply === c ZZ the array to be filtered === c KEY grid type [t-grid or u-grid] === c === c======================================================================= c #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,idown,iup,icolb,icole,iend,ii,iodev,ip,istr,j,jend,jj, * jrowb,jrowe,jstr,kord,m,n,nord,key #ifdef coast integer l # if defined pindex & defined shapnocoastflux integer iz,jz # endif #endif FLOAT * fac,g(0:ijmx+1,0:1),zz(nwds) c c----------------------------------------------------------------------- c Begin excutable code. c----------------------------------------------------------------------- c iodev=(nord+1)/2-nord/2 fac=cm1+c2*FLoaT(iodev) fac=fac/c2**(2*nord) c c Set default bdy condition c #ifdef cyclic jrowb=1 icolb=1 #else if(key.eq.zgrid) then jrowb=2 else jrowb=1 endif if(key.eq.zgrid) then icolb=2 else icolb=1 endif #endif jrowe=n icole=m iup=1 c c Filter by rows. c do 60 j=jrowb,jrowe ip=(j-1)*imt #ifndef coast if(key.eq.zgrid) then istr=icolb iend=icole else istr=icolb+1 iend=icole-1 endif #else # if defined pindex & defined shapnocoastflux jz=j if (j.eq.1) jz=2 if (j.eq.n) jz=n-1 # endif do 55 l=1,lseg if(key.eq.vgrid) then istr=isq(j,l) elseif(key.eq.tgrid) then #if defined pindex & defined shapnocoastflux istr=isp(jz,l) #else istr=isz(j,l) #endif else if(isz(j,l).eq.0)goto 60 istr=max(isz(j,l),2) endif if(istr.eq.0) goto 60 if(key.eq.vgrid) then iend=ieq(j,l) elseif(key.eq.tgrid) then #if defined pindex & defined shapnocoastflux iend=iep(jz,l) #else iend=iez(j,l) #endif else iend=min(iez(j,l),imtm1) endif #endif if(istr.ge.iend) goto 55 ii=0 if (key.eq.zgrid) then g(ii,iup)=zz(ip+istr+1) #if defined pindex & defined shapnocoastflux elseif(key.eq.tgrid)then g(ii,iup)=zz(ip+istr) #endif else g(ii,iup)=zz(ip+istr-1) endif do 10 i=istr,iend ii=ii+1 g(ii,iup)=zz(ip+i) 10 continue if (key.eq.zgrid) then g(ii+1,iup)=zz(ip+iend-1) #if defined pindex & defined shapnocoastflux elseif(key.eq.tgrid)then g(ii+1,iup)=zz(ip+iend) #endif else g(ii+1,iup)=zz(ip+iend+1) endif do 40 kord=1,nord idown=1-iup do 20 i=1,ii g(i,idown)=cm2*g(i,iup)+g(i+1,iup)+g(i-1,iup) 20 continue iup=1-iup #ifdef cyclic # ifndef coast g(0,iup)=g(2,iup) g(ii+1,iup)=g(ii-1,iup) # else if(istr.eq.2) then g(0,iup)=g(ii,iup) g(ii+1,iup)=g(1,iup) else g(0,iup)=g(2,iup) g(ii+1,iup)=g(ii-1,iup) endif # endif #else g(0,iup)=g(2,iup) g(ii+1,iup)=g(ii-1,iup) #endif 40 continue ii=0 do 50 i=istr,iend ii=ii+1 zz(ip+i)=zz(ip+i)+fac*g(ii,iup) 50 continue #ifdef cyclic # ifndef coast zz(ip+1)=zz(ip+icole-1) zz(ip+icole)=zz(ip+2) # else if(istr.eq.2) then zz(ip+1)=zz(ip+icole-1) zz(ip+icole)=zz(ip+2) endif # endif #endif 55 continue 60 continue c c Filter by columns. c do 120 i=icolb,icole #ifndef coast if(key.eq.zgrid) then jstr=jrowb jend=jrowe else jstr=jrowb+1 jend=jrowe-1 endif #else # if defined pindex & defined shapnocoastflux iz=i if (i.eq.1) iz=2 # endif do 115 l=1,lseg if(key.eq.vgrid) then jstr=jsq(i,l) elseif(key.eq.tgrid) then #if defined pindex & defined shapnocoastflux jstr=jsp(iz,l) #else jstr=jsz(i,l) #endif else if(jsz(i,l).eq.0) goto 120 jstr=max(jsz(i,l),2) endif if(jstr.eq.0) goto 120 if(key.eq.vgrid) then jend=jeq(i,l) elseif(key.eq.tgrid) then #if defined pindex & defined shapnocoastflux jend=jep(iz,l) #else jend=min(jez(i,l),jmtm2) #endif else jend=min(jez(i,l),jmtm1) endif #endif if(jstr.ge.jend) goto 115 jj=0 if(key.eq.zgrid) then g(jj,iup)=zz((jstr)*imt+i) #if defined pindex & defined shapnocoastflux elseif(key.eq.tgrid) then g(jj,iup)=zz((jstr-1)*imt+i) #endif else g(jj,iup)=zz((jstr-1-1)*imt+i) endif do 70 j=jstr,jend jj=jj+1 ip=(j-1)*imt+i g(jj,iup)=zz(ip) 70 continue if(key.eq.zgrid) then g(jj+1,iup)=zz((jend-2)*imt+i) #if defined pindex & defined shapnocoastflux elseif(key.eq.tgrid) then g(jj+1,iup)=zz((jend-1)*imt+i) #endif else g(jj+1,iup)=zz((jend+1-1)*imt+i) endif do 100 kord=1,nord idown=1-iup do 80 j=1,jj g(j,idown)=cm2*g(j,iup)+g(j+1,iup)+g(j-1,iup) 80 continue iup=1-iup g(0,iup)=g(2,iup) g(jj+1,iup)=g(jj-1,iup) 100 continue jj=0 do 110 j=jstr,jend jj=jj+1 ip=(j-1)*imt+i zz(ip)=zz(ip)+fac*g(jj,iup) 110 continue 115 continue 120 continue return end