subroutine shap_indx(npt,nxp,nyp,mask,isx,nfxk,nfpxk,nfyk,ifx,ifpx,ify) c------------------------------------------------------------------------- c nfx : number of contiguous latitudinal segments c ifx(1,i): starting point (x-compressed index) of i'th segment c ifx(2,i): length of i-th segment c ifx(3,i): highest order filter used in i'th segment c nfpx : number of segments with periodic overlap c ifx(1,i): starting point western portion of segment c ifx(2,i): length of western portion of segment c ifx(3,i): starting point eastern portion of segment c ifx(4,i): length of eastern portion of segment c ifx(5,i): highest order filter used in segment c nfy : number of contiguous longitudinal segments c ify(1,j): starting point (y-compressed index) of j'th segment c ify(2,j): length of j'th segment c ify(3,j): highest order filter used in j'th segment include "comm_new.h" include 'comm_para.h' dimension mask(nxp,1),isx(1),ifx(3,1),ifpx(5,1),ify(3,1) logical prev, curr, peri, ZER common /shap_c25/ cs25(10), ZER common /new_filt/ MAXFO, nbx, nby, nfx, nfpx, nfy do i = 1, MAXFO cs25(i) = 1./(4.**i) enddo nfxk = 0 nfpxk = 0 do irow = 1, nyp ista = 1 prev = (mask(1, irow) .eq. 0) peri = (.not. prev) .and. (iglob .ne. 0) do icol = 2, nxp curr = (mask(icol,irow) .eq. 0) if (curr .ne. prev) then if (peri .and. ista .eq. 1) then nfpxk = nfpxk + 1 ifpx(1,nfpxk) = mask(1,irow) ifpx(2,nfpxk) = icol - 1 elseif (curr) then nfxk = nfxk + 1 ifx(1,nfxk) = mask(ista,irow) ifx(2,nfxk) = icol-ista ifx(3,nfxk) = min((icol-ista-1)/2,MAXFO) endif ista = icol prev = curr endif enddo if (.not. curr) then if (peri) then if (ista .eq. 1) then nfpxk = nfpxk + 1 ifpx(1,nfpxk) = mask(1,irow) ifpx(2,nfpxk) = nxp ifpx(3,nfpxk) = 0 ifpx(4,nfpxk) = 0 ifpx(5,nfpxk) = MAXFO else ifpx(3,nfpxk) = mask(ista,irow) ifpx(4,nfpxk) = mask(nxp,irow) - mask(ista,irow) + 1 ifpx(5,nfpxk) = min((ifpx(2,nfpxk)+ifpx(4,nfpxk)-1)/2,MAXFO) endif else nfxk = nfxk + 1 ifx(1,nfxk) = mask(ista,irow) ifx(2,nfxk) = nxp-ista+1 ifx(3,nfxk) = min((nxp-ista)/2,MAXFO) endif endif enddo nfyk = 0 do icol = 1, nxp jsta = 1 prev = (mask(icol, 1) .eq. 0) do irow = 2, nyp curr = (mask(icol,irow) .eq. 0) if (curr .ne. prev) then if (curr) then nfyk = nfyk + 1 ify(1,nfyk) = isx(mask(icol,jsta)) ify(2,nfyk) = irow-jsta ify(3,nfyk) = min((irow-jsta-1)/2,MAXFO) endif jsta = irow prev = curr endif enddo if (.not. curr) then nfyk = nfyk + 1 ify(1,nfyk) = isx(mask(icol,jsta)) ify(2,nfyk) = nyp-jsta+1 ify(3,nfyk) = min((nyp-jsta)/2,MAXFO) endif enddo return end c----------------------------------------------------------------------- subroutine shap_vec(nstep,npt,nz,uc,vc,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) c----------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension uc(npt,1),vc(npt,1),tp(1) dimension lxxk(MXBDY,nz),lyxk(MAXNB,nz),isyk(npt,nz),isk(npt,nz), * ifxk(9*MAXSID,nz), ifpxk(5*MAXSID,nz), ifyk(9*MAXSID,nz) common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) logical NOSLIP common /new_shap/ nordu, nordh, mshx, mshy, mshh common/shapi/nord,nshapu,nshaph,nshapt,mshapu,mshaph,mshapt common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20), + itherm,mlc,limp common /new_filt/ MAXFO, nbx, nby, nfx, nfpx, nfy external shap_1do, shap_1dco if (nstep.eq.0 .or. (nshapu.ne.0 .and. mod(nstep,nshapu).eq.0)) then NOSLIP = (mbc .eq. 1) do k = 1, nz npk = nptk(k) nbx = nbxk(k) nby = nbyk(k) nfx = nfxk(k) nfpx = nfpxk(k) nfy = nfyk(k) call shap_2d (nordu, mshapu, .true.,NOSLIP, npk, * lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k), * ifxk(1,k),ifpxk(1,k),ifyk(1,k), uc(1,k),tp,shap_1do) call shap_2d (nordu, mshapu, NOSLIP,.true., npk, * lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k), * ifxk(1,k),ifpxk(1,k),ifyk(1,k), vc(1,k),tp,shap_1do) enddo endif return end c----------------------------------------------------------------------- subroutine shap_scl(nstep,npt,nz,tem,lxxk,lyxk,isyk,isk,ifxk,ifpxk,ifyk,tp) c----------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension tem(npt,nz), tp(1) dimension lxxk(MXBDY,nz),lyxk(MAXNB,nz),isyk(npt,nz),isk(npt,nz), * ifxk(9*MAXSID,nz), ifpxk(5*MAXSID,nz), ifyk(9*MAXSID,nz) common /new_shap/ nordu, nordh, mshx, mshy, mshh 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) common /new_filt/ MAXFO, nbx, nby, nfx, nfpx, nfy external shap_1do, shap_1dco if (nstep.eq.0 .or. (nshaph.ne.0 .and. mod(nstep,nshaph).eq.0)) then do k = 1, nz npk = nptk(k) nbx = nbxk(k) nby = nbyk(k) nfx = nfxk(k) nfpx = nfpxk(k) nfy = nfyk(k) call shap_3d (nordh, mshaph, .false.,.false.,npk, * lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k), * ifxk(1,k),ifpxk(1,k),ifyk(1,k), tem(1,k),tp,shap_1dco) enddo endif return end c-------------------------------------------------------------------------- subroutine shap_2d (nord,key,SETX,SETY,npt,lxx,lyx,isy,isk, * ifx,ifpx,ify, aa,abc,filter) c-------------------------------------------------------------------------- c.....if key.eq.0 then reduce order for all points on short segments c since ifx(3,*) = min((ifx(2,*)-1)/2,MAXFO) : c length=4 -> nshap = 1 c length=5,6 -> nshap = min(2,nord) c length=7,8 -> nshap = min(3,nord) c etc. implicit real(a-h,o-z),integer(i-n) dimension lxx(1),lyx(1),isy(1),isk(1), aa(1), bb(1), cc(1), abc(npt,1) dimension ifx(3,1), ifpx(5,1), ify(3,1) logical SETX, SETY, ZER, REDUC pointer (pbb, bb), (pcc, cc) external filter common /shap_c25/ cs25(10), ZER common /new_filt/ MAXFO, nbx, nby, nfx, nfpx, nfy c sizeof(abc) = 3*npt pbb = loc(abc(1,1)) pcc = loc(abc(1,2)) c REDUC = (key .eq. 0) c.......apply Shapiro filter in X direction do i = 1, npt bb(i) = aa(isk(i)) enddo ZER = SETX ! used only if filter = shap_1do nshap = MAXFO do k = 1, nfx if (REDUC) nshap = min(nord, ifx(3,k)) call filter (nshap, ifx(2,k), bb(ifx(1,k)), cc) enddo do k = 1, nfpx if (REDUC) nshap = min(nord, ifpx(5,k)) if (ifpx(4,k) .eq. 0) then call shap_1dp0 (nshap, ifpx(2,k), bb(ifpx(1,k)), cc) else call shap_1dper (nshap, ifpx(1,k), bb, cc, filter) endif enddo do i = 1, npt cc(isk(i)) = bb(i) enddo if (SETX) then do i = 1, nbx cc(lxx(i)) = 0. enddo endif if (SETY) then do i = 1, nby cc(lyx(i)) = 0. enddo endif do j = 1, npt i = isk(j) aa(i) = aa(i) - cc(i) enddo c.......apply Shapiro filter in Y direction do i = 1, npt bb(i) = aa(isy(i)) enddo ZER = SETY ! used only if filter = shap_1do do k = 1, nfy if (REDUC) nshap = min(nord, ify(3,k)) call filter (nshap, ify(2,k), bb(ify(1,k)), cc) enddo do i = 1, npt cc(isy(i)) = bb(i) enddo if (SETX) then do i = 1, nbx cc(lxx(i)) = 0. enddo endif if (SETY) then do i = 1, nby cc(lyx(i)) = 0. enddo endif do j = 1, npt i = isk(j) aa(i) = aa(i) - cc(i) enddo return end subroutine shap_1do (nshap, nn, f, tmp) c---------------------------------------------- c non-conservative, order reduced near boundaries c ZER => no change at boundary c (i.e. u=0 before filter -> u=0 after filter) implicit real(a-h,o-z),integer(i-n) dimension f(1), tmp(1), aa(1), bb(1), cc(1) pointer (paa, aa), (pbb, bb), (pcc, cc) logical ZER common /shap_c25/ cs25(10), ZER paa = loc(f(1)) pbb = loc(tmp(1)) pcc = loc(tmp(nn+1)) do n = 1, nshap const = cs25(n) iab = paa paa = pbb pbb = iab f0 = bb(1) fp = bb(2) do i = 2, nn-1 fm = f0 f0 = fp fp = bb(i+1) aa(i) = (f0 - fp) + (f0 - fm) enddo if (n .eq. 1) then if (ZER) then aa(1) = 0. aa(nn) = 0. else aa(1) = bb(1) - bb(2) aa(nn) = bb(nn) - bb(nn-1) endif endif cc(n) = const*aa(n) cc(nn-n+1) = const*aa(nn-n+1) enddo do i = nshap+1, nn-nshap f(i) = const*aa(i) enddo do i = 1, nshap f(i) = cc(i) f(nn-i+1) = cc(nn-i+1) enddo return end subroutine shap_1dper (nshap, id, f, tmp, filter) c------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension id(4), f(1), tmp(1), aa(1), bb(1) pointer (paa, aa), (pbb, bb) external filter i1 = id(1) n1 = id(2) i2 = id(3) n2 = id(4) n12 = n1 + n2 paa = loc(f(i2)) do i = 1, n2 tmp(i) = aa(i) enddo paa = loc(f(i1)) k = n2 do i = 1, n1 k = k + 1 tmp(k) = aa(i) enddo call filter (nshap, n12, tmp, tmp(n12+1)) k = n2 do i = 1, n1 k = k + 1 aa(i) = tmp(k) enddo paa = loc(f(i2)) do i = 1, n2 aa(i) = tmp(i) enddo return end subroutine shap_1dco (nshap, nn, f, tmp) c--------------------------------------------- c conservative, order NOT reduced near boundaries implicit real(a-h,o-z),integer(i-n) dimension f(1), tmp(1), aa(1), bb(1) pointer (paa, aa), (pbb, bb) logical ZER common /shap_c25/ cs25(10), ZER paa = loc(f(1)) pbb = loc(tmp(1)) do n = 1, nshap iab = paa paa = pbb pbb = iab f0 = bb(1) fp = bb(2) do i = 2, nn-1 fm = f0 f0 = fp fp = bb(i+1) aa(i) = (f0 - fp) + (f0 - fm) enddo aa(1) = bb(1) - bb(2) aa(nn) = bb(nn) - bb(nn-1) enddo const = cs25(nshap) do i = 1, nn f(i) = const * aa(i) enddo return end subroutine shap_1dp0 (nshap, nn, f, tmp) c--------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension f(1), tmp(1), aa(1), bb(1) pointer (paa, aa), (pbb, bb) logical ZER common /shap_c25/ cs25(10), ZER paa = loc(f(1)) pbb = loc(tmp(1)) do n = 1, nshap iab = paa paa = pbb pbb = iab f1 = bb(1) f2 = bb(2) f0 = f1 fp = f2 do i = 2, nn-1 fm = f0 f0 = fp fp = bb(i+1) aa(i) = (f0 - fp) + (f0 - fm) enddo aa(1) = (f1 - f2) + (f1 - fp) aa(nn) = (fp - f1) + (fp - f0) enddo const = cs25(nshap) do i = 1, nn f(i) = const * aa(i) enddo return end c-------------------------------------------------------------------------- subroutine shap_3d (nord,key,SETX,SETY,npt,lxx,lyx,isy,isk, * ifx,ifpx,ify, aa,abc,filter) c-------------------------------------------------------------------------- c.....reduces order for all points on short segments c length=4 -> nshap = 1 c length=5,6 -> nshap = min(2,nord) c length=7,8 -> nshap = min(3,nord) c etc. c.....but if nshap <= LIM_SHAP, then DON'T FILTER segment at all implicit real(a-h,o-z),integer(i-n) dimension lxx(1),lyx(1),isy(1),isk(1), aa(1), bb(1), cc(1), abc(npt,1) dimension ifx(3,1), ifpx(5,1), ify(3,1) logical SETX, SETY, ZER, REDUC pointer (pbb, bb), (pcc, cc) external filter common /shap_c25/ cs25(10), ZER common /new_filt/ MAXFO, nbx, nby, nfx, nfpx, nfy LIM_SHAP = key - 2 c sizeof(abc) = 3*npt pbb = loc(abc(1,1)) pcc = loc(abc(1,2)) c.....use a reduced order filter in a narrow passages. c c.......apply Shapiro filter in X direction do i = 1, npt bb(i) = aa(isk(i)) enddo ZER = SETX do k = 1, nfx nshap = min(nord, ifx(3,k)) if (nshap .gt. LIM_SHAP) then call filter (nshap, ifx(2,k), bb(ifx(1,k)), cc) else call zero_em (ifx(2,k), bb(ifx(1,k))) endif enddo do k = 1, nfpx nshap = min(nord, ifpx(5,k)) if (nshap .gt. LIM_SHAP) then if (ifpx(4,k) .eq. 0) then call shap_1dp0 (nshap, ifpx(2,k), bb(ifpx(1,k)), cc) else call shap_1dper (nshap, ifpx(1,k), bb, cc, filter) endif else call zero_em (ifpx(2,k), bb(ifpx(1,k))) call zero_em (ifpx(4,k), bb(ifpx(3,k))) endif enddo do i = 1, npt cc(isk(i)) = bb(i) enddo if (SETX) then do i = 1, nbx cc(lxx(i)) = 0. enddo endif if (SETY) then do i = 1, nby cc(lyx(i)) = 0. enddo endif do j = 1, npt i = isk(j) aa(i) = aa(i) - cc(i) enddo c.......apply Shapiro filter in Y direction do i = 1, npt bb(i) = aa(isy(i)) enddo ZER = SETY do k = 1, nfy nshap = min(nord, ify(3,k)) if (nshap .gt. LIM_SHAP) then call filter (nshap, ify(2,k), bb(ify(1,k)), cc) else call zero_em (ify(2,k), bb(ify(1,k))) endif enddo do i = 1, npt cc(isy(i)) = bb(i) enddo if (SETX) then do i = 1, nbx cc(lxx(i)) = 0. enddo endif if (SETY) then do i = 1, nby cc(lyx(i)) = 0. enddo endif do j = 1, npt i = isk(j) aa(i) = aa(i) - cc(i) enddo return end