subroutine lap_depth (j,fjm1,fj,fjp1,mskm,msk,mskp,lap_xx,lap_yy, * lap_ug) c c======================================================================= c === c This routine computes the Laplacian of the given field in the === c specified slab and evaluated at the velocity points (entry === c LAPV_DEPTH) or at the tracer points (entry LAPT_DEPTH). The === c finite differencing is done along constant depths. === c === c On input: === c === c FJ The field in the slab J (real array). === c FJM1 The field in the slab max(J-1,1) (real array). === c FJP1 The field in the slab min(J+1,JMT) (real array). === c J The current slab number (integer). === c MSK Land mask in the slab J (real array). === c MSKM Land mask in the slab max(J-1,1) (real array). === c MSKP Land mask in the slab min(J+1,JMT) (real array). === c NO_LEV The number of vertical levels (integer). === c === c On output: === c === c LAP_UG The zonal gradient of the given velocity component. === c Used for velocity metric terms (real array) === c LAP_XX The zonal component of the Laplacian of the given === c field in the current slab (real array) === c LAP_YY The merdional component of the Laplacian of the === c given field in the current slab (real array) === c === #ifdef cubspl c Calls: INIT_SRCH, SPLINESLAB, SPLINTSLAB === #else c Calls: INIT_SRCH, LINTSLAB === #endif c === c Common Blocks: (only relevant variables documented) === c === 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 DXT2R Reciprocal of twice X-width of tracer boxes (real === c array; input) === c DXU X-width of velocity boxes (real array; input). === c DXUR Reciprocal DXU (real array; input). === c DYTR Reciprocal of Y-width of tracer boxes (real array; === c input). === c DYT2R Reciprocal of twice Y-width of tracer boxes (real === c array; input). === c DYU Y-width of velocity boxes (real array; input). === c DYUR Reciprocal DYU (real array; input). === c LPMTGD Metric factors for gradient metric term in velocity. === c (real array; input) === c === c /VERTSLAB/ === c === c DZQZ Vertical thickness of tracer boxes (real array; === c input) === c DZ2RQZ Reciprocal of 2 times DZQZ (real array; input) === c DZVQZ Vertical thickness of tracer boxes (real array; === c input) === c DZV2RQZ Reciprocal of 2 times DZVQZ (real array; input) === c JRN Pointer: J+1 slab (tracer) J slab (velocity) === c JRNP1 Pointer: J+2 slab (tracer) J+1 slab (velocity) === c JRS Pointer: J slab (tracer) J-1 slab (velocity) === c JRSM1 Pointer: J-1 slab (tracer) J-2 slab (velocity) === c TDEPTH Depths at center of tracer vertical boxes (real === c array; input) === c VDEPTH Depths at center of velocity vertical boxes (real === c array; input) === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c logical xtrpin,xtrpout integer i,j,k integer khi(imt,13),klow(imt,13) FLOAT * dyinn,dyins,dyout,fb,fbe,fbn,fbs,fbw,fe,fje,fjn,fjs,fjw,fn, * fs,fw,z,zbot FLOAT * dxin(imt),dxout(imt),fbotx(imt),fboty(imt),fj(imt,km), * fjm1(imt,km),fjp1(imt,km),ftopx(imt),ftopy(imt), * lap_ug(imt,km),lap_xx(imt,km),lap_yy(imt,km),msk(imt,km), * mskm(imt,km),mskp(imt,km) #ifdef cubspl FLOAT * spval FLOAT * dfz1(imt),dfzkm(imt),f2jm1(imt,km),f2j(imt,km), * f2jp1(imt,km),wk(imt,km) * parameter (spval=c1e30) data dfz1,dfzkm /imt*spval,imt*spval/ #endif c c======================================================================= c Compute Laplacian of field at velocty points. c======================================================================= c entry lapv_depth(j,fjm1,fj,fjp1,lap_xx,lap_yy,lap_ug) c c----------------------------------------------------------------------- c Grab grid spacings. c----------------------------------------------------------------------- c do 10 i=1,imum1 dxin(i)=dxtr(i+1)*csr(j) dxout(i)=dxur(i)*csr(j) 10 continue dxout(imu)=dxur(imu)*csr(j) dyout=dyur(j)*csr(j) dyins=dytr(j)*cst(j) dyinn=dytr(j+1)*cst(j+1) #ifdef cubspl c c----------------------------------------------------------------------- c Compute second derivatives for cubic splines. c----------------------------------------------------------------------- c call splineslab(vdepth(1,1,jrs),fjm1,imt,km,dfz1,dfzkm,wk,f2jm1) call splineslab(vdepth(1,1,jrn),fj,imt,km,dfz1,dfzkm,wk,f2j) call splineslab(vdepth(1,1,jrnp1),fjp1,imt,km,dfz1,dfzkm,wk,f2jp1) #endif c c----------------------------------------------------------------------- c Initialize search indices and bottom fluxes. c----------------------------------------------------------------------- c call init_srch(khi,klow,imt,km,13) c do 20 i=2,imum1 fbotx(i)=c0 fboty(i)=c0 20 continue c c----------------------------------------------------------------------- c Compute the laplacian in the slab. Note: extrapolations and c boundaries are set to give a no-flux boundary condition. c----------------------------------------------------------------------- c do 30 k=1,km do 30 i=2,imum1 #ifdef cubspl z=p5*(tdepth(i,k,jrn)+tdepth(i+1,k,jrn)) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i, * klow(i,1),khi(i,1),z,fjn,xtrpin) call splintslab(vdepth(1,1,jrnp1),fjp1,f2jp1,imt,km,i, * klow(i,2),khi(i,2),z,fn,xtrpout) if (xtrpin.or.xtrpout) fjn=fn z=p5*(tdepth(i+1,k,jrn)+tdepth(i+1,k,jrs)) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i, * klow(i,3),khi(i,3),z,fje,xtrpin) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i+1, * klow(i,4),khi(i,4),z,fe,xtrpout) if (xtrpin.or.xtrpout) fje=fe z=p5*(tdepth(i,k,jrs)+tdepth(i+1,k,jrs)) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i, * klow(i,5),khi(i,5),z,fjs,xtrpin) call splintslab(vdepth(1,1,jrs),fjm1,f2jm1,imt,km,i, * klow(i,6),khi(i,6),z,fs,xtrpout) if (xtrpin.or.xtrpout) fjs=fs z=p5*(tdepth(i,k,jrn)+tdepth(i,k,jrs)) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i, * klow(i,7),khi(i,7),z,fjw,xtrpin) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i-1, * klow(i,8),khi(i,8),z,fw,xtrpout) if (xtrpin.or.xtrpout) fjw=fw zbot=vdepth(i,k,jrn)+p5*dzvqz(i,k,0) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i, * klow(i,9),khi(i,9),zbot,fb,xtrpin) call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i+1, * klow(i,10),khi(i,10),zbot,fbe,xtrpout) if (xtrpin.or.xtrpout) fbe=fb call splintslab(vdepth(1,1,jrn),fj,f2j,imt,km,i-1, * klow(i,11),khi(i,11),zbot,fbw,xtrpout) if (xtrpin.or.xtrpout) fbw=fb call splintslab(vdepth(1,1,jrs),fjm1,f2jm1,imt,km,i, * klow(i,12),khi(i,12),zbot,fbs,xtrpout) if (xtrpin.or.xtrpout) fbs=fb call splintslab(vdepth(1,1,jrnp1),fjp1,f2jp1,imt,km,i, * klow(i,13),khi(i,13),zbot,fbn,xtrpout) if (xtrpin.or.xtrpout) fbn=fb # else z=p5*(tdepth(i,k,jrn)+tdepth(i+1,k,jrn)) call lintslab(vdepth(1,1,jrn),fj,imt,km,i, * klow(i,1),khi(i,1),z,fjn,xtrpin) call lintslab(vdepth(1,1,jrnp1),fjp1,imt,km,i, * klow(i,2),khi(i,2),z,fn,xtrpout) if (xtrpin.or.xtrpout) fjn=fn z=p5*(tdepth(i+1,k,jrn)+tdepth(i+1,k,jrs)) call lintslab(vdepth(1,1,jrn),fj,imt,km,i, * klow(i,3),khi(i,3),z,fje,xtrpin) call lintslab(vdepth(1,1,jrn),fj,imt,km,i+1, * klow(i,4),khi(i,4),z,fe,xtrpout) if (xtrpin.or.xtrpout) fje=fe z=p5*(tdepth(i,k,jrs)+tdepth(i+1,k,jrs)) call lintslab(vdepth(1,1,jrn),fj,imt,km,i, * klow(i,5),khi(i,5),z,fjs,xtrpin) call lintslab(vdepth(1,1,jrs),fjm1,imt,km,i, * klow(i,6),khi(i,6),z,fs,xtrpout) if (xtrpin.or.xtrpout) fjs=fs z=p5*(tdepth(i,k,jrn)+tdepth(i,k,jrs)) call lintslab(vdepth(1,1,jrn),fj,imt,km,i, * klow(i,7),khi(i,7),z,fjw,xtrpin) call lintslab(vdepth(1,1,jrn),fj,imt,km,i-1, * klow(i,8),khi(i,8),z,fw,xtrpout) if (xtrpin.or.xtrpout) fjw=fw zbot=vdepth(i,k,jrn)+p5*dzvqz(i,k,0) call lintslab(vdepth(1,1,jrn),fj,imt,km,i, * klow(i,9),khi(i,9),zbot,fb,xtrpin) call lintslab(vdepth(1,1,jrn),fj,imt,km,i+1, * klow(i,10),khi(i,10),zbot,fbe,xtrpout) if (xtrpin.or.xtrpout) fbe=fb call lintslab(vdepth(1,1,jrn),fj,imt,km,i-1, * klow(i,11),khi(i,11),zbot,fbw,xtrpout) if (xtrpin.or.xtrpout) fbw=fb call lintslab(vdepth(1,1,jrs),fjm1,imt,km,i, * klow(i,12),khi(i,12),zbot,fbs,xtrpout) if (xtrpin.or.xtrpout) fbs=fb call lintslab(vdepth(1,1,jrnp1),fjp1,imt,km,i, * klow(i,13),khi(i,13),zbot,fbn,xtrpout) if (xtrpin.or.xtrpout) fbn=fb #endif ftopx(i)=fbotx(i) ftopy(i)=fboty(i) fbotx(i)=(fbe-fbw)*dxu2r(i)*csr(j)*dxur(i)* * (cstr(j)*((tdepth(i+1,k,jrs)+p5*dzqz(i+1,k,0))- * (tdepth(i,k,jrs)+p5*dzqz(i,k,0))) * +cstr(j+1)*((tdepth(i+1,k,jrn)+p5*dzqz(i+1,k,1))- * (tdepth(i,k,jrn)+p5*dzqz(i,k,1)))) fboty(i)=(fbn-fbs)*dyu2r(i)*dyur(i)* * (((tdepth(i+1,k,jrn)+p5*dzqz(i+1,k,1))- * (tdepth(i+1,k,jrs)+p5*dzqz(i+1,k,0))) * +((tdepth(i,k,jrn)+p5*dzqz(i,k,1))- * (tdepth(i,k,jrs)+p5*dzqz(i,k,0)))) lap_xx(i,k)=(dxout(i)* * (dxin(i )*(fe-fje)*(dzqz(i+1,k,0)+dzqz(i+1,k,1))- * dxin(i-1)*(fjw-fw)*(dzqz(i ,k,0)+dzqz(i, k,1))) * +(fbotx(i)-ftopx(i))) * *dzv2rqz(i,k) lap_yy(i,k)=(dyout* * (dyinn*(fn-fjn)*(dzqz(i,k,1)+dzqz(i+1,k,1))- * dyins*(fjs-fs)*(dzqz(i,k,0)+dzqz(i+1,k,0))) * +(fboty(i)-ftopy(i))) * *dzv2rqz(i,k) lap_ug(i,k)=lpmtgd(j)*dxu2r(i)*((fje+fe)-(fjw+fw)) 30 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 40 k=1,km lap_xx(1 ,k)=lap_xx(imum1,k) lap_yy(1 ,k)=lap_yy(imum1,k) lap_ug(1 ,k)=lap_ug(imum1,k) lap_xx(imu,k)=lap_xx(2 ,k) lap_yy(imu,k)=lap_yy(2 ,k) lap_ug(imu,k)=lap_ug(2 ,k) 40 continue #endif return c c======================================================================= c Compute Laplacian of field at tracer points using their mask. c======================================================================= c entry lapt_depth(j,fjm1,fj,fjp1,mskm,msk,mskp,lap_xx,lap_yy) c c----------------------------------------------------------------------- c Grab grid spacings. c----------------------------------------------------------------------- c do 100 i=1,imt dxin(i)=dxur(i)*cstr(j) dxout(i)=dxtr(i)*cstr(j) 100 continue dyout=dytr(j)*cstr(j) dyins=dyur(j-1)*cs(j-1) dyinn=dyur(j)*cs(j) #ifdef cubspl c c----------------------------------------------------------------------- c Compute second derivatives for cubic splines. c----------------------------------------------------------------------- c call splineslab(tdepth(1,1,jrsm1),fjm1,imt,km,dfz1,dfzkm,wk,f2jm1) call splineslab(tdepth(1,1,jrs),fj,imt,km,dfz1,dfzkm,wk,f2j) call splineslab(tdepth(1,1,jrn),fjp1,imt,km,dfz1,dfzkm,wk,f2jp1) #endif c c----------------------------------------------------------------------- c Initialize search indices and bottom fluxes. c----------------------------------------------------------------------- c call init_srch(khi,klow,imt,km,13) c do 120 i=2,imum1 fbotx(i)=c0 fboty(i)=c0 120 continue c c----------------------------------------------------------------------- c Compute the laplacian in the slab. Note: extrapolations and c boundaries are set to give a no-flux boundary condition. c----------------------------------------------------------------------- c do 130 k=1,km do 130 i=2,imtm1 #ifdef cubspl z=(vdepth(i-1,k,jrn)*dxu(i-1)+vdepth(i,k,jrn)*dxu(i))*dxt2r(i) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i, * klow(i,1),khi(i,1),z,fjn,xtrpin) call splintslab(tdepth(1,1,jrn),fjp1,f2jp1,imt,km,i, * klow(i,2),khi(i,2),z,fn,xtrpout) if(xtrpin.or.xtrpout) fjn=fn z=(vdepth(i,k,jrn)*dyu(j)+vdepth(i,k,jrs)*dyu(j-1))*dyt2r(j) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i, * klow(i,3),khi(i,3),z,fje,xtrpin) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i+1, * klow(i,4),khi(i,4),z,fe,xtrpout) if(xtrpin.or.xtrpout) fje=fe z=(vdepth(i-1,k,jrs)*dxu(i-1)+vdepth(i,k,jrs)*dxu(i))*dxt2r(i) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i, * klow(i,5),khi(i,5),z,fjs,xtrpin) call splintslab(tdepth(1,1,jrsm1),fjm1,f2jm1,imt,km,i, * klow(i,6),khi(i,6),z,fs,xtrpout) if(xtrpin.or.xtrpout) fjs=fs z=(vdepth(i-1,k,jrn)*dyu(j)+vdepth(i-1,k,jrs)*dyu(j-1))* * dyt2r(j) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i, * klow(i,7),khi(i,7),z,fjw,xtrpin) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i-1, * klow(i,8),khi(i,8),z,fw,xtrpout) if(xtrpin.or.xtrpout) fjw=fw zbot=tdepth(i,k,jrs)+p5*dzqz(i,k,0) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i, * klow(i,9),khi(i,9),zbot,fb,xtrpin) call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i-1, * klow(i,10),khi(i,10),zbot,fbw,xtrpout) if(xtrpin.or.xtrpout) fbw=fb call splintslab(tdepth(1,1,jrs),fj,f2j,imt,km,i+1, * klow(i,11),khi(i,11),zbot,fbe,xtrpout) if(xtrpin.or.xtrpout) fbe=fb call splintslab(tdepth(1,1,jrsm1),fjm1,f2jm1,imt,km,i, * klow(i,12),khi(i,12),zbot,fbs,xtrpout) if(xtrpin.or.xtrpout) fbs=fb call splintslab(tdepth(1,1,jrn),fjp1,f2jp1,imt,km,i, * klow(i,13),khi(i,13),zbot,fbn,xtrpout) if(xtrpin.or.xtrpout) fbn=fb # else z=(vdepth(i-1,k,jrn)*dxu(i-1)+vdepth(i,k,jrn)*dxu(i))*dxt2r(i) call lintslab(tdepth(1,1,jrs),fj,imt,km,i, * klow(i,1),khi(i,1),z,fjn,xtrpin) call lintslab(tdepth(1,1,jrn),fjp1,imt,km,i, * klow(i,2),khi(i,2),z,fn,xtrpout) if(xtrpin.or.xtrpout) fjn=fn z=(vdepth(i,k,jrn)*dyu(j)+vdepth(i,k,jrs)*dyu(j-1))*dyt2r(j) call lintslab(tdepth(1,1,jrs),fj,imt,km,i, * klow(i,3),khi(i,3),z,fje,xtrpin) call lintslab(tdepth(1,1,jrs),fj,imt,km,i+1, * klow(i,4),khi(i,4),z,fe,xtrpout) if(xtrpin.or.xtrpout) fje=fe z=(vdepth(i-1,k,jrs)*dxu(i-1)+vdepth(i,k,jrs)*dxu(i))*dxt2r(i) call lintslab(tdepth(1,1,jrs),fj,imt,km,i, * klow(i,5),khi(i,5),z,fjs,xtrpin) call lintslab(tdepth(1,1,jrsm1),fjm1,imt,km,i, * klow(i,6),khi(i,6),z,fs,xtrpout) if(xtrpin.or.xtrpout) fjs=fs z=(vdepth(i-1,k,jrn)*dyu(j)+vdepth(i-1,k,jrs)*dyu(j-1))* * dyt2r(j) call lintslab(tdepth(1,1,jrs),fj,imt,km,i, * klow(i,7),khi(i,7),z,fjw,xtrpin) call lintslab(tdepth(1,1,jrs),fj,imt,km,i-1, * klow(i,8),khi(i,8),z,fw,xtrpout) if(xtrpin.or.xtrpout) fjw=fw zbot=tdepth(i,k,jrs)+p5*dzqz(i,k,0) call lintslab(tdepth(1,1,jrs),fj,imt,km,i, * klow(i,9),khi(i,9),zbot,fb,xtrpin) call lintslab(tdepth(1,1,jrs),fj,imt,km,i-1, * klow(i,10),khi(i,10),zbot,fbw,xtrpout) if(xtrpin.or.xtrpout) fbw=fb call lintslab(tdepth(1,1,jrs),fj,imt,km,i+1, * klow(i,11),khi(i,11),zbot,fbe,xtrpout) if(xtrpin.or.xtrpout) fbe=fb call lintslab(tdepth(1,1,jrsm1),fjm1,imt,km,i, * klow(i,12),khi(i,12),zbot,fbs,xtrpout) if(xtrpin.or.xtrpout) fbs=fb call lintslab(tdepth(1,1,jrn),fjp1,imt,km,i, * klow(i,13),khi(i,13),zbot,fbn,xtrpout) if(xtrpin.or.xtrpout) fbs=fb #endif ftopx(i)=fbotx(i) ftopy(i)=fboty(i) fbotx(i)=(fbe-fbw)*dxt2r(i)*cstr(j)*dxtr(i)* * (csr(j-1)*((vdepth(i,k,jrs)+p5*dzvqz(i,k,1))- * (vdepth(i-1,k,jrs)+p5*dzvqz(i-1,k,1))) * +csr(j)*((vdepth(i,k,jrn)+p5*dzvqz(i,k,0))- * (vdepth(i-1,k,jrn)+p5*dzvqz(i-1,k,0)))) fboty(i)=(fbn-fbs)*dyt2r(i)*dytr(i)* * (((vdepth(i-1,k,jrn)+p5*dzvqz(i-1,k,0))- * (vdepth(i-1,k,jrs)+p5*dzvqz(i-1,k,1))) * +((vdepth(i,k,jrn)+p5*dzvqz(i,k,0))- * (vdepth(i,k,jrs)+p5*dzvqz(i,k,1)))) lap_xx(i,k)=(dxout(i)*dytr(j)* * (dxin(i )*msk(i+1,k)*(fe-fje)* * (dzvqz(i ,k,0)*dyu(j)+dzvqz(i ,k,1)*dyu(j-1)) * -dxin(i-1)*msk(i-1,k)*(fjw-fw)* * (dzvqz(i-1,k,0)*dyu(j)+dzvqz(i-1,k,1)*dyu(j-1))) * +(fbotx(i)-ftopx(i))) * *dz2rqz(i,k,0) lap_yy(i,k)=(dyout*dxtr(i)* * (dyinn*mskp(i,k)*(fn-fjn)* * (dzvqz(i,k,0)*dxu(i)+dzvqz(i-1,k,0)*dxu(i-1)) * -dyins*mskm(i,k)*(fjs-fs)* * (dzvqz(i,k,1)*dxu(i)+dzvqz(i-1,k,0)*dxu(i-1))) * +(fboty(i)-ftopy(i))) * *dz2rqz(i,k,0) 130 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 140 k=1,km lap_xx(1 ,k)=lap_xx(imtm1,k) lap_yy(1 ,k)=lap_yy(imtm1,k) lap_xx(imt,k)=lap_xx(2 ,k) lap_yy(imt,k)=lap_yy(2 ,k) 140 continue #endif return end