elcirc5_01_01c.f90

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

          do it=1,iths
            read(51,*) ttt,qq
          enddo !it
        endif

      	if(ntetype>0) then
          do it=1,iths
            read(52,*) ttt,te
          enddo !it
      	endif

      	if(nsatype>0) then
          do it=1,iths
            read(53,*) ttt,sal
          enddo !it
      	endif

!...  end hot start section
!...
      endif !ihot.ne.0

      if(nscreen.eq.1) write(*,*)'done initializing variables...'
      write(16,*)'done initializing variables...'

!                                                                             *
!******************************************************************************
!                                                                             *
!			open output files				      *
!                                                                             *
!******************************************************************************
!                                                                             *

!...  write global output headers
!...

      if(ihot<=1) then
   	ifile=1 !output file #
!       Convert it to a string
      	write(ifile_char,'(i12)') ifile
      endif
      call header

      if(nscreen.eq.1) write(*,*)'done initializing outputs'
      write(16,*)'done initializing outputs'

!                                                                             *
!******************************************************************************
!                                                                             *
!	    	          Time stepping 				      *
!                                                                             *
!******************************************************************************
!                                                                             *

      if(ihot.eq.0) iths=0
      if(nscreen.eq.1) write(*,*)'time stepping begins...',iths+1,nt
      write(16,*)'time stepping begins...',iths+1,nt

!...  Compute initial levels
      call levels(iths,iths,itur)
      if(nscreen.eq.1) write(*,*)'done computing initial levels...'
      write(16,*)'done computing initial levels...'

!...  Compute total # of active faces
      nfaces=0
      do i=1,ns
        if(kfs(i)/=0) nfaces=nfaces+kfs(i)-kbs(i)+1
      enddo !i=1,ns
      if(nscreen.eq.1) write(*,*)'Total # of faces=',nfaces
      write(16,*)'Total # of faces=',nfaces

!...  Compute nodal vel. 
      call nodalvel
      if(nscreen.eq.1) write(*,*)'done computing initial nodal vel...'
      write(16,*)'done computing initial nodal vel...'

!...  Compute initial density
      call eqstate
      if(nscreen.eq.1) write(*,*)'done computing initial density...'
      write(16,*)'done computing initial density...'

!...  Store initial elevation field (hot start) for relaxation in sponge layer
!...
      do i=1,ne
        etaic(i)=eta1(i)
      enddo

!...  For UB & MYG, compute vertical diffusivities at whole levels 
      if(itur.eq.3) then
        do i=1,ns
    	  do j=kbs(i),kfs(i) !wet sides
	    call asm(g,i,j,vd,td,qd,qd2)
            vdiff(i,j)=dmin1(diffmax(i),dmax1(bgiff,vd))
            tdiff(i,j)=dmin1(diffmax(i),dmax1(bgdiff,td))
            qdiff(i,j)=dmin1(diffmax(i),dmax1(bgdiff,qd))
            qdiff2(i,j)=dmin1(diffmax(i),dmax1(bgdiff,qd2))
          enddo !j=kbs(i),kfs(i)
        enddo !i=1,ns
      endif !itur=3

!...  Initialize heat budget model
      if(nws==2.and.ihconsv/=0) then
        call surf_fluxes(wtime1,windx1,windy1,pr1,airt1,shum1,&
     &srad,fluxsu,fluxlu,hradu,hradd,tauxz,tauyz)
      	do i=1,np
          sflux(i)=-fluxsu(i)-fluxlu(i)-(hradu(i)-hradd(i))
      	enddo
      	if(nscreen.eq.1) write(*,*)'heat budge model completes...'
      	write(16,*)'heat budge model completes...'
      endif

!...
!...  Begin time stepping
!...
      do it=iths+1,nt

      time=it*dt 

!...  define ramp function for boundary elevation forcing, wind and pressure
!...  forcing and tidal potential forcing
!...
      if(ibc.eq.0.and.nrampbc.ne.0) then
        rampbc=tanh(2*time/86400/drampbc)
      else
      	rampbc=1
      endif

      if(nws.gt.0.and.nrampwind.ne.0) then
        rampwind=tanh(2*time/86400/drampwind)
      else
        rampwind=1
      endif

      if(nramp.eq.1) then
 	ramp=tanh(2*time/86400/dramp)
      	ramptide=ramp
      else
        ramp=1
        ramptide=1
      endif

!...  Bottom friction apart from the vel. magnitude
!	Error: shall we change tk back to no-btrack?
      do i=1,ns
      	if(ntau.eq.0) then
          Cd(i)=Cd0
      	else if(ntau.eq.1) then
          if(dps(i).le.hestu) then
            bthick=bthick_estu !thickness of boundary layer (in meters)
          else if(dps(i).le.hcont) then
            bthick=bthick_cont
          else !dps(i) > hcont
            bthick=bthick_ocea
          endif
          gthick=dmax1(dz(i,kbs(i))/2,bthick_estu) !thickness of half bottom layer
          Cd(i)=1/(2.5*dlog(gthick/1.0d-2))**2 
          Cdmin=1/(2.5*dlog(bthick/1.0d-2))**2 
          Cd(i)=dmax1(Cd(i),Cdmin)
        else !ntau=2
          Cd(i)=bdragc(i)
        endif
        tk(i)=0 !initialize for some weird sides
      enddo !i=1,ns

      if(nscreen.eq.1) write(*,*)'done Smag. and bottom fric...'
      write(16,*)'done Smag. and bottom fric...'
      
!...  process new wind info 
!...
      if(nws.eq.1) then
        if(time.gt.wtime2) then
          wtime1=wtime2
          wtime2=wtime2+wtiminc
          read(22,*) wx2,wy2
          do i=1,np
            windx1(i)=windx2(i)
            windy1(i)=windy2(i)
            windx2(i)=wx2
            windy2(i)=wy2
          enddo
        endif

        wtratio=(time-wtime1)/wtiminc
!	Compute wind at nodes for output only
        do i=1,np
          windxp(i)=windx1(i)+wtratio*(windx2(i)-windx1(i))
          windyp(i)=windy1(i)+wtratio*(windy2(i)-windy1(i))
        enddo !i

        do i=1,ns
          n1=isidenode(i,1)
          n2=isidenode(i,2)
          windx1_s=(windx1(n1)+windx1(n2))/2
          windx2_s=(windx2(n1)+windx2(n2))/2
          windy1_s=(windy1(n1)+windy1(n2))/2
          windy2_s=(windy2(n1)+windy2(n2))/2
          windx(i)=windx1_s+wtratio*(windx2_s-windx1_s)
          windy(i)=windy1_s+wtratio*(windy2_s-windy1_s)
        enddo
      endif !nws=1

!	CORIE mode
      if(nws.eq.2) then
        if(time>=wtime2) then
!...      Heat budget & wind stresses
          if(ihconsv.ne.0) then
            call surf_fluxes(wtime2,windx2,windy2,pr2,airt2,shum2,&
     &srad,fluxsu,fluxlu,hradu,hradd,tauxz,tauyz)
            do i=1,np
     	      sflux(i)=-fluxsu(i)-fluxlu(i)-(hradu(i)-hradd(i))
            enddo
            if(nscreen.eq.1) write(*,*)'heat budge model completes...'
            write(16,*)'heat budge model completes...'
          endif !ihconsv.ne.0

          wtime1=wtime2
          wtime2=wtime2+wtiminc
          do i=1,np
            windx1(i)=windx2(i)
            windy1(i)=windy2(i)
	    pr1(i)=pr2(i)
            airt1(i)=airt2(i)
            shum1(i)=shum2(i)
          enddo
          call get_wind(wtime2,windx2,windy2,pr2,airt2,shum2)

        endif !time.ge.wtime2

        wtratio=(time-wtime1)/wtiminc

!	Compute wind at nodes for output only
      	do i=1,np
          windxp(i)=windx1(i)+wtratio*(windx2(i)-windx1(i))
          windyp(i)=windy1(i)+wtratio*(windy2(i)-windy1(i))
        enddo !i

        do i=1,ns
          n1=isidenode(i,1)
          n2=isidenode(i,2)
          windx1_s=(windx1(n1)+windx1(n2))/2
          windx2_s=(windx2(n1)+windx2(n2))/2
          windy1_s=(windy1(n1)+windy1(n2))/2
          windy2_s=(windy2(n1)+windy2(n2))/2
          windx(i)=windx1_s+wtratio*(windx2_s-windx1_s)
          windy(i)=windy1_s+wtratio*(windy2_s-windy1_s)
        enddo
      endif !nws=2

!...  compute wind stress components
      dragcmin=1.0d-3*(0.61+0.063*6)
      dragcmax=1.0d-3*(0.61+0.063*50)
      do i=1,ns
        if(nws.eq.0) then
	  tau_n(i)=0
	  tau_t(i)=0
	else if(nws.eq.1.or.nws.eq.2.and.ihconsv.eq.0) then
          wmag=dsqrt(windx(i)**2+windy(i)**2)
          windn=windx(i)*snx(i)+windy(i)*sny(i)
          windt=-windx(i)*sny(i)+windy(i)*snx(i)
          dragcoef=1.0d-3*(0.61+0.063*wmag)
          dragcoef=dmin1(dmax1(dragcoef,dragcmin),dragcmax)
          tau_n(i)=dragcoef*0.001293*wmag*windn*rampwind
          tau_t(i)=dragcoef*0.001293*wmag*windt*rampwind
	else !nws=2 and ihconsv !=0; tauxz and tauyz defined
	  n1=isidenode(i,1)
	  n2=isidenode(i,2)
	  if(kfp(n1).eq.0.or.kfp(n2).eq.0) then
	    tau_n(i)=0
	    tau_t(i)=0
	  else
	    tau_n(i)=(tauxz(n1)+tauxz(n2))/2*snx(i)+(tauyz(n1)+tauyz(n2))/2*sny(i)
	    tau_t(i)=-(tauxz(n1)+tauxz(n2))/2*sny(i)+(tauyz(n1)+tauyz(n2))/2*snx(i)
	    tau_n(i)=-tau_n(i)/rho0*rampwind !sign and scale difference between stresses tauxz and tau_n
	    tau_t(i)=-tau_t(i)/rho0*rampwind
	  endif
	endif !nws
      enddo !i=1,ns


      if(nscreen.eq.1) write(*,*)'done adjusting wind stress ...'
      write(16,*)'done adjusting wind stress ...'

!...  Get new t.h. values fort.5[0-3]
!...
      if(nettype.gt.0) then
        read(50,*) ttt,(ath(i),i=1,nettype)
        if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then
          write(11,*)'Starting time wrong for eta',it,ttt
          stop
        endif
      
        icount=0
        do k=1,nope
          if(iettype(k).eq.1) then
            icount=icount+1
            if(icount.gt.nettype) then
              write(11,*)'Wrong counting 1'
              stop
            endif
            eth(k)=ath(icount)
          endif
        enddo 
      endif !nettype.gt.0

      if(nfltype.gt.0) then
        read(51,*) ttt,(ath(i),i=1,nfltype)
        if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then
          write(11,*)'Starting time wrong for flux',it,ttt
          stop
        endif

        icount=0
      	do k=1,nope
          if(ifltype(k).eq.1) then
            icount=icount+1
            if(icount.gt.nfltype) then
              write(11,*)'Wrong counting 2'
              stop
            endif
            qth(k)=ath(icount)
          endif
        enddo !k
      endif

      if(ntetype.gt.0) then
      	read(52,*) ttt,(ath(i),i=1,ntetype)
      	if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then
          write(11,*)'Starting time wrong for temp',it,ttt
          stop
      	endif
        
      	icount=0
      	do k=1,nope
          if(itetype(k).eq.1) then
            icount=icount+1
            if(icount.gt.ntetype) then
              write(11,*)'Wrong counting 3'
              stop
            endif
            tth(k)=ath(icount)
          endif
      	enddo !k
      endif

      if(nsatype.gt.0) then
      	read(53,*) ttt,(ath(i),i=1,nsatype)
      	if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then
          write(11,*)'Starting time wrong for salt',it,ttt
          stop
        endif

        icount=0
        do k=1,nope
          if(isatype(k).eq.1) then
            icount=icount+1
            if(icount.gt.nsatype) then
              write(11,*)'Wrong counting 4'
              stop
            endif
            sth(k)=ath(icount)
          endif
      	enddo !k
      endif

!...  Compute new vel. for flow b.c.
!...
      do i=1,nope
        if(ifltype(i)>=1) then
          isd1=iosd(i,1)
          n1=isidenode(isd1,1)
          n2=isidenode(isd1,2)
          H2=dps(isd1)+(peta(n1)+peta(n2))/2
          if(H2.lt.h0) then
            write(11,*)'Dry bnd side',n1,n2,H2
            stop
          endif
          uth(i)=qth(i)*ramp/cwidth(i)/H2*snx(isd1)
          vth(i)=qth(i)*ramp/cwidth(i)/H2*sny(isd1)
      	endif !ifltype(i).ge.1
      enddo !i=1,nope

      if(nscreen.eq.1) write(*,*)'done flow b.c.'
      write(16,*)'done flow b.c.'

!
!************************************************************************
!									*
!			Backtracking 					*
!									*
!************************************************************************
!

!...  Compute # of subdivisions for nsubfl=2
!...
      if(nsubfl==2) then
      	do i=1,np
          do k=1,nvrt
            sum=0
            do j=1,nne(i)
              ie=ine(i,j)
	      id=iself(i,j)
	      devm=0
	      do l=1,2
		if(l==1) then
	          nd=nm(ie,nx(i34(ie),id,1))
		else
	          nd=nm(ie,nx(i34(ie),id,i34(ie)-1))
		endif
	        rl=dsqrt((x(i)-x(nd))**2+(y(i)-y(nd))**2)
	        uvel1=(uu2(i,k)+uu2(i,k-1))/2
	        uvel2=(uu2(nd,k)+uu2(nd,k-1))/2
	        vvel1=(vv2(i,k)+vv2(i,k-1))/2
	        vvel2=(vv2(nd,k)+vv2(nd,k-1))/2
                dev=dsqrt((uvel1-uvel2)**2+(vvel1-vvel2)**2)/rl*dt
		if(dev>devm) devm=dev
	      enddo !l
              wfc=(ww2(i,k)+ww2(i,k-1))/2
              iw=2*wfc*dt/delzmin
              idelta=max0(int(devm/1.e-1),iw)
              sum=sum+dfloat(idelta)/nne(i)
            enddo !j=1,nne(i)
            ndelt(i,k)=max0(ndelt_min,min0(ndelt_max,int(sum)))
          enddo !k=1,nvrt
      	enddo !i=1,np
      endif !nsubfl=2

      st=0 !secnds(0.0) !timing the process

!...  From centers of faces 
!...
      do 451 i=1,ns
        n1=isidenode(i,1)
      	n2=isidenode(i,2)
      	if(kfe(is(i,1)).ne.0) then
          ie0=is(i,1)
      	else if(is(i,2).ne.0.and.kfe(is(i,2)).ne.0)