!-----------------------------------------------------------------------
  subroutine handle_ncerr(istat,msg,ifatal)
    use netcdf
    use params,only: iulog
    implicit none
!
! Handle a netcdf lib error:
!
    integer,intent(in) :: istat,ifatal
    character(len=*),intent(in) :: msg
!
    write(iulog,"(/72('-'))")
    write(iulog,"('>>> Error from netcdf library:')")
    write(iulog,"(a)") trim(msg)
    write(iulog,"('istat=',i5)") istat
    write(iulog,"(a)") nf90_strerror(istat)
    write(iulog,"(72('-')/)")
    if (ifatal > 0) call shutdown('Fatal netcdf error')
  end subroutine handle_ncerr
!-----------------------------------------------------------------------
  subroutine shutdown(msg)
    use esmf
    use params,only: iulog
    implicit none
#include <mpif.h>
!
! Args:
    character(len=*),intent(in) :: msg
!
! Local:
    integer :: ier

    write(iulog,"(/,28('>'),' SHUTDOWN ',28('<'))")
    write(iulog,"('Shutdown stop message: ',a)") trim(msg)

    call ESMF_Finalize(endflag=ESMF_END_KEEPMPI,rc=ier)
    if (ier /= ESMF_SUCCESS) &
      write(iulog,"('>>> Shutdown: error return from ESMF_Finalize: ier=',i4)") ier

    write(iulog,"(72('>'))")
!
! Shutdown MPI:
    call MPI_Abort(MPI_COMM_WORLD,ier)
    if (ier /= 0) &
      write(iulog,"('>>> Shutdown: error return from MPI_Abort: ier=',i4)") ier

    stop 'shutdown'
  end subroutine shutdown
!-----------------------------------------------------------------------
  integer function julian_day(month,day)
!
! Given 2-digit month and day, return julian day (assuming no leap-years)
!
! Args:
  integer,intent(in) :: month,day
!
! Local:
  integer :: i
  integer :: days_per_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

    julian_day = 0
    do i=1,month-1
      julian_day = julian_day + days_per_month(i)
    enddo
    julian_day = julian_day+day
!   write(iulog,"('julian_day: month=',i3,' day=',i3,' julian_day=',i4)") &
!     month,day,julian_day
  end function julian_day
!-----------------------------------------------------------------------
  integer function ixfind(iarray,idim,itarget,icount)
!
! Search iarray(idim) for itarget, returning first index in iarray
! where iarray(idim)==target. Also return number of elements of
! iarray that == itarget in icount.
!
! Args:
    integer,intent(in) :: idim,itarget
    integer,intent(in) :: iarray(idim)
    integer,intent(out) :: icount
!
! Local:
    integer :: i
!
    ixfind = 0
    icount = 0
    if (.not.any(iarray==itarget)) return
    icount = count(iarray==itarget)
    do i=1,idim
      if (iarray(i)==itarget) then
        ixfind = i
        exit
      endif
    enddo
  end function ixfind
!-----------------------------------------------------------------------
  subroutine reverse_vec(vec,n)
    use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals
!
! Reverse order of elements in vector vec(n)
!
    integer,intent(in) :: n
    real(r8),intent(inout) :: vec(n)

    real(r8) :: tmp(n)
    integer :: i

    do i=1,n
      tmp(i) = vec(n-i+1)
    enddo
    vec = tmp

  end subroutine reverse_vec
!-----------------------------------------------------------------------
  subroutine shift_lon(f,nlon,lonseq,iscoord)
    use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals
    use params,only: iulog
!
! Shift longitude vector f(nlon) forward 180 degrees according to input 
! string lonseq. Input f can be either arbitrary field values or 
! the coordinate array itself. Shift f in the 'lonseq' manner, as follows:
!
! If lonseq='-180to180', then shift from 0->360 to -180->+180
! If lonseq='zeroto360', then shift from -180->+180 to 0->360
!
! WARNING: This routine works with WACCM-X history files, where nlon=144, 72, or 80
!          It has not been tested with other models or resolutions.
!          (e.g., there is no test for center point, its assumed to be nlon/2)
!
! Args:
    integer,intent(in) :: nlon
    real(r8),intent(inout) :: f(nlon)
    character(len=*),intent(in) :: lonseq 
    logical,intent(in) :: iscoord ! if true, f is a coordinate, otherwise it is data
!
! Local:
    character(len=80) :: msg
    integer :: ihalf

    if (lonseq /= '-180to180'.and.lonseq /= 'zeroto360') then
      write(msg,"('shift_lon: bad lonseq=',a,' must be either ''-180to180'' or ''zeroto360''')") &
        lonseq
      call shutdown(msg)
    endif

    ihalf = nlon/2
    if (lonseq == '-180to180') then ! shift to -180 -> +180
      f = cshift(f,ihalf)           ! cshift is circular shift intrinsic
      if (iscoord) then
        do i=1,ihalf
          f(i) = f(i)-360._r8
        enddo
      endif
    else                           ! shift to 0 -> 360
      f = cshift(f,ihalf)          ! cshift is circular shift intrinsic
      if (iscoord) then
        do i=ihalf+1,nlon
          f(i) = f(i)+360._r8
        enddo
      endif
    endif 

  end subroutine shift_lon
!-----------------------------------------------------------------------