!
      program tgcmproc
!
! tgcm post-model processor.
!
      use proc
      use input
      use fields 
      use hist 
      use get_flds 
      use mk_lats 
      use mk_lons 
      use mk_maps 
      use mk_xylocs
      use mk_utvert
      use mk_utlat
      use mk_utlon
      use mk_xyut
      use mk_satut
      use mk_v5d
      use mkcdf_module
      use plt,only: setplt,igks_cgm,igks_ps,igks_x11,
     +                     iwk_cgm, iwk_ps, iwk_x11
      implicit none
!
! Locals:
      integer :: 
     +  it,		! model time loop index
     +  i,j,k,		! msc loop indices
     +  ier,		! error stat
     +  appid,		! application id for ncarg
     +  nflds,		! number of fields to process
     +  iprhist,	! print flag to getflds
     +  ivolp=-1,	! index to histvols of current hist (pert if diffs)
     +  ivolc=-1,	! index to histvols of current hist (cntr if diffs)
     +  ivolp_prev=-1,	! index to histvols of previous hist (pert if diffs)
     +  ivolc_prev=-1,	! index to histvols of previous hist (cntr if diffs)
     +  lup=0,		! logical unit for hist (pert if diffs)
     +  luc=0		! logical unit for hist (cntr if diffs)
      type(history) :: h,h_cntr
      type (field), allocatable,save :: flds_prev(:)
      character(len=120) :: flnm_satcdf	! satellite netcdf file name
      real :: dum
      logical :: isopen
!
! Declarations for timing info:
      integer :: date_time(8)
      character :: date*8, time*10, zone*5
!
! Externals:
      integer,external :: nextlu,mssrcp
!
! Print starting date and time:
!
      call date_and_time(date,time,zone,date_time)
      write(6,"(/'tgcmproc: begin execution at ',
     +  i2.2,'/',i2.2,'/',i4,' ',i2.2,':',i2.2,':',i2.2)")
     +  date_time(2),date_time(3),date_time(1),	! mo,day,yr
     +  (date_time(i),i=5,7)			! hr,min,sec
      write(6,"('tgcmproc_f90 version ',a)") trim(proc_version)
!
! Initialize some stuff:
!
      pi=4.*atan(1.)	! in proc.f
      call set_ohalt
!
! Define geographic grid (note gcmlev is defined after gethist):
!
      do i=1,nlon
        gcmlon(i) = glon1+(i-1)*dlon
      enddo
      do j=1,nlat
        gcmlat(j) = glat1+(j-1)*dlat
      enddo
!
! Get user input parameters:
!
      call getcwd(cwd) ! get exec dir
      call getinp
!
! Define attributes of known fields (no data):
!
      call fset_known
      if (iprint_fknown > 0) then
        write(6,"(/'tgcmproc: ',i4,' known fields are defined:')")
     +    nfknown
        call printfields_table(flds_known,nfknown) 
      endif
!
! Allocate fields structure array for requested fields. 
! (data components are not allocated by allocf)
!
      call allocf
!
! Initialize HLU package, create application object, and needed
! workstations for cgm and/or ps: 
!
      if (iplot > 0) then
        call setplt

        write(6,"('tgcmproc: outplt=',3a8)") outplt

        if (len_trim(outplt(1))==0.and.len_trim(outplt(2))==0.and.
     +      len_trim(outplt(3))==0) iwk_cgm = 1
        do i=1,3
          if (trim(outplt(i))=='cgm') iwk_cgm = 1
          if (trim(outplt(i))=='ps' ) iwk_ps = 1
          if (trim(outplt(i))=='x11'.or.trim(outplt(i))=='X11') 
     +      iwk_x11 = 1
        enddo

        call opnwrk(appid,
     +    flnm_cgm,iwk_cgm,igks_cgm,
     +    flnm_ps, iwk_ps, igks_ps, psmode,
     +             iwk_x11,igks_x11,icolor)

!       call opnwrk(appid,iwk_ps,flnm_ps)
      endif
!
! History (model time) loop:
!
      iprhist = 1
      lup = nextlu()
      if (diffs) luc = nextlu()
      time_loop: do it=1,ntimes
!
! Get history and fields:
!
        nflds = size(flds)
!
! Update fields at previous time for satut:
        if (ipltsatut > 0.and.it /= 1) flds_prev=flds
!
        call getflds(histvols,nvols,mxdiskvols,ivolp,lup,mtimes(1:3,it),
     +    tmpdir,flds,nflds,h,modelhts,iden,ionvel,0,iprhist)
!
! Allocate previous-time fields for satut:
! (6/6/05 btf: added allocation of flds_prev(i)%data)
        if (ipltsatut > 0.and.it == 1) then
          allocate(flds_prev(nflds),stat=ier)
          if (ier /= 0) call allocerr(ier,"allocating flds_prev")
          do i=1,nflds
            allocate(flds_prev(i)%data(h%nlon-1,h%nlat,flds(i)%nlev),
     |        stat=ier)
            if (ier /= 0) call allocerr(ier,"allocating flds_prev data")
          enddo
          flds_prev = flds
        endif
!
! If history was not found, then start search for next history if
! isearch > 0. Otherwise exit history loop:
!
        if (ivolp==0) then		! history not found
          if (isearch > 0) then
            ivolp = -abs(ivolp_prev)
            cycle time_loop
          endif
          exit time_loop
        elseif (ivolp < 0) then
          write(6,"('>>> tgcmproc: fatal error from getflds.')")
          exit time_loop
        else
          ivolp_prev = ivolp
        endif
!
! Get control fields if doing diffs:
        if (diffs) then
            call getflds(histvols_cntr,nvols_cntr,mxdiskvols,ivolc,luc,
     +      mtimes_cntr(1:3,it),tmpdir,flds_cntr,nflds,h_cntr,
     +      modelhts,iden,ionvel,1,iprhist)
          if (ivolc == 0) then		! history not found
            if (isearch > 0) then
              ivolc = -abs(ivolc_prev)
              cycle time_loop
            endif
            exit time_loop
          elseif (ivolc < 0) then
            write(6,"('>>> tgcmproc: fatal error from getflds.')")
            exit time_loop
          else
            ivolc_prev = ivolc
          endif
!
! Check some diffs stuff and report to stdout:
          call checkdiffs(flds,flds_cntr,nflds,h,h_cntr,iprhist)
        endif	! diffs
!
! Define model pressure coordinates for this history:
! (defines npress, zpb, zpt, and gcmlev(npress) for midpoints, and 
!  equivalent variables for interfaces, see module proc)
!
        call setlev(h,it)
!
! Report 3d min,max if requested:
!
        nflds = size(flds)
        if (iprint_fldminmax>0) then
          if (.not.diffs) then
            write(6,"('Global 3d min,max:')")
            call fldminmax(flds)
          else
            write(6,"('Global 3d min,max for perturbed case:')")
            call fldminmax(flds)
            write(6,"('Global 3d min,max for control case:')")
            call fldminmax(flds_cntr)
            write(6,"('Global 3d min,max differences:')")
            call diffs3d(flds,flds_cntr,nflds)
          endif
        endif
!
! Set up f%cmin,cmax,cint if requested by user (input fmnmxint):
! Also set f%scalefac from optional user provided fscale.
!
        call setfmnmxint(flds,nflds,fmnmxint,fscale,mxfproc,iprhist)
!
! Define user-provided field units: 
!
        call setunits(flds,nflds,cunits,mxfproc,iprhist)
!
! Make netcdf file for current history:
! (1/98: not yet avail on sgi)
        if (len_trim(sendcdf) > 0 .or. len_trim(flnm_cdf) > 0) then

!         write(6,"('tgcmproc calling mkcdf: it=',i3,' nflds=',i3)")
!    |      it,nflds

          call mkcdf(flds,flds_cntr,nflds,h,flnm_cdf,it,ier)
          if (ier /= 0) then
            write(6,"('>>> tgcmproc: error return from mkcdf -- ',
     +        'will not make netcdf file.')")
            sendcdf = ' '
          endif
        endif
!
! Make vis5d file for current history:
! (1/98: not yet avail on sgi)
        if (len_trim(v5d%sendv5d)>0.or.len_trim(v5d%sendmsv5d)>0) then
#if defined(SUN) || defined (LINUX)
          write(6,"('>>> v5d not available on the Sun.')")
          stop 'v5d on Sun'
#else
          call mkv5d(flds,nflds,h,it,mtimes,ier)
          if (ier /= 0) then
            write(6,"('>>> tgcmproc: error return from mkv5d -- ',
     +        'will not make vis5d file.')")
            v5d%flnm = ' '
          endif
#endif
        endif
!
! Activate workstation(s) for LLU calls:
! (note this is inside time_loop because of deactivation before
!  hlu calls):
!
        if (iplot > 0) then
          if (igks_cgm > 0) call gacwk(igks_cgm)
          if (igks_ps  > 0) call gacwk(igks_ps )
          if (igks_x11 > 0) call gacwk(igks_x11)
        endif
!
! Qadglb makes "Quick-And-Dirty" global contours
!
!       call qadglb(flds,nflds,h)
!
! Make maps:
        if (ipltmaps > 0)call mkmaps(flds,flds_cntr,nflds,h,h_cntr)
!
! Make lon slices:
! (diffs not available for amp/phase)
        if (ipltlon > 0) then
          call mklons(flds,flds_cntr,nflds,h,h_cntr,it)
          if (amphase>0) call mkampha(flds,nflds,h,h_cntr)
!         if (istream>0) call mkstrm(flds,nflds,h,h_cntr)
        endif
!
! Make lat slices:
        if (ipltlat > 0) call mklats(flds,flds_cntr,nflds,h,h_cntr)
!
! Make vertical profiles at selected locations:
!
        if (ipltxyloc > 0) call mkxyloc(flds,flds_cntr,nflds,h,h_cntr)
!
! Contour ut vs zp/ht at selected locations:
!
        if (ipltutvert > 0) 
     +    call mkutvert(flds,flds_cntr,nflds,h,h_cntr,it)
!
! Contour ut vs zp/ht along satellite track:
!
        if (ipltsatut > 0) then
          call mksatut(flds,flds_cntr,flds_prev,nflds,
     +                 h,h_cntr,flnm_satcdf,it)
        endif
!
! Contour ut vs latitude at selected zp/ht, lons:
!
        if (ipltutlat > 0) 
     +    call mkutlat(flds,flds_cntr,nflds,h,h_cntr,it)
!
! Contour ut vs longitude at selected zp/ht, lats:
!
        if (ipltutlon > 0) 
     +    call mkutlon(flds,flds_cntr,nflds,h,h_cntr,it)
!
! Line plot ut vs field at selected locations and levels:
!
        if (ipltxyut > 0)
     |    call mkxyut(flds,flds_cntr,nflds,h,h_cntr,it)
!
! Release gcmlev (was allocated above after getflds):
!
        deallocate(gcmlev,stat=ier)
        if (iplot > 0) then
          if (igks_cgm > 0) call gdawk(igks_cgm)
          if (igks_ps  > 0) call gdawk(igks_ps)
          if (igks_x11 > 0) call gdawk(igks_x11)
        endif
!
! End history (model time) loop:
        iprhist = 0
!
! Release allocated memory from field structure data 
! pointers f(i)%data. These were allocated by allocfdat,
! called by getflds.
!
        call deallocfdat(flds,nflds)
      enddo time_loop
!      
! Retire HLU package:
!
      if (iplot > 0) then
        if (iwk_cgm > 0) call NhlFDestroy(iwk_cgm,ier)
        if (iwk_ps > 0) call NhlFDestroy(iwk_ps,ier)
        call NhlFDestroy(appid,ier)
        call NhlFClose
      endif
!
! Send netcdf file to remote:
      if (len_trim(sendcdf) > 0) 
     +  call scpfile(0,flnm_cdf,sendcdf)
!
! Send lon slice netcdf file to remote:
      if (len_trim(sendcdf_lons) > 0) 
     |  call scpfile(0,flnm_cdf_lons,sendcdf_lons)
!
! Send xyut netcdf file to remote.
      if (len_trim(sendcdf_xyut) > 0) 
     |  call scpfile(0,flnm_cdf_xyut,sendcdf_xyut)
!
! Send satellite netcdf file to remote:
      if (len_trim(sendsat) > 0) 
     +  call scpfile(0,flnm_satcdf,sendsat)
!
! Send vis5d file to mss and/or remote machine:
      if (len_trim(v5d%flnm) > 0) then
        if (len_trim(v5d%sendmsv5d) > 0) then
          ier=mssrcp(' ',trim(v5d%flnm),'mss:'//trim(v5d%sendmsv5d))
          if (ier/=0)
     |      write(6,"('>>> tgcmproc: Error copying disk file ',a,
     |        ' to mss path ',a)") trim(v5d%flnm),trim(v5d%sendmsv5d)
        endif
        if (len_trim(v5d%sendv5d) > 0) 
     +    call scpfile(0,v5d%flnm,v5d%sendv5d)
      endif
!
! Send cgm file to remote:
      write(6,"(' ')")
      if (len_trim(sendcgm) > 0.and.iwk_cgm > 0) 
     +  call scpfile(0,flnm_cgm,sendcgm)
!
! Send ps files to remote:
      if (len_trim(sendps) > 0.and.iwk_ps > 0) 
     +  call scpfile(0,flnm_ps,sendps)
!
! Close dat file:
      if (len_trim(flnm_dat) > 0) then
        inquire(ludat,opened=isopen)
        if (isopen) then
          write(6,"('tgcmproc: closing ascii data file ',a)") 
     |      trim(flnm_dat)
          close(ludat)
        endif
      endif
!
! Send ascii data file to remote:
      if (len_trim(senddat) > 0)
     +  call scpfile(ludat,flnm_dat,senddat)
!
! Send ascii data file in Bougher format to remote:
      if (len_trim(senddat_bf) > 0)
     +  call scpfile(ludat_bf,flnm_bf_dat,senddat_bf)
!
! Send xdr data file to remote:
      if (len_trim(sendxdr) > 0.or.len_trim(sendms_xdr) > 0) then

!       write(6,"('before wrxdr: len_trim(flnm_xdr)=',i3,
!    |    ' trim(flnm_xdr)=$',a,'$')") len_trim(flnm_xdr),
!    |    trim(flnm_xdr)

        call wrxdr(trim(flnm_xdr),
     |    dum,dum,dum,dum,dum,dum,dum,dum,dum,dum,dum,dum,1)

!       write(6,"('before scpfile: len_trim(flnm_xdr)=',i3,' flnm_xdr=',
!    |    '$',a,'$')") len_trim(flnm_xdr),trim(flnm_xdr)
!       write(6,"('before scpfile: len_trim(sendxdr)=',i3,' sendxdr=',
!    |    '$',a,'$')") len_trim(sendxdr),trim(sendxdr)

        if (len_trim(sendxdr) > 0) call scpfile(0,flnm_xdr,sendxdr)
        if (len_trim(sendms_xdr) > 0) then
          ier=mssrcp(' ',trim(flnm_xdr),'mss:'//trim(sendms_xdr))
          if (ier/=0)
     |      write(6,"('>>> tgcmproc: Error copying disk file ',a,
     |        ' to mss path ',a)") trim(flnm_xdr),trim(sendms_xdr)
        endif
      endif
!
! Print ending date and time:
!
      call date_and_time(date,time,zone,date_time)
      write(6,"(/'tgcmproc: end execution at ',
     |  i2.2,'/',i2.2,'/',i4,' ',i2.2,':',i2.2,':',i2.2)")
     |  date_time(2),date_time(3),date_time(1),	! mo,day,yr
     |  (date_time(i),i=5,7)			! hr,min,sec
      stop 'done'
      end
!-------------------------------------------------------------------
      subroutine setlev(h,it)
      use proc
      use hist,only: history
      implicit none
!
! Define model vertical pressure resolution for this history:
! (npress, zpb, zpt, gcmlev(:) in module proc)
! This is called from tgcmproc after rdhist in history loop. gcmlev
!  is deallocated by tgcmproc at end of each iter of history loop.
!
! Args:
      integer,intent(in) :: it
      type(history),intent(in) :: h
! Locals:
      integer :: k,ier
!
      npress = h%nzp
!
! Midpoints:
      zpb = h%zpb
      zpt = h%zpt
      dlev = (zpt-zpb)/float(npress-1)
      if (allocated(gcmlev)) deallocate(gcmlev)
      allocate(gcmlev(npress),stat=ier)
      if (ier /= 0) call allocerr(ier,"setlev allocating gcmlev")
      do k=1,npress
        gcmlev(k) = zpb+(k-1)*dlev
      enddo
!
! Interfaces:
      zpib = h%zpib
      zpit = h%zpit
      if (allocated(gcmilev)) deallocate(gcmilev)
      allocate(gcmilev(npress),stat=ier)
      if (ier /= 0) call allocerr(ier,"setlev allocating gcmilev")
      do k=1,npress
        gcmilev(k) = zpib+(k-1)*dlev
      enddo

!     if (it==1) then
!       write(6,"('setlev interfaces: npress=',i4,' zpib=',f8.3,
!    |    ' zpit=',f8.3,' dlev=',f8.3,' gcmilev=',/,(9f8.3))") 
!    |    npress,zpib,zpit,dlev,gcmilev
!       if (h%isnew)
!    |    write(6,"('setlev midpoints: npress=',i4,' zpb=',f8.3,
!    |      ' zpt=',f8.3,' dlev=',f8.3,' gcmlev=',/,(9f8.3))") 
!    |      npress,zpb,zpt,dlev,gcmlev
!     endif

      if (gcmlev(npress).ne.zpt) then
        write(6,"('>>> WARNING: gcmlev(npress) != zpt')")
        write(6,"('    npress=',i3,' zpb,zpt=',2f14.10)")npress,zpb,zpt
        write(6,"('    gcmlev=',/(5f14.10))") gcmlev
        stop "gcmlev"
      endif
      if (gcmilev(npress).ne.zpit) then
        write(6,"('>>> WARNING: gcmilev(npress) != zpit')")
        write(6,"('    npress=',i3,' zpib,zpit=',2f14.10)")npress,zpib,
     |    zpit
        write(6,"('    gcmilev=',/(5f14.10))") gcmilev
        stop "gcmlev"
      endif
      return
      end subroutine setlev
!-------------------------------------------------------------------
      subroutine setfmnmxint(f,nf,fmnmxint,fscale,mxf,iprint)
!
! If user provided fixed min,max,interval (by setting fmnmxint(:,ixf)),
!   transfer this info into f(ixf)%cmin,cmax,cint
!   (fmnmxint was partially validated in getinp (input.f))
! Also set scale factors f(ixf)%scalefac if fscale(2,mxf) was set by user
!
      use proc,only: spval
      use input,only: cmnmxint,cscale
      use fields,only: field
      implicit none
!
! Args:
      integer,intent(in) :: nf,mxf,iprint
      type(field),intent(inout) :: f(nf)
!     real,intent(in) :: fmnmxint(4,mxf),fscale(2,mxf)
      real,intent(in) :: fscale(2,mxf)
      real(kind=8),intent(in) :: fmnmxint(4,mxf)
!
! Locals:
      integer :: ixf,if
!
! Externals:
      character(len=16),external :: float_to_str,float8_to_str
!
      user_loop: do ixf=1,mxf			! loop through all fmnmxint
        if (fmnmxint(1,ixf)/=spval) then	! user provided limits
          fields_loop: do if=1,nf		! find which field
            if (trim(float8_to_str(fmnmxint(1,ixf)))==
     |        trim(f(if)%fname8)) then
              f(if)%cmin = fmnmxint(2,ixf)
              f(if)%cmax = fmnmxint(3,ixf)
              f(if)%cint = fmnmxint(4,ixf)
              cycle user_loop
            endif
          enddo fields_loop
          if (iprint > 0) then
            write(6,"('>>> WARNING: could not find field ',a,' to set',
     +        ' fmnmxint (cmin,cmax,cint)')") fmnmxint(1,ixf)
            write(6,"('    (This field was either not available ',
     +        'or not requested.)')")
          endif
        endif
      enddo user_loop
!
! Also check for cmnmxint:
      user_loop1: do ixf=1,mxf			! loop through all fmnmxint
        if (len_trim(cmnmxint(1,ixf))>0) then	! user provided limits
          fields_loop1: do if=1,nf		! find which field
            if (trim(cmnmxint(1,ixf))==trim(f(if)%fname8)) then
              read(cmnmxint(2,ixf),fmt=*) f(if)%cmin
              read(cmnmxint(3,ixf),fmt=*) f(if)%cmax
              read(cmnmxint(4,ixf),fmt=*) f(if)%cint
              cycle user_loop1
            endif
          enddo fields_loop1
          if (iprint > 0) then
            write(6,"('>>> WARNING: could not find field ',a,' to set',
     |        ' cmnmxint (cmin,cmax,cint)')") cmnmxint(1,ixf)
            write(6,"('    (This field was either not available ',
     |        'or not requested.)')")
          endif
        endif
      enddo user_loop1
!
! Also set scale factors:
      user_loop2: do ixf=1,mxf			! loop through all fmnmxint
        if (fscale(1,ixf)/=spval) then	! user provided limits
          fields_loop2: do if=1,nf		! find which field
            if (trim(float_to_str(fscale(1,ixf)))==trim(f(if)%fname8))
     +        then
              f(if)%scalefac = fscale(2,ixf)
              cycle user_loop2
            endif
          enddo fields_loop2
          if (iprint > 0) then
            write(6,"('>>> WARNING: could not find field ',a,' to set',
     +        ' fscale (scale factor)')") fscale(1,ixf)
            write(6,"('    (This field was either not available ',
     +        'or not requested.)')")
          endif
        endif
      enddo user_loop2
!
! Check for cscale:
      user_loop3: do ixf=1,mxf			! loop through all fmnmxint
        if (len_trim(cscale(1,ixf))>0) then	! user provided limits
          fields_loop3: do if=1,nf		! find which field
            if (trim(cscale(1,ixf))==trim(f(if)%fname8)) then
              read(cscale(2,ixf),fmt=*) f(if)%scalefac
              cycle user_loop3
            endif
          enddo fields_loop3
          if (iprint > 0) then
            write(6,"('>>> WARNING: could not find field ',a,' to set',
     +        ' cscale (scale factor)')") cscale(1,ixf)
            write(6,"('    (This field was either not available ',
     +        'or not requested.)')")
          endif
        endif
      enddo user_loop3
      return
      end subroutine setfmnmxint
!-------------------------------------------------------------------
      subroutine setunits(f,nf,cunits,mxf,iprint)
      use fields,only: field
!
! Set any user-provided units for fields (cunits(2,mxf) from input):
!
! Args:
      integer,intent(in) :: nf,mxf,iprint
      type(field),intent(inout) :: f(nf)
      character(len=*),intent(in) :: cunits(2,mxf)
!
! Locals:
      integer :: ixf,if
!
! Externals:
      character(len=16),external :: float_to_str
!
      user_loop: do ixf=1,mxf			! loop through all cunits
        if (len_trim(cunits(1,ixf))>0) then     ! user provided units
          fields_loop: do if=1,nf		! find which field
            if (cunits(1,ixf)==f(if)%fname8) then
              f(if)%units = cunits(2,ixf)
              if (iprint > 0) 
     +          write(6,"('Set user provided units: field = ',
     +            a8,' units = ',a)") cunits(:,ixf) 
              cycle user_loop
            endif
          enddo fields_loop
          if (iprint > 0) then
            write(6,"('>>> WARNING: could not find field ',a,' to set',
     +        ' units (from cunits)')") cunits(1,ixf)
            write(6,"('    (This field was either not available ',
     +        'or not requested.)')")
          endif
        endif
      enddo user_loop

      end subroutine setunits
!-------------------------------------------------------------------
      subroutine checkdiffs(fpert,fcntr,nf,hpert,hcntr,iprhist)
!
! Check that ut's are same for perturbed and control cases,
!   deallocate any fields not available on both histories,
!   set field difftype, and report to stdout.
!
      use fields,only: field
      use hist,only: history
      use input,only: idifpercent
      use proc,only: ispval
      implicit none
!
! Args:
      integer,intent(in) :: nf,iprhist
      type(field),intent(inout) :: fpert(nf),fcntr(nf)
      type(history),intent(in) :: hpert,hcntr
!
! Locals:
      integer :: i,nalloc
      logical :: palloc,calloc
!
      if (iprhist > 0) then
        write(6,"(/'DIFFERENCE FIELDS RUN:')")
        write(6,"('  Perturbed history: ',a,' (',3i4,')')")
     +    trim(hpert%mssvol),hpert%mtime
        write(6,"('  Control history:   ',a,' (',3i4,')')")
     +    trim(hcntr%mssvol),hcntr%mtime
      endif
!
! Confirm same model times for pert and control:
!
      if (hpert%mtime(2)/=hcntr%mtime(2).or.
     +    hpert%mtime(3)/=hcntr%mtime(3)) then
        write(6,"(/'>>> checkdiffs: model times (ut) of',
     +    ' perturbed and control cases differ:')")
        write(6,"(14x,' perturbed mtime=',3i4,
     +    ', control mtime=',3i4)") hpert%mtime,hcntr%mtime
        stop 'diffs'
      endif
!
! Deallocate any fields that were not available in both
! perturbed and control histories. Stop if no requested
! fields were available on both histories:
!
      nalloc = 0
      do i=1,nf
        palloc = associated(fpert(i)%data)
        calloc = associated(fcntr(i)%data)
        if (palloc.and..not.calloc) then
          if (iprhist > 0) 
     +      write(6,"('>>> WARNING: field ',a,' is not available from ',
     +      'the control history,',/,'  so will not be processed.')")
     +      fpert(i)%fname8
          deallocate(fpert(i)%data)
        elseif (calloc.and..not.palloc) then
          if (iprhist > 0) 
     +      write(6,"('>>> WARNING: field ',a,' is not available from ',
     +      'the perturbed history,',/,'  so will not be processed.')")
     +      fpert(i)%fname8
          deallocate(fcntr(i)%data)
        elseif (palloc.and.fpert(i)%requested.and.
     +          calloc.and.fcntr(i)%requested) then
          nalloc = nalloc+1
!
! Set difftype:
          if (idifpercent==0.or.(idifpercent==ispval.and.
     +        trim(fpert(i)%type)/='DENSITY')) then
            fpert(i)%difftype = "RAW"
            fcntr(i)%difftype = "RAW"
          elseif (idifpercent==1.or.(idifpercent==ispval.and.
     +            trim(fpert(i)%type)=='DENSITY')) then
            fpert(i)%difftype = "PERCENT"
            fcntr(i)%difftype = "PERCENT"
          endif
        endif
!
! Make sure difftype is set (this might include unrequested 
!   dependencies for derived fields)
        if (len_trim(fpert(i)%difftype)==0.or.
     |      len_trim(fcntr(i)%difftype)==0) then
          if (idifpercent==0.or.(idifpercent==ispval.and.
     |        trim(fpert(i)%type)/='DENSITY')) then
            fpert(i)%difftype = "RAW"
            fcntr(i)%difftype = "RAW"
          elseif (idifpercent==1.or.(idifpercent==ispval.and.
     |            trim(fpert(i)%type)=='DENSITY')) then
            fpert(i)%difftype = "PERCENT"
            fcntr(i)%difftype = "PERCENT"
          endif
        endif

!       write(6,"('checkdiffs: i=',i3,' nf=',i3,' field ',a,
!    |    ' idifpercent=',i3,' difftype=',a,' pert.req=',l1,
!    |    ' cntr.req=',l1)")
!    |    i,nf,fpert(i)%fname8,idifpercent,fpert(i)%difftype,
!    |    fpert(i)%requested,fcntr(i)%requested

      enddo
      if (nalloc==0) then
        write(6,"(/'>>> No requested fields were available in both ',
     +    'perturbed and control histories.'/)")
        stop 'diffs'
      endif
!
! idifpercent==ispval (user did not provide override) -> 
!   take percent diffs of species, raw diffs of all other fields.
! idifpercent==0 (user overrode w/ idifpercent<=0) ->
!   take raw diffs of all fields.
! idifpercent==1 (user overrode w/ idifpercent>0) ->
!   take percent diffs of all fields.
!
      if (iprhist > 0) then
        if (idifpercent==ispval) then
          write(6,"('  Will calculate percent diffs of ',
     +      'species, raw diffs of all other fields.')")
        elseif (idifpercent<=0) then
          write(6,"('  Will calculate raw diffs of all fields.')")
        else
          write(6,"('  Will calculate percent diffs of all fields.')")
        endif
        write(6,"(' ')")
      endif
      end subroutine checkdiffs
!-------------------------------------------------------------------
      subroutine diffs3d(pert,cntr,nf)
      use fields,only: field
      implicit none
!
! Calculate and report min,max of 3d differences of allocated 
! fields pert minus cntr:
!
! Args:
      integer,intent(in) :: nf
      type(field),intent(in) :: pert(nf),cntr(nf)
!
! Locals:
      integer :: i,n,ier
      real,allocatable :: diffs(:)
      real :: dmin,dmax
!
      floop: do i=1,nf
!
! Pert and cntr fields data must be associated and requested, 
! field names, data sizes, and difftypes must match:
!
        if (.not.associated(pert(i)%data).or.
     +      .not.associated(cntr(i)%data)) cycle floop
        if (.not.pert(i)%requested.or..not.cntr(i)%requested)
     +    cycle floop
        if (pert(i)%fname8 /= cntr(i)%fname8) then
          write(6,"('WARNING diffs3d: field names do not match:',
     +      ' i=',i2,' pert name=',a,' cntr name=',a)")
     +      i,pert(i)%fname8,cntr(i)%fname8
          cycle floop
        endif
        n = size(pert(i)%data)
        if (n /= size(cntr(i)%data)) then
          write(6,"('WARNING diffs3d: field data sizes do not match:',
     +      ' field ',a,' i=',i2,' pert size=',i8,' cntr size=',i8)")
     +      pert(i)%fname8,i,n,size(cntr(i)%data)
          cycle floop 
        endif
        if (pert(i)%difftype /= cntr(i)%difftype) then
          write(6,"('WARNING diffs3d: field ',a,' difftypes do not ',
     +      'match: pert difftype=',a,' cntr difftype=',a)")
     +      pert(i)%fname8,pert(i)%difftype,cntr(i)%difftype
          cycle floop 
        endif
        allocate(diffs(n),stat=ier)
        if (ier /= 0) call allocerr(ier,"diffs3d allocating diffs")
        call mkdiffs(pert(i)%data,cntr(i)%data,diffs,n,pert(i)%difftype)
        call fminmax(diffs,n,dmin,dmax)
        if (pert(i)%difftype=="PERCENT ") then 
          write(6,"('Field ',a,' 3d min,max percent differences=',
     +      2e12.4)") pert(i)%fname8,dmin,dmax
        elseif (pert(i)%difftype=="RAW     ") then 
          write(6,"('Field ',a,' 3d min,max raw differences=',4x,
     +      2e12.4)") pert(i)%fname8,dmin,dmax
        else ! should not happen
          write(6,"('Field ',a,' 3d min,max differences=',2e12.4)")
     +      pert(i)%fname8,dmin,dmax
        endif
        deallocate(diffs)
      enddo floop
      end subroutine diffs3d
!-----------------------------------------------------------------------
      subroutine deallocfdat(f,nf)
      use fields,only: field
      implicit none
!
! Args:
      integer,intent(in) :: nf
      type (field),intent(inout) :: f(nf)
!
! Local:
      integer :: i,istat,nalloc
!
      nalloc = 0
      do i=1,nf
        if (associated(f(i)%data)) then
          deallocate(f(i)%data,stat=istat)
          nalloc = nalloc+1
        endif
      enddo
!     write(6,"('Deallocated f(ix)%data for ',i3,
!    |  ' fields.')") nalloc
      end subroutine deallocfdat
