diff --git a/cmake/source_list_KKRhost.txt b/cmake/source_list_KKRhost.txt
index c80ef24f89fd07e03536c80acb960abe0ec9ec9f..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 "")
@@ -389,6 +391,7 @@ add_library(lib_kkrhost STATIC
     source/KKRhost/wrmoms.f90
     source/KKRhost/wunfiles.F90
     source/KKRhost/ylag.f90
+    source/KKRhost/write_gflle_npy.f90
 )
 # disable cmake auto add of 'lib' prefix to .so file
 SET_TARGET_PROPERTIES(lib_kkrhost PROPERTIES PREFIX "")
diff --git a/source/KKRhost/rhovalnew.F90 b/source/KKRhost/rhovalnew.F90
index 864896ebf18b56dc621b8890c2735ea0cafd122b..6bedba2072cbb5b26c723275c2638ea07b6df83b 100644
--- a/source/KKRhost/rhovalnew.F90
+++ b/source/KKRhost/rhovalnew.F90
@@ -41,7 +41,7 @@ contains
     use :: mod_save_wavefun, only: t_wavefunctions, read_wavefunc
     use :: mod_runoptions, only: calc_exchange_couplings, calc_gmat_lm_full, disable_tmat_sratrick, fix_nonco_angles, &
                                  use_qdos, write_complex_qdos, write_pkkr_operators, write_DOS_lm, set_cheby_nospeedup, &
-                                 decouple_spin_cheby, disable_print_serialnumber, set_gmat_to_zero, use_rllsll
+                                 decouple_spin_cheby, disable_print_serialnumber, gflle_to_npy, set_gmat_to_zero, use_rllsll
     use :: mod_version_info, only: version_print_header
     use :: global_variables, only: lmmaxd, iemxd, ncleb, lmxspd, irmd, ntotd, nrmaxd, lmpotd, nspotd, nfund, korbit, mmaxd, nspind, angles_cutoff
     use :: mod_constants, only: czero, cvlight, cone, pi, ci
@@ -62,6 +62,7 @@ contains
     use :: mod_wunfiles, only: t_params
     use :: mod_bfield, only: add_bfield
     use :: mod_torque, only: calc_torque
+    use mod_write_gflle, only: write_gflle_to_npy
 
     implicit none
 
@@ -409,7 +410,7 @@ contains
 #else
     i1_myrank = i1                 ! lmlm-dos ruess
 #endif
-    if ((calc_gmat_lm_full) .and. (i1_myrank==1)) then ! lmlm-dos ruess
+    if ((calc_gmat_lm_full) .and. (i1_myrank==1) .and..not.gflle_to_npy) then ! lmlm-dos ruess
       lrecgflle = nspin*(1+korbit)*lmmaxd*lmmaxd*ielast*nqdos ! lmlm-dos ruess
       open (91, access='direct', recl=lrecgflle, file='gflle' & ! lmlm-dos ruess
         , form='unformatted', status='replace', err=110, iostat=ierr) ! lmlm-dos ruess
@@ -1035,7 +1036,11 @@ contains
         if (t_inc%i_write>0) then  ! lmlm-dos
           write (1337, *) 'gflle:', shape(gflle), shape(gflle_part), lrecgflle ! lmlm-dos
         end if                     ! lmlm-dos
-        write (91, rec=i1) gflle   ! lmlm-dos
+        if (gflle_to_npy) then
+          call write_gflle_to_npy(lmmaxd, ielast, nqdos, i1, gflle)
+        else
+          write (91, rec=i1) gflle   ! lmlm-dos
+        end if
       end if                       ! lmlm-dos
 
       allocate (rhotemp(irmdnew,lmpotd), stat=i_stat)
diff --git a/source/KKRhost/write_gflle_npy.f90 b/source/KKRhost/write_gflle_npy.f90
new file mode 100644
index 0000000000000000000000000000000000000000..17afbec6ac7a3c9f8d47605ff3292efc3efaa5f1
--- /dev/null
+++ b/source/KKRhost/write_gflle_npy.f90
@@ -0,0 +1,31 @@
+module mod_write_gflle 
+
+  implicit none
+
+contains
+
+  !-------------------------------------------------------------------------------
+  !> Summary: Write gflle file out in npy format
+  !> Author: Philipp Rüßmann
+  !> Category: writeout
+  !> Deprecated: False 
+  !> Creates one file per atom and energy, otherwise file can be very large which
+  !> might be problematic in post-processing
+  !-------------------------------------------------------------------------------
+  subroutine write_gflle_to_npy(lmmaxd, ielast, nqdos, i1, gflle)
+
+    use mod_datatypes, only: dp
+    use m_npy, only: save_npy
+    implicit none
+    integer, intent(in) :: lmmaxd, ielast, nqdos, i1
+    complex (kind=dp) :: gflle(lmmaxd,lmmaxd,ielast,nqdos)
+    character (len=100) :: filename
+    integer :: ie
+    do ie = 1, ielast
+      write(filename, "(A,1I0.3,A,1I0.3,A)") "gllke.", I1, ".", IE, ".npy"
+      call save_npy(trim(filename), gflle(:,:,ie, :))
+    end do
+
+  end subroutine write_gflle_to_npy
+
+end module mod_write_gflle 
diff --git a/source/common/runoptions.F90 b/source/common/runoptions.F90
index 4d3e4ae03d0eecc5d8cb22766110b8a796ded479..ae760610d160b297fd1a3a6e3cdf7d603d23dcf6 100644
--- a/source/common/runoptions.F90
+++ b/source/common/runoptions.F90
@@ -34,6 +34,7 @@ module mod_runoptions
   logical :: fix_nonco_angles = .false.                !! fix direction of non-collinear magnetic moments (Chebychev solver) (former: 'FIXMOM')
   logical :: formatted_file = .false.                  !! write files ascii-format. only effective with some other write-options (former: 'fileverb')
   logical :: write_npy = .false.                       !! write files in npy format (python readable) (former: 'writenpy')
+  logical :: gflle_to_npy = .false.                    !! write gflle file to npy instead of to gflle file (contains G(k,e)LL'
   logical :: impurity_operator_only = .false.          !! only for `write_pkkr_operators`: disable costly recalculation of host operators (former: 'IMP_ONLY')
   logical :: modify_soc_Dirac = .false.                !! modify SOC for Dirac solver (former: 'SOC')
   logical :: no_madelung = .false.                     !! do not add some energy terms (coulomb, XC, eff. pot.) to total energy (former: 'NoMadel')
@@ -167,6 +168,7 @@ module mod_runoptions
     call set_runoption(write_lloyd_tralpha_file      , '<write_lloyd_tralpha_file>'      , '<wrttral>' )
     call set_runoption(write_lloyd_cdos_file         , '<write_lloyd_cdos_file>'         , '<wrtcdos>' )
     call set_runoption(calc_gmat_lm_full             , '<calc_gmat_lm_full>'             , '<lmlm-dos>')
+    call set_runoption(gflle_to_npy                  , '<gflle_to_npy>')
     call set_runoption(simulate_asa                  , '<simulate_asa>'                  , '<simulasa>')
     call set_runoption(use_readcpa                   , '<use_readcpa>'                   , '<readcpa>' )
     call set_runoption(print_kpoints                 , '<print_kpoints>'                 , '<BZKP>'    )
@@ -716,6 +718,7 @@ module mod_runoptions
     call mpi_bcast(write_lloyd_tralpha_file      , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(write_lloyd_cdos_file         , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(calc_gmat_lm_full             , 1, mpi_logical, master, mpi_comm_world, ierr)
+    call mpi_bcast(gflle_to_npy                  , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(simulate_asa                  , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(use_readcpa                   , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(print_kpoints                 , 1, mpi_logical, master, mpi_comm_world, ierr)
@@ -828,6 +831,7 @@ module mod_runoptions
     write(iounit, '(A35,1x,1L,3x,A)') '<calc_exchange_couplings>=', calc_exchange_couplings, "calculate magnetic exchange coupling parameters (former: 'XCPL')"
     write(iounit, '(A35,1x,1L,3x,A)') '<calc_exchange_couplings_energy>=', calc_exchange_couplings_energy, "write energy-resolved Jij-files also if npol/=0 (former: 'Jijenerg')"
     write(iounit, '(A35,1x,1L,3x,A)') '<calc_gmat_lm_full>=', calc_gmat_lm_full, "calculate all lm-lm components of systems greens function and store to file `gflle` (former: 'lmlm-dos')"
+    write(iounit, '(A35,1x,1L,3x,A)') '<gflle_to_npy>=', gflle_to_npy, "Write G_LL'(k,E) to npy files, one file per atom and energy"
     write(iounit, '(A35,1x,1L,3x,A)') '<dirac_scale_SpeefOfLight>=', dirac_scale_SpeefOfLight, "scale the speed of light for Dirac solver (former: 'CSCALE')"
     write(iounit, '(A35,1x,1L,3x,A)') '<disable_charge_neutrality>=', disable_charge_neutrality, "no charge neutrailty required: leaving Fermi level unaffected (former: 'no-neutr')"
     write(iounit, '(A35,1x,1L,3x,A)') '<disable_print_serialnumber>=', disable_print_serialnumber, "deactivate writing of serial number and version information to files (for backwards compatibility) (former: 'noserial')"
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