!//=========================================================================
!// Oslo CTM3
!//=========================================================================
!// Ole Amund Sovde, May 2015
!//=========================================================================
!// netCDF utilities.
!//=========================================================================
module ncutils
  !// ----------------------------------------------------------------------
  !// MODULE: ncutils
  !// DESCRIPTION: Routines for extracting different fields from
  !//              netcdf-files.
  !//
  !// Contains
  !//   subroutine readnc_3d_from4d
  !//   subroutine readnc_2d_from3d
  !//   subroutine readnc_1d
  !//   subroutine readnc_1d_from2d
  !//   subroutine get_netcdf_var_1d
  !//   subroutine get_netcdf_var_2d
  !//   subroutine get_netcdf_var_3d
  !//   subroutine get_netcdf_var_4d
  !//   subroutine get_netcdf_att_char
  !//   subroutine get_netcdf_var_dims
  !//   subroutine get_netcdf_dim
  !//   subroutine get_netcdf_r4var_1d_from_2d
  !//   subroutine get_netcdf_r4var_2d_from_3d
  !//   subroutine get_netcdf_r4var_3d_from_4d
  !//   subroutine handle_err
  !//   subroutine handle_error
  !//
  !// Ole Amund Sovde, March 2011
  !// ----------------------------------------------------------------------
  implicit none
  !// ----------------------------------------------------------------------
  character(len=*), parameter, private :: f90file = 'ncutils.f90'
  !//-----------------------------------------------------------------------
  public
  !// ----------------------------------------------------------------------

contains

  !// ----------------------------------------------------------------------
  subroutine readnc_3d_from4d(INFILE, DIMNAME1, DIM1, DIMNAME2, DIM2, &
       DIMNAME3, DIM3, DIMNAME4, GET4ENTRY, GETFIELDNAME, R8XYZ)
    !// --------------------------------------------------------------------
    !// Read netCDF file containing 4D fields (DIM1,DIM2,DIM3,DIM4), and
    !// returns 3D-field for DIM4-entry TSTEP. Typically, we have
    !//   DIM1: longitude
    !//   DIM2: latitude
    !//   DIM3: height
    !//   DIM4: time
    !//
    !// Ole Amund Sovde, March 2011
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in)   :: infile     !// Name of netCDFfile
    character(len=*), intent(in)   :: dimname1   !// Name of dimension1
    integer,intent(in)          :: DIM1       !// Dimension of dimname1
    character(len=*), intent(in)   :: dimname2   !// Name of dimension2
    integer,intent(in)          :: DIM2       !// Dimension of dimname2
    character(len=*), intent(in)   :: dimname3   !// Name of dimension3
    integer,intent(in)          :: DIM3       !// Dimension of dimname3
    character(len=*), intent(in)   :: dimname4   !// Name of dimension4
                                              !// No need for DIM4
    integer,intent(in)          :: GET4ENTRY  !// Entry to get from dimension4
    character(len=*), intent(in)   :: getfieldname  !// Name of field to get

    !// Output
    real(r8), intent(out)        :: R8XYZ(DIM1,DIM2,DIM3)

    !// Locals
    integer  :: field_id      !// Id for 1d-field
    integer  :: dim1_id       !// Id for dimension1
    integer  :: dim2_id       !// Id for dimension2
    integer  :: dim3_id       !// Id for dimension3
    integer  :: dim4_id       !// Id for dimension4
    integer  :: nsize1        !// Size of dim1
    integer  :: nsize2        !// Size of dim2
    integer  :: nsize3        !// Size of dim3
    integer  :: nsize4        !// Size of dim4

    integer  :: srt_entry(4)  !// Start point
    integer  :: cnt_entry(4)  !// Count indexes
    integer  :: status        !// status of process (0=OK)
    integer  :: ncid          !// file id
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='readnc_3d_from4d'
    !// --------------------------------------------------------------------

    status=nf90_noerr  !// Status is 0 and should be kept that way !!

    !// Array which tells you where to start picking your 3D field
    srt_entry = (/ 1, 1, 1, GET4ENTRY /)    !Start array
    !// Array which tells you how far to count when picking it
    cnt_entry = (/ DIM1, DIM2, DIM3, 1 /)         !Count array

    !// Open the existing file
    status = nf90_open(infile, nf90_nowrite, ncid)
    if (status /= NF90_NOERR ) call handle_error(status, &
         subr//':'//trim(infile)//':'//TRIM(getfieldname)//':open')

    !// Inquire dimension ids
    !// DIM1
    status = nf90_inq_dimid(ncid,dimname1,dim1_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 1')
    !// DIM2
    status = nf90_inq_dimid(ncid,dimname2,dim2_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 2')
    !// DIM3
    status = nf90_inq_dimid(ncid,dimname3,dim3_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 3')
    !// DIM4
    status = nf90_inq_dimid(ncid,dimname4,dim4_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 4')

    !// Inquire dimensions
    status = nf90_Inquire_Dimension(ncid,dim1_id,len=nsize1)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 1')
    if (nsize1/=DIM1) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM1 = ',nsize1
       write(6,'(a,i7)') 'your array has dimension ',DIM1
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim2_id,len=nsize2)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 2')
    if (nsize2/=DIM2) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM2 = ',nsize2
       write(6,'(a,i7)') 'your array has dimension',DIM2
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim3_id,len=nsize3)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 3')
    if (nsize3/=DIM3) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM3 = ',nsize3
       write(6,'(a,i7)') '    your array has dimension',DIM3
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim4_id,len=nsize4)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 4')
    !// Check GET4ENTRY against nsize4
    if (GET4ENTRY .gt. nsize4 .or. nsize4 .le.0) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM4 = ',nsize4
       write(6,'(a,i7)') 'you try to get entry ',GET4ENTRY
       stop
    endif

    !// Number of variables
    !status = nf90_Inquire(ncid,nDims,nVars)

    !// Get variable ID from name
    status=nf90_inq_varid(ncid,getfieldname,field_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_varid')

    !// Get the variable you want and put it in the threeDfield array
    status=nf90_get_var(ncid,field_id,R8XYZ, &
         start=srt_entry, count=cnt_entry )
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':get_var')

    write(6,'(a)') f90file//':'//subr//': Got variable '//trim(getfieldname)

    !// Close file
    status=nf90_close(ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':close')

    !// --------------------------------------------------------------------
  end subroutine readnc_3d_from4d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine readnc_2d_from3d(INFILE, DIMNAME1, DIM1, DIMNAME2, DIM2, &
       DIMNAME3, GET3ENTRY, GETFIELDNAME, R8XY)
    !// --------------------------------------------------------------------
    !// Read netCDF file containing 4D fields (DIM1,DIM2,DIM3), and
    !// returns 2D-field for DIM3-entry TSTEP. Typically, we have
    !//   DIM1: longitude
    !//   DIM2: latitude
    !//   DIM3: time
    !//
    !// Amund Sovde, November 2014
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in) :: infile     !// Name of netCDFfile
    character(len=*), intent(in) :: dimname1   !// Name of dimension1
    integer, intent(in)          :: DIM1       !// Dimension of dimname1
    character(len=*), intent(in) :: dimname2   !// Name of dimension2
    integer, intent(in)          :: DIM2       !// Dimension of dimname2
    character(len=*), intent(in) :: dimname3   !// Name of dimension3

    integer, intent(in)          :: GET3ENTRY  !// Entry to get from dimension3
    character(len=*), intent(in) :: getfieldname  !// Name of field to get

    !// Output
    real(r8), intent(out)        :: R8XY(DIM1,DIM2)

    !// Locals
    integer  :: field_id      !// Id for 1d-field
    integer  :: dim1_id       !// Id for dimension1
    integer  :: dim2_id       !// Id for dimension2
    integer  :: dim3_id       !// Id for dimension3
    integer  :: nsize1        !// Size of dim1
    integer  :: nsize2        !// Size of dim2
    integer  :: nsize3        !// Size of dim3

    integer  :: srt_entry(3)  !// Start point
    integer  :: cnt_entry(3)  !// Count indexes
    integer  :: status        !// status of process (0=OK)
    integer  :: ncid          !// file id
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='readnc_2d_from3d'
    !// --------------------------------------------------------------------

    status=nf90_noerr  !// Status is 0 and should be kept that way !!

    !// Array which tells you where to start picking your 3D field
    srt_entry = (/ 1, 1, GET3ENTRY /)    !Start array
    !// Array which tells you how far to count when picking it
    cnt_entry = (/ DIM1, DIM2, 1 /)      !Count array

    !// Open the existing file
    status=nf90_open(infile, nf90_nowrite, ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':open')

    !// Inquire dimension ids
    !// DIM1
    status = nf90_inq_dimid(ncid,dimname1,dim1_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 1')
    !// DIM2
    status = nf90_inq_dimid(ncid,dimname2,dim2_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 2')
    !// DIM3
    status = nf90_inq_dimid(ncid,dimname3,dim3_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 3')

    !// Inquire dimensions
    status = nf90_Inquire_Dimension(ncid,dim1_id,len=nsize1)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 1')
    if (nsize1/=DIM1) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM1 = ',nsize1
       write(6,'(a,i7)') 'your array has dimension ',DIM1
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim2_id,len=nsize2)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 2')
    if (nsize2/=DIM2) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM2 = ',nsize2
       write(6,'(a,i7)') 'your array has dimension',DIM2
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim3_id,len=nsize3)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 3')
    !// Check GET3ENTRY against nsize3
    if (GET3ENTRY .gt. nsize3 .or. nsize3 .le.0) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM3 = ',nsize3
       write(6,'(a,i7)') 'you try to get entry ',GET3ENTRY
       stop
    endif

    !// Number of variables
    !status = nf90_Inquire(ncid,nDims,nVars)

    !// Get variable ID from name
    status=nf90_inq_varid(ncid,getfieldname,field_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_varid')

    !// Get the variable you want and put it in the threeDfield array
    status=nf90_get_var(ncid,field_id,R8XY, &
         start=srt_entry, count=cnt_entry )
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':get_var')

    write(6,'(a)') f90file//':'//subr//': Got variable '//trim(getfieldname)

    !// Close file
    status=nf90_close(ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':close')

    !// --------------------------------------------------------------------
  end subroutine readnc_2d_from3d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine readnc_sub3d_from3d(INFILE, DIMNAME1, DIM1, DIMNAME2, DIM2, &
       DIMNAME3, DIM3, GETFIELDNAME, srt_entry, cnt_entry, R8XY)
    !// --------------------------------------------------------------------
    !// Read netCDF file containing 4D fields (DIM1,DIM2,DIM3), and
    !// returns 2D-field for DIM3-entry TSTEP. Typically, we have
    !//   DIM1: longitude
    !//   DIM2: latitude
    !//   DIM3: time
    !//
    !// Amund Sovde Haslerud, January 2017
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in) :: infile     !// Name of netCDFfile
    character(len=*), intent(in) :: dimname1   !// Name of dimension1
    integer, intent(in)          :: DIM1       !// Dimension of dimname1
    character(len=*), intent(in) :: dimname2   !// Name of dimension2
    integer, intent(in)          :: DIM2       !// Dimension of dimname2
    character(len=*), intent(in) :: dimname3   !// Name of dimension3
    integer, intent(in)          :: DIM3       !// Dimension of dimname3

    character(len=*), intent(in) :: getfieldname  !// Name of field to get
    integer, intent(in)  :: srt_entry(3)  !// Start point
    integer, intent(in)  :: cnt_entry(3)  !// Count indexes

    !// Output (note the third dimension)
    real(r8), intent(out)        :: R8XY(DIM1,DIM2,CNT_ENTRY(3))

    !// Locals
    integer  :: field_id      !// Id for 1d-field
    integer  :: dim1_id       !// Id for dimension1
    integer  :: dim2_id       !// Id for dimension2
    integer  :: dim3_id       !// Id for dimension3
    integer  :: nsize1        !// Size of dim1
    integer  :: nsize2        !// Size of dim2
    integer  :: nsize3        !// Size of dim3

    integer  :: status        !// status of process (0=OK)
    integer  :: ncid          !// file id
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='readnc_sub3d_from3d'
    !// --------------------------------------------------------------------

    status=nf90_noerr  !// Status is 0 and should be kept that way !!

    !// Open the existing file
    status=nf90_open(infile, nf90_nowrite, ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':open')

    !// Inquire dimension ids
    !// DIM1
    status = nf90_inq_dimid(ncid,dimname1,dim1_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 1')
    !// DIM2
    status = nf90_inq_dimid(ncid,dimname2,dim2_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 2')
    !// DIM3
    status = nf90_inq_dimid(ncid,dimname3,dim3_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_dimid 3')

    !// Inquire dimensions
    status = nf90_Inquire_Dimension(ncid,dim1_id,len=nsize1)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 1')
    if (nsize1/=DIM1) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM1 = ',nsize1
       write(6,'(a,i7)') 'your array has dimension ',DIM1
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim2_id,len=nsize2)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 2')
    if (nsize2/=DIM2) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM2 = ',nsize2
       write(6,'(a,i7)') 'your array has dimension',DIM2
       stop
    endif

    status = nf90_Inquire_Dimension(ncid,dim3_id,len=nsize3)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(infile)//':'//trim(getfieldname)//':inquire_dimension 3')
    !// Check DIM3 against size of srt_entry and cnt_entry
    !// Remember that cnt_entry(3) is 12 when reading 12 datasets,
    !// so that when checking against DIM3, we have to use
    !// srt_entry(3)+cnt_entry(3)-1:
    if (srt_entry(3) .gt. DIM3 .or. &
         (srt_entry(3)+cnt_entry(3)-1) .gt. DIM3) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(infile)
       write(6,'(a,i7)') '    reports DIM3 = ',DIM3
       write(6,'(a,i7,a3,i7)') 'you try to get entries ',&
            srt_entry(3), ' : ', (srt_entry(3)+cnt_entry(3)-1)
       stop
    endif

    !// Number of variables
    !status = nf90_Inquire(ncid,nDims,nVars)

    !// Get variable ID from name
    status=nf90_inq_varid(ncid,getfieldname,field_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':inq_varid')

    !// Get the variable you want and put it in the threeDfield array
    status=nf90_get_var(ncid,field_id,R8XY, &
         start=srt_entry, count=cnt_entry )
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':get_var')

    write(6,'(a)') f90file//':'//subr//': Got variable '//trim(getfieldname)

    !// Close file
    status=nf90_close(ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(infile)//':'//trim(getfieldname)//':close')

    !// --------------------------------------------------------------------
  end subroutine readnc_sub3d_from3d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine readnc_1d(oneDfield, getfieldname, dimname1, DIM1, filename)
    !// --------------------------------------------------------------------
    !// Open a netCDF file and get 1D field from the file.
    !// Reads _only_ 1D fields, i.e., not by extracting 1D from 2D fields
    !// (see readnc_1d_from2d for that).
    !//
    !// Ole Amund Sovde, March 2011
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in)       :: getfieldname !// Name of field to get
    character(len=*), intent(in)       :: filename     !// Name of netCDFfile
    character(len=*), intent(in)       :: dimname1     !// Name of dimension1
    integer,intent(in)              :: DIM1         !// Dimension of dimname1
                                                    !//    and getfieldname

    !// Output
    real(r8), intent(out)   :: oneDfield(DIM1)        !// 1D field

    !// Locals
    integer  :: dim_id       !// Id for dimension
    integer  :: field_id     !// Id for 1d-field
    integer  :: nsize        !// Size of 1d-field, read from file

    integer  :: srt_entry(1) !// Start point
    integer  :: cnt_entry(1) !// Count indexes
    integer  :: status       !// status of process (0=OK)
    integer  :: ncid         !// file id
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='readnc_1d'
    !// --------------------------------------------------------------------

    !// Array which tells you where to start picking your field
    srt_entry = (/ 1 /)    !Start array
    !// Array which tells you how far to count when picking it
    cnt_entry = (/ DIM1 /)   !Count array

    status=nf90_noerr  !Status is 0 and should be kept that way !!

    !// Open the existing file
    status=nf90_open(filename, nf90_nowrite, ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':open')

    !// Inquire dimension ID
    status = nf90_inq_dimid(ncid,dimname1,dim_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inq_dimid')

    !// Inquire dimension (i.e. size)
    status = nf90_Inquire_Dimension(ncid,dim_id,len=nsize)
    if (status /= NF90_NOERR) call handle_error(status, &
       subr//':'//trim(filename)//':'//trim(getfieldname)//':inquire_dimension')
    if (nsize /= DIM1) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(filename)
       write(6,'(a,i7)') '    reports DIM1 = ',nsize
       write(6,'(a,i7)') 'your array has dimension ',DIM1
       stop
    endif

    !// Get variable ID for 1D field
    status=nf90_inq_varid(ncid,getfieldname,field_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inq_varid')

    !// Get the variable you want and put it in the oneDfield array
    status=nf90_get_var(ncid,field_id,oneDfield, &
         start=srt_entry, count=cnt_entry )
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':get_var')

    write(6,'(a)') f90file//':'//subr//': Got variable ',trim(getfieldname)

    !Closing file
    status=nf90_close(ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':close')

    !// --------------------------------------------------------------------
  end subroutine readnc_1d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine readnc_1d_from2d(oneDfield, getfieldname, dimname1, DIM1, &
       dimname2, GETENTRY, filename )
    !// --------------------------------------------------------------------
    !// Open a netCDF file and get 1D field from a 2D field on the file.
    !// Reads a 1D field for specified entry of dimension 2 (DIM2).
    !// DIM2 does not have to be sent in to routine.
    !// In other words, we extract the whole DIM1 for a specified entry
    !// of DIM2. This entry is specified by GETENTRY.
    !//
    !// Ole Amund Sovde, March 2011
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in)  :: filename      !// Name of netCDFfile
    character(len=*), intent(in)  :: getfieldname  !// Name of field to get
    character(len=*), intent(in)  :: dimname1      !// Name of dimension1
    integer,intent(in)         :: DIM1          !// Dimension of dimname1
    character(len=*), intent(in)  :: dimname2      !// Name of dimension2
                                                !// No need for DIM2
    integer,intent(in)         :: GETENTRY      !// Entry to get from dimension2

    !// Output
    real(r8), intent(out)  :: oneDfield(DIM1)  !// 1D field

    !// Locals
    integer  :: dim1_id       !// Id for dimension1
    integer  :: dim2_id       !// Id for dimension2
    integer  :: field_id      !// Id for 1d-field
    integer  :: nsize1        !// Size of dim1
    integer  :: nsize2        !// Size of dim2

    integer  :: srt_entry(2)  !// Start point
    integer  :: cnt_entry(2)  !// Count indexes
    integer  :: status        !// status of process (0=OK)
    integer  :: ncid          !// file id
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='readnc_1d_from2d'
    !// --------------------------------------------------------------------

    !// Array which tells you where to start picking your field
    srt_entry = (/ 1, GETENTRY /)    !Start array
    !// Array which tells you how far to count when picking it
    cnt_entry = (/ DIM1, 1 /)   !Count array

    status=nf90_noerr  !Status is 0 and should be kept that way !!

    !// Open the existing file
    status=nf90_open(filename, nf90_nowrite, ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':open')

    !// Inquire dimension1 ID
    status = nf90_inq_dimid(ncid,dimname1,dim1_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inq_dimid 1')

    !// Inquire dimension1 (i.e. size)
    status = nf90_Inquire_Dimension(ncid,dim1_id,len=nsize1)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inquire_dimension 1')
    if (nsize1 /= DIM1) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(filename)
       write(6,'(a,i7)') '    reports DIM1 = ',nsize1
       write(6,'(a,i7)') 'your array has dimension ',DIM1
       stop
    endif

    !// Inquire dimension2 ID
    status = nf90_inq_dimid(ncid,dimname2,dim2_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inq_dimid 2')

    !// Inquire dimension1 (i.e. size)
    status = nf90_Inquire_Dimension(ncid,dim2_id,len=nsize2)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inquire_dimension 2')
    !// Test against GETENTRY
    if (GETENTRY .gt. nsize2 .or. nsize2 .le.0) then
       write(6,'(a)') f90file//':'//subr//': File '//trim(filename)
       write(6,'(a,i7)')'    reports DIM2 = ',nsize2
       write(6,'(a,i7)')'you try to get entry ',GETENTRY
       stop
    endif

    !// Get variable ID for 2D field
    status=nf90_inq_varid(ncid,getfieldname,field_id)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':inq_varid')

    !// Get the variable you want and put it in the oneDfield array
    status=nf90_get_var(ncid,field_id,oneDfield, &
         start=srt_entry, count=cnt_entry )
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':get_var')

    write(6,'(a)') f90file//':'//subr//': Got variable ',trim(getfieldname)

    !Closing file
    status=nf90_close(ncid)
    if (status /= NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//trim(getfieldname)//':close')

    !// --------------------------------------------------------------------
  end subroutine readnc_1d_from2d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_0d( filename, varname, data )
    !// --------------------------------------------------------------------
    !// Reads 0D variable (varname) from netCDF file.
    !//
    !// Amund Sovde Haslerud, January 2017
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    !// Output
    real(r8), intent(out):: data

    !// Local variables
    integer :: status, ncID, varID, dimLen, xtype, ndims
    integer :: dimIDs(1)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_0d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':inq_varid')

    !// Variable information
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs,xtype=xtype,ndims=ndims )

!    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
!         subr//':'//trim(filename)//':'//TRIM(varname)//':variable')

!    !// Dimension of the variable
!    status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(1), len=dimLen )
!    !IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
!    IF (status /= NF90_EBADDIM) CALL HANDLE_ERROR(status, &
!         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension')

!    if (dimLen .ne. 0) then
    if (ndims .ne. 0) then
       write(6,'(a,i3)') f90file//':'//subr//': Not zero dimension length: ',ndims
       stop
    end if

    !// Read the variable
    status = NF90_GET_VAR( ncid, varID, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':get_var')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_0d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_1d( filename, varname, data )
    !// --------------------------------------------------------------------
    !// Reads 1D variable (varname) from netCDF file.
    !// From Chris D. Holmes, UCI
    !//
    !// Ole Amund Sovde, March 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    !// Output
    real(r8),allocatable,intent(out):: data(:)

    !// Local variables
    integer :: status, ncID, varID, dimLen
    integer :: dimIDs(1)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_1d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':inq_varid')

    !// Variable information
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':variable')

    !// Dimension of the variable
    status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(1), len=dimLen )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 1')

    !// Allocate the variable
    ALLOCATE( data(dimLen), STAT=status )
    IF (status < 0 ) then
       write(6,'(a)') f90file//':'//subr//': data allocation failed'
       stop
    end if

    !// Read the variable
    status = NF90_GET_VAR( ncid, varID, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':get_var')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_1d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_2d( filename, varname, data, nlon,nlat )
    !// --------------------------------------------------------------------
    !// Reads 2D variable (varname) from netCDF file.
    !// From Chris D. Holmes, UCI
    !//
    !// Ole Amund Sovde, September 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    integer, intent(in)           :: nlon,nlat
    !// Output
    real(r8),intent(out)            :: data(nlon,nlat)

    !// Local variables
    integer :: status, ncID, varID
    integer :: dimIDs(3), dimLen(3)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_2d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':inq_varid')

    !// Variable information
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':variable')

    !// 1st Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(1),len=dimLen(1))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 1')

    !// 2nd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(2),len=dimLen(2))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 2')

    !// Read the variable
    status = NF90_GET_VAR( ncid, varID, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':get_var')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_2d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_3d( filename, varname, data, nlon,nlat,ntime )
    !// --------------------------------------------------------------------
    !// Reads 3D variable (varname) from netCDF file.
    !// From Chris D. Holmes, UCI
    !//
    !// Ole Amund Sovde, March 2012
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in)  :: filename, varname
    integer, intent(in)           :: nlon,nlat,ntime
    !// Output
    real(r8), intent(out)         :: data(nlon,nlat,ntime)

    !// Local variables
    integer :: status, ncID, varID
    integer :: dimIDs(3), dimLen(3)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_3d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':inq_varid')

    !// Variable information
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':variable')

    !// 1st Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(1),len=dimLen(1))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 1')

    !// 2nd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(2),len=dimLen(2))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 2')

    !// 3rd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(3),len=dimLen(3))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 3')

    !// Read the variable
    status = NF90_GET_VAR( ncid, varID, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':get_var')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_3d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine inq_netcdf_var(filename, varname, error)
    !// --------------------------------------------------------------------
    !// Inquires if a variable exist on file. Return error code.
    !//
    !// Amund Sovde Haslerud, August 2017
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*), intent(in) :: filename, varname
    !// Output
    integer, intent(out)         :: error

    !// Local variables
    integer :: status, ncID, varID
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='inq_nc_var'
    !// --------------------------------------------------------------------

    !// Open file
    status = nf90_open(trim(filename), NF90_NoWrite, ncid)
    if (status .ne. NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//': open')

    !// Inquire variable id
    error = nf90_inq_varid(ncid, trim(varname), varID)
    if (.not. (status .eq. NF90_NOERR .or. status .eq. NF90_ENOTVAR)) &
         call handle_error(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//': inq_varid')

    !// Close
    status = nf90_close(ncid)
    if (status .ne. NF90_NOERR) call handle_error(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//': close')

    !// --------------------------------------------------------------------
  end subroutine inq_netcdf_var
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_4d(filename,varname,data, nlon,nlat,nlev,ntime)
    !// --------------------------------------------------------------------
    !// Reads 4D variable (varname) from netCDF file.
    !//
    !// Ole Amund Sovde, April 2016
    !// --------------------------------------------------------------------
    use cmn_precision, only: r8
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    integer, intent(in)           :: nlon,nlat,nlev,ntime
    !// Output
    real(r8),intent(out)            :: data(nlon,nlat,nlev,ntime)

    !// Local variables
    integer :: status, ncID, varID
    integer :: dimIDs(4), dimLen(4)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_4d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':inq_varid')

    !// Variable information
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':variable')

    !// 1st Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(1),len=dimLen(1))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 1')

    !// 2nd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(2),len=dimLen(2))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 2')

    !// 3rd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(3),len=dimLen(3))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 3')

    !// 4rd Dimension of the variable
    status = NF90_INQUIRE_DIMENSION(ncid,dimIDs(4),len=dimLen(4))
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':dimension 4')

    !// Read the variable
    status = NF90_GET_VAR( ncid, varID, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':get_var')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, &
         subr//':'//trim(filename)//':'//TRIM(varname)//':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_4d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_att_char( filename, varname, attname, data )
    !// --------------------------------------------------------------------
    !// Reads attribute (attname) data for variable (varname) from netCDF file.
    !// From Chris D. Holmes, UCI
    !//
    !// Ole Amund Sovde, March 2012
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname, attname
    !// Output
    character(len=*),intent(out)  :: data

    !// Local variables
    integer :: status, ncID, varID, dimLen
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_att_char'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status,subr // ':open')

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status,subr // ':inq_varid')

    !// Read the attribute
    status = NF90_GET_ATT( ncid, varID, attname, data )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status,subr // ':get_att')

    !// Trim character string
    data = TRIM( data )

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, subr // ':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_att_char
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_var_dims( filename, varname, dims )
    !// --------------------------------------------------------------------
    !// Reads dimensions of variable (varname) from netCDF file.
    !//
    !// Ole Amund Sovde, July 2015
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    !// Output
    integer,allocatable,intent(out):: dims(:)

    !// Local variables
    integer :: status, ncID, varID, I, ND
    integer, parameter :: maxdims = 20 ! should be more than enough
    integer, dimension(maxdims) :: dimIDs, dimlen
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_dims'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr // ':open '//TRIM(filename))

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    if (status /= NF90_NOERR) call handle_error(status,subr // ':inq_varid')

    !// Variable information
    dimIDs(:) = 0
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    if (status /= NF90_NOERR) call handle_error(status,subr // ':variable')

    !// Dimension of the variable
    ND = 0 ! number of dimensions
    dimLen(:) = 0
    do I = 1, maxdims
       if (dimIDs(I) .gt. 0) then
          status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(I), len=dimLen(I) )
          IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status,subr//':dimension')
          ND = ND + 1
       end if
    end do

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, subr // ':close')

    !// Allocate the variable
    ALLOCATE( dims(ND), STAT=status )
    IF (status < 0 ) then
       write(6,'(a)') f90file//':'//subr//': data allocation failed'
       stop
    end if

    !// Return the dimensions
    dims(1:ND) = dimLen(1:ND)

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_var_dims
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_dim( filename, dimname, dimLen)
    !// --------------------------------------------------------------------
    !// Finds the size of a given dimension on netCDF file.
    !//
    !// Amund Sovde Haslerud, November 2017
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, dimname
    !// Output
    integer, intent(out):: dimLen

    !// Local variables
    integer :: status, ncID
    integer  :: dimID
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr='get_netcdf_var_dims'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr // ':open '//TRIM(filename))

    !// Locate dimension
    status = NF90_INQ_DIMID( ncid, TRIM(dimname), dimID )
    if (status /= NF90_NOERR) call handle_error(status,subr // ':inq_dimid')

    !// Dimension length
    status = NF90_INQUIRE_DIMENSION( ncid, dimID, len=dimLen )
    if (status /= NF90_NOERR) call handle_error(status,subr // ':inq_dimension')

    !// Close
    status = NF90_CLOSE( ncid )
    IF (status /= NF90_NOERR) CALL HANDLE_ERROR(status, subr // ':close')

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_dim
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_r4var_1d_from_2d(filename, varname, data, &
       n1, n2, nstep)
    !// --------------------------------------------------------------------
    !// Reads netCDF4 file.
    !// Retrieves 1D variable varname(n1) from a 2D field of size (n1,n2).
    !// The n2 step defined by nstep is returned.
    !//
    !// Ole Amund Sovde, July 2015
    !// --------------------------------------------------------------------
    use cmn_precision, only: r4
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    integer, intent(in)           :: n1, n2, nstep
    !// Output
    real(r4),intent(out)          :: data(n1)

    !// Local variables
    integer :: status, ncID, varID, I, ND
    integer, parameter :: maxdims = 20 ! should be more than enough
    integer, dimension(maxdims) :: dimIDs, dimLen
    !// --------------------------------------------------------------------
    character(len=*),parameter :: subr='get_netcdf_var_1d_from_2d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':open: '//trim(filename))

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':inq_varid: '//trim(varname))

    !// Variable information
    dimIDs(:) = 0
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':variable: '//trim(varname))

    !// Find the dimensions n1 and n2
    ND = 0 ! number of dimensions
    dimLen(:) = 0
    do I = 1, maxdims
       if (dimIDs(I) .gt. 0) then
          status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(I), len=dimLen(I) )
          if (status /= NF90_NOERR) &
               call handle_error(status,subr//':dimension: '//trim(varname))
          ND = ND + 1
       end if
    end do

    if (ND .ne. 2) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong dimension on field: '//trim(varname)
       write(6,'(a,i6)') ' Expected ND=2, found: ',ND
       stop
    end if

    !// Check dimensions of field against expected (n1, n2)
    if (dimLen(1) .ne. n1 .or. (dimLen(2) .ne. n2)) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong dimension on field: '//trim(varname)
       write(6,'(a,2i6)') ' Expected/found n1: ',n1,dimLen(1)
       write(6,'(a,2i6)') ' Expected/found n2: ',n2,dimLen(2)
       stop
    end if

    !// Read the variable for nstep
    if (nstep .gt. n2) then
       write(6,'(a,i7)') f90file//':'//subr//': nstep > n2',nstep
       stop
    end if
    status = NF90_GET_VAR( ncid, varID, data, &
         start = (/1, nstep/), count=(/n1, 1 /) )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':get_var: '//trim(varname))

    !// Close
    status = NF90_CLOSE( ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':close: '//trim(filename))

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_r4var_1d_from_2d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_r4var_2d_from_3d(filename, varname, data, &
       n1, n2, n3, nstep)
    !// --------------------------------------------------------------------
    !// Reads netCDF4 file.
    !// Retrieves 2D variable varname(n1,n2) from a 3D field of size
    !// (n1,n2,n3).
    !// The n3 step defined by nstep is returned.
    !//
    !// Ole Amund Sovde, July 2015
    !// --------------------------------------------------------------------
    use cmn_precision, only: r4
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    integer, intent(in)           :: n1, n2, n3, nstep
    !// Output
    real(r4),intent(out)            :: data(n1,n2)

    !// Local variables
    integer :: status, ncID, varID, I, ND
    integer, parameter :: maxdims = 20 ! should be more than enough
    integer, dimension(maxdims) :: dimIDs, dimLen
    !// --------------------------------------------------------------------
    character(len=*),parameter :: subr='get_netcdf_var_2d_from_3d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':open: '//trim(filename))

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':inq_varid: '//trim(varname))

    !// Variable information
    dimIDs(:) = 0
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':variable: '//trim(varname))

    !// Find the dimensions
    ND = 0 ! number of dimensions
    dimLen(:) = 0
    do I = 1, maxdims
       if (dimIDs(I) .gt. 0) then
          status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(I), len=dimLen(I) )
          if (status /= NF90_NOERR) &
               call handle_error(status,subr//':dimension: '//trim(varname))
          ND = ND + 1
       end if
    end do

    if (ND .ne. 3) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong dimension on field: '//trim(varname)
       write(6,'(a,i6)') ' Expected ND=3, found: ',ND
       stop
    end if

    !// Check dimensions of field against expected (n1, n2, n3)
    if (dimLen(1) .ne. n1 .or. (dimLen(2) .ne. n2) &
         .or. (dimLen(3) .ne. n3)) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong size on field: '//trim(varname)
       write(6,'(a,2i6)') ' Expected/found n1: ',n1,dimLen(1)
       write(6,'(a,2i6)') ' Expected/found n2: ',n2,dimLen(2)
       write(6,'(a,2i6)') ' Expected/found n3: ',n3,dimLen(3)
       stop
    end if

    !// Read the variable for nstep
    if (nstep .gt. n3) then
       write(6,'(a,i7)') f90file//':'//subr//': nstep > n3',nstep
       stop
    end if
    status = NF90_GET_VAR( ncid, varID, data, &
         start = (/1, 1, nstep/), count=(/n1, n2, 1 /) )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':get_var: '//trim(varname))

    !// Close
    status = NF90_CLOSE( ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':close: '//trim(filename))

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_r4var_2d_from_3d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_netcdf_r4var_3d_from_4d( filename, varname, data, &
       n1, n2, n3, n4, nstep )
    !// --------------------------------------------------------------------
    !// Reads netCDF4 file.
    !// Retrieves 3D variable varname(n1,n2,n3) from a 4D field of size
    !// (n1,n2,n3,n4).
    !// The n4 step defined by nstep is returned.
    !//
    !// Ole Amund Sovde, July 2015
    !// --------------------------------------------------------------------
    use cmn_precision, only: r4
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in)   :: filename, varname
    integer, intent(in)           :: n1, n2, n3, n4, nstep
    !// Output
    real(r4),intent(out)            :: data(n1,n2,n3)

    !// Local variables
    integer :: status, ncID, varID, I, ND
    integer, parameter :: maxdims = 20 ! should be more than enough
    integer, dimension(maxdims) :: dimIDs, dimLen
    !// --------------------------------------------------------------------
    character(len=*),parameter :: subr='get_netcdf_var_3d_from_4d'
    !// --------------------------------------------------------------------

    !// Open file
    status = NF90_OPEN( TRIM(filename), NF90_NoWrite, ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':open: '//trim(filename))

    !// Locate variable
    status = NF90_INQ_VARID( ncid, TRIM(varname), varID )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':inq_varid: '//trim(varname))

    !// Variable information
    dimIDs(:) = 0
    status = NF90_INQUIRE_VARIABLE( ncid, varID, dimIDs=dimIDs )
    if (status /= NF90_NOERR) &
         call handle_error(status,subr//':variable: '//trim(varname))

    !// Find the dimensions
    ND = 0 ! number of dimensions
    dimLen(:) = 0
    do I = 1, maxdims
       if (dimIDs(I) .gt. 0) then
          status = NF90_INQUIRE_DIMENSION( ncid, dimIDs(I), len=dimLen(I) )
          if (status /= NF90_NOERR) &
               call handle_error(status,subr//':dimension: '//trim(varname))
          ND = ND + 1
       end if
    end do

    if (ND .ne. 4) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong dimension on field: '//trim(varname)
       write(6,'(a,i6)') ' Expected ND=4, found: ',ND
       stop
    end if

    !// Check dimensions of field against expected (n1, n2, n3)
    if (dimLen(1) .ne. n1 .or. (dimLen(2) .ne. n2) &
         .or. (dimLen(3) .ne. n3) .or. (dimLen(4) .ne. n4)) then
       write(6,'(a)') f90file//':'//subr// &
            ': Wrong size on field: '//trim(varname)
       write(6,'(a,2i6)') ' Expected/found n1: ',n1,dimLen(1)
       write(6,'(a,2i6)') ' Expected/found n2: ',n2,dimLen(2)
       write(6,'(a,2i6)') ' Expected/found n3: ',n3,dimLen(3)
       write(6,'(a,2i6)') ' Expected/found n4: ',n4,dimLen(4)
       stop
    end if

    !// Read the variable for nstep
    if (nstep .gt. n4) then
       write(6,'(a,i7)') f90file//':'//subr//': nstep > n4',nstep
       stop
    end if
    status = NF90_GET_VAR( ncid, varID, data, &
         start = (/1, 1, 1, nstep/), count=(/n1, n2, n3, 1 /) )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':get_var: '//trim(varname))

    !// Close
    status = NF90_CLOSE( ncid )
    if (status /= NF90_NOERR) &
         call handle_error(status, subr//':close: '//trim(filename))

    !// --------------------------------------------------------------------
  end subroutine get_netcdf_r4var_3d_from_4d
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine handle_err(status)
    !// --------------------------------------------------------------------
    !// Handle netCDF errors. This routine does not print out much
    !// helpful information.
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    integer, intent(in)   :: status  !Error status
    !// --------------------------------------------------------------------
    character(len=*),parameter :: subr='handle_err'
    !// --------------------------------------------------------------------
    if (status/=nf90_noerr) then
       write(6,'(a,i7)') f90file//':'//subr// &
            ': ERROR: '//trim(nf90_strerror(status))
       stop "stopped"
    end if
    !// --------------------------------------------------------------------
  end subroutine handle_err
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine handle_error(ECODE,EMSG)
    !// --------------------------------------------------------------------
    !// Handle netCDF errors slightly better than only handle_err.
    !// Ole Amund Sovde, March 2011
    !// --------------------------------------------------------------------
    use netcdf
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    integer, intent(in) :: ECODE
    character(len=*) :: EMSG
    !// --------------------------------------------------------------------
    character(len=*),parameter :: subr='handle_error'
    !// --------------------------------------------------------------------
    if (ECODE /= nf90_noerr) then
       write(6,'(a)') f90file//':'//subr//': netCDF error!'
       write(*,'(a,i5)') '  Error code: ',ECODE
       write(*,'(a)')    '  Error message: '//EMSG
       write(*,'(a)')    '  nf90_strerror: '//trim(nf90_strerror(ECODE))
       stop "STOPPED"
    end if
    !// --------------------------------------------------------------------
  end subroutine handle_error
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
end module ncutils
!//=========================================================================


!!===================================================================
!!
!! Module:
!!
!!    misc          Last Modified by Meng Wu (2020)
!!
!! Routines:
!!
!! 1. checknorm()   By Meng Wu (2020)
!!
!!    Check the normalization condition of Bloch waves for one kpoint, one band
!!
!! 2. get_volume()       Originally By (SIB)     Last Modified 6/12/2008 (JRD)
!!
!!    This assumes that b is a symmetric matrix.  It sets
!!    vol = (2pi)^3 / square_root(|det(b)|)
!!    This makes sense if b is the matrix of dot products of the recip
!!    lattice vectors, so vol is the real space volume.
!!
!! 3. findvector()       Originally By (SIB)     Last Modified 6/12/2008 (JRD)
!!
!!    Given Gvectors as a 3D vector, return its index in the gvec%components(:,:) array
!!    If not found, return 0
!!
!! 4. M33INV()
!!
!!    Inverse of a 3x3 matrix
!!
!! 5. procmem()          Originally By (gsm)     Last Modified 4/14/2009 (gsm)
!!
!!    Determines the amount of free memory per processor
!!    from the proc file system
!!
!! 6. sizeof_scalar()    Originally By (DAS)     Last Modified 1/25/2011 (DAS)
!!
!!    Return the size of the SCALAR type, for memory estimation.
!!
!! 7. get_gumk3()       By Meng Wu (2020)
!!
!!    Given a kpoint (could be outside 1st BZ), return its Umklapp vector
!!
!! 8. get_gumk3()       By Meng Wu (2020)
!!
!!    Given a kpoint (could be outside 1st BZ), return its Umklapp vector
!!
!! 9. generate_full_WS_BZ()   By Meng Wu (2020)
!!
!!    Given kgrid, generate all the kpoints in the Wigner-Seitz BZ
!!
!! 10. prepare symmetry()     By Meng Wu (2020)
!!
!!    Generate point group symmetry matrices in the cartesian
!!    basis and reciprocal lattice basis
!!
!! 11. subgroup()             By Meng Wu (2020)
!!
!!    Generate a subgroup of syms that preserves dq(1:3)
!!
!!===================================================================

#include "f_defs.h"

module misc_m
  use global_m
  use blas_m
  use scalapack_m
  implicit none
  private
  public ::   checknorm, get_volume, findvector, procmem, sizeof_scalar, &
       get_gumk3, M33INV, generate_full_WS_BZ, prepare_syms, subgroup
contains

  subroutine checknorm(filename, iband, ik, nspin, wfn, ispin, tol)
    character (len=*), intent(in) :: filename
    integer, intent(in) :: iband,ik, nspin
    SCALAR, intent(in) :: wfn(:,:) !< (ng,nspin*nspinor)
    integer, optional, intent(in) :: ispin
    real(DP), optional, intent(in) :: tol
    real(DP) :: tol_
    SCALAR :: norm
    integer :: wfnsize
    SCALAR, allocatable :: wfn_1d(:)

    PUSH_SUB(checknorm)

    if (present(tol)) then
       tol_ = tol
    else
       tol_ = TOL_SMALL
    endif
    !> wfnsize = ng * nspinor
    wfnsize = SIZE(wfn)/nspin
    SAFE_ALLOCATE(wfn_1d, (wfnsize))

    if (nspin .eq. 1) then
       !> Do a dot product, make sure that the result has 1.0D0 as real part and 0.0D0 as imaginary part
       !> DOT(CONJG(wfn),wfn)

       !> Flatten the 2D wfn into 1D wfn_1d
       wfn_1d = PACK(wfn, .true.)
       norm = blas_dot(wfnsize,wfn_1d,1,wfn_1d,1)

       if (ABS(norm - 1.0D0) > tol_) then
          write(0,555) TRUNC(filename),ABS(norm-1.0D0),iband,1,ik
          call die("Incorrect normalization.")
       endif
    elseif (nspin .eq. 2) then
       if (present(ispin)) then
          if (ispin .eq. 1) then
             wfn_1d = PACK(wfn(:,1), .true.)
             norm = blas_dot(wfnsize,wfn_1d,1,wfn_1d,1)

             if (ABS(norm - 1.0D0) > tol_) then
                write(0,555) TRUNC(filename),ABS(norm-1.0D0),iband,1,ik
                call die("Incorrect normalization.")
             endif
          elseif (ispin .eq. 2) then
             wfn_1d = PACK(wfn(:,2), .true.)
             norm = blas_dot(wfnsize,wfn_1d,1,wfn_1d,1)

             if (ABS(norm - 1.0D0) > tol_) then
                write(0,555) TRUNC(filename),ABS(norm-1.0D0),iband,2,ik
                call die("Incorrect normalization.")
             endif
          else
             call die("wrong input of ispin.")
          endif
       else
          wfn_1d = PACK(wfn(:,1), .true.)
          norm = blas_dot(wfnsize,wfn_1d,1,wfn_1d,1)

          if (ABS(norm - 1.0D0) > tol_) then
             write(0,555) TRUNC(filename),ABS(norm-1.0D0),iband,1,ik
             call die("Incorrect normalization.")
          endif

          wfn_1d = PACK(wfn(:,2), .true.)
          norm = blas_dot(wfnsize,wfn_1d,1,wfn_1d,1)

          if (ABS(norm - 1.0D0) > tol_) then
             write(0,555) TRUNC(filename),ABS(norm-1.0D0),iband,2,ik
             call die("Incorrect normalization.")
          endif
       endif
    else
       call die("checknorm_2: Check nspin",only_root_writes=.true.)
    endif

    SAFE_DEALLOCATE(wfn_1d)

555 format(1x,'Wavefunction is not normalized in file',1x,a,/,&
         3x,'abs(norm - 1) =',f10.7,/,&
         3x,'iband =',i6,1x,'ispin =',i2,1x,'ik =',i6,/)

    POP_SUB(checknorm)

    return
  end subroutine checknorm

  subroutine get_volume(vol,b)
    real(DP), intent(out) :: vol
    real(DP), intent(in)  :: b(3,3)

    PUSH_SUB(get_volume)

    vol = b(1,1)*(b(2,2)*b(3,3) - b(2,3)**2) &
         + 2*b(1,2)*b(2,3)*b(3,1) &
         - b(2,2)*b(1,3)**2 - b(3,3)*b(1,2)**2
    vol = sqrt(abs(vol))
    vol = ((2.0d0*PI_D)**3)/vol

    POP_SUB(get_volume)

    return
  end subroutine get_volume

  !=====================================================================

  subroutine findvector(iout,kk,gvec)
    integer, intent(out) :: iout
    integer, intent(in)  :: kk(3)
    type (gspace), intent(in) :: gvec

    ! no push/pop since called too frequently

    iout=((kk(1)+gvec%FFTgrid(1)/2)*gvec%FFTgrid(2)+kk(2)+gvec%FFTgrid(2)/2)* &
         gvec%FFTgrid(3)+kk(3)+gvec%FFTgrid(3)/2+1
    if (iout .ge. 1 .and. iout .le. gvec%nFFTgridpts) then
       iout=gvec%index_vec(iout)
       if (iout .ge. 1 .and. iout .le. gvec%ng) then
          if (any(kk(1:3) /= gvec%components(1:3, iout))) iout = 0
       else
          iout = 0
       endif
    else
       iout = 0
    endif

    return
  end subroutine findvector

  SUBROUTINE M33INV (A, AINV, OK_FLAG)

    IMPLICIT NONE

    DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN)  :: A
    DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: AINV
    LOGICAL, INTENT(OUT) :: OK_FLAG

    DOUBLE PRECISION, PARAMETER :: EPS = 1.0D-10
    DOUBLE PRECISION :: DET
    DOUBLE PRECISION, DIMENSION(3,3) :: COFACTOR

    DET =   A(1,1)*A(2,2)*A(3,3)  &
         - A(1,1)*A(2,3)*A(3,2)  &
         - A(1,2)*A(2,1)*A(3,3)  &
         + A(1,2)*A(2,3)*A(3,1)  &
         + A(1,3)*A(2,1)*A(3,2)  &
         - A(1,3)*A(2,2)*A(3,1)

    IF (ABS(DET) .LE. EPS) THEN
       AINV = 0.0D0
       OK_FLAG = .FALSE.
       RETURN
    END IF

    COFACTOR(1,1) = +(A(2,2)*A(3,3)-A(2,3)*A(3,2))
    COFACTOR(1,2) = -(A(2,1)*A(3,3)-A(2,3)*A(3,1))
    COFACTOR(1,3) = +(A(2,1)*A(3,2)-A(2,2)*A(3,1))
    COFACTOR(2,1) = -(A(1,2)*A(3,3)-A(1,3)*A(3,2))
    COFACTOR(2,2) = +(A(1,1)*A(3,3)-A(1,3)*A(3,1))
    COFACTOR(2,3) = -(A(1,1)*A(3,2)-A(1,2)*A(3,1))
    COFACTOR(3,1) = +(A(1,2)*A(2,3)-A(1,3)*A(2,2))
    COFACTOR(3,2) = -(A(1,1)*A(2,3)-A(1,3)*A(2,1))
    COFACTOR(3,3) = +(A(1,1)*A(2,2)-A(1,2)*A(2,1))

    AINV = TRANSPOSE(COFACTOR) / DET

    OK_FLAG = .TRUE.

    RETURN

  END SUBROUTINE M33INV

  ! subroutine procmem(mem, ranks_per_node, nfreq_group)
  subroutine procmem(mem, ranks_per_node)
    real(DP), intent(out) :: mem !< Memory per MPI rank (lower bound)
    integer, intent(out) :: ranks_per_node !< # of MPI ranks per node (upper bound)
    ! integer, intent(in), optional :: nfreq_group

    integer, parameter :: host_len=64
    integer :: ierr,eof,info,iunit,m,n,i,pagesize
    real(DP) :: x,y,mac_m,mac_n
    !    integer :: ntot
    !    real(DP) :: xtot,ytot ! we do not use the total memory actually
    character(len=256) :: s, filename, my_hostname
    character(len=host_len), allocatable :: hostnames(:)

    PUSH_SUB(procmem)

    !-----------------------------------------------------
    !> determine the amount of free memory per node in kB

    m=0
    iunit=14
    x=0

    call open_file(unit=iunit,file='/proc/meminfo',form='formatted',iostat=ierr,status='old')
    if (ierr.eq.0) then
       eof=0
       do while(eof.eq.0)
          read(iunit,'(a)',iostat=eof)s
          if (s(1:7).eq."MemFree") then
             read(s(9:),*)n
             m=m+n
          endif
          !        if (s(1:8).eq."MemTotal") then
          !          read(s(10:),*)ntot
          !        endif
          if (s(1:6).eq."Cached") then
             read(s(8:),*)n
             m=m+n
          endif
       enddo
       x=dble(m)/dble(peinf%npes)
       call close_file(iunit)
    endif

    if(m == 0) then
       !> this is for Mac OS
       !! total memory is accessible instead from sysctl -n hw.usermem
       write(filename,'(a,i9.9)') 'vm_stat_', peinf%inode
       SYSTEMCALL("vm_stat > " + TRUNC(filename) + " 2> /dev/null")
       !> Fortran 2008 would use execute_command_line instead
       !! even if the command failed, still open file in order to delete it
       call open_file(unit=iunit,file=TRUNC(filename),form='formatted',iostat=ierr,status='old')
       if (ierr.eq.0) then
          eof=0
          do while(eof.eq.0)
             read(iunit,'(a)',iostat=eof)s
             if (s(1:45).eq."Mach Virtual Memory Statistics: (page size of") then
                read(s(46:),*)pagesize ! in bytes
             endif
             if (s(1:11).eq."Pages free:") then
                read(s(12:),*) mac_n
                mac_m = mac_m + mac_n
             endif
             if (s(1:18).eq."Pages speculative:") then
                read(s(19:),*) mac_n
                mac_m = mac_m + mac_n
             endif
          enddo
          call close_file(iunit, delete = .true.)
          x = mac_m * dble(pagesize) / dble(peinf%npes * 1024) ! to kB
       endif
    endif

    if(m == 0 .and. mac_m == 0) then ! BSD
       write(filename,'(a,i9.9)') 'sysctl_', peinf%inode
       SYSTEMCALL("sysctl -a > " + TRUNC(filename) + " 2> /dev/null")
       !> Fortran 2008 would use execute_command_line instead
       !! even if the command failed, still open file in order to delete it
       call open_file(unit=iunit,file=TRUNC(filename),form='formatted',iostat=ierr,status='old')
       if (ierr.eq.0) then
          eof=0
          do while(eof.eq.0)
             read(iunit,'(a)',iostat=eof)s
             if (s(1:12).eq."hw.pagesize:") then
                read(s(13:),*)pagesize ! in bytes
             endif
             if (s(1:25).eq."vm.stats.vm.v_free_count:") then
                read(s(26:),*) mac_n
                mac_m = mac_m + mac_n
             endif
             if (s(1:26).eq."vm.stats.vm.v_cache_count:") then
                read(s(27:),*) mac_n
                mac_m = mac_m + mac_n
             endif
          enddo
          call close_file(iunit, delete = .true.)
          x = mac_m * dble(pagesize) / dble(peinf%npes * 1024) ! to kB
       endif
    endif

    !    xtot=dble(ntot)/dble(peinf%npes)

#ifdef MPI
    call MPI_Allreduce(x,y,1,MPI_REAL_DP,MPI_SUM,MPI_COMM_WORLD,mpierr)
#else
    y=x
    !    ytot=xtot
#endif

    !----------------------------------------------
    ! Determine the number of processors per node
    ! FHJ: we do this in parallel. Each MPI rank compares with all other hostnames,
    ! and figures out the number of ranks running on the same node. At the end, we
    ! pick the upper bound for ranks_per_node.

    SAFE_ALLOCATE(hostnames, (peinf%npes))
    HOSTNAMECALL(my_hostname, info)
    hostnames(peinf%inode+1) = my_hostname(1:host_len)

#ifdef MPI
    call MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, hostnames(1), &
         host_len, MPI_BYTE, MPI_COMM_WORLD, mpierr)
    ranks_per_node = 0
    do i = 1, peinf%npes
       ! FHJ: no need to trim, all strings have the same length
       if (hostnames(peinf%inode+1)==hostnames(i)) ranks_per_node = ranks_per_node + 1
    enddo
    call MPI_Allreduce(MPI_IN_PLACE, ranks_per_node, 1, MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD, mpierr)
#else
    ranks_per_node = 1
#endif

    !-----------------------------------
    !> report the available memory in B

    if (ranks_per_node>1) y = y / dble(ranks_per_node)
    mem=y*1024.0d0

    !-----------------------------------
    !> warn if zero memory

    if (mem .lt. TOL_Small .and. peinf%inode .eq. 0) then
       write(0,666)
666    format(1x,'WARNING: estimation of memory available is zero, probably failed.',/)
    endif

    POP_SUB(procmem)

    return
  end subroutine procmem

  !================================================================================

  !> for memory estimation, tell what size of SCALAR type is
  integer function sizeof_scalar()

#ifdef NOSIZEOF
    !> if no sizeof, take a guess

#ifdef CPLX
    sizeof_scalar = 16
#else
    sizeof_scalar = 8
#endif

#else

    SCALAR :: dummy
    sizeof_scalar = sizeof(dummy)

#endif

  end function sizeof_scalar

  !> (kvector(:) - dble(gumk_out(:))) falls into first WS BZ, having the shortest distance to Gamma
  subroutine get_gumk3(bdot, kvector, gumk_out)
    real(DP), intent(in) :: bdot(3,3)
    real(DP), intent(in) :: kvector(3)
    integer, intent(out) :: gumk_out(3)
    real(DP) :: len2_kvector_temp, len2_kvector_temp_min, kvector_temp(3), kvector_(3)
    integer :: b1, b2, b3, ib1, ib2, ib3, b_out(3), b1list(7), b2list(7), b3list(7), nb1, nb2, nb3, nb(3), id
    logical :: dim_flag(3)

    !> If kvector(id) = 0 and bdot(id,id) is block-diagonal, e.g. bdot(1,3)=DOT[b1,b3] = 0, bdot(2,3)=DOT[b2,b3]=0
    !> do not apply integer shift!
    if ( ABS(kvector(1)) < TOL_ZERO .and. ABS(bdot(1,2)) < TOL_ZERO .and. ABS(bdot(1,3)) < TOL_ZERO ) then
       nb1 = 1
       b1list = 0
    else
       nb1 = 7
       b1list(:) = (/0, -1, 1, -2, 2, -3, 3/)
    endif

    if ( ABS(kvector(2)) < TOL_ZERO .and. ABS(bdot(1,2)) < TOL_ZERO .and. ABS(bdot(2,3)) < TOL_ZERO ) then
       nb2 = 1
       b2list = 0
    else
       nb2 = 7
       b2list(:) = (/0, -1, 1, -2, 2, -3, 3/)
    endif

    if ( ABS(kvector(3)) < TOL_ZERO .and. ABS(bdot(1,3)) < TOL_ZERO .and. ABS(bdot(2,3)) < TOL_ZERO ) then
       nb3 = 1
       b3list = 0
    else
       nb3 = 7
       b3list(:) = (/0, -1, 1, -2, 2, -3, 3/)
    endif

    len2_kvector_temp_min = INF
    kvector_ = 0.0D0
    b_out = 0

    nb = (/ nb1, nb2, nb3 /)
    dim_flag = (nb > 1)

    do ib1 = 1, nb1
       b1 = b1list(ib1)
       do ib2 = 1, nb2
          b2 = b2list(ib2)
          do ib3 = 1, nb3
             b3 = b3list(ib3)

             kvector_temp(:) = kvector(:) - (/DBLE(b1), DBLE(b2), DBLE(b3)/)
             len2_kvector_temp = DOT_PRODUCT(kvector_temp, MATMUL(bdot, kvector_temp))

             !> If T, it means the kvector is at the edge of BZ
             if (ABS(len2_kvector_temp - len2_kvector_temp_min) < TOL_SMALL) then
                !> Determine whether we should replace b_out
                !> We try to make the first meaningful fractional coordinates >= 0 and as large as possible
                !> Since BZ has inversion symmetry, we can always find an BZ edge point with the first meaningful fractional coordinate >= 0
                loop_id: do id = 1, 3
                   if (.not. dim_flag(id)) cycle
                   !> kvector_temp(id)
                   !> Pick the larger fractional coordinate
                   if (kvector_temp(id) > kvector_(id)) then
                      len2_kvector_temp_min = len2_kvector_temp
                      kvector_(:) = kvector_temp(:)
                      b_out(:) = (/ b1, b2, b3 /)
                      !> If this fractional coordinate is the same, cycle to the next one
                   elseif ( ABS(kvector_temp(id) - kvector_(id)) < TOL_ZERO) then
                      cycle
                   endif
                   exit loop_id
                enddo loop_id
             elseif ( len2_kvector_temp < len2_kvector_temp_min ) then
                !> Determine whether this is at edge of BZ
                len2_kvector_temp_min = len2_kvector_temp
                kvector_(:) = kvector_temp(:)
                b_out(:) = (/ b1, b2, b3 /)
             endif
          enddo
       enddo
    enddo
    gumk_out(:) = b_out(:)

  end subroutine get_gumk3

  subroutine generate_full_WS_BZ(kgrid, qshift, bdot, nfk, fk)
    integer, intent(in) :: kgrid(3)
    real(DP), intent(in) :: qshift(3) !> in reciprocal space crystal coordinates [b1,b2,b3]
    real(DP), intent(in) :: bdot(3,3)
    integer, intent(out) :: nfk
    real(DP), allocatable, intent(out) :: fk(:,:)
    integer :: ifk, gumk(3), i1, i2, i3

    nfk = PRODUCT(kgrid(:))
    if (nfk .le. 0) then
       call die("kgrid mistake.", only_root_writes=.true.)
    endif
    SAFE_ALLOCATE(fk, (3,nfk))

    !> Generate kmesh + qshift points
    ifk = 0
    do i1 = 0, kgrid(1)-1
       do i2 = 0, kgrid(2)-1
          do i3 = 0, kgrid(3)-1
             ifk = ifk + 1
             !> kpoint in crystal coordinates [0,1]
             fk(1,ifk) = dble(i1) / dble(kgrid(1)) + qshift(1)
             fk(2,ifk) = dble(i2) / dble(kgrid(2)) + qshift(2)
             fk(3,ifk) = dble(i3) / dble(kgrid(3)) + qshift(3)
             !> Move kpoint to Wigner-Seitz BZ
             call get_gumk3(bdot, fk(:,ifk), gumk)
             !> move kpoint to 1st WS BZ
             fk(:,ifk) = fk(:,ifk) - dble(gumk(:))
          enddo
       enddo
    enddo

    if (ifk .ne. nfk) then
       call die("nfk mistach.", only_root_writes=.true.)
    endif
  end subroutine generate_full_WS_BZ

  subroutine prepare_syms(crys, syms)
    type(crystal), intent(in) :: crys
    type(symmetry), intent(inout) :: syms
    real(DP), dimension(3,3) :: atTat
    real(DP), dimension(3,3) :: inv_atTat
    real(DP), dimension(3,3) :: inv_at
    integer :: itran
    logical :: flag_inv

    ! crys%avec(:,:) = at(:,:)
    ! at(:,:) = [a1, a2, a3]
    atTat = MATMUL(TRANSPOSE(crys%avec),crys%avec)
    ! -------------------------------
    ! [a1]
    ! ------adot = [a1 a2 a3] \cdot [a2]
    ! [a3]
    ! -------------------------------
    call M33INV(atTat, inv_atTat, flag_inv)
    if (.not. flag_inv) then
       call die("prepare_syms: atTat not invertable", only_root_writes=.true.)
    endif
    call M33INV(crys%avec,inv_at, flag_inv)
    if (.not. flag_inv) then
       call die("prepare_syms: crys%avec not invertable", only_root_writes=.true.)
    endif
    syms%mtrx_reci(:,:,:) = 0.0D0
    syms%mtrx_cart(:,:,:) = 0.0D0
    do itran = 1, syms%ntran
       syms%mtrx_reci(1:3,1:3,itran) = nint(MATMUL(MATMUL(atTat, dble(syms%mtrx(1:3,1:3,itran))), inv_atTat))
       syms%mtrx_cart(1:3,1:3,itran) = MATMUL(MATMUL(crys%avec, dble(syms%mtrx(1:3,1:3,itran))), inv_at)
    enddo

  end subroutine prepare_syms

  !> Generate a subgroup of syms that preserves dq(1:3)
  subroutine subgroup(syms_sub, syms, crys, qshift)
    type (symmetry), intent(inout) :: syms_sub
    type (crystal), intent(in) :: crys
    type (symmetry), intent(in) :: syms
    real(DP), intent(in) :: qshift(3)
    integer :: gumk(3), itran
    real(DP) :: dq(3)

    syms_sub%ntran = 0
    syms_sub%mtrx(:,:,:) = 0
    syms_sub%mtrx_reci(:,:,:) = 0.0D0
    syms_sub%mtrx_cart(:,:,:) = 0.0D0
    syms_sub%kgzero(:,:) = 0.0D0
    syms_sub%tnp(:,:) = 0.0D0

    do itran = 1, syms%ntran
       dq(1:3) = MATMUL(DBLE(syms%mtrx_reci(1:3, 1:3, itran)), qshift(1:3)) - qshift(1:3)
       call get_gumk3(crys%bdot, dq, gumk)
       dq(:) = dq(:) - DBLE(gumk(:))
       if (all(abs(dq(1:3)) .lt. TOL_Small)) then
          !> Store index of element of subgroup
          syms_sub%ntran = syms_sub%ntran + 1
          syms_sub%mtrx(1:3,1:3,syms_sub%ntran) = syms%mtrx(1:3,1:3,itran)
          syms_sub%mtrx_reci(1:3,1:3,syms_sub%ntran) = syms%mtrx_reci(1:3,1:3,itran)
          syms_sub%mtrx_cart(1:3,1:3,syms_sub%ntran) = syms%mtrx_cart(1:3,1:3,itran)
          syms_sub%tnp(1:3,syms_sub%ntran) = syms%tnp(1:3,itran)
          !> matmul(syms_sub%mtrx_reci(1:3, 1:3, itran), qshift(1:3)) + syms_sub%kgzero(1:3,itran) = qshift(1:3), itran = 1, syms_sub%ntran
          syms_sub%kgzero(1:3, syms_sub%ntran) = - gumk(1:3)
       endif
    enddo

  end subroutine subgroup

end module misc_m


text:
      subroutine ZZ_Hgaga_g(p,msq)
      implicit none
c--- Weak Boson Fusion by Z-Z exchange only
c---Matrix element squared averaged over initial colors and spins
c
c     q(-p1)+q(-p2) -->  H(p3,p4)+q(p5)+q(p6)+g(p7)
c                           |
c                           |
c                           |
c                           ---> b(p3)+bbar(p4)
c---- Extension to photon decay contributed by Fabian Stoeckli
      include 'constants.f'
      include 'masses.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'zcouple.f'
      include 'sprods_com.f'
      integer j,k,nup,ndo
      double precision p(mxpart,4),facqq,facqg
      double precision msq(-nf:nf,-nf:nf),hdecay,
     . ud_udg_LL,udb_udbg_LL,ud_udg_LR,udb_udbg_LR,
     . ug_uddb_LL,ug_uddb_LR,gu_ddbu_LL,gu_ddbu_LR
      double precision msqgamgam

      parameter (nup=2,ndo=3)

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo

      call dotem(7,p,s)

      hdecay=msqgamgam(hmass)/((s(3,4)-hmass**2)**2+(hmass*hwidth)**2)

      facqq=0.25d0*gsq*Cf*gwsq**3*hdecay
      facqg=-facqq*(aveqg/aveqq)
C Extra factor of gsq*Cf compared to lowest order

C q-q and qbar-qbar
c--- u(1)+d(2) -> u(5)+d(6)+g(p7)
c--- ub(1)+db(2) -> ub(5)+db(6)+g(p7)
      call msq_gpieces_gaga_zz(1,2,5,6,7,ud_udg_LL,ud_udg_LR)

C q-qbar and qbar-q
c--- u(1)+db(2) -> u(5)+db(6)+g(p7)
c--- ub(1)+d(2) -> ub(5)+d(6)+g(p7)
      call msq_gpieces_gaga_zz(1,6,5,2,7,udb_udbg_LL,udb_udbg_LR)

C q-g and qbar-g
c--- u(1)+g(2) -> u(5)+[q](6)+[qb](7)
c--- ub(1)+g(2) -> ub(5)+[qb](6)+[q](7)
      call msq_gpieces_gaga_zz(1,7,5,6,2,ug_uddb_LL,ug_uddb_LR)

C g-q and g-qbar
c--- g(1)+u(2) -> [q](5)+[qb](6)+u(7)
c--- g(1)+ub(2) -> [qb](5)+[q](6)+ub(7)
      call msq_gpieces_gaga_zz(6,2,5,7,1,gu_ddbu_LL,gu_ddbu_LR)

      do j=-nf,nf
      do k=-nf,nf
        if     ((j .gt. 0) .and. (k .lt. 0)) then
          msq(j,k)=facqq*(
     .     +udb_udbg_LL*((L(+j)*L(-k))**2+(R(+j)*R(-k))**2)
     .     +udb_udbg_LR*((L(+j)*R(-k))**2+(R(+j)*L(-k))**2))
        elseif ((j .lt. 0) .and. (k .gt. 0)) then
          msq(j,k)=facqq*(
     .     +udb_udbg_LL*((L(-j)*L(k))**2+(R(-j)*R(k))**2)
     .     +udb_udbg_LR*((L(-j)*R(k))**2+(R(-j)*L(k))**2))
        elseif ((j .gt. 0) .and. (k .gt. 0)) then
          msq(j,k)=facqq*(
     .     +ud_udg_LL*((L(+j)*L(+k))**2+(R(+j)*R(+k))**2)
     .     +ud_udg_LR*((L(+j)*R(+k))**2+(R(+j)*L(+k))**2))
        elseif ((j .lt. 0) .and. (k .lt. 0)) then
          msq(j,k)=facqq*(
     .     +ud_udg_LL*((L(-j)*L(-k))**2+(R(-j)*R(-k))**2)
     .     +ud_udg_LR*((L(-j)*R(-k))**2+(R(-j)*L(-k))**2))
        elseif ((j .ne. 0) .and. (k .eq. 0)) then
          msq(j,k)=facqg*dfloat(nup)*(
     .     +ug_uddb_LL*((L(abs(j))*L(+2))**2+(R(abs(j))*R(+2))**2)
     .     +ug_uddb_LR*((L(abs(j))*R(+2))**2+(R(abs(j))*L(+2))**2))
     .            +facqg*dfloat(ndo)*(
     .     +ug_uddb_LL*((L(abs(j))*L(+1))**2+(R(abs(j))*R(+1))**2)
     .     +ug_uddb_LR*((L(abs(j))*R(+1))**2+(R(abs(j))*L(+1))**2))
        elseif ((j .eq. 0) .and. (k .ne. 0)) then
          msq(j,k)=facqg*dfloat(nup)*(
     .     +gu_ddbu_LL*((L(abs(k))*L(+2))**2+(R(abs(k))*R(+2))**2)
     .     +gu_ddbu_LR*((L(abs(k))*R(+2))**2+(R(abs(k))*L(+2))**2))
     .            +facqg*dfloat(ndo)*(
     .     +gu_ddbu_LL*((L(abs(k))*L(+1))**2+(R(abs(k))*R(+1))**2)
     .     +gu_ddbu_LR*((L(abs(k))*R(+1))**2+(R(abs(k))*L(+1))**2))
        endif
      enddo
      enddo

      return
      end

      subroutine msq_gpieces_gaga_zz(i1,i2,i5,i6,i7,zll,zlr)
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'ewcouple.f'
      double precision zll,zlr,htheta
      double precision msqll_a,msqll_b,msqlr_a,msqlr_b,
     . msq_gsamehel,msq_gopphel
      double precision propz,x
      double precision s15,s26,s157,s267
      integer i1,i2,i5,i6,i7
C--- define Heaviside theta function (=1 for x>0) and (0 for x < 0)
      htheta(x)=half+sign(half,x)
      propz(s15)=sign(one,(s15-zmass**2))
     . *dsqrt(dsqrt(1d0-xw)/xw/2d0/zmass
     . *((s15-zmass**2)**2+htheta(s15)*(zmass*zwidth)**2))

c--- Calculate some invariants
      s15=s(i1,i5)
      s26=s(i2,i6)
      s157=s15+s(i1,i7)+s(i5,i7)
      s267=s26+s(i2,i7)+s(i6,i7)

      msqll_a=msq_gsamehel(i1,i2,i5,i6,i7)
      msqll_b=msq_gsamehel(i2,i1,i6,i5,i7)
      msqlr_a=msq_gopphel(i1,i2,i5,i6,i7)
      msqlr_b=msq_gopphel(i2,i1,i6,i5,i7)

c--- catch the unwanted diagrams for the gluon-quark processes
      if (i7 .eq. 1) then
        msqll_b=0d0
        msqlr_b=0d0
      endif
      if (i7 .eq. 2) then
        msqll_a=0d0
        msqlr_a=0d0
      endif

      zll=msqll_a/(propz(s157)*propz(s26))**2
     .   +msqll_b/(propz(s267)*propz(s15))**2

      zlr=msqlr_a/(propz(s157)*propz(s26))**2
     .   +msqlr_b/(propz(s267)*propz(s15))**2

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C--------------------------------------------------------------------------
        subroutine ftkshf(iunit,colmin,colmax,incre,status)

C       shift the index value on any existing column keywords
C       This routine will modify the name of any keyword that begins with 'T'
C       and has an index number in the range COLMIN - COLMAX, inclusive.

C       if incre is positive, then the index values will be incremented.
C       if incre is negative, then the kewords with index = COLMIN
C       will be deleted and the index of higher numbered keywords will
C       be decremented.

C       iunit   i  Fortran I/O unit number
C       colmin  i  starting column number to be incremented
C       colmax  i  maximum column number to be increment
C       incre   i  amount by which the index value should be shifted
C       status  i  returned error status (0=ok)

        integer iunit,colmin,colmax,incre,status

C       COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nf,nb,ne
        parameter (nb = 20)
        parameter (nf = 3000)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
        integer tfield,tstart,tbcol,rowlen,tdtype,trept,tnull,scount
        integer theap,nxheap
        double precision tscale,tzero
        common/ft0002/tfield(nb),tstart(nb),tbcol(nf),rowlen(nb),
     &  tdtype(nf),trept(nf),tscale(nf),tzero(nf),tnull(nf),scount(nb)
     &  ,theap(nb),nxheap(nb)
C       END OF COMMON BLOCK DEFINITIONS-----------------------------------

        integer ibuff,typhdu,tflds,nkeys,nmore,nrec,ival,tstat,i1
        character rec*80,newkey*8,q*4

C       define the number of the buffer used for this file
        ibuff=bufnum(iunit)

C       test that the CHDU is an ASCII table or BINTABLE
        typhdu=hdutyp(ibuff)
        if (typhdu .ne. 1 .and. typhdu .ne. 2)then
                status=235
                call ftpmsg('Can only operate on TABLE or '//
     &          'BINTABLE extension (FTKSHF)')
                return
        end if

C       test column number limits
        tflds=tfield(ibuff)
        if (colmin .lt. 1 .or. colmax .lt. 1)then
             status=302
             return
        else if (colmin .gt. colmax .or. colmin .gt. tflds)then
             return
        end if

C       get the number of keywords in the header
        call ftghsp(iunit,nkeys,nmore,status)

C       go thru header starting with the 9th keyword looking for 'TxxxxNNN'

        nrec=9
100     call ftgrec(iunit,nrec,rec,status)

        if (rec(1:1) .eq. 'T')then
            q=rec(2:5)
            i1=6

C           search list of 5-character 'official' indexed keywords
            if ( q .eq. 'BCOL' .or. q .eq. 'FORM' .or. q .eq. 'TYPE'
     &      .or. q .eq. 'UNIT' .or. q .eq. 'NULL' .or. q .eq. 'SCAL'
     &      .or. q .eq. 'ZERO' .or. q .eq. 'DISP')go to 20

C           search list of 5-character 'local' indexed keywords
            if ( q .eq. 'LMIN' .or. q .eq. 'LMAX' .or. q .eq. 'DMIN'
     &      .or. q .eq. 'DMAX' .or. q .eq. 'CTYP' .or. q .eq. 'CRPX'
     &      .or. q .eq. 'CRVL' .or. q .eq. 'CDLT' .or. q .eq. 'CROT'
     &      .or. q .eq. 'CUNI')go to 20

            q=rec(1:4)
            i1=5
C           search list of 4-character 'official' indexed keywords
            if (q .eq. 'TDIM')go to 20

C           no match so go on to next keyword
            go to 90

20          continue
C           try reading the index number suffix
            tstat=0
            call ftc2ii(rec(i1:8),ival,tstat)
            if (tstat .eq. 0 .and. ival .ge. colmin .and.
     &          ival .le. colmax)then
                if (incre .le. 0 .and. ival .eq. colmin)then
C                   delete keyword related to this column
                    call ftdrec(iunit,nrec,status)
                    nkeys=nkeys-1
                    nrec=nrec-1
                else
                    ival=ival+incre
                    i1=i1-1
                    call ftkeyn(rec(1:i1),ival,newkey,status)
                    rec(1:8)=newkey
C                   modify the index number of this keyword
                    call ftmrec(iunit,nrec,rec,status)
                end if
            end if
        end if

90      nrec=nrec+1
        if (nrec .le. nkeys)go to 100
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf