subroutine bigrid(depth) 1,8 c c --- set loop bounds for irregular basin in c-grid configuration c --- q,u,v,p are vorticity, u-velocity, v-velocity, and mass points, resp. c --- 'depth' = basin depth array, zero values indicate land c c --- this version works for both cyclic and noncyclic domains. c --- land barrier at i=ii and/or j=jj signals closed-basin conditions c implicit none include 'dimensions.h' include 'dimension2.h' c real depth(idm,jdm) integer nfill,nzero,jsec,jfrst,jlast character fmt*12,char2*2 data fmt/'(i4,1x,75i1)'/ c c$OMP PARALLEL DO do 17 j=1,jj do 17 i=1,ii ip(i,j)=0 iq(i,j)=0 iu(i,j)=0 17 iv(i,j)=0 c$OMP END PARALLEL DO c go to 888 ! no change in topo after weights are done c --- fill single-width inlets 16 nfill=0 c$OMP PARALLEL DO PRIVATE(ja,jb,ia,ib,nzero) REDUCTION(+:nfill) do 15 j=1,jj ja=mod(j-2+jj,jj)+1 jb=mod(j ,jj)+1 do 15 i=1,ii ia=mod(i-2+ii,ii)+1 ib=mod(i ,ii)+1 nzero=0 if (depth(i,j).gt.0.) then if (depth(ia,j).le.0.) nzero=nzero+1 if (depth(ib,j).le.0.) nzero=nzero+1 if (depth(i,ja).le.0.) nzero=nzero+1 if (depth(i,jb).le.0.) nzero=nzero+1 if (nzero.ge.3) then write (lp,'(a,i4,a,i4,a)') ' depth(',i,',',j,') set to zero' stop 'pre-process depth' depth(i,j)=0. nfill=nfill+1 end if end if 15 continue c$OMP END PARALLEL DO if (nfill.gt.0) go to 16 c 888 continue c --- mass points are defined where water depth is greater than zero c$OMP PARALLEL DO do 2 j=1,jj do 2 i=1,ii if (depth(i,j).gt.0.) ip(i,j)=1 2 continue c$OMP END PARALLEL DO c c --- u,v points are located halfway between any 2 adjoining mass points c$OMP PARALLEL DO PRIVATE(ia) do 3 j=1,jj do 3 i=1,ii ia=mod(i-2+ii,ii)+1 if (ip(ia,j).gt.0.and.ip(i,j).gt.0) iu(i,j)=1 3 continue c$OMP END PARALLEL DO c$OMP PARALLEL DO PRIVATE(ja) do 4 j=1,jj ja=mod(j-2+jj,jj)+1 do 4 i=1,ii if (ip(i,ja).gt.0.and.ip(i,j).gt.0) iv(i,j)=1 4 continue c$OMP END PARALLEL DO c c --- 'interior' q points require water on all 4 sides. c$OMP PARALLEL DO PRIVATE(ja,ia) do 5 j=1,jj ja=mod(j-2+jj,jj)+1 do 5 i=1,ii ia=mod(i-2+ii,ii)+1 if (min(ip(i,j),ip(ia,j),ip(i,ja),ip(ia,ja)).gt.0) iq(i,j)=1 5 continue c$OMP END PARALLEL DO c c --- 'promontory' q points require water on 3 (or at least 2 diametrically c --- opposed) sides c$OMP PARALLEL DO PRIVATE(ja,ia) do 10 j=1,jj ja=mod(j-2+jj,jj)+1 do 10 i=1,ii ia=mod(i-2+ii,ii)+1 if ((ip(i ,j).gt.0.and.ip(ia,ja).gt.0).or. . (ip(ia,j).gt.0.and.ip(i ,ja).gt.0)) iq(i,j)=1 10 continue c$OMP END PARALLEL DO c c --- determine loop bounds for vorticity points, including interior and c --- promontory points call indxi(iq,ifq,ilq,isq) call indxj(iq,jfq,jlq,jsq) c c --- determine loop indices for mass and velocity points call indxi(ip,ifp,ilp,isp) call indxj(ip,jfp,jlp,jsp) call indxi(iu,ifu,ilu,isu) call indxj(iu,jfu,jlu,jsu) call indxi(iv,ifv,ilv,isv) call indxj(iv,jfv,jlv,jsv) c c --- write out -ip- array c --- data are written in strips 75 points wide jsec=(jj-1)/75 do 9 jfrst=0,75*jsec,75 jlast=min(jj,jfrst+75) write (char2,'(i2)') jlast-jfrst fmt(8:9)=char2 write (lp,'(''ip array, cols'',i5,'' --'',i5)') jfrst+1,jlast 9 write (lp,fmt) (i,(10*ip(i,j),j=jfrst+1,jlast),i=1,ii) c return end c c subroutine indxi(ipt,if,il,is) 4 c c --- input array ipt contains 1 at grid point locations, 0 elsewhere c --- output is arrays if, il, is where c --- if(j,k) gives row index of first point in column j for k-th section c --- il(j,k) gives row index of last point c --- is(j) gives number of sections in column j (maximum: ms) c implicit none include 'dimensions.h' include 'dimension2.h' c integer ipt(idm,jdm),if(jdm,ms),il(jdm,ms),is(jdm) do 1 j=1,jj is(j)=0 do 4 k=1,ms if(j,k)=0 4 il(j,k)=0 i=1 k=1 3 if (ipt(i,j).ne.0) go to 2 i=i+1 if (i.le.ii) go to 3 go to 1 2 if (k.gt.ms) then write (lp,'('' error in indxi - ms too small at i,j ='',2i5)') i,j write (lp,'('' j-th line of ipt array:'',/(7(1x,10i1)))') . (ipt(l,j),l=1,ii) stop '(indxi)' end if if(j,k)=i 6 i=i+1 if (i.le.ii) go to 5 il(j,k)=ii is(j)=k go to 1 5 if (ipt(i,j).ne.0) go to 6 il(j,k)=i-1 is(j)=k k=k+1 go to 3 1 continue return end c c subroutine indxj(jpt,jf,jl,js) 4 c c --- input array jpt contains 1 at grid point locations, 0 elsewhere c --- output is arrays jf, jl, js where c --- jf(i,k) gives column index of first point in row i for k-th section c --- jl(i,k) gives column index of last point c --- js(i) gives number of sections in row i (maximum: ms) c implicit none include 'dimensions.h' include 'dimension2.h' c integer jpt(idm,jdm),jf(idm,ms),jl(idm,ms),js(idm) do 1 i=1,ii js(i)=0 do 4 k=1,ms jf(i,k)=0 4 jl(i,k)=0 j=1 k=1 3 if (jpt(i,j).ne.0) go to 2 j=j+1 if (j.le.jj) go to 3 go to 1 2 if (k.gt.ms) then write (lp,'('' error in indxj - ms too small at i,j ='',2i5)') i,j write (lp,'('' i-th line of jpt array:'',/(7(1x,10i1)))') . (jpt(i,l),l=1,jj) stop '(indxj)' end if jf(i,k)=j 6 j=j+1 if (j.le.jj) go to 5 jl(i,k)=jj js(i)=k go to 1 5 if (jpt(i,j).ne.0) go to 6 jl(i,k)=j-1 js(i)=k k=k+1 go to 3 1 continue return end c c c> Revision history: c> c> May 2000 - routine generalized to accomodate cyclic & noncyclic b.c. c> Mar. 2006 - added bering strait exchange logic