elcirc5_01_01c.f90

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                                	!
!                                ELCIRC:                                    		!
!         A Three-Dimensional Baroclinic Model for Unstructured Grids         		!
!		        Version 5.01.01c (Sept. 05, 2003)            	      		!
!                                                                             		!
!                 Center for Coastal and Land-Margin Research                 		!
!             Department of Environmental Science and Engineering             		!
!                   OGI School of Science and Engineering,		      		!
!	              Oregon Health & Science University		      		!
!                       Beaverton, Oregon 97006, USA                            	!
!								 	      		!
!                   Scientific direction: Antonio Baptista                    		!
!                   Code development: Joseph Zhang (since April 2001).	      		!
!		                      Ed Myers (1999-May 2001);               		!
!											!
!	        Copyright 1999-2003 Oregon Health and Science University		!
!  		               All Rights Reserved					!
!                                                                             		!
!     	Formulation of Elcirc was inspired by Casulli and Zanolli (1998).       	!
!	Sparse matrix solver dsrc2c.f90 was adapted from ITPACK:			!
!	D. Young and D. Kincaid. ``The ITPACK Package for Large Sparse Linear Systems,''!
!	in Elliptic Problem Solvers, (M. Schultz, ed.), Academic Press, 		!
!	New York, 1981, pp.163-185. 							!
!    									      		!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             		!

!...  Data type consts
      module kind_par
        implicit none
        integer, parameter :: sng_kind1=4
        integer, parameter :: dbl_kind1=8
      	real(kind=dbl_kind1), parameter :: small1=1.e-6 !small non-negative number; must be identical to that in global
      end module kind_par


!...  definition of variables
!...
!
!************************************************************************
!     			mnp < mne < mns					*
!************************************************************************
!
      module global
      	implicit none
	integer, parameter :: sng_kind=4
        integer, parameter :: dbl_kind=8

!...  	Dimensioning parameters
        integer, parameter :: mnp=38000
      	integer, parameter :: mne=56000
      	integer, parameter :: mns=85000
      	integer, parameter :: mnv=62
      	integer, parameter :: mnei=30 !neighbor
      	integer, parameter :: mnope=6 !# of open bnd segements
       	integer, parameter :: mnond=20000 !max. # of open-bnd nodes on each segment
      	integer, parameter :: mnland=50 !# of land bnd segements
      	integer, parameter :: mnlnd=10000 !max. # of land nodes on each segment
      	integer, parameter :: mnoe=20000 !max. # of open-bnd elements on each segment
      	integer, parameter :: mnosd=20000 !max # of open-bnd sides on each segment
      	integer, parameter :: mnbfr=15 !# of forcing freqs.
      	integer, parameter :: itmax=5000 !# of iteration for itpack solvers used for dimensioning
      	integer, parameter :: nwksp=6*mne+4*itmax !available work space for itpack solvers
      	integer, parameter :: nbyte=4
      	integer, parameter :: mnout=100 !max. # of output files
      	integer, parameter :: mirec=1109000000 !997000000) !max. record # to prevent output ~> 4GB
      	real(kind=dbl_kind), parameter :: small1=1.e-6 !small non-negative number; must be identical to that in kind_par

!...  	Important variables
      	integer :: np,ne,ns,nvrt
      	real(kind=dbl_kind) :: zmsl,h0,denom,q2min,rho0,dt
!	Consts. used in Canuto's ASM
        character(len=2) :: mid
        real(kind=dbl_kind) :: ubd0,ubd1,ubd2,ubd3,ubd4,ubd5,ubs0,ubs1,ubs2,ubs4,ubs5,ubs6, &
     &a2_cm03,schk,schpsi

!...    Output handles
        character(len=48) :: start_time,version,data_format='DataFormat v3.0'
        character(len=12) :: ifile_char
        character(len=48), dimension(mnout) :: outfile,variable_nm,variable_dim
        integer :: ihot,nrec,nspool,igmp,noutgm,ifile,noutput
        integer, dimension(mnout) :: ichan,irec,iof
        real(kind=dbl_kind), dimension(mnout) :: vpos
        
! evm   character buffers for binary type conversions
	integer :: iwrite
        character(len=48)   :: a_48
        character(len=16)   :: a_16
        character(len=8)    :: a_8
        character(len=4)    :: a_4


!...  1D arrays
	integer :: nne(mnp),nnp(mnp),kfs(mns),kbs(mns),kbp(mnp),kfp(mnp), &
	   &kfe(mne),kbe(mne)

	real(kind=dbl_kind), dimension(mnp) :: x,y,dp,peta
	real(kind=dbl_kind), dimension(mne) :: areas,radiel,xctr,yctr,dpe,eta1,eta2
	real(kind=dbl_kind), dimension(mns) :: snx,sny,distj,xcj,ycj,delj,dps,xlmin1,xlmin2

	real(kind=dbl_kind) :: delz(mnv),ztot(0:mnv),tem0(mnp,mnv),sal0(mnp,mnv)

!...  2D arrays
	integer :: i34(mne),nm(mne,4),nx(4,4,3),ic3(mne,4),ine(mnp,mnei),js(mne,4),is(mns,2), &
	   &isidenode(mns,2),inp(mnp,mnei),iself(mnp,mnei)

	real(kind=dbl_kind) :: ssign(mne,4),dz(mns,mnv),dzhalf(mns,0:mnv),dc(mne,mnv), &
	   &dchalf(mne,mnv),dzp(mnp,mnv),dzphalf(mnp,mnv),we(mne,0:mnv), &
	   &vn2(mns,0:mnv),vt2(mns,0:mnv),tsd(mns,mnv),ssd(mns,mnv),tnd(mnp,mnv),snd(mnp,mnv), &
	   &srho(mns,mnv),erho(mne,mnv),prho(mnp,mnv),q2(mns,0:mnv),xl(mns,0:mnv)

	real(kind=dbl_kind), dimension(mnp,0:mnv) :: uu1,vv1,ww1, uu2,vv2,ww2

      end module global

!...  Main program
      program elcirc
      use global
      implicit real(kind=dbl_kind)(a-h,o-z),integer(i-n)

!...  Output handles
      real(kind=sng_kind) :: floatout,floatout2,st,en
      character(len=12) :: it_char
      character(len=40) :: date,timestamp

!... Geometry
      dimension xlon(mnp),ylat(mnp),x1j(mns),y1j(mns),x2j(mns),y2j(mns)
      dimension xlon_e(mne),ylat_e(mne),dl(2,3) 

!...  Boundary forcings
      character(len=5) :: bountag(mnbfr)
      character(len=10) :: fq_nm(mnbfr)
      dimension isbnd(mnp),sarea(mnp),ibad(mnp),idrynode(mnp)
      dimension isbnode(mnp,2),isbe(mne),isbs(mns)
      dimension nond(mnope),iond(mnope,mnond)
      dimension nlnd(mnland),ilnd(mnland,mnlnd)
      dimension nosd(mnope),iosd(mnope,mnosd),noe(mnope),ioe(mnope,mnoe)
      dimension isflowside(mns),isflowside2(mns),cwidth(mnope)
      dimension iettype(mnope),ifltype(mnope)
      dimension itetype(mnope),isatype(mnope)
      dimension tamp(mnbfr),tnf(mnbfr),tfreq(mnbfr),jspc(mnbfr),tear(mnbfr)
      dimension fun_lat(mnp,0:2),fun_lat_e(mne,0:2)
      dimension amig(mnbfr),ff(mnbfr),face(mnbfr)
      dimension emo(mnope,mnbfr,mnoe),efa(mnope,mnbfr,mnoe)
      dimension vmo(mnope,mnbfr,mnosd),vfa(mnope,mnbfr,mnosd)
      dimension eth(mnope),qth(mnope),tth(mnope),sth(mnope)
      dimension uth(mnope),vth(mnope),ath(mnope)

!...  Flow arrays
      dimension tsdbt(mns,mnv),ssdbt(mns,mnv),vnbt(mns,mnv)
      dimension tndbt(mnp,mnv),sndbt(mnp,mnv),vtbt(mns,mnv)
      dimension zhat1(mns,mnv),ghat(mns,mnv),gvec(mnv)
      dimension zhat2(mns,mnv),fhat(mns,mnv),fvec(mnv),ndelt(mnp,mnv)
      dimension icoef(mne+1),jcoef(mne*5),e2coef(mne*5),qel(mne)
      dimension sparsem(mne,0:4),elbc(mne),imape(mne),qel2(mne),eta3(mne)
      dimension windx1(mnp),windy1(mnp),windx2(mnp),windy2(mnp)
      dimension windx(mns),windy(mns),tau_n(mns),tau_t(mns),advt(mns)
      dimension pr1(mnp),airt1(mnp),shum1(mnp),pr2(mnp),airt2(mnp),shum2(mnp)
      dimension sflux(mnp),srad(mnp)
      dimension fluxsu(mnp),fluxlu(mnp),hradu(mnp),hradd(mnp)
      dimension tauxz(mnp),tauyz(mnp),windxp(mnp),windyp(mnp)
      dimension tk(mns),cori(mns),bdragc(mns),Cd(mns)
      dimension vdiff(mns,mnv),tdiff(mns,mnv)
      dimension tdiffp(mnp,mnv),etaic(mne),relax(mnp)
      dimension rich(mns,mnv)
      dimension hor_t(mne),hor_s(mne)

!...  Wild-card arrays
      dimension nwild(mnp),nwild2(mne),swild(mnp) !,slope(mne)

!... Solver arrays for TRIDAG
      dimension alow(mnv),bdia(mnv),cupp(mnv),rrhs(mnv,5)
      dimension soln(mnv,5),gam(mnv)
      dimension trhs(mns,mnv),srhs(mns,mnv)

!     MY-G turbulence closure arrays
      dimension qdiff(mns,mnv),qdiff2(mns,mnv),q2tmp(mnv),xltmp(mnv)
      dimension q2bt(mns,mnv),xlbt(mns,mnv)
      dimension rzbt(mns,mnv),shearbt(mns,mnv),xlmax(mnv),diffmax(mns)
      dimension q2nd(mnp,mnv),xlnd(mnp,mnv)
!     Test arrays
      dimension testa(mns)

!...  variables used by the itpack solvers
      dimension iwksp(3*mne),wksp(nwksp),iparm(12),rparm(12)


!...
!...  First executible statement of Elcirc
!...

!... Initialize arrays
      idrynode=0 !not dry
      isbnd=0
      isbe=0
      isbs=0
      pr1=0; pr2=0 !uniform pressure (the const. is unimportant)
!	for output
      airt1=0; shum1=0;  airt2=0; shum2=0; srad=0; fluxsu=0; fluxlu=0 
      hradu=0; hradd=0; sflux=0; windxp=0; windyp=0

!...  define some constants and initial values
!...
      omega=7.29d-5 !angular freq. of earth rotation 
      rearth=6378206.4 !earth radius
      g=9.81
      rho0=1025. !ref. density for S=33 and T=10C
      pi=3.141592653589793d0
      deg2rad=pi/180
      rad2deg=180/pi
      shw=4184  !specific heat of pure water

      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


!...  Boundary layer characteristics
      hestu=60. !threshold depth for estuary (m)
      hcont=100.  !threshold depth for continuental shelf
      bthick_estu=0.25 !char. BL thickness for estuary
      bthick_cont=1.00 !char. BL thickness for continuental shelf
      bthick_ocea=29.8 !char. BL thickness for open ocean
      if(hestu.gt.hcont) then
        write(*,*)'Check depth scales',hestu,hcont
      stop
      endif

!...  set the maximum number of digits of precision the grid can be expected to have
!...  note: if the grid was build on a 32 bit computer, it should be accurate to about
!...  7 digits.  Thus, if the nodal spacing requires more than 5 digits of precision,
!...  it is unlikely that the model results will be trustworthy.

      nprec=5



!                                                                             *
!******************************************************************************
!                                                                             *
!			open input files				      *
!                                                                             *
!******************************************************************************
!                                                                             *

      open(14,file='hgrid.gr3',status='old')
      open(15,file='param.in',status='old')
      open(19,file='vgrid.in',status='old')

      open(11,file='fort.11') !fatal error message output
      open(12,file='fort.12') !non-fatal error message output
      open(16,file='mirror.out')
      rewind(11)
      rewind(16)

      call date_and_time(date,timestamp)
      write(16,*)'Run begins at ',date,timestamp


!...  read the vertical layers information from vgrid.in
!...
      read(19,*) nvrt,zmsl
      if(nvrt.gt.mnv) then
 	write(11,*)'nvrt > mnv'
      	stop
      endif

      ztot(0)=0 !0th level
      delzmin=1.0d15 !min. delz(i)
      do i=1,nvrt
        read(19,*) jki,delz(i),ztot(i)
      	if(delz(i).lt.delzmin) delzmin=delz(i)
      	if(dabs(ztot(i)-zmsl).lt.1.0d-5) then
          write(11,*)'zmsl is too close to level line',i
          stop
      	endif
      enddo
      close(19)

!     Check delz and ztot
      do i=1,nvrt
      	if(dabs(ztot(i-1)+delz(i)-ztot(i)).gt.1.e-6) then
          write(11,*)'Wrong input fort.19 at level',i,ztot(i-1),delz(i)
          stop
     	endif
      enddo
      if(zmsl.gt.ztot(nvrt).or.zmsl.lt.0) then
      	write(11,*)'MSL above/below  all levels'
      	stop
      endif

!...  read unit 15 file
!...

      read(15,'(a48)') version
      read(15,'(a48)') start_time
      read(15,*) ipre !pre-processing flag (to output obe.out)
      read(15,*) nscreen
      read(15,*) iwrite !write mode

      read(15,*) ihot
      if(ihot.lt.0.or.ihot.gt.2) then
      	write(11,*)'Unknown ihot',ihot
      	stop
      endif

      read(15,*) ics
      if(ics.ne.1.and.ics.ne.2) then
        write(11,*)'Unknown ics',ics
      	stop
      endif

!...  Center of projection in degrees
      read(15,*) slam0,sfea0
      slam0=slam0*deg2rad
      sfea0=sfea0*deg2rad

!... Enough info to read unit 14 grid file
!...
      read(14,*) 
      read(14,*) ne,np
      if(ne.gt.mne.or.np.gt.mnp) then
        write(11,*)'Increase mne/mnp',mne,mnp,ne,np
        stop
      endif

      do i=1,np
        if(ics.eq.1) then
          read(14,*) j,x(i),y(i),dp(i)
        else if(ics.eq.2) then
          read(14,*) j,xlon(i),ylat(i),dp(i)
          ylat(i)=ylat(i)*deg2rad
          xlon(i)=xlon(i)*deg2rad
          call cpp(x(i),y(i),xlon(i),ylat(i),slam0,sfea0)
        endif
	if(dp(i).le.0) idrynode(i)=1
      enddo !i=1,np

      do i=1,ne
        read(14,*) j,i34(i),(nm(i,k),k=1,i34(i))
	if(i34(i)/=3.and.i34(i)/=4) then
	  write(*,*)'Unknown type of element',i
	  stop
	endif

	n1=nm(i,1)
	n2=nm(i,2)
	n3=nm(i,3)
	if(i34(i)==3) then
          areas(i)=signa(x(n1),x(n2),x(n3),y(n1),y(n2),y(n3))
	else !quad
	  n4=nm(i,4)
!         Check convexity
          ar1=signa(x(n1),x(n2),x(n3),y(n1),y(n2),y(n3))
          ar2=signa(x(n1),x(n3),x(n4),y(n1),y(n3),y(n4))
          ar3=signa(x(n1),x(n2),x(n4),y(n1),y(n2),y(n4))
          ar4=signa(x(n2),x(n3),x(n4),y(n2),y(n3),y(n4))
          if(ar1<=0.or.ar2<=0.or.ar3<=0.or.ar4<=0) then
            write(*,*)'Concave quadrangle',i,ar1,ar2,ar3,ar4
            stop
          endif

          areas(i)=ar1+ar2
	endif
        if(areas(i).le.0.0) then
          write(11,*)'Negative area at',i
          stop
        endif
        radiel(i)=dsqrt(areas(i)/pi)  !equivalent radius
      enddo !i=1,ne

!     Open bnds
      read(14,*) nope
      if(nope.gt.mnope) then
      	write(11,*) 'nope > mnope' 
        stop
      endif

      read(14,*) neta
      jnmm=0
      do k=1,nope
        read(14,*) nond(k)
        if(nond(k).gt.mnond) then
  	  write(11,*) 'nond(k) > mnond'
          stop
        endif
        do i=1,nond(k)
          read(14,*) iond(k,i)
          isbnd(iond(k,i))=k
        enddo
        jnmm=jnmm+nond(k)
      enddo

      if(neta.ne.jnmm) then
        write(11,*)'neta /= total # of open bnd nodes',neta,jnmm
        stop
      endif

!     Land bnds
      read(14,*) nland
      if(nland.gt.mnland) then
      	write(11,*) 'nland > mnland'
        stop
      endif

      read(14,*) nvel
      do k=1,nland
        read(14,*) nlnd(k)
        if(nlnd(k).gt.mnlnd) then
          write(11,*)'nlnd(k) > mnlnd',nlnd(k),mnlnd
          stop
        endif
        do i=1,nlnd(k)
          read(14,*) ilnd(k,i)
          if(isbnd(ilnd(k,i)).eq.0) isbnd(ilnd(k,i))=-1 !overlap of open bnd
        enddo
      enddo !k=1,nland
      close(14)
!...  End fort.14

!                                                                             *
!                                                                             *
!******************************************************************************
!                                                                             *
!     			Compute geometry 				      *
!                                                                             *
!******************************************************************************
!                                                                             *
!                                                                             *