module newhist_module
      use readhist_module,only: history,fminmaxspv,field_3d,
     |  fminmax,mxflds
      use grid_module
      implicit none
      contains
!-----------------------------------------------------------------------
      subroutine newhist(h_old,h_new,ihist)
!
! Define new geographic grid on history structure h_new, 
!   from newgrid structure (grid.F)
!
! Define magnetic grid on h_new from h_old (vertical mag grid may
!   be interpolated from 0.5 to 0.25, but the horizontal magnetic 
!   grid is unchanged from old history)
!
      implicit none
!
! Args:
      type(history),intent(inout)  :: h_old
      type(history),intent(inout) :: h_new 
      integer,intent(in) :: ihist
!
! Local:
      integer :: ier,n
      real :: dlat,dlon,dlev
      character(len=8) :: create_date,create_time
      real :: glat(h_old%nlat),glon(h_old%nlon)
      real :: fmin,fmax
!
! Copy h_old into h_new, excluding pointers:
      call copyhist(h_old,h_new)
!
! Define geographic grid on new history structure:
!
! Set latitudes for new history (see grid.F):
      h_new%nlat = newgrid%nlat
      if (associated(h_new%glat)) deallocate(h_new%glat)
      allocate(h_new%glat(h_new%nlat),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new glat:',
     |    ' nlat=',i3)") h_new%nlat
      h_new%glat(:) = newgrid%glat(:)
      dlat = h_new%glat(2)-h_new%glat(1) ! assume regular spacing
!
! Set longitudes for new history (see grid.F):
      h_new%nlon = newgrid%nlon
      if (associated(h_new%glon)) deallocate(h_new%glon)
      allocate(h_new%glon(h_new%nlon),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new glon:',
     |    ' nlon=',i3)") h_new%nlon
      h_new%glon(:) = newgrid%glon(:)
      dlon = h_new%glon(2)-h_new%glon(1) ! assume regular spacing
!
! Set midpoint levels for new history (see grid.F):
      h_new%nlev = newgrid%nlev
      if (associated(h_new%lev)) deallocate(h_new%lev)
      allocate(h_new%lev(h_new%nlev),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new lev:',
     |    ' nlev=',i3)") h_new%nlev
      h_new%lev(:) = newgrid%zlev(:)
      dlev = h_new%lev(2)-h_new%lev(1) ! assume regular spacing
!
! Set interface levels for new history (see grid.F):
      h_new%nilev = newgrid%nilev
      if (associated(h_new%ilev)) deallocate(h_new%ilev)
      allocate(h_new%ilev(h_new%nilev),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new ilev:',
     |    ' nilev=',i3)") h_new%nilev
      h_new%ilev(:) = newgrid%zilev(:)
      dlev = h_new%ilev(2)-h_new%ilev(1) ! assume regular spacing
!
! Define magnetic grid on new history structure:
!
! Note: This code does not change the magnetic grid, 
!       it is only echoed from the input history.
!
! mlon,mlat,mlev,    ! magnetic grid dimensions
! imlev,             ! number of mag interface levels
! real,pointer :: gmlat(:),gmlon(:)     ! geomagnetic lat lon grid
! real,pointer :: gmlev(:),gimlev(:)    ! midpoint and interface mag levels
!
! Set mag latitudes for new history (see grid.F):
      h_new%mlat = newgrid%mlat
      if (associated(h_new%gmlat)) deallocate(h_new%gmlat)
      allocate(h_new%gmlat(h_new%mlat),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new gmlat:',
     |    ' mlat=',i3)") h_new%mlat
      h_new%gmlat(:) = newgrid%gmlat(:)
!
! Set mag longitudes for new history (see grid.F):
      h_new%mlon = newgrid%mlon
      if (associated(h_new%gmlon)) deallocate(h_new%gmlon)
      allocate(h_new%gmlon(h_new%mlon),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new gmlon:',
     |    ' mlon=',i3)") h_new%mlon
      h_new%gmlon(:) = newgrid%gmlon(:)
!
! Set midpoint mag midpoint levels for new history (see grid.F):
      h_new%mlev = newgrid%mlev
      if (associated(h_new%gmlev)) deallocate(h_new%gmlev)
      allocate(h_new%gmlev(h_new%mlev),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new gmlev:',
     |    ' gmlev=',i3)") h_new%gmlev
      h_new%gmlev(:) = newgrid%zmlev(:)
!
! Set midpoint mag interface levels for new history (see grid.F):
      h_new%imlev = newgrid%imlev
      if (associated(h_new%gimlev)) deallocate(h_new%gimlev)
      allocate(h_new%gimlev(h_new%imlev),stat=ier)
      if (ier /= 0)
     |  write(6,"('newhist: error allocating new gimlev:',
     |    ' gimlev=',i3)") h_new%gimlev
      h_new%gimlev(:) = newgrid%zimlev(:)
!
! Report to stdout:
      write(6,"(/,'newhist: nlat=',i3,' dlat=',f6.2,' glat=',/,
     |  (9f8.2))") h_new%nlat,dlat,h_new%glat
      write(6,"(/,'newhist: nlon=',i3,' dlon=',f6.2,' glon=',/,
     |  (9f8.2))") h_new%nlon,dlon,h_new%glon
      write(6,"(/,'newhist: nlev=',i3,' dlev=',f6.2,' midpoint levs=',/,
     |  (9f8.3))") h_new%nlev,dlev,h_new%lev
      write(6,"(/,'newhist: nilev=',i3,' dlev=',f6.2,
     |  ' interface levs=',/,(9f8.3))") h_new%nilev,dlev,h_new%ilev

      write(6,"(/,'newhist: mlat=',i3,' gmlat =',/,(9f8.2))") 
     |  h_new%mlat, h_new%gmlat
      write(6,"(/,'newhist: mlon=',i3,' gmlon =',/,(9f8.2))") 
     |  h_new%mlon, h_new%gmlon
      write(6,"(/,'newhist: mlev=',i3,' gmlev =',/,(9f8.3))") 
     |  h_new%mlev, h_new%gmlev
      write(6,"(/,'newhist: imlev=',i3,' gimlev=',/,(9f8.3))")
     |  h_new%imlev, h_new%gimlev

      end subroutine newhist
!-----------------------------------------------------------------------
      subroutine horizontal_interp(h_old,h_new)
      implicit none
!
! Interpolate fields to new horizontal resolution 
! (vertical res is same for both input and output)
! Regridding routine rgrd2 is from SCD's regridpack lib,
!   http://www.scd.ucar.edu/softlib/REGRIDPACK.html
!
! Args:
      type(history),intent(in)  :: h_old
      type(history),intent(inout) :: h_new 
!
! Local:
      integer :: n,ier,lw,liw,isign,iprint,intpol(2)
      integer,allocatable :: iw(:)
      integer :: ilen,ilen1
      real,allocatable :: w(:)
      real :: tmplon_old(h_old%nlon+4),tmpold(h_old%nlon+4,2),
     |  tmplon_olddata(h_old%nlon+4),tmplat2(2),dlat_old
      real :: fmin,fmax
!
! Allocate workspace for regrid routine rgrd2.f:
      lw = h_new%nlon+h_new%nlat+2*h_new%nlon
      allocate(w(lw),stat=ier)
      if (ier /= 0) write(6,"('>>> horizontal_interp: error allocating',
     |  ' w(lw): lw=',i6,' ier=',i4)") lw,ier
      liw = h_new%nlon+h_new%nlat
      allocate(iw(liw),stat=ier)
      if (ier /= 0) write(6,"('>>> horzontal_interp: error allocating',
     |  ' iw(liw): liw=',i6,' ier=',i4)") liw,ier
      intpol(:) = 1 ! linear (not cubic) interp in both dimensions
!
! Transform to new horizontal:
      do n=1,h_new%nf2d_geo
!
! Regrid does bilinear interpolation (source from local ncar lib regrid)
!   subroutine rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)
!   subroutine rgrd2 interpolates the values p(i,j) on the orthogonal
!   grid (x(i),y(j)) for i=1,...,nx and j=1,...,ny onto q(ii,jj) on the
!   orthogonal grid (xx(ii),yy(jj)) for ii=1,...,mx and jj=1,...,my.
!
! subroutine rgrd2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,w,lw,iw,liw,ier)
!
        call rgrd2(
     |    h_old%nlon, h_old%nlat, h_old%glon, h_old%glat,
     |    h_old%f2d_geo(n)%data(:,:),
     |    h_new%nlon-1, h_new%nlat-2, h_new%glon(1:h_new%nlon-1),
     |    h_new%glat(2:h_new%nlat-1),
     |    h_new%f2d_geo(n)%data(1:h_new%nlon-1,2:h_new%nlat-1),
     |    intpol,w,lw,iw,liw,ier) 
        if (ier /= 0)
     |    write(6,"('>>> WARNING horzontal_interp: error from rgrd2 ',
     |      'interpolating 2d: ier=',i3)") ier
!
! Interpolate between i=1 and i=nlon-1 to get i=nlon:
        h_new%f2d_geo(n)%data(h_new%nlon,:) = 
     |    0.5*(h_new%f2d_geo(n)%data(1,:)+
     |         h_new%f2d_geo(n)%data(h_new%nlon-1,:))
!
! Extrapolate to latitude boundaries:
        h_new%f2d_geo(n)%data(:,1) = 
     |    1.5*h_new%f2d_geo(n)%data(:,2)-
     |    0.5*h_new%f2d_geo(n)%data(:,3)

        h_new%f2d_geo(n)%data(:,h_new%nlat) = 
     |    1.5*h_new%f2d_geo(n)%data(:,h_new%nlat-1)-
     |    0.5*h_new%f2d_geo(n)%data(:,h_new%nlat-2)
      enddo

      end subroutine horizontal_interp
!-----------------------------------------------------------------------
      subroutine vertical_interp(h_old,h_new)
      use grid_module,only: oldgrid,newgrid
      implicit none
!
! Interpolate fields to new vertical resolution (horizontal resolutions
! are the same). (rgrd1 is source from scd's regridpack lib).
!
! Args:
      type(history),intent(in)  :: h_old
      type(history),intent(inout) :: h_new 

      end subroutine vertical_interp
!-----------------------------------------------------------------------
      subroutine grid3d_interp(h_old,h_new)
      implicit none
#include "mksrc.h"
!
! Interpolate fields in 3d to new horizontal and vertical resolution:
! (rgrd3 is source from scd's regridpack lib).
!
! Args:
      type(history),intent(in)  :: h_old
      type(history),intent(inout) :: h_new 
!
! Local:
      integer :: ier,n,k,kbot,i,j,isign,iprint
      integer :: lw,liw,intpol(3)
      integer,allocatable :: iw(:)
      real,allocatable :: w(:)
      real :: fmin,fmax
      type(field_3d) :: fold,fnew
!
! Interpolate 2d geo fields for lower boundaries (TLBC, ULBC, VLBC, etc)
      call horizontal_interp(h_old,h_new) 
!
! Set up work arrays for the regridding routine rgrd3:
! lw = at least mx + my+2*mx + 2*mx*my+mz
!
      lw = h_new%nlon+(h_new%nlat+2*h_new%nlon)+
     |  (2*h_new%nlon*h_new%nlat+h_new%nlev)
      if (allocated(w)) deallocate(w)
      allocate(w(lw),stat=ier)
      if (ier /= 0) write(6,"('>>> grid3d_interp: error allocating',
     |  ' w(lw): lw=',i6,' ier=',i4)") lw,ier
!
! liw = at least mx+my+mz
      liw = h_new%nlon+h_new%nlat+h_new%nlev
      if (allocated(iw)) deallocate(iw)
      allocate(iw(liw),stat=ier)
      if (ier /= 0) write(6,"('>>> grid3d_interp: error allocating',
     |  ' iw(liw): liw=',i6,' ier=',i4)") liw,ier
      intpol(:) = 1 ! linear (not cubic) interp in both dimensions
!
! Do 3d interpolation.
! Fields loop:
      do n=1,h_new%nf3d_geo
        fold = h_old%f3d_geo(n)
        fnew = h_new%f3d_geo(n)
!
! Interpolate in 3-d. Regridding routine rgrd3 does not allow extrapolation,
! so this call excludes j=1 and j=nlat, and i=nlon, and bottom boundary. 
! After the rgrd3 call, these points are dealt with separately.  
!
!       write(6,"('call rgrd3: n=',i3,' new field ',a,' (',a,')',
!    |    ' h_new nlon,nlat,nlev=',3i4)") n,
!    |    h_new%f3d_geo(n)%short_name(1:8),h_new%f3d_geo(n)%levdimname,
!    |    h_new%nlon,h_new%nlat,h_new%nlev
!       write(6,"('call rgrd3: n=',i3,' old field ',a,' (',a,')',
!    |    ' h_old nlon,nlat,nlev=',3i4)") n,
!    |    h_old%f3d_geo(n)%short_name(1:8),h_old%f3d_geo(n)%levdimname,
!    |    h_old%nlon,h_old%nlat,h_old%nlev

!       call fminmax(fold%data,size(fold%data),fmin,fmax)
!       write(6,"(/,'Before rgrd3: old      ',a,' (',a,') ',
!    |    'min,max=',2e12.4)") fold%short_name,fold%levdimname,
!    |    fmin,fmax

        fnew%data = 0.
!
! Interpolate fields at midpoints (exclude bottom boundary):
!     subroutine rgrd3(nx,ny,nz,x,y,z, p,mx,my,mz,xx,yy,zz,q,intpol,
!    +                 w,lw,iw,liw,ier)
! where (nx,ny,nz,x,y,z,p) is old grid, (mx,my,mz,xx,yy,zz,q) is new grid.
!
!     ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
!                yy(1) < y(1) or y(ny) < yy(my) (or)
!                zz(1) < z(1) or z(nz) < zz(mz)
! xx(1) = h_new%glon(1)
!
        if (fnew%levdimname == 'lev') then             ! midpoints
          call rgrd3(
     |      h_old%nlon, h_old%nlat, h_old%nlev,                     ! nx,ny,nz
     |      h_old%glon, h_old%glat, h_old%lev, fold%data,           ! x,y,z,p,
     |      h_new%nlon-1, h_new%nlat-2, h_new%nlev-1,               ! mx,my,mz
     |      h_new%glon(1:h_new%nlon-1), h_new%glat(2:h_new%nlat-1), ! xx,yy,
     |      h_new%lev(2:h_new%nlev),                                ! zz
     |      fnew%data(1:h_new%nlon-1,2:h_new%nlat-1,2:h_new%nlev),  ! q
     |      intpol,w,lw,iw,liw,ier)

          if (ier /= 0) then
            write(6,"('>>> grid3d_interp: error from rgrd3: ier=',i3,
     |        ' n=',i3,' field ',a,' at midpoints')") ier,n,
     |        fnew%short_name
            if (ier == 3) then
            write(6,"('ier=3 if xx(1) < x(1) or x(nx) < xx(mx) (or)',/,
     |                '         yy(1) < y(1) or y(ny) < yy(my) (or)',/,
     |                '         zz(1) < z(1) or z(nz) < zz(mz)')")
            write(6,"('xx(1)  = h_new%glon(1) = ',f8.2)") h_new%glon(1)
            write(6,"(' x(1)  = h_old%glon(1) = ',f8.2)") h_old%glon(1)
            write(6,"('yy(1)  = h_new%glat(2) = ',f8.2)") h_new%glat(2)
            write(6,"(' y(1)  = h_old%glat(1) = ',f8.2)") h_old%glat(1)
            write(6,"('zz(1)  = h_new%lev(2)  = ',f8.2)") h_new%lev(2)
            write(6,"(' z(1)  = h_old%lev(1)  = ',f8.2)") h_old%lev(1)
            write(6,"(' x(nx) = h_old%glon(nlon)   = ',f8.2)") 
     |        h_old%glon(h_old%nlon)
            write(6,"('xx(mx) = h_new%glon(nlon-1) = ',f8.2)") 
     |        h_new%glon(h_new%nlon-1)
            write(6,"(' y(ny) = h_old%glat(nlat)   = ',f8.2)")
     |        h_old%glat(h_old%nlat)
            write(6,"('yy(my) = h_new%glat(nlat-1) = ',f8.2)")
     |        h_new%glat(h_new%nlat-1)
            write(6,"(' z(nz) = h_old%lev(nlev)    = ',f8.2)")
     |        h_old%lev(h_old%nlev)
            write(6,"('zz(mz) = h_new%lev(nlev)    = ',f8.2)")
     |        h_new%lev(h_new%nlev)
            endif
            stop 'rgrd3'
          endif

!         write(6,"('grid3d_interp for midpoints: field ',a,
!    |      ' h_new%lev(nlev)=',f8.3,' fnew%data min,max at ',
!    |      ' top lev=',2e12.4)")
!    |      fnew%short_name(1:8),h_new%lev(h_new%nlev),
!    |      minval(fnew%data(:,:,h_new%nlev)),
!    |      maxval(fnew%data(:,:,h_new%nlev))

          if (all(fold%data(:,:,h_old%nlev)==spval)) then
            write(6,"('grid3d_interp: old field ',a,' has all spval in '
     |        ,'top midpoint lev ',f8.3)") fold%short_name,
     |        h_old%lev(h_old%nlev)
            fnew%data(:,:,h_new%nlev) = spval
            fnew%data(:,:,h_new%nlev-1) = fnew%data(:,:,h_new%nlev-2)
          endif
          kbot = 2
!
! Interpolate fields at interfaces (including top and bottom boundaries):
!
        else                                           ! interfaces
          call rgrd3(
     |      h_old%nlon, h_old%nlat, h_old%nilev,
     |      h_old%glon, h_old%glat, h_old%ilev, fold%data,
     |      h_new%nlon-1, h_new%nlat-2, h_new%nilev,
     |      h_new%glon(1:h_new%nlon-1), h_new%glat(2:h_new%nlat-1),
     |      h_new%ilev,
     |      fnew%data(1:h_new%nlon-1,2:h_new%nlat-1,1:h_new%nilev), 
     |      intpol,w,lw,iw,liw,ier)
          if (ier /= 0) then
            write(6,"('>>> grid3d_interp: error from rgrd3: ier=',i3,
     |        ' n=',i3,' field ',a,' at interfaces')") ier,n,
     |        h_new%f3d_geo(n)%short_name
            stop 'rgrd3'
          endif

!         write(6,"('grid3d_interp for interfaces: field ',a,
!    |      ' h_new%ilev(nilev)=',f8.3,' fnew%data min,max at ',
!    |      ' top ilev=',2e12.4)")
!    |      fnew%short_name(1:8),h_new%ilev(h_new%nilev),
!    |      minval(fnew%data(:,:,h_new%nilev)),
!    |      maxval(fnew%data(:,:,h_new%nilev))

          if (all(fold%data(:,:,h_old%nilev)==spval)) then
            write(6,"('grid3d_interp: field ',a,' has all spval in ',
     |        'top interface ilev ',f8.3)") fnew%short_name,
     |        h_old%ilev(h_old%nilev)
            fnew%data(:,:,h_new%nilev) = spval
            fnew%data(:,:,h_new%nilev-1) = fnew%data(:,:,h_new%nilev-2)
          endif

          kbot = 1
        endif ! midpoints or interfaces

!       call fminmaxspv(fnew%data,fnew%id1*fnew%id2*fnew%id3,
!    |   fmin,fmax,fnew%missing_value)
!       write(6,"('After rgrd3: field     ',a,' (',a,') min,max=',
!    |    2e12.4)") 
!    |    fnew%short_name,fnew%levdimname,fmin,fmax
!
! Vertical loop to handle horizontal boundaries (exclude bottom midpoint):
        do k=kbot,h_new%nlev
!
! Interpolate between i=1 and i=nlon-1 to get i=nlon:
          fnew%data(h_new%nlon,:,k) = 
     |      0.5*(fnew%data(1,:,k)+fnew%data(h_new%nlon-1,:,k))
!
! Extrapolate to latitude boundaries:
          fnew%data(:,1,k) = 1.5*fnew%data(:,2,k)-0.5*fnew%data(:,3,k)

          fnew%data(:,h_new%nlat,k) = 1.5*fnew%data(:,h_new%nlat-1,k)-
     |                                0.5*fnew%data(:,h_new%nlat-2,k)
        enddo ! k=2,h_new%nlev

!       call fminmaxspv(fnew%data,fnew%id1*fnew%id2*fnew%id3,
!    |   fmin,fmax,fnew%missing_value)
!       write(6,"('After horiz boundaries:',a,' (',a,') min,max=',
!    |    2e12.4)") fnew%short_name,fnew%levdimname,fmin,fmax
!
! Extrapolate to bottom boundary midpoints:
        if (kbot==2) then
          do j=1,h_new%nlat
            fnew%data(:,j,1) = 1.5*fnew%data(:,j,2)- 
     |                         0.5*fnew%data(:,j,3)
            if (trim(h_new%f3d_geo(n)%type) == 'DENSITY') then
              do i=1,h_new%nlon
                if (fnew%data(i,j,1) <= 0.) 
     |            fnew%data(i,j,1) = fnew%data(i,j,2)
              enddo          
            endif
          enddo ! j=1,h_new%nlat
        endif
 
!       call fminmaxspv(fnew%data,fnew%id1*fnew%id2*fnew%id3,
!    |   fmin,fmax,fnew%missing_value)
!       write(6,"('Final 3d field:        ',a,' (',a,') min,max=',
!    |    2e12.4)") 
!    |    fnew%short_name,fnew%levdimname,fmin,fmax

        h_new%f3d_geo(n) = fnew

      enddo ! n=1,h_new%nf3d_geo

      end subroutine grid3d_interp
!-----------------------------------------------------------------------
      subroutine allocdata(h)
      implicit none
!
! Args:
      type(history),intent(inout) :: h    ! h_new
!
! Local:
      integer :: i,id1,id2,id3,ier,n
      real :: fmin,fmax
!
! 3d geo data:
      if (h%nf3d_geo > 0) then
        h%f3d_geo%id1 = newgrid%nlon
        h%f3d_geo%id2 = newgrid%nlat
        h%f3d_geo%id3 = newgrid%nlev
        do i=1,h%nf3d_geo
          id1 = h%f3d_geo(i)%id1
          id2 = h%f3d_geo(i)%id2
          id3 = h%f3d_geo(i)%id3

!         write(6,"('Allocating 3d data for field ',
!    |      a,' id1,2,3=',3i4,' (',a,')')") 
!    |      h%f3d_geo(i)%short_name(1:8),id1,id2,id3,
!    |      h%f3d_geo(i)%levdimname(1:4)

          if (h%f3d_geo(i)%levdimname == 'ilev')
     |      h%f3d_geo(i)%id3 = newgrid%nilev

          if (associated(h%f3d_geo(i)%data)) 
     |      deallocate(h%f3d_geo(i)%data)

          allocate(h%f3d_geo(i)%data(id1,id2,id3),stat=ier)

          if (ier /= 0) then
            write(6,"('>>> Error allocating h%f3d_geo(i)%data:',
     |        ' id1,2,3=',3i4)") id1,id2,id3
            stop 'allocation error'
          else
!           write(6,"('i=',i4,' allocated h%f3d_geo(i)%data for field ',
!    |        a,' id1,2,3=',3i4,' (',a,')')") i,
!    |        h%f3d_geo(i)%short_name(1:8),id1,id2,id3,
!    |        h%f3d_geo(i)%levdimname(1:4)
          endif
        enddo ! i=1,nf3d_geo
      endif ! nf3d_geo > 0
!
! 3d mag data:
!     if (h%nf3d_mag > 0) then
!       h%f3d_mag%id1 = newgrid%nlon
!       h%f3d_mag%id2 = newgrid%nlat
!       h%f3d_mag%id3 = newgrid%nlev
!       do i=1,h%nf3d_mag
!         id1 = h%f3d_mag(i)%id1
!         id2 = h%f3d_mag(i)%id2
!         id3 = h%f3d_mag(i)%id3
!         if (h%f3d_mag(i)%levdimname == 'ilev')
!    |      h%f3d_mag%id3 = newgrid%nilev
!         if (associated(h%f3d_mag(i)%data)) 
!    |      deallocate(h%f3d_mag(i)%data)
!         allocate(h%f3d_mag(i)%data(id1,id2,id3),stat=ier)
!         if (ier /= 0) then
!           write(6,"('>>> Error allocating h%f3d_mag(i)%data:',
!    |        ' id1,2,3=',3i4)") id1,id2,id3
!           stop 'allocation error'
!         else
!           write(6,"('allocated h%f3d_mag(i)%data for field ',
!    |        a,' id1,2,3=',3i4)") h%f3d_mag(i)%short_name(1:8),
!    |        id1,id2,id3
!         endif
!       enddo ! i=1,nf3d_mag
!     endif ! nf3d_mag > 0
!
! 2d geo data:
      if (h%nf2d_geo > 0) then
        h%f2d_geo%id1 = newgrid%nlon
        h%f2d_geo%id2 = newgrid%nlat
        do i=1,h%nf2d_geo
          id1 = h%f2d_geo(i)%id1
          id2 = h%f2d_geo(i)%id2
          if (associated(h%f2d_geo(i)%data)) 
     |      deallocate(h%f2d_geo(i)%data)
          allocate(h%f2d_geo(i)%data(id1,id2),stat=ier)
          if (ier /= 0) then
            write(6,"('>>> Error allocating h%f2d_geo(i)%data:',
     |        ' id1,2=',2i4)") id1,id2
            stop 'allocation error'
          else
!           write(6,"('allocated h%f2d_geo(i)%data for field ',
!    |        a,' id1,2=',2i4)") h%f2d_geo(i)%short_name(1:8),
!    |        id1,id2
          endif
        enddo ! i=1,nf2d_geo
      endif ! nf2d_geo > 0
!
! 2d mag data:
!     if (h%nf2d_mag > 0) then
!       h%f2d_mag%id1 = newgrid%nlon
!       h%f2d_mag%id2 = newgrid%nlat
!       do i=1,h%nf2d_mag
!         id1 = h%f2d_mag(i)%id1
!         id2 = h%f2d_mag(i)%id2
!         if (associated(h%f2d_mag(i)%data)) 
!    |      deallocate(h%f2d_mag(i)%data)
!         allocate(h%f2d_mag(i)%data(id1,id2),stat=ier)
!         if (ier /= 0) then
!           write(6,"('>>> Error allocating h%f2d_mag(i)%data:',
!    |        ' id1,2=',2i4)") id1,id2
!           stop 'allocation error'
!         else
!           write(6,"('allocated h%f2d_mag(i)%data for field ',
!    |        a,' id1,2=',2i4)") h%f2d_mag(i)%short_name(1:8),
!    |        id1,id2
!         endif
!       enddo ! i=1,nf2d_mag
!     endif ! nf2d_mag > 0
      end subroutine allocdata
!-----------------------------------------------------------------------
      subroutine copyhist(hsrc,hdest)
!
! Copy hsrc into hdest.
!
      implicit none
!
! Args:
      type(history),intent(in) :: hsrc 
      type(history),intent(out)  :: hdest
!
! Local:
      integer :: i
!
! Copy entire structure, including pointers:
      hdest = hsrc
!
! Allocate new pointers the size of the source pointers: 
!
!     allocate(hdest%fnames3d_geo(size(hsrc%fnames3d_geo)))
!     allocate(hdest%fnames3d_mag(size(hsrc%fnames3d_mag)))
!     allocate(hdest%fnames2d_geo(size(hsrc%fnames2d_geo)))
!     allocate(hdest%fnames2d_mag(size(hsrc%fnames2d_mag)))
      allocate(hdest%glat(size(hsrc%glat)))
      allocate(hdest%glon(size(hsrc%glon)))
      allocate(hdest%lev(size(hsrc%lev)))
      allocate(hdest%ilev(size(hsrc%ilev)))
      allocate(hdest%gmlat(size(hsrc%gmlat)))
      allocate(hdest%gmlon(size(hsrc%gmlon)))
      allocate(hdest%mtimes(size(hsrc%mtimes,1),
     |                      size(hsrc%mtimes,2)))
!
! Assign from source structure:
!
!     hdest%fnames3d_geo = hsrc%fnames3d_geo
!     hdest%fnames3d_mag = hsrc%fnames3d_mag
!     hdest%fnames2d_geo = hsrc%fnames2d_geo
!     hdest%fnames2d_mag = hsrc%fnames2d_mag
      hdest%glat   = hsrc%glat
      hdest%glon   = hsrc%glon
      hdest%lev    = hsrc%lev
      hdest%ilev   = hsrc%ilev
      hdest%gmlat  = hsrc%gmlat
      hdest%gmlon  = hsrc%gmlon
      hdest%mtimes = hsrc%mtimes
!
! 3d and 2d fields:
!
      do i=1,hsrc%nf3d_geo
        allocate(hdest%f3d_geo(i)%data(size(hsrc%f3d_geo(i)%data,1),
     |                                 size(hsrc%f3d_geo(i)%data,2),
     |                                 size(hsrc%f3d_geo(i)%data,3)))
        hdest%f3d_geo(i)%data = hsrc%f3d_geo(i)%data
      enddo
      do i=1,hsrc%nf3d_mag
        allocate(hdest%f3d_mag(i)%data(size(hsrc%f3d_mag(i)%data,1),
     |                                 size(hsrc%f3d_mag(i)%data,2),
     |                                 size(hsrc%f3d_mag(i)%data,3)))
        hdest%f3d_mag(i)%data = hsrc%f3d_mag(i)%data
      enddo
      do i=1,hsrc%nf2d_geo
        allocate(hdest%f2d_geo(i)%data(size(hsrc%f2d_geo(i)%data,1),
     |                                 size(hsrc%f2d_geo(i)%data,2)))
        hdest%f2d_geo(i)%data = hsrc%f2d_geo(i)%data
      enddo
      do i=1,hsrc%nf2d_mag
        allocate(hdest%f2d_mag(i)%data(size(hsrc%f2d_mag(i)%data,1),
     |                                 size(hsrc%f2d_mag(i)%data,2)))
        hdest%f2d_mag(i)%data = hsrc%f2d_mag(i)%data
      enddo
      end subroutine copyhist
!-----------------------------------------------------------------------
      end module newhist_module