elcirc5_01_01c.f90

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

      if(nadv.eq.0) then
      	open(42,file='adv.gr3',status='old')
        read(42,*)
        read(42,*) !ne,np
        do i=1,np
          read(42,*)j,xtmp,ytmp,tmp
          nwild(i)=tmp
      	enddo
        do i=1,ns
          n1=isidenode(i,1)
          n2=isidenode(i,2)
          advt(i)=max0(nwild(n1),nwild(n2))
        enddo
        close(42)
      endif

!...  Minimum depth allowed
      read(15,*) h0 
      if(h0.le.0) then
        write(11,*)'h0 must be positive'
        stop
      endif
      denom=dmin1(h0*(0.5-1.0e-4),delzmin*(0.5-1.e-4)) !min. denominator for [a,s,t]mat 

!...  Bottom friction
      read(15,*) ntau
      if(ntau.lt.0.or.ntau.gt.2) then
      	write(11,*)'Unknown ntau', ntau
      	stop
      endif
      if(ntau.eq.0) then
        read(15,*) Cd0
      else if(ntau.eq.2) then
!	drag.bp is a build pt file because gredit cannot keep enough precision
!	for gredit depth
      	open(21,file='drag.bp',status='old')
      	read(21,*)
      	read(21,*) npbp
      	if(npbp.ne.np) then
          write(11,*)'drag.bp is hgrid.gr3 based'
          stop
      	endif
        do i=1,np
          read(21,*)j,xtmp,ytmp,swild(i)
      	enddo
      	do i=1,ns
          n1=isidenode(i,1)
          n2=isidenode(i,2)
          bdragc(i)=(swild(n1)+swild(n2))/2
      	enddo
      	close(21)
      endif

!     Coriolis
      read(15,*) ncor
      if(abs(ncor)>1) then
        write(11,*)'Unknown ncor',ncor
        stop
      endif
      if(ncor==-1) then
        read(15,*) tmp
 	coricoef=2*omega*dsin(tmp/180*pi)
      else if(ncor==0) then
	read(15,*) coricoef
      else !ncor=1
      	write(*,*)'Check slam0 and sfea0 as variable Coriolis is used'
      	write(16,*)'Check slam0 and sfea0 as variable Coriolis is used'
      	open(32,file='hgrid.ll',status='old')
      	read(32,*)
      	read(32,*) !ne,np
      	do i=1,np
          read(32,*)j,xlon(i),ylat(i)
          xlon(i)=xlon(i)*deg2rad
          ylat(i)=ylat(i)*deg2rad
      	enddo !i
      	close(32)
      endif

!     Wind
      read(15,*) nws,wtiminc
      if(nws.lt.0.or.nws.gt.2) then
        write(11,*)'Unknown nws',nws
        stop
      endif
      if(nws.gt.0.and.dt.gt.wtiminc) then
      	write(11,*)'wtiminc < dt'
      	stop
      endif

      if(nws.gt.0) read(15,*) nrampwind,drampwind

!	Heat and salt conservation flags
      read(15,*) ihconsv,isconsv
      if(ihconsv.lt.0.or.ihconsv.gt.1.or.isconsv.lt.0.or.isconsv.gt.1) then
        write(11,*)'Unknown ihconsv or isconsv',ihconsv,isconsv
      	stop
      endif
      if(ihconsv+isconsv.ne.0.and.nws.ne.2) then
        write(11,*)'Heat/salt budge model must have nws=2'
      	stop
      endif
      if(ihconsv.ne.0) then
        write(16,*)'Warning: you have chosen a heat conservation model'
      	write(16,*)'which assumes start time at 0:00 PST!'
      endif

!...  Turbulence closure options
      read(15,*) itur
      if(itur.lt.-2.or.itur.gt.3) then
      	write(11,*)'Unknown turbulence closure model',itur
      	stop
      endif
      if(ibc.eq.1.and.ibtp.eq.0.and.itur.gt.0) then
        write(11,*)'Barotropic model cannot have turbulence closure'
        stop
      endif

      if(itur.eq.0) then
        read(15,*) vdiffcon,tdiffcon
        do j=1,nvrt
          do i=1,ns
            vdiff(i,j)=vdiffcon
            tdiff(i,j)=tdiffcon
          enddo
        enddo !j=1,nvrt
      else if(itur.eq.-1) then !VVD
        open(43,file='vvd.dat',status='old')
      	read(43,*) !nvrt
      	do j=1,nvrt
          read(43,*)k,vdiffcon,tdiffcon
          do i=1,mns
            vdiff(i,j)=vdiffcon
            tdiff(i,j)=tdiffcon
          enddo
      	enddo !j=1,nvrt
        close(43)
      else if(itur.eq.-2) then !HVD
        open(43,file='hvd.mom',status='old')
        open(44,file='hvd.tran',status='old')
      	read(43,*)
      	read(43,*) !nsbp
      	read(44,*)
      	read(44,*) !nsbp
      	do i=1,ns
          read(43,*)k,xtmp,ytmp,vdiffcon
          read(44,*)k,xtmp,ytmp,tdiffcon
          do j=1,nvrt
            vdiff(i,j)=vdiffcon
            tdiff(i,j)=tdiffcon
          enddo
        enddo !i=1,ns
      	close(43)
      	close(44)
      else if(itur.eq.2) then !read in P&P coefficients
        read(15,*) tdmin_pp,h1_pp,vdmax_pp1,vdmin_pp1,h2_pp,vdmax_pp2,vdmin_pp2
	if(h1_pp.ge.h2_pp) then
	  write(11,*)'h1_pp >= h2_pp in P&P'
	  stop
	endif
	if(vdmax_pp1.lt.vdmin_pp1.or.vdmax_pp2.lt.vdmin_pp2) then
	  write(11,*)'Wrong limits in P&P:',vdmax_pp1,vdmin_pp1,vdmax_pp2,vdmin_pp2
	  stop
	endif
      else if(itur.eq.3) then !read in const. (cf. Umlauf and Burchard 2003)
	read(15,*) mid
        read(15,*) bgdiff,hestu_my,diffmax_est,xlmin_est,hcont_my,diffmax_sea,xlmin_sea
      	if(hestu_my.gt.hcont_my) then
          write(11,*)'hestu_my > hcont_my'
          stop
        endif
	if(mid.eq.'UB') then
	  read(15,*) rmub,rnub,cpsi2,cpsi3,schk,schpsi !m,n,c_psi[2-3], and Schmidt numbers
	  if(rnub==0.or.schk<=0.or.schpsi<=0) then
	    write(11,*)'Wrong input for rnub etc:',rnub,schk,schpsi
	    stop
	  endif
	  cpsi1=rmub
!	  Consts. used in Canuto's ASM
   	  ubl1=0.1070
   	  ubl2=0.0032
   	  ubl3=0.0864
   	  ubl4=0.12
   	  ubl5=11.9
   	  ubl6=0.4
   	  ubl7=0
   	  ubl8=0.48
   	  ubs0=1.5*ubl1*ubl5**2
   	  ubs1=-ubl4*(ubl6+ubl7)+2*ubl4*ubl5*(ubl1-ubl2/3-ubl3)+1.5*ubl1*ubl5*ubl8
   	  ubs2=-0.375*ubl1*(ubl6**2-ubl7**2)
   	  ubs4=2*ubl5
   	  ubs5=2*ubl4
   	  ubs6=2*ubl5/3*(3*ubl3**2-ubl2**2)-0.5*ubl5*ubl1*(3*ubl3-ubl2)+0.75*ubl1*(ubl6-ubl7)
   	  ubd0=3*ubl5**2
   	  ubd1=ubl5*(7*ubl4+3*ubl8)
   	  ubd2=ubl5**2*(3*ubl3**2-ubl2**2)-0.75*(ubl6**2-ubl7**2)
   	  ubd3=ubl4*(4*ubl4+3*ubl8)
   	  ubd4=ubl4*(ubl2*ubl6-3*ubl3*ubl7-ubl5*(ubl2**2-ubl3**2))+ubl5*ubl8*(3*ubl3**2-ubl2**2)
   	  ubd5=0.25*(ubl2**2-3*ubl3**2)*(ubl6**2-ubl7**2)

!	  print*, 'ubd2=',ubd2,',ubd4=',ubd4,',ubd2/ubd4=',ubd2/ubd4

	else if(mid.eq.'MY') then
	  rmub=1; rnub=1; cpsi1=0.9; cpsi3=0.9 !others not used
	else
	  write(11,*)'Unknown closure model id',mid
	  stop
	endif
	cmiu0=dsqrt(0.3d0)
	a2_cm03=2/cmiu0**3
        q2min=1.e-9/2 !min. kinetic energy
        if(diffmax_est<bgdiff.or.diffmax_sea<bgdiff) then
          write(11,*)'Bg. diffusivity > diffmax'
          stop
        endif

        do i=1,ns
	  xlmin2(i)=2*q2min*0.1*dmax1(h0,dps(i)) !floor for non-surface layers
          if(dps(i)<=hestu_my) then
            xlmin1(i)=xlmin_est
	    diffmax(i)=diffmax_est
          else if(dps(i)<=hcont_my) then
            xlmin1(i)=xlmin_est+(xlmin_sea-xlmin_est)*(dps(i)-hestu_my)/(hcont_my-hestu_my)
	    diffmax(i)=diffmax_est+(diffmax_sea-diffmax_est)*(dps(i)-hestu_my)/(hcont_my-hestu_my)
          else !dps > hcont_my
            xlmin1(i)=xlmin_sea
	    diffmax(i)=diffmax_sea
          endif
          xlmin1(i)=dmax1(xlmin1(i),xlmin2(i))
        enddo !i

      endif !itur

!...  HORCON parameter
      read(15,*) ihorcon,smagcoef

!      if(ihorcon.eq.0) then
!	read(15,*) hor 
!	do i=1,mns
!	  smagcoef(i)=hor
!	enddo
!      else !ihorcon=1
!	open(99,file='fort.99',status='old')
!	read(99,*)
!	read(99,*) nsbp
!	do i=1,nsbp
!	  read(99,*)j,xtmp,ytmp,hor
!	  smagcoef(i)=hor
!	enddo
!	close(99)
!      endif

      read(15,*) ihorcon_st
      if(ihorcon_st.eq.0) then
        read(15,*) hor_t0,hor_s0 
        do i=1,ne
          hor_t(i)=hor_t0
          hor_s(i)=hor_s0
        enddo
      else !ihorcon_st=1
      	open(97,file='fort.99_T',status='old')
      	open(98,file='fort.99_S',status='old')
      	read(97,*)
      	read(97,*) !ne
      	read(98,*)
      	read(98,*) !ne
      	do i=1,ne
          read(97,*)j,xtmp,ytmp,hor_t0
          read(98,*)j,xtmp,ytmp,hor_s0
          hor_t(i)=hor_t0
          hor_s(i)=hor_s0
        enddo
        close(97)
        close(98)
      endif

      read(15,*)ictemp,icsalt
      if(ictemp.ne.1.and.ictemp.ne.2.or.icsalt.ne.1.and.icsalt.ne.2) then
      	write(11,*)'Unknown i.c. flag',ictemp,icsalt
      	stop
      endif

!...  Sponge layer
      read(15,*) isponge
      if(isponge.ne.0) then
        open(96,file='sponge.gr3',status='old')
      	read(96,*)
      	read(96,*) !ne,np
      	do i=1,np
          read(96,*)j,xtmp,ytmp,relax(i)
      	enddo !i
      	close(96)
      endif

!...  Earth tidal potential
      read(15,*) ntip,tip_dp !cut-off depth for applying tidal potential
      if(ntip.gt.mnbfr) then
        write(11,*)'ntip > mnbfr',ntip,mnbfr
        stop
      endif
      
      if(ntip.gt.0) then
	open(32,file='hgrid.ll',status='old')
        read(32,*)
        read(32,*) !ne,np
        do i=1,np
          read(32,*)j,xlon(i),ylat(i)
          xlon(i)=xlon(i)*deg2rad
          ylat(i)=ylat(i)*deg2rad
!...      Pre-compute species function to save time
	  fun_lat(i,0)=3*dsin(ylat(i))**2-1
	  fun_lat(i,1)=dsin(2*ylat(i))
	  fun_lat(i,2)=dcos(ylat(i))**2
        enddo !i
        close(32)

	do i=1,ne
	  xlon_e(i)=0
	  ylat_e(i)=0
	  do j=1,i34(i)
	    xlon_e(i)=xlon_e(i)+xlon(nm(i,j))/i34(i)
	    ylat_e(i)=ylat_e(i)+ylat(nm(i,j))/i34(i)
	  enddo !j
	  fun_lat_e(i,0)=3*dsin(ylat_e(i))**2-1
	  fun_lat_e(i,1)=dsin(2*ylat_e(i))
	  fun_lat_e(i,2)=dcos(ylat_e(i))**2
	enddo !i
      endif !ntip.gt.0
	
      do i=1,ntip
        read(15,*) !tag
	read(15,*) jspc(i),tamp(i),tfreq(i),tnf(i),tear(i)
	if(jspc(i).lt.0.or.jspc(i).gt.2) then
	  write(11,*)'Illegal tidal species #',jspc(i)
	  stop
	endif
        tear(i)=tear(i)*deg2rad
      enddo !i

!...  Boundary forcing freqs.
      read(15,*) nbfr
      if(nbfr.gt.mnbfr) then
        write(11,*)'nbfr > mnbfr',nbfr,mnbfr
        stop
      endif

      do i=1,nbfr
        read(15,'(a5)') bountag(i)
        read(15,*) amig(i),ff(i),face(i) !freq., nodal factor and earth equil.
        face(i)=face(i)*deg2rad
      enddo

      read(15,*) nope1
      if(nope1.ne.nope) then
        write(11,*)'Inconsistent # of open bnds',nope1,nope
      	stop
      endif

      nettype=0 !total # of type I eta bnds
      nfltype=0
      ntetype=0
      nsatype=0
      do k=1,nope
        read(15,*) ntmp,iettype(k),ifltype(k),itetype(k),isatype(k)
      	if(ntmp.ne.noe(k)) then
          write(11,*)'Inconsistent # of elements at open boundary',k
          write(11,*)ntmp,noe(k)
          stop
        endif

        if(iettype(k).eq.1) then
          nettype=nettype+1
!	  Mock reading
          open(50,file='elev.th',status='old')
          do j=1,nt
            read(50,*) ttt,et
          enddo !j
          rewind(50)
        else if(iettype(k).eq.2) then
          read(15,*) eth(k)
        else if(iettype(k).eq.3) then
          do i=1,nbfr
            read(15,'(a10)') fq_nm(i) !freq. name
            do j=1,noe(k)
              read(15,*) emo(k,i,j),efa(k,i,j) !amp. and phase
              efa(k,i,j)=efa(k,i,j)*deg2rad
            enddo
          enddo
        else if(iettype(k).ne.-1) then
          write(11,*)'INVALID VALUE FOR IETTYPE'
          stop
        endif

        if(ifltype(k).eq.1) then
          nfltype=nfltype+1
          open(51,file='flux.th',status='old')
          do j=1,nt
            read(51,*) ttt,qq
          enddo
          rewind(51)
      	else if(ifltype(k).eq.2) then
          read(15,*) qth(k)
        else if(ifltype(k).ne.-1.and.ifltype(k).ne.0) then
          write(11,*) 'INVALID VALUE FOR IFLTYPE'
          stop
        endif

        if(itetype(k).eq.1) then
          ntetype=ntetype+1
          open(52,file='temp.th',status='old')
          do j=1,nt
            read(52,*) ttt,temp
          enddo
          rewind(52)
        else if(itetype(k).eq.2) then
          read(15,*) tth(k)
      	else if(itetype(k).eq.-1) then
          read(15,*) tth(k)
        else if(itetype(k).lt.-1.or.itetype(k).gt.3) then
          write(11,*) 'INVALID VALUE FOR ITETYPE'
          stop
        endif

        if(isatype(k).eq.1) then
          nsatype=nsatype+1
          open(53,file='salt.th',status='old')
          do j=1,nt
            read(53,*) ttt,sal
          enddo
          rewind(53)
        else if(isatype(k).eq.2) then
          read(15,*) sth(k)
        else if(isatype(k).eq.-1) then
          read(15,*) sth(k)
        else if(isatype(k).lt.-1.or.isatype(k).gt.3) then
          write(11,*) 'INVALID VALUE FOR ISATYPE'
          stop
        endif
      enddo !k=1,nope