elcirc5_01_01c.f90

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

!     Flow sides 
      do i=1,ns
      	isflowside(i)=0
      	isflowside2(i)=0
      enddo

      do i=1,nope
        if(ifltype(i)>=1) then 
          do j=1,noe(i)
            iel=ioe(i,j)
            do k=1,i34(iel)
              isd=js(iel,k)
              isflowside(isd)=i
              isflowside2(isd)=i
            enddo !k
      	  enddo !j
        endif
      enddo

!     Global output parameters
      noutput=20
      if(noutput.gt.mnout) then
        write(11,*)'Increase mnout in the header to',noutput
        stop
      endif
      outfile(1)='elev.61'
      outfile(2)='pres.61'
      outfile(3)='airt.61'
      outfile(4)='shum.61'
      outfile(5)='srad.61'
      outfile(6)='flsu.61'
      outfile(7)='fllu.61'
      outfile(8)='radu.61'
      outfile(9)='radd.61'
      outfile(10)='flux.61'
      outfile(11)='wind.62'
      outfile(12)='wist.62'
      outfile(13)='hvel.64'
      outfile(14)='vert.63'
      outfile(15)='temp.63'
      outfile(16)='salt.63'
      outfile(17)='conc.63'
      outfile(18)='tdff.63'
      outfile(19)='kine.63'
      outfile(20)='mixl.63'
      variable_nm(1)='surface elevation'
      variable_nm(2)='atmopheric pressure'
      variable_nm(3)='air temperature'
      variable_nm(4)='specific humidity'
      variable_nm(5)='solar radiation'
      variable_nm(6)='fluxsu'
      variable_nm(7)='fluxlu'
      variable_nm(8)='hradu'
      variable_nm(9)='hradd'
      variable_nm(10)='total flux'
      variable_nm(11)='wind speed'
      variable_nm(12)='wind stress (m^2/s^2)'
      variable_nm(13)='horizontal velocity'
      variable_nm(14)='vertical velocity'
      variable_nm(15)='temperature in C'
      variable_nm(16)='salinity in psu'
      variable_nm(17)='density in kg/m^3'
      variable_nm(18)='diffusivity for transport'
      variable_nm(19)='turbulent kinetic energy'
      variable_nm(20)='turbulent mixing length'
      do i=1,10
        variable_dim(i)='2D scalar'
      enddo
      variable_dim(11)='2D vector'
      variable_dim(12)='2D vector'
      variable_dim(13)='3D vector'
      do i=14,noutput
        variable_dim(i)='3D scalar'
      enddo
      do i=1,noutput
        if(i.le.12) then
          vpos(i)=0
      	else if(i.ge.15.and.i.le.17) then
          vpos(i)=0.5
      	else 
          vpos(i)=1
      	endif
      enddo

      read(15,*) nspool,ihfskip !output and file spools
      if(nspool.eq.0.or.ihfskip.eq.0) then
        write(11,*)'Zero nspool'
        stop
      endif
      if(mod(ihfskip,nspool).ne.0) then
	write(11,*)'ihfskip/nspool != integer'
	stop
      endif
      nrec=min(nt,ihfskip)/nspool

      do i=1,noutput
        read(15,*) iof(i)
        if(iof(i).ne.0.and.iof(i).ne.1) then
          write(11,*)'Unknown output option',i,iof(i)
          stop
        endif
      enddo !i=1,noutput

!...  Test output parameters
      read(15,*) noutgm 
      if(noutgm.ne.1.and.noutgm.ne.0) then
        write(11,*)'Unknown noutgm',noutgm
        stop
      endif
      
!...  input information about hot start output
!...
      read(15,*) nhstar
      if(nhstar.ne.0.and.nhstar.ne.1) then
        write(11,*)'Unknown nhstar',nhstar
        stop
      endif

!... Itpack solver info
      read(15,*) isolver,itmax1,iremove,zeta,tol
      if(itmax1.gt.itmax) then
        write(11,*)'Increase itmax in header file'
      	stop
      endif
      if(isolver.lt.1.or.isolver.gt.4) then
        write(11,*)'Unknown solver',isolver
        stop
      endif

!...  Compute flux flag
      read(15,*) iflux,ihcheck

!.... W switch; inactive
      read(15,*) !iwmode

!.... Mode splitting factor
      read(15,*) nsplit

!...  Check last parameter read in from fort.15
      write(*,*)'Last parameter in fort.15 is nsplit=',nsplit
      close(15)
!     End reading fort.15

      if(nscreen.eq.1) write(*,*)'done reading fort.14, 15, and 19...'
      write(16,*)'done reading fort.14, 15, and 19...'


!								   *
!*******************************************************************
!								   *
!	Initialization for cold and hot start			   *
!								   *
!*******************************************************************
!								   *

!...  Coriolis parameter
!...
      if(ncor<=0) then
      	do i=1,ns
          cori(i)=coricoef
      	enddo
      else !ncor=1
        open(31,file='coriolis.out')
      	fc=2*omega*dsin(sfea0)
      	beta=2*omega*dcos(sfea0)
      	do i=1,ns
          id1=isidenode(i,1)
          id2=isidenode(i,2)
          sphi=(ylat(id1)+ylat(id2))/2
          cori(i)=fc+beta*(sphi-sfea0)
	  if(iwrite.eq.0) then
            write(31,*)i,xcj(i),ycj(i),cori(i)
	  else !evm
	    write(31,"(a,i6,a,f16.9,a,f16.9,a,es22.14e3,a)",advance="no") &
     &           " ",i," ",xcj(i)," ",ycj(i)," ",cori(i),"\n"
	  endif
      	enddo !i=1,ns
      	close(31)
      endif


!								   *
!*******************************************************************
!								   *
!	Initialization for cold start alone			   *
!								   *
!*******************************************************************
!								   *

      if(ihot.eq.0) then
!------------------------------------------------------------------
!...  read the initial salinity and temperature values from 
!...  salinity.bp and temperature.bp files. Initial S,T fields may vary
!...  either horizontally (and vertically homogeneous) or vertically 
!...  (horizontally homogeneous). For more general 3D case, use hot start.
!...
      if(ibc.eq.1.and.ibtp.eq.0) then
!	Reset ictemp and icsalt
      	ictemp=1
      	icsalt=1
        if(10.lt.tempmin.or.10.gt.tempmax.or.33.lt.saltmin.or.33.gt.saltmax) then
          write(11,*)'Pls reset ST range to include S=33 and T=10C'
          stop
        endif
      	do i=1,np
	  do k=1,nvrt
            tem0(i,k)=10
            sal0(i,k)=33
	  enddo !k
      	enddo !i
      else !read in S,T
        open(24,file='temp.ic',status='old')
        open(25,file='salt.ic',status='old')
        if(ictemp.eq.1) then
          read(24,*) 
          read(24,*) !np
          do i=1,np
            read(24,*) num,xtmp,ytmp,te
            if(te.lt.tempmin.or.te.gt.tempmax) then
              write(11,*)'Initial invalid T at',i,te
              stop
            endif
	    do k=1,nvrt
              tem0(i,k)=te
	    enddo !k
          enddo !i
        else !ictemp=2
          read(24,*) !nvrt
          do k=1,nvrt
            read(24,*)num,te
            if(te.lt.tempmin.or.te.gt.tempmax) then
              write(11,*)'Initial invalid T at',k,te
              stop
            endif
	    do i=1,np
              tem0(i,k)=te
	    enddo !i
          enddo !k
        endif

        if(icsalt.eq.1) then
          read(25,*) 
          read(25,*) !np
          do i=1,np
            read(25,*) num,xtmp,ytmp,sa
            if(sa.lt.saltmin.or.sa.gt.saltmax) then
              write(11,*)'Initial invalid S at',i,sa
              stop
            endif
	    do k=1,nvrt
              sal0(i,k)=sa
	    enddo !k
          enddo
        else !icsalt=2
          read(25,*) !nvrt
          do k=1,nvrt
            read(25,*)num,sa
            if(sa.lt.saltmin.or.sa.gt.saltmax) then
              write(11,*)'Initial invalid S at',k,sa
              stop
            endif
	    do i=1,np
              sal0(i,k)=sa
	    enddo !i
          enddo !i
        endif
        close(24)
        close(25)
      endif !ibc.eq.1.and.ibtp.eq.0

!...  initialize S,T
      do i=1,np
        do j=1,nvrt
          tnd(i,j)=tem0(i,j)
          snd(i,j)=sal0(i,j)
      	enddo
      enddo

      do i=1,ns
	n1=isidenode(i,1)
	n2=isidenode(i,2)
        do j=1,nvrt
          tsd(i,j)=(tem0(n1,j)+tem0(n2,j))/2
          ssd(i,j)=(sal0(n1,j)+sal0(n2,j))/2
      	enddo
      enddo

!...  initialize elevations and vel. 
!...
      do i=1,ne
        eta1(i)=0
        eta2(i)=0
        do j=0,nvrt
          we(i,j)=0
        enddo
      enddo

      do i=1,np
        peta(i)=0
      	ibad(i)=0
        do k=0,nvrt
          uu1(i,k)=0.!only for internal mode
          vv1(i,k)=0
          ww1(i,k)=0
      	enddo
      enddo 

      do i=1,ns
        do j=0,nvrt
          vn2(i,j)=0 !0.1*snx(i)
          vt2(i,j)=0 !-0.1*sny(i)
        enddo
      enddo

!... initialize q2 and xl in MY-G scheme
!...
      if(itur.eq.3) then
        do i=1,ns
      	  do j=0,nvrt
            xl(i,j)=dmax1(xlmin2(i),0.1*dmax1(h0,dps(i))) 
            q2(i,j)=q2min
          enddo !j
        enddo !i=1,ns	
      endif !itur=3

!...  initialize wind for nws=1,2 (first two lines)
!...
      if(nws.eq.1) then
        open(22,file='wind.th',status='old')
        read(22,*) wx1,wy1
        read(22,*) wx2,wy2
        do i=1,np
          windx1(i)=wx1
          windy1(i)=wy1
          windx2(i)=wx2
          windy2(i)=wy2
        enddo
        wtime1=0
      	wtime2=wtiminc 
      endif

!	CORIE mode
      if(nws.eq.2) then
      	wtime1=0
      	wtime2=wtiminc 
      	call get_wind(wtime1,windx1,windy1,pr1,airt1,shum1)
      	call get_wind(wtime2,windx2,windy2,pr2,airt2,shum2)
      endif !nws=2
!------------------------------------------------------------------
      endif !ihot=0

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

!                                                                             
!******************************************************************************
!                                                                             *
!		hot start setup of the program				      *
!                                                                             *
!******************************************************************************
!
!	Record length for hot start files (double precision for all reals)
      ihot_len=nbyte*(3+4*ne+2*ne*(nvrt+1)+4*ns*(nvrt+1)+4*ns*nvrt+ &
     &	       3*np+7*np*(nvrt+1)+8*np*nvrt+1)+12
      if(itur.eq.3) ihot_len=ihot_len+nbyte*4*ns*(nvrt+1)

      if(ihot.ne.0) then
      	open(36,file='hotstart.in',access='direct',recl=ihot_len)
      	if(itur==3) then
          read(36,rec=1)time,iths,(eta1(i),eta2(i), &
     &(we(i,j),j=0,nvrt),i=1,ne), &
     &((vn2(i,j),vt2(i,j),j=0,nvrt),(tsd(i,j),ssd(i,j),j=1,nvrt),i=1,ns) &
     &,(peta(i),ibad(i),(uu1(i,j),vv1(i,j),ww1(i,j),junk,j=0,nvrt), &
     &(tnd(i,j),snd(i,j),tem0(i,j),sal0(i,j),j=1,nvrt),i=1,np), &
     &((q2(i,j),xl(i,j),j=0,nvrt),i=1,ns), &
     &ifile,ifile_char

	  do i=1,ns
            do j=0,nvrt
              q2(i,j)=dmax1(q2min,q2(i,j))
              xl(i,j)=dmax1(xlmin2(i),xl(i,j))
            enddo
          enddo
      	else !itur.ne.3
          read(36,rec=1)time,iths,(eta1(i),eta2(i), &
     &(we(i,j),j=0,nvrt),i=1,ne), &
     &((vn2(i,j),vt2(i,j),j=0,nvrt),(tsd(i,j),ssd(i,j),j=1,nvrt),i=1,ns) &
     &,(peta(i),ibad(i),(uu1(i,j),vv1(i,j),ww1(i,j),junk,j=0,nvrt), &
     &(tnd(i,j),snd(i,j),tem0(i,j),sal0(i,j),j=1,nvrt),i=1,np), &
     &ifile,ifile_char
      	endif
        close(36)

!...    change time and iteration for forecast mode
!...    Causion: this affects all t.h. files (fort.5[0-3]) and wind files
        if(ihot==1) then
          time=0
          iths=0
        endif

        write(*,*)'hot start at time=',time,iths
        write(16,*)'hot start at time=',time,iths

!...  find position in the wind input file for nws=1,2, and read in wind[x,y][1,2]
!...
        if(nws.eq.1) then
          open(22,file='wind.th',status='old')
          rewind(22)
          ninv=time/wtiminc
          wtime1=ninv*wtiminc 
          wtime2=(ninv+1)*wtiminc 
          do it=0,ninv
            read(22,*)wx1,wy1
          enddo
          read(22,*)wx2,wy2
          do i=1,np
            windx1(i)=wx1
            windy1(i)=wy1
            windx2(i)=wx2
            windy2(i)=wy2
          enddo
        endif

        if(nws.eq.2) then
          ninv=time/wtiminc
          wtime1=ninv*wtiminc 
          wtime2=(ninv+1)*wtiminc 
          call get_wind(wtime1,windx1,windy1,pr1,airt1,shum1)
          call get_wind(wtime2,windx2,windy2,pr2,airt2,shum2)
        endif !nws=2

!...  Find positions in t.h. files fort.5[0-3] 
        if(nettype>0) then
          do it=1,iths
            read(50,*) ttt,et
          enddo !it
        endif

        if(nfltype>0) then