subroutine new_topo (nxp,nyp,nz,npt,zin,dzin,hin,nsig,sigma,dept,h,nptk,nzi) dimension zin(nz+1),dzin(nz+1),dept(npt),nptk(nz),h(npt,nz) * ,nptk(nz),nzi(npt),hin(nz),sigma(nz) include 'comm_new.h' do k = 1, nz do i = 1, npt c h(i,k) = -98765432. h(i,k) = 0. enddo enddo do i = 1, npt dep = dept(i) n = 1 do k = 2, nz if (dep.lt.zin(k)) goto 10 n = k enddo 10 nzi(i) = n enddo if (initbt.eq.0) then do i = 1, npt mz = nzi(i) do k = 1, mz h(i,k) = dzin(k+1) + dzin(k) enddo enddo elseif (initbt.eq.3) then mmz = 0 do i = 1, npt mmz = max(nzi(i),mmz) enddo do i = 1, npt mz = nzi(i) do k = 1, mz - 1 h(i,k) = dzin(k+1) + dzin(k) enddo if (mz.eq.mmz) then h(i,mz) = dzin(mz) + (dept(i)-zin(mz)) else h(i,mz) = dzin(mz+1) + dzin(mz) endif enddo else do i = 1, npt dep = dept(i) mz = nzi(i) do k = 1, mz - 1 h(i,k) = dzin(k+1) + dzin(k) enddo h(i,mz) = dzin(mz) + (dep-zin(mz)) enddo endif hin(1) = (zin(2) + zin(1))/2. do i = 2, nz - 1 hin(i) = (zin(i+1) - zin(i-1))/2. enddo hin(nz) = dept(i) - (zin(nz) + zin(nz-1))/2. return end c --------------------------------------------------------------------- subroutine data_init (npt,nptk,nz,isk,u,v,uc,vc,fu,fv,ft,fsal,bdiv,ubar,vbar,use_salt) c --------------------------------------------------------------------- dimension nptk(nz),isk(npt,nz),bdiv(npt),ubar(npt),vbar(npt), * u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz), * fu(npt,nz), fv(npt,nz), ft(npt,nz), fsal(npt,nz) logical use_salt do i = 1, npt bdiv(i) = 0. ubar(i) = 0. vbar(i) = 0. enddo do k = 1, nz do j = 1, nptk(k) i = isk(j,k) u(i,k) = 0. v(i,k) = 0. uc(i,k) = 0. vc(i,k) = 0. fu(i,k) = 0. fv(i,k) = 0. ft(i,k) = 0. enddo enddo if (.not.use_salt) return do k = 1, nz do j = 1, nptk(k) i = isk(j,k) fsal(i,k) = 0. enddo enddo return end c --------------------------------------------------------------------- subroutine dfdxk (f,df,npt,npk,nbu,nbx,nby,ncs,lxx,lyx,snx * ,npbc,lpbcw,lpbce,isk) c --------------------------------------------------------------------- c............ a version with Periodic B.C. (Senya Basin) implicit real(a-h,o-z),integer(i-n) dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4), lyx(nby) dimension lpbcw(npbc), lpbce(npbc),isk(npk) c compute centered first derivative for entire grid. #ifdef fourth_order do i = 3, npk - 2 j = isk(i) df(j) = (8.*(f(j+1)-f(j-1)) - (f(j+2)-f(j-2)))/12. enddo #else !second_order do i = 2, npk - 1 j = isk(i) df(j) = (f(j+1)-f(j-1))/2. enddo #endif c....................... periodic B.C. do i = 1, npbc i2 = lpbce(i) f4 = f(i2) f3 = f(i2-1) f2 = f(i2-2) f1 = f(i2-3) i1 = lpbcw(i) f5 = f(i1) f6 = f(i1+1) f7 = f(i1+2) f8 = f(i1+3) #ifdef fourth_order df(i1) = (8.*(f6 - f4) - (f7 - f3))/12. df(i1+1) = (8.*(f7 - f5) - (f8 - f4))/12. df(i2) = (8.*(f5 - f3) - (f6 - f2))/12. df(i2-1) = (8.*(f4 - f2) - (f5 - f1))/12. #else !second_order df(i1) = (f6 - f4)/2. df(i1+1) = (f7 - f5)/2. df(i2) = (f5 - f3)/2. df(i2-1) = (f4 - f2)/2. #endif enddo nb = nbx if (nbu .eq. -1 .or. nbu .eq. 2) then do i = 1, nb i1 = lxx(i,1) i2 = lxx(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lxx(i,3)) f4 = f(lxx(i,4)) df(i1) = 0. #ifdef fourth_order df(i2) = snx(i)*( 2.*(f2-f1) + 5.*(f3-f2) - (f4-f3))/6. #endif enddo else do i = 1, nb i1 = lxx(i,1) i2 = lxx(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lxx(i,3)) f4 = f(lxx(i,4)) #ifdef fourth_order df(i1) = snx(i)*(11.*(f2-f1) + 7.*(f2 - f3) + 2.*(f4-f3))/6. df(i2) = snx(i)*( 2.*(f2-f1) + 5.*(f3-f2) - (f4-f3))/6. #else df(i1) = snx(i)*(3.*(f2 - f1) + (f2 - f3))/2. #endif enddo endif if (nbu .eq. -1 .or. nbu.eq.1) then c..................set the derivative along the boundary equal to zero. do i = 1, nby df(lyx(i)) = 0. enddo endif return c end of dfdxk. end c ------------------------------------------------------------------ subroutine dfdyk(f,df,npt,npk,nbv,nby,nbx,ncs,lyy,lxy,sny,isy) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4),lxy(nbx) * ,isy(npk) c note, isy: k-y-comp -> x-comp, c lyy: k-y-comp-bound -> x-comp , etc. c #ifdef fourth_order do i = 3, npk-2 j = isy(i) df(j)=(8.*(f(isy(i+1))-f(isy(i-1))) - (f(isy(i+2))-f(isy(i-2))))/12. enddo #else !second_order do i = 2, npk-1 j = isy(i) df(j)=(f(isy(i+1))-f(isy(i-1)))/2. enddo #endif nb = nby if (nbv.eq.-1 .or. nbv.eq.2) then do i = 1, nb i1 = lyy(i,1) i2 = lyy(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lyy(i,3)) f4 = f(lyy(i,4)) df(i1) = 0. #ifdef fourth_order df(i2) = sny(i)*(2.*(f2-f1) + 5.*(f3-f2) - (f4-f3))/6. #endif enddo else do i = 1, nb i1 = lyy(i,1) i2 = lyy(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lyy(i,3)) f4 = f(lyy(i,4)) #ifdef fourth_order df(i1) = sny(i)*(11.*(f2-f1) + 7.*(f2-f3) + 2.*(f4-f3))/6. df(i2) = sny(i)*( 2.*(f2-f1) + 5.*(f3-f2) - (f4-f3))/6. #else !second_order df(i1) = sny(i)*(3.*(f2-f1) + (f2-f3))/2. #endif enddo endif if(nbv.eq.-1 .or. nbv.eq.1) then c.....................set the derivative along the x boundary equal to zero. do i = 1, nbx df(lxy(i)) = 0. enddo endif return c end of dfdyk. end c --------------------------------------------------------------------- subroutine dfdx1 (f,df,npt,nbu,nbx,nby,ncs,lxx,lyx,snx * ,npbc,lpbcw,lpbce) c --------------------------------------------------------------------- c............ a version with Periodic B.C. (Senya Basin) implicit real(a-h,o-z),integer(i-n) dimension f(npt),df(npt),snx(nbx),lxx(nbx+ncs,4), lyx(nby) dimension lpbcw(npbc), lpbce(npbc) c compute fourth order centered first derivative for entire grid. #ifdef fourth_order do j = 3, npt - 2 df(j) = (8.*(f(j+1)-f(j-1)) - (f(j+2)-f(j-2)))/12. enddo #else !second_order do j = 2, npt - 1 df(j) = (f(j+1)-f(j-1))/2. enddo #endif c....................... periodic B.C. do i = 1, npbc i2 = lpbce(i) f4 = f(i2) f3 = f(i2-1) f2 = f(i2-2) f1 = f(i2-3) i1 = lpbcw(i) f5 = f(i1) f6 = f(i1+1) f7 = f(i1+2) f8 = f(i1+3) #ifdef fourth_order df(i1) = (8.*(f6 - f4) - (f7 - f3))/12. df(i1+1) = (8.*(f7 - f5) - (f8 - f4))/12. df(i2) = (8.*(f5 - f3) - (f6 - f2))/12. df(i2-1) = (8.*(f4 - f2) - (f5 - f1))/12. #else !second_order df(i1) = (f6 - f4)/2. df(i1+1) = (f7 - f5)/2. df(i2) = (f5 - f3)/2. df(i2-1) = (f4 - f2)/2. #endif enddo do i = 1, nbx i1 = lxx(i,1) i2 = lxx(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lxx(i,3)) f4 = f(lxx(i,4)) #ifdef fourth_order df(i1) = snx(i)*(3.*(f2-f1) + (f2-f3))/2. c df(i1) = snx(i)*(11.*(f2-f1) + 7.*(f2 - f3) + 2.*(f4-f3))/6. df(i2) = snx(i)*(2.*(f2-f1) + 5.*(f3-f2) + (f3-f4))/6. #else df(i1) = snx(i)*(f2 - f1) #endif enddo if (nbu.eq.1) then c..................set the derivative along the boundary equal to zero. do i = 1, nby df(lyx(i)) = 0. enddo endif return end c ------------------------------------------------------------------ subroutine dfdy1(f,df,npt,nbv,nby,nbx,ncs,lyy,lxy,sny,isy) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension f(npt),df(npt),sny(nby),lyy(nby+ncs,4),lxy(nbx) * ,isy(npt) c note, isy: k-y-comp -> x-comp, c lyy: k-y-comp-bound -> x-comp , etc. c #ifdef fourth_order do i = 3, npt-2 j = isy(i) df(j)=(8.*(f(isy(i+1))-f(isy(i-1))) - (f(isy(i+2))-f(isy(i-2))))/12. enddo #else !second_order do i = 2, npt-1 j = isy(i) df(j)=(f(isy(i+1))-f(isy(i-1)))/2. enddo #endif do i = 1, nby i1 = lyy(i,1) i2 = lyy(i,2) f1 = f(i1) f2 = f(i2) f3 = f(lyy(i,3)) f4 = f(lyy(i,4)) #ifdef fourth_order df(i1) = sny(i)*(3.*(f2-f1) + (f2-f3))/2. c df(i1) = sny(i)*(11.*(f2-f1) + 7.*(f2-f3) + 2.*(f4-f3))/6. df(i2) = sny(i)*(2.*(f2-f1) + 5.*(f3-f2) + (f3-f4))/6. #else df(i1) = sny(i)*(f2 - f1) #endif enddo if(nbv.eq.1) then c.....................set the derivative along the x boundary equal to zero. do i = 1, nbx df(lxy(i)) = 0. enddo endif return c end of dfdy1. end subroutine baro_dept(npt,nz,nzi,nzi_b,h,lxxk,lyyk,mbc,dept,nz_x,nz_y) c --------------------------------------------------------------------- c...... requires no-slip boundary condition c...... if point is an x or y boundary point on the k-th grid, but not c...... the first grid, then it must be treated differently in baro_shap c...... and baro_updat in order to satisfy boundary conditions implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension h(npt,nz),dept(npt), * lxxk(MXBDY,nz),lyyk(MXBDY,nz),nzi(npt),nzi_b(npt) * ,nz_x(npt),nz_y(npt) common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) do i = 1, npt nz_x(i) = nzi(i) nz_y(i) = nzi(i) enddo c find x-direction restrictions on depth do k = nz, 2, -1 do j = 1, nbxk(k) i = lxxk(j,k) nz_x(i) = k-1 enddo enddo k = 1 do j = 1, nbxk(k) i = lxxk(j,k) nz_x(i) = nzi(i) enddo do i = 1, npt nzi_b(i) = nz_x(i) enddo #ifdef fourth_order do k = nz, 2, -1 nxck = nbxk(k) + ncsk(k) do i = 1, nbxk(k) j = lxxk(i+nxck,k) nz_x(j) = k-1 enddo enddo k = 1 nxck = nbxk(k) + ncsk(k) do i = 1, nbxk(k) j = lxxk(i+nxck,k) nz_x(j) = nzi_b(j) enddo #endif c find y-direction restrictions on depth do k = nz, 2, -1 do j = 1, nbyk(k) i = lyyk(j,k) nz_y(i) = k-1 enddo enddo k = 1 do j = 1, nbyk(k) i = lyyk(j,k) nz_y(i) = nzi(i) enddo do i = 1, npt nzi_b(i) = nz_y(i) enddo #ifdef fourth_order do k = nz, 2, -1 nyck = nbyk(k) + ncsk(k) do i = 1, nbyk(k) j = lyyk(i+nyck,k) nz_y(j) = k-1 enddo enddo k = 1 nyck = nbyk(k) + ncsk(k) do i = 1, nbyk(k) j = lyyk(i+nyck,k) nz_y(j) = nzi_b(j) enddo #endif do i = 1, npt nzi_b(i) = min(nz_x(i),nz_y(i)) enddo do i = 1, npt dept(i) = h(i,1) do k = 2, nzi_b(i) dept(i) = dept(i) + h(i,k) enddo enddo return end subroutine baro_shap (nstep,npt,nz,nzi,nzi_b,dept,h,uc,vc,ubar,vbar,u,v, * lxxk,lyyk) c--------------------------------------------------------------------------- c This subroutine is responsible for preserving the rigid lid c assumption. It must be called immediately before ddiv to c work properly c Note that in the non-constant depth scenario, we need c div(sum(uc_k)) = sum(div(uc_k)). c - We need the same divergence operators on both sides. c We accomplish this by using the k=1 divergence operator c on all levels, assuming uc(i,k)=0 at all 'mud' (non-ocean) points. c - In order for the divergence and summation (over k) operators c to commute, we need also need the summation to be independent of c horizontal position, hence the sum must be over all 'nz' levels c but no normal vertical flow through the bottom requires that c div(uc_k) = 0 at all 'mud' points. c Thus we have two constraints while computing divergences in ddiv: c uc_k=0 and div(uc_k)=0 at mud points c NOTE: A distinction is made between 'mud' points and 'land' c points. A 'mud' point is by definitions a point which is c not water on some level k>1, but is water on the k=1 level. c These 'fixes' are not done on the land points, only mud points. c Zero mudpoint transport (uc_k = 0) is done in ddiv. c Zero mudpoint divergence (div(uc_k) = 0) is accomplished by c enforcing the normal component of velocity to be zero at c one(two) grid point(s) adjacent to land in the case of c second(fourth) order approximations to the derivatives. This c is done here and in bcset, but in order for the total transport to c remain fixed, we must redistribute the zero-ed out barotropic c transport in this routine. The redistribution of ubar(vbar) is c determined by nzi_b. This restores the correct barotropic transport c (which was spoiled by bcset and the shapiro filter). c----NHN:Dec.21,95 implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension h(npt,nz),uc(npt,nz),vc(npt,nz),ubar(npt),vbar(npt) dimension u(npt,nz),v(npt,nz) dimension dept(npt),nzi(npt),nzi_b(npt) dimension lxxk(MXBDY,nz), lyyk(MXBDY,nz) common/shapi/nord,nshapu,nshaph,nshapt,mshapu,mshaph,mshapt common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) do i = 1, npt do k = 1, nzi_b(i) hi = h(i,k) uc(i,k) = u(i,k)*hi vc(i,k) = v(i,k)*hi enddo enddo do i = 1, npt do k = nzi_b(i)+1, nzi(i) uc(i,k) = 0. vc(i,k) = 0. enddo enddo call baro_sum (npt, nz, nzi_b, uc, vc, u, v) call baro_scale (npt, u, v, dept) do i = 1, npt do k = 1, nzi_b(i) uc(i,k) = uc(i,k) + h(i,k)*(ubar(i) - u(i,1)) enddo do k = 1, nzi_b(i) vc(i,k) = vc(i,k) + h(i,k)*(vbar(i) - v(i,1)) enddo enddo call decap (npt, nz, nzi, u,v,uc,vc,h) return end