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