!
      module mk_v5d
!
! Make vis5d output file:
!
      use proc
      use hist,only: history,hdr
      use fields,only: field
      implicit none
!
      type vis5d
        character(len=120) :: 
     +    flnm,		! local disk name for v5d output file
     +    sendv5d,	! host:path to send v5d output (user: sendv5d)
     +    sendmsv5d	! mss path to send v5d output (user: sendms_v5d)
        character(len=8) :: 
     +    vtype		! Vertical coord type
          		!   user: v5d_vtype= 'reg_zp' or 'ireg_mb'
        real :: 
     +    fix_slt,	! opt slt at which to fix map data (user: v5d_slt)
     +    zprange(2)	! opt zp range (bot,top) (user: v5d_zprange)
      end type vis5d
      type (vis5d) v5d	! init by input and setflev in getflds
!
#if defined(SUN) || defined(LINUX)
      end module mk_v5d
#else
!
      contains
!-------------------------------------------------------------------
      subroutine mkv5d(f,nf,h,itime,mtimes,ier)
!
! Make v5d file containing 3d fields for each history for use w/
!   vis5d application.
! As in mkcdf.f, this routine is called for each history.
!   When called for 1st history (itime==1), the vis5d file is created, 
!   and the header is written (v5dcreate). For itime > 1, fields are 
!   added to the open file (v5dwrite). When called for the last
!   history (itime==ntimes), the file is closed after writing the
!   fields.
! Because the datestamps for all times are required at file creation
!   time, the yyddd of the 1st history is used, then is incremented
!   by 1 at day boundaries (i.e., when timestamp(i)<=timestamp(i-1).
! It may be desirable to limit the vertical extent of the fields.
!   This may be done in 2 ways: First, setting user input v5d%zprange(2) 
!   (bot,top) will limit all fields to within this range (see global 
!   vis5d file header var vert_args). Second, can set f(i)%vregion(2)
!   to further limit vertical range of selected fields. (see setflev in 
!   getflds.f). This is accomplished be setting field values to spval
!   when outside the vregion.
! The data may be optionally rotated to a fixed local time v5d%fix_slt
!   at center of map longitude axis (see function rotslt). Note because 
!   of the model 5 deg grid resolution, this works only when delta time 
!   divides evenly by 20 minutes (earth rotates 5 deg every 20 mins).
!   (v5d_slt is optional user input). Also note that vis5d will not
!   rotate map projections w/ time, so maps should not be displayed
!   with files in which v5d_slt was set.
!
! v5dcreate and v5dwrite, and other related functions are in v5d.c
!   and binio.c. As of 9/97, these are vis5d version 4.3 for SunOS 5.x, 
!   obtained from U.Wisconsin (http://www.ssec.wisc.edu/~billh/vis5d.html) 
!   (Principle author Bill Hibbard whibbard@macc.wisc.edu)
!   (see /e/foster/vis5d/vis5d-4.3 for source, doc, etc)
! Note this code executes on the Cray and does cray_to_ieee conversion
!   so the .v5d files may be read w/ vis5d on SUNS or SGI's.
!
      implicit none
!
! Args:
      integer,intent(in) :: nf,itime,mtimes(:,:)
      type(field),intent(in) :: f(nf)
      type(history),intent(in) :: h
      integer,intent(out) :: ier
!
! Externals:
      integer,external :: v5dcreate,v5dcreatesimple,v5dwrite,
     +                    v5dsetunits,v5dclose,v5dsetlowlev
      integer,external :: ixfind
!
! Interface block for external functions that return arrays:
      interface
        function rotslt(forig,slt,ut,gcmlon,nlon,dlon,spval)
          real :: rotslt(nlon)	! function result
          integer,intent(in) :: nlon
          real,intent(in) :: forig(nlon),slt,ut,gcmlon(nlon),dlon,spval
        end function rotslt
        function mkpmb(zp,nzp,p0,zprange,idimmb,spval,v5dhts)
          real :: mkpmb(idimmb)	! function result
          integer,intent(in) :: nzp,idimmb
          real,intent(in) :: zp(nzp),p0,zprange(2),spval
          real,optional,intent(out) :: v5dhts(idimmb)
        end function mkpmb
      end interface
!
! 5-D grid limits, must match those in v5d.h (also v5df.h, if used)
      integer,parameter :: 
     +  MAXVARS=30, 
     +  MAXTIMES=400, 
     +  MAXROWS=200,
     +  MAXCOLUMNS=200, 
     +  MAXLEVELS=100, 
     +  MAXVERTARGS=MAXLEVELS+1,
     +  MAXPROJARGS=100,
     +  IMISSING=-987654
      real,parameter :: MISSING=1.0E35
!
! Locals:
      character(len=8)   :: charmtime,charmtime1
      character(len=10) :: fnames(MAXVARS)
      integer :: stat,i,ii,ixf,j,k,kk,iier,maxnl,nlev,nflds=0
      integer :: compress=1
      integer,allocatable :: timestamp(:),datestamp(:)
      integer,allocatable,save :: nlevs(:),lowlevs(:)
      integer :: projection=1	! assume rectilinear CE for now
      real :: proj_args(MAXPROJARGS)
!
! vertical is int vertical coord type passed to v5d_create:
!   For v5d%vtype='reg_zp',  vertical=0, equally spaced generic units (zp)
!   For v5d%vtype='ireg_mb', vertical=3, unequally spaced millibars
! (vert_args are different for different vertical coord types)
!
      integer :: vertical
      real :: vert_args(MAXVERTARGS),
     +  v5d_mbs(MAXVERTARGS),v5d_hts(MAXVERTARGS)
      real :: flon(nlon)
      real,allocatable :: ftmp(:,:,:),vrange(:)
!
! Begin exec:
      ier = 0
!
! If this is not the 1st time, just add fields to existing v5d file,
      if (itime > 1) goto 100
!
! Otherwise, create the file and write the header:
!
! Construct a local v5d disk file name:
!   If only 1 time, v5d%flnm = "1sthistvol_mtime.v5d"
!   If multiple times, v5d%flnm = "1sthistvol_1stmtime-lastmtime.v5d",
!   For example: CCMS135_0751200-0770000.v5d
!
      call tail(h%mssvol,v5d%flnm)
      write(charmtime,"(i3.3,i2.2,i2.2)") h%mtime
      if (ntimes == 1) then
        v5d%flnm = trim(v5d%flnm)//'_'//trim(charmtime)//'.v5d'
      else
        write(charmtime1,"(i3.3,i2.2,i2.2)") mtimes(1:3,ntimes)
        v5d%flnm = trim(v5d%flnm)//'_'//trim(charmtime)//'-'
     +  //trim(charmtime1)//'.v5d'
      endif
!
! nflds = number of fields to write to v5d file:
! (note only ZP field are written (no height-only fields))
      nflds = 0
      do i=1,nf
        if ((f(i)%requested.and.associated(f(i)%data).and.
     +    trim(f(i)%vtype)=='ZP').or.trim(f(i)%fname8)== 'Z') 
     +  nflds = nflds+1
      enddo
      if (nflds > MAXVARS) then
        write(6,"('>>> mkv5d: too many fields for vis5d file: nflds=',
     +    i3)") nflds
        write(6,"('    Need to increase MAXVARS (currently MAXVARS=',
     +    i2)") MAXVARS
        ier = 1
        return
      endif
!
! allocate local space for per field info names, nlevs, lowlevs:
! (lowlevs is not used at this time, but may be useful later when
!  combining ccm/tgcm coupled histories) 
!
      fnames = ' '
      if (.not.allocated(nlevs)) then
        allocate(nlevs(nflds),stat=iier)
        if (iier /= 0) call allocerr(iier,"mkv5d allocate nlevs")
      endif 
      if (.not.allocated(lowlevs)) then
        allocate(lowlevs(nflds),stat=iier)
        if (iier /= 0) call allocerr(iier,"mkv5d allocate lowlevs")
      endif 
      ii = 0
!
! Define field names and number of levels for each var.
! (Use 'U' and 'V' rather than 'UN' and 'VN' to match vis5d default
!  wind names)
!
      do i=1,nf
        if ((f(i)%requested.and.associated(f(i)%data).and.
     +    trim(f(i)%vtype)=='ZP').or.trim(f(i)%fname8)== 'Z') then
          ii = ii+1
          if (trim(f(i)%fname8)=='UN') then
            fnames(ii) = 'U         '
          elseif (trim(f(i)%fname8)=='VN') then
            fnames(ii) = 'V         '
          else
            fnames(ii) = f(i)%fname8
          endif
!
! mkvrange returns array of vertical coords (dim nlevs), optionally 
! limited by v5d%zprange. It is called here only to get the number
! of vertical levels:
!
          allocate(vrange(f(i)%nlev),stat=iier)
          if (iier /= 0) call allocerr(iier,"mkv5d allocate vrange")
          vrange = mkvrange(f(i),v5d%zprange,f(i)%nlev,nlevs(ii))
          deallocate(vrange)
          lowlevs(ii) = 0
        endif 
      enddo
!
! int timestamp(ntimes) is time in HHMMSS format
! int datestamp(ntimes) is date in YYDDD  format
!
      allocate(timestamp(ntimes),stat=iier)
      if (iier /= 0) call allocerr(iier,"mkv5d allocate timestamp")
      allocate(datestamp(ntimes),stat=iier)
      if (iier /= 0) call allocerr(iier,"mkv5d allocate datestamp")
      datestamp(1) = h%iyd	! init to 1st date
      do i=1,ntimes
        timestamp(i) = mtimes(2,i)*10000+mtimes(3,i)*100
        if (i > 1) then		! increment date if necessary
          if (timestamp(i) <= timestamp(i-1)) then
            datestamp(i) = datestamp(i-1)+1
          else
            datestamp(i) = datestamp(i-1)
          endif
        endif
      enddo
      proj_args = 0.
      proj_args(1) = gcmlat(nlat)	! North boundary
      proj_args(2) = gcmlon(nlon)	! West boundary (vis5d uses +West)
      proj_args(3) = dlat
      proj_args(4) = dlon
!
      maxnl = nlevs(1)
      do i=1,nflds
        if (nlevs(i) > maxnl) maxnl = nlevs(i)
      enddo
!
! For v5d%vtype='reg_zp',  vertical=0, equally spaced generic units (zp)
! For v5d%vtype='ireg_mb', vertical=3, unequally spaced millibars
! (vert_args are different for different v5d%vtype)
!
      vert_args = 0.
      if (trim(v5d%vtype)=='reg_zp') then	! vtype is regular zp
        vertical = 0
        vert_args(1) = v5d%zprange(1)		! bottom boundary
        vert_args(2) = dlev			! delta vertical
        write(6,"('mkv5d: vertical=',i2,' vert_args(1) (bottom)=',f8.2, 
     +    ' vert_args(2) (dlev)=',f8.2)") vertical,vert_args(1:2)
      else					! vtype is irregular mb
!
! real function mkpmb(zp,nzp,p0,zprange,idimmb,spval,v5dhts=[optional])
        vertical = 3
        v5d_mbs = mkpmb(gcmlev,npress,p0,v5d%zprange,MAXVERTARGS,
     +    spval,v5dhts=v5d_hts)
!       vert_args = v5d_hts
        vert_args = v5d_mbs
        write(6,"('mkv5d: vertical=',i2,' maxnl=',i2,
     +    ' vert_args(1:maxnl) (mb)=',/(6e13.4))") vertical,maxnl,
     +    vert_args(1:maxnl)
      endif
!
! Create v5d file and write header:
      stat = v5dcreate(v5d%flnm,ntimes,nflds,nlat,nlon,nlevs,fnames,
     +  timestamp,datestamp,compress,projection,proj_args,vertical,
     +  vert_args)
      if (stat==0) then
        write(6,"('>>> mkv5d: error return from v5dcreate.')")
        ier = 1
        return
      else
        write(6,"('Created v5d file ',a)") trim(v5d%flnm)
        write(6,"('  nflds=',i3,' fnames=',/(8a8))") 
     +    nflds,(fnames(i),i=1,nflds)
      endif
      deallocate(timestamp)
      deallocate(datestamp)
!
! Add units for each var:
      ii = 0
      do i=1,nf
        if ((f(i)%requested.and.associated(f(i)%data).and.
     +    trim(f(i)%vtype)=='ZP').or.trim(f(i)%fname8)== 'Z') then
          ii = ii+1
          stat = v5dsetunits(ii,trim(f(i)%units)//char(0))
          if (stat==0) then
            write(6,"('>>> mkv5d: error return from v5dsetunits',
     +        ' setting units for field ',a,' units=',a)")
     +        f(i)%fname8,f(i)%units
            ier = 1
          endif
        endif 
      enddo
!
! Set low level offset (if any) for each field (lowlevs(nflds)):
! (not used at this time, 9/97)
      stat = v5dsetlowlev(lowlevs)
      if (stat==0) then
        write(6,"('>>> mkv5d: error return from v5dsetlowlev',
     +    ' setting low level offset for each field:',
     +    ' lowlevs(:)=',/(10i3))") lowlevs
        ier = 1
      endif
!
! Check error status:
      if (ier /= 0) return
!
! End itime==1
 100  continue
!
! Add fields to existing v5d file for current time:
! f(i)%data(nlon,nlat,nlev) must be xferred to ftmp(nlat,nlon,nlev)
!   for v5dwrite, with lat dimension reversed to be N->S.
!
      maxnl = nlevs(1)
      do i=1,nflds
        if (nlevs(i) > maxnl) maxnl = nlevs(i)
      enddo
      allocate(ftmp(nlat,nlon,maxnl),stat=iier)
      if (iier /= 0) call allocerr(iier,"mkv5d allocate ftmp")
!
! Field loop:
      nflds = 0
      floop: do ixf=1,nf
        if ((.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or.
     +    trim(f(ixf)%vtype)/='ZP').and.trim(f(ixf)%fname8)/= 'Z')
     +    cycle floop
        nflds = nflds+1
        ftmp = spval	! init to spval
!
! Calculate vertical array of field (vrange), which may be restricted
!   by v5d%zprange (from optional user input v5d_zprange):
!
        allocate(vrange(f(ixf)%nlev),stat=iier)
        if (iier /= 0) call allocerr(iier,"mkv5d allocate vrange")
        vrange = mkvrange(f(ixf),v5d%zprange,f(ixf)%nlev,nlev)
!
! Define field over its optionally restricted vertical range:
!
        do k=1,nlevs(nflds)
          kk = ixfind(gcmlev,npress,vrange(k),f(ixf)%dlev)
          if (kk <= 0) write(6,"('>>> mkv5d: error finding ',
     +      'vrange(k)=',f8.2,' in gcmlev=',/(9f8.2))") vrange(k),
     +      gcmlev
!
! f(ixf)%vregion(2) is bottom and top of "range of interest" for field.
!   vregion was optionally set for selected fields by setflev in getflds.f.
!   (for fields in which vregion is not used, vregion(1:2)==spval)
! vrange(f(ixf)%nlev) is array of total vertical range of the field,
!   which has been optionally restricted by v5d%zprange.
! If current level vrange(k) is within desired vregion, define
!   field, otherwise leave spval (as init above).
!
          if (f(ixf)%vregion(1) >= f(ixf)%vregion(2).or.	! no vregion
     +       (f(ixf)%vregion(1) <  f(ixf)%vregion(2).and.	! vregion
     +        vrange(k) >= f(ixf)%vregion(1).and.
     +        vrange(k) <= f(ixf)%vregion(2))) then
            do j=1,nlat
!
! Optionally rotate data to center local time v5d%fix_slt (user input):
              if (v5d%fix_slt/=spval) then
                flon = rotslt(f(ixf)%data(:,nlat-j+1,kk),v5d%fix_slt,
     +            h%ut,gcmlon,nlon,dlon,spval)
                ftmp(j,:,k) = flon(:)
              else
                ftmp(j,:,k) = f(ixf)%data(:,nlat-j+1,kk)
              endif	! v5d%fix_slt
!
! Multiply by scale factor, if set:
              if (f(ixf)%scalefac/=1.)
     +          where(ftmp(j,:,k)/=spval) 
     +            ftmp(j,:,k)=ftmp(j,:,k)*f(ixf)%scalefac
            enddo	! j=1,nlat
          endif		! vregion
        enddo		! k=1,nlevs(ixf)
        stat = v5dwrite(itime,nflds,ftmp)
        if (stat==0) then
          write(6,"('>>> mkv5d: error return from v5dwrite: itime=',
     +      i2,' ixf=',i2)") itime,ixf
          ier = 1
          return
        else
          if (f(ixf)%scalefac==1.) then
            write(6,"('Wrote field ',a,' to vis5d file ',a,
     +        ' for itime=',i3)") f(ixf)%fname8,trim(v5d%flnm),itime
          else
            write(6,"('Wrote field ',a,' to vis5d file ',a,
     +        ' for itime=',i3,/'  (scale factor for ',a,' = ',
     +        1pe12.4,')')") 
     +        f(ixf)%fname8,trim(v5d%flnm),itime,
     +        f(ixf)%fname8,f(ixf)%scalefac
          endif
        endif
        deallocate(vrange)
      enddo floop ! ixf=1,nf
      deallocate(ftmp)
!
! Close v5d file if last time:
      if (itime==ntimes) then
        stat = v5dclose()
        if (stat==0) then
          write(6,"('>>> mkv5d: error closing closed v5d file ',a)") 
     +      trim(v5d%flnm)
          ier = 1
          return
        else
          write(6,"('Closed v5d file ',a)") trim(v5d%flnm)
        endif
        if (allocated(nlevs)) deallocate(nlevs)
        if (allocated(lowlevs)) deallocate(lowlevs)
      endif
      return
      end subroutine mkv5d
!-------------------------------------------------------------------
      function mkvrange(f,zprange,mxlev,nlev)
!
! Return array of vertical coords for field f, optionally restricted
!   by zprange(2) (bot,top):
! Also return nlev as number of defined elements in returned array.
! (because this function returns an array, the calling routine must
!  include an interface block for this function)
!
      use fields,only: field
      use proc,only: spval
      implicit none
!
      real :: mkvrange(mxlev)	! function result
!
! Args:
      integer,intent(in) :: mxlev
      type(field),intent(in) :: f
      real,intent(in) :: zprange(2)
      integer,intent(out) :: nlev
!
! Locals:
      integer :: i
      real :: blev,tlev,lev,mxrange(mxlev)
      integer,external :: ixfind
!
      blev = f%lev(1)
      tlev = blev+(f%nlev-1)*f%dlev
      do i=1,mxlev
        mxrange(i) = blev+(i-1)*f%dlev
      enddo
      if (zprange(1) < zprange(2)) then	! zprange is in effect
        if (zprange(1) > blev .and. zprange(1) < tlev) then
          i = ixfind(mxrange,mxlev,zprange(1),f%dlev)
          if (i>0) blev = mxrange(i)
        endif
        if (zprange(2) > blev .and. zprange(2) < tlev) then
          i = ixfind(mxrange,mxlev,zprange(2),f%dlev)
          if (i>0) tlev = mxrange(i)
        endif
      endif
      mkvrange = spval
      nlev = 0 
      do
        nlev = nlev+1
        lev = blev+(nlev-1)*f%dlev
        if (lev > tlev) then
          nlev = nlev-1
          exit
        endif
        mkvrange(nlev) = lev
      enddo
      return
      end function mkvrange
      end module mk_v5d
#endif