! Generate bnd list for GNOME: nbe(ne,3) (elem ball); bnd(:,4) - 2 end nodes of a bnd seg, ! followed by 2 flags. First flag indicates if the seg is on exterior ! bnd (0) or island (1,2,...). The second flag indicates if it's an open ocean ! bnd (0: no (inclduing river bnd); 1: yes). ! GNOME only works for pure triangular grids. ! ! ifort -O2 -Bstatic -CB -g -traceback -o gen_GNOME_bnd gen_GNOME_bnd.f90 ! Input: ! (1) hgrid.gr3 with bnd info, generated with gredit. No open bnd's ! on islands; (2) screen input ! Output: nbe,out; GNOME_bnd_ext.out (exterior bnd) and GNOME_bnd_island.out (island); debug: fort.98 ! Then: cat GNOME_bnd_ext.out GNOME_bnd_island.out > GNOME_bnd.out (i.e. bnd(:,4)) ! Use gen_GNOME_current.m to generate nc input to GNOME ! implicit none ! integer, parameter :: debug=1 integer, parameter :: mnp=500000 integer, parameter :: mne=1000000 integer, parameter :: mns=1500000 integer, parameter :: mnope=10 !max # of open bnd segments integer, parameter :: mnond=1000 !max # of open bnd nodes in each segment integer, parameter :: mnei=20 !neighbor integer, parameter :: nbyte=4 real, parameter :: small1=1.e-3 !used to check area ratios ! Vertical postion, salinity, and temperature integer, dimension(:,:), allocatable :: kbp integer :: ier ! allocate error return. dimension xnd(mnp),ynd(mnp),nm(mne,4),dp(mnp),i34(mne) dimension iest(mnp),ixy(mnp,2),arco(3) dimension wild(100),wild2(100,2) dimension nne(mnp),ine(mnp,mnei),ic3(mne,4),nx(4,4,3),js(mne,4),is(mns,2),isidenode(mns,2) dimension xcj(mns),ycj(mns),nond(mnope),iond(mnope,mnond),isbnd(mnp) dimension ie_land(mne),icolor(mnp),list_ext(mnp) print*, 'Input open bnd segment # corresponding to open ocean bnd:' !' read*, iocean ! Read in hgrid and vgrid open(13,file='GNOME_bnd_ext.out',status='replace') open(15,file='GNOME_bnd_island.out',status='replace') open(17,file='nbe.out',status='replace') open(14,file='hgrid.gr3',status='old') read(14,*) read(14,*)ne,np if(np.gt.mnp.or.ne.gt.mne) then write(*,*)'Increase mnp/mne' stop endif do i=1,np read(14,*)j,xnd(i),ynd(i),dp(i) enddo !i do i=1,ne read(14,*)j,i34(i),(nm(i,l),l=1,i34(i)) n1=nm(i,1) n2=nm(i,2) n3=nm(i,3) if(i34(i)==3) then area=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3)) else !quad n4=nm(i,4) ! Check convexity ar1=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3)) ar2=signa(xnd(n1),xnd(n3),xnd(n4),ynd(n1),ynd(n3),ynd(n4)) ar3=signa(xnd(n1),xnd(n2),xnd(n4),ynd(n1),ynd(n2),ynd(n4)) ar4=signa(xnd(n2),xnd(n3),xnd(n4),ynd(n2),ynd(n3),ynd(n4)) if(ar1.le.0.or.ar2.le.0.or.ar3.le.0.or.ar4.le.0) then write(*,*)'Concave quadrangle',i,ar1,ar2,ar3,ar4 stop endif area=ar1+ar2 endif if(area<=0) then write(*,*)'Negative area at',i stop endif write(12,*)i,area enddo !i=1,ne close(12) ! Open bnds: isbnd() ! 0: internal; >0: open bnd segmt #; -1: exterior land bnd (not on ! any open bnd); -2: island bnd nodes isbnd=0 read(14,*) nope read(14,*) neta ntot=0 if(nope>mnope) stop 'Increase mnope (2)' do k=1,nope read(14,*) nond(k) if(nond(k)>mnond) stop 'Increase mnond' do i=1,nond(k) read(14,*) iond(k,i) isbnd(iond(k,i))=k enddo !i enddo !k !Init list of exterior bnd nodes nlist_ext=0 if(nond(1)>0) then nlist_ext=nond(1) !# of exterior bnd nodes in the list now list_ext(1:nlist_ext)=iond(1,1:nlist_ext) endif ! Land bnds read(14,*) nland read(14,*) nvel !total # nseg=0 !final # of bnd segments (exterior counted as a single bnd) do k=1,nland read(14,*)nlnd,ifl do i=1,nlnd read(14,*)ilnd icolor(i)=ilnd !temp save if(isbnd(ilnd)==0) then if(ifl==0) then !exterior isbnd(ilnd)=-1 else !island isbnd(ilnd)=-2 endif endif enddo !i if(nlist_ext==0.and.ifl==0) then nlist_ext=nlnd list_ext(1:nlist_ext)=icolor(1:nlnd) endif !Output islands if(ifl==1) then icolor(nlnd+1)=icolor(1) !cyclic nseg=nseg+1 do i=1,nlnd write(15,*)icolor(i),icolor(i+1),nseg,0 enddo !i endif !ifl enddo !k close(14) print*, 'List of exterior bnd init:',nlist_ext,list_ext(1:nlist_ext) ! Compute geometry do k=3,4 do i=1,k do j=1,k-1 nx(k,i,j)=i+j if(nx(k,i,j)>k) nx(k,i,j)=nx(k,i,j)-k if(nx(k,i,j)<1.or.nx(k,i,j)>k) then write(*,*)'nx wrong',i,j,k,nx(k,i,j) stop endif enddo !j enddo !i enddo !k do i=1,np nne(i)=0 enddo do i=1,ne do j=1,i34(i) nd=nm(i,j) nne(nd)=nne(nd)+1 if(nne(nd)>mnei) then write(*,*)'Too many neighbors',nd stop endif ine(nd,nne(nd))=i enddo enddo ! Compute ball info; this won't be affected by re-arrangement below do i=1,ne do j=1,i34(i) ic3(i,j)=0 !index for bnd sides nd1=nm(i,nx(i34(i),j,1)) nd2=nm(i,nx(i34(i),j,2)) do k=1,nne(nd1) ie=ine(nd1,k) if(ie/=i.and.(nm(ie,1)==nd2.or.nm(ie,2)==nd2.or.nm(ie,3)==nd2.or.(i34(ie)==4.and.nm(ie,4)==nd2))) ic3(i,j)=ie enddo !k enddo !j write(17,*)i,ic3(i,1:3) enddo !i ns=0 !# of sides do i=1,ne do j=1,i34(i) nd1=nm(i,nx(i34(i),j,1)) nd2=nm(i,nx(i34(i),j,2)) if(ic3(i,j)==0.or.imns) then write(*,*)'Too many sides' stop endif js(i,j)=ns is(ns,1)=i isidenode(ns,1)=nd1 isidenode(ns,2)=nd2 xcj(ns)=(xnd(nd1)+xnd(nd2))/2 ycj(ns)=(ynd(nd1)+ynd(nd2))/2 is(ns,2)=ic3(i,j) !bnd element => bnd side ! Corresponding side in element ic3(i,j) if(ic3(i,j)/=0) then !old internal side iel=ic3(i,j) index=0 do k=1,i34(iel) if(ic3(iel,k)==i) then index=k exit endif enddo !k if(index==0) then write(*,*)'Wrong ball info',i,j stop endif js(iel,index)=ns endif !ic3(i,j).ne.0 endif !ic3(i,j)==0.or.inp) stop 'overflow(0)' list_ext(nlist_ext)=nd exit loop1 endif if(icolor(nd)==0.and.(isbnd(nd)>0.or.isbnd(nd)==-1)) then !new front on exterior lflound=1 icolor(nd)=1 nlist_ext=nlist_ext+1 if(nlist_ext>np) stop 'overflow(1)' list_ext(nlist_ext)=nd ifront=nd write(98,*)'new front node:',nd cycle loop1 endif enddo !j enddo !m if(lflound==0) then write(*,*)'Unable to find next front node:',ifront stop endif enddo loop1 write(98,*)nlist_ext,' nodes on exterior bnd' write(98,*)list_ext(1:nlist_ext) do i=1,nlist_ext-1 n1=list_ext(i); n2=list_ext(i+1) if(isbnd(n1)==iocean.and.isbnd(n2)==iocean) then itmp=1 else itmp=0 endif write(13,*)n1,n2,0,itmp enddo !i stop end function signa(x1,x2,x3,y1,y2,y3) !... Compute signed area formed by pts 1,2,3 ! implicit real*8(a-h,o-z) signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2 return end