LCOV - code coverage report
Current view: top level - shared/common/src/27_toolbox_oop - m_nctk.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 74.6 % 614 458
Test Date: 2026-09-19 15:24:51 Functions: 81.8 % 44 36

            Line data    Source code
       1              : !!****m* ABINIT/m_nctk
       2              : !! NAME
       3              : !! m_nctk
       4              : !!
       5              : !! FUNCTION
       6              : !!  Tools and wrappers for NETCDF-IO.
       7              : !!
       8              : !! COPYRIGHT
       9              : !!  Copyright (C) 2008-2026 ABINIT group (MG)
      10              : !!  This file is distributed under the terms of the
      11              : !!  GNU General Public License, see ~abinit/COPYING
      12              : !!  or http://www.gnu.org/copyleft/gpl.txt .
      13              : !!
      14              : !! TODO
      15              : !!   Remove create_nc_file, write_var_netcdf, the output of OUT.nc is dangereous
      16              : !!     because we can create too many dimensions and get
      17              : !!    nf90_def_dim - NetCDF library returned:   NetCDF: NC_MAX_DIMS exceeded
      18              : !!   Moreover the multiple calls to redef render the IO very inefficient
      19              : !!   That part should be rationalized!
      20              : !!
      21              : !! SOURCE
      22              : 
      23              : #if defined HAVE_CONFIG_H
      24              : #include "config.h"
      25              : #endif
      26              : 
      27              : #include "abi_common.h"
      28              : 
      29              : MODULE m_nctk
      30              : 
      31              :  use, intrinsic :: iso_c_binding
      32              :  use defs_basis
      33              :  use m_abicore
      34              :  use m_build_info
      35              :  use m_errors
      36              :  use m_xmpi
      37              :  use netcdf
      38              : 
      39              :  use m_fstrings,  only : itoa, sjoin, lstrip, char_count, strcat, endswith, startswith, ltoa
      40              :  use m_io_tools,  only : pick_aname, delete_file, file_exists
      41              :  use m_yaml,      only : DTSET_IDX
      42              : 
      43              :  implicit none
      44              : 
      45              :  private
      46              : !!***
      47              : 
      48              :  integer,public,parameter :: nctk_noid = -huge(1)
      49              :  ! This value is used to signal to procedures that IO should not be performed.
      50              : 
      51              :  ! Basic variables
      52              :  character(len=*),public,parameter :: etsfio_file_format = "ETSF Nanoquanta"
      53              :  character(len=*),public,parameter :: etsfio_conventions = "http://www.etsf.eu/fileformats/"
      54              : 
      55              :  integer,public,parameter :: etsfio_charlen = abi_slen
      56              :  ! The value corresponding to character_string_len
      57              : 
      58              :  real,public,parameter :: etsfio_version = 3.3
      59              :  ! This is clearly wrong because one should use strings instad of floats that cannot be represented exactly
      60              :  ! but, unfortunately, it's in the specifications and we have to live with it!
      61              : 
      62              :  integer,public,parameter :: nctk_max_dims = NF90_MAX_DIMS
      63              :  ! Maximum number of dimensions
      64              : 
      65              :  integer,public,parameter :: nctk_slen = NF90_MAX_NAME
      66              :  ! String length used for the names of dimensions and variables
      67              :  ! The maximum allowable number of characters
      68              : 
      69              :  ! netcdf4-hdf5 is the default
      70              :  integer,save,private :: def_cmode_for_seq_create = ior(ior(nf90_clobber, nf90_netcdf4), nf90_write)
      71              :  ! netcdf4 classic
      72              :  !integer,save,private :: def_cmode_for_seq_create = ior(nf90_clobber, nf90_write)
      73              : 
      74              :  character(len=5),private,parameter :: NCTK_IMPLICIT_DIMS(10) = [ &
      75              :    "one  ", "two  ", "three", "four ", "five ", "six  ", "seven", "eight", "nine ", "ten  "]
      76              : 
      77              : !!****t* m_nctk/nctkerr_t
      78              : !! NAME
      79              : !! nctkerr_t
      80              : !!
      81              : !! FUNCTION
      82              : !!
      83              : !! SOURCE
      84              : !!
      85              :  type,private :: nctkerr_t
      86              :    integer :: ncerr = nf90_noerr
      87              :    integer :: line = 0
      88              :    character(len=fnlen) :: file = "Dummy File"
      89              :    character(len=2048) :: msg="No error detected"
      90              :  end type nctkerr_t
      91              : !!***
      92              : 
      93              :  type(nctkerr_t),private,save :: einfo
      94              : 
      95              : !!****t* m_nctk/ncfdim_t
      96              : !! NAME
      97              : !! nctkdim_t
      98              : !!
      99              : !! FUNCTION
     100              : !!  Stores the name and the value of a netcdf dimension
     101              : !!
     102              : !! SOURCE
     103              : 
     104              :  type,public :: nctkdim_t
     105              :    character(len=nctk_slen) :: name   ! name of the dimension.
     106              :    integer :: value                   ! value of the dimension.
     107              :    !integer :: id
     108              :  end type nctkdim_t
     109              : !!***
     110              : 
     111              : !!****t* m_nctk/nctkarr_t
     112              : !! NAME
     113              : !! nctkarr_t
     114              : !!
     115              : !! FUNCTION
     116              : !!  Stores the name and the shape of a netcdf array
     117              : !!
     118              : !! SOURCE
     119              : 
     120              :  type,public :: nctkarr_t
     121              :    character(len=nctk_slen) :: name        ! name of the array.
     122              :    character(len=4) :: dtype               ! string specifying the type.
     123              :    character(len=nctk_slen) :: shape_str   ! string with the shape. e.g. "dim1, dim2" for [dim1, dim2] array.
     124              :  end type nctkarr_t
     125              : !!***
     126              : 
     127              : !!****s* m_nctk/nctkvar_t
     128              : !! NAME
     129              : !!  nctkvar_t
     130              : !!
     131              : !! FUNCTION
     132              : !!  This structure stores variable information, such as
     133              : !!  name, NetCDF id, type, shape and dimensions. It contains the following elements:
     134              : !!
     135              : !! SOURCE
     136              : 
     137              :  type nctkvar_t
     138              : 
     139              :    integer :: id
     140              :    ! the id used by NetCDF to access this variable.
     141              : 
     142              :    integer :: xtype
     143              :    ! the type of the variable
     144              : 
     145              :    integer :: ndims
     146              :    ! the number of dimensions (0 for scalar variable).
     147              : 
     148              :    integer :: natts
     149              :    ! The number of attributes associated to the variable
     150              : 
     151              :    character(len=nctk_slen) :: name
     152              :    ! the variable name.
     153              : 
     154              :    character(len=nctk_slen) :: dimnames(nctk_max_dims)
     155              :    ! the name corresponding to each dimension
     156              :    ! Only if array variable, only (1:ndims) are relevent
     157              : 
     158              :    integer :: dimids(nctk_max_dims) = -1
     159              :    ! The id of the dimensions. only (1:ndims) are relevent
     160              : 
     161              :    integer :: dimlens(nctk_max_dims) = 0
     162              :    ! the size for each dimension if array variable, only (1:ndims) are relevent
     163              : 
     164              :    !character(len=nctk_slen) :: dimnames(nctk_max_dims)
     165              :    !character(len=nctk_slen), pointer :: ncattnames(:)
     166              :    ! * ncattnames: the name corresponding to all associated attributes
     167              : 
     168              :  end type nctkvar_t
     169              :  !!***
     170              : 
     171              :  public :: nctk_use_classic_for_seq ! Use netcdf-classic for files that are used in sequential.
     172              :                                     ! instead of the default that is netcdf4/hdf5.
     173              :  public :: nctk_idname              ! Return the nc identifier from the name of the variable.
     174              :  public :: nctk_idgroup             ! Return the nc identifier from the name of a group.
     175              :  public :: nctk_ncify               ! Append ".nc" to ipath if ipath does not end with ".nc"
     176              :  public :: nctk_string_from_occopt  ! Return human-readable string with the smearing scheme.
     177              :  public :: nctk_fort_or_ncfile      ! Test wheter a path exists (fortran or nc file) and
     178              :                                     ! select iomode depending on file extension.
     179              :  public :: nctk_try_fort_or_ncfile  ! Return fortran or netcdf filename depending on the existence of the file.
     180              :  public :: nctk_test_mpiio          ! Test at run-time whether the netcdf library supports parallel IO.
     181              : 
     182              :  public :: ab_define_var            ! Helper function used to declare a netcdf variable.
     183              : 
     184              :  ! Helper functions
     185              :  public :: nctk_open_read           ! Open a file in read-only mode.
     186              :  public :: nctk_open_create         ! Create a netcdf file for modifications.
     187              :  public :: nctk_open_modify         ! Open an already existent file for modifications.
     188              :  public :: nctk_add_etsf_header     ! Add the ETSF-IO header.
     189              :  public :: nctk_set_defmode         ! Set the file in define mode (metadata)
     190              :  public :: nctk_set_datamode        ! Set the file in datamode (IO)
     191              :  public :: nctk_set_collective      ! Set collective access for a netcdf variable
     192              : 
     193              :  public :: nctk_def_dims            ! Define dimensions in a netcdf file.
     194              :  interface nctk_def_dims
     195              :    module procedure nctk_def_one_dim
     196              :    module procedure nctk_def_dim_list
     197              :  end interface nctk_def_dims
     198              : 
     199              :  public :: nctk_set_atomic_units    ! Set the value of the attributes "units" and "scale_to_atomic_units".
     200              :  public :: nctk_def_basedims        ! Define the basic dimensions used in ETSF-IO files.
     201              :  public :: nctk_def_scalars_type    ! Declare a list of scalars of the given type.
     202              :  public :: nctk_def_iscalars        ! Declare a list of integer scalars.
     203              :  public :: nctk_def_dpscalars       ! Declare a list of double precision scalars.
     204              :  public :: nctk_write_iscalars      ! Write a list of integer scalars.
     205              :  public :: nctk_write_dpscalars     ! Write a list of double precision scalars.
     206              :  public :: nctk_defnwrite_ivars     ! Declare and write a list of integer scalars.
     207              :  public :: nctk_defnwrite_dpvars    ! Declare and write a list of double precisions scalars.
     208              :  public :: nctk_write_ibz           ! Write k-points in the IBZ and corresponding weights.
     209              : 
     210              :  public :: nctk_def_one_array
     211              :  public :: nctk_def_arrays          ! Define netcdf arrays.
     212              : 
     213              :  interface nctk_def_arrays
     214              :    module procedure nctk_def_one_array
     215              :    module procedure nctk_def_array_list
     216              :  end interface nctk_def_arrays
     217              : 
     218              :  public :: nctk_get_dim
     219              : 
     220              :  public :: nctk_write_datar
     221              :  public :: nctk_read_datar
     222              :  public :: nctk_prepare_mpiio
     223              :  ! This function appears to be required to prevent deadlocks during I/O operations in single mode.
     224              :  ! It's called automatically when using nctk_open_modify and nctk_open_read
     225              : 
     226              :  public :: create_nc_file              ! FIXME: Deprecated
     227              :  public :: write_var_netcdf            ! FIXME: Deprecated
     228              :  public :: write_eig                   ! FIXME: Deprecated
     229              : 
     230              :  !integer,save ABI_PROTECTED, public :: nctk_cache_size = 32000000
     231              :  ! If the cache_size is provided when opening a netCDF-4/HDF5 file, it will be used instead
     232              :  ! of the default (32000000) as the size, in bytes, of the HDF5 chunk cache.
     233              : 
     234              :  !integer,save ABI_PROTECTED, public :: nctk_cache_nelems = 1000
     235              :  ! If cache_nelems is provided when opening a netCDF-4/HDF5 file, it will be used instead
     236              :  ! of the default (1000) as the maximum number of elements in the HDF5 chunk cache.
     237              : 
     238              :  !real,save ABI_PROTECTED, public :: nctk_cache_premtion = 0.75
     239              :  ! If cache_preemption is provided when opening a netCDF-4/HDF5 file, it will be used
     240              :  ! instead of the default (0.75) as the preemption value for the HDF5 chunk cache.
     241              : 
     242              :  logical, save ABI_PROTECTED, public :: nctk_has_mpiio = .false.
     243              :  ! This flag is set to true if the netcdf library supports parallel IO.
     244              :  ! Cannot use CPP flags because nf90_open_par and other similar functions are always
     245              :  ! exported by netcdf. As a consequence we have to check at run-time if we can
     246              :  ! perform parallel IO and we use nctk_has_mpiio to select the IO algorithms.
     247              : 
     248              : CONTAINS
     249              : 
     250              : 
     251              : !!****f* m_nctk/nctk_set_default_for_seq
     252              : !! NAME
     253              : !!  nctk_set_default_for_seq
     254              : !!
     255              : !! FUNCTION
     256              : !!  Use netcdf classic mode for new files when only sequential-IO needs to be performed
     257              : !!
     258              : !! SOURCE
     259              : 
     260            0 : subroutine nctk_use_classic_for_seq()
     261              : 
     262              :  ! Use netcdf classic mode.
     263            0 :  def_cmode_for_seq_create = ior(nf90_clobber, nf90_write)
     264            0 :  ABI_COMMENT("Using netcdf-classic mode")
     265              : 
     266            0 : end subroutine nctk_use_classic_for_seq
     267              : !!***
     268              : 
     269              : !----------------------------------------------------------------------
     270              : 
     271              : !!****f* m_nctk/nctk_idname
     272              : !! NAME
     273              : !!  nctk_idname
     274              : !!
     275              : !! FUNCTION
     276              : !!  Return the nc identifier from the name of the variable
     277              : !!
     278              : !! SOURCE
     279              : 
     280      2964823 : integer function nctk_idname(ncid, varname) result(varid)
     281              : 
     282              : !Arguments ------------------------------------
     283              :  integer,intent(in) :: ncid
     284              :  character(len=*),intent(in) :: varname
     285              : 
     286              : !Local variables-------------------------------
     287              : !scalars
     288              :  integer :: ncerr
     289              :  character(len=1000) :: msg
     290              : ! *********************************************************************
     291              : 
     292      2964823 :  ncerr = nf90_inq_varid(ncid, varname, varid)
     293              : 
     294      2964823 :  if (ncerr /= nf90_noerr) then
     295              :    write(msg,'(6a)')&
     296            0 :      "NetCDF library returned: `",trim(nf90_strerror(ncerr)), "`", ch10,&
     297            0 :      "while trying to get the ncid of variable: ",trim(varname)
     298            0 :    ABI_ERROR(msg)
     299              :  end if
     300              : 
     301      2964823 : end function nctk_idname
     302              : !!***
     303              : 
     304              : !----------------------------------------------------------------------
     305              : 
     306              : !!****f* m_nctk/nctk_idgroup
     307              : !! NAME
     308              : !!  nctk_idgroup
     309              : !!
     310              : !! FUNCTION
     311              : !!  Return the nc identifier from the name of a group
     312              : !!
     313              : !! SOURCE
     314              : 
     315        17105 : integer function nctk_idgroup(ncid, grpname) result(grpid)
     316              : 
     317              : !Arguments ------------------------------------
     318              :  integer,intent(in) :: ncid
     319              :  character(len=*),intent(in) :: grpname
     320              : 
     321              : !Local variables-------------------------------
     322              :  integer :: ncerr
     323              :  character(len=1000) :: msg
     324              : ! *********************************************************************
     325              : 
     326        17105 :  ncerr = nf90_inq_ncid(ncid, grpname, grpid)
     327              : 
     328        17105 :  if (ncerr /= nf90_noerr) then
     329              :    write(msg,'(6a)')&
     330            0 :      "NetCDF library returned: `",trim(nf90_strerror(ncerr)), "`", ch10,&
     331            0 :      "while trying to get the ncid of group: ",trim(grpname)
     332            0 :    ABI_ERROR(msg)
     333              :  end if
     334              : 
     335        17105 : end function nctk_idgroup
     336              : !!***
     337              : 
     338              : !----------------------------------------------------------------------
     339              : !!****f* m_nctk/nctk_ncify
     340              : !! NAME
     341              : !!  nctk_ncify
     342              : !!
     343              : !! FUNCTION
     344              : !!  Append ".nc" to ipath if ipath does not end with ".nc"
     345              : !!
     346              : !! SOURCE
     347              : 
     348        56817 : function nctk_ncify(ipath) result(opath)
     349              : 
     350              :  character(len=*),intent(in) :: ipath
     351              :  character(len=fnlen) :: opath
     352              : ! *********************************************************************
     353              : 
     354        56817 :  if (.not. endswith(ipath, ".nc")) then
     355        54520 :    opath = trim(ipath)//'.nc'
     356              :  else
     357         2297 :    opath = ipath
     358              :  end if
     359              : 
     360        56817 : end function nctk_ncify
     361              : !!***
     362              : 
     363              : 
     364              : !----------------------------------------------------------------------
     365              : 
     366              : !!****f* m_nctk/nctk_string_from_occopt
     367              : !! NAME
     368              : !!  nctk_string_from_occopt
     369              : !!
     370              : !! FUNCTION
     371              : !!
     372              : !! SOURCE
     373              : 
     374        27700 : pure function nctk_string_from_occopt(occopt) result(smearing)
     375              : 
     376              :  integer,intent(in) :: occopt
     377              :  character(len=etsfio_charlen) :: smearing
     378              : ! *********************************************************************
     379              : 
     380        29191 :  select case (occopt)
     381              :  case (3)
     382         1491 :    smearing = "Fermi-Dirac"
     383              :  case (4)
     384         1802 :    smearing = "cold smearing of N. Marzari with minimization of the bump"
     385              :  case (5)
     386           42 :    smearing = "cold smearing of N. Marzari with monotonic function in the tail"
     387              :  case (6)
     388           30 :    smearing = "Methfessel and Paxton"
     389              :  case (7)
     390         3573 :    smearing = "gaussian"
     391              :  case (8)
     392            0 :    smearing = "uniform"
     393              :  case default
     394        27700 :    smearing = "none"
     395              :  end select
     396              : 
     397        27700 : end function nctk_string_from_occopt
     398              : !!***
     399              : 
     400              : !----------------------------------------------------------------------
     401              : 
     402              : !!****f* m_nctk/nctk_fort_or_ncfile
     403              : !! NAME
     404              : !!  nctk_fort_or_ncfile
     405              : !!
     406              : !! FUNCTION
     407              : !!  Return the iomode used to perform IO operations on filename.
     408              : !!  If filename does not exist, a similar file with extension `.nc` is tried
     409              : !!  and iomode is set to IO_MODE_ETSF if the file exists.
     410              : !!  This trick is used to run the Abinit test suite in netcdf mode without changing the input files.
     411              : !!  The modification (if any) is logged to std_out.
     412              : !!
     413              : !! SIDE EFFECTS
     414              : !!  filename=Tentative filename in input. Changed to netcdf file if input filename does not exist
     415              : !!   and a file with netcdf extension is found.
     416              : !!
     417              : !! OUTPUT
     418              : !!  iomode=Flag selecting the IO library. Set to IO_MODE_ETSF if netcdf file, else IO_MODE_MPI
     419              : !!    if MPI supports it, finally IO_MODE_FORTRAN
     420              : !!  errmsg=String with error message. Use `if (len_trim(errmsg) /= 0) ABI_ERROR(errmsg)`
     421              : !!    to handle possible errors in the caller.
     422              : !!
     423              : !! SOURCE
     424              : 
     425           65 : subroutine nctk_fort_or_ncfile(filename, iomode, errmsg)
     426              : 
     427              :  character(len=*),intent(inout) :: filename
     428              :  character(len=*),intent(out) :: errmsg
     429              :  integer,intent(out) :: iomode
     430              : ! *********************************************************************
     431           65 :   errmsg = ""
     432              : 
     433              :  ! Default value
     434              : #ifdef HAVE_MPI_IO
     435           65 :  iomode = IO_MODE_MPI
     436              : #else
     437              :  iomode = IO_MODE_FORTRAN
     438              : #endif
     439              : 
     440              :  !  Checking the existence of data file
     441           65 :  if (.not.file_exists(filename)) then
     442              :    ! Trick needed to run Abinit test suite in netcdf mode.
     443           50 :    if (file_exists(nctk_ncify(filename))) then
     444           50 :      write(std_out,"(3a)")"- File: ",trim(filename)," does not exist but found netcdf file with similar name."
     445          100 :      filename = nctk_ncify(filename); iomode = IO_MODE_ETSF
     446              :    end if
     447           50 :    if (.not. file_exists(filename)) then
     448            0 :      errmsg = 'Missing file: '//trim(filename)
     449              :    end if
     450              :  end if
     451              : 
     452           65 : end subroutine nctk_fort_or_ncfile
     453              : !!***
     454              : 
     455              : !----------------------------------------------------------------------
     456              : 
     457              : !!****f* m_nctk/nctk_try_fort_or_ncfile
     458              : !! NAME
     459              : !!  nctk_try_fort_or_ncfile
     460              : !!
     461              : !! FUNCTION
     462              : !!  If filename does not exist, a similar file with extension `.nc` is tried
     463              : !!  This trick is used to run the Abinit test suite in netcdf mode without changing the input files.
     464              : !!  The modification (if any) is logged to unit (Default: std_out)
     465              : !!
     466              : !! SIDE EFFECTS
     467              : !!  filename=Tentative filename in input. Changed to netcdf file if input filename does not exist
     468              : !!   and a file with netcdf extension is found.
     469              : !!
     470              : !! OUTPUT
     471              : !!  errmsg=String with error message if return value /= 0
     472              : !!
     473              : !! SOURCE
     474              : 
     475        11823 : integer function nctk_try_fort_or_ncfile(filename, errmsg, unit) result(ierr)
     476              : 
     477              : !Arguments ------------------------------------
     478              :  character(len=*),intent(inout) :: filename
     479              :  character(len=*),intent(out) :: errmsg
     480              :  integer,optional,intent(in) :: unit
     481              : 
     482              : !Local variables-------------------------------
     483              :  integer :: unt
     484              : ! *********************************************************************
     485              : 
     486            0 :  unt = std_out; if (present(unit)) unt = unit
     487        11823 :  ierr = 0; errmsg = ""
     488              : 
     489        11823 :  if (.not.file_exists(filename)) then
     490              :    ! Try netcdf exists.
     491         2510 :    if (file_exists(nctk_ncify(filename))) then
     492         2510 :      if (unt /= dev_null) then
     493         2510 :        write(unt,"(3a)")"- File: ",trim(filename)," does not exist but found netcdf file with similar name."
     494              :      end if
     495         5020 :      filename = nctk_ncify(filename)
     496              :    end if
     497         2510 :    if (.not. file_exists(filename)) then
     498            0 :      ierr = 1; errmsg = 'Cannot find file: '//trim(filename)
     499              :    end if
     500              :  end if
     501              : 
     502        11823 : end function nctk_try_fort_or_ncfile
     503              : !!***
     504              : 
     505              : !----------------------------------------------------------------------
     506              : 
     507              : !!****f* m_nctk/nctk_test_mpiio
     508              : !! NAME
     509              : !!  nctk_test_mpiio
     510              : !!
     511              : !! FUNCTION
     512              : !!  Test at run-time whether the netcdf library supports parallel IO and
     513              : !!  set the value of the module variable `nctk_has_mpiio`.
     514              : !!  This is a COLLECTIVE routine that should be called by all processors
     515              : !!  in MPI_COMM_WORLD at the beginning of the calculation
     516              : !!
     517              : !! INPUTS
     518              : !!  [print_warning]=TRUE if a warning about paral_kgb use has to be printed
     519              : !!                  Optional, default=yes
     520              : !!
     521              : !! SOURCE
     522              : 
     523         1447 : subroutine nctk_test_mpiio(print_warning)
     524              : 
     525              :  logical,intent(in),optional :: print_warning
     526              : 
     527              : !Local variables-------------------------------
     528              : !scalars
     529              :  logical :: my_print_warning
     530              : #ifdef HAVE_NETCDF_MPI
     531              :  integer,parameter :: master=0
     532              :  integer :: ierr,ncid,ncerr
     533              :  character(len=500) :: msg
     534              :  character(len=fnlen) :: apath
     535              : #endif
     536              : ! *********************************************************************
     537              : 
     538         1447 :  nctk_has_mpiio = .False.
     539         1447 :  my_print_warning=.true. ; if (present(print_warning)) my_print_warning=print_warning
     540              : 
     541              : !FIXME nf90create fails when using NVHPC
     542              : ! This might be due to my environment, maybe not, need to investigate this...
     543              : !!#ifndef FC_NVHPC
     544              : #ifdef HAVE_NETCDF_MPI
     545         1447 :  if (xmpi_comm_rank(xmpi_world) == master) then
     546              :    ! Try to open a file with hdf5.
     547         1169 :    apath = pick_aname()
     548              :    ncerr = nf90_create(apath, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), ncid=ncid, &
     549         1169 :      comm=xmpi_comm_self, info=xmpio_info)
     550              : 
     551         1169 :    if (ncerr == nf90_noerr) then
     552         1169 :      nctk_has_mpiio = .True.
     553         1169 :      call wrtout(std_out," Netcdf library supports MPI-IO", "COLL")
     554            0 :    else if (ncerr == nf90_enopar) then
     555              :      ! This is the value returned by the C function ifndef USE_PARALLEL
     556            0 :      ABI_WARNING(sjoin("Netcdf lib does not support MPI-IO and: ", nf90_strerror(ncerr)))
     557            0 :      nctk_has_mpiio = .False.
     558              :    else
     559              :      ! Maybe something wrong in the low-level layer!
     560            0 :      ABI_WARNING(sjoin("Strange, netcdf seems to support MPI-IO but: ", nf90_strerror(ncerr)))
     561            0 :      nctk_has_mpiio = .False.
     562              :    end if
     563              : 
     564         1169 :    ncerr = nf90_close(ncid)
     565         1169 :    call delete_file(apath, ierr)
     566              :  end if
     567              : 
     568              :  ! Master broadcast nctk_has_mpiio
     569         1447 :  call xmpi_bcast(nctk_has_mpiio,master,xmpi_world,ierr)
     570              : 
     571         1447 :  if ((.not. nctk_has_mpiio) .and. my_print_warning) then
     572              :    write(msg,"(5a)") &
     573            0 :       "The netcdf library does not support parallel IO, see message above",ch10,&
     574            0 :       "Abinit won't be able to produce files in parallel when e.g. paral_kgb == 1 is used.",ch10,&
     575            0 :       "Action: install a netcdf4+HDF5 library with MPI-IO support."
     576            0 :    ABI_WARNING(msg)
     577              :  end if
     578              : #endif
     579              : !!#endif
     580              : 
     581              : #ifdef HAVE_NETCDF_DEFAULT
     582         1447 :  if (.not. nctk_has_mpiio) then
     583            0 :    ABI_ERROR("--netcdf-default is on but netcdf library does not support MPI-IO. Aborting now")
     584              :  end if
     585              : #endif
     586              : 
     587         1447 : end subroutine nctk_test_mpiio
     588              : !!***
     589              : 
     590              : !!****f* m_nctk/str2xtype
     591              : !! NAME
     592              : !!  str2xtype
     593              : !!
     594              : !! FUNCTION
     595              : !!  Return the netcdf type from a string. Possible values:
     596              : !!    c or ch   for NF90_CHAR
     597              : !!    i or int  for NF90_INTtrue
     598              : !!   sp         for NF90_FLOAT
     599              : !!   dp         for NF90_DOUBLE
     600              : !!
     601              : !! SOURCE
     602              : 
     603      1200687 : integer function str2xtype(string) result(xtype)
     604              : 
     605              : !Arguments ------------------------------------
     606              :  character(len=*),intent(in) :: string
     607              : ! *********************************************************************
     608              : 
     609              :  !Type  FORTRAN API Mnemonic    Bits
     610              :  !byte      NF90_BYTE           8
     611              :  !char      NF90_CHAR           8
     612              :  !short     NF90_SHORT          16
     613              :  !int       NF90_INT            32
     614              :  !float     NF90_FLOAT          32
     615              :  !double    NF90_DOUBLE         64
     616              : 
     617              :  select case (string)
     618              :  case ("c", "ch", "char")
     619              :    xtype = NF90_CHAR
     620              :  case ("i", "int")
     621              :    xtype = NF90_INT
     622              :  case ("sp")
     623              :    xtype = NF90_FLOAT
     624              :  case ("dp")
     625            0 :    xtype = NF90_DOUBLE
     626              :  case default
     627      1200687 :    ABI_ERROR(sjoin("Invalid string type:", string))
     628              :  end select
     629              : 
     630      1200687 : end function str2xtype
     631              : !!***
     632              : 
     633              : !!****f* m_nctk/bail_if_ncerr
     634              : !! NAME
     635              : !!  bail_if_ncerr
     636              : !!
     637              : !! FUNCTION
     638              : !!
     639              : !! INPUTS
     640              : !!  ncerr=Netcdf error
     641              : !!  line=line number of the file where problem occurred
     642              : !!  file=name of the f90 file containing the caller
     643              : !!
     644              : !! SOURCE
     645              : 
     646              : logical function bail_if_ncerr(ncerr, file, line) result(bail)
     647              : 
     648              : !Arguments ------------------------------------
     649              :  integer,intent(in) :: ncerr
     650              :  character(len=*),optional,intent(in) :: file
     651              :  integer,optional,intent(in) :: line
     652              : ! *********************************************************************
     653              : 
     654              :  bail = (ncerr /= nf90_noerr)
     655              : 
     656              :  if (bail) then
     657              :    einfo%ncerr = ncerr
     658              :    einfo%file = "Subroutine Unknown"; if (present(file)) einfo%file = file
     659              :    einfo%line = 0; if (present(line)) einfo%line = line
     660              :    ! Append Netcdf string to user-defined message.
     661              :    write(einfo%msg,'(2a)')'NetCDF library raised: ',trim(nf90_strerror(ncerr))
     662              :  end if
     663              : 
     664              : end function bail_if_ncerr
     665              : !!***
     666              : 
     667              : !----------------------------------------------------------------------
     668              : 
     669              : !!****f* m_nctk/nctk_open_read
     670              : !! NAME
     671              : !!  nctk_open_read
     672              : !!
     673              : !! FUNCTION
     674              : !!   Open a netcdf file in read-only mode. Return exit status.
     675              : !!
     676              : !! INPUTS
     677              : !!  ncid=Netcdf identifier.
     678              : !!  comm=MPI communicator.
     679              : !!
     680              : !! SOURCE
     681              : 
     682        22031 : integer function nctk_open_read(ncid, path, comm) result(ncerr)
     683              : 
     684              : !Arguments ------------------------------------
     685              :  integer,intent(out) :: ncid
     686              :  integer,intent(in) :: comm
     687              :  character(len=*),intent(in) :: path
     688              : 
     689              : !Local variables-------------------------------
     690              :  integer :: nprocs
     691              : ! *********************************************************************
     692              : 
     693        22031 :  nprocs = xmpi_comm_size(comm)
     694              : 
     695              :  ! Enforce netcdf4 only if the communicator contains more than one processor.
     696        22031 :  if (nctk_has_mpiio .and. nprocs > 1) then
     697              : #ifdef HAVE_NETCDF_MPI
     698              :    ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_nowrite),&
     699           48 :                      comm=comm, info=xmpio_info, ncid=ncid)
     700              : #else
     701              :    ncerr = nf90_einval
     702              :    ABI_WARNING("Netcdf without MPI support. Cannot open file, will abort in caller")
     703              : #endif
     704           48 :    NCF_CHECK_MSG(ncerr, sjoin("opening file:", path))
     705              :  else
     706        21983 :    ncerr = nf90_open(path, mode=nf90_nowrite, ncid=ncid)
     707        21983 :    NCF_CHECK_MSG(ncerr, sjoin("Opening file:", path))
     708              :    !if (ncerr /= NC_EHDFERR) then
     709              :    !  ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4), nf90_nowrite), comm=comm, info=xmpio_info, ncid=ncid)
     710              :    !end if
     711        21983 :    if (nprocs > 1) then
     712            0 :      ncerr = nf90_einval
     713            0 :      ABI_WARNING("netcdf without MPI-IO support with nprocs > 1! Will abort in the caller")
     714              :    end if
     715              :  end if
     716              : 
     717        22031 : end function nctk_open_read
     718              : !!***
     719              : 
     720              : !----------------------------------------------------------------------
     721              : 
     722              : !!****f* m_nctk/nctk_open_create
     723              : !! NAME
     724              : !!  nctk_open_create
     725              : !!
     726              : !! FUNCTION
     727              : !!  Create and open the netcdf file. Return exis status.
     728              : !!
     729              : !! INPUTS
     730              : !!  path=Name of the file
     731              : !!  comm=MPI communicator.
     732              : !!
     733              : !! OUTPUT
     734              : !!  ncid=Netcdf identifier.
     735              : !!
     736              : !! SOURCE
     737              : 
     738        15885 : integer function nctk_open_create(ncid, path, comm) result(ncerr)
     739              : 
     740              : !Arguments ------------------------------------
     741              :  integer,intent(out) :: ncid
     742              :  integer,intent(in) :: comm
     743              :  character(len=*),intent(in) :: path
     744              : 
     745              : !Local variables-------------------------------
     746              :  integer :: input_len, cmode
     747              :  character(len=strlen) :: my_string
     748              : ! *********************************************************************
     749              : 
     750              :  ! Always use mpiio mode (i.e. hdf5) if available so that one can perform parallel IO
     751        15885 :  if (nctk_has_mpiio) then
     752        15649 :    ncerr = nf90_einval
     753              : #ifdef HAVE_NETCDF_MPI
     754        15649 :    write(my_string,'(2a)') "- Creating HDf5 file with MPI-IO support: ",trim(path)
     755        15649 :    call wrtout(std_out, my_string)
     756              :    ! Believe it or not, I have to use xmpi_comm_self even in sequential to avoid weird SIGSEV in the MPI layer!
     757        15649 :    ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), ncid=ncid, comm=comm, info=xmpio_info)
     758              : #endif
     759              :  else
     760              :    ! Note that here we don't enforce nf90_netcdf4 hence the netcdf file with be in classic model.
     761          236 :    write(my_string,'(2a)') "- Creating HDf5 file WITHOUT MPI-IO support: ",trim(path)
     762          236 :    call wrtout(std_out, my_string)
     763              :    !ncerr = nf90_create(path, ior(nf90_clobber, nf90_write), ncid)
     764          236 :    cmode = def_cmode_for_seq_create
     765          236 :    ncerr = nf90_create(path, cmode=cmode, ncid=ncid)
     766          236 :    if (xmpi_comm_size(comm) > 1) then
     767            0 :      ABI_WARNING("netcdf without MPI-IO support with nprocs > 1!")
     768              :    end if
     769              :  end if
     770        15885 :  NCF_CHECK(ncerr)
     771              : 
     772              :  ! Write etsf_header: file format, version and conventions.
     773        15885 :  NCF_CHECK(nf90_put_att(ncid, NF90_GLOBAL, "file_format", etsfio_file_format))
     774        15885 :  NCF_CHECK(nf90_put_att(ncid, NF90_GLOBAL, "file_format_version", etsfio_version))
     775        15885 :  NCF_CHECK(nf90_put_att(ncid, NF90_GLOBAL, "Conventions", etsfio_conventions))
     776              : 
     777              :  ! Add info on the code that produced this file. These are extensions not in the standard.
     778        15885 :  NCF_CHECK(nf90_put_att(ncid, NF90_GLOBAL, "code", "Abinit"))
     779        15885 :  NCF_CHECK(nf90_put_att(ncid, NF90_GLOBAL, "abinit_version", abinit_version))
     780              : 
     781              :  ! Define the basic dimensions used in ETSF-IO files.
     782        15885 :  NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
     783              : 
     784              :  ! INPUT_STRING is allocated and initialized in parsefile
     785        15885 :  if (allocated(INPUT_STRING)) then
     786              :    ! Write string with input.
     787        15808 :    my_string = trim(INPUT_STRING)
     788        15808 :    if (DTSET_IDX /= -1 .and. index(INPUT_STRING, "jdtset ") == 0) then
     789        14021 :      my_string = "jdtset " // trim(itoa(DTSET_IDX)) // "  " // trim(INPUT_STRING)
     790              :    end if
     791              : 
     792        15808 :    input_len = len_trim(my_string)
     793        15808 :    NCF_CHECK(nctk_def_dims(ncid, nctkdim_t("input_len", input_len)))
     794        15808 :    NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("input_string", "c", "input_len")))
     795              :    !print *, "input_len, strlen:", input_len, strlen
     796              : 
     797        15808 :    if (xmpi_comm_rank(comm) == 0) then
     798        15808 :      NCF_CHECK(nctk_set_datamode(ncid))
     799              :      ! Pass my_string(1:input_len)) instead from trim(string) to avoid SIGSEV on higgs_intel_19.0_serial
     800        15808 :      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "input_string"), my_string(1:input_len)))
     801        15808 :      NCF_CHECK(nctk_set_defmode(ncid))
     802              :    end if
     803              :  end if
     804              : 
     805        15885 : end function nctk_open_create
     806              : !!***
     807              : 
     808              : !----------------------------------------------------------------------
     809              : 
     810              : !!****f* m_nctk/nctk_open_modify
     811              : !! NAME
     812              : !!  nctk_open_modfy
     813              : !!
     814              : !! FUNCTION
     815              : !!   Open an already existent netcdf file for modifications. Return exit status.
     816              : !!
     817              : !! INPUTS
     818              : !!  path=File name.
     819              : !!  comm=MPI communicator.
     820              : !!
     821              : !! OUTPUT
     822              : !!  ncid=Netcdf identifier.
     823              : !!
     824              : !! SOURCE
     825              : 
     826        11859 : integer function nctk_open_modify(ncid, path, comm) result(ncerr)
     827              : 
     828              : !Arguments ------------------------------------
     829              :  integer,intent(out) :: ncid
     830              :  character(len=*),intent(in) :: path
     831              :  integer,intent(in) :: comm
     832              : ! *********************************************************************
     833              : 
     834        11859 :  if (.not. nctk_has_mpiio .and. xmpi_comm_size(comm) > 1) then
     835            0 :    ABI_ERROR("netcdf without MPI-IO support and nprocs > 1!")
     836              :  end if
     837              : 
     838        11859 :  if (xmpi_comm_size(comm) > 1 .or. nctk_has_mpiio) then
     839        11858 :    call wrtout(std_out, sjoin(" nctk_open_modify: Opening HDf5 file with MPI-IO support:", path))
     840              : #ifdef HAVE_NETCDF_MPI
     841              :    ncerr = nf90_open_par(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
     842        11858 :                          comm=comm, info=xmpio_info, ncid=ncid)
     843        11858 :    NCF_CHECK_MSG(ncerr, sjoin("nf90_open_par: ", path))
     844              : #else
     845              :    ABI_ERROR("nprocs > 1 but netcdf does not support MPI-IO")
     846              : #endif
     847              :  else
     848            1 :    call wrtout(std_out, sjoin(" nctk_open_modify: Opening netcdf file without MPI-IO support:", path))
     849            1 :    ncerr = nf90_open(path, nf90_write, ncid)
     850            1 :    NCF_CHECK_MSG(ncerr, sjoin("nf90_open: ", path))
     851              :  end if
     852              : 
     853              :  ! Set file in define mode.
     854        11859 :  NCF_CHECK(nctk_set_defmode(ncid))
     855              :  !call wrtout(std_out, "- Returning from nctk_open_modify")
     856              : 
     857        11859 : end function nctk_open_modify
     858              : !!***
     859              : 
     860              : !----------------------------------------------------------------------
     861              : 
     862              : !!****f* m_nctk/nctk_add_etsf_header
     863              : !! NAME
     864              : !!  nctk_add_etsf_header
     865              : !!
     866              : !! FUNCTION
     867              : !!   Add the etsf-io header to a file associated to a netcdf file handler.
     868              : !!
     869              : !! INPUTS
     870              : !!  ncid=Netcdf file identifier.
     871              : !!  * version = the number of version to be created.
     872              : !!  * title = (optional) a title for the file (80 characters max).
     873              : !!  * history = (optional) the first line of history (1024 characters max).
     874              : !!  * with_etsf_header = (optional) if true, will create a header
     875              : !!                       as defined in the ETSF specifications (default is .true.).
     876              : !!                       When value is .false., arguments title, history and version
     877              : !!                       are ignored.
     878              : !!
     879              : !! SOURCE
     880              : 
     881         4561 : integer function nctk_add_etsf_header(ncid, title, history) result(ncerr)
     882              : 
     883              : !Arguments ------------------------------------
     884              :  integer,intent(in) :: ncid
     885              :  character(len=*),optional,intent(in) :: title,history
     886              : 
     887              : !Local variables-------------------------------
     888              :  !integer :: ncerr
     889              :  character(len=*), parameter :: file_format = "ETSF Nanoquanta"
     890              :  character(len=*), parameter :: conventions = "http://www.etsf.eu/fileformats/"
     891              :  real :: format_version = 3.3 ! Real is not a good choice for a version!
     892              : ! *********************************************************************
     893              : 
     894         4561 :  ncerr = nctk_set_defmode(ncid)
     895         4561 :  if (ncerr /= nf90_noerr) return
     896              : 
     897              :  ! The file format
     898         4561 :  ncerr = nf90_put_att(ncid, NF90_GLOBAL, "file_format", file_format)
     899         4561 :  if (ncerr /= nf90_noerr) return
     900              : 
     901              :  ! The version
     902         4561 :  ncerr = nf90_put_att(ncid, NF90_GLOBAL, "file_format_version", format_version)
     903         4561 :  if (ncerr /= nf90_noerr) return
     904              : 
     905              :  ! The conventions
     906         4561 :  ncerr = nf90_put_att(ncid, NF90_GLOBAL, "Conventions", conventions)
     907         4561 :  if (ncerr /= nf90_noerr) return
     908              : 
     909              :  ! The history if present
     910         4561 :  if (present(history)) then
     911         4561 :    ncerr = nf90_put_att(ncid, NF90_GLOBAL, "history", history(1:min(1024, len(history))))
     912         4561 :    if (ncerr /= nf90_noerr) return
     913              :  end if
     914              : 
     915              :  ! The title if present
     916         4561 :  if (present(title)) then
     917         4561 :    ncerr = nf90_put_att(ncid, NF90_GLOBAL, "title", title(1:min(80, len(title))))
     918         4561 :    if (ncerr /= nf90_noerr) return
     919              :  end if
     920              : 
     921              :  ! These are extensions not in the standard.
     922              :  ! Add info on the code that produced this file
     923         4561 :  ncerr = nf90_put_att(ncid, NF90_GLOBAL, "code", "Abinit")
     924         4561 :  if (ncerr /= nf90_noerr) return
     925              : 
     926         4561 :  ncerr = nf90_put_att(ncid, NF90_GLOBAL, "code_version", ABINIT_VERSION)
     927              :  if (ncerr /= nf90_noerr) return
     928              : 
     929              : end function nctk_add_etsf_header
     930              : !!***
     931              : 
     932              : !----------------------------------------------------------------------
     933              : 
     934              : !!****f* m_nctk/nctk_set_defmode
     935              : !! NAME
     936              : !!  nctk_set_defmode
     937              : !!
     938              : !! FUNCTION
     939              : !!   Set the file in define mode, return exit status.
     940              : !!
     941              : !! INPUTS
     942              : !!  ncid=Netcdf identifier.
     943              : !!
     944              : !! SOURCE
     945              : 
     946       214743 : integer function nctk_set_defmode(ncid) result(ncerr)
     947              : 
     948              : !Arguments ------------------------------------
     949              :  integer,intent(in) :: ncid
     950              : ! *********************************************************************
     951              : 
     952       214743 :  ncerr = nf90_redef(ncid)
     953              :  ! Use same trick as in etsf_io
     954       214743 :  if (ncerr /= nf90_noerr .and. ncerr /= -39) then
     955            0 :    NCF_CHECK(ncerr)
     956              :  else
     957       214743 :    ncerr = nf90_noerr
     958              :  end if
     959              : 
     960       214743 : end function nctk_set_defmode
     961              : !!***
     962              : 
     963              : !----------------------------------------------------------------------
     964              : 
     965              : !!****f* m_nctk/nctk_set_datamode
     966              : !! NAME
     967              : !!  nctk_set_datamode
     968              : !!
     969              : !! FUNCTION
     970              : !!  Set the file in data mode. Return exit status
     971              : !!
     972              : !! INPUTS
     973              : !!  ncid=Netcdf identifier.
     974              : !!  [reserve]
     975              : !!
     976              : !! OUTPUT
     977              : !!  ncerr=Exit status
     978              : !!
     979              : !! SOURCE
     980              : 
     981       141658 : integer function nctk_set_datamode(ncid, reserve) result(ncerr)
     982              : 
     983              : !Arguments ------------------------------------
     984              :  integer,intent(in) :: ncid
     985              :  logical,intent(in),optional :: reserve
     986              : 
     987              : !Local variables-------------------------------
     988              : !scalars
     989              :  logical :: do_reserve
     990              : ! *********************************************************************
     991              : 
     992       141658 :  do_reserve = .False.; if (present(reserve)) do_reserve = reserve
     993              : 
     994       141658 :  ncerr = nf90_enddef(ncid)
     995              : 
     996              :  ! Use same trick as in etsf_io
     997              :  ! needed otherwise netcdf complains if the file is already in def mode.
     998       141658 :  if (ncerr /= nf90_noerr .and. ncerr /= -38) then
     999            0 :    NCF_CHECK(ncerr)
    1000              :  else
    1001       141658 :    ncerr = nf90_noerr
    1002              :  end if
    1003              : 
    1004       141658 :  return
    1005              : 
    1006              :  ! TODO
    1007              :  if (do_reserve) then
    1008              :    ncerr = nf90_enddef(ncid)
    1009              :    !ncerr = nf90__enddef(ncid)
    1010              :  else
    1011              :    ncerr = nf90_enddef(ncid)
    1012              :  end if
    1013              : 
    1014              : end function nctk_set_datamode
    1015              : !!***
    1016              : 
    1017              : !----------------------------------------------------------------------
    1018              : 
    1019              : !!****f* m_nctk/nctk_set_collective
    1020              : !! NAME
    1021              : !!  nctk_set_collective
    1022              : !!
    1023              : !! FUNCTION
    1024              : !!  Use collective IO for the given netcdf variable. Return exit status.
    1025              : !!
    1026              : !! INPUTS
    1027              : !!  ncid=Netcdf file identifier.
    1028              : !!  varid=Netcdf variable identifier.
    1029              : !!  [independent]=True of use indepedent mode.
    1030              : !!
    1031              : !! SOURCE
    1032              : 
    1033       102560 : integer function nctk_set_collective(ncid, varid, independent) result(ncerr)
    1034              : 
    1035              : !Arguments ------------------------------------
    1036              :  integer,intent(in) :: ncid, varid
    1037              :  logical,optional,intent(in) :: independent
    1038              : 
    1039              : !Local variables-------------------------------
    1040              :  logical :: independent__
    1041              : ! *********************************************************************
    1042              : 
    1043       102560 :   ncerr = nf90_einval
    1044       102302 :   independent__ = .false.
    1045              : #ifdef HAVE_NETCDF_MPI
    1046       102302 :   if (present(independent)) independent__ = independent
    1047            0 :   if (independent__) then
    1048            0 :     ncerr = nf90_var_par_access(ncid, varid, nf90_independent)
    1049              :   else
    1050       102302 :     ncerr = nf90_var_par_access(ncid, varid, nf90_collective)
    1051              :   end if
    1052              : #else
    1053              :   ABI_ERROR("nctk_set_collective should not be called if NETCDF does not support MPI-IO")
    1054              :   ABI_UNUSED((/ncid, varid/))
    1055              : #endif
    1056              : 
    1057       102302 : end function nctk_set_collective
    1058              : !!***
    1059              : 
    1060              : !!****f* m_nctk/nctk_def_one_dim
    1061              : !! NAME
    1062              : !!  nctk_def_one_dim
    1063              : !!
    1064              : !! FUNCTION
    1065              : !!  Define list of dimensions variables and write their values.
    1066              : !!  Return immediately if error
    1067              : !!
    1068              : !! INPUTS
    1069              : !!  ncid=Netcdf identifier.
    1070              : !!  dimnames(:)=List of strings with the name of the dimensions
    1071              : !!  values(:)=List of integer scalars
    1072              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1073              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1074              : !!
    1075              : !! SOURCE
    1076              : 
    1077      1343108 : integer function nctk_def_one_dim(ncid, nctkdim, defmode, prefix) result(ncerr)
    1078              : 
    1079              : !Arguments ------------------------------------
    1080              : !scalars
    1081              :  integer,intent(in) :: ncid
    1082              :  logical,optional,intent(in) :: defmode
    1083              :  character(len=*),optional,intent(in) :: prefix
    1084              : !arrays
    1085              :  type(nctkdim_t),intent(in) :: nctkdim
    1086              : 
    1087              : !Local variables-------------------------------
    1088              :  integer :: dimid,dimlen
    1089              :  character(len=nctk_slen) :: dname
    1090              :  character(len=500) :: msg
    1091              : ! *********************************************************************
    1092              : 
    1093      1343108 :  ncerr = nf90_noerr
    1094              : 
    1095      1343108 :  if (present(defmode)) then
    1096            0 :    if (defmode) then
    1097            0 :      NCF_CHECK(nctk_set_defmode(ncid))
    1098              :    end if
    1099              :  end if
    1100              : 
    1101      1343108 :  if (present(prefix)) then
    1102          396 :    if (any(nctkdim%name == NCTK_IMPLICIT_DIMS)) then
    1103            0 :      dname = nctkdim%name
    1104              :    else
    1105           36 :      dname = strcat(prefix, nctkdim%name)
    1106              :    end if
    1107              :  else
    1108      1343072 :    dname = nctkdim%name
    1109              :  end if
    1110              : 
    1111              :  ! if dimension already exists, test whether it has the same value else define new dim.
    1112      1343108 :  ncerr = nf90_inq_dimid(ncid, dname, dimid)
    1113              : 
    1114      1343108 :  if (ncerr == nf90_noerr) then
    1115       511852 :    NCF_CHECK(nf90_inquire_dimension(ncid, dimid, len=dimlen))
    1116       511852 :    if (dimlen /= nctkdim%value) then
    1117              :      write(msg, "(4a,2(a,i0))")&
    1118            0 :         "dimension ", trim(dname)," already exists but with a different value",ch10,&
    1119            0 :         "from file: ", dimlen, "; about to write: ", nctkdim%value
    1120            0 :      ABI_ERROR(msg)
    1121              :    end if
    1122              :  else
    1123       831256 :    ncerr = nf90_def_dim(ncid, dname, nctkdim%value, dimid)
    1124       831256 :    NCF_CHECK(ncerr)
    1125              :  end if
    1126              : 
    1127      1343108 : end function nctk_def_one_dim
    1128              : !!***
    1129              : 
    1130              : !----------------------------------------------------------------------
    1131              : 
    1132              : !!****f* m_nctk/nctk_def_dim_list
    1133              : !! NAME
    1134              : !!  nctk_def_dim_list
    1135              : !!
    1136              : !! FUNCTION
    1137              : !!  Define list of dimensions variables and write their values.
    1138              : !!  Return immediately if error
    1139              : !!
    1140              : !! INPUTS
    1141              : !!  ncid=Netcdf identifier.
    1142              : !!  dimnames(:)=List of strings with the name of the dimensions
    1143              : !!  values(:)=List of integer scalars
    1144              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1145              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1146              : !!
    1147              : !! SOURCE
    1148              : 
    1149       219257 : integer function nctk_def_dim_list(ncid, nctkdims, defmode, prefix) result(ncerr)
    1150              : 
    1151              : !Arguments ------------------------------------
    1152              : !scalars
    1153              :  integer,intent(in) :: ncid
    1154              :  logical,optional,intent(in) :: defmode
    1155              :  character(len=*),optional,intent(in) :: prefix
    1156              : !arrays
    1157              :  type(nctkdim_t),intent(in) :: nctkdims(:)
    1158              : 
    1159              : !Local variables-------------------------------
    1160              :  integer :: ii
    1161              : ! *********************************************************************
    1162              : 
    1163       219257 :  ncerr = nf90_noerr
    1164       219257 :  if (present(defmode)) then
    1165       105791 :    if (defmode) then
    1166       105791 :      NCF_CHECK(nctk_set_defmode(ncid))
    1167              :    end if
    1168              :  end if
    1169              : 
    1170      1538000 :  do ii=1,size(nctkdims)
    1171      1318743 :    if (present(prefix)) then
    1172           36 :      ncerr = nctk_def_one_dim(ncid, nctkdims(ii), prefix=prefix)
    1173              :    else
    1174      1318707 :      ncerr = nctk_def_one_dim(ncid, nctkdims(ii))
    1175              :    end if
    1176      1538000 :    if (ncerr /= nf90_noerr) return
    1177              :  end do
    1178              : 
    1179       219257 : end function nctk_def_dim_list
    1180              : !!***
    1181              : 
    1182              : !----------------------------------------------------------------------
    1183              : 
    1184              : !!****f* m_nctk/nctk_set_atomic_units
    1185              : !! NAME
    1186              : !!  nctk_set_atomic_units
    1187              : !!
    1188              : !! FUNCTION
    1189              : !!  Set the attributes "units" to "atomic units" and "scale_to_atomic_units" to one.
    1190              : !!
    1191              : !! INPUTS
    1192              : !!  ncid=Netcdf identifier.
    1193              : !!  varname=Name of the variable
    1194              : !!
    1195              : !! SOURCE
    1196              : 
    1197        97604 : integer function nctk_set_atomic_units(ncid, varname) result(ncerr)
    1198              : 
    1199              : !Arguments ------------------------------------
    1200              :  integer,intent(in) :: ncid
    1201              :  character(len=*),intent(in) :: varname
    1202              : 
    1203              : !Local variables-------------------------------
    1204              :  integer :: varid
    1205              : ! *********************************************************************
    1206              : 
    1207        97604 :  ncerr = nf90_noerr
    1208              : 
    1209        97604 :  varid = nctk_idname(ncid, varname)
    1210        97604 :  NCF_CHECK(nf90_put_att(ncid, varid, "units", "atomic units"))
    1211        97604 :  NCF_CHECK(nf90_put_att(ncid, varid, "scale_to_atomic_units", one))
    1212              : 
    1213        97604 : end function nctk_set_atomic_units
    1214              : !!***
    1215              : 
    1216              : !----------------------------------------------------------------------
    1217              : 
    1218              : !!****f* m_nctk/nctk_def_basedims
    1219              : !! NAME
    1220              : !!  nctk_def_basedims
    1221              : !!
    1222              : !! FUNCTION
    1223              : !!  Define the basic dimensions used in ETSF-IO files.
    1224              : !!
    1225              : !! INPUTS
    1226              : !!  ncid=Netcdf identifier.
    1227              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1228              : !!
    1229              : !! SOURCE
    1230              : 
    1231        34526 : integer function nctk_def_basedims(ncid, defmode) result(ncerr)
    1232              : 
    1233              : !Arguments ------------------------------------
    1234              :  integer,intent(in) :: ncid
    1235              :  logical,optional,intent(in) :: defmode
    1236              : ! *********************************************************************
    1237              : 
    1238        34526 :  ncerr = nf90_noerr
    1239              : 
    1240        34526 :  if (present(defmode)) then
    1241        32593 :    if (defmode) then
    1242        32593 :      NCF_CHECK(nctk_set_defmode(ncid))
    1243              :    end if
    1244              :  end if
    1245              : 
    1246              :  ! Basic ETSF-IO dimensions that should be always present in the file.
    1247              :  ncerr = nctk_def_dims(ncid, [&
    1248              :    nctkdim_t("complex", 2), nctkdim_t("symbol_length", 2), nctkdim_t("character_string_length", etsfio_charlen),&
    1249              :    nctkdim_t("number_of_cartesian_directions", 3), nctkdim_t("number_of_reduced_dimensions", 3),&
    1250              :    nctkdim_t("number_of_vectors", 3) &
    1251       241682 :  ])
    1252        34526 :  NCF_CHECK(ncerr)
    1253              : 
    1254              :  ! Useful integers.
    1255              :  ncerr = nctk_def_dims(ncid, [ &
    1256              :    nctkdim_t("one", 1), nctkdim_t("two", 2), nctkdim_t("three", 3), &
    1257              :    nctkdim_t("four", 4), nctkdim_t("five", 5), nctkdim_t("six", 6), &
    1258              :    nctkdim_t("seven", 7), nctkdim_t("eight", 8), nctkdim_t("nine", 9), nctkdim_t("ten", 10), &
    1259              :    nctkdim_t("fnlen", fnlen + 1) &
    1260       414312 :  ])
    1261        34526 :  NCF_CHECK(ncerr)
    1262              : 
    1263        34526 : end function nctk_def_basedims
    1264              : !!***
    1265              : 
    1266              : !!****f* m_nctk/ab_define_var
    1267              : !!
    1268              : !! NAME
    1269              : !! ab_define_var
    1270              : !!
    1271              : !! FUNCTION
    1272              : !! Write the definition of a variable, including units and mnemonics
    1273              : !!
    1274              : !! INPUTS
    1275              : !! ncid = Identifier of the netcdf dataset
    1276              : !! var_dim_id = Identifier of the Dimensions
    1277              : !! var_id     = Identifier of the variable
    1278              : !! var_mnemo  = String of mnemonics
    1279              : !! var_name   = String with the name of the variable
    1280              : !! var_type   = NetCDF type of variable (NF90_DOUBLE, etc)
    1281              : !! var_units  = String of units
    1282              : !!
    1283              : !! OUTPUT
    1284              : !!  (only writing)
    1285              : !!
    1286              : !! SOURCE
    1287              : 
    1288        17265 : subroutine ab_define_var(ncid, var_dim_id, var_id, var_type, var_name, var_mnemo, var_units)
    1289              : 
    1290              : !Arguments ------------------------------------
    1291              : !scalars
    1292              :  integer, intent(in) :: ncid
    1293              :  integer, intent(out) :: var_id
    1294              :  character(len=*), intent(in) :: var_mnemo,var_units,var_name
    1295              :  integer,intent(in) :: var_type
    1296              : !arrays
    1297              :  integer,intent(in) :: var_dim_id(:)
    1298              : 
    1299              : !Local variables-------------------------------
    1300              :  integer :: ncerr
    1301              : ! *************************************************************************
    1302              : 
    1303        17265 :  ncerr = nf90_def_var(ncid, trim(var_name), var_type, var_dim_id, var_id)
    1304        17265 :  NCF_CHECK_MSG(ncerr," define variable "//trim(var_name))
    1305              : 
    1306        17265 :  ncerr = nf90_put_att(ncid, var_id,  "units",trim(var_units))
    1307        17265 :  NCF_CHECK_MSG(ncerr," define attribute for "//trim(var_name))
    1308              : 
    1309        17265 :  ncerr = nf90_put_att(ncid, var_id,  "mnemonics", trim(var_mnemo))
    1310        17265 :  NCF_CHECK_MSG(ncerr," define attribute for "//trim(var_name))
    1311              : 
    1312        17265 : end subroutine ab_define_var
    1313              : !!***
    1314              : 
    1315              : !!****f* m_nctk/nctk_def_scalars_type
    1316              : !! NAME
    1317              : !!  nctk_def_scalars_type
    1318              : !!
    1319              : !! FUNCTION
    1320              : !!  Define list of **scalar** variables with a given type, Return immediately if error
    1321              : !!
    1322              : !! INPUTS
    1323              : !!  ncid=Netcdf identifier.
    1324              : !!  varnames(:)=List of strings with the name of the variables
    1325              : !!  xtype=Type of the variables
    1326              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1327              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1328              : !!
    1329              : !! SOURCE
    1330              : 
    1331       216039 : integer function nctk_def_scalars_type(ncid, varnames, xtype, defmode, prefix) result(ncerr)
    1332              : 
    1333              : !Arguments ------------------------------------
    1334              : !scalars
    1335              :  integer,intent(in) :: ncid,xtype
    1336              :  logical,optional,intent(in) :: defmode
    1337              :  character(len=*),optional,intent(in) :: prefix
    1338              : !arrays
    1339              :  character(len=*),intent(in) :: varnames(:)
    1340              : 
    1341              : !Local variables-------------------------------
    1342              : !scalars
    1343              :  integer :: ii,varid
    1344              :  character(len=nctk_slen) :: vname
    1345              :  type(nctkvar_t) :: var
    1346              : ! *********************************************************************
    1347              : 
    1348       216039 :  ncerr = nf90_noerr
    1349       216039 :  if (present(defmode)) then
    1350        39754 :    if (defmode) then
    1351        39754 :      NCF_CHECK(nctk_set_defmode(ncid))
    1352              :    end if
    1353              :  end if
    1354              : 
    1355              :  ! Special case where dimension is null
    1356      1153814 :  do ii=1,size(varnames)
    1357       937775 :    vname = varnames(ii)
    1358       937775 :    if (present(prefix)) vname = strcat(prefix, vname)
    1359              : 
    1360      1153814 :    if (nf90_inq_varid(ncid, vname, varid) == nf90_noerr)  then
    1361              :        ! Variable already exists. Check if type and dimensions agree
    1362        98998 :        call var_from_id(ncid, varid, var)
    1363              : 
    1364        98998 :        if (.not. (var%xtype == xtype .and. var%ndims == 0)) then
    1365            0 :          ABI_ERROR("variable already exists with a different definition.")
    1366              :        else
    1367              :          cycle ! Dimension matches, skip definition.
    1368              :        end if
    1369              : 
    1370              :    else
    1371              :      ! Define variable since it doesn't exist.
    1372       838777 :      ncerr = nf90_def_var(ncid, vname, xtype, varid)
    1373       838777 :      NCF_CHECK(ncerr)
    1374              :    end if
    1375              :  end do
    1376              : 
    1377    442663911 : end function nctk_def_scalars_type
    1378              : !!***
    1379              : 
    1380              : !!****f* m_nctk/nctk_def_iscalars
    1381              : !! NAME
    1382              : !!  nctk_def_iscalars
    1383              : !!
    1384              : !! FUNCTION
    1385              : !!  Define list of integer **scalar** variables. Return immediately if error
    1386              : !!
    1387              : !! INPUTS
    1388              : !!  ncid=Netcdf identifier.
    1389              : !!  varnames(:)=List of strings with the name of the variables
    1390              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1391              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1392              : !!
    1393              : !! SOURCE
    1394              : 
    1395       112130 : integer function nctk_def_iscalars(ncid, varnames, defmode, prefix) result(ncerr)
    1396              : 
    1397              : !Arguments ------------------------------------
    1398              : !scalars
    1399              :  integer,intent(in) :: ncid
    1400              :  logical,optional,intent(in) :: defmode
    1401              :  character(len=*),optional,intent(in) :: prefix
    1402              : !arrays
    1403              :  character(len=*),intent(in) :: varnames(:)
    1404              : ! *********************************************************************
    1405              : 
    1406       112130 :  if (present(defmode)) then
    1407        28868 :    ncerr = nctk_def_scalars_type(ncid, varnames, nf90_int, defmode=defmode)
    1408              :  else
    1409        83262 :    if (present(prefix)) then
    1410           36 :      ncerr = nctk_def_scalars_type(ncid, varnames, nf90_int, prefix=prefix)
    1411              :    else
    1412        83226 :      ncerr = nctk_def_scalars_type(ncid, varnames, nf90_int)
    1413              :    end if
    1414              :  end if
    1415              : 
    1416       112130 : end function nctk_def_iscalars
    1417              : !!***
    1418              : 
    1419              : !!****f* m_nctk/nctk_def_dpscalars
    1420              : !! NAME
    1421              : !!  nctk_def_dpscalars
    1422              : !!
    1423              : !! FUNCTION
    1424              : !!  Define list of double precision **scalar** variables. Return immediately if error
    1425              : !!
    1426              : !! INPUTS
    1427              : !!  ncid=Netcdf identifier.
    1428              : !!  varnames(:)=List of strings with the name of the variables
    1429              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1430              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1431              : !!
    1432              : !! SOURCE
    1433              : 
    1434       207818 : integer function nctk_def_dpscalars(ncid, varnames, defmode, prefix) result(ncerr)
    1435              : 
    1436              : !Arguments ------------------------------------
    1437              : !scalars
    1438              :  integer,intent(in) :: ncid
    1439              :  logical,optional,intent(in) :: defmode
    1440              :  character(len=*),optional,intent(in) :: prefix
    1441              : !arrays
    1442              :  character(len=*),intent(in) :: varnames(:)
    1443              : 
    1444              : !Local variables-------------------------------
    1445              :  character(len=nctk_slen) :: prefix_
    1446              : ! *********************************************************************
    1447       103909 :  prefix_ = ""; if (present(prefix)) prefix_ = prefix
    1448              : 
    1449       103909 :  if (present(defmode)) then
    1450        10886 :    ncerr = nctk_def_scalars_type(ncid, varnames, nf90_double, defmode=defmode, prefix=prefix_)
    1451              :  else
    1452        93023 :    ncerr = nctk_def_scalars_type(ncid, varnames, nf90_double, prefix=prefix_)
    1453              :  end if
    1454              : 
    1455       103909 : end function nctk_def_dpscalars
    1456              : !!***
    1457              : 
    1458              : !!****f* m_nctk/nctk_def_one_array
    1459              : !! NAME
    1460              : !!  nctk_def_one_array
    1461              : !!
    1462              : !! FUNCTION
    1463              : !!  Define list of arrays with a given type, Return immediately if error
    1464              : !!
    1465              : !! INPUTS
    1466              : !!  ncid=Netcdf identifier.
    1467              : !!  nctk_array=Array descriptor.
    1468              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1469              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1470              : !!
    1471              : !! SOURCE
    1472              : 
    1473              : 
    1474      1200687 : integer function nctk_def_one_array(ncid, nctk_array, defmode, varid, prefix) result(ncerr)
    1475              : 
    1476              : !Arguments ------------------------------------
    1477              : !scalars
    1478              :  integer,intent(in) :: ncid
    1479              :  logical,optional,intent(in) :: defmode
    1480              :  integer,optional,intent(out) :: varid
    1481              :  type(nctkarr_t),intent(in) :: nctk_array
    1482              :  character(len=*),optional,intent(in) :: prefix
    1483              : 
    1484              : !Local variables-------------------------------
    1485              : !scalars
    1486              :  integer :: ii,xtype,prev,cnt,nn,vid
    1487              :  character(len=500) :: msg
    1488              : !arrays
    1489              :  integer :: dimids(NF90_MAX_DIMS),dimvals(NF90_MAX_DIMS)
    1490              :  character(len=nctk_slen) :: sarr(NF90_MAX_DIMS), string, pre, vname, dimname
    1491              :  type(nctkvar_t) :: var
    1492              : ! *********************************************************************
    1493              : 
    1494           90 :  pre = ""; if (present(prefix)) pre = prefix
    1495              : 
    1496      1200687 :  ncerr = nf90_noerr
    1497      1200687 :  if (present(defmode)) then
    1498            0 :    if (defmode) then
    1499            0 :      NCF_CHECK(nctk_set_defmode(ncid))
    1500              :    end if
    1501              :  end if
    1502              : 
    1503      1200687 :  xtype = str2xtype(nctk_array%dtype)
    1504      1200687 :  vname = strcat(pre, nctk_array%name)
    1505              : 
    1506              :  ! Build array of strings with the dimensions.
    1507      1200687 :  string = nctk_array%shape_str
    1508      1200687 :  nn = char_count(string, ",")
    1509      1200687 :  ABI_CHECK(nn <= NF90_MAX_DIMS, "Too many dimensions!")
    1510              : 
    1511              :  ! Parse dimension names and add prefix (if any).
    1512      1200687 :  if (nn == 0) then
    1513       557427 :    cnt = 1
    1514       557427 :    dimname = lstrip(string)
    1515      6054447 :    if (any(dimname == NCTK_IMPLICIT_DIMS)) then
    1516        10936 :      sarr(1) = dimname
    1517              :    else
    1518       546491 :      sarr(1) = strcat(pre, dimname)
    1519              :    end if
    1520              :  else
    1521       643260 :    prev = 0; cnt = 0
    1522     34811649 :    do ii=1,len_trim(string)
    1523     34811649 :      if (string(ii:ii) == ",") then
    1524       981316 :        cnt = cnt + 1
    1525       981316 :        dimname = lstrip(string(prev+1:ii-1))
    1526     10062342 :        if (any(dimname == NCTK_IMPLICIT_DIMS)) then
    1527        93627 :          sarr(cnt) = dimname
    1528              :        else
    1529       887689 :          sarr(cnt) = strcat(pre, dimname)
    1530              :        end if
    1531              :        prev = ii
    1532              :      end if
    1533              :    end do
    1534       643260 :    cnt = cnt + 1
    1535       643260 :    dimname = lstrip(string(prev+1:ii-1))
    1536      7073031 :    if (any(dimname == NCTK_IMPLICIT_DIMS)) then
    1537          385 :      sarr(cnt) = dimname
    1538              :    else
    1539       642875 :      sarr(cnt) = strcat(pre, dimname)
    1540              :    end if
    1541              :  end if
    1542              : 
    1543              :  ! Get dimids
    1544      3382690 :  do ii=1,cnt
    1545      2182003 :    NCF_CHECK_MSG(nf90_inq_dimid(ncid, sarr(ii), dimids(ii)), sarr(ii))
    1546      3382690 :    NCF_CHECK(nf90_inquire_dimension(ncid, dimids(ii), len=dimvals(ii)))
    1547              :  end do
    1548              : 
    1549              :  ! Check if dimension already exists.
    1550              :  ! Variable already exists. Check if type and dimensions agree
    1551      1200687 :  if (nf90_inq_varid(ncid, vname, vid) == nf90_noerr)  then
    1552       225639 :    call var_from_id(ncid, vid, var)
    1553       225639 :    if (.not. (var%xtype == xtype .and. var%ndims == cnt)) then
    1554              :       write(msg,"(4a,2(2(a,i0),a))")&
    1555            0 :         "variable ",trim(vname)," already exists with a different definition:",ch10,&
    1556            0 :         "In file:     xtype = ",var%xtype,", ndims = ",var%ndims,ch10,&
    1557            0 :         "From caller: xtype = ",xtype,", ndims = ",cnt,ch10
    1558            0 :       ABI_ERROR(msg)
    1559              :    end if
    1560       641401 :    if (any(dimvals(1:cnt) /= var%dimlens(1:var%ndims))) then
    1561              :       write(msg,"(4a,2(3a))")&
    1562            0 :         "variable ",trim(vname)," already exists but with different shape.",ch10,&
    1563            0 :         "In file:     dims = ",trim(ltoa(var%dimlens(:var%ndims))),ch10,&
    1564            0 :         "From caller  dims = ",trim(ltoa(dimvals(:cnt))),ch10
    1565            0 :       ABI_ERROR(msg)
    1566              :    end if
    1567       225639 :    if (present(varid)) varid = vid
    1568       225639 :    return
    1569              :  end if
    1570              : 
    1571              :  ! Define variable since it doesn't exist.
    1572       975048 :  ncerr = nf90_def_var(ncid, vname, xtype, dimids(1:cnt), vid)
    1573       975048 :  NCF_CHECK(ncerr)
    1574              : 
    1575       975048 :  if (present(varid)) varid = vid
    1576              : 
    1577   2461182711 : end function nctk_def_one_array
    1578              : !!***
    1579              : 
    1580              : !!****f* m_nctk/nctk_def_array_list
    1581              : !! NAME
    1582              : !!  nctk_def_array_list
    1583              : !!
    1584              : !! FUNCTION
    1585              : !!  Define list of arrays with a given type, Return immediately if error
    1586              : !!
    1587              : !! INPUTS
    1588              : !!  ncid=Netcdf identifier.
    1589              : !!  nctk_arrays(:)=List of array descriptors.
    1590              : !!  [defmode]=If True, the nc file is set in define mode (default=False)
    1591              : !!  [prefix]=Prefix added to varnames and dimensions. Empty string if not specified.
    1592              : !!
    1593              : !! SOURCE
    1594              : 
    1595              : 
    1596       247507 : integer function nctk_def_array_list(ncid, nctk_arrays, defmode, prefix) result(ncerr)
    1597              : 
    1598              : !Arguments ------------------------------------
    1599              : !scalars
    1600              :  integer,intent(in) :: ncid
    1601              :  logical,optional,intent(in) :: defmode
    1602              :  character(len=*),optional,intent(in) :: prefix
    1603              : !arrays
    1604              :  type(nctkarr_t),intent(in) :: nctk_arrays(:)
    1605              : 
    1606              : !Local variables-------------------------------
    1607              :  integer :: ia
    1608              : ! *********************************************************************
    1609              : 
    1610       247507 :  ncerr = nf90_noerr
    1611       247507 :  if (present(defmode)) then
    1612         4364 :    if (defmode) then
    1613         4364 :      NCF_CHECK(nctk_set_defmode(ncid))
    1614              :    end if
    1615              :  end if
    1616              : 
    1617      1334745 :  do ia=1,size(nctk_arrays)
    1618      1334745 :    if (present(prefix)) then
    1619           90 :      NCF_CHECK(nctk_def_one_array(ncid, nctk_arrays(ia), prefix=prefix))
    1620              :    else
    1621      1087148 :      NCF_CHECK(nctk_def_one_array(ncid, nctk_arrays(ia)))
    1622              :    end if
    1623              :  end do
    1624              : 
    1625       247507 : end function nctk_def_array_list
    1626              : !!***
    1627              : 
    1628              : !!****f* m_nctk/nctk_write_iscalars
    1629              : !! NAME
    1630              : !!  nctk_write_iscalars
    1631              : !!
    1632              : !! FUNCTION
    1633              : !!  Write a list of **scalar** integer variables. Return immediately if error
    1634              : !!
    1635              : !! INPUTS
    1636              : !!  ncid=Netcdf identifier.
    1637              : !!  varnames(:)=List of strings with the name of the variables
    1638              : !!  values(:)=List of integer scalars
    1639              : !!  [datamode]=If True, the nc file is set in data mode (default=False)
    1640              : !!
    1641              : !! OUTPUT
    1642              : !!  ncerr=Exit status
    1643              : !!
    1644              : !! SOURCE
    1645              : 
    1646        40733 : integer function nctk_write_iscalars(ncid, varnames, values, datamode) result(ncerr)
    1647              : 
    1648              : !Arguments ------------------------------------
    1649              : !scalars
    1650              :  integer,intent(in) :: ncid
    1651              :  logical,optional,intent(in) :: datamode
    1652              : !arrays
    1653              :  integer,intent(in) :: values(:)
    1654              :  character(len=*),intent(in) :: varnames(:)
    1655              : 
    1656              : !Local variables-------------------------------
    1657              :  integer :: ii,varid
    1658              : ! *********************************************************************
    1659              : 
    1660        40733 :  ABI_CHECK(size(varnames) == size(values), "Different size in varnames, values")
    1661              : 
    1662        40733 :  ncerr = nf90_noerr
    1663        40733 :  if (present(datamode)) then
    1664           14 :    if (datamode) then
    1665           14 :      NCF_CHECK(nctk_set_datamode(ncid))
    1666              :    end if
    1667              :  end if
    1668              : 
    1669       245006 :  do ii=1,size(varnames)
    1670       204273 :    NCF_CHECK_MSG(nf90_inq_varid(ncid, varnames(ii), varid), sjoin("Inquiring: ", varnames(ii)))
    1671       245006 :    NCF_CHECK(nf90_put_var(ncid, varid, values(ii)))
    1672              :  end do
    1673              : 
    1674        40733 : end function nctk_write_iscalars
    1675              : !!***
    1676              : 
    1677              : !!****f* m_nctk/nctk_write_dpscalars
    1678              : !! NAME
    1679              : !!  nctk_write_dpscalars
    1680              : !!
    1681              : !! FUNCTION
    1682              : !!  Write a list of **scalar** real(dp) variables. Return immediately if error
    1683              : !!
    1684              : !! INPUTS
    1685              : !!  ncid=Netcdf identifier.
    1686              : !!  varnames(:)=List of strings with the name of the variables
    1687              : !!  values(:)=List of real(dp) scalars
    1688              : !!  [datamode]=If True, the nc file is set in data mode (default=False)
    1689              : !!
    1690              : !! SOURCE
    1691              : 
    1692        37328 : integer function nctk_write_dpscalars(ncid, varnames, values, datamode) result(ncerr)
    1693              : 
    1694              : !Arguments ------------------------------------
    1695              : !scalars
    1696              :  integer,intent(in) :: ncid
    1697              :  logical,optional,intent(in) :: datamode
    1698              : !arrays
    1699              :  real(dp),intent(in) :: values(:)
    1700              :  character(len=*),intent(in) :: varnames(:)
    1701              : 
    1702              : !Local variables-------------------------------
    1703              :  integer :: ii,varid
    1704              : ! *********************************************************************
    1705              : 
    1706        37328 :  ncerr = nf90_noerr
    1707              : 
    1708        37328 :  ABI_CHECK(size(varnames) == size(values), "Different size in varnames, values")
    1709              : 
    1710        37328 :  if (present(datamode)) then
    1711         4147 :    if (datamode) then
    1712         4147 :      NCF_CHECK(nctk_set_datamode(ncid))
    1713              :    end if
    1714              :  end if
    1715              : 
    1716       269719 :  do ii=1,size(varnames)
    1717       232391 :    NCF_CHECK(nf90_inq_varid(ncid, varnames(ii), varid))
    1718       269719 :    NCF_CHECK(nf90_put_var(ncid, varid, values(ii)))
    1719              :  end do
    1720              : 
    1721        37328 : end function nctk_write_dpscalars
    1722              : !!***
    1723              : 
    1724              : !!****f* m_nctk/nctk_defnwrite_ivars
    1725              : !! NAME
    1726              : !!  nctk_defnwrite_ivars
    1727              : !!
    1728              : !! FUNCTION
    1729              : !!  Define list of integer **scalar** variables and write their values.
    1730              : !!  Return immediately if error
    1731              : !!
    1732              : !! INPUTS
    1733              : !!  ncid=Netcdf identifier.
    1734              : !!  varnames(:)=List of strings with the name of the variables
    1735              : !!  values(:)=List of integer scalars
    1736              : !!
    1737              : !! SOURCE
    1738              : 
    1739        13554 : integer function nctk_defnwrite_ivars(ncid, varnames, values) result(ncerr)
    1740              : 
    1741              : !Arguments ------------------------------------
    1742              : !scalars
    1743              :  integer,intent(in) :: ncid
    1744              : !arrays
    1745              :  integer,intent(in) :: values(:)
    1746              :  character(len=*),intent(in) :: varnames(:)
    1747              : 
    1748              : !Local variables-------------------------------
    1749              : !scalars
    1750              :  integer :: ii,varid
    1751              : ! *********************************************************************
    1752              : 
    1753        13554 :  ABI_CHECK(size(varnames) == size(values), "Different size in varnames, values")
    1754              : 
    1755        13554 :  ncerr = nctk_def_iscalars(ncid, varnames, defmode=.True.)
    1756        13554 :  NCF_CHECK(ncerr)
    1757              : 
    1758        13554 :  NCF_CHECK(nctk_set_datamode(ncid))
    1759        39259 :  do ii=1,size(varnames)
    1760        25705 :    varid = nctk_idname(ncid, varnames(ii))
    1761        39259 :    NCF_CHECK(nf90_put_var(ncid, varid, values(ii)))
    1762              :  end do
    1763              : 
    1764        13554 : end function nctk_defnwrite_ivars
    1765              : !!***
    1766              : 
    1767              : !!****f* m_nctk/nctk_defnwrite_dpvars
    1768              : !! NAME
    1769              : !!  nctk_defnwrite_dpvars
    1770              : !!
    1771              : !! FUNCTION
    1772              : !!  Define list of real(dp) **scalar** variables and write their values.
    1773              : !!  Return immediately if error
    1774              : !!
    1775              : !! INPUTS
    1776              : !!  ncid=Netcdf identifier.
    1777              : !!  varnames(:)=List of strings with the name of the variables
    1778              : !!  values(:)=List of integer scalars
    1779              : !!
    1780              : !! SOURCE
    1781              : 
    1782        10886 : integer function nctk_defnwrite_dpvars(ncid, varnames, values) result(ncerr)
    1783              : 
    1784              : !Arguments ------------------------------------
    1785              : !scalars
    1786              :  integer,intent(in) :: ncid
    1787              : !arrays
    1788              :  real(dp),intent(in) :: values(:)
    1789              :  character(len=*),intent(in) :: varnames(:)
    1790              : 
    1791              : !Local variables-------------------------------
    1792              : !scalars
    1793              :  integer :: ii,varid
    1794              : !arrays
    1795              : ! *********************************************************************
    1796              :  ncerr = nf90_noerr
    1797              : 
    1798        10886 :  ABI_CHECK(size(varnames) == size(values), "Different size in varnames, values")
    1799              : 
    1800        10886 :  ncerr = nctk_def_dpscalars(ncid, varnames, defmode=.True.)
    1801        10886 :  NCF_CHECK(ncerr)
    1802              : 
    1803        10886 :  NCF_CHECK(nctk_set_datamode(ncid))
    1804       244368 :  do ii=1,size(varnames)
    1805              :    !write(std_out,*)varnames(ii)
    1806       233482 :    varid = nctk_idname(ncid, varnames(ii))
    1807       244368 :    NCF_CHECK(nf90_put_var(ncid, varid, values(ii)))
    1808              :  end do
    1809              : 
    1810        10886 : end function nctk_defnwrite_dpvars
    1811              : !!***
    1812              : 
    1813              : !!****f* m_nctk/nctk_write_ibz
    1814              : !! NAME
    1815              : !!  nctk_write_ibz
    1816              : !!
    1817              : !! FUNCTION
    1818              : !!  Write the list of the k-points in the IBZ with the corresponding weights in Netcdf format.
    1819              : !!  Mainly used for passing data to AbiPy. This routine should be called by master only.
    1820              : !!
    1821              : !! INPUTS
    1822              : !!  fname=File name
    1823              : !!  kpoints(:,:)=List of k-points
    1824              : !!  weights(:)=K-point weights
    1825              : !!
    1826              : !! OUTPUT
    1827              : !!  ncerr=Exit status
    1828              : !!
    1829              : !! SOURCE
    1830              : 
    1831            0 : integer function nctk_write_ibz(fname, kpoints, weights) result(ncerr)
    1832              : 
    1833              : !Arguments ------------------------------------
    1834              : !scalars
    1835              :  character(len=*),intent(in) :: fname
    1836              : !arrays
    1837              :  real(dp),intent(in) :: kpoints(:,:),weights(:)
    1838              : 
    1839              : !Local variables-------------------------------
    1840              :  integer :: nkpts,ncid
    1841              : ! *********************************************************************
    1842              : 
    1843            0 :  ABI_CHECK(size(kpoints, dim=2) == size(weights), "size(kpoints, dim=2) != size(weights)")
    1844            0 :  nkpts = size(kpoints, dim=2)
    1845              : 
    1846            0 :  NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), sjoin("Creating:", fname))
    1847              : 
    1848              :  ncerr = nctk_def_dims(ncid, [ &
    1849            0 :    nctkdim_t("number_of_reduced_dimensions",3), nctkdim_t("number_of_kpoints", nkpts)], defmode=.True.)
    1850            0 :  NCF_CHECK(ncerr)
    1851              : 
    1852              :  ncerr = nctk_def_array_list(ncid, [&
    1853              :    nctkarr_t('reduced_coordinates_of_kpoints', "dp", "number_of_reduced_dimensions, number_of_kpoints"),&
    1854            0 :    nctkarr_t('kpoint_weights', "dp", "number_of_kpoints")])
    1855              :  NCF_CHECK(ncerr)
    1856              : 
    1857            0 :  NCF_CHECK(nctk_set_datamode(ncid))
    1858              : 
    1859            0 :  ncerr = nf90_put_var(ncid, nctk_idname(ncid, 'reduced_coordinates_of_kpoints'), kpoints)
    1860            0 :  NCF_CHECK(ncerr)
    1861              : 
    1862            0 :  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'kpoint_weights'), weights))
    1863              : 
    1864            0 :  NCF_CHECK(nf90_close(ncid))
    1865              : 
    1866            0 : end function nctk_write_ibz
    1867              : !!***
    1868              : 
    1869              : !!****f* m_nctk/nctk_get_dim
    1870              : !! NAME
    1871              : !!  nctk_get_dim
    1872              : !!
    1873              : !! FUNCTION
    1874              : !!  Get the value of a dimension from its name.
    1875              : !!
    1876              : !! INPUTS
    1877              : !!  ncid=Netcdf identifier.
    1878              : !!  dimname=Name of the dimension.
    1879              : !!  [datamode]=If True, the nc file is set in data mode (default=False)
    1880              : !!
    1881              : !! OUTPUT
    1882              : !!  dimlen=Value of the dimension.
    1883              : !!
    1884              : !! SOURCE
    1885              : 
    1886       311234 : integer function nctk_get_dim(ncid, dimname, dimlen, datamode) result(ncerr)
    1887              : 
    1888              : !Arguments ------------------------------------
    1889              : !scalars
    1890              :  integer,intent(in) :: ncid
    1891              :  character(len=*),intent(in) :: dimname
    1892              :  integer,intent(out) :: dimlen
    1893              :  logical,optional,intent(in) :: datamode
    1894              : 
    1895              : !Local variables-------------------------------
    1896              :  integer :: dimid
    1897              : ! *********************************************************************
    1898              : 
    1899       311234 :  ncerr = nf90_noerr
    1900              : 
    1901       311234 :  if (present(datamode)) then
    1902            2 :    if (datamode) then
    1903            2 :      NCF_CHECK(nctk_set_datamode(ncid))
    1904              :    end if
    1905              :  end if
    1906              : 
    1907       311234 :  ncerr = nf90_inq_dimid(ncid, dimname, dimid)
    1908       311234 :  if (ncerr == nf90_noerr) then
    1909       311234 :    ncerr = nf90_inquire_dimension(ncid, dimid, len=dimlen)
    1910              :  end if
    1911              : 
    1912       311234 : end function nctk_get_dim
    1913              : !!***
    1914              : 
    1915              : !----------------------------------------------------------------------
    1916              : 
    1917              : !!****f* m_nctk/nctk_write_datar
    1918              : !! NAME
    1919              : !! nctk_write_datar
    1920              : !!
    1921              : !! FUNCTION
    1922              : !!  Write an array in real space in netcdf format
    1923              : !!
    1924              : !! INPUTS
    1925              : !!  path=Filename
    1926              : !!  varname=Name of the variable to write.
    1927              : !!  ngfft(18)=information about 3D FFT
    1928              : !!  cplex=1 for real arrays (e.g. GS rhor), 2 for complex array.
    1929              : !!  nfft=number of points in the real space FFT mesh treated by this MPI proc
    1930              : !!  nspden=number of spin-density components
    1931              : !!  comm_fft=MPI communicator (used only if MPI-FFT).
    1932              : !!  fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
    1933              : !!  ffti3_local(n3)=local index for 3d dimension
    1934              : !!  datar(cplex*nfft,nspden)= array in real space.
    1935              : !!  [action]
    1936              : !!
    1937              : !! OUTPUT
    1938              : !!  Only writing
    1939              : !!
    1940              : !! SOURCE
    1941              : 
    1942         7672 : integer function nctk_write_datar(varname,path,ngfft,cplex,nfft,nspden,&
    1943         7672 :    comm_fft,fftn3_distrib,ffti3_local,datar,action) result(ncerr)
    1944              : 
    1945              : !Arguments ------------------------------------
    1946              : !scalars
    1947              :  integer,intent(in) :: cplex,nfft,nspden,comm_fft
    1948              :  character(len=*),intent(in) :: path,varname
    1949              :  character(len=*),optional,intent(in) :: action
    1950              : !arrays
    1951              :  integer,intent(in) :: ngfft(18),fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
    1952              :  real(dp),target,intent(in) :: datar(cplex*nfft,nspden)
    1953              : 
    1954              : !Local variables-------------------------------
    1955              : !scalars
    1956              :  integer,parameter :: master=0
    1957              :  integer :: ncid,varid,i3,nproc_fft,me_fft,i3_glob,n1,n2,n3,ispden,cmode
    1958              :  logical :: ionode
    1959              :  character(len=nctk_slen) :: cplex_name,my_action
    1960              :  !character(len=500) :: msg
    1961              : !arrays
    1962         7672 :  real(dp),allocatable :: glob_datar(:,:)
    1963              : ! *************************************************************************
    1964              : 
    1965              :  ! FIXME: Default should be open but this enters into conflict with the abi_estf stuff!
    1966              :  ! if we are using MPI-IO since file is not open with HDF5.
    1967              :  !my_action = "create"; if (present(action)) my_action = action
    1968         7672 :  my_action = "open"; if (present(action)) my_action = action
    1969              : 
    1970         7672 :  nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
    1971         7672 :  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
    1972              : 
    1973              :  ! TODO: Be careful here because we should always create with HDF5 if available
    1974              :  ! to avoid problems if we have to reread with nproc_fft > 1 and MPI-IO
    1975         7672 :  ionode = .True.; ncerr = nf90_noerr
    1976         7672 :  if (nproc_fft == 1) then
    1977            0 :    select case(my_action)
    1978              :    case ("open")
    1979            0 :      ncerr = nf90_open(path, mode=nf90_write, ncid=ncid)
    1980              :    case ("create")
    1981         7414 :      ncerr = nctk_open_create(ncid, path, comm_fft)
    1982              :    case default
    1983         7414 :      ABI_ERROR(sjoin("Wrong action: ", my_action))
    1984              :    end select
    1985              : 
    1986              :  else
    1987          258 :    if (nctk_has_mpiio) then
    1988          258 :      call wrtout(std_out, strcat("nctk_write_datar: using MPI-IO to write ", varname, path), "COLL")
    1989              : 
    1990          258 :      ncerr = nf90_einval
    1991              : #ifdef HAVE_NETCDF_MPI
    1992            0 :      select case(my_action)
    1993              :      case ("open")
    1994            0 :        ncerr = nf90_open(path, mode=nf90_write, comm=comm_fft, info=xmpio_info, ncid=ncid)
    1995              :      case ("create")
    1996              :        ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
    1997          258 :          comm=comm_fft, info=xmpio_info, ncid=ncid)
    1998              :      case default
    1999          258 :        ABI_ERROR(strcat("Wrong action:", my_action))
    2000              :      end select
    2001              : #endif
    2002              :    else
    2003              :      ! MPI-FFT without MPI-support. Only master performs IO
    2004            0 :      ionode = (me_fft == master)
    2005            0 :      if (ionode) then
    2006            0 :        select case(my_action)
    2007              :        case ("open")
    2008            0 :          ncerr = nf90_open(path, mode=nf90_write, ncid=ncid)
    2009              :        case ("create")
    2010              :          !ncerr = nf90_create(path, cmode=nf90_clobber, ncid=ncid)
    2011            0 :          cmode = def_cmode_for_seq_create
    2012            0 :          ncerr = nf90_create(path, cmode=cmode, ncid=ncid)
    2013              :        case default
    2014            0 :          ABI_ERROR(strcat("Wrong action:", my_action))
    2015              :        end select
    2016              :      end if
    2017              :    end if
    2018              :  end if
    2019         7672 :  NCF_CHECK_MSG(ncerr, sjoin("opening file:", path))
    2020              : 
    2021         7672 :  if (ionode) then
    2022              :    ! Define dims and variables
    2023              :    !write(std_out,*)"defing dims",trim(varname)," in file: ",path
    2024         7672 :    cplex_name = strcat("real_or_complex_", varname)
    2025              :    ncerr = nctk_def_dims(ncid, [&
    2026              :      nctkdim_t(cplex_name, cplex),&
    2027              :      nctkdim_t("number_of_grid_points_vector1", n1),&
    2028              :      nctkdim_t("number_of_grid_points_vector2", n2),&
    2029              :      nctkdim_t("number_of_grid_points_vector3", n3),&
    2030        46032 :      nctkdim_t("number_of_components", nspden)], defmode=.True.)
    2031         7672 :    NCF_CHECK(ncerr)
    2032              : 
    2033              :    ncerr = nctk_def_one_array(ncid, nctkarr_t(name=varname, dtype="dp", shape_str=strcat(cplex_name, &
    2034              : ", number_of_grid_points_vector1, number_of_grid_points_vector2, number_of_grid_points_vector3, number_of_components")),&
    2035         7672 :    varid=varid)
    2036              : 
    2037              :    ! Add attributes
    2038         7672 :    varid = nctk_idname(ncid, varname)
    2039         7672 :    NCF_CHECK(nf90_put_att(ncid, varid, "units", "atomic units"))
    2040         7672 :    NCF_CHECK(nf90_put_att(ncid, varid, "scale_to_atomic_units", one))
    2041              :  end if
    2042              : 
    2043         7672 :  if (nproc_fft == 1) then
    2044              :    ! no MPI-FFT --> write data directly.
    2045         7414 :    varid = nctk_idname(ncid, varname)
    2046         7414 :    NCF_CHECK(nctk_set_datamode(ncid))
    2047        44484 :    NCF_CHECK(nf90_put_var(ncid, varid, datar, start=[1,1,1,1,1], count=[cplex, n1, n2, n3, nspden]))
    2048         7414 :    NCF_CHECK(nf90_close(ncid))
    2049              : 
    2050              :  else
    2051              :    ! Must handle data distribution.
    2052          258 :    ABI_CHECK(mod(n3, nproc_fft) == 0, "assuming mod(n3, nproc_fft) == 0")
    2053              : 
    2054          258 :    i3_glob = -1
    2055         2814 :    do i3=1,ngfft(3)
    2056         2814 :      if (fftn3_distrib(i3) == me_fft) then
    2057              :         i3_glob = i3
    2058              :         exit
    2059              :      end if
    2060              :    end do
    2061          258 :    ABI_CHECK(i3_glob > 0, "negative i3_glob")
    2062              : 
    2063              :    !do i3=i3_glob,ngfft(3)
    2064              :    !  if (fftn3_distrib(i3) /= me_fft) exit
    2065              :    !  !assert all(ffn3_distrib(i3_glob:i3_glob -1 + ngfft(3) / nproc_fft) == me_fft)
    2066              :    !end do
    2067              :    !i3_glob = i3- 1
    2068              :    !print*,"i3_glob",i3_glob
    2069              : 
    2070          258 :    if (ionode) then
    2071          258 :      NCF_CHECK(nf90_enddef(ncid))
    2072              :    end if
    2073              : 
    2074              :    ! Array on disk has shape [cplex, n1, n2, n3, nspden]
    2075          258 :    if (nctk_has_mpiio) then
    2076              :      ! Use collective IO.
    2077          258 :      ncerr = nf90_einval
    2078          258 :      NCF_CHECK(nctk_set_collective(ncid, varid))
    2079              : 
    2080          518 :      do ispden=1,nspden
    2081              :        ncerr = nf90_put_var(ncid, varid, datar(:,ispden), start=[1,1,1,i3_glob,ispden], &
    2082         2860 :                 count=[cplex,n1,n2,n3/nproc_fft,1])
    2083          518 :        NCF_CHECK(ncerr)
    2084              :      end do
    2085              :    else
    2086              :      ! MPI-FFT without MPI-IO. Collect data (requires more memory and communication)
    2087            0 :      ABI_MALLOC(glob_datar, (cplex*product(ngfft(1:3)), nspden))
    2088            0 :      call collect_datar(ngfft,cplex,nfft,nspden,datar,comm_fft,fftn3_distrib,ffti3_local,glob_datar,master=master)
    2089              : 
    2090            0 :      if (ionode) then
    2091              :        ! Write global array.
    2092            0 :        NCF_CHECK(nf90_put_var(ncid, varid, glob_datar, start=[1,1,1,1,1], count=[cplex,n1,n2,n3,nspden]))
    2093              :      end if
    2094            0 :      ABI_FREE(glob_datar)
    2095              :    end if
    2096              : 
    2097          258 :    if (ionode) then
    2098          258 :      NCF_CHECK(nf90_close(ncid))
    2099              :    end if
    2100              :  end if
    2101              : 
    2102              :  !ok = .True.
    2103              :  ! Sequential IO
    2104              :  !do rank=0,nproc_fft-1
    2105              :  !  if (rank == me_fft) then
    2106              :  !     ncerr = nf90_open(path, mode=nf90_write, ncid=ncid)
    2107              :  !     do ispden=1,nspden
    2108              :  !       ncerr = nf90_put_var(ncid, varid, datar(1:,ispden), start=[1,1,1,i3_glob,ispden], &
    2109              :  !                count=[cplex,ngfft(1),ngfft(2),ngfft(3)/nproc_fft,1])
    2110              :  !     end do
    2111              :  !     ncerr = nf90_close(ncid)
    2112              :  !  end if
    2113              :  !  call xmpi_barrier(comm_fft)
    2114              :  !end do
    2115              : 
    2116        15344 : end function nctk_write_datar
    2117              : !!***
    2118              : 
    2119              : !!****f* m_nctk/nctk_read_datar
    2120              : !! NAME
    2121              : !! nctk_read_datar
    2122              : !!
    2123              : !! FUNCTION
    2124              : !!  Read an array in real space in netcdf format
    2125              : !!
    2126              : !! INPUTS
    2127              : !!  path=Filename
    2128              : !!  varname=Name of the variable to read.
    2129              : !!  ngfft(18)=information about 3D FFT
    2130              : !!  cplex=1 for real arrays (e.g. GS rhor), 2 for complex array.
    2131              : !!  nfft=number of points in the real space FFT mesh treated by this MPI proc
    2132              : !!  nspden=number of spin-density components
    2133              : !!  comm_fft=MPI communicator (used only if MPI-FFT).
    2134              : !!  fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
    2135              : !!  ffti3_local(n3)=local index for 3d dimension
    2136              : !!  datar(cplex*nfft,nspden)= array in real space.
    2137              : !!
    2138              : !! OUTPUT
    2139              : !!  Only writing
    2140              : !!
    2141              : !! SOURCE
    2142              : 
    2143            0 : integer function nctk_read_datar(path,varname,ngfft,cplex,nfft,nspden,&
    2144            0 :    comm_fft,fftn3_distrib,ffti3_local,datar) result(ncerr)
    2145              : 
    2146              : !Arguments ------------------------------------
    2147              : !scalars
    2148              :  integer,intent(in) :: cplex,nfft,nspden,comm_fft
    2149              :  character(len=*),intent(in) :: path,varname
    2150              : !arrays
    2151              :  integer,intent(in) :: ngfft(18),fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
    2152              :  real(dp),intent(out) :: datar(cplex*nfft,nspden)
    2153              : 
    2154              : !Local variables-------------------------------
    2155              : !scalars
    2156              :  integer,parameter :: master=0
    2157              :  integer :: ncid,varid,i3,nproc_fft,me_fft,i3_glob,n1,n2,n3,ispden
    2158              :  logical :: ionode
    2159              : !arrays
    2160            0 :  real(dp),allocatable :: glob_datar(:,:)
    2161              : ! *************************************************************************
    2162              : 
    2163            0 :  nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
    2164            0 :  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
    2165              : 
    2166              :  ! TODO: Be careful here because we should always create with HDF5 if available
    2167              :  ! to avoid problems if we have to reread with nproc_fft > 1 and MPI-IO
    2168            0 :  ionode = .True.
    2169            0 :  if (nproc_fft == 1) then
    2170            0 :    ncerr = nf90_open(path, mode=nf90_nowrite, ncid=ncid)
    2171              :  else
    2172            0 :    if (nctk_has_mpiio) then
    2173              :      !write(std_out,*)"open_par: ",trim(path)
    2174              :      !ncerr = nf90_open_par(path, nf90_nowrite,
    2175              :      ! Don't know why but the format is not autodected!
    2176              :      ncerr = nf90_einval
    2177              : #ifdef HAVE_NETCDF_MPI
    2178              :      ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_nowrite),&
    2179            0 :                        comm=comm_fft, info=xmpio_info, ncid=ncid)
    2180              : #endif
    2181              :    else
    2182              :      ! MPI-FFT without MPI-support. Only master does IO
    2183            0 :      ionode = (me_fft == master); ncerr = nf90_noerr
    2184            0 :      if (ionode) ncerr = nf90_open(path, nf90_nowrite, ncid)
    2185              :    end if
    2186              :  end if
    2187            0 :  NCF_CHECK_MSG(ncerr, sjoin("opening file: ",path))
    2188              : 
    2189            0 :  NCF_CHECK(nf90_inq_varid(ncid, varname, varid))
    2190              :  !write(std_out,*)"about to read varname, ngfft, cplex, nfft, nspden:", trim(varname), ngfft(:3), cplex,nfft,nspden
    2191              : 
    2192              :  ! netcdf array has shape [cplex, n1, n2, n3, nspden]
    2193            0 :  if (nproc_fft == 1) then
    2194              :    ! No MPI-FFT --> easy
    2195            0 :    NCF_CHECK(nf90_get_var(ncid, varid, datar, start=[1,1,1,1], count=[cplex, n1, n2, n3, nspden]))
    2196            0 :    NCF_CHECK(nf90_close(ncid))
    2197              : 
    2198              :  else
    2199              :    ! Handle data distribution.
    2200            0 :    ABI_CHECK(mod(ngfft(3), nproc_fft) == 0, "assuming mod(n3, nproc_fft) == 0")
    2201              : 
    2202            0 :    i3_glob = -1
    2203            0 :    do i3=1,ngfft(3)
    2204            0 :      if (fftn3_distrib(i3) == me_fft) then
    2205              :         i3_glob = i3
    2206              :         exit
    2207              :      end if
    2208              :    end do
    2209            0 :    ABI_CHECK(i3_glob > 0, "negative i3_glob")
    2210              : 
    2211            0 :    if (nctk_has_mpiio) then
    2212              :      ! Use parallel IO with collective calls.
    2213            0 :      ncerr = nf90_einval
    2214            0 :      NCF_CHECK(nctk_set_collective(ncid, varid))
    2215              : 
    2216            0 :      do ispden=1,nspden
    2217              :        ncerr = nf90_get_var(ncid, varid, datar(:,ispden), start=[1,1,1,i3_glob,ispden], &
    2218            0 :                 count=[cplex,ngfft(1),ngfft(2),ngfft(3)/nproc_fft,1])
    2219            0 :        NCF_CHECK(ncerr)
    2220              :      end do
    2221              :    else
    2222              :      ! MPI-FFT without MPI-IO. Master read and broadcast (requires more memory and communication)
    2223            0 :      ABI_MALLOC(glob_datar, (cplex*product(ngfft(1:3)), nspden))
    2224            0 :      if (ionode) then
    2225            0 :        NCF_CHECK(nf90_get_var(ncid, varid, glob_datar, start=[1,1,1,1,1], count=[cplex,n1,n2,n3,nspden]))
    2226              :      end if
    2227              : 
    2228            0 :      call distrib_datar(ngfft,cplex,nfft,nspden,glob_datar,master,comm_fft,fftn3_distrib,ffti3_local,datar)
    2229            0 :      ABI_FREE(glob_datar)
    2230              :    end if
    2231              : 
    2232            0 :    if (ionode) then
    2233            0 :      NCF_CHECK(nf90_close(ncid))
    2234              :    end if
    2235              :  end if
    2236              : 
    2237            0 : end function nctk_read_datar
    2238              : !!***
    2239              : 
    2240              : !!****f* m_nctk/nctk_prepare_mpiio
    2241              : !! NAME
    2242              : !! nctk_mpiio
    2243              : !!
    2244              : !! FUNCTION
    2245              : !! This function appears to be required to prevent deadlocks during I/O operations in single mode.
    2246              : !! Although single mode is the default, on some architectures or compilers, nf90_put_var
    2247              : !! can deadlock if not all processors in the communicator invoke the function.
    2248              : !! This solution was proposed by Hsiao-Yi Tsai.
    2249              : 
    2250           43 : integer function nctk_prepare_mpiio(ncid, varname) result(ncerr)
    2251              : 
    2252              : !Arguments ------------------------------------
    2253              :  integer,intent(in) :: ncid
    2254              :  character(len=*),intent(in) :: varname
    2255              : 
    2256              : !Local variables-------------------------------
    2257              :  integer :: vid
    2258              :  character(len=nctk_slen) :: out_varname
    2259              : ! *************************************************************************
    2260              : 
    2261           43 :  vid = nctk_idname(ncid, varname)
    2262           43 :  ncerr = nf90_inquire_variable(ncid, vid, out_varname)
    2263              : 
    2264           43 : end function nctk_prepare_mpiio
    2265              : !!***
    2266              : 
    2267              : !----------------------------------------------------------------------
    2268              : 
    2269              : !!****f* m_nctk/collect_datar
    2270              : !! NAME
    2271              : !!  collect_datar
    2272              : !!
    2273              : !! FUNCTION
    2274              : !! Collect a real-space MPI-FFT distributed array on each proc.
    2275              : !!
    2276              : !! INPUTS
    2277              : !!  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
    2278              : !!  cplex=1 if real array, 2 for complex
    2279              : !!  nfft=Number of FFT points treated by this MPI proc
    2280              : !!  nspden=Second dimension of rhor
    2281              : !!  rhor(cplex*nfft,nspden)=Array in real space (MPI-FFT distributed)
    2282              : !!  fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
    2283              : !!  fftn3_local(n3)=local i3 indices
    2284              : !!  comm_fft=MPI-FFT communicator
    2285              : !!  [master]=MPI rank, Optional. If present, the global array is available only on master node.
    2286              : !!
    2287              : !! OUTPUT
    2288              : !!   rhor_glob(cplex*nfft_tot,nspden)=Global array
    2289              : !!
    2290              : !! SOURCE
    2291              : 
    2292            0 : subroutine collect_datar(ngfft,cplex,nfft,nspden,rhor,comm_fft,fftn3_distrib,ffti3_local,rhor_glob,master)
    2293              : 
    2294              : !Arguments ------------------------------------
    2295              : !scalars
    2296              :  integer,intent(in) :: cplex,nfft,nspden,comm_fft
    2297              :  integer,optional,intent(in) :: master
    2298              : !arrays
    2299              :  integer,intent(in) :: ngfft(18)
    2300              :  integer,intent(in) :: fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
    2301              :  real(dp),intent(in) :: rhor(cplex*nfft,nspden)
    2302              :  real(dp),intent(out) :: rhor_glob(cplex*product(ngfft(1:3)),nspden)
    2303              : 
    2304              : !Local variables-------------------------------
    2305              :  integer :: ispden,i1,i2,i3,me_fft,i3_local,my_fftbase,glob_fftbase
    2306              :  integer :: n1,n2,n3,ierr,nfft_tot
    2307              : ! *************************************************************************
    2308              : 
    2309            0 :  nfft_tot = product(ngfft(1:3)); me_fft = xmpi_comm_rank(comm_fft)
    2310            0 :  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
    2311              : 
    2312            0 :  if (nfft_tot == nfft) then
    2313              :    ! full rhor on each node, just do a copy
    2314            0 :    rhor_glob = rhor
    2315              :  else
    2316              :    ! if MPI-FFT we have to gather the global array on each node.
    2317            0 :    rhor_glob = zero
    2318            0 :    do ispden=1,nspden
    2319            0 :      do i3=1,n3
    2320            0 :        if (me_fft == fftn3_distrib(i3)) then
    2321            0 :          i3_local = ffti3_local(i3)
    2322            0 :          do i2=1,n2
    2323            0 :            my_fftbase =   cplex * ( (i2-1)*n1 + (i3_local-1)*n1*n2 )
    2324            0 :            glob_fftbase = cplex * ( (i2-1)*n1 + (i3-1)*n1*n2 )
    2325            0 :            do i1=1,cplex * n1
    2326            0 :              rhor_glob(i1+glob_fftbase,ispden) = rhor(i1+my_fftbase,ispden)
    2327              :            end do
    2328              :          end do
    2329              :        end if
    2330              :      end do
    2331              :    end do
    2332            0 :    if (present(master)) then
    2333            0 :      call xmpi_sum_master(rhor_glob,master,comm_fft,ierr)
    2334              :    else
    2335            0 :      call xmpi_sum(rhor_glob,comm_fft,ierr)
    2336              :    end if
    2337              :  end if
    2338              : 
    2339            0 : end subroutine collect_datar
    2340              : !!***
    2341              : 
    2342              : !----------------------------------------------------------------------
    2343              : 
    2344              : !!****f* m_nctk/distrib_datar
    2345              : !! NAME
    2346              : !!  distrib_datar
    2347              : !!
    2348              : !! FUNCTION
    2349              : !! distribute a real-space MPI-FFT
    2350              : !!
    2351              : !! INPUTS
    2352              : !!  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
    2353              : !!  cplex=1 if real array, 2 for complex
    2354              : !!  nfft=Number of FFT points treated by this MPI proc
    2355              : !!  nspden=Second dimension of rhor
    2356              : !!  rhor_glob(cplex*nfft_tot,nspden)=Global array
    2357              : !!  master=The rank of the node that owns the global array.
    2358              : !!  comm_fft=MPI-FFT communicator
    2359              : !!  fftn3_distrib(n3)=rank of the processors which own fft planes in 3rd dimension.
    2360              : !!  fftn3_local(n3)=local i3 indices
    2361              : !!
    2362              : !! OUTPUT
    2363              : !!  rhor(cplex*nfft,nspden)=Array in real space (MPI-FFT distributed)
    2364              : !!
    2365              : !! SOURCE
    2366              : 
    2367            0 : subroutine distrib_datar(ngfft,cplex,nfft,nspden,rhor_glob,master,comm_fft,fftn3_distrib,ffti3_local,rhor)
    2368              : 
    2369              : !Arguments ------------------------------------
    2370              : !scalars
    2371              :  integer,intent(in) :: cplex,nfft,nspden,comm_fft,master
    2372              : !arrays
    2373              :  integer,intent(in) :: ngfft(18)
    2374              :  integer,intent(in) :: fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
    2375              :  real(dp),intent(out) :: rhor(cplex*nfft,nspden)
    2376              :  real(dp),intent(inout) :: rhor_glob(cplex*product(ngfft(1:3)),nspden)
    2377              : 
    2378              : !Local variables-------------------------------
    2379              :  integer :: ispden,i1,i2,i3,me_fft,i3_local,my_fftbase,glob_fftbase
    2380              :  integer :: n1,n2,n3,ierr,nfft_tot
    2381              : ! *************************************************************************
    2382              : 
    2383            0 :  nfft_tot = product(ngfft(1:3)); me_fft = xmpi_comm_rank(comm_fft)
    2384            0 :  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
    2385              : 
    2386            0 :  if (nfft_tot == nfft) then
    2387              :    ! full rhor on each node, just do a copy
    2388            0 :    rhor = rhor_glob
    2389              :  else
    2390              :    ! if MPI-FFT we have to gather the global array on each node.
    2391            0 :    call xmpi_bcast(rhor_glob,master,comm_fft,ierr)
    2392            0 :    do ispden=1,nspden
    2393            0 :      do i3=1,n3
    2394            0 :        if (me_fft == fftn3_distrib(i3)) then
    2395            0 :          i3_local = ffti3_local(i3)
    2396            0 :          do i2=1,n2
    2397            0 :            my_fftbase =   cplex * ( (i2-1)*n1 + (i3_local-1)*n1*n2 )
    2398            0 :            glob_fftbase = cplex * ( (i2-1)*n1 + (i3-1)*n1*n2 )
    2399            0 :            do i1=1,cplex * n1
    2400            0 :              rhor(i1+my_fftbase,ispden) = rhor_glob(i1+glob_fftbase,ispden)
    2401              :            end do
    2402              :          end do
    2403              :        end if
    2404              :      end do
    2405              :    end do
    2406              :  end if
    2407              : 
    2408            0 : end subroutine distrib_datar
    2409              : !!***
    2410              : 
    2411              : !----------------------------------------------------------------------
    2412              : 
    2413              : !!****f* m_nctk/var_from_id
    2414              : !! NAME
    2415              : !!  var_from_id
    2416              : !!
    2417              : !! FUNCTION
    2418              : !!  Initialize a nctkvar_t object from the variable id
    2419              : !!
    2420              : !! INPUTS
    2421              : !!  ncid=NC file handle
    2422              : !!  varid=Variable ID
    2423              : !!
    2424              : !! OUTPUT
    2425              : !!  var<nctkvar_t>=Info on the variable.
    2426              : !!
    2427              : !! SOURCE
    2428              : 
    2429    665505850 : subroutine var_from_id(ncid, varid, var)
    2430              : 
    2431              : !Arguments ------------------------------------
    2432              :  integer, intent(in) :: ncid, varid
    2433              :  type(nctkvar_t), intent(out) :: var
    2434              : 
    2435              : !Local variables-------------------------------
    2436              : !scalars
    2437              :  integer :: ii, ncerr
    2438              :  !character(len=NF90_MAX_NAME) :: ncname
    2439              : ! *********************************************************************
    2440              : 
    2441              :  ! Get info about the variable.
    2442       324637 :  var%id = varid
    2443              :  ncerr = nf90_inquire_variable(ncid, var%id, &
    2444       324637 :    name=var%name, xtype=var%xtype, ndims=var%ndims, dimids=var%dimids, natts=var%natts)
    2445       324637 :  NCF_CHECK(ncerr)
    2446              : 
    2447              :  ! Get info about dimensions.
    2448       324637 :  if (var%ndims > 0) then
    2449       641401 :    do ii=1,var%ndims
    2450       415762 :      ncerr = nf90_inquire_dimension(ncid, var%dimids(ii), len=var%dimlens(ii), name=var%dimnames(ii))
    2451       641401 :      NCF_CHECK(ncerr)
    2452              :    end do
    2453              :  end if
    2454              : 
    2455              :  ! Get the number of attributes and their names.
    2456              :  !if (var%natts > 0) then
    2457              :  !   do ii=1,var%natts
    2458              :  !      ncerr = nf90_inq_attname(ncid, var%id, ii, var%attnames(ii))
    2459              :  !   end do
    2460              :  !end if
    2461              : 
    2462       324637 : end subroutine var_from_id
    2463              : !!***
    2464              : 
    2465              : !----------------------------------------------------------------------
    2466              : 
    2467              : !!****f* m_nctk/var_from_name
    2468              : !! NAME
    2469              : !!  var_from_name
    2470              : !!
    2471              : !! FUNCTION
    2472              : !!  Initialize a nctkvar_t object from the variable name
    2473              : !!
    2474              : !! INPUTS
    2475              : !!  ncid=NC file handle
    2476              : !!  name=Variable name
    2477              : !!
    2478              : !! OUTPUT
    2479              : !!  var<nctkvar_t>=Info on the variable.
    2480              : !!
    2481              : !! SOURCE
    2482              : 
    2483              : subroutine var_from_name(ncid, name, var)
    2484              : 
    2485              : !Arguments ------------------------------------
    2486              :  integer, intent(in) :: ncid
    2487              :  character(len=*),intent(in) :: name
    2488              :  type(nctkvar_t), intent(out) :: var
    2489              : 
    2490              : !Local variables-------------------------------
    2491              : !scalars
    2492              :  integer :: varid
    2493              : ! *********************************************************************
    2494              : 
    2495              :  varid = nctk_idname(ncid, name)
    2496              :  call var_from_id(ncid, varid, var)
    2497              : 
    2498              : end subroutine var_from_name
    2499              : !!***
    2500              : 
    2501              : !!****f* m_nctk/create_nc_file
    2502              : !! NAME
    2503              : !! create_nc_file
    2504              : !!
    2505              : !! FUNCTION
    2506              : !! Create an NetCDF file including a dimension one definition
    2507              : !!
    2508              : !! INPUTS
    2509              : !!
    2510              : !! OUTPUT
    2511              : !!
    2512              : !! TODO:
    2513              : !!  Remove
    2514              : !!
    2515              : !! SOURCE
    2516              : 
    2517         4162 : subroutine create_nc_file(filename, ncid)
    2518              : 
    2519              : !Arguments ------------------------------------
    2520              : !scalars
    2521              :  integer,intent(out) :: ncid
    2522              : !arrays
    2523              : character(len=*),intent(in) :: filename
    2524              : 
    2525              : !Local variables-------------------------------
    2526              : integer :: one_id, ncerr, cmode
    2527              : ! *************************************************************************
    2528              : 
    2529              :  ncid = 0
    2530              :  ! Create the NetCDF file
    2531              :  !ncerr = nf90_create(path=filename,cmode=NF90_CLOBBER,ncid=ncid)
    2532         4162 :  cmode = def_cmode_for_seq_create
    2533         4162 :  ncerr = nf90_create(path=filename, cmode=cmode, ncid=ncid)
    2534         4162 :  NCF_CHECK_MSG(ncerr, sjoin('Error while creating:', filename))
    2535         4162 :  ncerr=nf90_def_dim(ncid,'one',1,one_id)
    2536         4162 :  NCF_CHECK_MSG(ncerr,'nf90_def_dim')
    2537              : 
    2538         4162 :  end subroutine create_nc_file
    2539              : !!***
    2540              : 
    2541              : !!****f* m_nctk/write_var_netcdf
    2542              : !!
    2543              : !! NAME
    2544              : !! write_var_netcdf
    2545              : !!
    2546              : !! FUNCTION
    2547              : !! Write variable into a netcdf dataset
    2548              : !!
    2549              : !! TODO:
    2550              : !!  Remove!
    2551              : !!
    2552              : !! INPUTS
    2553              : !! arr_int
    2554              : !! arr_real
    2555              : !! marr
    2556              : !! narr
    2557              : !! typevar
    2558              : !! varname
    2559              : !!
    2560              : !! OUTPUT
    2561              : !!  (only writing)
    2562              : !!
    2563              : !! SOURCE
    2564              : 
    2565       421844 : subroutine write_var_netcdf(arr_int,arr_real,marr,narr,ncid,typevar,varname)
    2566              : 
    2567              : !Arguments ------------------------------------
    2568              : !scalars
    2569              :  integer,intent(in) :: narr,marr,ncid
    2570              :  character(len=*),intent(in) :: varname
    2571              :  character(len=3),intent(in) :: typevar
    2572              : !arrays
    2573              :  integer,intent(in) :: arr_int(marr)
    2574              :  real(dp),intent(in) :: arr_real(marr)
    2575              : 
    2576              : !Local variables-------------------------------
    2577              : !scalars
    2578              :  integer :: var_id,var_type,vardim_id,ncerr
    2579              :  !character(len=500) :: msg
    2580              : ! *************************************************************************
    2581              : 
    2582              :  !write(std_out,*)"about to write varname: ",trim(varname)
    2583              : 
    2584       421844 :  if (ncid>0) then
    2585              : !  ### Put the file in definition mode
    2586       290756 :    ncerr=nf90_redef(ncid)
    2587       290756 :    if (ncerr/=NF90_NOERR.and.ncerr/=NF90_EINDEFINE) then
    2588            0 :      NCF_CHECK_MSG(ncerr,'nf90_redef')
    2589              :    end if
    2590              : !  ### Define the dimensions
    2591       290756 :    if (narr==1)then
    2592       206846 :      ncerr=nf90_inq_dimid(ncid,'one',vardim_id)
    2593       206846 :      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
    2594              :    else
    2595        83910 :      ncerr=nf90_def_dim(ncid,trim(varname),narr,vardim_id)
    2596        83910 :      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
    2597              :    end if
    2598              : !  ### Define the variables
    2599       290756 :    if (typevar=='INT') then
    2600       176894 :      var_type=NF90_INT
    2601       113862 :    else if (typevar=='DPR') then
    2602       113862 :      var_type=NF90_DOUBLE
    2603              :    end if
    2604       290756 :    ncerr=nf90_def_var(ncid, trim(varname), var_type, vardim_id, var_id)
    2605       290756 :    NCF_CHECK_MSG(ncerr,'nf90_def_var')
    2606              : !  ### Put the file in data mode
    2607       290756 :    ncerr=nf90_enddef(ncid)
    2608       290756 :    if (ncerr/=NF90_NOERR.and.ncerr/=NF90_ENOTINDEFINE) then
    2609            0 :      NCF_CHECK_MSG(ncerr,'nf90_enddef')
    2610              :    end if
    2611              : !  ### Write variables into the dataset
    2612       290756 :    if (typevar=='INT') then
    2613       353788 :      ncerr=nf90_put_var(ncid,var_id,arr_int,start=(/1/),count=(/narr/))
    2614       113862 :    else if (typevar=='DPR') then
    2615       227724 :      ncerr=nf90_put_var(ncid,var_id,arr_real,start=(/1/),count=(/narr/))
    2616              :    end if
    2617       290756 :    NCF_CHECK_MSG(ncerr,'nf90_put_var')
    2618              :  end if
    2619              : 
    2620       421844 : end subroutine write_var_netcdf
    2621              : !!***
    2622              : 
    2623              : !!****f* ABINIT/write_eig
    2624              : !!
    2625              : !! NAME
    2626              : !! write_eig
    2627              : !!
    2628              : !! FUNCTION
    2629              : !! Write the eigenvalues band by band and k point by k point
    2630              : !! in a NetCDF file format
    2631              : !!
    2632              : !! INPUTS
    2633              : !! filname = Filename of the file where the history will be stored
    2634              : !!
    2635              : !! OUTPUT
    2636              : !!  (only writing)
    2637              : !!
    2638              : !! SOURCE
    2639              : 
    2640         2694 : subroutine write_eig(eigen,fermie,filename,kptns,mband,nband,nkpt,nsppol,&
    2641              : & extfpmd_eshift) ! Optional arguments
    2642              : 
    2643              : !Arguments ------------------------------------
    2644              : !scalars
    2645              :  character(len=fnlen),intent(in) :: filename
    2646              :  integer,intent(in) :: nkpt,nsppol,mband
    2647              :  real(dp),intent(in) :: fermie
    2648              :  real(dp),optional,intent(in) :: extfpmd_eshift
    2649              : !arrays
    2650              :  integer,intent(in) :: nband(nkpt*nsppol)
    2651              :  real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
    2652              :  real(dp),intent(in) :: kptns(3,nkpt)
    2653              : 
    2654              : !Local variables-------------------------------
    2655              : !scalars
    2656              :  integer :: ncerr,ncid,ii, cmode
    2657              :  integer :: xyz_id,nkpt_id,mband_id,nsppol_id
    2658              :  integer :: eig_id,fermie_id,kpt_id,nbk_id,nbk
    2659              :  integer :: extfpmd_eshift_id
    2660              :  integer :: ikpt,isppol,nband_k,band_index
    2661              :  real(dp):: convrt
    2662              : !arrays
    2663              :  integer :: dimEIG(3),dimKPT(2),dimNBK(2)
    2664              :  integer :: count2(2),start2(2)
    2665              :  integer :: count3(3),start3(3)
    2666              :  integer :: dim0(0)
    2667         5388 :  real(dp):: band(mband)
    2668              : ! *********************************************************************
    2669              : 
    2670         2694 :  convrt=1.0_dp
    2671              : 
    2672              : !1. Create netCDF file
    2673              :  !ncerr = nf90_create(path=trim(filename),cmode=NF90_CLOBBER, ncid=ncid)
    2674         2694 :  cmode = def_cmode_for_seq_create
    2675         2694 :  ncerr = nf90_create(path=trim(filename), cmode=cmode, ncid=ncid)
    2676         2694 :  NCF_CHECK_MSG(ncerr," create netcdf EIG file")
    2677              : 
    2678              : !2. Define dimensions
    2679         2694 :  ncerr = nf90_def_dim(ncid,"xyz",3,xyz_id)
    2680         2694 :  NCF_CHECK_MSG(ncerr," define dimension xyz")
    2681              : 
    2682         2694 :  ncerr = nf90_def_dim(ncid,"mband",mband,mband_id)
    2683         2694 :  NCF_CHECK_MSG(ncerr," define dimension mband")
    2684              : 
    2685         2694 :  ncerr = nf90_def_dim(ncid,"nkpt",nkpt,nkpt_id)
    2686         2694 :  NCF_CHECK_MSG(ncerr," define dimension nkpt")
    2687              : 
    2688         2694 :  ncerr = nf90_def_dim(ncid,"nsppol",nsppol,nsppol_id)
    2689         2694 :  NCF_CHECK_MSG(ncerr," define dimension nsppol")
    2690              : 
    2691              : !Dimensions for EIGENVALUES
    2692        10776 :  dimEIG = (/ mband_id, nkpt_id, nsppol_id /)
    2693              : !Dimensions for kpoint positions
    2694         8082 :  dimKPT = (/ xyz_id, nkpt_id /)
    2695              : !Dimensions for number kpoints per band and spin
    2696              : !dimNBK = (/ nkpt_id, nsppol_id /)
    2697         8082 :  dimNBK = (/ nkpt_id, nsppol_id /)
    2698              : 
    2699              : !3. Define variables
    2700              :  call ab_define_var(ncid,dim0,fermie_id,NF90_DOUBLE,&
    2701         2694 : & "fermie","Chemical potential","Hartree")
    2702              :  call ab_define_var(ncid, dimEIG, eig_id, NF90_DOUBLE,&
    2703              : & "Eigenvalues",&
    2704              : & "Values of eigenvalues",&
    2705         2694 : & "Hartree")
    2706              :  call ab_define_var(ncid, dimKPT, kpt_id, NF90_DOUBLE,"Kptns",&
    2707              : & "Positions of K-points in reciprocal space",&
    2708         2694 : & "Dimensionless")
    2709              :  call ab_define_var(ncid, dimNBK, nbk_id, NF90_INT,"NBandK",&
    2710              : & "Number of bands per kpoint and Spin",&
    2711         2694 : & "Dimensionless")
    2712         2694 :  if(present(extfpmd_eshift)) then
    2713              :     call ab_define_var(ncid,dim0,extfpmd_eshift_id,NF90_DOUBLE,&
    2714            0 : &    "extfpmd_eshift","Extended FPMD energy shift","Hartree")
    2715              :  end if
    2716              : 
    2717              : !4. End define mode
    2718         2694 :  ncerr = nf90_enddef(ncid)
    2719         2694 :  NCF_CHECK_MSG(ncerr," end define mode")
    2720              : 
    2721              : !5 Write kpoint positions
    2722        30647 :  do ikpt=1,nkpt
    2723        83859 :    start2 = (/ 1, ikpt /)
    2724        27953 :    count2 = (/ 3, 1 /)
    2725              :    ncerr = nf90_put_var(ncid, kpt_id,&
    2726              : &   kptns(1:3,ikpt),&
    2727              : &   start = start2,&
    2728        27953 : &   count = count2)
    2729        30647 :    NCF_CHECK_MSG(ncerr," write variable kptns")
    2730              :  end do
    2731              : 
    2732              : !6.1 Write chemical potential
    2733         2694 :  ncerr = nf90_put_var(ncid, fermie_id, fermie)
    2734         2694 :  NCF_CHECK_MSG(ncerr," write variable fermie")
    2735              : 
    2736              : !6.2 Write extfpmd shiftfactor
    2737         2694 :  if(present(extfpmd_eshift)) then
    2738            0 :    ncerr = nf90_put_var(ncid, extfpmd_eshift_id, extfpmd_eshift)
    2739            0 :    NCF_CHECK_MSG(ncerr," write variable extfpmd_eshift")
    2740              :  end if
    2741              : 
    2742              : !6.3 Write eigenvalues
    2743              :  band_index=0
    2744         5652 :  do isppol=1,nsppol
    2745        34650 :    do ikpt=1,nkpt
    2746        28998 :      nband_k=nband(ikpt+(isppol-1)*nkpt)
    2747       115992 :      start3 = (/ 1, ikpt, isppol /)
    2748       115992 :      count3 = (/ mband, 1, 1 /)
    2749       291590 :      band(:)=zero
    2750       291514 :      do ii=1,nband_k
    2751       291514 :        band(ii)=eigen(band_index+ii)
    2752              :      end do
    2753              :      ncerr = nf90_put_var(ncid, eig_id,&
    2754              : &     band,&
    2755              : &     start = start3,&
    2756        28998 : &     count = count3)
    2757        28998 :      NCF_CHECK_MSG(ncerr," write variable band")
    2758              : 
    2759        31956 :      band_index=band_index+nband_k
    2760              :    end do
    2761              :  end do
    2762              : 
    2763              : !6.4 Write Number of bands per kpoint and Spin
    2764              : 
    2765         5652 :  do isppol=1,nsppol
    2766        34650 :    do ikpt=1,nkpt
    2767        86994 :      start2 = (/ ikpt, 1 /)
    2768        28998 :      count2 = (/ 1, 1 /)
    2769        28998 :      nbk=nband(ikpt+(isppol-1)*nkpt)
    2770              :      ncerr = nf90_put_var(ncid, nbk_id,&
    2771              : &     nbk,&
    2772        28998 : &     start = start2)
    2773        31956 :      NCF_CHECK_MSG(ncerr," write variable nband")
    2774              :    end do
    2775              :  end do
    2776              : 
    2777              : !7 Close file
    2778              : 
    2779         2694 :  ncerr = nf90_close(ncid)
    2780         2694 :  NCF_CHECK_MSG(ncerr," close netcdf EIG file")
    2781              : 
    2782         2694 : end subroutine write_eig
    2783              : !!***
    2784              : 
    2785            0 : END MODULE m_nctk
    2786              : !!***
        

Generated by: LCOV version 2.3-1