!This script is used to post-process the NWM coupling files. !It combines the sinks with adjacent sources if there is any. !Two options are available: distance-based and neighboring element based. !Inputs files: ! inbr: 0 (distance-based); 1 (neighboring element based) ! dist: distance if inbr; number of neighboring tiers ! !Input files ! hgrid.cpp : hgrid in cpp projection ! source_sink.in : contains the element ID for each ! intersection of the NWM streams and the land boundary. ! msource.th : contains the salinity and temprature ! of the source element along the land boundary. ! Salinity is set to be 0, temp = -9999. ! vsource.th : input of the stream flows of source elements. ! vsink.th : input of the stream flows of sink elements. ! !Output files ! source_sink.in.1, msource.th.1, vsource.th.1, vsink.th.1 ! !serial: !ifort -CB -O2 -o combine_sink_source combine_sink_source.F90 !pgf90 -O2 -o combine_sink_source combine_sink_source.F90 ! !openmp: !ifort -CB -O2 -qopenmp -o combine_sink_source combine_sink_source.F90 !pgf90 -O2 -mp -o combine_sink_source combine_sink_source.F90 !export OMP_NUM_TREADS=8 !setenv OMP_NUM_TREADS 8 program coupling_nwm use omp_lib implicit none real(8) :: start_time,end_time integer,parameter::max_nei=1000 real(8) :: tmp,distance,dist,max_dist real(8), allocatable :: temp(:),xx(:),yy(:),xel(:),yel(:) integer, allocatable :: elnode(:,:),i34(:),indel(:,:),n_sink_source(:), & & i_sink_source(:,:),ntier1(:),ele_source(:),ele_sink(:),imap_ele2source(:), & & tier1(:,:),tier_n(:,:),ntier_n(:),nne(:),nxq(:,:,:),ic3(:,:),isbnd(:), & & ncount(:),nlbnd(:),lbnd(:,:),nobnd(:),obnd(:,:),i_island(:) integer :: mnei,inbr,i,j,k,nsource,nsink,itmp,nt,istat,nd,ntracer,n_tier_dist integer :: ne,np,id,id1,ie,ie_nbr,isource,kk,k0,nope,nland,nn,n1,n2,l, & & ip,je,new,ii,icount,tid,nthreads,ntmp real(8), allocatable :: vsource(:,:),vsink(:,:),msource(:,:,:), & & time_stamp(:) real(8) :: t,s t = -9999 s = 0 !$omp parallel !tid = omp_get_thread_num() !if (tid==0) then ! nthreads= omp_get_num_threads() ! print*,'Number 0f threads=',nthreads ! allocate(ncount(0:nthreads-1),stat=istat) ! if (istat/=0) stop 'Failed to alloc. ncount' !endif !$omp end parallel !Read in cmd line inputs print*, 'Input search option (0: distance-based; 1: neighboring element based):' read*, inbr !inbr=1 print*, 'Input search radius (if inbr=0: distance; if inbr=1: neighboring tiers):' if (inbr==0) then read*, distance elseif (inbr==1) then read*, n_tier_dist else print*, 'wrong inbr ', inbr stop endif !distance=1 !read hgrid !call cpu_time(start_time) !start_time=omp_get_wtime(); print*, 'reading inputs ...' open(14,file='hgrid.cpp',status='old')!lambert projection, NWM has shorter precision read(14,*) read(14,*)ne,np write(*,*)'# of elements',ne,'# of nodes',np allocate(xx(np),yy(np),nne(np),stat=istat) if(istat/=0) stop 'Failed to alloc. xx, yy, nne' allocate(i34(ne),elnode(4,ne),stat=istat) if(istat/=0) stop 'Failed to alloc. i34, elnode' do i=1,np read(14,*)j,xx(i),yy(i),tmp enddo do i=1,ne read(14,*) j,i34(i),elnode(1:i34(i),i) enddo if (inbr==1) then allocate(isbnd(np),stat=istat) if (istat/=0) stop 'Failed to alloc. isbnd' isbnd=0 read(14,*) nope read(14,*) ntmp allocate(obnd(ntmp,nope),nobnd(nope),stat=istat) if (istat/=0) stop 'Failed to alloc. obnd, nobnd' do k=1,nope read(14,*) nobnd(k) do i=1,nobnd(k) read(14,*) ip obnd(i,k)=ip isbnd(ip)=1 !open bnd enddo !i enddo !k read(14,*) nland read(14,*) ntmp allocate(lbnd(ntmp,nland),nlbnd(nland),i_island(nland),stat=istat) if (istat/=0) stop 'Failed to alloc. lbnd, nlbnd, i_island' do k=1,nland read(14,*) nlbnd(k), i_island(k) do i=1,nlbnd(k) read(14,*) ip lbnd(i,k)=ip if (isbnd(ip)==1) then !open bnd isbnd(ip)=-1 !open and land bnd else isbnd(ip)=1 !land bnd endif enddo !i enddo !k endif close(14) ! read existing source/sink files open(13,file='source_sink.in',status='old') read(13,*) nsource allocate(ele_source(nsource),stat=istat) if (istat/=0) stop 'Failed to alloc. ele_source' do i=1,nsource read(13,*)ele_source(i) enddo read(13,*) read(13,*) nsink allocate(ele_sink(nsink),stat=istat) if (istat/=0) stop 'Failed to alloc. ele_sink' do i=1,nsink read(13,*)ele_sink(i) enddo close(13) nt=0 open(15,file='vsource.th',status='old') do while (.true.) read (15,*,end=999) nt=nt+1 enddo 999 continue close(15) allocate(vsource(nsource,nt),time_stamp(nt),stat=istat) if (istat/=0) stop 'Failed to alloc. time_stamp' open(15,file='vsource.th',status='old') do i=1,nt read (15,*) time_stamp(i),vsource(:,i) enddo close(15) !ntracer=2 !allocate(msource(nsource,ntracer,nt),stat=istat) !if (istat/=0) stop 'Failed to alloc. msource' !open(15,file='msource.th',status='old') !do i=1,nt ! read (15,*) tmp,msource(:,:,i) !enddo !close(15) allocate(vsink(nsink,nt),stat=istat) if (istat/=0) stop 'Failed to alloc. vsink' open(15,file='vsink.th',status='old') do i=1,nt read (15,*) tmp,vsink(:,i) enddo close(15) !call cpu_time(end_time) !end_time=omp_get_wtime(); !print*, 'reading inputs takes ',(end_time-start_time)/60d0,' minutes' !------done reading inputs----------------------------------- !call cpu_time(start_time) !start_time=omp_get_wtime(); print*, 'calculating geometry ...' if (inbr==1) then ! elem ball nne=0 do i=1,ne do j=1,i34(i) nd=elnode(j,i) nne(nd)=nne(nd)+1 enddo enddo mnei=maxval(nne) allocate(indel(mnei,np),stat=istat) if(istat/=0) stop 'Failed to alloc. indel' nne=0 do i=1,ne do j=1,i34(i) nd=elnode(j,i) nne(nd)=nne(nd)+1 if(nne(nd)>mnei) then write(*,*)'Too many neighbors',nd stop endif indel(nne(nd),nd)=i enddo enddo !i allocate(nxq(3,4,4),ic3(4,ne),stat=istat) if (istat/=0) stop 'Failed to alloc. nxq, ic3' !re-order indel(:,nd) in counter-clockwise order !setup nxq (cyclic node index) do k=3,4 !elem. type do i=1,k !local index do j=1,k-1 !offset nxq(j,i,k)=i+j if(nxq(j,i,k)>k) nxq(j,i,k)=nxq(j,i,k)-k if(nxq(j,i,k)<1.or.nxq(j,i,k)>k) then write(*,*)'INIT: nx wrong',i,j,k,nxq(j,i,k) stop endif enddo !j enddo !i enddo !k !ele-side-ele table do ie=1,ne do k=1,i34(ie) ic3(k,ie)=0 !index for boundary sides n1=elnode(nxq(1,k,i34(ie)),ie) n2=elnode(nxq(2,k,i34(ie)),ie) do l=1,nne(n1) je=indel(l,n1) if(je/=ie.and.(elnode(1,je)==n2.or.elnode(2,je)==n2.or. & &elnode(3,je)==n2.or.(i34(je)==4.and.elnode(4,je)==n2))) ic3(k,ie)=je enddo !l je=ic3(k,ie) if(je/=0) then do l=1,i34(je) if(elnode(nxq(1,l,i34(je)),je)==n1.and.elnode(nxq(2,l,i34(je)),je)==n2) then write(*,*) 'Elem ', ie, ' and ', je, ' have opposite orientation' stop endif end do !l endif enddo !k enddo !ie do i=1,np if(isbnd(i)/=0) then !bnd ball ! Look for starting bnd element icount=0 do j=1,nne(i) ie=indel(j,i) ii=0 !local index do l=1,i34(ie) if(elnode(l,ie)==i) then ii=l; exit endif enddo !l if(ii==0) stop('AQUIRE_HGRID: bomb (1)') if(ic3(nxq(i34(ie)-1,ii,i34(ie)),ie)==0) then icount=icount+1 indel(1,i)=ie endif enddo !j=1,nne if(icount/=1) then write(*,*)'Illegal bnd node',i,isbnd(i),icount stop endif endif !bnd ball ! For internal balls, starting elem. is not altered ! Sequential search for the rest of elements do j=2,nne(i) ie=indel(j-1,i) ii=0 !local index do l=1,i34(ie) if(elnode(l,ie)==i) then ii=l; exit endif enddo !l if(ii==0) stop('AQUIRE_HGRID: bomb (2)') new=ic3(nxq(i34(ie)-2,ii,i34(ie)),ie) if(new==0) then write(*,*)'Incomplete ball',i stop endif indel(j,i)=new enddo !j=2,nne(i) enddo !i=1,np endif !--------------------------------------------------------------- ! find neighboring sources for each sink !--------------------------------------------------------------- allocate(n_sink_source(nsink),i_sink_source(max_nei,nsink),stat=istat) if (istat/=0) stop 'Failed to alloc. n_sink_source, i_sink_source' n_sink_source=0; i_sink_source=0 !--------------------------------------------------------------- if (inbr==0) then !distance-based !--------------------------------------------------------------- ! elem coordinates allocate(xel(ne),yel(ne),stat=istat) if (istat/=0) stop 'Failed to alloc. xel, yel' xel=0d0; yel=0d0 do i=1,ne xel(i)=sum(xx(elnode(1:i34(i),i)))/dble(i34(i)) yel(i)=sum(yy(elnode(1:i34(i),i)))/dble(i34(i)) enddo !<<<<<<< Output sources/sinks as bpfiles<<<<<<<<< open(24,file='sources_sinks.bp',status='replace') write(24,*) write(24,*) nsource+nsink do i=1,nsource ie=ele_source(i) write(24,*) i,xel(ie),yel(ie),ie enddo do i=1,nsink ie=ele_sink(i) write(24,*) i,xel(ie),yel(ie),-ie enddo close(24) open(24,file='sources.bp',status='replace') write(24,*) write(24,*) nsource do i=1,nsource ie=ele_source(i) write(24,*) i,xel(ie),yel(ie),ie enddo close(24) open(24,file='sinks.bp',status='replace') write(24,*) write(24,*) nsink do i=1,nsink ie=ele_sink(i) write(24,*) i,xel(ie),yel(ie),ie enddo close(24) print*, 'done outputing bpfiles for sources/sinks' !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> max_dist=distance**2 do i=1,nsink ie=ele_sink(i) do j=1,nsource ie_nbr=ele_source(j) dist=(xel(ie)-xel(ie_nbr))**2+(yel(ie)-yel(ie_nbr))**2 if (dist0) then !self has both sink and source n_sink_source(i)=n_sink_source(i)+1 i_sink_source(n_sink_source(i),i)=imap_ele2source(ie) endif call mark_bnd_neighbors(i,ie,0,n_tier_dist) endif enddo !$omp end parallel do !--------------------------------------------------------------- elseif (inbr==2) then !neighboring element tier in 2D !--------------------------------------------------------------- allocate(tier1(0:mnei*4,ne),ntier1(ne),stat=istat) !tier1(0) should be 0 at all time if (istat/=0) stop 'Failed to alloc. tier1, ntier1' allocate(tier_n(0:mnei*400,nsink),ntier_n(nsink),stat=istat) if (istat/=0) stop 'Failed to alloc. tier_n, ntier_n' ncount=0 !$omp parallel do private (ie,j,k,n1,k0,kk,tid) do ie=1,ne !tid = omp_get_thread_num() !ncount(tid)=ncount(tid)+1 ntier1(ie)=0 !number of tier 1 elements do j=1,i34(ie) n1=elnode(j,ie) !find position of ie in the nodal ball do k=1,nne(n1) if (indel(k,n1).eq.ie) then k0=k exit endif enddo !k do k=1,nne(n1)-1 !counter-clockwise: from the first non-ie element to the last non-ie element !thus, the possible overlapping of tier1 elements (from next node of ie) can only !occur at the first or last element in the current tier1 kk=k+k0; if (kk>nne(n1)) kk=kk-nne(n1) if (indel(kk,n1).ne.0 .and. indel(kk,n1).ne.tier1(1,ie) .and. indel(kk,n1).ne.tier1(ntier1(ie),ie)) then ntier1(ie)=ntier1(ie)+1 if (ntier1(ie)<=mnei*4) then tier1(ntier1(ie),ie)=indel(kk,n1) else write(*,*) 'tier 1 > mnei*4' stop endif endif enddo !k enddo enddo !ie !$omp end parallel do print*, ncount !mark sources allocate(imap_ele2source(ne),stat=istat) if (istat/=0) stop 'Failed to alloc. imap_ele2source' imap_ele2source=0 do i=1,nsource imap_ele2source(ele_source(i))=i enddo tier_n=0; ntier_n=0 !$omp parallel do private (i,ie) do i=1,nsink ie=ele_sink(i) tier_n(0,i)=ie !self call mark_neighbors(i,ie,0,n_tier_dist) enddo !$omp end parallel do !call cpu_time(end_time) !end_time=omp_get_wtime(); !print*, 'calculating geometry takes ',(end_time-start_time)/60d0,' minutes' !--------------------------------------------------------------- endif !inbr !--------------------------------------------------------------- print*, 'redistributing sinks ...' !debug open(33,file='i_sink_source.txt',status='replace') do i=1,nsink write(33,*) 'Sink #',i,' at Ele #',ele_sink(i), ' has ',n_sink_source(i),'neighboring sources' write(33,*) i_sink_source(1:n_sink_source(i),i) write(33,*) enddo ! Redistribute vsink to neighboring vsources ! Do not use omp on this loop, serial mode assumes smaller i is ! treated (distributed to nearby sources) first, which is not true for omp. ! This step takes little time anyway. do i=1,nsink if (n_sink_source(i)>0) then do k=1,nt do j=1,n_sink_source(i) isource=i_sink_source(j,i) if (vsource(isource,k)>abs(vsink(i,k))) then vsource(isource,k)=vsource(isource,k)+vsink(i,k) vsink(i,k)=0d0 exit !no need for finding next neighboring source else vsink(i,k)=vsink(i,k)+vsource(isource,k) vsource(isource,k)=0.0d0 endif enddo ! j=1,nsource enddo !k=1,nt endif enddo !i=1,nsink close(33) do i=1,nsink if (n_sink_source(i)>0) then do k=1,nt do j=1,n_sink_source(i) isource=i_sink_source(j,i) if (vsource(isource,k)>abs(vsink(i,k))) then vsource(isource,k)=vsource(isource,k)+vsink(i,k) vsink(i,k)=0d0 exit !no need for finding next neighboring source else vsink(i,k)=vsink(i,k)+vsource(isource,k) vsource(isource,k)=0.0d0 endif enddo ! j=1,nsource enddo !k=1,nt endif enddo !i=1,nsink !call cpu_time(end_time) print*, 'Done redistributing sinks' ! write vsource.th file open(16,file='vsource.th.1',status='replace') do i=1,nt write(16,'(10000F15.3)')time_stamp(i),(vsource(k,i),k=1,nsource) end do ! write vsink.th open(17,file='vsink.th.1') do i=1,nt write(17,'(10000F15.3)')time_stamp(i),(vsink(k,i),k=1,nsink) end do ! write msource.th open(18,file='msource.th.1') do i=1,nt write(18,'(10000F15.3)')time_stamp(i),(t,k=1,nsource),(s,k=1,nsource) end do ! write source_sink.in open(19,file='source_sink.in.1') write(19,*)nsource do i=1,nsource write(19,*)ele_source(i) end do write(19,*) write(19,*)nsink do i=1,nsink write(19,*)ele_sink(i) end do contains recursive subroutine mark_bnd_neighbors(isink,ie_ctr,i_depth,n_tier_dist) implicit none integer, intent(in) :: isink, ie_ctr, i_depth,n_tier_dist integer :: iter,ie000,kk,ie1,ie2,nd,nd1,nd2,ie3,ie_new(2),nn,i logical :: inew if (i_depth >= n_tier_dist) then return endif !find tier 1 neighbors, at most 2, at least 1 (when self is quad) !the following searching procedure didn't stop at non-land bnd elements. !This is okay since non-land bnd elements do not have sources anyway nn=0 do i=1,i34(ie_ctr) nd = elnode(i,ie_ctr) if (isbnd(nd)==0) cycle !internal nodes ie1=indel(1,nd) ie2=indel(nne(nd),nd) if (ie1==ie_ctr .and. ie2==ie_ctr) then cycle !tip of a small stream elseif (ie1==ie_ctr .and. ie2.ne.ie_ctr) then nn=nn+1 if (nn>2) then print*, 'tier 1 neighbors > 2: ', isink,ie_ctr,i_depth,ie1,ie2 endif ie_new(nn)=ie2 elseif (ie1.ne.ie_ctr .and. ie2==ie_ctr) then nn=nn+1 if (nn>2) then print*, 'tier 1 neighbors > 2: ', isink,ie_ctr,i_depth,ie1,ie2 endif ie_new(nn)=ie1 else !could happen in narrow 2DV channel, where the other node is also a bnd !node but on the other side of the channel cycle endif enddo !push new elements in tier n into record loop1: do iter=1,nn ie000=ie_new(iter) inew=.true. do kk=0,ntier_n(isink) if (ie000==tier_n(kk,isink)) then inew=.false. !avoid duplicated points, but its neighbor may be new exit endif enddo !register new neighbor if (inew) then ntier_n(isink)=ntier_n(isink)+1 tier_n(ntier_n(isink),isink)=ie000 !debug !write(*,*) 'ie_ctr,i_depth,ntier_n(isink),ie000:',ie_ctr,i_depth,ntier_n(isink),ie000 if (imap_ele2source(ie000)>0) then n_sink_source(isink)=n_sink_source(isink)+1 i_sink_source(n_sink_source(isink),isink)=imap_ele2source(ie000) endif endif call mark_bnd_neighbors(isink,ie000,i_depth+1,n_tier_dist) enddo loop1 end subroutine mark_bnd_neighbors recursive subroutine mark_neighbors(isink,ie_ctr,i_depth,n_tier_dist) implicit none integer, intent(in) :: isink, ie_ctr, i_depth,n_tier_dist integer :: iter,ie000,kk logical :: inew if (i_depth >= n_tier_dist) then return endif loop1: do iter=1,ntier1(ie_ctr) ie000=tier1(iter,ie_ctr) inew=.true. do kk=0,ntier_n(isink) if (ie000==tier_n(kk,isink)) then inew=.false. !avoid duplicated points, but its neighbor may be new exit endif enddo !register new neighbor if (inew) then ntier_n(isink)=ntier_n(isink)+1 tier_n(ntier_n(isink),isink)=ie000 !debug !write(*,*) 'ie_ctr,i_depth,ntier_n(isink),ie000:',ie_ctr,i_depth,ntier_n(isink),ie000 if (imap_ele2source(ie000)>0) then n_sink_source(isink)=n_sink_source(isink)+1 i_sink_source(n_sink_source(isink),isink)=imap_ele2source(ie000) endif endif call mark_neighbors(isink,ie000,i_depth+1,n_tier_dist) enddo loop1 end subroutine mark_neighbors end