c    @(#)clinic.f	1.6
      program clinic
c
c --- this is a package of routines for a finite-element solution of
c --- the quasi-geostrophic equations,as described in "a baroclinic
c --- quasi-geostrophic open ocean model" (1981) by r.n. miller,
c --- a.r. robinson, and d.b. haidvogel of harvard.
c
c --- bottom topography (for slopes less than o(r0)) and top vertical
c --- velocity have been implemented (10/15/82).
c
c --- a driving package is required to provide initial and boundary
c --- values of the streamfunction, vorticity, and top and bottom
c --- density variations; as well as the topography. this package
c --- is called infld.
c
c --- a variety of options are available:
c
c        ifdiff = 1/0     finite difference/collocation in depth
c        ifpert = 1/0    calculate perturbation fields (yes/no)
c        ifrst = 1/0     restart at time rdt (yes/no)
c        iftop = 1/0     use top density information (yes/no)
c        ifbot = 1/0     use bottom density information (yes/no)
c        iftvv = 1/0     use top vertical velocity (yes/no)
c        ifrel = 1/0     use bottom relief (yes/no)
c        ifsbl = 1/0     use surface boundary layer
c        ifeva = 1/0     execute EVA preprocessor
c        ifsdo = 1/0     write the data for small domain 
c
      include 'qg_par.inc'
c
      character*4 restid,runid
      character*70 idstring

c
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva
      common/solver/vam(m0),vbm(m0),van(n0),vbn(n0)
      common/edges/icorn(4),ncorn(4),mast(4),nast(4)
      common/philtr/iltord,iltfrq,iltcnt
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)
      common/vert/rho0,tsurf,xc(k01),salsur,eta(k01)
      common/user/pram(16)
      common/dynam/alpha,beta,rf
      common/unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
      common/cnn/n,nm1,nm2,np1
      common/extpar/rlat0,rlng0,t0,v0,dhor,ht,time0,
     +                r0,f0star
      common/inform/ifdiff,ifpert,titl(20,2),runid(6),restid(9),
     +              idstring,ident
      common/control/ddt(4),adt(4),nalev(4),ialev(4,k02),iprnt(k02),
     +      lprnt(k0)
      common/rstinfo/ifrst,fmlrrs,rdt,imlcnt
      common/asspar/isd,ibeg,iasday,iastep,ifsdo,imf,icit
      common/assfield/ch3(mn0),ch4(mn0,k0),ch5(mn0,k0),bnd3(mn0),
     +     bnd4(mn0,k0),bnd5(mn0,k0)
      common/reffield/ref3(mn0),ref4(mn0,k0),ref5(mn0,k0)
      common/pref/irefs
      common/unit1/isdo,iref,istc,istr,ivoc,ivor,isbc,isbr,ivbc,ivbr
      character*80 ident,fniin,fniorr,fname,fmlrrs
      dimension pdt(4),nplev(4),iplev(4,k02),iplot(k02),lplot(k0)
      dimension uqg(k0),vqg(k0),tqg(k0)
c
      equivalence (nprnt, nplot)
c
c --- rdt = time at which to restart
c --- iin = unit from which to read job specifications
c --- iout = unit on which to write job results
c --- idff = unit on which to write full field diagnostics
c --- idif = unit on which to write interior field diagnostics
c --- iorr = unit from which to read restart data
c --- iowr = unit on which to write restart data
c --- iplt = unit on which to write plot data
c --- iced = no of times diagnostic routine has been called
c --- iwrite = write title for plot data (0=yes/1=no)
c --- isdo = unit on which to write interior fields
c
      data rdt/1000./,runid/'    ','    ','    ','    ','    ','    '/
      data restid/'    ','    ','    ','    ','    ','    ','    ',
     +            '    ','    '/
      data iin/3/,iout/4/,idff/7/,idif/8/,iorr/9/,iowr/10/,iplt/11/
      data iuvel/13/,isdo/23/,iref/25/,istc/26/,istr/27/,ivoc/28/
      data ivor/29/,isbc/30/,isbr/31/,ivbc/32/,ivbr/33/
      data iced/0/, iwrite/0/, idiag/0/, irest/0/, ifndd/0/
      data icnt,itcnt,ibcnt/0,0,0/
c
c ---   open units for input/output and read * file names
c
c      ident    used to generate the file names for the
c	        listing file (*.out), full field diagnostics
c	        (*.dff), interior field diagnostics (*.dif),
c	        restart output file (*.wrs), store interior fields
c               (*.sdo) and plot output
c		file (*.plt)
c
c	fniin	the parameter input file
c
c	fniorr  the input restart file
c
c
        write(idstring,
     +  '(''@(#)Version 1.13 of the Harvard Open Ocean Model'')')
        write(*,'('' Enter identifier for file names    '',$)')
	read(*,'(a)')ident
        write(*,'('' Enter name of input parameter file '',$)')
	read(*,'(a)')fniin
        write(*,'('' Enter name of restart file         '',$)')
	read(*,'(a)')fniorr
c
      iflen=lnblk(ident,80)
      open(unit=iin,file=fniin,form='formatted',status='old')
      fname=ident(:iflen)//'.out'
      open(unit=iout,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.dff'
      open(unit=idff,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.dif'
      open(unit=idif,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.plt'
      open(unit=iplt,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.sdo'
      open(unit=isdo,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.ref'
      open(unit=iref,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.stc'
      open(unit=istc,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.str'
      open(unit=istr,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.voc'
      open(unit=ivoc,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.vor'
      open(unit=ivor,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.sbc'
      open(unit=isbc,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.sbr'
      open(unit=isbr,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.vbc'
      open(unit=ivbc,file=fname,form='formatted',status='new')
      fname=ident(:iflen)//'.vbr'
      open(unit=ivbr,file=fname,form='formatted',status='new')
c

c
c
c --- initialize timing and get specs for this run
c --- an end code read in estart terminates processing
c
   10 continue
      icalc=0
      isd=0
      ibeg=0
      irefs=0
      iasday=50
      iastep=50
      imf=0
      icit=0
      ntitl=0
      call timstamp(runid)
c
c --- read and write info on data cards
c
      write(iout, *)idstring
      write(iout, 600)runid
      write(iout,'('' File identifier:    '',a,/,
     +             '' Parameter file :    '',a,/,
     +             '' Restart file   :    '',a,/,)')ident,fniin,fniorr
c
c
    1 continue
      read(iin, *, end=120)icard
      if(icard.le.0 .or. icard.gt.19) go to 132
      go to(101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 
     +  111, 112, 113, 114, 115, 116, 117, 118, 119), icard
c
  101 continue
      ntitl=ntitl+1
      if(ntitl.gt.2)ntitl=2
      read(iin, 501, end=130)(titl(ititl,ntitl),ititl=1,20)
      write(iout, 601)(titl(ititl,ntitl),ititl=1,20)
      go to 1
c
  102 continue
      read(iin, *, end=130)alpha, beta, rf
      write(iout, 602)alpha, beta, rf
      go to 1
c
  103 continue
      read(iin, *, end=130)m, n, kz
      write(iout, 603)m, n, kz
      kzp1=kz+1
      go to 1
c
  104 continue
      read(iin, *, end=130)(hz(i), i=1, kz)
      write(iout, 604)kz
      write(iout,'(10x,4(f9.3,2x))') (hz(i), i=1, kz)
      go to 1
c
c --- note that sigz(1) is no longer scsig and the
c --- addition of sigz(1) and sigz(kz+1)
c
  105 continue
      read(iin, *, end=130)scsig, (sigz(i), i=1, kzp1)
      write(iout,605)scsig
      write(iout,'(10x,4(f9.3,2x))') (sigz(i),i=1,kzp1)
      go to 1
c
  106 continue
      read(iin, *,end=130)iltord,iltfrq,iltcnt
      write(iout,606)iltord,iltfrq,iltcnt
      go to 1
c
c --- real array for driver
c
  107 continue
      read(iin,*,end=130)(pram(i),i=1,16)
      write(iout,607)(pram(i),i=1,16)
      go to 1
  108 continue
      read(iin,*,end=130)dt,tstart,tmax,xbasin
      write(iout,608)dt,tstart,tmax,xbasin
      go to 1
c
c --- switches for finite diff./colloc., f-/beta-plane, diagnostics,
c --- top and bottom density, top vertical velocity, and bottom
c --- topography.
c
  109 continue
      read(iin,*,end=130)ifdiff,ifpert,ifrst,iftop,ifbot,
     +                   iftvv,ifrel,ifsbl,ifeva,ifsdo
      if(iftvv.ne.0)iftop=1
      if(ifrel.ne.0)ifbot=1
      write(iout,609)ifdiff,ifpert,ifrst,iftop,ifbot,
     +                   iftvv,ifrel,ifsbl,ifeva,ifsdo
      if(ifrst.eq.0)go to 1
      open(unit=iorr,file=fniorr,form='unformatted',status='old')
      read(iorr,end=190)restid,titl
      write(iout,618)restid,titl
      go to 1
c
c --- if collocation used in depth, read in eigenstuff:
c
  110 continue
      read(iin,*,end=130)(eigval(i),i=1,kz)
      do 1110 i=1,kz
        read(iin,*,end=130)(amat(j1,i),j1=1,kz)
 1110 continue
      go to 1
c
c --- parameters for diagnostic options
c
  111 continue
      read(iin,*,end=130)nddt,(ddt(i),i=1,iabs(nddt))
      if(nddt.lt.0)ifndd=1
      nddt=iabs(nddt)
      if(nddt.gt.0)go to 1111
      write(iout,1611)
      goto 1
 1111 continue
      write(iout,610)
      do 2111 i=1,nddt
        write(iout,611)ddt(i)
 2111 continue
      go to 1
c
c --- parameters for array printing options
c
  112 continue
      read(iin,*,end=130)nadt,(adt(i),i=1,nadt)
      if(nadt.gt.0)go to 1112
      write(iout,1612)
      go to 1
 1112 continue
      write(iout,612)
      read(iin,*,end=130)(nalev(i),(ialev(i,j),j=1,nalev(i)),i=1,nadt)
      do 2112 i=1,nadt
        write(iout,613)adt(i),(ialev(i,j),j=1,nalev(i))
 2112 continue
      go to 1
c
c --- parameters for array plotting options
c
  113 continue
      read(iin,*,end=130)npdt,(pdt(i),i=1,npdt)
      if(npdt.gt.0)go to 1113
      write(iout,1614)
      go to 1
 1113 continue
      write(iout,614)
      read(iin,*,end=130)(nplev(i),(iplev(i,j),j=1,nplev(i)),i=1,npdt)
      do 2113 i=1,npdt
        write(iout,613)pdt(i),(iplev(i,j),j=1,nplev(i))
 2113 continue
      go to 1
c
c --- rdt = time interval at which to print restart information
c
  114 continue
      read(iin,*,end=130)rdt
      write(iout,616)rdt
      go to 1
c
  115 continue
      read(iin,*,end=130)theta
      write(iout,617)theta
      theta=theta*(3.1415926/180.)
      goto 1
c
  116 continue
      read(iin,*)rlat0,rlng0,t0,v0,dhor,ht,time0,rho0
      go to 1
c
  117 continue
      read(iin,*)tsurf, (xc(i),i=1,kzp1)
      write(iout,623)tsurf
      write(iout,"(5(e13.6,',  '),/)" ) (xc(i),i=1,kzp1)
      go to 1
  118 continue
      read(iin,*)salsur, (eta(i),i=1,kzp1)
      write(iout,624)salsur
      write(iout,"(5(e13.6,',  '),/)") (xc(i),i=1,kzp1)
      go to 1
  119 continue
      read(iin,*)cmdt
      if (cmdt .gt. 0.) then
         iflen=lnblk(ident,80)
         fname=ident(:iflen)//'.vel'
         open(unit=iuvel,file=fname,form='formatted',status='new')
        endif
      go to 1
c
c
  120 continue
      write(iout,121)
  121 format(' EOF encountered reading icard')
      stop
c
  130 continue
      write(iout,131)
  131 format(' EOF encountered reading input variable')
      stop
c
c --- come here after reading input parameters
c
  132 continue
      if(rdt.le.(tmax-tstart+dt/2))then
        fname=ident(:iflen)//'.wrs'
        open(unit=iowr,file=fname,form='unformatted',status='new')
      endif
c
c
      if(t0.eq.0.0)go to 136
      write(iout,625)rlat0,rlng0,t0,v0,dhor,ht,time0,rho0
c
      pi=acos(-1.0)
      f0star=4*pi*sin(rlat0*pi/180.)/86164.
      betstr=4*pi*cos(rlat0*pi/180.)/86164./6356912.
      r0=1.0/(t0*f0star)
      astar=t0*v0/dhor
      bstar=t0*betstr*dhor
      write(iout,626)f0star,betstr,r0,astar,bstar
c
      dt=dt*86400./t0
      tstart=(tstart-time0)*86400./t0
      tmax=(tmax-time0)*86400./t0


      do 135 i=1,4
        ddt(i)=ddt(i)*86400./t0
        adt(i)=adt(i)*86400./t0
        pdt(i)=pdt(i)*86400./t0
  135 continue
c
      cmdt=cmdt*86400./t0
c
      rdt=rdt*86400./t0
      go to 140
c
c --- translation for printed output if not dimensional input
c
  136 continue
      time0=0.0
      t0=86400.0
      ht=1.0
      dhor=1.0
      v0=1.0
c
  140 continue
c
c --- set up the stratification and associated matrices and constants
c
      call set3d(kz,ifdiff)
c
c --- constants for common/parm/
c
      mm1=m-1
      mm2=m-2
      mp1=m+1
      msq=m*n
c
      nm1=n-1
      nm2=n-2
      np1=n+1
c
      hy=xbasin/float(m-1)
      rea=hy*float(nm1*mm1)
      kzm1=kz-1
      kzp2=kz+2
c
c --- bottom topography and surface density variation is only
c --- implemented for finite difference method
c
      if(iftop.eq.0 .and. ifbot.eq.0)goto 146
      if(ifdiff.eq.0)goto 142
      do 141 i=1,kz
        topfct(i)=0.
        botfct(i)=0.
  141 continue
      topfct(1)=1./hz(1)
      botfct(kz)=1./hz(kz)
      goto 146
  142 continue
      write(iout,143)
  143 format('  Non-uniform surface or bottom density implemented',
     +  /,'  only for finite differences.')
      stop
  146 continue
c
c --- level printing and plotting information. allows for printing
c --- and plotting of density fields.
c
      do 151 na=1,nadt
        do 148 i=1,kzp2
          iprnt(i)=0
  148   continue
        do 149 i=1,nalev(na)
          lp=ialev(na,i)
          if(lp.eq.kzp1)lp=kzp2
          if(lp.eq.0)lp=kzp1
          iprnt(lp)=1
  149   continue
        do 150 i=1,kzp2
          ialev(na,i)=iprnt(i)
  150   continue
  151 continue
c
      do 155 np=1,npdt
        do 152 i=1,kzp2
          iplot(i)=0
  152   continue
        do 153 i=1,nplev(np)
          lp=iplev(np,i)
          if(lp.eq.kzp1)lp=kzp2
          if(lp.eq.0)lp=kzp1
          iplot(lp)=1
  153   continue
        do 154 i=1,kzp2
          iplev(np,i)=iplot(i)
  154   continue
  155 continue
c
c --- array 'mast'
c
      mast(1)=1
      mast(2)=m
      mast(3)=-1
      mast(4)=-m
c
c --- array 'nast'
c
      nast(1)=m
      nast(2)=-1
      nast(3)=-m
      nast(4)=1
c
c --- array 'icorn'
c
      icorn(1)=1
      icorn(2)=m
      icorn(3)=msq
      icorn(4)=msq-m+1
c
c --- array 'ncorn'
c
      ncorn(1)=m+2
      ncorn(2)=2*m-1
      ncorn(3)=msq-m-1
      ncorn(4)=msq-2*m+2
c
c --- coordinate values for 'latitude'
c
      ml=m
      if(n.gt.m)ml=n
      do 159 j=1,ml
        yy(j)=hy*float(j-1)-.5*xbasin
  159 continue
c
c --- define the meridional coordinate under model rotation
c
      ip=0
      do 1592 j=1,n
        do 1591 i=1,m
          ip=ip+1
          pvort(ip)=(yy(j)*cos(theta)+yy(i)*sin(theta))*beta
 1591   continue
 1592 continue
c
c --- constants for common/solver/
c
      vbm(1)=3.5
      vbn(1)=3.5
      do 160 i=3,mm2
        vam(i-1)=1./vbm(i-2)
        vbm(i-1)=4.-vam(i-1)
  160 continue
      vam(mm2)=1./vbm(mm2-1)
      vbm(mm2)=vbm(1)-vam(mm2)
      do 165 i=3,nm2
	van(i-1)=1./vbn(i-2)
	vbn(i-1)=4.-van(i-1)
  165 continue
      van(nm2)=1./vbn(nm2-1)
      vbn(nm2)=vbn(1)-van(nm2)
c
c --- initialize the streamfunction and vorticity arrays
c
      if(ifrst.ne.0)go to 185
c
c --- if this is an original calculation come here.
c
      tim=tstart
      tm1=tim-dt
c --- read ref data
c
      if(irefs.eq.1)then
      open(unit=51,file='tpersd.dat',form='formatted',status='old')
      endif

c ---read vertical profile
c
      open(unit=50,file='eof_dh3.dat',form='formatted',status='old')
c --- read input data
c
      open(unit=24,file='tpersf_13.dat',form='formatted',status='old')
      open(unit=52,file='dh3v1015d.dat',form='formatted',
     1 status='old')
c
c --- read  fields: step 1
c

        call meshob(0) 

          do 167 k=1,kz
           do 166 j=1,msq
             ch4(j,k)=bnd4(j,k)
             ch5(j,k)=bnd5(j,k)
  166  continue
  167  continue
           do 170 j=1,msq
            ch3(j)=bnd3(j)
  170   continue

      if(iftop.eq.1)call infld(topden,1,tm1,7,0)
      if(ifbot.eq.1)call infld(botden,kz,tm1,8,0)
      if(ifrel.eq.1)call infld(relief,kz,tm1,6,0)
c
c --- initialize for time = tim -dt
c
      do 175 k=1,kz
        call infld(stream(1,k),k,tm1,1,0)
        call infld(vort(1,k),k,tm1,2,0)
        call vorteq(k)
  175 continue
c
c --- initialize sbl model
c
      if( ifsbl .gt. 0 ) call infld(r,k,tm1,10,0)

c ---      call infld(r,k,tm1,9,0)
c
c --- initialize for time tim 

       continue

      if(iftop.eq.1)call infld(topden,1,tim,7,0)
      if(ifbot.eq.1)call infld(botden,kz,tim,8,0)
c

      do 180 k=1,kz

        call infld(stream(1,k),k,tim,1,0)
        call infld(vort(1,k),k,tim,2,0)

  180 continue
c
c --- read step 2: 

            call meshob(0)
c
c --- initialize sbl model
c
      if( ifsbl .gt. 0 ) call infld(r,k,tim,10,0)

c -----      call infld(r,k,tim,9,0)

      go to 199
c -------------------------------------------------------------
c
c --- if this is a restart of a previous calculation come here.
c
  185 continue
c
c --- read the dt used in the previous calculation and insure that it
c --- varies from the present dt by less than .01%.
c
      read(iorr,end=190)odt,otstart
      if(abs((odt-dt)/dt).le.0.0001)go to 186
      write(iout,619)odt,dt
      stop
  186 continue
c
c --- read in the relief if bottom topography is to be used.
c
      if(ifrel.ne.0)read(iorr,end=190)(relief(ip),ip=1,msq)
c
c --- tlx is the time corresponding to the last full set of restart
c --- information read.  if tlx=0.0 then no full sets have been read.
c
      tlx=0.0
c
c --- read the time of the first set of information.
c
      read(iorr,end=190)tx
  187 continue
c
c --- now read in the set of information corresponding to tx.
  
c --- if end of file is detected, then go to section which determines
c --- whether this is an error or this indicates that the last set 
c --- should be used.
c
      read(iorr,end=189)((stream(ip,k),ip=1,msq),k=1,kz),
     +                  ((vort(ip,k)  ,ip=1,msq),k=1,kz),
     +                  ((prvjac(ip,k),ip=1,msq),k=1,kz)
      if(iftop.ne.0)read(iorr,end=189)(topden(ip),ip=1,msq),
     +                                (prvtop(ip),ip=1,msq)
      if(ifbot.ne.0)read(iorr,end=189)(botden(ip),ip=1,msq),
     +                                (prvbot(ip),ip=1,msq)
c
c --- a full set of information has been read, hence update tlx.
c
      tlx=tx
c
c --- if time = otstart (tstart of the restart data.), then 
c --- we are searching for the last full set of
c --- of information, while if it is not equal, we are searching for a
c --- set within + or - dt of tstart.
c
      if(tstart.ne.otstart)go to 188
c
c --- we are looking for the last full set, hence if the end of
c --- file is found when reading the next time, all of the information
c --- has been stored properly (except the time) hence begin the 
c --- calculation.
c
      read(iorr,end=198)tx
c
c --- if no end of file, read in the next set.
c
      go to 187
  188 continue
c
c --- we are looking for a set witin + or - dt of time, hence check tx.
c --- also, if tx is greater than time+dt, there is no hope of
c --- satisfying this condition.
c
      if(tx.gt.tstart+(1.5*dt))go to 191
      if(abs(tstart-tx).le.abs(1.5*dt))go to 198
      read(iorr,end=191)tx
      go to 187
  189 continue
c
c --- here if eof in the middle of a set of restart information.
c --- if looking for a specific time (time not=0) then this is an error
c --- since a set of information corresponding to this time does not 
c --- exist (no complete set).
c --- if tlx=0 then there are no complete sets of data in the file.
c
      if(tlx.eq.0.0)go to 190
      if(tstart.ne.otstart)go to 191
c
c --- if looking for the last complete set, then this corresponds to
c --- time=tlx, hence start over and find this specific time.
c
      tstart=tlx
      rewind(iorr)
      read(iorr)titl
      go to 185
  190 continue
c
c --- here if no complete sets of restart information.
c
      write(iout,621)
      stop
  191 continue
c
c --- here if unable to find a data set corresponding to tstart.
      write(iout,622)tstart,tx
      print *,'tstart,tx,tim,dt,tlx,otstart,odt'
      print *,tstart,tx,tim,dt,tlx,otstart,odt
      stop
c
c --- restart information has been read in at time tlx, so set time and
c --- enter main computation loop.
c
  198 continue
      tim=tlx
      tstart=otstart
c
c --- initialize sbl model if necessary
c
      if (ifsbl.gt.0) then
         call infld(r,k,tim,10,0)
         call mlrrs
       endif
c
c --- end of  restart initialization
c
c -------------------------------------------------------------
c
c -------------------------------------------------------------
c --- top of main computation loop
c


  199 continue
c ------------------------------------------------------------
c     output section
c
         if(isd.ne.iasday) go to 200

          do 1991 j=1,msq
           ch3(j)=bnd3(j)
 1991   continue

            do 1993 i=1,kz
            do 1992 j=1,msq
                ch4(j,i)=bnd4(j,i)
                ch5(j,i)=bnd5(j,i)
 1992        continue
 1993     continue

  200  CONTINUE
c
c --- if desired, calculate the diagnostics for this time step.
c
      if(nddt.le.0)go to 210
      if(idiag.ne.0)go to 204
      idiag=1
      if(ifndd.ne.0)go to 201
      write(idff,801)runid,restid,
     .                   ((titl(ititl,jtitl),ititl=1,20),jtitl=1,2)
      write(idff,800)
      write(idif,802)runid,restid,
     .                   ((titl(ititl,jtitl),ititl=1,20),jtitl=1,2)
      write(idif,800)
      go to 202
c
  201 continue
      write(idff)'full field integrated quantities    '
      call ndinfo(idff)
      write(idif)'interior field integrated quantities'
      call ndinfo(idif)
c
  202 continue
  204 continue
      do 205 i=1,nddt
c --        if(amod(tim-tstart+dt/2,ddt(i)).ge.dt)go to 205
        if(icalc.eq.0) go to 2041

        if(imf.ne.iastep) go to 205
          imf=0
 2041  continue

        
        rz=stream(16,1)-ch4(16,1)

        write(isbr,'(8E15.6)')tim*t0/86400.0,rz,(ch4(16,k),k=1,kz)

        write(istc,'(7E15.6)') tim*t0/86400.0,(stream(416,k),
     +  k=1,kz)
        write(ivoc,'(775f10.4)') (vort(j,1),j=1,msq)

        write(isbc,'(7E15.6)') tim*t0/86400.0,(stream(16,k),k=1,kz)

         call ediag
         write(ivbc,'(775f10.4)')(stream(j,1),j=1,msq)
         if(irefs.eq.1)then
         write(ivbr,'(775f10.4)')(ref4(j,1),j=1,msq)
         write(istr,'(7E15.6)') tim*t0/86400.0,(ref4(416,k),k=1,kz)
         write(ivor,'(775f10.4)')(ref5(j,1),j=1,msq)
         else
         write(ivbr,'(775f10.4)')(ch4(j,1),j=1,msq)
         write(istr,'(7E15.6)') tim*t0/86400.0,(ch4(447,k),k=1,kz)
         write(ivor,'(775f10.4)')(ch5(j,1),j=1,msq)
c         write(ivor,'(775f10.4)')((ch5(j,k),j=1,msq),k=1,kz)
         endif 
        
        iced=iced+1
        go to 210
  205 continue
  210 continue
c
c --- if desired, print out the arrays
c
      if(nadt.le.0)go to 320
      ifprnt=0
      do 220 j=1,kzp2
        iprnt(j)=0
  220 continue
      do 240 i=1,nadt
        if(amod(tim-tstart+dt/2,adt(i)).ge.dt)go to 240
        if(nalev(i).le.0)go to 240
        ifprnt=1
        do 230 j=1,kzp2
          iprnt(j)=iprnt(j)+ialev(i,j)
  230   continue
  240 continue
      if(ifprnt.ne.1)go to 320
      nprnt=0
      do 250 i=1,kz
        if(iprnt(i).eq.0)go to 250
        nprnt=nprnt+1
        lprnt(nprnt)=i
  250 continue
c
c --- streamfunction
c
      if(icalc.ne.0)go to 259
      rmax=0.0
      do 256 k=1,kz
        do 255 ip=1,msq
          if(abs(stream(ip,k)).gt.rmax)rmax=abs(stream(ip,k))
          if(abs(vort(ip,k)).gt.rmax)rmax=abs(vort(ip,k))
  255   continue
  256 continue
c
      fac=1000.
      if(rmax.eq.0.0)go to 258
      ipower=5-log10(rmax)
      fac=10.**ipower
  258 continue
      write(iout,707)fac
c
c --- if there are no levels to print(nprnt=0) then go to density print.
c
      if(nprnt.eq.0)go to 311
c
  259 continue
c
      write(iout,701)tim*t0/86400.+time0,(lprnt(i),i=1,nprnt)
      do 260 i=1,nprnt
        write(iout,700)lprnt(i)
        call harvpt(stream(1,lprnt(i)),m,n,fac)
  260 continue
c
c --- vorticity
c
      write(iout,702)tim*t0/86400.+time0,(lprnt(i),i=1,nprnt)
      do 270 i=1,nprnt
        write(iout,700)lprnt(i)
        call harvpt(vort(1,lprnt(i)),m,n,fac)
  270 continue
c
c --- perturbation streamfunction
c
      if(ifpert.eq.0)goto 311
      write(iout,703)tim*t0/86400.+time0,(lprnt(i),i=1,nprnt)
      do 290 i=1,nprnt
        k=lprnt(i)
        call infld(r,k,tim,3,0)
        do 280 ip=1,msq
          r(ip)=stream(ip,k)-r(ip)
  280   continue
        write(iout,700)k
        call harvpt(r,m,n,fac)
  290 continue
c
c --- perturbation vorticity
c
      write(iout,704)tim*t0/86400.+time0,(lprnt(i),i=1,nprnt)
      do 310 i=1,nprnt
        k=lprnt(i)
        call infld(r1,k,tim,4,0)
        do 300 ip=1,msq
          r1(ip)=vort(ip,k)-r1(ip)
  300   continue
        write(iout,700)k
        call harvpt(r1,m,n,fac)
  310 continue
c
c --- now print density fields if required.
c
  311 continue
      if(iprnt(kzp1).eq.0)go to 312
      write(iout,705)tim*t0/86400.+time0
      call harvpt(topden,m,n,fac)
  312 continue
      if(iprnt(kzp2).eq.0)go to 320
      write(iout,706)tim*t0/86400.+time0
      call harvpt(botden,m,n,fac)
c
c --- come to line 320 if no array printing or when finished.
c
  320 continue
      if(npdt.le.0)go to 370
      ifplot=0
      do 330 j=1,kzp2
        iplot(j)=0
  330 continue
      do 350 i=1,npdt
c --        if(amod(tim-tstart+dt/2,pdt(i)).ge.dt)go to 350
           if(icalc.eq.0) go to 335
       if(isd.ne.iasday) go to 350

  335      if(nplev(i).le.0)go to 350
        ifplot=1
        do 340 j=1,kzp2
          iplot(j)=iplot(j)+iplev(i,j)
  340   continue
  350 continue
      if(ifplot.ne.1)go to 370
      nplot=0
      do 360 i=1,kz
      if(iplot(i).eq.0)go to 360
        nplot=nplot+1
        lplot(nplot)=i
  360 continue
c
c --- insert calls to plot fields here.
c
c --- the variable nplot specifies the number of levels to be plotted,
c --- while lplot(i),i=1,nplot specify the levels to be plotted.
c
      if(iwrite.ne.0)go to 365
      iwrite=1
      call ndinfo(isdo)
      call ndinfo(iplt)
      call ndinfo(iref)
      if(ifrel.ne.0)call xwrite(relief,msq,6,6,6,6,6,iplt)
c
  365 continue
      if(nplot.eq.0)go to 367
      do 366 i=1,nplot
        k=lplot(i)
        call xwrite(stream(1,k),msq,tim,1,k,k,k,iplt)
        call xwrite(vort(1,k),msq,tim,2,k,k,k,iplt)
c ----
        if(irefs.eq.1)then
         do 3650 j=1,msq
          r1(j)=ref4(j,k)
 3650   continue
         else
        call infld(r1,k,tim,1,0)
         endif
c ----
        call xwrite(r1,msq,tim,1,k,k,k,iref)
         do 3651 j=1,msq
           r1(j)=stream(j,k)-r1(j)
 3651  continue
        call xwrite(r1,msq,tim,1,k,k,k,isdo)
c ----
       if(irefs.eq.1)then
         do 3652 j=1,msq
          r1(j)=ref5(j,k)
 3652  continue
        else
        call infld(r1,k,tim,2,0)
         endif
c --------
        call xwrite(r1,msq,tim,2,k,k,k,iref)
         do 3653 j=1,msq
           r1(j)=vort(j,k)-r1(j)
 3653   continue
        call xwrite(r1,msq,tim,2,k,k,k,isdo)
  366 continue
c
c --- now plot the top and/or bottom density fields.
c
  367 continue
      if(iplot(kzp1).eq.0)go to 368
      call xwrite(topden,msq,tim,7,7,7,7,iplt)
c ----
       if(irefs.eq.1)then
        do 3670 j=1,msq
         r1(j)=ref3(j)
 3670   continue
       else
      call infld(r1,k,tim,7,0)
       endif
c ----
      call xwrite(r1,msq,tim,7,7,7,7,iref)
       do 3671 j=1,msq
         r1(j)=topden(j)-r1(j)
 3671  continue       
      call xwrite(r1,msq,tim,7,7,7,7,isdo)
  368 continue
      if(iplot(kzp2).eq.0)go to 370
      call xwrite(botden,msq,tim,8,8,8,8,iplt)
c
c --- come here if no plots or if finished with plots.
c
  370 continue
c
c --- velocity record printing
c
      if (cmdt .gt. 0.) then
         if (amod(tim-tstart+dt/2,cmdt).lt.dt) then
            do 372 k=1,kz
               call qgvel(stream(1,k),m,n,hy,v0,uqg(k),vqg(k))
  372        continue
            call qgtemp(tqg)
            write(iuvel,373)tim*t0/86400.,(uqg(k),vqg(k),tqg(k),k=1,kz)
           endif
        endif
  373 format(99e12.4)
c
c --- restart information is printed here.
c
      if(rdt.gt.999.)go to 380
      if(amod(tim-tstart+dt/2,rdt).ge.dt)go to 380
      if(irest.ne.0)go to 375
      irest=1
      restid(1)=' res'
      restid(2)='tart'
      restid(3)=' of '
      do 374 i=1,6
        restid(i+3)=runid(i)
  374 continue
      write(iowr)restid,titl
      write(iowr)dt,tstart
      if(ifrel.ne.0)write(iowr)(relief(ip),ip=1,msq)
      if(icalc.le.1)go to 380
  375 continue
      write(iowr)tim
      write(iowr)((stream(ip,k),ip=1,msq),k=1,kz), 
     +         ((vort(ip,k),  ip=1,msq),k=1,kz),
     +         ((prvjac(ip,k),ip=1,msq),k=1,kz)
      if(iftop.ne.0)write(iowr)(topden(ip),ip=1,msq),
     +                         (prvtop(ip),ip=1,msq)
      if(ifbot.ne.0)write(iowr)(botden(ip),ip=1,msq),
     +                         (prvbot(ip),ip=1,msq)
      if (ifsbl .gt. 0) call mlwrs
  380 continue
c
c     end of output section
c ---------------------------------------------------------------
       if(icalc.eq.0) go to 440

c -- if(amod(tim-tstart+dt/2,asday).ge.dt) go to 440
c
c --- read new boundary conditions
c
      if(isd.ne.iasday) go to 440

          isd=0

  410  continue
           call meshob(0)

  440  CONTINUE

      if(tim.ge.tmax) go to 9999
c --------------------      tim=tim+dt
      icalc=icalc+1
      tim=dt*icalc
      isd=isd+1
      imf=imf+1
c
c --- test for divergence
c
      midlev=kz/2
      if(midlev.lt.1)midlev=1
      if(abs(vort(msq/2,midlev)).ge.1.e6) go to 998
c
c --- call loop3d to perform calculations
c
      call loop3d
c
c --- compute sbl if needed
c
      if (ifsbl.gt.0) call infld(r,k,tim,10,0)
c
c --- return to top of main loop
c
      go to 199
c
c --- come to line 998 if vorticity blows up
c --- come to line 9999 for normal exit
c
  998 continue
      write(iout,999)tim*t0/86400.+time0
  999 format(1h ,' Error during run: Vorticity field too big:TIME=',
     +e10.3,/)
 9999 continue
c -----------------------------------------------------------------
c
c --- input format statements for data cards
c
  500 format(40i2)
  501 format(20a4)
  502 format(8f10.0,/,8f10.0)
  503 format(i5,4f10.0)
  504 format(11f8.0)
  505 format(16i5)
  506 format(f10.2)
  507 format(12f10.0)
c
c --- output format statements for data cards
c
  600 format(' ',6a4)
  601 format(/,20a4)
  602 format(/'Ratio of characteristic to advective time scales=',f8.4,
     + //,'Ratio of characteristic to rossby wave time scales=',f8.4,//,
     +    'Bottom friction=',f8.4)
  603 format(/'Number of horizontal grid points, m=',i3,/,
     +        ' Number of merdional grid points, n=',i3,/,
     +        '      Number of vertical levels, kz=',i3)
  604 format(/'No of layers= ',i3,'  Layer Depths : ')

  605 format(/'Stratification scale = ',1pe11.4,/,' Sigma values:')

  606 format(/,'Shapiro filter:  Order=',i3,5x,'Frequency=',i3,
     +       5x,'every ',i3,' time steps.')
  607 format(/,'Array PRAM:  ',4(1pg13.6,1x),/,13x,4(g13.6,1x),/,
     +       13x,4(g13.6,1x),/,13x,4(g13.6,1x))
  608 format(/'Dt= ',f8.5,'   Tstart= ',f12.5,'   Tmax= ',f12.5,
     +         //'Xbasin= ',f8.5)
  609 format(/,'Switches:     ',
     +       /,'vertical: finite diff(1)/colloc(0) ifdiff: ',i1,
     +       /,'     compute perturbation field(1) ifpert: ',i1,
     +       /,'            restart at time rdt(1)  ifrst: ',i1,
     +       /,'        non-uniform top density(1)  iftop: ',i1,
     +       /,'     non-uniform bottom density(1)  ifbot: ',i1,
     +       /,'     use  top vertical velocity(1)  iftvv: ',i1,
     +       /,'        use  bottom topography(1)   ifrel: ',i1,
     +       /,' compute surface boundary layer(1)  ifsbl: ',i1,
     +       /,'         compute fields for EVA(1)  ifeva: ',i1,
     +       /,'          write internal fields(1)  ifsdo: ',i1)

  610 format(/'Diag printing:'//,4x,'Interval')
  611 format(3x,f10.6)
 1611 format(/'No diagnostic printing')
  612 format(/'Array printing:'//,4x,'Interval',6x,'levels')
 1612 format(/'No array printing')
  613 format(3x,f10.6,4x,9i3)
  614 format(/'Array plotting:'//,4x,'Interval',6x,'levels')
 1614 format(/'No array plotting')
  615 format(/'       Plotting interval    No of levels   ',40i3)
  616 format(/'Interval at which to write restart information:',1pe10.3)
  617 format(/' Domain is rotated',f12.3,' degrees cck')
  618 format(/,9a4,/,' ',20a4)
  619 format(/' DT specified by input and specified by restart data',
     +     /,'  File differ by more than .01%.  odt,dt:',1pe15.8,', ',
     +       e15.8)
  621 format(/' No complete sets of data.')
  622 format(/' Unable to find set of data, whose time is within DT',
     + /,'  of the time specified.  TSTART, TX:',1pe15.8,', ',e15.8)
  623 format(/' Vertical temperature profile ',//,
     +       ' Surface temperature: ',1pe13.6,' deg c.',
     +      /' Vertical derivative of temperature: ( deg c. / meter)')
  624 format(/' Vertical salinity profile ',//,
     +       ' Surface salinity: ',1pe13.6,' .001',
     +      /' Vertical derivative of salinity: ( .001 / meter)')
  625 format(/' **** Run includes dimensional external parameters ***',
     +     //,' The units of time specifed in the input data file are',
     +     /,'  in days.  ',
     +     /,'  Warning : The input ALPHA and  BETA should be checked',
     +     /,'  with those computed from the dimensional input.',//,
     +     //,' Center grid Latitude, longitude: ',2(f12.4),
     +             '(deg north and west)',
     +     //,' T0 ',e14.6,' s, V0',e14.6,' m/s,',/,
     +        '  L ',e14.6,' m,  H',e14.6,' m',
     +     //,' Time offset ',0pf10.2,' (day), Mean density ',
     +       1pe14.6,' g/cm/cm/cm')
  626 format(/' Parameters determined from the dimensional input',/,
     +       /'     F0* ',1pe14.6,'(/s), BETA* ',e14.6,'(/s/m)',
     +       /'     R0* ',e14.6,
     +       /'     A*  ',e14.6,'         B* ',e14.6)
c
c --- format statements for array printing
c
  700 format(/'Level ',i2)
  701 format(/,30x,'Streamfunction',5x,5hTIME=,f10.4,
     +' Levels',9(1x,i2))
  702 format(/,30x,'Vorticity',5x,5hTIME=,f10.4,' Levels',9(1x,i2))
  703 format(/,30x,'Perturbation streamfunction',5x,5hTIME=,f10.4,
     +' Levels',9(1x,i2))
  704 format(/,30x,'Perturbation vorticity',5x,5hTIME=,f10.4,
     +' Levels',9(1x,i2))
  705 format(/,30x,'Top density',8x,5hTIME=,f10.4)
  706 format(/,30x,'Bottom density'5x,5hTIME=,f10.4)
  707 format('Scale = ',1pe6.0)
c
c --- formats for diagnostic printing.
c
  800 format(/'Level ',
     +             40hICALC   T    I(z)      I(q)    I(nrg)   ,
     +             36hMAX(v)   DI(nrg)  MAXD(s)  RMSD(s)  ,
     +             43hRMSD(z)  STD(s)   RATIO   ABS CORR TEN CORR)
  801 format(6a4,9a4,/,' ',20a4,/,' ',20a4,//,
     .                                ' Diagnostics for the full field')
  802 format(6a4,9a4,/,' ',20a4,/,' ',20a4,//,
     .                            ' diagnostics for the interior field')
c
      stop
      end
