do it=1,iths            read(51,*) ttt,qq          enddo !it        endif      	if(ntetype>0) then          do it=1,iths            read(52,*) ttt,te          enddo !it      	endif      	if(nsatype>0) then          do it=1,iths            read(53,*) ttt,sal          enddo !it      	endif!...  end hot start section!...      endif !ihot.ne.0      if(nscreen.eq.1) write(*,*)'done initializing variables...'      write(16,*)'done initializing variables...'!                                                                             *!******************************************************************************!                                                                             *!			open output files				      *!                                                                             *!******************************************************************************!                                                                             *!...  write global output headers!...      if(ihot<=1) then   	ifile=1 !output file #!       Convert it to a string      	write(ifile_char,'(i12)') ifile      endif      call header      if(nscreen.eq.1) write(*,*)'done initializing outputs'      write(16,*)'done initializing outputs'!                                                                             *!******************************************************************************!                                                                             *!	    	          Time stepping 				      *!                                                                             *!******************************************************************************!                                                                             *      if(ihot.eq.0) iths=0      if(nscreen.eq.1) write(*,*)'time stepping begins...',iths+1,nt      write(16,*)'time stepping begins...',iths+1,nt!...  Compute initial levels      call levels(iths,iths,itur)      if(nscreen.eq.1) write(*,*)'done computing initial levels...'      write(16,*)'done computing initial levels...'!...  Compute total # of active faces      nfaces=0      do i=1,ns        if(kfs(i)/=0) nfaces=nfaces+kfs(i)-kbs(i)+1      enddo !i=1,ns      if(nscreen.eq.1) write(*,*)'Total # of faces=',nfaces      write(16,*)'Total # of faces=',nfaces!...  Compute nodal vel.       call nodalvel      if(nscreen.eq.1) write(*,*)'done computing initial nodal vel...'      write(16,*)'done computing initial nodal vel...'!...  Compute initial density      call eqstate      if(nscreen.eq.1) write(*,*)'done computing initial density...'      write(16,*)'done computing initial density...'!...  Store initial elevation field (hot start) for relaxation in sponge layer!...      do i=1,ne        etaic(i)=eta1(i)      enddo!...  For UB & MYG, compute vertical diffusivities at whole levels       if(itur.eq.3) then        do i=1,ns    	  do j=kbs(i),kfs(i) !wet sides	    call asm(g,i,j,vd,td,qd,qd2)            vdiff(i,j)=dmin1(diffmax(i),dmax1(bgiff,vd))            tdiff(i,j)=dmin1(diffmax(i),dmax1(bgdiff,td))            qdiff(i,j)=dmin1(diffmax(i),dmax1(bgdiff,qd))            qdiff2(i,j)=dmin1(diffmax(i),dmax1(bgdiff,qd2))          enddo !j=kbs(i),kfs(i)        enddo !i=1,ns      endif !itur=3!...  Initialize heat budget model      if(nws==2.and.ihconsv/=0) then        call surf_fluxes(wtime1,windx1,windy1,pr1,airt1,shum1,&     &srad,fluxsu,fluxlu,hradu,hradd,tauxz,tauyz)      	do i=1,np          sflux(i)=-fluxsu(i)-fluxlu(i)-(hradu(i)-hradd(i))      	enddo      	if(nscreen.eq.1) write(*,*)'heat budge model completes...'      	write(16,*)'heat budge model completes...'      endif!...!...  Begin time stepping!...      do it=iths+1,nt      time=it*dt !...  define ramp function for boundary elevation forcing, wind and pressure!...  forcing and tidal potential forcing!...      if(ibc.eq.0.and.nrampbc.ne.0) then        rampbc=tanh(2*time/86400/drampbc)      else      	rampbc=1      endif      if(nws.gt.0.and.nrampwind.ne.0) then        rampwind=tanh(2*time/86400/drampwind)      else        rampwind=1      endif      if(nramp.eq.1) then 	ramp=tanh(2*time/86400/dramp)      	ramptide=ramp      else        ramp=1        ramptide=1      endif!...  Bottom friction apart from the vel. magnitude!	Error: shall we change tk back to no-btrack?      do i=1,ns      	if(ntau.eq.0) then          Cd(i)=Cd0      	else if(ntau.eq.1) then          if(dps(i).le.hestu) then            bthick=bthick_estu !thickness of boundary layer (in meters)          else if(dps(i).le.hcont) then            bthick=bthick_cont          else !dps(i) > hcont            bthick=bthick_ocea          endif          gthick=dmax1(dz(i,kbs(i))/2,bthick_estu) !thickness of half bottom layer          Cd(i)=1/(2.5*dlog(gthick/1.0d-2))**2           Cdmin=1/(2.5*dlog(bthick/1.0d-2))**2           Cd(i)=dmax1(Cd(i),Cdmin)        else !ntau=2          Cd(i)=bdragc(i)        endif        tk(i)=0 !initialize for some weird sides      enddo !i=1,ns      if(nscreen.eq.1) write(*,*)'done Smag. and bottom fric...'      write(16,*)'done Smag. and bottom fric...'      !...  process new wind info !...      if(nws.eq.1) then        if(time.gt.wtime2) then          wtime1=wtime2          wtime2=wtime2+wtiminc          read(22,*) wx2,wy2          do i=1,np            windx1(i)=windx2(i)            windy1(i)=windy2(i)            windx2(i)=wx2            windy2(i)=wy2          enddo        endif        wtratio=(time-wtime1)/wtiminc!	Compute wind at nodes for output only        do i=1,np          windxp(i)=windx1(i)+wtratio*(windx2(i)-windx1(i))          windyp(i)=windy1(i)+wtratio*(windy2(i)-windy1(i))        enddo !i        do i=1,ns          n1=isidenode(i,1)          n2=isidenode(i,2)          windx1_s=(windx1(n1)+windx1(n2))/2          windx2_s=(windx2(n1)+windx2(n2))/2          windy1_s=(windy1(n1)+windy1(n2))/2          windy2_s=(windy2(n1)+windy2(n2))/2          windx(i)=windx1_s+wtratio*(windx2_s-windx1_s)          windy(i)=windy1_s+wtratio*(windy2_s-windy1_s)        enddo      endif !nws=1!	CORIE mode      if(nws.eq.2) then        if(time>=wtime2) then!...      Heat budget & wind stresses          if(ihconsv.ne.0) then            call surf_fluxes(wtime2,windx2,windy2,pr2,airt2,shum2,&     &srad,fluxsu,fluxlu,hradu,hradd,tauxz,tauyz)            do i=1,np     	      sflux(i)=-fluxsu(i)-fluxlu(i)-(hradu(i)-hradd(i))            enddo            if(nscreen.eq.1) write(*,*)'heat budge model completes...'            write(16,*)'heat budge model completes...'          endif !ihconsv.ne.0          wtime1=wtime2          wtime2=wtime2+wtiminc          do i=1,np            windx1(i)=windx2(i)            windy1(i)=windy2(i)	    pr1(i)=pr2(i)            airt1(i)=airt2(i)            shum1(i)=shum2(i)          enddo          call get_wind(wtime2,windx2,windy2,pr2,airt2,shum2)        endif !time.ge.wtime2        wtratio=(time-wtime1)/wtiminc!	Compute wind at nodes for output only      	do i=1,np          windxp(i)=windx1(i)+wtratio*(windx2(i)-windx1(i))          windyp(i)=windy1(i)+wtratio*(windy2(i)-windy1(i))        enddo !i        do i=1,ns          n1=isidenode(i,1)          n2=isidenode(i,2)          windx1_s=(windx1(n1)+windx1(n2))/2          windx2_s=(windx2(n1)+windx2(n2))/2          windy1_s=(windy1(n1)+windy1(n2))/2          windy2_s=(windy2(n1)+windy2(n2))/2          windx(i)=windx1_s+wtratio*(windx2_s-windx1_s)          windy(i)=windy1_s+wtratio*(windy2_s-windy1_s)        enddo      endif !nws=2!...  compute wind stress components      dragcmin=1.0d-3*(0.61+0.063*6)      dragcmax=1.0d-3*(0.61+0.063*50)      do i=1,ns        if(nws.eq.0) then	  tau_n(i)=0	  tau_t(i)=0	else if(nws.eq.1.or.nws.eq.2.and.ihconsv.eq.0) then          wmag=dsqrt(windx(i)**2+windy(i)**2)          windn=windx(i)*snx(i)+windy(i)*sny(i)          windt=-windx(i)*sny(i)+windy(i)*snx(i)          dragcoef=1.0d-3*(0.61+0.063*wmag)          dragcoef=dmin1(dmax1(dragcoef,dragcmin),dragcmax)          tau_n(i)=dragcoef*0.001293*wmag*windn*rampwind          tau_t(i)=dragcoef*0.001293*wmag*windt*rampwind	else !nws=2 and ihconsv !=0; tauxz and tauyz defined	  n1=isidenode(i,1)	  n2=isidenode(i,2)	  if(kfp(n1).eq.0.or.kfp(n2).eq.0) then	    tau_n(i)=0	    tau_t(i)=0	  else	    tau_n(i)=(tauxz(n1)+tauxz(n2))/2*snx(i)+(tauyz(n1)+tauyz(n2))/2*sny(i)	    tau_t(i)=-(tauxz(n1)+tauxz(n2))/2*sny(i)+(tauyz(n1)+tauyz(n2))/2*snx(i)	    tau_n(i)=-tau_n(i)/rho0*rampwind !sign and scale difference between stresses tauxz and tau_n	    tau_t(i)=-tau_t(i)/rho0*rampwind	  endif	endif !nws      enddo !i=1,ns      if(nscreen.eq.1) write(*,*)'done adjusting wind stress ...'      write(16,*)'done adjusting wind stress ...'!...  Get new t.h. values fort.5[0-3]!...      if(nettype.gt.0) then        read(50,*) ttt,(ath(i),i=1,nettype)        if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then          write(11,*)'Starting time wrong for eta',it,ttt          stop        endif              icount=0        do k=1,nope          if(iettype(k).eq.1) then            icount=icount+1            if(icount.gt.nettype) then              write(11,*)'Wrong counting 1'              stop            endif            eth(k)=ath(icount)          endif        enddo       endif !nettype.gt.0      if(nfltype.gt.0) then        read(51,*) ttt,(ath(i),i=1,nfltype)        if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then          write(11,*)'Starting time wrong for flux',it,ttt          stop        endif        icount=0      	do k=1,nope          if(ifltype(k).eq.1) then            icount=icount+1            if(icount.gt.nfltype) then              write(11,*)'Wrong counting 2'              stop            endif            qth(k)=ath(icount)          endif        enddo !k      endif      if(ntetype.gt.0) then      	read(52,*) ttt,(ath(i),i=1,ntetype)      	if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then          write(11,*)'Starting time wrong for temp',it,ttt          stop      	endif              	icount=0      	do k=1,nope          if(itetype(k).eq.1) then            icount=icount+1            if(icount.gt.ntetype) then              write(11,*)'Wrong counting 3'              stop            endif            tth(k)=ath(icount)          endif      	enddo !k      endif      if(nsatype.gt.0) then      	read(53,*) ttt,(ath(i),i=1,nsatype)      	if(it.eq.iths+1.and.abs(ttt-time).gt.1.e-4) then          write(11,*)'Starting time wrong for salt',it,ttt          stop        endif        icount=0        do k=1,nope          if(isatype(k).eq.1) then            icount=icount+1            if(icount.gt.nsatype) then              write(11,*)'Wrong counting 4'              stop            endif            sth(k)=ath(icount)          endif      	enddo !k      endif!...  Compute new vel. for flow b.c.!...      do i=1,nope        if(ifltype(i)>=1) then          isd1=iosd(i,1)          n1=isidenode(isd1,1)          n2=isidenode(isd1,2)          H2=dps(isd1)+(peta(n1)+peta(n2))/2          if(H2.lt.h0) then            write(11,*)'Dry bnd side',n1,n2,H2            stop          endif          uth(i)=qth(i)*ramp/cwidth(i)/H2*snx(isd1)          vth(i)=qth(i)*ramp/cwidth(i)/H2*sny(isd1)      	endif !ifltype(i).ge.1      enddo !i=1,nope      if(nscreen.eq.1) write(*,*)'done flow b.c.'      write(16,*)'done flow b.c.'!!************************************************************************!									*!			Backtracking 					*!									*!************************************************************************!!...  Compute # of subdivisions for nsubfl=2!...      if(nsubfl==2) then      	do i=1,np          do k=1,nvrt            sum=0            do j=1,nne(i)              ie=ine(i,j)	      id=iself(i,j)	      devm=0	      do l=1,2		if(l==1) then	          nd=nm(ie,nx(i34(ie),id,1))		else	          nd=nm(ie,nx(i34(ie),id,i34(ie)-1))		endif	        rl=dsqrt((x(i)-x(nd))**2+(y(i)-y(nd))**2)	        uvel1=(uu2(i,k)+uu2(i,k-1))/2	        uvel2=(uu2(nd,k)+uu2(nd,k-1))/2	        vvel1=(vv2(i,k)+vv2(i,k-1))/2	        vvel2=(vv2(nd,k)+vv2(nd,k-1))/2                dev=dsqrt((uvel1-uvel2)**2+(vvel1-vvel2)**2)/rl*dt		if(dev>devm) devm=dev	      enddo !l              wfc=(ww2(i,k)+ww2(i,k-1))/2              iw=2*wfc*dt/delzmin              idelta=max0(int(devm/1.e-1),iw)              sum=sum+dfloat(idelta)/nne(i)            enddo !j=1,nne(i)            ndelt(i,k)=max0(ndelt_min,min0(ndelt_max,int(sum)))          enddo !k=1,nvrt      	enddo !i=1,np      endif !nsubfl=2      st=0 !secnds(0.0) !timing the process!...  From centers of faces !...      do 451 i=1,ns        n1=isidenode(i,1)      	n2=isidenode(i,2)      	if(kfe(is(i,1)).ne.0) then          ie0=is(i,1)      	else if(is(i,2).ne.0.and.kfe(is(i,2)).ne.0)