c*********************************************************************** c c routines for irregular geometry. c c*********************************************************************** c ------------------------------------------------------------------ subroutine bndrys(npt,iox,ioy,isxk,isyk,mask,h,nzi, * ixk,iyk,lxxk,lyyk,lxyk,lyxk,snxk,snyk,lok, * lpbcwk,lpbcek,ifxk,ifpxk,ifyk,dept) c ------------------------------------------------------------------ c finds the indices of the land-ocean boundary grid points. c c from the common block grid: c nxp,nyp = (input) # of grid points in the x and y directions. c nxyc = (input) # of ocean grid points. c mxbdy = (input) maximum storage space for each array lxx, lyy. c maxnb = (input) max storage space for each array lxy, lyx, snx, sny. c iox = (output) nxyc indices of the x-sorted ocean grid points. c ioy = (output) nxyc indices of the ocean points for a y-sort. c isx = (output) nxyc indices to gather data to the compressed c x-sort from the compressed y-sort. c isy = nxyc indices to gather data to the compressed y-sort c from the compressed x-sort. c lxx = MINSEG*(nbx+ncs) indices of points on and next to x-bndrys c for the compressed x-sort. c lyy = MINSEG*(nby+ncs) x-sort indices of points on and next c to y-bndrys for the compressed y-sort. c lxy = nbx+ncs indices of x-boundary points for the compressed c y-sort. c lyx = nby+ncs indices of y-boundary points for the compressed c x-sort. c c snx = (output) nbx+ncs signs (+1. or -1.) associated with the c x-boundary indices in lxx. if snx(i)=1., then the ocean c grid point with index lxx(i) is to the east (pos. x-dir.) c of the adjacent land point. if snx(i)=-1., c then the ocean pt. lxx(i) is to the west (neg. x-dir) of c the adjacent land grid point. c sny = (output) nby+ncs signs (+1. or -1.) associated with the c y-boundary indices in lyy in the same since as snx and lxx. c nbx = (output) # of x-boundary grid points. c nby = (output) # of y-boundary grid points. c ncs = (output) # of interior corner boundary grid points. c c ***************************************************************** c a few words about the different data sorts and bndry indices: c c the boundary grid points are the ocean grid points which are c adjacent to one or more land grid points. an ocean point is an c x-boundary point if the adjacent land point is in the pos. or neg. c x-direction, and an ocean point is a y-boundary pt if the adjacent c land point is in the pos. or neg. y-direction. thus a point c can be both an x-boundary and a y-boundary point. there are nbx c and nby x and y boundary points, respectively. a special case is c an ocean point which is neither an x nor a y boundary point, c but is adjacent to a land point in a diagonal direction. these are c referred to as interior corner boundary points (there are always c some on a peninsula). there are ncs ocean points of this type. c c assume that a data field is stored in the matrix u(i,j) whose c first dimension has been declared as nxp elements. the index i c corresponds to increasing x, and j corresponds to increasing y. in c fortran, the elements of u are stored sequentially with the index c i increasing before j, ie. u(1,1),u(2,1),...,u(nxp,1),u(1,2),... c this is the x or regular x-sort. the y or regular y-sort is simply c a sequential storage of u with the j index increasing before i, c u(1,1),u(1,2),...,u(1,nyp),u(2,1),u(2,2),.... the compressed c x-sort is the regular x-sort, excluding land points. the elements c representing ocean points are simply shifted towards u(1,1) to c fill the gaps left by the land elements, so that all ocean points c are stored consecutively. similarly, the compressed y-sort is the c regular y-sort, excluding land points. c for each type of sort the data u can obviously be identified by c a single index, say k. if k is simply the sequential position of c the element as it is stored, then for the regular sorts, k can be c expressed in terms of the original indices i and j: c k = i + (j-1)*nxp regular x-sort c k = j + (i-1)*nyp regular y-sort c the indices iox are the k-indices of the ocean pts for an x-sort c and the indices ioy are the k-indices of the ocean pts for a c regular y-sort. some relations between sorts and indices are: c uxc(i) = ux(iox(i)) compressed x-sort from regular x-sort c uyc(i) = uy(ioy(i)) compressed y-sort from regular y-sort c uyc(i) = uxc(isy(i)) compressed y-sort from compressed x-sort, c where isy is a columnwise ordering of the c compressed x-sort indices c uxc(i) = uyc(isx(i)) compressed x-sort from compressed y-sort, c where isx is a rowwise ordering of the c compressed y-sort indices c isx(isy(k)) = k, and isy(isx(k)) = k c c the k-index for the regular x-sort can be expressed in terms c of the k-index for the regular y-sort, and vice versa. so the c indices of ocean grid points in the two sorts are related as: c ioy(i) = ((iox(i)-1)/nxp)*(1-nxp*nyp) + (iox(i)-1)*nyp + 1 c and c iox(i) = ((ioy(i)-1)/nyp)*(1-nxp*nyp) + (ioy(i)-1)*nxp + 1 c c (the divisions must operate on integers) c c ***************************************************************** implicit real(a-h,o-z),integer(i-n) c include 'comm_para.h' dimension isx(1),isy(1),iox(1),ioy(1),h(1) dimension ixk(npt,nz),iyk(npt,nz),isxk(npt,nz),isyk(npt,nz), + lxxk(MXBDY,nz),lyyk(MXBDY,nz),lxyk(MAXNB,nz),lyxk(MAXNB,nz), + snxk(MAXNB,nz),snyk(MAXNB,nz), mask(nxp*nyp,nz), dept(npt), + ifxk(9*MAXSID,nz),ifpxk(5*MAXSID,nz),ifyk(9*MAXSID,nz), + lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz),lok(4*MAXSID,nz),nzi(npt) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) common /new_filt/ MAXFO, nxk, nyk, nfx, nfpx, nfy c if (MINSEG.le.0 .or. nxp.le.0 .or. nyp.le.0) return nptk(1) = npt c.....iox and mask have already been created for top level, here we c.....initialize mask for rest of levels call maskk (npt,nxp,nyp,nz,nzi,iox,mask) c.....get the nbx x-bndry and the nby y-bndry indices call bndxy (npt,iox,ioy,lxxk,lyyk,lxyk,lyxk,snxk,snyk, * isxk,isyk,lpbcwk,lpbcek,mask,nzi,h,ixk,iyk,lok,dept) c.....compute indices for shapiro filters do k = 1, nz call shap_indx (nptk(k),nxp,nyp,mask(1,k),isxk(1,k), * nfxk(k),nfpxk(k),nfyk(k),ifxk(1,k),ifpxk(1,k),ifyk(1,k)) enddo c.....find the ncs interior-corner indices for the compressed x-sort. c.....store them at the end of lxx, and store the signs in snx. do k = 1, nz ncx = nbxk(k) ncy = nbyk(k) call newcorn (nxp,nyp,ncsk(k),lxxk(ncx+1,k),snxk(ncx+1,k), * lyyk(ncy+1,k),snyk(ncy+1,k),mask(1,k),isxk(1,k)) enddo need = MINSEG*(max0(nbx,nby)+ncs) if (need .gt. MXBDY) call wspace('MXBDY', need) c.....store copies of the corner indices for both sorts in lyx and lxy. do k = 1, nz ncx = nbxk(k) ncy = nbyk(k) do i = 1, ncsk(k) lyxk(ncy+i,k) = lxxk(ncx+i,k) lxyk(ncx+i,k) = lyyk(ncy+i,k) enddo enddo c do k = 1, nz nc = ncsk(k) ncx = nbxk(k) ncy = nbyk(k) nbxc = ncx + nc nbyc = ncy + nc do m = 1, MINSEG - 1 do i = 1, ncx lxxk(m*nbxc+i,k) = lxxk(i,k) + isign(m, int(snxk(i,k))) enddo do i = 1, ncy lyyk(m*nbyc+i,k) = lyyk(i,k) + isign(m, int(snyk(i,k))) enddo enddo do i = 1, nc lxxk(nbxc+i+ncx,k) = lxxk(i+ncx,k) + isign(1, int(snxk(i+ncx,k))) lyyk(nbyc+i+ncy,k) = lyyk(i+ncy,k) + isign(1, int(snyk(i+ncy,k))) enddo enddo c now combine a few of the mappings: do k = 1, nz npk = nptk(k) do i = 1, npk if (isyk(i,k).lt.1.or.isyk(i,k).gt.npt) print*,'trouble',i,k isyk(i,k) = ixk(isyk(i,k),k) enddo nbyck = nbyk(k) + ncsk(k) nbxck = nbxk(k) + ncsk(k) do i = 1, MINSEG*nbyck ly = lyyk(i,k) if (ly.le.npk.and.ly.ge.1) lyyk(i,k) = isyk(ly,k) enddo do i = 1, nbxck ly = lxyk(i,k) if (ly.le.npk.and.ly.ge.1) lxyk(i,k) = isyk(ly,k) enddo do i = 1, MINSEG*nbxck ly = lxxk(i,k) if (ly.le.npk.and.ly.ge.1) lxxk(i,k) = ixk(ly,k) enddo do i = 1, nbyck ly = lyxk(i,k) if (ly.le.npk.and.ly.ge.1) lyxk(i,k) = ixk(ly,k) enddo do i = 1, npbck(k) ly = lpbcek(i,k) if (ly.le.npk.and.ly.ge.1) lpbcek(i,k) = ixk(ly,k) ly = lpbcwk(i,k) if (ly.le.npk.and.ly.ge.1) lpbcwk(i,k) = ixk(ly,k) enddo do i = 1, nlok(k) ly = lok(i,k) if (ly.le.npk.and.ly.ge.1) lok(i,k) = ixk(ly,k) enddo enddo return c end of bndrys. end c------------------------------------------------------------------ subroutine bndxy(npt,iox,ioy,lxx,lyy,lxy,lyx,snx,sny, + isxk,isyk, lpbcw, lpbce, mask, nzi,h,ixk, iyk, lok, dept) c------------------------------------------------------------------ c Get the compressed x-sort boundary indices. c iox = (input) nxyc indices of the ocean points for the x-sort. c maxnb = (input) max. storage space for lxx, lyy, lxy, or lyx. c = default maximum allowable value for max0(nbx,nby). c minseg = (input) required minimum # of consecutive ocean points c interior to and including each boundary point. c lxx = (output) x-bndry indices for the compressed x-sort. c lyy = (output) y-bndry indices for the compressed y-sort. c lxy = (output) x-bndry indices for the compressed y-sort. c lyx = (output) y-bndry indices for the compressed x-sort. c snx = (output) nbx signs for the x-boundaries. c sny = (output) nby signs for the y-boundaries. c ioy = (output) nxyc indices of the ocean points for the y-sort. c isx = (output) nxyc indices to gather the compressed x-sort c from the compressed y-sort. c isy = (output) nxyc indices to gather the compressed y-sort c from the compressed x-sort. c implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' include 'comm_para.h' dimension iox(1),ioy(1),isxk(npt,nz),isyk(npt,nz), + lxx(MXBDY,nz),lyy(MXBDY,nz),lxy(MAXNB,nz),lyx(MAXNB,nz), + snx(MAXNB,nz),sny(MAXNB,nz), mask(nxp*nyp,nz), nzi(npt), + h(1),lpbcw(MAXSID,nz),lpbce(MAXSID,nz),lok(4*MAXSID,nz), + ixk(npt,nz),iyk(npt,nz), dept(npt) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) c.....convert the x-sort indices in iox to y-sort indices ioy. nxy1 = 1-nxp*nyp c c.....convert cumulation location of ocean points from rowwise to c.....colomnwise. do i=1,nxyc ioy(i) = ((iox(i)-1)/nxp)*nxy1 + (iox(i)-1)*nyp + 1 enddo c sort by increasing value so that ioy(1) is the bottom of the c first ocean column and not the beginning of the first ocean c row. the sort order is isy. call sorti(nxyc,ioy,1,isyk) c form isxk from isyk, i.e. find the rowwise ordering of the c compressed y-sort indices. do i=1,nxyc isxk(isyk(i,1),1) = i enddo c.....find x-bndry indices for the compressed x-sort. do i=1,nxyc ixk(i,1) = i iyk(i,1) = i enddo if (iglob .eq. 0) then call bound(nxyc,iox,ixk,nxp,maxnb,minseg,nbxk,lxx,snx) else call set_pbc (nxp,nyp, npbck, lpbcw, lpbce, mask) call set_bpx (nxp,nyp,mask,maxnb,minseg,nbxk,lxx,snx) endif call reset_mask (nxyc,nxp,nyp,nz,nzi,h,mask,MINSEG,nptk, * ixk,iyk,isxk,isyk,dept) do k = 2, nz do i=1,nptk(k) isxk(isyk(i,k),k) = i enddo enddo call make_lok(npt,nxp,nyp,nz,iox,mask,nlok,lok) do k = 2, nz npbck(k) = 0 if (iglob .eq. 0) then call bound(nptk(k),iox,ixk(1,k),nxp,maxnb,minseg, * nbxk(k),lxx(1,k),snx(1,k)) else call set_pbck (nxp,nyp,npbck(k),lpbcw(1,k),lpbce(1,k), * mask(1,k)) call set_bpxk (nxp,nyp,mask(1,k),nbxk(k),lxx(1,k),snx(1,k)) endif enddo do k = 1, nz c.....find y-bndry indices for the compressed y-sort. call bound(nptk(k),ioy,iyk(1,k),nyp,maxnb,minseg, * nbyk(k),lyy(1,k),sny(1,k)) c find the y-bndry indices for the compressed x-sort using the c columnwise ordering of the compressed x-sort. do i=1,nbyk(k) lyx(i,k) = isyk(lyy(i,k),k) enddo c sort in sequential order. call sorti(nbyk(k),lyx(1,k),0) enddo c find the x-bndry indices for the compressed y-sort using the c rowwise ordering of the compressed x-sort. do k = 1, nz do i=1,nbxk(k) lxy(i,k) = isxk(lxx(i,k),k) enddo call sorti(nbxk(k),lxy(1,k),0) enddo return c end of bndxy. end c ------------------------------------------------------------------ subroutine bound(npk,ioc,ixk,nsid,maxnb,minseg,nb,lbn,sgn) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) character*72 msg dimension ioc(1),lbn(1),sgn(1),ixk(1) c if (npk.eq.0) then nb = 0 return endif nb = 1 lbn(1) = 1 sgn(1) = 1. c c cycle through all ocean points and note the location of all c nonconsecutive indices and points at the absolute extreme of c the rectangular grid. c do 10 i=2,npk ii = ixk(i) im = ixk(i-1) if(ioc(ii)-ioc(im).gt.1 .or. mod(ioc(ii)-1,nsid).eq.0) then if(nb+3 .gt. maxnb) goto 20 c c store location and set sign for the previous end of row c or top of column. c nb = nb + 1 lbn(nb) = i - 1 sgn(nb) = -1. c c make sure the boundary geometry will permit a 4th order c differencing. c if(lbn(nb)-lbn(nb-1)+1 .lt. minseg) goto 30 c c store location and set sign for the start of row or bottom c of column. c nb = nb + 1 lbn(nb) = i sgn(nb) = 1. endif 10 continue c c register last point. c nb = nb + 1 lbn(nb) = npk sgn(nb) = -1. if(lbn(nb)-lbn(nb-1)+1 .lt. minseg) goto 30 return 20 write(msg,21) maxnb 21 format('bound: insufficient space for lb. maxnb=',i10,'$') call perror1(msg,1) 30 j = (ioc(lbn(nb))-1)/nsid + 1 i = ioc(lbn(nb)) - (j-1)*nsid write(msg,31) lbn(nb)-lbn(nb-1)+1,i,j 31 format('bound: only',i3,' consecutive ocean grid pts', + ' near (i,j or j,i)= ',2i5,'$') call perror1(msg,1) c end of bound. end c ------------------------------------------------------------------ subroutine gridxy(nxp,nyp,x1,x2,y1,y2,nsx,nsy,nystrch,xs,alpha,beta, * x,y,xp,yp,xpp,ypp) c ------------------------------------------------------------------ c compute the x and y grid point coordinates. c c nxp,nyp = (input) # of grid points in the x and y directions. c x1,x2 = (input) minimum and maximum x-coordinate. c y1,y2 = (input) minimum and maximum y-coordinate. c nsx,nsy = (input) # of atan's composing the stretched grid c transformation function for the x and y directions. c xs = (input) locations of the atan's for x: xs(1 to nsx), c and for y: xs(nsx+1 to nsx+nsy). c alpha = (input) parameters that will determine the # of grid c points in a stretched region; c x: alpha(1 to nsx); and for y: alpha(nsx+1 to nsx+nsy). c beta = (input) parameters that determines the scale width of c each stretched region; c x: beta(1 to nsx); and for y: beta(nsx+1 to nsx+nsy). c x,y = (output) nxp x and nyp y grid coordinates. c xp = (output) nxp derivatives of the x-transformation c function: d(psi1(x))/d(x). c yp = (output) nyp derivatives of the y-transformation c function: d(psi2(y))/d(y). c dimension xs(1),alpha(1),beta(1),x(1),y(1),xp(1),yp(1),xpp(1),ypp(1) c if(nsx.le.0) then c no stretching of the x-coordinates. delx = (x2-x1)/float(nxp-1) do i = 1, nxp x(i) = x1 + real(i-1)*delx xp(i) = 1. xpp(i) = 0. enddo else c coordinate stretching in the x-direction. call stretch(nxp,x1,x2,nsx,xs,alpha,beta,x,xp,xpp) endif if (nystrch.eq.1) then if(nsy.le.0) then c no stretching of the y-coordinates. dely = (y2-y1)/float(nyp-1) do i = 1, nyp y(i) = y1 + (i-1)*dely yp(i) = 1. ypp(i) = 0. enddo else c coordinate stretching in the y-direction. i = nsx + 1 call stretch(nyp,y1,y2,nsy,xs(i),alpha(i),beta(i),y,yp,ypp) endif elseif (nystrch.eq.2) then dely = (y2-y1)/float(nyp-1) b = (y2-y1)/(sind(y2)-sind(y1)) a = y1 - b* sind(y1) to_pie = atan(1.)/45. do i = 1, nyp yl = y1 + (i-1)*dely y(i) = a + b * sind(yl) yp(i) = to_pie* b * cosd(yl) ypp(i)= - to_pie*to_pie* b * sind(yl) enddo endif return c end of gridxy. end c ------------------------------------------------------------------ subroutine sorti(n,ix,key,ist) c ------------------------------------------------------------------ c sort ix by increasing values and return the sort order in ist. c c n = (input) length of ix. c ix = (input) array to be sorted. c = (output) array sorted by increasing values. c key = (input) flag: c = 1; return the sort order in array ist. c = otherwise; do not use array ist. (call sorti(n,ix,0) is ok.) c ist = (output) array containing the sort order if key=1; c ie. ix output(i) = ix input(ist(i)), i=1,n c implicit real(a-h,o-z),integer(i-n) dimension ix(0:n-1),ist(0:n-1) c c see "the c programming language", kernighan and ritchie, page 58 c for this algorithm. c igap = n/2 if(key.eq.0) then 10 if(igap .le. 0) return do 30 i=igap,n-1 do 20 j=i-igap,0,-igap if(ix(j) .lt. ix(j+igap)) goto 30 itemp = ix(j) ix(j) = ix(j+igap) ix(j+igap) = itemp 20 continue 30 continue igap = igap/2 goto 10 else c do 40 i=0,n-1 40 ist(i) = i+1 50 if(igap .le. 0) return do 70 i=igap,n-1 do 60 j=i-igap,0,-igap if(ix(j) .lt. ix(j+igap)) goto 70 itemp = ix(j) ix(j) = ix(j+igap) ix(j+igap) = itemp itemp = ist(j) ist(j) = ist(j+igap) ist(j+igap) = itemp 60 continue 70 continue igap = igap/2 goto 50 endif c end of sorti. end subroutine newcorn (nxp,nyp,nc,lxx,snx,lyy,sny,mask,isx) c------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) dimension in(4), lxx(1), lyy(1), mask(1), snx(1), sny(1), isx(1) do j = 1, nyp-1 do i = 1, nxp-1 k = nxp*(j-1) + i in(1) = k in(2) = k + 1 in(3) = k + nxp + 1 in(4) = k + nxp land = 0 do m = 1, 4 if (mask(in(m)) .eq. 0) then land = land + 1 m0 = m endif enddo if (land .eq. 1) then nc = nc + 1 if (m0 .le. 2) then n0 = m0 + 2 sy = 1. sx = 3-2*m0 else n0 = m0 - 2 sy = -1. sx = 2*m0-7 endif it = in(n0) ma = mask(it) lxx(nc) = ma snx(nc) = sx c lyy(nc) = (1-nxp*nyp)*((it-1)/nxp) + (it-1)*nyp + 1 lyy(nc) = isx(ma) sny(nc) = sy endif enddo enddo return end subroutine make_iox (nx, ny, mask, iox, nlo, lo, nsponge, lsponge, * nrelax, lrelax, iglob, inx, iny) c------------------------------------------------------------- dimension mask(1), iox(1), lo(1), lsponge(1), lrelax(1), inx(1), iny(1) nlo = 0 nsponge = 0 nrelax = 0 k = 0 do iyy = 1, ny do ixx = 1, nx i = (iyy-1)*nx + ixx ma = mask(i) if (ma .ne. 0) then k = k + 1 iox(k) = i inx(k) = ixx iny(k) = iyy if (ma .eq. 2) then nlo = nlo + 1 lo(nlo) = k elseif (ma .eq. 3) then nsponge = nsponge + 1 lsponge(nsponge) = i elseif (ma .eq. 4) then nrelax = nrelax + 1 lrelax(nrelax) = i elseif (ma .eq. 5) then nsponge = nsponge + 1 lsponge(nsponge) = i nrelax = nrelax + 1 lrelax(nrelax) = i endif mask(i) = k endif enddo enddo if (iglob .eq. 1) then do j = 1, ny j1 = 1 + (j-1)*nx jnx = j*nx if ( (mask(j1) .eq. 0 .and. mask(jnx) .ne. 0).or. * (mask(jnx) .eq. 0 .and. mask(j1) .ne. 0)) then print*,'ocean/land periodic boundary not allowed' stop endif enddo endif return end c ------------------------------------------------------------------ subroutine maskk (npt,nxp,nyp,nz,nzi,iox,mask) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension mask(nxp*nyp,nz), nzi(npt), iox(npt) do k = 2, nz do i = 1, nxp*nyp mask(i,k) = 0 enddo enddo do i = 1, npt do k = 2, nzi(i) ixy = iox(i) mask(ixy,k) = mask(ixy,1) enddo enddo return end c ------------------------------------------------------------------ subroutine reset_mask(npt,nxp,nyp,nz,nzi,h,mask,MINSEG,nptk, * ixk,iyk,isxk,isyk,dept) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_new.h' dimension mask(nxp,nyp,nz), h(npt,nz), nptk(nz), nzi(npt), * ixk(npt,nz),iyk(npt,nz),isxk(npt,nz),isyk(npt,nz),dept(npt) logical prev, curr c do k = nz, 2, -1 10 continue iflag = 0 do irow = 1, nyp prev = (mask(1, irow, k) .eq. 0) ista = 1 do icol = 2, nxp ixy = mask(icol, irow, k) curr = (ixy .eq. 0) if ( curr .ne. prev ) then if ( prev ) then ista = icol else if (icol-ista .lt. MINSEG) then do i = ista, icol-1 ixym = mask(i,irow,1) mask(i,irow,k) = 0 if (initbt.eq.2)h(ixym,k-1)=h(ixym,k-1)+h(ixym,k) c h(ixym, k) = -98765432. h(ixym, k) = 0. nzi(ixym) = nzi(ixym) - 1 iflag = 1 enddo endif endif prev = curr endif enddo c check last segment if (.not.curr.and.nxp-ista+1.lt.MINSEG) then do i = ista, nxp ixym = mask(i,irow,1) mask(i,irow,k) = 0 if (initbt.eq.2)h(ixym,k-1)=h(ixym,k-1)+h(ixym,k) c h(ixym, k) = -98765432. h(ixym, k) = 0. nzi(ixym) = nzi(ixym) - 1 iflag = 1 enddo endif enddo c now make sure periodic boundaries and land boundaries do not coincide if (iglob .eq. 1) then do irow = 1, nyp ixy = mask(1,irow,1) ixyp = mask(nxp,irow,1) ixyk = mask(1,irow,k) ixypk = mask(nxp,irow,k) if (ixyk.ne.0.and.ixypk.eq.0) then c--------------- turn (1,irow,k) into land point mask(1,irow,k) = 0 h(ixy,k) = 0. nzi(ixy) = nzi(ixy) - 1 iflag = 1 endif if (ixypk.ne.0.and.ixyk.eq.0) then c--------------- turn (nxp,irow,k) into land point mask(nxp,irow,k) = 0 h(ixyp,k) = 0. nzi(ixyp) = nzi(ixyp) - 1 iflag = 1 endif enddo endif c now same in y-direction do icol = 1, nxp prev = (mask(icol, 1, k) .eq. 0) ista = 1 do irow = 2, nyp ixy = mask(icol, irow, k) curr = (ixy .eq. 0) if ( curr .ne. prev ) then if ( prev ) then ista = irow else if (irow-ista .lt. MINSEG) then do i = ista, irow-1 ixym = mask(icol,i,1) mask(icol,i,k) = 0 if (initbt.eq.2)h(ixym,k-1)=h(ixym,k-1)+h(ixym,k) c h(ixym, k) = -98765432. h(ixym, k) = 0. nzi(ixym) = nzi(ixym) - 1 iflag = 1 enddo endif endif prev = curr endif enddo c check last segment if (.not.curr.and.nyp-ista+1.lt.MINSEG) then do i = ista, nyp ixym = mask(icol,i,1) mask(icol,i,k) = 0 if (initbt.eq.2)h(ixym,k-1)=h(ixym,k-1)+h(ixym,k) c h(ixym, k) = -98765432. h(ixym, k) = 0. nzi(ixym) = nzi(ixym) - 1 iflag = 1 enddo endif enddo c repeat until iflag = 0 if (iflag.ne.0) goto 10 c now count the number of points on the k-th grid: npk = 0 do icol = 1, nxp do irow = 1, nyp ma = mask(icol,irow,k) if (ma.ne.0) then npk = npk + 1 iyk(npk,k) = isxk(ma,1) endif enddo enddo npk = 0 do irow = 1, nyp do icol = 1, nxp ma = mask(icol,irow,k) if (ma.ne.0) then npk = npk + 1 ixk(npk,k) = ma mask(icol,irow,k) = npk endif enddo enddo nptk(k) = npk npk = 0 do icol = 1, nxp do irow = 1, nyp ma = mask(icol,irow,k) if (ma.ne.0) then npk = npk + 1 isyk(npk,k) = ma endif enddo enddo enddo nzt = nzi(1) nzi(1) = 0 write(80,*)'number of vertical layers: (n.b., lower left value should be ',nzt write(80,*)'to make a postscript version of this page,' write(80,*)'use: a2ps -nP -1 -F4 fort.80 > file.ps' write(80,*) do j = nyp,1,-1 write(80,101)(nzi(max(1,mask(i,j,1))),i=1,nxp) enddo write(80,*) write(80,*) write(80,*) do i = 1,nxp write(80,101)(nzi(max(1,mask(i,j,1))),j=1,nyp) enddo nzi(1) = nzt close(80) do i = 1, npt dept(i) = h(i,1) if (nzi(i).eq.1) then print*,'only one grid level at grid point',i print*,'aborting your run. try changing Bath_min, or Z_profile' stop endif do k = 2, nzi(i) dept(i) = dept(i) + h(i,k) enddo enddo do j = nyp,1,-1 write(91,102)(mod(mask(i,j,1),100),i=1,nxp) enddo close(91) 101 format(300i2) 102 format(300i2) return end subroutine set_pbck (nxp, nyp, npbc, lpbcw, lpbce, mask) c---------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension lpbcw(1), lpbce(1), mask(nxp*nyp) do j = 1, nyp j1 = 1 + (j-1)*nxp jnx = j*nxp ierr = 0 if ( mask(j1) .ne. 0 ) then if ( mask(jnx) .eq. 0 ) then ierr = nxp+1 else do i = 1, MINSEG - 1 if (ierr.eq.0 .and. mask(j1+i).eq.0) then ierr = i elseif (ierr.eq.0 .and. mask(jnx-i).eq.0) then ierr = nxp+1-i endif enddo endif if (ierr .ne. 0) then write (6, *) 'A land is too close for PBC [i,j]=',ierr,j else npbc = npbc + 1 lpbcw(npbc) = mask(j1) lpbce(npbc) = mask(jnx) endif endif enddo return end c ------------------------------------------------------------------ subroutine set_bpxk (nxp, nyp, mask, nbx,lxx,snx) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) dimension mask(nxp,1), lxx(1),snx(1) logical prev, curr c nbx = 0 do irow = 1, nyp prev = (mask(1, irow) .eq. 0) ista = 1 do icol = 2, nxp ixy = mask(icol, irow) curr = (ixy .eq. 0) if ( curr .ne. prev ) then nbx = nbx + 1 if ( prev ) then lxx(nbx) = ixy snx(nbx) = 1. ista = icol else lxx(nbx) = mask(icol-1, irow) snx(nbx) = -1. endif prev = curr endif enddo enddo end c ------------------------------------------------------------------ subroutine make_lok(npt,nxp,nyp,nz,iox,mask,nlok,lok) c ------------------------------------------------------------------ implicit real(a-h,o-z),integer(i-n) include 'comm_para.h' dimension mask(nxp*nyp,nz), iox(npt), lok(4*MAXSID,nz), nlok(nz) do k = 2, nz nl = 0 do i= 1, nlok(1) ma = mask( iox( lok(i,1)), k) if (ma.ne.0) then nl = nl+1 lok(nl,k) = ma endif enddo enddo return end c------------------------------------------------------------------------ subroutine comp_rotma (al, be, ga) c------------------------------------------------------------------------ c Computes the tranformation matrices A() and A'() which define the c transformations between Geographical [G] and Rotated [R] coordinate c systems: c [XG] = [A] * [XR]; and: [XR] = [A'] * [XG] c (reference: A.Korn,M.Korn "Mathematical Handbook", 1968) c Senya Basin, 1996 c------------------------------------------------------------------------- common /pole_rotm/ a(3,3), ap(3,3) call rmat32 ( al, be, ga, a ) call rmat32 (-ga, -be, -al, ap ) return end c------------------------------------------------------------------------ subroutine rot_g2r (n, x, y) c------------------------------------------------------------------------ c Converts a vector of latitude longitude points [X(N),Y(N)] c from Geo to Rotated System. c Senya Basin, 1996 c------------------------------------------------------------------------- common /pole_rotm/ a(3,3), ap(3,3) call rot_coor180(n, x, y, ap) return end c------------------------------------------------------------------------ subroutine rot_r2g (n, x, y) c------------------------------------------------------------------------ c Converts a vector of latitude longitude points [X(N),Y(N)] c from Rotated to Geo System. c Senya Basin, 1996 c------------------------------------------------------------------------- common /pole_rotm/ a(3,3), ap(3,3) call rot_coor180(n, x, y, a) return end subroutine rmat32 (al, be, ga, a) c------------------------------------------------------------------------ c Computes the A32() transformation matrix:: c A32() = R(Z,alpha) * R(Y,beta) * R(Z,gamma); c (reference: A.Korn,M.Korn "Mathematical Handbook", 1968) c Senya Basin, 1996 c------------------------------------------------------------------------- dimension a(3,3) sin_a = sind(al) cos_a = cosd(al) sin_b = sind(be) cos_b = cosd(be) sin_g = sind(ga) cos_g = cosd(ga) co_ab = cos_a * cos_b sa_cb = sin_a * cos_b c.....transformation from Geo to a Rotated system: a(1,1)= co_ab * cos_g - sin_a * sin_g a(1,2)= - co_ab * sin_g - sin_a * cos_g a(1,3)= cos_a * sin_b a(2,1)= sa_cb * cos_g + cos_a * sin_g a(2,2)= - sa_cb * sin_g + cos_a * cos_g a(2,3)= sin_a * sin_b a(3,1)= - sin_b * cos_g a(3,2)= sin_b * sin_g a(3,3)= cos_b return end subroutine rot_coor(n, x, y, a) c------------------------------------------------------------------------ c Purpose: To rotate a given sequence of points according to matrix A() c x(output) : [0:360] Senya Basin, 1996 c------------------------------------------------------------------------- dimension x(1), y(1), a(3,3) parameter (PI = 3.14159265358979323846) parameter (D2R = PI/180., R2D = 180./PI) do i = 1, n z1 = SIN(D2R * y(i)) cos_t1 = COS(D2R * y(i)) y1 = cos_t1 * SIN(D2R * x(i)) x1 = cos_t1 * COS(D2R * x(i)) x2 = a(1,1)*x1 + a(1,2)*y1 + a(1,3)*z1 y2 = a(2,1)*x1 + a(2,2)*y1 + a(2,3)*z1 z2 = a(3,1)*x1 + a(3,2)*y1 + a(3,3)*z1 y(i) = R2D * ASIN(z2) x(i) = R2D * ATAN2(y2, x2) + 180. * (1.- SIGN(1., y2)) enddo return end subroutine rot_coor180(n, x, y, a) c------------------------------------------------------------------------ c Purpose: To rotate a given sequence of points according to matrix A() c x(output) : [-180:180] Senya Basin, 1996 c------------------------------------------------------------------------- dimension x(1), y(1), a(3,3) parameter (PI = 3.14159265358979323846) parameter (D2R = PI/180., R2D = 180./PI) do i = 1, n z1 = SIN(D2R * y(i)) cos_t1 = COS(D2R * y(i)) y1 = cos_t1 * SIN(D2R * x(i)) x1 = cos_t1 * COS(D2R * x(i)) write(93,*)x(i),y(i) x2 = a(1,1)*x1 + a(1,2)*y1 + a(1,3)*z1 y2 = a(2,1)*x1 + a(2,2)*y1 + a(2,3)*z1 z2 = a(3,1)*x1 + a(3,2)*y1 + a(3,3)*z1 write(94,120)x1,y1,z1,x2,y2,z2 y(i) = R2D * ASIN(z2) x(i) = R2D * ATAN2(y2, x2) write(95,*)x(i),y(i) enddo close(93) close(94) close(95) 120 format(6f13.9) return end function rot_fcr2g(x, y) c------------------------------------------------------------------------ c A function for a computation of the Coriolis Term in a c a System with a Rotated Pole: c (input): in degrees, (output):sin(theta) Senya Basin, 1996 c------------------------------------------------------------------------- parameter (PI = 3.14159265358979323846) parameter (D2R = PI/180., R2D = 180./PI) common /pole_rotm/ a(3,3), ap(3,3) z1 = SIN(D2R * y) cos_t1 = COS(D2R * y) y1 = cos_t1 * SIN(D2R * x) x1 = cos_t1 * COS(D2R * x) rot_fcr2g = a(3,1)*x1 + a(3,2)*y1 + a(3,3)*z1 return end subroutine rot_point180(xin, yin, x, y) c------------------------------------------------------------------------ c Purpose: To rotate a point according to matrix A() c x : [-180:180] Naomi, 1997 c------------------------------------------------------------------------- parameter (PI = 3.14159265358979323846) parameter (D2R = PI/180., R2D = 180./PI) common /pole_rotm/ a(3,3), ap(3,3) z1 = SIN(D2R * yin) cos_t1 = COS(D2R * yin) y1 = cos_t1 * SIN(D2R * xin) x1 = cos_t1 * COS(D2R * xin) x2 = a(1,1)*x1 + a(1,2)*y1 + a(1,3)*z1 y2 = a(2,1)*x1 + a(2,2)*y1 + a(2,3)*z1 z2 = a(3,1)*x1 + a(3,2)*y1 + a(3,3)*z1 y = R2D * ASIN(z2) x = R2D * ATAN2(y2, x2) return end c------------------------------------------------------------------------ subroutine rot_vec (nx,ny,n, iox, x, y, f, g) c------------------------------------------------------------------------ c Converts an unrotated vector [f(n),g(n)] at c rotated lat,long points [X(nx),Y(ny)] c to a rotated vector at the rotated lat,long points c Naomi, 1998 c------------------------------------------------------------------------- parameter (PI = 3.14159265358979323846) parameter (D2R = PI/180., R2D = 180./PI) dimension x(nx), y(ny), f(n), g(n), iox(n) common /pole_rotm/ a(3,3), ap(3,3) do k = 1, n j = (iox(k)-1)/nx + 1 i = iox(k) - (j-1)*nx x1 = x(i) y1 = y(j) shift = 0. if (x1.gt.180) shift = 180. c if (x1.lt.-180) shift = -180. zz = SIN(D2R * y1) cos_t1 = COS(D2R * y1) yy = cos_t1 * SIN(D2R * x1) xx = cos_t1 * COS(D2R * x1) x0 = a(1,1)*xx + a(1,2)*yy + a(1,3)*zz y0 = a(2,1)*xx + a(2,2)*yy + a(2,3)*zz z0 = a(3,1)*xx + a(3,2)*yy + a(3,3)*zz y1 = R2D * ASIN(z0) x1 = R2D * ATAN2(y0, x0) + shift * (1.- SIGN(1., y0)) z1 = sqrt(f(k)**2 + g(k)**2) if (z1.ne.0) then x1 = x1 + f(k)/z1 y1 = y1 + g(k)/z1 endif zz = SIN(D2R * y1) cos_t1 = COS(D2R * y1) yy = cos_t1 * SIN(D2R * x1) xx = cos_t1 * COS(D2R * x1) x0 = a(1,1)*xx + a(2,1)*yy + a(3,1)*zz y0 = a(1,2)*xx + a(2,2)*yy + a(3,2)*zz z0 = a(1,3)*xx + a(2,3)*yy + a(3,3)*zz y1 = R2D * ASIN(z0) x1 = R2D * ATAN2(y0, x0) + shift * (1.- SIGN(1., y0)) f(k) = x1 - x(i) g(k) = y1 - y(j) if (f(k).gt.180) f(k) = f(k) + 360. if (f(k).lt.-180) f(k) = f(k) + 360. s = sqrt(f(k)**2 + g(k)**2) if (s.ne.0) then f(k) = z1*f(k)/s g(k) = z1*g(k)/s endif enddo return end