!     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