diff --git a/cmake/source_list_KKRhost.txt b/cmake/source_list_KKRhost.txt
index 4b6cea7bf9e4a36444cdb5fb3ad65f03e2b6e265..c693b7e31c3cd9b38606b94b9412ba839ebbf6bf 100644
--- a/cmake/source_list_KKRhost.txt
+++ b/cmake/source_list_KKRhost.txt
@@ -110,6 +110,8 @@ add_library(lib_common STATIC
     source/common/mathtools.f90
     source/common/bfield.f90
     source/common/torque.f90
+    source/external/NPY-for-Fortran/src/npy.F90
+    source/external/NPY-for-Fortran/src/endian_swap.f90
 )
 # disable cmake auto add of 'lib' prefix to .so file
 SET_TARGET_PROPERTIES(lib_common PROPERTIES PREFIX "")
diff --git a/source/KKRhost/write_gflle_npy.f90 b/source/KKRhost/write_gflle_npy.f90
index e12d0bd6a49810ceeece86a84a1ddcb94d0080a9..17afbec6ac7a3c9f8d47605ff3292efc3efaa5f1 100644
--- a/source/KKRhost/write_gflle_npy.f90
+++ b/source/KKRhost/write_gflle_npy.f90
@@ -1,962 +1,3 @@
-module m_npy
-   implicit none
-
-   integer(4), parameter               :: p_un = 23
-   character, parameter                :: magic_num = achar(147) ! x93
-   character, parameter                :: major = achar(2)   !major *.npy version
-   character, parameter                :: minor = achar(0)   !minor *.npy version
-   character(len=*), parameter         :: zip_flag = "-q0"
-   character(len=*), parameter         :: magic_str = "NUMPY"
-
-   interface save_npy
-      module procedure write_int64_vec, write_int64_mtx, &
-         write_int32_vec, write_int32_mtx, write_int32_3d, &
-         write_int16_vec, write_int16_mtx, &
-         write_int8_vec, write_int8_mtx, write_int8_3d, &
-         write_dbl_vec, write_dbl_mtx, &
-         write_sng_vec, write_sng_mtx, &
-         write_cmplx_sgn_vec, write_cmplx_sgn_mtx, &
-         write_cmplx_dbl_vec, write_cmplx_dbl_mtx, &
-         write_sng_3dT, write_dbl_3dT, &
-         write_sng_4dT, write_dbl_4dT, &
-         write_dbl_5dT, &
-         write_cmplx_dbl_3dT, &
-         write_cmplx_dbl_4dT, &
-         write_cmplx_dbl_5dT, &
-         write_cmplx_dbl_6dT
-
-   end interface save_npy
-   interface add_npz
-      module procedure addrpl_int8_vec, addrpl_int8_mtx, &
-         addrpl_int16_vec, addrpl_int16_mtx, &
-         addrpl_int32_vec, addrpl_int32_mtx, &
-         addrpl_int64_vec, addrpl_int64_mtx, &
-         addrpl_sng_vec, addrpl_sng_mtx, &
-         addrpl_dbl_vec, addrpl_dbl_mtx, &
-         addrpl_cmplx_dbl_vec, addrpl_cmplx_dbl_mtx, &
-         addrpl_cmplx_sng_vec, addrpl_cmplx_sng_mtx
-   end interface add_npz
-
-contains
-   subroutine run_sys(cmd, stat)
-      implicit none
-      character(len=*), intent(in)     :: cmd
-      integer(4), intent(out)          :: stat
-
-      call execute_command_line(cmd, wait=.True., exitstat=stat)
-   end subroutine run_sys
-
-   subroutine addrpl_cmplx_sng_vec(zipfile, var_name, vec)
-      implicit none
-      complex(4), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_cmplx_sng_vec
-
-   subroutine addrpl_cmplx_sng_mtx(zipfile, var_name, mtx)
-      implicit none
-      complex(4), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_cmplx_sng_mtx
-
-   subroutine addrpl_cmplx_dbl_vec(zipfile, var_name, vec)
-      implicit none
-      complex(8), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_cmplx_dbl_vec
-
-   subroutine addrpl_cmplx_dbl_mtx(zipfile, var_name, mtx)
-      implicit none
-      complex(8), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_cmplx_dbl_mtx
-
-   subroutine addrpl_dbl_vec(zipfile, var_name, vec)
-      implicit none
-      real(8), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_dbl_vec
-
-   subroutine addrpl_dbl_mtx(zipfile, var_name, mtx)
-      implicit none
-      real(8), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_dbl_mtx
-
-   subroutine addrpl_sng_vec(zipfile, var_name, vec)
-      implicit none
-      real(4), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_sng_vec
-
-   subroutine addrpl_sng_mtx(zipfile, var_name, mtx)
-      implicit none
-      real(4), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_sng_mtx
-
-   subroutine addrpl_int8_vec(zipfile, var_name, vec)
-      implicit none
-      integer(1), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int8_vec
-
-   subroutine addrpl_int8_mtx(zipfile, var_name, mtx)
-      implicit none
-      integer(1), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int8_mtx
-
-   subroutine addrpl_int16_vec(zipfile, var_name, vec)
-      implicit none
-      integer(2), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int16_vec
-
-   subroutine addrpl_int16_mtx(zipfile, var_name, mtx)
-      implicit none
-      integer(2), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int16_mtx
-
-   subroutine addrpl_int32_vec(zipfile, var_name, vec)
-      implicit none
-      integer(4), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int32_vec
-
-   subroutine addrpl_int32_mtx(zipfile, var_name, mtx)
-      implicit none
-      integer(4), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int32_mtx
-
-   subroutine addrpl_int64_vec(zipfile, var_name, vec)
-      implicit none
-      integer(8), intent(in)           :: vec(:)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, vec)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int64_vec
-
-   subroutine addrpl_int64_mtx(zipfile, var_name, mtx)
-      implicit none
-      integer(8), intent(in)           :: mtx(:, :)
-      character(len=*), intent(in)     :: zipfile, var_name
-      character(len=:), allocatable    :: npy_name
-      integer(4)                       :: succ
-
-      npy_name = var_name//".npy"
-
-      call save_npy(npy_name, mtx)
-      ! just store and be quite while zipping
-      call run_sys("zip "//zip_flag//" "//zipfile &
-                   //" "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute zip command"
-      endif
-
-      call run_sys("rm "//npy_name, succ)
-      if (succ /= 0) then
-         write (*, *) "Can't execute rm command"
-      endif
-   end subroutine addrpl_int64_mtx
-
-   Subroutine write_cmplx_sgn_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(4), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<c8"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-      write (p_un) header_len
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_cmplx_sgn_mtx
-
-   Subroutine write_cmplx_sgn_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(4), intent(in)           :: vec(:)
-      character(len=*), parameter      :: var_type = "<c8"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_cmplx_sgn_vec
-
-   Subroutine write_cmplx_dbl_6dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: tensor(:, :, :, :, :, :)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_6dT
-
-   Subroutine write_cmplx_dbl_5dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: tensor(:, :, :, :, :)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_5dT
-
-   Subroutine write_cmplx_dbl_4dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: tensor(:, :, :, :)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_4dT
-
-   Subroutine write_cmplx_dbl_3dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: tensor(:, :, :)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_3dT
-
-   Subroutine write_cmplx_dbl_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_mtx
-
-   Subroutine write_cmplx_dbl_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      complex(8), intent(in)           :: vec(:)
-      character(len=*), parameter      :: var_type = "<c16"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_cmplx_dbl_vec
-
-   Subroutine write_sng_3dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(4), intent(in)              :: tensor(:, :, :)
-      character(len=*), parameter      :: var_type = "<f4"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_sng_3dT
-
-   Subroutine write_sng_4dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(4), intent(in)              :: tensor(:, :, :, :)
-      character(len=*), parameter      :: var_type = "<f4"
-      integer(4)                       :: header_len
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_sng_4dT
-
-   Subroutine write_sng_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(4), intent(in)              :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<f4"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_sng_mtx
-
-   Subroutine write_sng_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(4), intent(in)              :: vec(:)
-      character(len=*), parameter      :: var_type = "<f4"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_sng_vec
-
-   Subroutine write_dbl_3dT(filename, tensor)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(8), intent(in)              :: tensor(:, :, :)
-      character(len=*), parameter      :: var_type = "<f8"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor))
-      write (p_un) tensor
-      close (unit=p_un)
-   End Subroutine write_dbl_3dT
-
-   Subroutine write_dbl_4dT(filename, tensor4)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(8), intent(in)              :: tensor4(:, :, :, :)
-      character(len=*), parameter      :: var_type = "<f8"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor4)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor4))
-      write (p_un) tensor4
-      close (unit=p_un)
-   End Subroutine write_dbl_4dT
-
-   Subroutine write_dbl_5dT(filename, tensor5)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(8), intent(in)              :: tensor5(:, :, :, :, :)
-      character(len=*), parameter      :: var_type = "<f8"
-      integer(4)                       :: header_len, i, j, k
-
-      header_len = len(dict_str(var_type, shape(tensor5)))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, shape(tensor5))
-      write (p_un) tensor5
-      close (unit=p_un)
-   End Subroutine write_dbl_5dT
-
-   Subroutine write_dbl_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(8), intent(in)              :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<f8"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_dbl_mtx
-
-   Subroutine write_dbl_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      real(8), intent(in)              :: vec(:)
-      character(len=*), parameter      :: var_type = "<f8"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_dbl_vec
-
-   Subroutine write_int64_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(8), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<i8"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_int64_mtx
-
-   Subroutine write_int64_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(8), intent(in)           :: vec(:)
-      character(len=*), parameter      :: var_type = "<i8"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_int64_vec
-
-   Subroutine write_int32_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(4), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<i4"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_int32_mtx
-
-   Subroutine write_int32_3d(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(4), intent(in)           :: mtx(:,:,:)
-      character(len=*), parameter      :: var_type = "<i4"
-      integer(4)                       :: header_len, s_mtx(3), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_int32_3d
-
-   Subroutine write_int32_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(4), intent(in)           :: vec(:)
-      character(len=*), parameter      :: var_type = "<i4"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_int32_vec
-
-   Subroutine write_int16_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(2), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<i2"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-      close (unit=p_un)
-   End Subroutine write_int16_mtx
-
-   Subroutine write_int16_vec(filename, vec)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(2), intent(in)           :: vec(:)
-      character(len=*), parameter      :: var_type = "<i2"
-      integer(4)                       :: header_len, s_vec(1), i
-
-      s_vec = shape(vec)
-      header_len = len(dict_str(var_type, s_vec))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_vec)
-
-      write (p_un) vec
-
-      close (unit=p_un)
-   End Subroutine write_int16_vec
-
-   Subroutine write_int8_mtx(filename, mtx)
-      Implicit None
-      character(len=*), intent(in)     :: filename
-      integer(1), intent(in)           :: mtx(:, :)
-      character(len=*), parameter      :: var_type = "<i1"
-      integer(4)                       :: header_len, s_mtx(2), i, j
-
-      s_mtx = shape(mtx)
-      header_len = len(dict_str(var_type, s_mtx))
-
-      open (unit=p_un, file=filename, form="unformatted", &
-            access="stream")
-      write (p_un) magic_num, magic_str, major, minor
-
-      write (p_un) header_len
-
-      write (p_un) dict_str(var_type, s_mtx)
-
-      write (p_un) mtx
-
-
 module mod_write_gflle 
 
   implicit none
diff --git a/source/external/NPY-for-Fortran/.gitignore b/source/external/NPY-for-Fortran/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..442d1b0776e7b36fdc6a7fc967f6219d7d46d179
--- /dev/null
+++ b/source/external/NPY-for-Fortran/.gitignore
@@ -0,0 +1,6 @@
+CMakeCache.txt
+CMakeFiles/*
+Makefile
+cmake_install.cmake
+*.x
+*.mod
diff --git a/source/external/NPY-for-Fortran/CMakeLists.txt b/source/external/NPY-for-Fortran/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..33389d1dcedbb7c035ef0f8ec6f2f3015390e232
--- /dev/null
+++ b/source/external/NPY-for-Fortran/CMakeLists.txt
@@ -0,0 +1,35 @@
+cmake_minimum_required(VERSION 2.8)
+project(STB)
+
+file(GLOB_RECURSE sources  src/*.f90
+                           src/*.F90
+                           src/*.h )
+
+add_executable(npy.x ${sources})
+
+enable_language(Fortran)
+set(CMAKE_Fortran_COMPILER ifort)
+set(CMAKE_Fortran_COMPILER_ID "Intel")
+
+
+if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
+    set(dialect "-ffree-form -cpp -std=gnu -fimplicit-none")
+	#set(bounds "-fbounds-check")
+endif()
+if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel")
+	#set(bounds "-check bounds")
+endif()
+if(CMAKE_Fortran_COMPILER_ID MATCHES "PGI")
+    set(dialect "-Mfreeform -Mdclchk -Mstandard -Mallocatable=03")
+	#set(bounds "-C")
+endif()
+
+set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${bounds}")
+set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${dialect}")
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
+MESSAGE( STATUS "cmake_module_path:    " ${CMAKE_MODULE_PATH})
+
+#
+# Compile.
+#
diff --git a/source/external/NPY-for-Fortran/LICENSE b/source/external/NPY-for-Fortran/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..e5165484155e63e0730ed35778e2202aed30a23a
--- /dev/null
+++ b/source/external/NPY-for-Fortran/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2017 Matthias Redies
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/source/external/NPY-for-Fortran/README.md b/source/external/NPY-for-Fortran/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..9d1c39dbb72cc3750a7ca811fa1c3a9c150e8b6b
--- /dev/null
+++ b/source/external/NPY-for-Fortran/README.md
@@ -0,0 +1,73 @@
+# NPY for Fortran
+This Fortran module allows to save numerical Fortran arrays in Numpy's .npy or .npz format. Currently supported are:
+```fortran
+1. integer(1), integer(2), integer(4), integer(8)
+2. real(4), real(8)
+3. complex(4), complex(8)
+```
+### *.npy files
+Saving an array into a .npy-file is simply done by calling:
+```fortran
+call save_npy("filename.npy", array)
+```
+
+
+### *.npz files
+In order to save .npz-files the commandline tool 'zip' has to be installed. By calling 
+```fortran
+call add_npz("example.npz", "temperature", data_array)
+```
+one creates an .npz-file containing data_array, with the name "temperature". If example.npz already exists the field "temperature" is added to it. If the field temperature already exsits in example.npz it will be overwritten.
+
+
+Reading .npy and .npz files isn't currently supported. (Maybe someone can give me ideas on dynamic typing in Fortran...)
+
+### Compiling using ifort
+
+The code uses the somewhat out-dated 
+```fortran
+succ = system(...)
+```
+command for which ifort needs the the IFPORT library:
+
+```fortran
+#ifdef INTEL_COMPILER_USED
+    USE IFPORT
+#endif
+```
+
+Therefore, Intel users need to add the flag:
+```
+-DINTEL_COMPILER_USED
+```
+
+### Compiling using gfortran
+Since the function 'system' is not standard Fortran one cannot use flags such as
+```
+-std=f2008
+```
+but instead can (not must) use the equivalent
+```
+-std=gnu
+```
+### Control Flags
+
+In the source file the parameter
+
+```fortran
+    character(len=*), parameter         :: zip_flag  = "-q0"  
+```
+
+can be used to control the compression. "-q0" tells the zip tool to have no output and no compression. 
+
+
+
+### Endianess
+
+Using the flag
+
+```fortran
+    logical, parameter                  :: use_big_endian = .False.
+```
+
+the output [endianess](https://en.wikipedia.org/wiki/Least_significant_bit) can be set. This works independent of the compiler. Little endian (.False.) is a reasonable default for most personal computers.
\ No newline at end of file
diff --git a/source/external/NPY-for-Fortran/_config.yml b/source/external/NPY-for-Fortran/_config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2f7efbeab578c8042531ea7908ee8ffd7589fe46
--- /dev/null
+++ b/source/external/NPY-for-Fortran/_config.yml
@@ -0,0 +1 @@
+theme: jekyll-theme-minimal
\ No newline at end of file
diff --git a/source/external/NPY-for-Fortran/src/endian_swap.f90 b/source/external/NPY-for-Fortran/src/endian_swap.f90
new file mode 100644
index 0000000000000000000000000000000000000000..4d2cf4f8cc9e2a825c419cb318915007ad253fde
--- /dev/null
+++ b/source/external/NPY-for-Fortran/src/endian_swap.f90
@@ -0,0 +1,183 @@
+module endian_swap
+    implicit none
+
+    PRIVATE
+    PUBLIC :: Big_Endian
+    PUBLIC :: Swap_Endian
+
+    INTERFACE Swap_Endian
+        module procedure SWAP_I1
+        module procedure SWAP_I2 
+        module procedure SWAP_I4
+        module procedure SWAP_I8
+        module procedure SWAP_F4
+        module procedure SWAP_F8 
+        module procedure SWAP_F16
+        module procedure SWAP_C4
+        module procedure SWAP_C8
+    END INTERFACE Swap_Endian
+
+
+    CONTAINS
+
+        FUNCTION Big_Endian()
+
+            LOGICAL :: Big_Endian
+
+            Big_Endian = ichar(transfer(1,'a')) == 0
+
+        END FUNCTION Big_Endian
+
+        function SWAP_I4(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz = 4
+            integer(b_sz), intent(in) :: input 
+            integer(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_I4 
+
+        function SWAP_I2(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  2
+            integer(b_sz), intent(in) :: input 
+            integer(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_I2
+        
+        function SWAP_I1(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  1
+            integer(b_sz), intent(in) :: input 
+            integer(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_I1
+        
+        function SWAP_I8(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  8
+            integer(b_sz), intent(in) :: input 
+            integer(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_I8
+
+
+        function SWAP_F4(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  4
+            real(b_sz), intent(in) :: input 
+            real(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_F4
+        
+        function SWAP_F8(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  8
+            real(b_sz), intent(in) :: input 
+            real(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_F8 
+
+        function SWAP_F16(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz =  16
+            real(b_sz), intent(in) :: input 
+            real(b_sz)             :: output 
+
+            integer(1) :: byte_arr(b_sz), byte_arr_tmp(b_sz)
+            integer(1) :: i
+
+            byte_arr_tmp =  transfer(input, byte_arr_tmp)
+
+            do i = 1,b_sz
+                byte_arr(i) =  byte_arr_tmp(1 + b_sz - i)
+            enddo
+
+            output =  transfer(byte_arr, output)
+        end function SWAP_F16
+        
+        function SWAP_C8(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz = 8
+            complex(b_sz), intent(in) :: input
+            complex(b_sz)             :: output
+            real(b_sz)                :: re, im
+
+            re = Swap_Endian(real(input))
+            im = Swap_Endian(aimag(input))
+
+            output = cmplx(re, im)
+        end function swap_c8
+
+        function SWAP_C4(input) result(output)
+            implicit none
+            integer(4), parameter  :: b_sz = 4
+            complex(b_sz), intent(in) :: input
+            complex(b_sz)             :: output
+            real(b_sz)                :: re, im
+
+            re = Swap_Endian(real(input))
+            im = Swap_Endian(aimag(input))
+
+            output = cmplx(re, im)
+        end function swap_c4
+END module endian_swap
diff --git a/source/external/NPY-for-Fortran/src/example.f90 b/source/external/NPY-for-Fortran/src/example.f90
new file mode 100644
index 0000000000000000000000000000000000000000..bc6aad92f7d9197f65a96e5891f745d5f439239f
--- /dev/null
+++ b/source/external/NPY-for-Fortran/src/example.f90
@@ -0,0 +1,32 @@
+program main
+    use m_npy
+    use endian_swap
+
+    complex(8)       :: a(10,20), b(10), c(2,3,4)
+    integer(4)       :: i, j, k
+    real(4)  :: test1
+    real(8)  :: test2
+
+    do i =  1,size(a,1)
+        do j =  1,size(a,2)
+            a(i,j) = i * j
+        enddo
+    enddo
+
+    do i =  1,size(b,1)
+        b(i) = 2 *  i
+    enddo
+
+    do i = 1,size(c,1)
+        do j = 1,size(c,2)
+            do k =1,size(c,3)
+                c(i,j,k) = cmplx(1, 2)
+            enddo
+        enddo
+    enddo
+    
+    c(1,1,1) = cmplx(3,4)
+    c(2,3,4) = cmplx(0,1)
+    c(2,2,2) = cmplx(1,0)
+    call save_npy("c.npy", c)
+end program main
diff --git a/source/external/NPY-for-Fortran/src/npy.F90 b/source/external/NPY-for-Fortran/src/npy.F90
new file mode 100644
index 0000000000000000000000000000000000000000..f51f992e6708e17770ef9fa11fe462bd3786aa13
--- /dev/null
+++ b/source/external/NPY-for-Fortran/src/npy.F90
@@ -0,0 +1,1184 @@
+module  m_npy
+    use endian_swap
+    implicit none
+
+    integer(4), parameter               :: p_un      = 23
+    character, parameter                :: magic_num = achar(147) ! x93
+    character, parameter                :: major     = achar(2)   !major *.npy version 
+    character, parameter                :: minor     = achar(0)   !minor *.npy version
+    logical, parameter                  :: use_big_endian = .False.
+    character(len=*), parameter         :: zip_flag  = "-q0"  
+    character(len=*), parameter         :: magic_str = "NUMPY"
+
+    interface save_npy
+        module procedure write_int64_vec,     write_int64_mtx, &
+                         write_int32_vec,     write_int32_mtx, &
+                         write_int16_vec,     write_int16_mtx, &
+                         write_int8_vec,      write_int8_mtx, &
+                         write_dbl_vec,       write_dbl_mtx,   &
+                         write_sng_vec,       write_sng_mtx,   &
+                         write_cmplx_sgn_vec, write_cmplx_sgn_mtx, &
+                         write_cmplx_dbl_vec, write_cmplx_dbl_mtx, &
+                         write_sng_3dT,       write_dbl_3dT, &
+                         write_dbl_4dT,&
+                         write_dbl_5dT,&
+                         write_cmplx_dbl_3dT,&
+                         write_cmplx_dbl_4dT,&
+                         write_cmplx_dbl_5dT,&
+                         write_cmplx_dbl_6dT
+
+
+    end interface save_npy
+    interface add_npz
+        module procedure addrpl_int8_vec,      addrpl_int8_mtx,      &
+                         addrpl_int16_vec,     addrpl_int16_mtx,     &
+                         addrpl_int32_vec,     addrpl_int32_mtx,     &
+                         addrpl_int64_vec,     addrpl_int64_mtx,     &
+                         addrpl_sng_vec,       addrpl_sng_mtx,       &
+                         addrpl_dbl_vec,       addrpl_dbl_mtx,       &
+                         addrpl_cmplx_dbl_vec, addrpl_cmplx_dbl_mtx, &
+                         addrpl_cmplx_sng_vec, addrpl_cmplx_sng_mtx
+    end interface add_npz
+
+contains
+    subroutine run_sys(cmd, stat)
+        implicit none
+        character(len=*), intent(in)     :: cmd
+        integer(4), intent(out)          :: stat
+        
+        call execute_command_line(cmd, wait=.True., exitstat=stat)
+    end subroutine run_sys
+    
+    subroutine addrpl_cmplx_sng_vec(zipfile, var_name, vec)
+        implicit none
+        complex(4), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_cmplx_sng_vec
+    
+    subroutine addrpl_cmplx_sng_mtx(zipfile, var_name, mtx)
+        implicit none
+        complex(4), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_cmplx_sng_mtx
+    
+    subroutine addrpl_cmplx_dbl_vec(zipfile, var_name, vec)
+        implicit none
+        complex(8), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_cmplx_dbl_vec
+    
+    subroutine addrpl_cmplx_dbl_mtx(zipfile, var_name, mtx)
+        implicit none
+        complex(8), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_cmplx_dbl_mtx
+    
+    subroutine addrpl_dbl_vec(zipfile, var_name, vec)
+        implicit none
+        real(8), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_dbl_vec
+    
+    subroutine addrpl_dbl_mtx(zipfile, var_name, mtx)
+        implicit none
+        real(8), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_dbl_mtx
+    
+    subroutine addrpl_sng_vec(zipfile, var_name, vec)
+        implicit none
+        real(4), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_sng_vec
+    
+    subroutine addrpl_sng_mtx(zipfile, var_name, mtx)
+        implicit none
+        real(4), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_sng_mtx
+    
+    subroutine addrpl_int8_vec(zipfile, var_name, vec)
+        implicit none
+        integer(1), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int8_vec
+    
+    subroutine addrpl_int8_mtx(zipfile, var_name, mtx)
+        implicit none
+        integer(1), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int8_mtx
+    
+    subroutine addrpl_int16_vec(zipfile, var_name, vec)
+        implicit none
+        integer(2), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int16_vec
+    
+    subroutine addrpl_int16_mtx(zipfile, var_name, mtx)
+        implicit none
+        integer(2), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int16_mtx
+    
+    subroutine addrpl_int32_vec(zipfile, var_name, vec)
+        implicit none
+        integer(4), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int32_vec
+    
+    subroutine addrpl_int32_mtx(zipfile, var_name, mtx)
+        implicit none
+        integer(4), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int32_mtx
+    
+    subroutine addrpl_int64_vec(zipfile, var_name, vec)
+        implicit none
+        integer(8), intent(in)           :: vec(:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, vec)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int64_vec
+    
+    subroutine addrpl_int64_mtx(zipfile, var_name, mtx)
+        implicit none
+        integer(8), intent(in)           :: mtx(:,:)
+        character(len=*), intent(in)     :: zipfile, var_name
+        character(len=:), allocatable    :: npy_name
+        integer(4)                       :: succ
+
+        npy_name =  var_name // ".npy"
+
+        call save_npy(npy_name, mtx)
+        ! just store and be quite while zipping
+        call run_sys("zip " // zip_flag // " " // zipfile &
+                                 // " " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute zip command"
+        endif
+
+        call run_sys("rm " // npy_name, succ)
+        if(succ /=  0) then
+            write (*,*) "Can't execute rm command"
+        endif
+    end subroutine addrpl_int64_mtx
+
+    Subroutine write_cmplx_sgn_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(4), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<c8"
+        integer(4)                       :: header_len, s_mtx(2), i, j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_mtx)
+
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+
+        close(unit=p_un)
+    End Subroutine write_cmplx_sgn_mtx
+    
+    Subroutine  write_cmplx_sgn_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(4), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<c8"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_cmplx_sgn_vec
+
+
+    Subroutine  write_cmplx_dbl_6dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: tensor(:,:,:,:,:,:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "6D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_6dT
+
+    Subroutine  write_cmplx_dbl_5dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: tensor(:,:,:,:,:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "5D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_5dT
+
+     Subroutine  write_cmplx_dbl_4dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: tensor(:,:,:,:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "4D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_4dT
+
+    Subroutine  write_cmplx_dbl_3dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: tensor(:,:,:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "3D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_3dT
+
+    Subroutine  write_cmplx_dbl_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, s_mtx(2),i,j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_mtx
+
+
+    
+    Subroutine  write_cmplx_dbl_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        complex(8), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<c16"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_cmplx_dbl_vec
+
+    Subroutine  write_sng_3dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(4), intent(in)              :: tensor(:,:,:)
+        character(len=*), parameter      :: var_type =  "<f4"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "3D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_sng_3dT
+
+
+    Subroutine  write_sng_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(4), intent(in)              :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<f4"
+        integer(4)                       :: header_len, s_mtx(2), i, j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+    
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_sng_mtx
+    
+    Subroutine  write_sng_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(4), intent(in)              :: vec(:)
+        character(len=*), parameter      :: var_type =  "<f4"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_sng_vec
+    
+    Subroutine  write_dbl_3dT(filename, tensor)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(8), intent(in)              :: tensor(:,:,:)
+        character(len=*), parameter      :: var_type =  "<f8"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "3D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor))
+        write (p_un) tensor
+        close(unit=p_un)
+    End Subroutine write_dbl_3dT
+
+
+    Subroutine  write_dbl_4dT(filename, tensor4)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(8), intent(in)              :: tensor4(:,:,:,:)
+        character(len=*), parameter      :: var_type =  "<f8"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor4)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "4D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor4))
+        write (p_un) tensor4
+        close(unit=p_un)
+    End Subroutine write_dbl_4dT
+
+
+    Subroutine  write_dbl_5dT(filename, tensor5)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(8), intent(in)              :: tensor5(:,:,:,:,:)
+        character(len=*), parameter      :: var_type =  "<f8"
+        integer(4)                       :: header_len, i,j, k
+
+        header_len =  len(dict_str(var_type, shape(tensor5)))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (*,*) "5D tensors not implemented on BigEndian"
+            write (*,*) "write in issue if you need it"
+            stop 7
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, shape(tensor5))
+        write (p_un) tensor5
+        close(unit=p_un)
+    End Subroutine write_dbl_5dT
+
+    Subroutine  write_dbl_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(8), intent(in)              :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<f8"
+        integer(4)                       :: header_len, s_mtx(2), i,j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_dbl_mtx
+    
+    Subroutine  write_dbl_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        real(8), intent(in)              :: vec(:)
+        character(len=*), parameter      :: var_type =  "<f8"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_dbl_vec
+    
+    Subroutine  write_int64_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(8), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<i8"
+        integer(4)                       :: header_len, s_mtx(2), i, j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int64_mtx
+    
+    Subroutine  write_int64_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(8), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<i8"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int64_vec
+
+
+    Subroutine  write_int32_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(4), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<i4"
+        integer(4)                       :: header_len, s_mtx(2), i, j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int32_mtx
+    
+    Subroutine  write_int32_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(4), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<i4"
+        integer(4)                       :: header_len, s_vec(1), i 
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int32_vec
+    
+    Subroutine  write_int16_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(2), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<i2"
+        integer(4)                       :: header_len, s_mtx(2), i, j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int16_mtx
+    
+    Subroutine  write_int16_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(2), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<i2"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int16_vec
+    
+    Subroutine  write_int8_mtx(filename, mtx)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(1), intent(in)           :: mtx(:,:)
+        character(len=*), parameter      :: var_type =  "<i1"
+        integer(4)                       :: header_len, s_mtx(2), i,j
+
+        s_mtx = shape(mtx)
+        header_len =  len(dict_str(var_type, s_mtx))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_mtx)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) mtx
+        else
+            do j = 1,size(mtx,2)
+                do i =  1,size(mtx,1)
+                    write (p_un) Swap_Endian(mtx(i,j))
+                enddo
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int8_mtx
+    
+    Subroutine  write_int8_vec(filename, vec)
+        Implicit None
+        character(len=*), intent(in)     :: filename
+        integer(1), intent(in)           :: vec(:)
+        character(len=*), parameter      :: var_type =  "<i1"
+        integer(4)                       :: header_len, s_vec(1), i
+
+        s_vec = shape(vec)
+        header_len =  len(dict_str(var_type, s_vec))
+        
+        open(unit=p_un, file=filename, form="unformatted",&
+             access="stream")
+        write (p_un) magic_num, magic_str, major, minor
+        if(Big_Endian()) then
+            write (p_un) Swap_Endian(header_len)
+        else
+            write (p_un) header_len
+        endif
+        write (p_un) dict_str(var_type, s_vec)
+        
+        if(use_big_endian .eqv. Big_Endian()) then  
+            write (p_un) vec
+        else
+            do i =  1,size(vec)
+                write (p_un) Swap_Endian(vec(i))
+            enddo
+        endif
+        close(unit=p_un)
+    End Subroutine write_int8_vec
+
+    function dict_str(var_type, var_shape) result(str)
+        implicit none
+        character(len=*), intent(in)   :: var_type 
+        integer(4), intent(in)         :: var_shape(:)
+        character(len=:), allocatable  :: str
+        integer(4)                     :: cnt
+
+        cnt =  len("{'descr': '")
+        cnt =  cnt + len(var_type)
+        cnt =  cnt +  len("', 'fortran_order': True, 'shape': (")
+        cnt =  cnt +  len(shape_str(var_shape))
+        cnt =  cnt +  len(",), }")
+        do while(mod(cnt +  10, 16) /= 0)
+            cnt =  cnt +  1
+        enddo
+
+        allocate(character(cnt) :: str)
+
+        str = "{'descr': '" // var_type // &
+              "', 'fortran_order': True, 'shape': (" // &
+              shape_str(var_shape) //  "), }"
+
+        do while(mod(len(str) + 11, 16) /= 0)
+            str = str // " "
+        enddo
+
+        str = str // achar(10)
+
+    end function dict_str
+
+    function shape_str(var_shape) result(fin_str)
+        implicit none
+        integer(4), intent(in)        :: var_shape(:)
+        character(len=:), allocatable :: str, small_str, fin_str
+        integer(4)                    :: i, length, start, halt
+
+        length = 14 * size(var_shape)
+        allocate(character(length) :: str)
+        allocate(character(14)     :: small_str)
+        str =  " "
+        
+        do i =  1, size(var_shape)
+            start = (i-1) * length + 1
+            halt  = i     * length +  1
+            write (small_str, "(I13,A)") var_shape(i), ","
+            str =  trim(str) // adjustl(small_str)
+        enddo
+        
+        fin_str =  trim(str)
+    end function shape_str 
+end module  m_npy