diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5d26a228..564b3b0d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,7 @@ include(FetchContent)
 include(cmake/options.cmake)
 include(cmake/compilers.cmake)
 
+include(cmake/functional.cmake)
 include(cmake/h5fortran.cmake)
 include(cmake/json.cmake)
 
@@ -62,7 +63,13 @@ add_library(neural
   src/nf/io/nf_io_hdf5.f90
   src/nf/io/nf_io_hdf5_submodule.f90
 )
-target_link_libraries(neural PRIVATE h5fortran::h5fortran HDF5::HDF5 jsonfortran::jsonfortran)
+
+target_link_libraries(neural PRIVATE
+  functional::functional
+  h5fortran::h5fortran
+  HDF5::HDF5
+  jsonfortran::jsonfortran
+)
 
 install(TARGETS neural)
 
diff --git a/README.md b/README.md
index 085ddd96..b66561de 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 [![GitHub issues](https://img.shields.io/github/issues/modern-fortran/neural-fortran.svg)](https://github.com/modern-fortran/neural-fortran/issues)
 
-A parallel neural net microframework. 
+A parallel framework for deep learning.
 Read the paper [here](https://arxiv.org/abs/1902.06714).
 
 * [Features](https://github.com/modern-fortran/neural-fortran#features)
@@ -18,9 +18,11 @@ Read the paper [here](https://arxiv.org/abs/1902.06714).
 
 * Dense, fully connected neural layers
 * Convolutional and max-pooling layers (experimental, forward propagation only)
+* Flatten layers (forward and backward pass)
+* Loading dense and convolutional models from Keras h5 files
 * Stochastic and mini-batch gradient descent for back-propagation
 * Data-based parallelism
-* Several activation functions
+* Several activation functions and their derivatives
 
 ### Available layer types
 
@@ -48,16 +50,18 @@ Required dependencies are:
 * A Fortran compiler
 * [HDF5](https://www.hdfgroup.org/downloads/hdf5/)
   (must be provided by the OS package manager or your own build from source)
-* [h5fortran](https://github.com/geospace-code/h5fortran),
+* [functional-fortran](https://github.com/wavebitscientific/functional-fortran),
+  [h5fortran](https://github.com/geospace-code/h5fortran),
   [json-fortran](https://github.com/jacobwilliams/json-fortran)
-  (both handled by neural-fortran's build systems, no need for a manual install)
+  (all handled by neural-fortran's build systems, no need for a manual install)
 * [fpm](https://github.com/fortran-lang/fpm) or
   [CMake](https://cmake.org) for building the code
 
 Optional dependencies are:
 
 * OpenCoarrays (for parallel execution with GFortran)
-* BLAS, MKL (optional)
+* BLAS, MKL, or similar (for offloading `matmul` and `dot_product` calls)
+* curl (for downloading testing and example datasets)
 
 Compilers tested include:
 
@@ -200,13 +204,15 @@ examples, in increasing level of complexity:
   dataset
 4. [cnn](example/cnn.f90): Creating and running forward a simple CNN using
   `input`, `conv2d`, `maxpool2d`, `flatten`, and `dense` layers.
-5. [mnist_from_keras](example/mnist_from_keras.f90): Creating a pre-trained
-  model from a Keras HDF5 file.
+5. [dense_from_keras](example/dense_from_keras.f90): Creating a pre-trained
+  dense model from a Keras HDF5 file and running the inference.
+6. [cnn_from_keras](example/cnn_from_keras.f90): Creating a pre-trained
+  convolutional model from a Keras HDF5 file and running the inference.
 
 The examples also show you the extent of the public API that's meant to be
 used in applications, i.e. anything from the `nf` module.
 
-The MNIST example uses [curl](https://curl.se/) to download the dataset,
+Examples 3-6 rely on [curl](https://curl.se/) to download the needed datasets,
 so make sure you have it installed on your system.
 Most Linux OSs have it out of the box.
 The dataset will be downloaded only the first time you run the example in any
diff --git a/cmake/functional.cmake b/cmake/functional.cmake
new file mode 100644
index 00000000..d8f21200
--- /dev/null
+++ b/cmake/functional.cmake
@@ -0,0 +1,18 @@
+FetchContent_Declare(functional
+  GIT_REPOSITORY https://github.com/wavebitscientific/functional-fortran
+  GIT_TAG 0.6.1
+  GIT_SHALLOW true
+)
+
+FetchContent_Populate(functional)
+
+add_library(functional ${functional_SOURCE_DIR}/src/functional.f90)
+target_include_directories(functional PUBLIC
+$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/include>
+$<INSTALL_INTERFACE:include>
+)
+
+add_library(functional::functional INTERFACE IMPORTED GLOBAL)
+target_link_libraries(functional::functional INTERFACE functional)
+
+install(TARGETS functional)
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index b8442417..7cb24732 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -1,4 +1,16 @@
-foreach(execid cnn mnist mnist_from_keras simple sine)
+foreach(execid
+  cnn
+  cnn_from_keras
+  dense_from_keras
+  mnist
+  simple
+  sine
+)
   add_executable(${execid} ${execid}.f90)
-  target_link_libraries(${execid} PRIVATE neural h5fortran::h5fortran jsonfortran::jsonfortran ${LIBS})
+  target_link_libraries(${execid} PRIVATE
+    neural
+    h5fortran::h5fortran
+    jsonfortran::jsonfortran
+    ${LIBS}
+  )
 endforeach()
diff --git a/example/cnn_from_keras.f90 b/example/cnn_from_keras.f90
new file mode 100644
index 00000000..5be47854
--- /dev/null
+++ b/example/cnn_from_keras.f90
@@ -0,0 +1,58 @@
+program cnn_from_keras
+
+  ! This example demonstrates loading a convolutional model
+  ! pre-trained on the MNIST dataset from a Keras HDF5
+  ! file and running an inferrence on the testing dataset.
+
+  use nf, only: network, label_digits, load_mnist
+  use nf_datasets, only: download_and_unpack, keras_cnn_mnist_url
+
+  implicit none
+
+  type(network) :: net
+  real, allocatable :: training_images(:,:), training_labels(:)
+  real, allocatable :: validation_images(:,:), validation_labels(:)
+  real, allocatable :: testing_images(:,:), testing_labels(:)
+  character(*), parameter :: keras_cnn_path = 'keras_cnn_mnist.h5'
+  logical :: file_exists
+  real :: acc
+
+  inquire(file=keras_cnn_path, exist=file_exists)
+  if (.not. file_exists) call download_and_unpack(keras_cnn_mnist_url)
+
+  call load_mnist(training_images, training_labels, &
+                  validation_images, validation_labels, &
+                  testing_images, testing_labels)
+
+  print '("Loading a pre-trained CNN model from Keras")'
+  print '(60("="))'
+
+  net = network(keras_cnn_path)
+
+  call net % print_info()
+
+  if (this_image() == 1) then
+    acc = accuracy( &
+      net, &
+      reshape(testing_images(:,:), shape=[1,28,28,size(testing_images,2)]), &
+      label_digits(testing_labels) &
+    )
+    print '(a,f5.2,a)', 'Accuracy: ', acc * 100, ' %'
+  end if
+
+contains
+
+  real function accuracy(net, x, y)
+    type(network), intent(in out) :: net
+    real, intent(in) :: x(:,:,:,:), y(:,:)
+    integer :: i, good
+    good = 0
+    do i = 1, size(x, dim=4)
+      if (all(maxloc(net % output(x(:,:,:,i))) == maxloc(y(:,i)))) then
+        good = good + 1
+      end if
+    end do
+    accuracy = real(good) / size(x, dim=4)
+  end function accuracy
+
+end program cnn_from_keras
diff --git a/example/mnist_from_keras.f90 b/example/dense_from_keras.f90
similarity index 64%
rename from example/mnist_from_keras.f90
rename to example/dense_from_keras.f90
index b5a67e41..cf473d6d 100644
--- a/example/mnist_from_keras.f90
+++ b/example/dense_from_keras.f90
@@ -1,10 +1,11 @@
-program mnist_from_keras
+program dense_from_keras
 
-  ! This example demonstrates loading a pre-trained MNIST model from Keras
-  ! from an HDF5 file and running an inferrence on the testing dataset.
+  ! This example demonstrates loading a dense model
+  ! pre-trained on the MNIST dataset from a Keras HDF5
+  ! file and running an inferrence on the testing dataset.
 
   use nf, only: network, label_digits, load_mnist
-  use nf_datasets, only: download_and_unpack, keras_model_dense_mnist_url
+  use nf_datasets, only: download_and_unpack, keras_dense_mnist_url
 
   implicit none
 
@@ -12,20 +13,20 @@ program mnist_from_keras
   real, allocatable :: training_images(:,:), training_labels(:)
   real, allocatable :: validation_images(:,:), validation_labels(:)
   real, allocatable :: testing_images(:,:), testing_labels(:)
-  character(*), parameter :: test_data_path = 'keras_dense_mnist.h5'
+  character(*), parameter :: keras_dense_path = 'keras_dense_mnist.h5'
   logical :: file_exists
 
-  inquire(file=test_data_path, exist=file_exists)
-  if (.not. file_exists) call download_and_unpack(keras_model_dense_mnist_url)
+  inquire(file=keras_dense_path, exist=file_exists)
+  if (.not. file_exists) call download_and_unpack(keras_dense_mnist_url)
 
   call load_mnist(training_images, training_labels, &
                   validation_images, validation_labels, &
                   testing_images, testing_labels)
 
-  print '("Loading a pre-trained MNIST model from Keras")'
+  print '("Loading a pre-trained dense model from Keras")'
   print '(60("="))'
 
-  net = network(test_data_path)
+  net = network(keras_dense_path)
 
   call net % print_info()
 
@@ -48,4 +49,4 @@ real function accuracy(net, x, y)
     accuracy = real(good) / size(x, dim=2)
   end function accuracy
 
-end program mnist_from_keras
+end program dense_from_keras
diff --git a/fpm.toml b/fpm.toml
index 90ee5d9d..0cabfcac 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -1,5 +1,5 @@
 name = "neural-fortran"
-version = "0.5.0"
+version = "0.6.0"
 license = "MIT"
 author = "Milan Curcic"
 maintainer = "milancurcic@hey.com"
@@ -10,5 +10,6 @@ external-modules = "hdf5"
 link = ["hdf5", "hdf5_fortran"]
 
 [dependencies]
+functional = { git = "https://github.com/wavebitscientific/functional-fortran" }
 h5fortran = { git = "https://github.com/geospace-code/h5fortran" }
 json-fortran = { git = "https://github.com/jacobwilliams/json-fortran" }
diff --git a/src/nf/io/nf_io_hdf5.f90 b/src/nf/io/nf_io_hdf5.f90
index aa6e20f1..ac74524e 100644
--- a/src/nf/io/nf_io_hdf5.f90
+++ b/src/nf/io/nf_io_hdf5.f90
@@ -31,7 +31,7 @@ module subroutine get_hdf5_dataset_real32_1d(filename, object_name, values)
         !! HDF5 file name
       character(*), intent(in) :: object_name
         !! Object (dataset) name
-      real(real32), allocatable, intent(in out) :: values(:)
+      real(real32), allocatable, intent(out) :: values(:)
         !! Array to store the dataset values into
     end subroutine get_hdf5_dataset_real32_1d
 
@@ -41,10 +41,20 @@ module subroutine get_hdf5_dataset_real32_2d(filename, object_name, values)
         !! HDF5 file name
       character(*), intent(in) :: object_name
         !! Object (dataset) name
-      real(real32), allocatable, intent(in out) :: values(:,:)
+      real(real32), allocatable, intent(out) :: values(:,:)
         !! Array to store the dataset values into
     end subroutine get_hdf5_dataset_real32_2d
   
+    module subroutine get_hdf5_dataset_real32_4d(filename, object_name, values)
+      !! Read a 4-d real32 array from an HDF5 dataset.
+      character(*), intent(in) :: filename
+        !! HDF5 file name
+      character(*), intent(in) :: object_name
+        !! Object (dataset) name
+      real(real32), allocatable, intent(out) :: values(:,:,:,:)
+        !! Array to store the dataset values into
+    end subroutine get_hdf5_dataset_real32_4d
+
   end interface get_hdf5_dataset
 
 end module nf_io_hdf5
diff --git a/src/nf/io/nf_io_hdf5_submodule.f90 b/src/nf/io/nf_io_hdf5_submodule.f90
index 263f8384..8cb6f760 100644
--- a/src/nf/io/nf_io_hdf5_submodule.f90
+++ b/src/nf/io/nf_io_hdf5_submodule.f90
@@ -50,7 +50,7 @@ module subroutine get_hdf5_dataset_real32_1d(filename, object_name, values)
 
     character(*), intent(in) :: filename
     character(*), intent(in) :: object_name
-    real(real32), allocatable, intent(in out) :: values(:)
+    real(real32), allocatable, intent(out) :: values(:)
 
     type(hdf5_file) :: f
     integer(int64), allocatable :: dims(:)
@@ -58,15 +58,7 @@ module subroutine get_hdf5_dataset_real32_1d(filename, object_name, values)
     call f % open(filename, 'r')
     call f % shape(object_name, dims)
 
-    ! If values is already allocated, re-allocate only if incorrect shape
-    if (allocated(values)) then
-      if (.not. all(shape(values) == dims)) then
-        deallocate(values)
-        allocate(values(dims(1)))
-      end if
-    else
-      allocate(values(dims(1)))
-    end if
+    allocate(values(dims(1)))
 
     call f % read(object_name, values)
     call f % close()
@@ -78,7 +70,7 @@ module subroutine get_hdf5_dataset_real32_2d(filename, object_name, values)
 
     character(*), intent(in) :: filename
     character(*), intent(in) :: object_name
-    real(real32), allocatable, intent(in out) :: values(:,:)
+    real(real32), allocatable, intent(out) :: values(:,:)
 
     type(hdf5_file) :: f
     integer(int64), allocatable :: dims(:)
@@ -86,19 +78,34 @@ module subroutine get_hdf5_dataset_real32_2d(filename, object_name, values)
     call f % open(filename, 'r')
     call f % shape(object_name, dims)
 
-    ! If values is already allocated, re-allocate only if incorrect shape
-    if (allocated(values)) then
-      if (.not. all(shape(values) == dims)) then
-        deallocate(values)
-        allocate(values(dims(1), dims(2)))
-      end if
-    else
-      allocate(values(dims(1), dims(2)))
-    end if
+    allocate(values(dims(1), dims(2)))
 
     call f % read(object_name, values)
     call f % close()
 
+    ! Transpose the array to respect Keras's storage order
+    values = transpose(values)
+
   end subroutine get_hdf5_dataset_real32_2d
 
+
+  module subroutine get_hdf5_dataset_real32_4d(filename, object_name, values)
+
+    character(*), intent(in) :: filename
+    character(*), intent(in) :: object_name
+    real(real32), allocatable, intent(out) :: values(:,:,:,:)
+
+    type(hdf5_file) :: f
+    integer(int64), allocatable :: dims(:)
+
+    call f % open(filename, 'r')
+    call f % shape(object_name, dims)
+
+    allocate(values(dims(1), dims(2), dims(3), dims(4)))
+
+    call f % read(object_name, values)
+    call f % close()
+
+  end subroutine get_hdf5_dataset_real32_4d
+
 end submodule nf_io_hdf5_submodule
diff --git a/src/nf/nf_datasets.f90 b/src/nf/nf_datasets.f90
index 3a19b1cf..2df35cb9 100644
--- a/src/nf/nf_datasets.f90
+++ b/src/nf/nf_datasets.f90
@@ -8,13 +8,19 @@ module nf_datasets
 
   private
 
-  public :: download_and_unpack, keras_model_dense_mnist_url, mnist_url
+  public :: &
+    download_and_unpack, &
+    keras_cnn_mnist_url, &
+    keras_dense_mnist_url, &
+    mnist_url
 
   character(*), parameter :: keras_snippets_baseurl = &
     'https://github.com/neural-fortran/keras-snippets/files'
   character(*), parameter :: neural_fortran_baseurl = &
     'https://github.com/modern-fortran/neural-fortran/files'
-  character(*), parameter :: keras_model_dense_mnist_url = &
+  character(*), parameter :: keras_cnn_mnist_url = &
+    keras_snippets_baseurl // '/8892585/keras_cnn_mnist.tar.gz'
+  character(*), parameter :: keras_dense_mnist_url = &
     keras_snippets_baseurl // '/8788739/keras_dense_mnist.tar.gz'
   character(*), parameter :: mnist_url = &
     neural_fortran_baseurl // '/8498876/mnist.tar.gz'
diff --git a/src/nf/nf_keras.f90 b/src/nf/nf_keras.f90
index 5d5ec2d9..04c9385b 100644
--- a/src/nf/nf_keras.f90
+++ b/src/nf/nf_keras.f90
@@ -9,12 +9,26 @@ module nf_keras
   public :: get_keras_h5_layers, keras_layer
 
   type :: keras_layer
-    !! Intermediate container to convey the Keras layer information
-    !! to neural-fortran layer constructors.
+
+    !! Intermediate container to convey the Keras layer
+    !! information to neural-fortran layer constructors.
+
+    ! General metadata that applies to any (or most) layers
     character(:), allocatable :: class
     character(:), allocatable :: name
     character(:), allocatable :: activation
-    integer, allocatable :: num_elements(:)
+
+    ! Dense
+    integer, allocatable :: units(:)
+
+    ! Conv2D
+    integer :: filters
+    integer, allocatable :: kernel_size(:)
+
+    ! MaxPooling2D
+    integer, allocatable :: pool_size(:)
+    integer, allocatable :: strides(:)
+
   end type keras_layer
 
   interface
diff --git a/src/nf/nf_keras_submodule.f90 b/src/nf/nf_keras_submodule.f90
index b97c12da..161b38be 100644
--- a/src/nf/nf_keras_submodule.f90
+++ b/src/nf/nf_keras_submodule.f90
@@ -1,5 +1,6 @@
 submodule(nf_keras) nf_keras_submodule
 
+  use functional, only: reverse
   use json_module, only: json_core, json_value
   use nf_io_hdf5, only: hdf5_attribute_string
 
@@ -18,7 +19,7 @@ module function get_keras_h5_layers(filename) result(res)
       model_config_json, layers_json, layer_json, layer_config_json
 
     real, allocatable :: tmp_array(:)
-    integer :: n, num_layers, num_elements
+    integer :: n, num_layers, units
     logical :: found
 
     model_config_string = hdf5_attribute_string(filename, '.', 'model_config')
@@ -51,13 +52,36 @@ module function get_keras_h5_layers(filename) result(res)
 
         case('InputLayer')
           call json % get(layer_config_json, 'batch_input_shape', tmp_array)
-          res(n) % num_elements = [tmp_array(2)]
+          res(n) % units = reverse(tmp_array(2:)) ! skip the 1st (batch) dim
       
         case('Dense')
-          call json % get(layer_config_json, 'units',  num_elements, found)
-          res(n) % num_elements = [num_elements]
+          call json % get(layer_config_json, 'units', units, found)
+          res(n) % units = [units]
           call json % get(layer_config_json, 'activation', res(n) % activation)
-      
+
+        case('Flatten')
+          ! Nothing to read here; merely a placeholder.
+          continue
+
+        case('Conv2D')
+          call json % get(layer_config_json, &
+            'filters', res(n) % filters, found)
+          call json % get(layer_config_json, &
+            'kernel_size', res(n) % kernel_size, found)
+          call json % get(layer_config_json, &
+            'activation', res(n) % activation)
+          ! Reverse to account for C -> Fortran order
+          res(n) % kernel_size = reverse(res(n) % kernel_size)
+
+        case('MaxPooling2D')
+          call json % get(layer_config_json, &
+            'pool_size',  res(n) % pool_size, found)
+          call json % get(layer_config_json, &
+            'strides',  res(n) % strides, found)
+          ! Reverse to account for C -> Fortran order
+          res(n) % pool_size = reverse(res(n) % pool_size)
+          res(n) % strides = reverse(res(n) % strides)
+
         case default
           error stop 'This Keras layer is not supported'
 
diff --git a/src/nf/nf_network_submodule.f90 b/src/nf/nf_network_submodule.f90
index 6333e5eb..c9f865ff 100644
--- a/src/nf/nf_network_submodule.f90
+++ b/src/nf/nf_network_submodule.f90
@@ -1,13 +1,15 @@
 submodule(nf_network) nf_network_submodule
 
+  use nf_conv2d_layer, only: conv2d_layer
   use nf_dense_layer, only: dense_layer
   use nf_flatten_layer, only: flatten_layer
   use nf_input1d_layer, only: input1d_layer
   use nf_input3d_layer, only: input3d_layer
+  use nf_maxpool2d_layer, only: maxpool2d_layer
   use nf_io_hdf5, only: get_hdf5_dataset
   use nf_keras, only: get_keras_h5_layers, keras_layer
   use nf_layer, only: layer
-  use nf_layer_constructors, only: dense, input
+  use nf_layer_constructors, only: conv2d, dense, flatten, input, maxpool2d
   use nf_loss, only: quadratic_derivative
   use nf_optimizers, only: sgd
   use nf_parallel, only: tile_indices
@@ -68,19 +70,51 @@ module function network_from_keras(filename) result(res)
 
       select case(keras_layers(n) % class)
 
+        case('Conv2D')
+
+          if (keras_layers(n) % kernel_size(1) &
+            /= keras_layers(n) % kernel_size(2)) &
+            error stop 'Non-square kernel in conv2d layer not supported.'
+
+          layers(n) = conv2d( &
+            keras_layers(n) % filters, &
+            !FIXME add support for non-square kernel
+            keras_layers(n) % kernel_size(1), &
+            keras_layers(n) % activation &
+          )
+
+        case('Dense')
+          layers(n) = dense( &
+            keras_layers(n) % units(1), &
+            keras_layers(n) % activation &
+          )
+
+        case('Flatten')
+          layers(n) = flatten()
+
         case('InputLayer')
-          if (size(keras_layers(n) % num_elements) == 1) then
+          if (size(keras_layers(n) % units) == 1) then
             ! input1d
-            layers(n) = input(keras_layers(n) % num_elements(1))
+            layers(n) = input(keras_layers(n) % units(1))
           else
             ! input3d
-            layers(n) = input(keras_layers(n) % num_elements)
+            layers(n) = input(keras_layers(n) % units)
           end if
 
-        case('Dense')
-          layers(n) = dense( &
-            keras_layers(n) % num_elements(1), &
-            keras_layers(n) % activation &
+        case('MaxPooling2D')
+
+          if (keras_layers(n) % pool_size(1) &
+            /= keras_layers(n) % pool_size(2)) &
+            error stop 'Non-square pool in maxpool2d layer not supported.'
+
+          if (keras_layers(n) % strides(1) &
+            /= keras_layers(n) % strides(2)) &
+            error stop 'Unequal strides in maxpool2d layer are not supported.'
+
+          layers(n) = maxpool2d( &
+            !FIXME add support for non-square pool and stride
+            keras_layers(n) % pool_size(1), &
+            keras_layers(n) % strides(1) &
           )
 
         case default
@@ -98,36 +132,45 @@ module function network_from_keras(filename) result(res)
 
       layer_name = keras_layers(n) % name
 
-      if (keras_layers(n) % class == 'Dense') then
-        select type(this_layer => res % layers(n) % p)
+      select type(this_layer => res % layers(n) % p)
 
-          type is(dense_layer)
+        type is(conv2d_layer)
+          ! Read biases from file
+          object_name = '/model_weights/' // layer_name // '/' &
+            // layer_name // '/bias:0'
+          call get_hdf5_dataset(filename, object_name, this_layer % biases)
 
-            ! Read biases from file
-            object_name = '/model_weights/' // layer_name // '/' &
-              // layer_name // '/bias:0'
-            call get_hdf5_dataset(filename, object_name, this_layer % biases)
+          ! Read weights from file
+          object_name = '/model_weights/' // layer_name // '/' &
+            // layer_name // '/kernel:0'
+          call get_hdf5_dataset(filename, object_name, this_layer % kernel)
 
-            ! Read weights from file
-            object_name = '/model_weights/' // layer_name // '/' &
-              // layer_name // '/kernel:0'
-            call get_hdf5_dataset(filename, object_name, this_layer % weights)
+        type is(dense_layer)
 
-            ! TODO Multidimensional arrays are stored in HDF5 in C-order.
-            ! TODO Here we transpose the array to get to the Fortran order.
-            ! TODO There may be a way to do this without re-allocating.
-            ! TODO It probably doesn't matter much since we do this once.
-            ! TODO Figure it out later.
-            this_layer % weights = transpose(this_layer % weights)
+          ! Read biases from file
+          object_name = '/model_weights/' // layer_name // '/' &
+            // layer_name // '/bias:0'
+          call get_hdf5_dataset(filename, object_name, this_layer % biases)
 
-          class default
-            error stop 'Internal error in network_from_keras(); ' &
-              // 'mismatch in layer types between the Keras and ' &
-              // 'neural-fortran model layers.'
+          ! Read weights from file
+          object_name = '/model_weights/' // layer_name // '/' &
+            // layer_name // '/kernel:0'
+          call get_hdf5_dataset(filename, object_name, this_layer % weights)
 
-        end select
+        type is(flatten_layer)
+          ! Nothing to do
+          continue
+
+        type is(maxpool2d_layer)
+          ! Nothing to do
+          continue
 
-    end if
+        class default
+          error stop 'Internal error in network_from_keras(); ' &
+            // 'mismatch in layer types between the Keras and ' &
+            // 'neural-fortran model layers.'
+
+      end select
 
   end do
 
@@ -216,6 +259,8 @@ module function output_1d(self, input) result(res)
         res = output_layer % output
       type is(flatten_layer)
         res = output_layer % output
+      class default
+        error stop 'network % output not implemented for this output layer'
     end select
 
   end function output_1d
@@ -232,10 +277,15 @@ module function output_3d(self, input) result(res)
     call self % forward(input)
 
     select type(output_layer => self % layers(num_layers) % p)
+      type is(conv2d_layer)
+        !FIXME flatten the result for now; find a better solution
+        res = pack(output_layer % output, .true.)
       type is(dense_layer)
         res = output_layer % output
       type is(flatten_layer)
         res = output_layer % output
+      class default
+        error stop 'network % output not implemented for this output layer'
     end select
 
   end function output_3d
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f30c6d79..5538d198 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,4 +1,17 @@
-foreach(execid input1d_layer input3d_layer dense_layer conv2d_layer maxpool2d_layer flatten_layer dense_network dense_network_from_keras conv2d_network io_hdf5 keras_read_model)
+foreach(execid
+  input1d_layer
+  input3d_layer
+  dense_layer
+  conv2d_layer
+  maxpool2d_layer
+  flatten_layer
+  dense_network
+  io_hdf5
+  keras_read_model
+  dense_network_from_keras
+  cnn_from_keras
+  conv2d_network
+  )
   add_executable(test_${execid} test_${execid}.f90)
   target_link_libraries(test_${execid} PRIVATE neural h5fortran::h5fortran jsonfortran::jsonfortran ${LIBS})
 
diff --git a/test/test_cnn_from_keras.f90 b/test/test_cnn_from_keras.f90
new file mode 100644
index 00000000..1c47325e
--- /dev/null
+++ b/test/test_cnn_from_keras.f90
@@ -0,0 +1,69 @@
+program test_cnn_from_keras
+
+  use iso_fortran_env, only: stderr => error_unit
+  use nf, only: network
+  use nf_datasets, only: download_and_unpack, keras_cnn_mnist_url
+
+  implicit none
+
+  type(network) :: net
+  character(*), parameter :: test_data_path = 'keras_cnn_mnist.h5'
+  logical :: file_exists
+  logical :: ok = .true.
+
+  inquire(file=test_data_path, exist=file_exists)
+  if (.not. file_exists) call download_and_unpack(keras_cnn_mnist_url)
+
+  net = network(test_data_path)
+
+  block
+
+    use nf, only: load_mnist, label_digits
+
+    real, allocatable :: training_images(:,:), training_labels(:)
+    real, allocatable :: validation_images(:,:), validation_labels(:)
+    real, allocatable :: testing_images(:,:), testing_labels(:)
+    real, allocatable :: input_reshaped(:,:,:,:)
+    real :: acc
+
+    call load_mnist(training_images, training_labels, &
+                    validation_images, validation_labels, &
+                    testing_images, testing_labels)
+
+    ! Use only the first 1000 images to make the test short
+    input_reshaped = reshape(testing_images(:,:1000), shape=[1,28,28,1000])
+
+    acc = accuracy(net, input_reshaped, label_digits(testing_labels(:1000)))
+
+    if (acc < 0.97) then
+      write(stderr, '(a)') &
+        'Pre-trained network accuracy should be > 0.97.. failed'
+      ok = .false.
+    end if
+
+  end block
+
+  if (ok) then
+    print '(a)', 'test_cnn_from_keras: All tests passed.'
+  else
+    write(stderr, '(a)') &
+      'test_cnn_from_keras: One or more tests failed.'
+    stop 1
+  end if
+
+contains
+
+  real function accuracy(net, x, y)
+    type(network), intent(in out) :: net
+    real, intent(in) :: x(:,:,:,:), y(:,:)
+    integer :: i, good
+    good = 0
+    do i = 1, size(x, dim=4)
+      if (all(maxloc(net % output(x(:,:,:,i))) == maxloc(y(:,i)))) then
+        good = good + 1
+      end if
+    end do
+    accuracy = real(good) / size(x, dim=4)
+  end function accuracy
+
+end program test_cnn_from_keras
diff --git a/test/test_dense_network_from_keras.f90 b/test/test_dense_network_from_keras.f90
index 435efcec..ebd0b501 100644
--- a/test/test_dense_network_from_keras.f90
+++ b/test/test_dense_network_from_keras.f90
@@ -2,7 +2,7 @@ program test_dense_network_from_keras
 
   use iso_fortran_env, only: stderr => error_unit
   use nf, only: network
-  use nf_datasets, only: download_and_unpack, keras_model_dense_mnist_url
+  use nf_datasets, only: download_and_unpack, keras_dense_mnist_url
 
   implicit none
 
@@ -12,7 +12,7 @@ program test_dense_network_from_keras
   logical :: ok = .true.
 
   inquire(file=test_data_path, exist=file_exists)
-  if (.not. file_exists) call download_and_unpack(keras_model_dense_mnist_url)
+  if (.not. file_exists) call download_and_unpack(keras_dense_mnist_url)
 
   net = network(test_data_path)
 
diff --git a/test/test_io_hdf5.f90 b/test/test_io_hdf5.f90
index dc8940af..e7e0fc88 100644
--- a/test/test_io_hdf5.f90
+++ b/test/test_io_hdf5.f90
@@ -1,7 +1,7 @@
 program test_io_hdf5
 
   use iso_fortran_env, only: stderr => error_unit
-  use nf_datasets, only: download_and_unpack, keras_model_dense_mnist_url
+  use nf_datasets, only: download_and_unpack, keras_dense_mnist_url
   use nf_io_hdf5, only: hdf5_attribute_string, get_hdf5_dataset
 
   implicit none
@@ -14,7 +14,7 @@ program test_io_hdf5
   logical :: ok = .true.
 
   inquire(file=test_data_path, exist=file_exists)
-  if (.not. file_exists) call download_and_unpack(keras_model_dense_mnist_url)
+  if (.not. file_exists) call download_and_unpack(keras_dense_mnist_url)
 
   attr = hdf5_attribute_string(test_data_path, '.', 'backend')
 
@@ -42,7 +42,7 @@ program test_io_hdf5
     write(stderr, '(a)') 'HDF5 1-d dataset dims inquiry is correct.. failed'
   end if
 
-  if (.not. all(shape(weights) == [30, 784])) then
+  if (.not. all(shape(weights) == [784, 30])) then
     ok = .false.
     write(stderr, '(a)') 'HDF5 2-d dataset dims inquiry is correct.. failed'
   end if
diff --git a/test/test_keras_read_model.f90 b/test/test_keras_read_model.f90
index 87d81ac4..ba021a5a 100644
--- a/test/test_keras_read_model.f90
+++ b/test/test_keras_read_model.f90
@@ -1,28 +1,26 @@
 program test_keras_read_model
 
   use iso_fortran_env, only: stderr => error_unit
-  use nf_datasets, only: download_and_unpack, keras_model_dense_mnist_url
+  use nf_datasets, only: download_and_unpack, keras_dense_mnist_url, &
+    keras_cnn_mnist_url
   use nf_keras, only: get_keras_h5_layers, keras_layer
-  use nf, only: layer, network, dense, input
 
   implicit none
 
-  character(:), allocatable :: model_config_string
-  character(*), parameter :: test_data_path = 'keras_dense_mnist.h5'
+  character(*), parameter :: keras_dense_path = 'keras_dense_mnist.h5'
+  character(*), parameter :: keras_cnn_path = 'keras_cnn_mnist.h5'
 
   type(keras_layer), allocatable :: keras_layers(:)
 
-  type(layer), allocatable :: layers(:)
-  type(network) :: net
-
-  integer :: n
   logical :: file_exists
   logical :: ok = .true.
 
-  inquire(file=test_data_path, exist=file_exists)
-  if (.not. file_exists) call download_and_unpack(keras_model_dense_mnist_url)
+  ! First test the dense model
+
+  inquire(file=keras_dense_path, exist=file_exists)
+  if (.not. file_exists) call download_and_unpack(keras_dense_mnist_url)
 
-  keras_layers = get_keras_h5_layers(test_data_path)
+  keras_layers = get_keras_h5_layers(keras_dense_path)
 
   if (size(keras_layers) /= 3) then
     ok = .false.
@@ -34,7 +32,7 @@ program test_keras_read_model
     write(stderr, '(a)') 'Keras first layer should be InputLayer.. failed'
   end if
 
-  if (.not. all(keras_layers(1) % num_elements == [784])) then
+  if (.not. all(keras_layers(1) % units == [784])) then
     ok = .false.
     write(stderr, '(a)') 'Keras first layer should have 784 elements.. failed'
   end if
@@ -51,6 +49,55 @@ program test_keras_read_model
       'Keras second and third layers should be dense.. failed'
   end if
 
+  ! Now testing for the CNN model
+
+  inquire(file=keras_cnn_path, exist=file_exists)
+  if (.not. file_exists) call download_and_unpack(keras_cnn_mnist_url)
+
+  keras_layers = get_keras_h5_layers(keras_cnn_path)
+
+  if (.not. all(keras_layers(1) % units == [1, 28, 28])) then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN input layer shape is expected.. failed'
+  end if
+
+  if (.not. keras_layers(2) % class == 'Conv2D') then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN second layer is Conv2D.. failed'
+  end if
+
+  if (.not. keras_layers(2) % filters == 8) then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN second layer number of filters is expected.. failed'
+  end if
+
+  if (.not. all(keras_layers(2) % kernel_size == [3, 3])) then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN second layer kernel_size is expected.. failed'
+  end if
+
+  if (.not. keras_layers(3) % class == 'MaxPooling2D') then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN third layer is MaxPooling2D.. failed'
+  end if
+
+  if (.not. all(keras_layers(3) % pool_size == [2, 2])) then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN second layer pool_size is expected.. failed'
+  end if
+
+  if (.not. all(keras_layers(3) % strides == [2, 2])) then
+    ok = .false.
+    write(stderr, '(a)') &
+      'Keras CNN second layer strides are expected.. failed'
+  end if
+
   if (ok) then
     print '(a)', 'test_keras_read_model: All tests passed.'
   else