elcirc5_01_01c.f90

河口模型 使用模拟盐水入侵、热量扩散等等 河口模型 使用模拟盐水入侵、热量扩散

!...  compute the elements connected to each node 
      do i=1,np
        nne(i)=0
      	nnp(i)=0
      	sarea(i)=0
      enddo

      do i=1,ne
        do j=1,i34(i)
          nd=nm(i,j)
          nne(nd)=nne(nd)+1
          if(nne(nd).gt.mnei) then
            write(11,*)'Too many neighbors',nd
            stop
          endif
          ine(nd,nne(nd))=i
          iself(nd,nne(nd))=j
          sarea(nd)=sarea(nd)+areas(i)
        enddo
      enddo

!...  Check hanging nodes
!...
      ihang=0
      do i=1,np
        if(nne(i).eq.0) then
          ihang=1
          write(11,*)'Hanging node',i
        endif
      enddo
      if(ihang.eq.1) then
        write(11,*)'Check fort.11 for hanging nodes'
        stop
      endif
      
!...  Compute surrounding nodes
!...
      do i=1,np
      	do j=1,nne(i)
          ie=ine(i,j)
          loopk: do k=1,i34(ie)
            nd=nm(ie,k)
            do l=1,nnp(i)
              if(nd.eq.i.or.nd.eq.inp(i,l)) cycle loopk
            enddo !l
            nnp(i)=nnp(i)+1
            if (nnp(i).gt.mnei) then
              write(11,*)'Too many surrounding nodes',i
              stop
            endif
            inp(i,nnp(i))=nd
	  end do loopk !k
      	enddo !j
      enddo !i=1,np

!     Compute ball info
      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)
            iel=ine(nd1,k)
            if(iel/=i.and.(nm(iel,1)==nd2.or.nm(iel,2)==nd2.or.&
     &nm(iel,3)==nd2.or.(i34(iel)==4.and.nm(iel,4)==nd2))) ic3(i,j)=iel
          enddo !k
        enddo !j
      enddo !i

!...  compute the sides information
!...
      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.i<ic3(i,j)) then !new sides
            ns=ns+1 
            if(ns.gt.mns) then
              write(11,*)'Too many sides'
              stop
            endif
            js(i,j)=ns
            is(ns,1)=i
            isidenode(ns,1)=nd1
            isidenode(ns,2)=nd2
            x1j(ns)=x(nd1)
            y1j(ns)=y(nd1)
            x2j(ns)=x(nd2)
            y2j(ns)=y(nd2)
            xcj(ns)=(x1j(ns)+x2j(ns))/2
            ycj(ns)=(y1j(ns)+y2j(ns))/2
            dps(ns)=(dp(nd1)+dp(nd2))/2
            distj(ns)=dsqrt((x2j(ns)-x1j(ns))**2+(y2j(ns)-y1j(ns))**2)
            if(distj(ns).eq.0) then
              write(11,*)'Zero side',ns
              stop
            endif
            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.i<ic3(i,j)
        enddo !j=1,i34(i)
      enddo !i=1,ne

      if(ns.lt.ne.or.ns.lt.np) then
        write(11,*)'Weird grid with ns < ne or ns < np',np,ne,ns
        stop
      endif

      do i=1,ne
        do j=1,i34(i)
          jsj=js(i,j)
          ssign(i,j)=(is(jsj,2)-2*i+is(jsj,1))/(is(jsj,2)-is(jsj,1))
        enddo
      enddo

      if(nscreen.eq.1) write(*,*) 'There are',ns,' sides in the grid...'
      write(16,*) 'There are',ns,' sides in the grid...'

!...  Compute open bnd elements and sides
!...  Each bnd node is assigned 2 flags to account for overlap of 2
!...  open bnds.
      do i=1,np
        isbnode(i,1)=0
        isbnode(i,2)=0
      enddo
      do i=1,nope
        do j=1,nond(i)
          nd=iond(i,j)
          if(isbnode(nd,1).eq.0) then
            isbnode(nd,1)=i
          else if(isbnode(nd,2).eq.0) then
            isbnode(nd,2)=i
          else
            write(11,*)'Illegal open bnd at node',nd
            stop
          endif
        enddo !j
      enddo !i=1,nope

      do j=1,nope
        nosd(j)=0 !initialize
        noe(j)=0
      enddo

      loop1: do i=1,ns
      	if(is(i,2)==0) then
          nd1=isidenode(i,1)
          nd2=isidenode(i,2)
          do j=1,nope
            if((isbnode(nd1,1).eq.j.or.isbnode(nd1,2).eq.j).and.&
     &         (isbnode(nd2,1).eq.j.or.isbnode(nd2,2).eq.j)) then
              nosd(j)=nosd(j)+1
              if(nosd(j).gt.mnosd) then
                write(11,*)'Too many open sides',j
                stop
              endif
              iosd(j,nosd(j))=i
              cycle loop1 !i can only be in one open bnd segment
            endif
          enddo !j=1,nope
        endif !is(i,2)=0
      end do loop1 !i=1,ns

      loop2: do i=1,ne
        do j=1,i34(i)
          nd=nm(i,j)
          do k=1,nope
            if(isbnode(nd,1)==k.or.isbnode(nd,2)==k) then
              noe(k)=noe(k)+1
              if(noe(k).gt.mnoe) then
                write(11,*)'Too many open elements',k
                stop
              endif
              ioe(k,noe(k))=i
!	In the case of dual role, choose any one of the segment
              cycle loop2
            endif
          enddo !k
        enddo !j
      end do loop2 !i=1,ne

!...  compute centers of each element & dpe
!...
      do i=1,ne
	xctr(i)=0
	yctr(i)=0
	dpe(i)=-5.e8
	do j=1,i34(i)
          xctr(i)=xctr(i)+x(nm(i,j))/i34(i)
          yctr(i)=yctr(i)+y(nm(i,j))/i34(i)
	  if(dps(js(i,j))>dpe(i)) dpe(i)=dps(js(i,j))
	enddo !j
      enddo !i=1,ne

!...  Output obe.out, centers.bp and sidecenters.bp if ipre.ne.0
      if(ipre.ne.0) then
        open(32,file='obe.out')
        if(iwrite.eq.0) then
          write(32,*) nope,' # of open bnd'
          write(32,*) 'Element list:'
        else !evm
          write(32,"(i12,a)",advance="no") nope,' # of open bnd\n'
          write(32,"(a)",advance="no") 'Element list:\n'
        endif
        do i=1,nope
          if(iwrite.eq.0) then
            write(32,*) noe(i),' bnd #',i 
          else !evm
            write(32,"(i12,a,i12,a)",advance="no") noe(i),' bnd #',i,'\n'
          endif
          do j=1,noe(i)
            if(iwrite.eq.0) then
              write(32,*) j,ioe(i,j)
            else !evm
              write(32,"(2i12,a)",advance="no") j,ioe(i,j),'\n'
            endif
          enddo !j
        enddo !i
        if(iwrite.eq.0) then
          write(32,*) 'Side list:'
        else !evm
          write(32,"(a)",advance="no") 'Side list:\n'
        endif
        do i=1,nope
          if(iwrite.eq.0) then
            write(32,*) nosd(i),' bnd #',i 
          else !evm
            write(32,"(i12,a,i12,a)",advance="no") nosd(i),' bnd #',i,'\n'
          endif
          do j=1,nosd(i)
            if(iwrite.eq.0) then
              write(32,*) j,iosd(i,j)
            else !evm
              write(32,"(2i12,a)",advance="no") j,iosd(i,j),'\n'
            endif
          enddo
        enddo !i
        close(32)

        open(32,file='sidecenters.bp')
        if(iwrite.eq.0) then
          write(32,*) 'Sidegrid'
          write(32,*) ns
        else !evm
          write(32,"(a)",advance="no") 'Sidegrid\n'
          write(32,"(i12,a)",advance="no") ns,'\n'
        endif
        do i=1,ns
          if(iwrite.eq.0) then
            write(32,*) i,xcj(i),ycj(i),real(dps(i))
          else !evm
            write(32,"(i12,a,f19.9,a,f19.9,a,f12.6,a)",advance="no") i, &
     &" ",xcj(i),"      ",ycj(i),"    ",real(dps(i)),'\n'
          endif
        enddo !i
        close(32)

        open(32,file='centers.bp')
        if(iwrite.eq.0) then
          write(32,*) 'centers pts'
          write(32,*) ne
        else !evm
          write(32,"(a)",advance="no") 'centers pts\n'
          write(32,"(i12,a)",advance="no") ne,'\n'
        endif
        do i=1,ne
          if(iwrite.eq.0) then
            write(32,*) i,xctr(i),yctr(i),real(dpe(i))
          else !evm
            write(32,"(i12,a,f19.9,a,f19.9,a,f12.6,a)",advance="no") i, &
     &" ",xctr(i),"      ",yctr(i),"    ",real(dpe(i)),'\n'
          endif
        enddo !i
        close(32)

        write(*,*)'obe.out successfully generated; bye'
        write(16,*)'obe.out successfully generated; bye'
        stop
      endif !ipre.ne.0

!...  compute angles (orientation of local x) and perpendicular distances for sides 
!...
      do j=1,ns
        if(is(j,2)/=0) then
          delj(j)=dsqrt((xctr(is(j,2))-xctr(is(j,1)))**2+(yctr(is(j,2))-yctr(is(j,1)))**2)
          if(delj(j).eq.0) then
            write(*,*)'Zero distance between centers',is(j,1),is(j,2)
            stop
          endif
        endif
        thetan=datan2(x1j(j)-x2j(j),y2j(j)-y1j(j))
      	snx(j)=dcos(thetan)
      	sny(j)=dsin(thetan)
      enddo

      if(nscreen.eq.1) write(*,*)'done computing geometry...'
      write(16,*)'done computing geometry...'

!...  Compute channel widths for type I b.c.
      do i=1,nope
        cwidth(i)=0
        do j=1,nosd(i)
          cwidth(i)=cwidth(i)+distj(iosd(i,j))
        enddo !j
      enddo !i

!...  classify boundary elements & sides
!...
      do i=1,nope
        do j=1,noe(i)
          isbe(ioe(i,j))=i
        enddo !j
        do j=1,nosd(i)
          isbs(iosd(i,j))=i
        enddo !j=1,nosd(i)
      enddo !i=1,nope

      do i=1,ns
      	if(is(i,2)==0.and.isbs(i)==0) isbs(i)=-1
      enddo

      if(nscreen.eq.1) write(*,*)'done classifying boundaries...'
      write(16,*)'done classifying boundaries...'

!...  Continue to read fort.15
!...

!... Implicitness for momentum
      read(15,*)thetai

!... Baroclinic flags
      read(15,*) ibc,ibtp
      if(ibc.ne.0.and.ibc.ne.1) then
      	write(*,*)'Unknown ibc'
      	stop
      endif
      if(ibtp.ne.0.and.ibtp.ne.1) then
      	write(*,*)'Unknown ibtp'
     	stop
      endif

      if(ibc.eq.0) then
        write(*,*)'You are using baroclinic model'
        write(16,*)'You are using baroclinic model'
      	read(15,*) nrampbc,drampbc 
      else !ibc=1
        if(ibtp.eq.0) then
          write(*,*)'Barotropic model without ST calculation'
          write(16,*)'Barotropic model without ST calculation'
        else !ibtp=1
          write(*,*)'Barotropic model with ST calculation'
          write(16,*)'Barotropic model with ST calculation'
        endif
      endif

      read(15,*) tempmin,tempmax,saltmin,saltmax
      if(tempmin.lt.0.or.tempmax.gt.40.or.saltmin.lt.0.or.saltmax.gt.42) then
  	write(*,*)'Specified ST range invalid'
      	stop
      endif

      read(15,*) rnday

!... dramp not used if nramp=0
      read(15,*) nramp,dramp
      if(nramp.ne.0.and.nramp.ne.1) then
      	write(*,*)'Unknown nramp',nramp
      	write(11,*)'Unknown nramp',nramp
        stop
      endif

      read(15,*) dt

!...  compute total number of time steps nt
      nt=rnday*86400.d0/dt+0.5

!...  input info on backtracking
!...
      read(15,*) nsubfl !flag
      if(nsubfl.lt.0.or.nsubfl.gt.2) then
        write(11,*)'Unknown subdivision flag',nsubfl
      	stop
      endif

      if(nsubfl.eq.0) then
        read(15,*) nsub
      	do i=1,mnp
          do j=1,nvrt
            ndelt(i,j)=nsub
          enddo
      	enddo
      else if(nsubfl.eq.1) then
      	open(41,file='ndelt.gr3',status='old')
        read(41,*)
        read(41,*) ntmp,npbp
        do i=1,npbp
          read(41,*)j,xtmp,ytmp,sub
          if(sub.le.0) then
            write(11,*)'# of subdivisions in btrack <=0',sub,i
            stop
          endif
          do j=1,nvrt
            ndelt(i,j)=sub
          enddo
        enddo !i
      	close(41)
      else !nsubfl=2
        read(15,*) ndelt_min,ndelt_max
      	if(ndelt_min.le.0.or.ndelt_min.gt.ndelt_max) then
          write(11,*)'Incorrect input for variable tracking!'
          write(11,*)ndelt_min,ndelt_max
          stop
        endif
      endif

!... Advection flag for momentum eq.
      read(15,*) nadv !flag
      if(nadv.ne.0.and.nadv.ne.1) then
        write(11,*)'Unknown advection flag',nadv
      	stop
      endif