!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                                                                                	!!                                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 				      *!                                                                             *!******************************************************************************!                                                                             *!                                                                             *