subroutine lap_filt (f,keytv) c c======================================================================= c === c This subroutine filters 2-D (barotropic) fields. === c === c Arguments: === c === c F The field to be filtered (real array; input/output) === c KEYTV Flag to indicate on which grid F is defined === c (integer; input): === c KEYTV = 0 => tracer grid. === c KEYTV = 1 => velocity grid. === c === c Calls: LAPT_LEV. === c === c Common Blocks: (only relevant variables documented) === c === #ifdef coast c /FULLWD/ === c === c LANDT Land mask at tracer points (integer array; input). === c LANDV Land mask at velocity points (integer array; input). === c === #endif c /ONEDIM/ === c === c CS Cosine metric factors at velocity points (real array;=== c input). === c CSR Reciprocal of CS (real array; input). === c CST Cosine metric factors at tracer points (real array; === c input). === c CSTR Reciprocal of CST (real array; input) === c DXTR Reciprocal of X-width of tracer boxes (real array; === c input). === c DXUR Reciprocal of X-width of velocity boxes (real array; === c input). === c DYTR Reciprocal of X-width of tracer boxes (real array; === c input). === c DYUR Reciprocal of X-width of velocity boxes (real array; === c input). === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #ifdef coast # include #endif c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,idn,j,jdn,keytv,imt3,jptr,jptrm,jptrp FLOAT * dxin,dxout,dxx,dyin,dyout,dyy FLOAT * f(imt,jmt),lap_ug(imt,jmt),lap_xx(imt,jmt),lap_yy(imt,jmt), * msk(imt,0:2) parameter(imt3=3*imt) save msk,jptrm,jptr,jptrp data msk,jptrm,jptr,jptrp/imt3*c1,0,1,2/ c c======================================================================= c Begin executable code. c======================================================================= c c----------------------------------------------------------------------- c Get Laplacian of supplied field. c----------------------------------------------------------------------- c #ifdef coast if (keytv.eq.tgrid) then do 10 i=1,imt if (landt(i,1).gt.0) then msk(i,jptr)=c1 msk(i,jptrm)=c1 else msk(i,jptr)=c0 msk(i,jptrm)=c0 end if if (landt(i,2).gt.0) then msk(i,jptrp)=c1 else msk(i,jptrp)=c0 end if 10 continue else do 20 i=1,imt if (landv(i,1).gt.0) then msk(i,jptr)=c1 msk(i,jptrm)=c1 else msk(i,jptr)=c0 msk(i,jptrm)=c0 end if if (landv(i,2).gt.0) then msk(i,jptrp)=c1 else msk(i,jptrp)=c0 end if 20 continue end if #endif if (keytv.eq.tgrid) then call lapt_lev(1,1,f(1,1),f(1,1),f(1,2),msk(1,jptrm),msk(1,jptr), * msk(1,jptrp),lap_xx(1,1),lap_yy(1,1)) else call lapv_lev(1,1,f(1,1),f(1,1),f(1,2),lap_xx(1,1),lap_yy(1,1), * lap_ug(1,1)) endif do 50 j=2,jmtm1 #ifdef coast jptr=mod(jptr+1,3) jptrm=mod(jptrm+1,3) jptrp=mod(jptrp+1,3) if (keytv.eq.tgrid) then do 30 i=1,imt if (landt(i,j+1).gt.0) then msk(i,jptrp)=c1 else msk(i,jptrp)=c0 end if 30 continue else do 40 i=1,imt if (landv(i,j+1).gt.0) then msk(i,jptrp)=c1 else msk(i,jptrp)=c0 end if 40 continue end if #endif if (keytv.eq.tgrid) then call lapt_lev(j,1,f(1,j-1),f(1,j),f(1,j+1),msk(1,jptrm), * msk(1,jptr),msk(1,jptrp),lap_xx(1,j),lap_yy(1,j)) else call lapv_lev(j,1,f(1,j-1),f(1,j),f(1,j+1),lap_xx(1,j), * lap_yy(1,j),lap_ug(1,j)) endif 50 continue if (keytv.eq.tgrid) then call lapt_lev(jmt,1,f(1,jmtm1),f(1,jmt),f(1,jmt),msk(1,jptr), * msk(1,jptrp),msk(1,jptrp),lap_xx(1,jmt),lap_yy(1,jmt)) else call lapv_lev(jmt,1,f(1,jmtm1),f(1,jmt),f(1,jmt), * lap_xx(1,jmt),lap_yy(1,jmt),lap_ug(1,jmt)) endif c c----------------------------------------------------------------------- c Recover grid spacing and use laplacian to filter supplied field. c----------------------------------------------------------------------- c if(keytv.eq.tgrid) then do 60 j=1,jmt jdn=max(j-1,1) dyin=p5*(csr(jdn)*dyu(jdn)+csr(j)*dyu(j)) dyout=cst(j)*dyt(j) dyy=dyin*dyout do 60 i=1,imt idn=max(i-1,1) dxin=cst(j)*p5*(dxu(idn)+dxu(i)) dxout=cst(j)*dxt(i) dxx=dxout*dxin f(i,j)=f(i,j)+dxx*lap_xx(i,j)+dyy*lap_yy(i,j) 60 continue else do 70 j=1,jmt jdn=min(j+1,jmt) dyin=p5*(cstr(jdn)*dyt(jdn)+cstr(j)*dyt(j)) dyout=cs(j)*dyu(j) dyy=dyin*dyout do 70 i=1,imt idn=min(i+1,imt) dxin=cs(j)*p5*(dxt(idn)+dxt(i)) dxout=cs(j)*dxu(i) dxx=dxout*dxin f(i,j)=f(i,j)+dxx*lap_xx(i,j)+dyy*lap_yy(i,j) 70 continue endif return end