c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timegcm/mksrc/mksrc.f
c------------------------------------------------------------------
c
      program mksrc
c
c Modified from ~/timegcm/mksrc/mksrc.f for zp=-13 to +5 (kmx=37)
c Does not use tigcm (glbmean only, with u,v,w = 0.)
c
      include 'mksrc.h'
      include 'glbm.h'
      dimension iasf(13),iflnm(3),tndzp(ntndown)
      data iasf/13*1/
      data ihistvol/'hist'/
      character*80 errmsg
      data ludata/88/, luglbm/9/
c
c Get user input:
c
      call getinp
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,glon1,dlon,imx,
     +  gcmzp,zp1,dzp,kmx)
      call cpsetup(cpspval)
      call setag
c
c Set up gks:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c
c Open new file for history (assigned to lutime): 
c
      call openhist
      if (idatvol.gt.0) call opendatv
c
c History loop (ntms should = 1):
c
      do it=1,ntms
c
c Get global mean history from which to assign fields to timegcm:
c
        call getglbm(luglbm)
c
c Get bottom boundaries for data volume:
c
        if (idatvol.gt.0) then
          ieqnx = 0
          if (iday.eq.80) ieqnx = 1
c
c Get bottom boundaries for tn at ground, and tn, ht, and un at zp -13
c   from cira88:
c
          do k=1,ntndown
            tndzp(k) = tndzp1+(k-1)*tnddel
          enddo
          call getcira(lucira,ieqnx,glbmzp(1),gcmlat,jmx,tng,tnlb,htlb,
     +                 unlb,rholb,tndown,ntndown,tndzp,ipltcira,iframe)
c
c 9/27/91: remake tn, tng, ht, and un from Jeff Forbes code:
c
          call getforbe(ieqnx,glbmzp(1),gcmlat,jmx,tng,tnlb,htlb,unlb,
     +      ipltzatm,iframe)
c         call getforbe(ieqnx,-7.,gcmlat,jmx,tng,tnlb,htlb,unlb)
c
c Write data statements for Ray:
c         write(ludata,"('C')")
c         if (ieqnx.gt.0) then
c           write(ludata,"('C EQUINOX')")
c         else
c           write(ludata,"('C SOLSTICE')")
c         endif
c         write(ludata,"('C'/'C TN (DEG K)'/'C')")
c         write(ludata,"(8x,6f10.4)") tnlb
c         write(ludata,"('C'/'C HEIGHT (KM)'/'C')")
c         write(ludata,"(8x,6f10.4)") htlb
c         write(ludata,"('C'/'C UN (M/S)'/'C')")
c         write(ludata,"(8x,6f10.4)") unlb
c
c 9/91: lower boundaries for h2o,h2,co,co2,h,oh,o,o3,ho2,no2,no
c       are found from acd 2d model (see ~/acd2d), and are
c       stored on the data volume.
c acd2d() returns lower boundaries at hgtlb height and at tgcm
c         latitude grid in speclb(jmx,nspec) These are written
c         onto the end of the data volume (not every lat slice)
c
          call acd2d(nspec,ispecno,acdname,acdhist,luacd,gcmlat,jmx,
     +               hgtlb,speclb,rholb,ispecpsi,ipltacd,ieqnx,iframe) 
c
c Fill in ox, noz, hox:
c ox = o + o3; noz = no + no2; hox = oh + ho2 + h
c In speclb:  h2o, h2, co, co2,  h, oh,  o, o3, ho2, no2, no ch4
c               1,  2,  3,   4,  5,  6,  7,  8,   9,  10, 11, 12
c
          do j=1,jmx
            oxlb(j) = speclb(j,7) + speclb(j,8)    
            rnozlb(j) = speclb(j,11) + speclb(j,10)
            hoxlb(j) = speclb(j,6) + speclb(j,9) + speclb(j,5)
          enddo
c
c Plot lower boundaries:
          if (ipltlb.gt.0) call pltlb(botbound,jmx,nbot,gcmlat,glbmzp(1)
     +      ,hgtlb,lbname,ieqnx,tndown,ntndown,tndzp,iframe)
        endif ! datavol
c
c Latitude loop:
c
        do j=1,jmx
c       do j=1,2
c
c Define f() at i=3->74 for current j:
c 
          call settime(j)
c
c Set periodic points in f:
c 1,2 = 73,74 and 75,76 = 3,4:
c
          do ip=1,nfldsd
            do i=1,2
              f(i,:,ip) = f(i+imxm1,:,ip)
              f(i+imxp1,:,ip) = f(i+2,:,ip)
            enddo
          enddo
c
c Define fdatv for data vol:
c
          if (idatvol.gt.0) call mkdatv(j,it)
c
c Fields at 'time step minus 1' are same as just defined:
c   nfldsw = 43 = number of fields written to history file
c   nfldst = 44 = total number of fields (for dimensions)
c   nfldsd = 29 = number of fields defined by this code
c (fields 30-43 are same as 1-14)
c
          do ip=nfldsd+1,nfldsw
            f(:,:,ip) = f(:,:,ip-nfldsd)
          enddo
c
c Write header and summary (only once per history):
c
          if (idispose.gt.0) then
            if (j.eq.1) call hdrsum
c
c Write latitude slice, history and data volume:
c
            call wrhist(j)
            if (idatvol.gt.0) call wrdatv(j)
          endif
c
c End latitude loop:
c
          write(6,"('mksrc: processed lat slice ',i3,' lat=',f8.2)")
     +      j,gcmlat(j)
        enddo
        if (idispose.gt.0) call disptp(itimeout,ihistvol,1,lutime)
c
c Dispose data volume if necessary ('datvol' must match assign
c                                    in opendatv)
      if (idatvol.gt.0) then
        close(ludatv)
        lendat = lenstr(datvnm)
        if (idispose.gt.0) then
          call mswrite(ier,'datvol',datvnm(1:lendat),',ECRIDLEY',
     +      367,' ')
          if (ier.ne.0) then
            write(6,"('mksrc: error ',i3,' from mswrite for data vol')")
     +        ier
            call mserror(ier,errmsg)
            write(6,"('errmsg=',a)") errmsg
          else
            write(6,"('mksrc: wrote data vol ',a)") datvnm(1:lendat)
          endif
        endif
      endif
c
c End time loop:
      enddo
c
      call clsgks
      stop 'done'
      end