Merge lp:~fluidity-core/fluidity/adjoint into lp:fluidity

[Patrick Farrell]

Overview of the major changes:
------------------------------

* The burgers equation adjoint now works, including nonlinear terms (and so their derivatives). The differentiation of the CG discretisation of the advection term is achieved using tapenade.

* The development of a program (scripts/optimality.py) to solve optimisation problems with fluidity/fluidity-like models.

* More (and more rigorous) test cases. The fundamental tool for verifying the discrete adjoint is the derivative test, which is implemented in optimality.py. Suppose you can compute J(m), and have a candidate for grad(J)(m). Then

|| J(m + dm) - J(m) || should converge at O(||dm||), and
|| J(m + dm) - J(m) - grad(J)(m) . dm || should converge at O(||dm||^2).

* Much stronger tools for debugging adjoint model development, including
  - forward replay verification
  - transpose action correctness
  - nonlinear derivative consistency
  - the derivative test in optimality.py

[lp:~fluidity-core/fluidity/adjoint]

3589. By Patrick Farrell on 2011-07-12

Removing shallow_water_mms_sphere: it never worked.

Simon: I checked out the commit which added this, and the
test is broken there. If you want this test, revive it,
fix it, and add it back.

3590. By Patrick Farrell on 2011-07-12

In the shallow water model, we now write_diagnostics before the timeloop
begins. (Simon: why? It was r3521. Something to do with the optimiser?)
So we have to change the test to look at stat[1] instead of stat[0] so that
it is looking at the same thing.

3591. By Patrick Farrell on 2011-07-12

Loosen some tolerances slightly, and activate two new assertions that had been
broken

3592. By Patrick Farrell on 2011-07-12

Try to get some useful information out of this test; it is failing on
the buildbot, but works when I run it by hand (on my workstation,
and on the buildbot builder)

3593. By Patrick Farrell on 2011-07-12

Implement David's suggestions during the code review.

[Patrick Farrell (pefarrell)]

For posterity, David's review requests were:

- get rid of the boundary_condition control from the schema -- the code currently doesn't honour it.
- remove the commented-out comparison code in advection_action in the burgers adjoint callbacks
- Makefile targets for ADing simple advection (but I don't think I can do this one)
- Rather than checking if gfortran version is >= 4.5 myself, use the handy env variable $GFORTRAN_4_5_OR_NEWER that David made
- the testharness should filter out long tests, by default
- have the testharness print out the criteria it's using to search for tests.

These were fixed in r3593.

[Patrick Farrell (pefarrell)]

All tests pass with --with-adjoint=no:

http://buildbot-ocean.ese.ic.ac.uk:8080/builders/x86_64-branch-adjoint-gcc4/builds/8

(Ignore the colour and look at the stdio; the buildbot is interpreting it incorrectly.)

Just waiting on the --with-adjoint=yes case.

[David Ham (david-ham)]

review: Approve

[Patrick Farrell (pefarrell)]

The --with-adjoint case also passes:

http://buildbot-ocean.ese.ic.ac.uk:8080/builders/adjoint/builds/9/steps/shell_9/logs/stdio

The only exception is the channel-flow-dg case, which is known broken in the trunk.

=== modified file 'Makefile.in'
--- Makefile.in	2011-06-27 14:01:01 +0000
+++ Makefile.in	2011-07-12 15:56:32 +0000
@@ -424,7 +424,7 @@
 	@echo "    CLEAN tests"
 	@cd tests; ../tools/testharness.py --clean >/dev/null
 
-test: fltools bin/@FLUIDITY@ bin/shallow_water serialtest
+test: fltools bin/@FLUIDITY@ bin/shallow_water serialtest bin/burgers_equation
 
 serialtest:
 ifeq (@MBA2D@,yes)
@@ -441,7 +441,7 @@
         endif
 endif
 
-mediumtest:
+mediumtest: bin/burgers_equation
 	@cd tests; ../tools/testharness.py -l medium $(EXCLUDE_TAGS) -n $(THREADS)
 
 mediumzoltantest:
 
=== modified file 'aclocal.m4'
--- aclocal.m4	2011-02-01 18:21:13 +0000
+++ aclocal.m4	2011-07-12 15:56:32 +0000
@@ -1030,10 +1030,11 @@
 	[adjoint="$withval"],
     [])
 
+bakLIBS=$LIBS
 tmpLIBS=$LIBS
 tmpCPPFLAGS=$CPPFLAGS
-if test $adjoint != no; then
-  if test $adjoint != yes; then
+if test "$adjoint" != "no"; then
+  if test "$adjoint" != "yes"; then
     adjoint_LIBS_PATH="$adjoint/lib"
     adjoint_INCLUDES_PATH="$adjoint/include"
     # Ensure the comiler finds the library...
@@ -1051,10 +1052,9 @@
 AC_CHECK_LIB(
 	[adjoint],
 	[adj_register_equation],
-	[AC_DEFINE(HAVE_ADJOINT,1,[Define if you have libadjoint.])],
-	[AC_MSG_ERROR( [Could not link in libadjoint ... exiting] )] )
+	[AC_DEFINE(HAVE_ADJOINT,1,[Define if you have libadjoint.])HAVE_ADJOINT=yes],
+	[AC_MSG_WARN( [Could not link in libadjoint ... ] );HAVE_ADJOINT=no;LIBS=$bakLIBS] )
 # Save variables...
 AC_LANG_RESTORE
 
-echo $LIBS
 ])dnl ACX_adjoint
 
=== added file 'adjoint/Adjoint_Controls.F90'
--- adjoint/Adjoint_Controls.F90	1970-01-01 00:00:00 +0000
+++ adjoint/Adjoint_Controls.F90	2011-07-12 15:56:32 +0000
@@ -0,0 +1,611 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+
+#include "fdebug.h"
+
+module adjoint_controls
+  use fields
+  use global_parameters, only : PYTHON_FUNC_LEN, OPTION_PATH_LEN
+  use boundary_conditions
+  use spud
+  use state_module
+  use python_state
+
+    implicit none
+
+    private
+
+    public :: adjoint_load_controls, adjoint_write_controls, adjoint_write_control_derivatives, allocate_and_insert_control_derivative_fields
+
+    contains 
+
+    ! Retrieves the control details from control number control_no in the optin tree. The outputs are:
+    subroutine get_control_details(states, timestep, control_no, control_name, state_id, field_name, field_type, &
+                                 & active, have_lb, lb_state_id, lb_material_phase_name, lb_field_name, have_ub, &
+                                 & ub_state_id, ub_material_phase_name, ub_field_name)
+
+       ! Inputs:
+      type(state_type), dimension(:), intent(in) :: states
+      integer, intent(in) :: control_no, timestep ! control_no: The index of the control in the option tree
+                                                  ! timestep: The timestep
+
+      ! Outputs:
+      character(len=OPTION_PATH_LEN), intent(out) :: control_name, field_type, field_name ! control_name: the name of the control 
+                                                                                          ! field_name: the name of the associated field
+                                                                                          ! field_type: the type of the associated field
+      integer, intent(out) :: state_id ! the state number in which the associated control field exists
+      logical, intent(out) :: active ! false if the control is active in this timestep (e.g. InitialConditions for timesteps >0)
+      logical, intent(out), optional :: have_ub, have_lb
+      integer, intent(out), optional :: lb_state_id, ub_state_id
+      character(len=OPTION_PATH_LEN), intent(out), optional :: lb_material_phase_name, ub_material_phase_name
+      character(len=OPTION_PATH_LEN), intent(out), optional :: ub_field_name, lb_field_name
+
+      integer :: s_idx
+      character(len=OPTION_PATH_LEN) :: control_type, material_phase_name, name 
+      type(scalar_field), pointer :: sfield => null(), sfield2 => null()
+      type(vector_field), pointer :: vfield => null(), vfield2 => null()
+      type(tensor_field), pointer :: tfield => null(), tfield2 => null()
+
+      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/name", control_name)
+      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/type/name", control_type)
+      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/type::" // trim(control_type) // "/field_name", name)
+
+      s_idx = scan(trim(name), "::")
+      if (s_idx == 0) then 
+        FLExit("The control " // trim(control_name) // " uses an invalid field_name. It should be of the form Materialphase::Field")
+      end if
+      material_phase_name = name(1:s_idx - 1)
+      field_name = name(s_idx + 2:len_trim(name))
+
+      active = .true.
+
+      ! Find state associated with the material phase
+      do state_id = 1, size(states) 
+        if (trim(states(state_id)%name) == trim(material_phase_name)) then
+          exit
+        end if
+      end do
+      if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
+        FLExit("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_name) // ".")
+      end if
+      
+      select case (trim(control_type))
+        case ("initial_condition")
+          if (timestep > 0) then
+            active = .false.
+            field_type = ''
+            field_name = ''
+            return
+          endif
+          assert(timestep == 0)
+          if (has_scalar_field(states(state_id), field_name)) then
+            field_type = "scalar"
+          elseif (has_vector_field(states(state_id), field_name)) then
+            field_type = "vector"
+          elseif (has_tensor_field(states(state_id), field_name)) then
+            field_type = "tensor"
+          else
+            ewrite(0, *) "The control field " // trim(field_name) // " specified in control " // trim(control_name) // " is not a field in the state."
+            ewrite(0, *) "The current state is: "
+            call print_state(states(state_id))
+            FLExit("Check your control's field settings.")
+          end if
+
+        case ("source_term")
+          if (timestep == 0) then
+            active = .false.
+            field_type = ''
+            field_name = ''
+            return
+          end if
+
+          if (has_scalar_field(states(state_id), field_name)) then
+            field_type = "scalar"
+          elseif (has_vector_field(states(state_id), field_name)) then
+            field_type = "vector"
+          elseif (has_tensor_field(states(state_id), field_name)) then
+            field_type = "tensor"
+          else
+            ewrite(0, *) "The control field " // trim(name) // " specified in control " // trim(control_name) // & 
+                       & " is not a field in state " // trim(material_phase_name) // "."
+            ewrite(0, *) "The current state is: "
+            call print_state(states(state_id))
+            FLExit("Check your control's field settings.")
+          end if
+
+        case ("boundary_condition")
+          FLAbort("Boundary condition control not implemented yet.")
+      end select
+
+      !!!! Check bounds settings !!!!
+      ! We only execute this when the optional files are present
+      if (.not. present(have_lb) .or. .not. present(have_ub)) then
+        return
+      endif
+
+      have_ub = have_option("/adjoint/controls/control[0]/bounds/upper_bound")
+      have_ub = have_option("/adjoint/controls/control[" // int2str(control_no) //"]/bounds/upper_bound")
+      have_lb = have_option("/adjoint/controls/control[" // int2str(control_no) //"]/bounds/lower_bound")
+
+      ! Find material phase for the lower bound field
+      if (have_lb) then 
+        call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/bounds/lower_bound/field_name", name)
+        if (have_lb) then
+          s_idx = scan(trim(name), "::")
+          if (s_idx == 0) then 
+            FLExit("The control lower bound for control " // trim(control_name) // " uses an invalid field_name. It should be of the form Materialphase::Field")
+          end if
+          lb_material_phase_name = name(1:s_idx - 1)
+          lb_field_name = name(s_idx + 2:len_trim(name))
+          ! Find state associated with the lower bound material phase
+          do lb_state_id = 1, size(states) 
+            if (trim(states(lb_state_id)%name) == trim(lb_material_phase_name)) then
+              exit
+            end if
+          end do
+          if (.not. trim(states(lb_state_id)%name) == trim(lb_material_phase_name)) then
+            FLExit("Could not find state " // trim(lb_material_phase_name) // " as specified in the lower bound of control " // trim(control_name) // ".")
+          end if
+        end if
+      end if
+
+      ! Find material phase for the upper bound field
+      if (have_ub) then
+        call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/bounds/upper_bound/field_name", name)
+        if (have_ub) then
+          s_idx = scan(trim(name), "::")
+          if (s_idx == 0) then 
+            FLExit("The control upper bound for control " // trim(control_name) // " uses an invalid field_name. It should be of the form Materialphase::Field")
+          end if
+          ub_material_phase_name = name(1:s_idx - 1)
+          ub_field_name = name(s_idx + 2:len_trim(name))
+          ! Find state associated with the upper bound material phase
+          do ub_state_id = 1, size(states) 
+            if (trim(states(ub_state_id)%name) == trim(ub_material_phase_name)) then
+              exit
+            end if
+          end do
+          if (.not. trim(states(ub_state_id)%name) == trim(ub_material_phase_name)) then
+            FLExit("Could not find state " // trim(ub_material_phase_name) // " as specified in the upper bound of control " // trim(control_name) // ".")
+          end if
+        end if
+      end if
+
+      ! Check that the lower/upper bound fields indeed exist and have the same mesh than the control
+      select case (trim(field_type))
+
+        case ("scalar")
+          ! Upper bound
+          ! Check that the specified field exist in the state
+          if (have_ub) then
+            if (.not. has_scalar_field(states(ub_state_id), ub_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(ub_field_name), " for the upper bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the upper bound use the same mesh
+            sfield => extract_scalar_field(states(state_id), field_name)
+            sfield2 => extract_scalar_field(states(ub_state_id), ub_field_name)
+            if (.not. sfield%mesh == sfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The upper bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+          ! Lower bound
+          if (have_lb) then
+            if (.not. has_scalar_field(states(state_id), lb_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(lb_field_name), " for the lower bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the lower bound use the same mesh
+            sfield => extract_scalar_field(states(state_id), field_name)
+            sfield2 => extract_scalar_field(states(lb_state_id), lb_field_name)
+            if (.not. sfield%mesh == sfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The lower bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+        case ("vector")
+          ! Upper bound
+          ! Check that the specified field exist in the state
+          if (have_ub) then
+            if (.not. has_vector_field(states(ub_state_id), ub_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(ub_field_name), " for the upper bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the upper bound use the same mesh
+            vfield => extract_vector_field(states(state_id), field_name)
+            vfield2 => extract_vector_field(states(ub_state_id), ub_field_name)
+            if (.not. vfield%mesh == vfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The upper bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+          ! Lower bound
+          if (have_lb) then
+            if (.not. has_vector_field(states(state_id), lb_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(lb_field_name), " for the lower bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the lower bound use the same mesh
+            vfield => extract_vector_field(states(state_id), field_name)
+            vfield2 => extract_vector_field(states(lb_state_id), lb_field_name)
+            if (.not. vfield%mesh == vfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The lower bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+        case ("tensor")
+          ! Upper bound
+          ! Check that the specified field exist in the state
+          if (have_ub) then
+            if (.not. has_tensor_field(states(ub_state_id), ub_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(ub_field_name), " for the upper bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the upper bound use the same mesh
+            tfield => extract_tensor_field(states(state_id), field_name)
+            tfield2 => extract_tensor_field(states(ub_state_id), ub_field_name)
+            if (.not. tfield%mesh == tfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The lower bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+          ! Lower bound
+          if (have_lb) then
+            if (.not. has_tensor_field(states(state_id), lb_field_name)) then
+              ewrite(0, *) "Can not find the specified field ", trim(lb_field_name), " for the lower bound of control ", trim(control_name)
+              FLExit("Check your control's field settings.")
+            end if
+
+            ! Check that the control and the lower bound use the same mesh
+            tfield => extract_tensor_field(states(state_id), field_name)
+            tfield2 => extract_tensor_field(states(lb_state_id), lb_field_name)
+            if (.not. sfield%mesh == sfield2%mesh) then
+              ewrite(0, *) "Error in control ", trim(control_name), " definition. The lower bound and the control's mesh have to be the same."
+              FLExit("Check your control's field settings.")
+            end if
+          end if
+
+      end select
+    end subroutine get_control_details
+
+    ! Loop through all controls in the option file and store them in files.
+    subroutine adjoint_write_controls(timestep, dt, states)
+      integer, intent(in) :: timestep
+      real, intent(in) :: dt
+      type(state_type), dimension(:), intent(in) :: states
+#ifdef HAVE_ADJOINT
+      integer :: nb_controls, i, state_id
+      logical :: active
+      character(len=OPTION_PATH_LEN) :: simulation_name, field_name, field_type, control_name 
+      logical :: have_ub, have_lb
+      integer :: lb_state_id, ub_state_id
+      character(len=OPTION_PATH_LEN) :: lb_material_phase_name, ub_material_phase_name
+      character(len=OPTION_PATH_LEN) :: ub_field_name, lb_field_name
+
+      if (.not. have_option("/adjoint/controls")) then
+        return
+      end if
+
+      nb_controls = option_count("/adjoint/controls/control")
+      call get_option("/simulation_name", simulation_name)
+      ! Now loop over the controls and save them to disk
+      do i = 0, nb_controls-1
+        call get_control_details(states, timestep, i, control_name, state_id, field_name, field_type, active, have_lb, lb_state_id, lb_material_phase_name, lb_field_name, have_ub, ub_state_id, ub_material_phase_name, ub_field_name)
+        if (active) then
+
+          ! Save the control parameter to disk
+          call python_reset()
+          call python_add_state(states(state_id))
+          call python_run_string("import pickle;" // &
+                              &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl';" // & 
+                              &  "f = open(fname, 'wb');" // &
+                              &  "field = state." // trim(field_type) // "_fields['" // trim(field_name) // "'];" // &
+                              &  "pickle.dump(field.val[:], f);" // &
+                              &  "f.close();")
+
+          ! And its bounds
+          call python_reset()
+          if (have_lb) then
+            call python_add_state(states(lb_state_id))
+            call python_run_string("import pickle;" // &
+                                &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_name) // "_lower_bound_" // int2str(timestep) // ".pkl';" // & 
+                                &  "f = open(fname, 'wb');" // &
+                                &  "field = state." // trim(field_type) // "_fields['" // trim(lb_field_name) // "'];" // &
+                                &  "pickle.dump(field.val[:], f);" // &
+                                &  "f.close();")
+            call python_reset()
+          end if
+
+          if (have_ub) then
+            call python_add_state(states(ub_state_id))
+            call python_run_string("import pickle;" // &
+                                &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_name) // "_upper_bound_" // int2str(timestep) // ".pkl';" // & 
+                                &  "f = open(fname, 'wb');" // &
+                                &  "field = state." // trim(field_type) // "_fields['" // trim(ub_field_name) // "'];" // &
+                                &  "pickle.dump(field.val[:], f);" // &
+                                &  "f.close();")
+            call python_reset()
+          endif
+
+        end if
+      end do
+#endif
+    end subroutine adjoint_write_controls
+
+    ! The counterpart of adjoint_write_controls: Loop through all controls in the option file and fill the fields from the control files.
+    subroutine adjoint_load_controls(timestep, dt, states)
+      integer, intent(in) :: timestep
+      real, intent(in) :: dt
+      type(state_type), dimension(:), intent(in) :: states
+#ifdef HAVE_ADJOINT
+      integer :: nb_controls, i, state_id, s_idx
+      logical :: active, file_exists
+      character(len=OPTION_PATH_LEN) :: field_name, field_type, control_name, simulation_name, control_type
+
+      ! Do not read anything if we are supposed to write the controls
+      if (.not. have_option("/adjoint/controls/load_controls")) then
+        return 
+      end if
+
+      call get_option("/simulation_name", simulation_name)
+      nb_controls = option_count("/adjoint/controls/control")
+
+      ! Now loop over the controls, read their controls files and fill the associated fields
+      do i = 0, nb_controls-1
+        call get_control_details(states, timestep, i, control_name, state_id, field_name, field_type, active)
+        call get_option("/adjoint/controls/control[" // int2str(i) // "]/type/name", control_type)
+
+        ! We load all OTHER controls, then apply our strong boundary conditions, then load BC controls
+        if (trim(control_type) == "boundary_condition") cycle
+
+        if (active) then
+          ! Check that the control parameter file extists
+          inquire(file="control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl", exist=file_exists)
+          if (.not. file_exists) then
+            ewrite (0, *) "File 'control_" // trim(simulation_name) // "_" // trim(control_name) // &
+                        & "_" // int2str(timestep) // ".pkl' does not exist."
+            FLExit("Error while loading control files. You might want to remove '/adjoint/controls/load_controls'?")
+          end if
+          ! Read the control parameter from disk
+          call python_reset()
+          call python_add_state(states(state_id))
+          call python_run_string("import pickle;" // &
+                              &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl';" // & 
+                              &  "f = open(fname, 'rb');" // &
+                              &  "field = state." // trim(field_type) // "_fields['" // trim(field_name) // "'];" // &
+                              &  "field.val[:] = pickle.load(f);" // &
+                              &  "f.close();")
+          call python_reset()
+        end if
+      end do
+
+      ! Now apply the BCs, because the BCs should overwrite (say) an initial_condition control
+      call set_dirichlet_consistent(states)
+
+      do i = 0, nb_controls-1
+        call get_control_details(states, timestep, i, control_name, state_id, field_name, field_type, active)
+        call get_option("/adjoint/controls/control[" // int2str(i) // "]/type/name", control_type)
+
+        ! We load all OTHER controls, then apply our strong boundary conditions, then load BC controls
+        if (trim(control_type) /= "boundary_condition") cycle
+
+        if (active) then
+          ! Check that the control parameter file extists
+          inquire(file="control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl", exist=file_exists)
+          if (.not. file_exists) then
+            ewrite (0, *) "File 'control_" // trim(simulation_name) // "_" // trim(control_name) // &
+                        & "_" // int2str(timestep) // ".pkl' does not exist."
+            FLExit("Error while loading control files. You might want to remove '/adjoint/controls/load_controls'?")
+          end if
+          ! Read the control parameter from disk
+          call python_reset()
+          call python_add_state(states(state_id))
+          call python_run_string("import pickle;" // &
+                              &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl';" // & 
+                              &  "f = open(fname, 'rb');" // &
+                              &  "field = state." // trim(field_type) // "_fields['" // trim(field_name) // "'];" // &
+                              &  "field.val[:] = pickle.load(f);" // &
+                              &  "f.close();")
+          call python_reset()
+        end if
+      end do
+#endif
+    end subroutine adjoint_load_controls
+
+    ! Extracts the total derivatives of the functional and saves them to disk.
+    subroutine adjoint_write_control_derivatives(states, timestep, functional_name)
+      type(state_type), dimension(:), intent(inout) :: states
+      integer, intent(in) :: timestep
+      character(len=*), intent(in) :: functional_name
+      character(len=OPTION_PATH_LEN) :: field_name, control_type, control_deriv_name, field_deriv_name, name, material_phase_name, simulation_name, field_type
+      integer :: nb_controls
+      integer :: i, state_id, s_idx
+      type(scalar_field), pointer :: sfield => null()
+      type(vector_field), pointer :: vfield => null()
+      type(tensor_field), pointer :: tfield => null()
+     
+      if (.not. have_option("/adjoint/controls")) then
+        return
+      end if
+      call get_option("/simulation_name", simulation_name)
+      ! Remove the suffix _adjoint from the simulation_name
+
+      nb_controls = option_count("/adjoint/controls/control")
+      do i = 0, nb_controls-1
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/name", control_deriv_name)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type/name", control_type)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type::" // trim(control_type) // "/field_name", name)
+        s_idx = scan(trim(name), "::")
+        if (s_idx == 0) then 
+          FLExit("The control " // trim(control_deriv_name) // " uses an invalid field_name. It should be of the form Materialphase::Fieldname")
+        end if
+        material_phase_name = name(1:s_idx - 1)
+        field_name = name(s_idx + 2:len_trim(name))
+        ! Find state associated with the material phase
+        do state_id = 1, size(states) 
+          if (trim(states(state_id)%name) == trim(material_phase_name)) then
+            exit
+          end if
+        end do
+        ! Make sure we found the state
+        if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
+          FLExit("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_deriv_name) // ".")
+        end if
+
+        field_name = "Adjoint" // trim(field_name)
+        control_deriv_name = trim(control_deriv_name) // "_TotalDerivative"
+        if (trim(control_type) == "initial_condition" .or. trim(control_type) == "source_term") then
+
+          if (trim(control_type) == "initial_condition" .and. timestep > 0) then
+            cycle
+          end if
+
+          if (trim(control_type) == "source_term" .and. timestep == 0) then
+            cycle
+          end if
+
+          field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+          if (has_scalar_field(states(state_id), field_name)) then
+            field_type = "scalar"
+          elseif (has_vector_field(states(state_id), field_name)) then
+            field_type = "vector"
+          elseif (has_tensor_field(states(state_id), field_name)) then
+            field_type = "tensor"
+          else
+            ewrite(0, *) "The control field " // trim(field_deriv_name) // " specified in control " // trim(control_deriv_name) // " is not a field in the state."
+            ewrite(0, *) "The current state is: "
+            call print_state(states(state_id))
+            FLExit("Check your control's field settings.")
+          end if
+        elseif (trim(control_type) == "boundary_condition") then
+          FLAbort("Boundary condition control not implemented yet.")
+        end if
+        ! Save the control parameter to disk
+        call python_reset()
+        call python_add_state(states(state_id))
+        call python_run_string("import pickle;" // &
+                            &  "fname = 'control_" // trim(simulation_name) // "_" // trim(control_deriv_name) // "_" // int2str(timestep) // ".pkl';" // & 
+                            &  "f = open(fname, 'wb');" // &
+                            &  "field = state." // trim(field_type) // "_fields['" // trim(field_deriv_name) // "'];" // &
+                            &  "pickle.dump(field.val[:], f);" // &
+                            &  "f.close();")
+        call python_reset()
+      end do
+    end subroutine adjoint_write_control_derivatives
+
+    subroutine allocate_and_insert_control_derivative_fields(states)
+      type(state_type), dimension(:), intent(inout) :: states
+      integer :: nb_controls, nb_functionals, i, state_id, functional, s_idx
+      character(len=OPTION_PATH_LEN) :: field_name, control_type, control_deriv_name, functional_name, field_deriv_name, material_phase_name, name
+      type(scalar_field), pointer :: sfield => null()
+      type(vector_field), pointer :: vfield => null()
+      type(tensor_field), pointer :: tfield => null()
+      type(scalar_field) :: sfield_deriv
+      type(vector_field) :: vfield_deriv
+      type(tensor_field) :: tfield_deriv
+
+      nb_controls = option_count("/adjoint/controls/control")
+      nb_functionals = option_count("/adjoint/functional")
+      ! Loop over the functionals
+      do functional = 0, nb_functionals-1
+        call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+        ! Now loop over the controls 
+        do i = 0, nb_controls-1
+          call get_option("/adjoint/controls/control[" // int2str(i) //"]/name", control_deriv_name)
+          call get_option("/adjoint/controls/control[" // int2str(i) //"]/type/name", control_type)
+          call get_option("/adjoint/controls/control[" // int2str(i) //"]/type::" // trim(control_type) // "/field_name", name)
+          s_idx = scan(trim(name), "::")
+          if (s_idx == 0) then 
+            FLExit("The control " // trim(control_deriv_name) // " uses an invalid field_name. It should be of the form Materialphase::Fieldname")
+          end if
+          material_phase_name = name(1:s_idx - 1)
+          field_name = name(s_idx + 2:len_trim(name))
+          ! Find state associated with the material phase
+          do state_id = 1, size(states) 
+            if (trim(states(state_id)%name) == trim(material_phase_name)) then
+              exit
+            end if
+          end do
+          ! Make sure we found the state
+          if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
+            FLExit("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_deriv_name) // ".")
+          end if
+
+          field_name = "Adjoint" // trim(field_name)
+          control_deriv_name = trim(control_deriv_name) // "_TotalDerivative"
+
+          if (trim(control_type) == "initial_condition" .or. trim(control_type) == "source_term") then
+            field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+
+            if (has_scalar_field(states(state_id), field_name)) then
+              sfield => extract_scalar_field(states(state_id), field_name)
+              call allocate(sfield_deriv, sfield%mesh, field_deriv_name)
+              call zero(sfield_deriv)
+              call insert(states(state_id), sfield_deriv, field_deriv_name)
+              call deallocate(sfield_deriv)
+
+            elseif (has_vector_field(states(state_id), field_name)) then
+              vfield => extract_vector_field(states(state_id), field_name)
+              call allocate(vfield_deriv, vfield%dim, vfield%mesh, field_deriv_name)
+              call zero(vfield_deriv)
+              call insert(states(state_id), vfield_deriv, field_deriv_name)
+              call deallocate(vfield_deriv)
+
+            elseif (has_tensor_field(states(state_id), field_name)) then
+              tfield => extract_tensor_field(states(state_id), field_name)
+              call allocate(tfield_deriv, tfield%mesh, field_deriv_name, tfield%field_type, tfield%dim)
+              call zero(tfield_deriv)
+              call insert(states(state_id), tfield_deriv, field_deriv_name)
+              call deallocate(tfield_deriv)
+
+            else
+              ewrite(0, *) "The control field " // trim(field_name) // " specified in control " // trim(control_deriv_name) // " is not a field in the state."
+              ewrite(0, *) "The current state is: "
+              call print_state(states(state_id))
+              FLExit("Check your control's field settings.")
+            end if
+
+          elseif (trim(control_type) == "boundary_condition") then
+            FLAbort("Boundary condition control not implemented yet.")
+          end if
+
+        end do
+      end do
+    end subroutine allocate_and_insert_control_derivative_fields
+
+end module adjoint_controls
 
=== modified file 'adjoint/Adjoint_Functional_Evaluation.F90'
--- adjoint/Adjoint_Functional_Evaluation.F90	2011-05-20 17:42:50 +0000
+++ adjoint/Adjoint_Functional_Evaluation.F90	2011-07-12 15:56:32 +0000
@@ -48,7 +48,7 @@
   implicit none
 
   private
-  public :: libadjoint_functional_derivative, adj_record_anything_necessary
+  public :: libadjoint_evaluate_functional, libadjoint_functional_derivative, adj_record_anything_necessary
 
   contains
 
@@ -250,32 +250,6 @@
     end if
 
     if (has_func) then
-      if (.not. functional_computed) then
-        ierr = adj_timestep_get_times(adjointer, timelevel, start_time, end_time)
-        call adj_chkierr(ierr)
-        assert(start_time < end_time)
-
-        write(buffer,*) timelevel + 1
-        call python_run_string("n = " // trim(buffer))
-
-        write(buffer,*) start_time
-        call python_run_string("time = " // trim(buffer))
-
-        write(buffer, *) abs(start_time - end_time)
-        call python_run_string("dt = " // trim(buffer))
-
-        call python_run_string(trim(code_func))
-        J = python_fetch_real("J")
-
-        ! So we've computed the component of the functional associated with this timestep.
-        ! We also want to sum them all up ...
-
-        call set_diagnostic(name=trim(functional_name_f) // "_component", statistic="value", value=(/J/))
-        fn_value => get_diagnostic(name=trim(functional_name_f), statistic="value")
-        J = J + fn_value(1)
-        call set_diagnostic(name=trim(functional_name_f), statistic="value", value=(/J/))
-        functional_computed = .true.
-      end if
 
       if (.not. has_deriv) then
         check_uncertainties = "try: " // achar(10) // &
@@ -341,6 +315,111 @@
 
   end subroutine libadjoint_functional_derivative
 
+  subroutine libadjoint_evaluate_functional(adjointer, timelevel, ndepends, dependencies, values, functional_name_c, output) bind(c) 
+    use iso_c_binding
+    use libadjoint_data_structures
+    type(adj_adjointer), intent(in) :: adjointer
+    integer(kind=c_int), intent(in), value :: timelevel
+    integer(kind=c_int), intent(in), value :: ndepends
+    type(adj_variable), dimension(ndepends), intent(in) :: dependencies
+    type(adj_vector), dimension(ndepends), intent(in) :: values
+    character(kind=c_char), dimension(ADJ_NAME_LEN), intent(in) :: functional_name_c
+    adj_scalar_f, intent(out) :: output
+
+    type(state_type), dimension(:,:), pointer :: states ! material_phases x timelevels
+    character(len=ADJ_NAME_LEN) :: functional_name_f
+    logical :: has_func
+    character(len=PYTHON_FUNC_LEN) :: code_func
+    integer :: ierr, max_timelevel, tmp_timelevel, min_timelevel, nmaterial_phases
+    character(len=OPTION_PATH_LEN) :: material_phase_name
+    integer :: no_timesteps, timestep, i
+    real :: start_time, end_time
+    character(len = 30) :: buffer, buffer2
+
+    call python_reset
+
+    functional_name_f = ' '
+    i = 1
+    do while (i <= ADJ_NAME_LEN .and. functional_name_c(i) /= c_null_char)
+      functional_name_f(i:i) = functional_name_c(i)
+      i = i + 1
+    end do
+
+    ! Fetch the python code the user has helpfully supplied.
+    has_func  = have_option("/adjoint/functional::" // trim(functional_name_f) // "/functional_value/algorithm")
+    if (.not. has_func) then
+      FLAbort("You need to supply code for the functional to evaluate it!")
+    endif
+
+    if (has_func) then
+      call get_option("/adjoint/functional::" // trim(functional_name_f) // "/functional_value/algorithm", code_func)
+    endif
+
+    if (ndepends==0) then
+      max_timelevel = 0
+      min_timelevel = 0
+    else
+      ! Convert from the list of adj_values/adj_vectors to actual states
+      ierr = adj_variable_get_timestep(dependencies(1), max_timelevel)
+      min_timelevel = max_timelevel
+      do i=1,ndepends
+        ierr = adj_variable_get_timestep(dependencies(i), tmp_timelevel)
+        max_timelevel = max(max_timelevel, tmp_timelevel)
+        min_timelevel = min(min_timelevel, tmp_timelevel)
+      end do
+    end if
+
+    nmaterial_phases = option_count("/material_phase")
+    allocate(states(nmaterial_phases, min_timelevel:max_timelevel))
+
+    ! Set the names of the material phases, so that we can look them up and know
+    ! which material_phase to put stuff into in convert_adj_vectors_to_states
+    do i=0,nmaterial_phases-1
+      call get_option("/material_phase[" // int2str(i) // "]/name", material_phase_name)
+      states(i+1,:)%name = trim(material_phase_name)
+    end do
+    call convert_adj_vectors_to_states(dependencies, values, states)
+    call python_add_states_time(states)
+
+    ! Also set up some useful variables for the user to use
+    ierr = adj_timestep_count(adjointer, no_timesteps)
+    call adj_chkierr(ierr)
+    write(buffer, *) no_timesteps
+    call python_run_string("times = [0] * " // trim(buffer))
+    do timestep=0,no_timesteps-2
+      ierr = adj_timestep_get_times(adjointer, timestep, start_time, end_time)
+      call adj_chkierr(ierr)
+      write(buffer, '(i0)') timestep
+      write(buffer2, *) start_time
+      call python_run_string("times[" // trim(buffer) // "] = " // trim(buffer2))
+      write(buffer, '(i0)') timestep+1
+      write(buffer2, *) end_time
+      call python_run_string("times[" // trim(buffer) // "] = " // trim(buffer2))
+    end do
+
+    ierr = adj_timestep_get_times(adjointer, timelevel, start_time, end_time)
+    call adj_chkierr(ierr)
+    assert(start_time < end_time)
+
+    write(buffer,*) timelevel + 1
+    call python_run_string("n = " // trim(buffer))
+
+    write(buffer,*) start_time
+    call python_run_string("time = " // trim(buffer))
+
+    write(buffer, *) abs(start_time - end_time)
+    call python_run_string("dt = " // trim(buffer))
+
+    call python_run_string(trim(code_func))
+    output = python_fetch_real("J")
+    
+    if (max_timelevel >= 0) then
+      call deallocate(states)
+      deallocate(states)
+    endif
+
+  end subroutine libadjoint_evaluate_functional
+
   subroutine convert_adj_vectors_to_states(dependencies, values, states)
     type(adj_variable), dimension(:), intent(in) :: dependencies
     type(adj_vector), dimension(:), intent(in) :: values
@@ -432,7 +511,7 @@
     call get_option("/timestepping/finish_time", finish_time)
 
     call get_option("/adjoint/functional::" // trim(functional) // "/functional_dependencies/algorithm", buf)
-    call adj_variables_from_python(buf, current_time, finish_time, python_timestep, vars)
+    call adj_variables_from_python(adjointer, buf, current_time, finish_time, python_timestep, vars)
 
     variable_loop: do j=1,size(vars)
       ierr = adj_variable_get_timestep(vars(j), var_timestep)
@@ -522,5 +601,6 @@
 
     deallocate(vars)
   end subroutine adj_record_anything_necessary
+
 #endif
 end module adjoint_functional_evaluation
 
=== modified file 'adjoint/Adjoint_Global_Variables.F90'
--- adjoint/Adjoint_Global_Variables.F90	2011-04-01 17:40:34 +0000
+++ adjoint/Adjoint_Global_Variables.F90	2011-07-12 15:56:32 +0000
@@ -35,6 +35,5 @@
 
   type(adj_dictionary) :: adj_path_lookup, adj_solver_path_lookup
   type(adj_adjointer)  :: adjointer
-  logical :: functional_computed = .false.
 #endif
 end module adjoint_global_variables
 
=== modified file 'adjoint/Adjoint_Main_Loop.F90'
--- adjoint/Adjoint_Main_Loop.F90	2011-05-11 17:15:33 +0000
+++ adjoint/Adjoint_Main_Loop.F90	2011-07-12 15:56:32 +0000
@@ -39,6 +39,7 @@
     use global_parameters, only: OPTION_PATH_LEN, running_adjoint
     use adjoint_global_variables
     use adjoint_functional_evaluation
+    use adjoint_controls
     use populate_state_module
     use signal_vars
     use mangle_options_tree
@@ -53,7 +54,6 @@
     implicit none
 
     private
-    public :: adjoint_main_loop_register_diagnostic
 #ifdef HAVE_ADJOINT
     public :: compute_adjoint
 #endif
@@ -61,9 +61,28 @@
     contains
 
 #ifdef HAVE_ADJOINT
-    subroutine compute_adjoint(state, dump_no)
+    ! Computes the adjoint equation
+    ! The optional adjoint timestep callback is called for each functional after every timestep in the adjoint loop.
+    ! It is commonly used to compute the diangostics from the adjoint solution, e.g. the total derivative of the functional.
+    subroutine compute_adjoint(state, dump_no, adjoint_timestep_callback, adjoint_timestep_callback_data)
+      use iso_c_binding, only: c_ptr
       type(state_type), dimension(:), intent(inout) :: state
       integer, intent(inout) :: dump_no
+      optional :: adjoint_timestep_callback
+      ! Data is a void pointer used to pass variables into the callback
+      type(c_ptr), intent(in), optional :: adjoint_timestep_callback_data
+      interface 
+        subroutine adjoint_timestep_callback(state, timestep, functional_name, data)
+          use iso_c_binding, only: c_ptr
+          use global_parameters, only: OPTION_PATH_LEN
+          use state_module
+          type(state_type), dimension(:), intent(inout) :: state
+          integer, intent(in) :: timestep
+          character(len=OPTION_PATH_LEN), intent(in) :: functional_name
+          ! Data is a void pointer used to pass variables into the callback
+          type(c_ptr), intent(in) :: data
+        end subroutine adjoint_timestep_callback
+      end interface
 
       type(adj_vector) :: rhs
       type(adj_vector) :: soln
@@ -80,10 +99,12 @@
       integer :: end_timestep, start_timestep, no_timesteps, timestep
       real :: start_time, end_time
 
-      character(len=OPTION_PATH_LEN) :: simulation_base_name, functional_name
+      character(len=OPTION_PATH_LEN) :: simulation_base_name
       type(stat_type), dimension(:), allocatable :: functional_stats
+      real :: J
+      real, dimension(:), pointer :: fn_value
 
-      character(len=ADJ_NAME_LEN) :: variable_name, field_name, material_phase_name
+      character(len=ADJ_NAME_LEN) :: variable_name, field_name, material_phase_name, functional_name
       type(scalar_field) :: sfield_soln, sfield_rhs
       type(vector_field) :: vfield_soln, vfield_rhs
       type(csr_matrix) :: csr_mat
@@ -112,17 +133,29 @@
       allocate(functional_stats(no_functionals))
 
       do functional=0,no_functionals-1
+        call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+        if (have_option("/adjoint/functional::" // trim(functional_name) // "/disable_adjoint_run")) then
+         cycle
+        end if 
         default_stat = functional_stats(functional + 1)
         call initialise_walltime
-        call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
         call initialise_diagnostics(trim(simulation_base_name) // '_' // trim(functional_name), state)
         functional_stats(functional + 1) = default_stat
-
+        
+        ! Register the callback to compute J
+        ierr = adj_register_functional_callback(adjointer, trim(functional_name), c_funloc(libadjoint_evaluate_functional))
+        call adj_chkierr(ierr)
+       
         ! Register the callback to compute delJ/delu
         ierr = adj_register_functional_derivative_callback(adjointer, trim(functional_name), c_funloc(libadjoint_functional_derivative))
         call adj_chkierr(ierr)
       end do
 
+      ! Insert the fields for storing the derivatives of the functionals with respect to the controls into the state
+       if (have_option("/adjoint/controls")) then
+        call allocate_and_insert_control_derivative_fields(state)
+      end if
+
       ierr = adj_timestep_count(adjointer, no_timesteps)
       call adj_chkierr(ierr)
 
@@ -145,18 +178,6 @@
           call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
           call set_option("/simulation_name", trim(simulation_base_name) // "_" // trim(functional_name))
           default_stat = functional_stats(functional + 1)
-          ! We don't want to compute a functional for the dummy timestep added at the end to act
-          ! as a container for the last equation
-          if (start_time == end_time) then
-            assert(timestep == no_timesteps-1)
-            if (have_option("/adjoint/functional[" // int2str(functional) // "]/functional_value")) then
-              call set_diagnostic(name=trim(functional_name), statistic="value", value=(/0.0/))
-              call set_diagnostic(name=trim(functional_name) // "_component", statistic="value", value=(/0.0/))
-            end if
-            functional_computed = .true.
-          else
-            functional_computed = .false.
-          end if
 
           do equation=end_timestep,start_timestep,-1
             ierr = adj_get_adjoint_equation(adjointer, equation, trim(functional_name), lhs, rhs, adj_var)
@@ -201,7 +222,7 @@
                       end if
 
                       call petsc_solve(sfield_soln, csr_mat, sfield_rhs, option_path=path)
-                      !call compute_inactive_rows(sfield_soln, csr_mat, sfield_rhs)
+                      call compute_inactive_rows(sfield_soln, csr_mat, sfield_rhs)
                     endif
                   case(ADJ_BLOCK_CSR_MATRIX)
                     FLAbort("Cannot map between scalar fields with a block_csr_matrix .. ")
@@ -240,7 +261,7 @@
                       end if
 
                       call petsc_solve(vfield_soln, csr_mat, vfield_rhs, option_path=path)
-                      call compute_inactive_rows(vfield_soln, csr_mat, vfield_rhs)
+                      !call compute_inactive_rows(vfield_soln, csr_mat, vfield_rhs)
                     endif
                   case(ADJ_BLOCK_CSR_MATRIX)
                     call matrix_from_adj_matrix(lhs, block_csr_mat)
@@ -282,7 +303,16 @@
               call adjoint_cleanup
               return
             end if
-          end do
+          end do ! End of equation loop
+          
+          ! Callback for the computation of the functional total derivatives
+          if (present(adjoint_timestep_callback)) then
+            call adjoint_timestep_callback(state, timestep, trim(functional_name), data=adjoint_timestep_callback_data)
+          end if
+          if (have_option("/adjoint/controls")) then
+            ! Write the functional's total derivatives to file
+            call adjoint_write_control_derivatives(state, timestep, trim(functional_name))
+          end if
 
           call calculate_diagnostic_variables(state, exclude_nonrecalculated = .true.)
           call calculate_diagnostic_variables_new(state, exclude_nonrecalculated = .true.)
@@ -296,13 +326,13 @@
           endif
 
           functional_stats(functional + 1) = default_stat
-        end do
+        end do ! End of functional loop
 
         ! Now forget
         ierr = adj_forget_adjoint_equation(adjointer, start_timestep)
         call adj_chkierr(ierr)
         current_time = start_time
-      end do
+      end do ! End of timestep loop
 
       call get_option("/timestepping/finish_time", finish_time)
       assert(current_time == finish_time)
@@ -314,6 +344,10 @@
       subroutine adjoint_cleanup
         ! Clean up stat files
         do functional=0,no_functionals-1
+          call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+          if (have_option("/adjoint/functional::" // trim(functional_name) // "/disable_adjoint_run")) then
+           cycle
+          end if 
           default_stat = functional_stats(functional + 1)
           call close_diagnostic_files
           call uninitialise_diagnostics
@@ -325,23 +359,7 @@
         call adj_chkierr(ierr)
       end subroutine adjoint_cleanup
     end subroutine compute_adjoint
+
 #endif
 
-    subroutine adjoint_main_loop_register_diagnostic
-      integer :: functional, no_functionals
-      character(len=OPTION_PATH_LEN) :: functional_name
-
-      no_functionals = option_count("/adjoint/functional")
-      
-      do functional=0,no_functionals-1
-        ! Register a diagnostic for each functional 
-        if (have_option("/adjoint/functional[" // int2str(functional) // "]/functional_value")) then
-          call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
-          call register_diagnostic(dim=1, name=trim(functional_name) // "_component", statistic="value")
-          call register_diagnostic(dim=1, name=trim(functional_name), statistic="value")
-        end if
-      end do
-
-       end subroutine adjoint_main_loop_register_diagnostic
-
 end module adjoint_main_loop
 
=== modified file 'adjoint/Adjoint_Python_Fortran.F90'
--- adjoint/Adjoint_Python_Fortran.F90	2011-05-01 19:08:55 +0000
+++ adjoint/Adjoint_Python_Fortran.F90	2011-07-12 15:56:32 +0000
@@ -40,9 +40,11 @@
   public :: adj_variables_from_python
 
   interface
-    subroutine adj_variables_from_python_c(fn, fnlen, start_time, end_time, timestep, result, result_len, stat) &
+    subroutine adj_variables_from_python_c(adjointer, fn, fnlen, start_time, end_time, timestep, result, result_len, stat) &
                                          & bind(c, name='adj_variables_from_python')
+        use libadjoint
         use iso_c_binding
+        type(adj_adjointer), intent(in) :: adjointer
         integer(kind=c_int), intent(in), value :: fnlen
         character(kind=c_char), dimension(fnlen), intent(in) :: fn
         real(kind=c_double), intent(in), value :: start_time, end_time
@@ -60,7 +62,8 @@
 
   contains
 
-  subroutine adj_variables_from_python(fn, start_time, end_time, timestep, result, extras, stat)
+  subroutine adj_variables_from_python(adjointer, fn, start_time, end_time, timestep, result, extras, stat)
+    type(adj_adjointer), intent(in) :: adjointer
     character(len=*), intent(in) :: fn
     real, intent(in) :: start_time, end_time
     integer, intent(in) :: timestep
@@ -88,7 +91,7 @@
 
     timestep_long = timestep
 
-    call adj_variables_from_python_c(fn_c, len_trim(fn), start_time, end_time, timestep_long, result_c_ptr, result_len_c, stat_c)
+    call adj_variables_from_python_c(adjointer, fn_c, len_trim(fn), start_time, end_time, timestep_long, result_c_ptr, result_len_c, stat_c)
 
     if (stat_c /= 0) then
       if (present(stat)) then
 
=== modified file 'adjoint/Burgers_Adjoint_Callbacks.F90'
--- adjoint/Burgers_Adjoint_Callbacks.F90	2011-04-01 17:40:34 +0000
+++ adjoint/Burgers_Adjoint_Callbacks.F90	2011-07-12 15:56:32 +0000
@@ -40,6 +40,9 @@
     use adjoint_global_variables, only: adj_path_lookup
     use mangle_options_tree, only: adjoint_field_path
     use mangle_dirichlet_rows_module
+    use populate_state_module
+    use burgers_assembly
+    use global_parameters, only: OPTION_PATH_LEN, running_adjoint
     implicit none
 
     private
@@ -61,6 +64,15 @@
       call adj_chkierr(ierr)
       ierr = adj_register_operator_callback(adjointer, ADJ_BLOCK_ACTION_CB, "TimesteppingOperator", c_funloc(timestepping_operator_action_callback))
       call adj_chkierr(ierr)
+
+      ierr = adj_register_operator_callback(adjointer, ADJ_NBLOCK_DERIVATIVE_ACTION_CB, "AdvectionOperator", c_funloc(advection_derivative_action_proc))
+      call adj_chkierr(ierr)
+
+      ierr = adj_register_operator_callback(adjointer, ADJ_NBLOCK_ACTION_CB, "AdvectionOperator", c_funloc(advection_action_proc))
+      call adj_chkierr(ierr)
+
+      ierr = adj_register_forward_source_callback(adjointer, c_funloc(burgers_equation_forward_source))
+      call adj_chkierr(ierr)
     end subroutine register_burgers_operator_callbacks
 
     subroutine velocity_identity_assembly_callback(nvar, variables, dependencies, hermitian, coefficient, context, output, rhs) bind(c)
@@ -114,9 +126,11 @@
       type(scalar_field) :: empty_velocity_field, previous_u, iter_u
       type(scalar_field), dimension(2) :: deps
       type(csr_matrix), pointer :: mass_mat, diffusion_mat
-      type(csr_matrix) :: burgers_mat, burgers_mat_T
+      type(csr_matrix) :: burgers_mat, burgers_mat_T, advection_mat
+      type(vector_field), pointer :: positions
       real :: dt, theta
       integer :: i
+      character(len=FIELD_NAME_LEN) :: path
 
       if (coefficient /= 1.0) then
         FLAbort("The coefficient in burgers_operator_assembly_callback has to be 1.0")
@@ -139,13 +153,20 @@
         FLAbort("Unknown number of dependencies in burgers_operator_assembly_callback")
       end if
 
-      call get_option(trim(adjoint_field_path("/material_phase::Fluid/scalar_field::Velocity/prognostic/temporal_discretisation/theta")), theta)
+      path = "/material_phase::Fluid/scalar_field::Velocity"
+      if (running_adjoint) then
+        path = adjoint_field_path("/material_phase::Fluid/scalar_field::Velocity")
+      end if
+
+      call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta, default=0.5)
       call get_option("/timestepping/timestep", dt)
+      dt = abs(dt)
 
       call c_f_pointer(context, matrices)
       u_mesh => extract_mesh(matrices, "VelocityMesh")
       diffusion_mat => extract_csr_matrix(matrices, "DiffusionMatrix")
       mass_mat => extract_csr_matrix(matrices, "MassMatrix")
+      positions => extract_vector_field(matrices, "Coordinate")
 
       call allocate(empty_velocity_field, u_mesh, trim("AdjointPressureRhs"))
       call zero(empty_velocity_field)
@@ -156,10 +177,18 @@
       call allocate(burgers_mat, mass_mat%sparsity, name="BurgersAdjointOutput")
       call zero(burgers_mat)
 
-      if (.not. have_option(trim(adjoint_field_path('/material_phase::Fluid/scalar_field::Velocity/prognostic/temporal_discretisation/remove_time_term')))) then
+      if (.not. have_option(trim(path) // '/prognostic/temporal_discretisation/remove_time_term')) then
         call addto(burgers_mat, mass_mat, 1.0/abs(dt))
       end if
-      !call addto(burgers_mat, advection_mat, theta)
+
+      if (.not. have_option(trim(path) // '/prognostic/remove_advection_term')) then
+        call allocate(advection_mat, mass_mat%sparsity, name="BurgersCallbackAdvectionMatrix")
+        call zero(advection_mat)
+        call assemble_advection_matrix(advection_mat, positions, previous_u, iter_u)
+        call addto(burgers_mat, advection_mat, theta)
+        call deallocate(advection_mat)
+      end if
+
       call addto(burgers_mat, diffusion_mat, theta)
 
       ! Ensure that the input velocity has indeed dirichlet bc attached
@@ -223,12 +252,15 @@
 
       type(state_type), pointer :: matrices
       type(csr_matrix), pointer :: diffusion_mat, mass_mat
-      type(csr_matrix) :: timestepping_mat
       type(mesh_type), pointer :: u_mesh
       ! Input/output variables
       type(scalar_field) :: u_input, previous_u, iter_u
-      type(scalar_field) :: u_output
+      type(scalar_field) :: u_output, tmp_u
       real :: theta, dt
+      character(len=FIELD_NAME_LEN) :: path
+      type(csr_matrix) :: advection_mat
+      type(vector_field), pointer :: positions
+      type(csr_matrix) :: timestepping_mat
 
       call c_f_pointer(context, matrices)
       
@@ -240,6 +272,7 @@
       u_mesh => extract_mesh(matrices, "VelocityMesh")
       diffusion_mat => extract_csr_matrix(matrices, "DiffusionMatrix")
       mass_mat => extract_csr_matrix(matrices, "MassMatrix")
+      positions => extract_vector_field(matrices, "Coordinate")
 
       if (nvar==2) then
         if (variables(1)%timestep==variables(2)%timestep-1) then
@@ -259,29 +292,45 @@
       end if
 
       ! Ensure that the input velocity has indeed dirichlet bc attached
-      ! Ensure that the input velocity has indeed dirichlet bc attached
       if (get_boundary_condition_count(previous_u)==0) then
         FLAbort("Velocity in callback must have dirichlet conditions attached")
       end if
 
       call field_from_adj_vector(input, u_input)
-      
+
       call allocate(u_output, u_mesh, "TimeSteppingAdjointVelocityOutput")
       call zero(u_output)
+      call allocate(tmp_u, u_mesh, "TimeSteppingAdjointVelocityOutput")
+      call zero(tmp_u)
 
       call allocate(timestepping_mat, mass_mat%sparsity, name="TimesteppingMatrix")
       call zero(timestepping_mat)
 
-      if (.not. have_option(trim(adjoint_field_path('/material_phase::Fluid/scalar_field::Velocity/prognostic/temporal_discretisation/remove_time_term')))) then
-        call addto(timestepping_mat, mass_mat, scale=1.0/dt)
-      end if
-
-!        call mult(tmp, advection_mat, u_input)
-!        call scale(tmp, theta - 1.0)
-!        call addto(u_output, tmp)
-
-      call addto(timestepping_mat, diffusion_mat, scale=theta - 1.0)
-
+      if (running_adjoint) then
+        path = adjoint_field_path("/material_phase::Fluid/scalar_field::Velocity")
+      else
+        path = "/material_phase::Fluid/scalar_field::Velocity"
+      end if
+
+      if (.not. have_option(trim(path) // '/prognostic/temporal_discretisation/remove_time_term')) then
+        call addto(timestepping_mat, mass_mat, scale=1.0/abs(dt))
+      end if
+
+      if (.not. have_option(trim(path) // '/prognostic/remove_advection_term')) then
+        call allocate(advection_mat, mass_mat%sparsity, name="BurgersCallbackAdvectionMatrix")
+        call zero(advection_mat)
+        call assemble_advection_matrix(advection_mat, positions, previous_u, iter_u)
+        call addto(timestepping_mat, advection_mat, scale=theta-1.0)
+        call deallocate(advection_mat)
+      end if
+
+      call addto(timestepping_mat, diffusion_mat, scale=theta-1.0)
+
+      call scale(timestepping_mat, -1.0)
+      if (get_boundary_condition_count(previous_u)==0) then
+        FLAbort("Velocity in callback must have dirichlet conditions attached")
+      end if
+      ! Mangle the dirichlet rows
       call mangle_dirichlet_rows(timestepping_mat, previous_u, keep_diag=.false.)
 
       if (hermitian==ADJ_FALSE) then
@@ -289,12 +338,283 @@
       else
         call mult_T(u_output, timestepping_mat, u_input)
       end if
+
+      call deallocate(timestepping_mat)
+
       call scale(u_output, coefficient)
       output = field_to_adj_vector(u_output)
       call deallocate(u_output)
-      call deallocate(timestepping_mat)
+      call deallocate(tmp_u)
     end subroutine timestepping_operator_action_callback
 
-
+    subroutine advection_derivative_action_proc(nvar, variables, dependencies, derivative, contraction, hermitian, &
+                                                  & input, coefficient, context, output) bind(c)
+      use iso_c_binding
+      use libadjoint_data_structures
+      use simple_advection_d
+      use simple_advection_b
+      integer(kind=c_int), intent(in), value :: nvar
+      type(adj_variable), dimension(nvar), intent(in) :: variables
+      type(adj_vector), dimension(nvar), intent(in) :: dependencies
+      type(adj_variable), intent(in), value :: derivative
+      type(adj_vector), intent(in), value :: contraction
+      integer(kind=c_int), intent(in), value :: hermitian
+      type(adj_vector), intent(in), value :: input
+      adj_scalar_f, intent(in), value :: coefficient
+      type(c_ptr), intent(in), value :: context
+      type(adj_vector), intent(out) :: output
+
+      type(state_type), pointer :: matrices
+      type(vector_field), pointer :: positions
+      type(scalar_field) :: u_left, u_right
+      type(scalar_field) :: contraction_field, udot, Acbar
+      type(scalar_field) :: output_field, tmp_field, u
+      character(len=OPTION_PATH_LEN) :: path
+
+      real :: itheta
+
+      assert(nvar == 1 .or. nvar == 2)
+      call c_f_pointer(context, matrices)
+      positions => extract_vector_field(matrices, "Coordinate")
+
+      call field_from_adj_vector(contraction, contraction_field)
+      call allocate(output_field, contraction_field%mesh, "NonlinearDerivativeOutput")
+      call zero(output_field)
+
+      if (have_option("/material_phase::Fluid/scalar_field::Velocity")) then
+        path = "/material_phase::Fluid/scalar_field::Velocity"
+      else
+        path = "/material_phase::Fluid/scalar_field::Velocity"
+        path = adjoint_field_path(path)
+      end if
+
+      if (hermitian == ADJ_FALSE) then
+        call field_from_adj_vector(input, udot)
+        call get_option(trim(path) // "/prognostic/temporal_discretisation/relaxation", itheta)
+        if (have_option(trim(path) // "/prognostic/remove_advection_term")) then
+          output = field_to_adj_vector(output_field)
+          call deallocate(output_field)
+          return
+        end if
+      else
+        call get_option(trim(path) // "/prognostic/temporal_discretisation/relaxation", itheta)
+        call field_from_adj_vector(input, Acbar)
+        if (have_option(trim(path) // "/prognostic/remove_advection_term")) then
+          output = field_to_adj_vector(output_field)
+          call deallocate(output_field)
+          return
+        end if
+      end if
+
+      call allocate(tmp_field, contraction_field%mesh, "TmpOutput")
+      call zero(tmp_field)
+
+      if (nvar == 1) then
+        call field_from_adj_vector(dependencies(1), u_left)
+        if (hermitian == ADJ_FALSE) then
+          call advection_action_d(positions%val(1,:), u_left%val, udot%val, contraction_field%val, tmp_field%val, output_field%val)
+        else
+          call advection_action_b(positions%val(1,:), u_left%val, output_field%val, contraction_field%val, tmp_field%val, Acbar%val)
+        end if
+      else if (nvar == 2) then
+        call field_from_adj_vector(dependencies(1), u_left)
+        call field_from_adj_vector(dependencies(2), u_right)
+        call allocate(u, u_left%mesh, "AdvectingVelocity")
+        call set(u, u_left)
+        call scale(u, (1.0-itheta))
+        call addto(u, u_right, scale=itheta)
+
+        if (hermitian == ADJ_FALSE) then
+          call advection_action_d(positions%val(1,:), u%val, udot%val, contraction_field%val, tmp_field%val, output_field%val)
+        else
+          call advection_action_b(positions%val(1,:), u%val, output_field%val, contraction_field%val, tmp_field%val, Acbar%val)
+        end if
+
+        call deallocate(u)
+
+        if (derivative == variables(1)) then
+          call scale(output_field, (1.0-itheta))
+        else if (derivative == variables(2)) then
+          call scale(output_field, itheta)
+        end if
+      end if
+
+      call deallocate(tmp_field)
+      call scale(output_field, coefficient)
+      output = field_to_adj_vector(output_field)
+      call deallocate(output_field)
+    end subroutine advection_derivative_action_proc
+
+    subroutine burgers_equation_forward_source(adjointer, var, ndepends, dependencies, values, context, output, has_output) bind(c)
+      type(adj_adjointer), intent(in) :: adjointer
+      type(adj_variable), intent(in), value :: var
+      integer(kind=c_int), intent(in), value :: ndepends
+      type(adj_variable), dimension(ndepends), intent(in) :: dependencies
+      type(adj_vector), dimension(ndepends), intent(in) :: values
+      type(c_ptr), intent(in), value :: context
+      type(adj_vector), intent(out) :: output
+      integer(kind=c_int), intent(out) :: has_output
+
+      character(len=ADJ_DICT_LEN) :: path
+      character(len=ADJ_NAME_LEN) :: name
+      integer :: timestep
+      real :: time, dt, end_time
+      real :: theta
+      type(state_type), pointer :: matrices
+      logical :: has_source
+      type(mesh_type), pointer :: u_mesh
+      type(state_type), dimension(1) :: dummy_state
+      type(scalar_field) :: u_output, tmp_u_output, dummy_u
+      type(csr_matrix), pointer :: mass_matrix
+      type(vector_field), pointer :: positions
+      type(csr_matrix) :: mangled_mass_matrix
+
+      integer :: ierr
+
+      ierr = adj_variable_get_name(var, name)
+      call adj_chkierr(ierr)
+
+      ierr = adj_variable_get_timestep(var, timestep)
+      call adj_chkierr(ierr)
+
+      call c_f_pointer(context, matrices)
+      assert(associated(matrices))
+      u_mesh => extract_mesh(matrices, "VelocityMesh")
+      positions => extract_vector_field(matrices, "Coordinate")
+      mass_matrix => extract_csr_matrix(matrices, "MassMatrix")
+
+      path = "/material_phase::Fluid/scalar_field::Velocity"
+
+      has_source = have_option(trim(path) // "/prognostic/scalar_field::Source")
+
+      if (timestep == 0) then ! initial condition
+        call allocate(u_output, u_mesh, "VelocityInitialCondition")
+        call zero(u_output)
+        u_output%option_path = trim(path)
+        call insert(dummy_state(1), positions, "Coordinate")
+        call insert(dummy_state(1), u_output, "Velocity")
+        call initialise_prognostic_fields(dummy_state)
+        
+        ! Set initial condition
+        call populate_boundary_conditions(dummy_state)
+        call set_boundary_conditions_values(dummy_state, shift_time=.false.)
+        call set_dirichlet_consistent(dummy_state)
+        call deallocate(dummy_state(1))
+
+        output = field_to_adj_vector(u_output)
+        has_output = ADJ_TRUE
+        call deallocate(u_output)
+      else if (timestep > 0) then
+        if (.not. has_source) then
+          has_output = ADJ_FALSE
+        else
+          call allocate(tmp_u_output, u_mesh, "VelocitySource")
+          call allocate(u_output, u_mesh, "VelocitySource")
+          call allocate(dummy_u, u_mesh, "DummyVelocity")
+          call zero(tmp_u_output)
+          call zero(u_output)
+          tmp_u_output%option_path = trim(path) // "/prognostic/scalar_field::Source"
+          dummy_u%option_path = trim(path)
+
+          call insert(dummy_state(1), positions, "Coordinate")
+          call insert(dummy_state(1), tmp_u_output, "VelocitySource")
+          call insert(dummy_state(1), dummy_u, "DummyVelocity")
+
+          call populate_boundary_conditions(dummy_state)
+
+          ierr = adj_timestep_get_times(adjointer, timestep-1, time, end_time)
+          call adj_chkierr(ierr)
+          dt = end_time - time
+
+          call get_option(trim(path) // "/prognostic/temporal_discretisation/theta", theta, default=0.5)
+          call set_prescribed_field_values(dummy_state, time=time + theta*abs(dt))
+          call deallocate(dummy_state(1))
+
+          call allocate(mangled_mass_matrix, mass_matrix%sparsity, name="MangledMassMatrix")
+          call set(mangled_mass_matrix, mass_matrix)
+          call mangle_dirichlet_rows(mangled_mass_matrix, dummy_u, keep_diag=.false.)
+          call mult(u_output, mangled_mass_matrix, tmp_u_output)
+          call deallocate(mangled_mass_matrix)
+
+          call deallocate(tmp_u_output)
+          call deallocate(dummy_u)
+
+          output = field_to_adj_vector(u_output)
+          has_output = ADJ_TRUE
+          call deallocate(u_output)
+        end if
+      end if
+
+    end subroutine burgers_equation_forward_source
+
+    subroutine advection_action_proc(nvar, variables, dependencies, input, context, output) bind(c)
+      use iso_c_binding
+      use libadjoint_data_structures
+      use simple_advection
+      integer(kind=c_int), intent(in), value :: nvar
+      type(adj_variable), dimension(nvar), intent(in) :: variables
+      type(adj_vector), dimension(nvar), intent(in) :: dependencies
+      type(adj_vector), intent(in), value :: input
+      type(c_ptr), intent(in), value :: context
+      type(adj_vector), intent(out) :: output
+
+      type(scalar_field) :: u_output, u_check, nonlinear_u
+      type(scalar_field) :: previous_u, iter_u, u_input
+      type(csr_matrix) :: advection_mat
+      type(csr_matrix), pointer :: diffusion_mat
+      type(vector_field), pointer :: positions
+      type(mesh_type), pointer :: u_mesh
+      type(state_type), pointer :: matrices
+      character(len=OPTION_PATH_LEN) :: path
+      real :: itheta
+
+      path = "/material_phase::Fluid/scalar_field::Velocity"
+      if (running_adjoint) then
+        path = adjoint_field_path(path)
+      end if
+
+      if (nvar==2) then
+        if (variables(1)%timestep==variables(2)%timestep-1) then
+          call field_from_adj_vector(dependencies(1), previous_u)
+          call field_from_adj_vector(dependencies(2), iter_u)
+        else if(variables(2)%timestep==variables(1)%timestep-1) then
+          call field_from_adj_vector(dependencies(2), previous_u)
+          call field_from_adj_vector(dependencies(1), iter_u)
+        else
+          FLAbort("Dependencies have no contiguous timesteps.")
+        end if
+      else if (nvar==1) then
+        call field_from_adj_vector(dependencies(1), previous_u)
+        call field_from_adj_vector(dependencies(1), iter_u)
+      else
+        FLAbort("Unknown number of dependencies in burgers_operator_assembly_callback")
+      end if
+
+      call field_from_adj_vector(input, u_input)
+
+      call c_f_pointer(context, matrices)
+      u_mesh => extract_mesh(matrices, "VelocityMesh")
+      diffusion_mat => extract_csr_matrix(matrices, "DiffusionMatrix")
+      positions => extract_vector_field(matrices, "Coordinate")
+
+      call allocate(advection_mat, diffusion_mat%sparsity, name="AdvectionActionMatrix")
+      call zero(advection_mat)
+      if (.not. have_option(trim(path) // "/prognostic/remove_advection_term")) then
+        call assemble_advection_matrix(advection_mat, positions, previous_u, iter_u)
+      end if
+
+      if (get_boundary_condition_count(previous_u) /= 1) then
+        FLAbort("Need to supply boundary conditions on previous_u")
+      end if
+
+      call mangle_dirichlet_rows(advection_mat, previous_u, keep_diag=.true.)
+
+      call allocate(u_output, u_mesh, "AdvectionActionOutput")
+      call mult(u_output, advection_mat, u_input)
+      call deallocate(advection_mat)
+
+      output = field_to_adj_vector(u_output)
+      call deallocate(u_output)
+    end subroutine advection_action_proc
 #endif
 end module burgers_adjoint_callbacks
 
=== added file 'adjoint/Burgers_Control_Callbacks.F90'
--- adjoint/Burgers_Control_Callbacks.F90	1970-01-01 00:00:00 +0000
+++ adjoint/Burgers_Control_Callbacks.F90	2011-07-12 15:56:32 +0000
@@ -0,0 +1,183 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+
+#include "fdebug.h"
+
+module burgers_adjoint_controls
+  use fields
+  use global_parameters, only : PYTHON_FUNC_LEN, OPTION_PATH_LEN
+  use spud
+  use state_module
+  use python_state
+  use sparse_matrices_fields
+  use manifold_projections
+  use mangle_dirichlet_rows_module
+
+    implicit none
+
+    private
+
+    public :: burgers_adjoint_timestep_callback
+
+    contains 
+    
+    ! Data is a void pointer used to pass variables into the callback
+    subroutine burgers_adjoint_timestep_callback(states, timestep, functional_name, data)
+      use iso_c_binding, only: c_ptr
+      use global_parameters, only: OPTION_PATH_LEN
+      use state_module
+      type(state_type), dimension(:), intent(inout) :: states
+      integer, intent(in) :: timestep
+      character(len=*), intent(in) :: functional_name
+      type(c_ptr), intent(in) :: data
+      
+      type(state_type), pointer :: matrices
+      character(len=OPTION_PATH_LEN) :: field_name, control_type, control_deriv_name, field_deriv_name, name, material_phase_name
+      integer :: nb_controls
+      integer :: i, state_id, s_idx
+      type(scalar_field), pointer :: sfield, adj_sfield, adj_sfield2
+      type(vector_field), pointer :: vfield, adj_vfield, positions
+      type(tensor_field), pointer :: tfield, adj_tfield
+      type(block_csr_matrix), pointer :: big_mat, div_mat, u_mass_mat
+      type(csr_matrix), pointer :: h_mass_mat
+      type(csr_matrix) :: h_mangled_mass_mat
+      type(state_type), dimension(1) :: dummy_state
+      
+      type(scalar_field) :: dummy_u
+      type(vector_field) :: local_src_tmp, local_src_tmp2, src_tmp
+      real :: theta, d0, dt 
+      integer :: stat
+     
+      ! Cast the data to the matrices state 
+      call c_f_pointer(data, matrices)
+      
+
+      nb_controls = option_count("/adjoint/controls/control")
+      do i = 0, nb_controls-1
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/name", control_deriv_name)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type/name", control_type)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type::" // trim(control_type) // "/field_name", name)
+        s_idx = scan(trim(name), "::")
+        if (s_idx == 0) then 
+          FLAbort("The control " // trim(control_deriv_name) // " uses an invalid field_name. It should be of the form Materialphase::Fieldname")
+        end if
+        material_phase_name = name(1:s_idx - 1)
+        field_name = name(s_idx + 2:len_trim(name))
+        ! Find state associated with the material phase
+        do state_id = 1, size(states) 
+          if (trim(states(state_id)%name) == trim(material_phase_name)) then
+            exit
+          end if
+        end do
+        ! Make sure we found the state
+        if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
+          FLAbort("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_deriv_name) // ".")
+        end if
+
+        control_deriv_name = trim(control_deriv_name) // "_TotalDerivative"
+        select case(trim(control_type))
+
+          !!!!!!!!!!!!! Initial condition !!!!!!!!!!!!
+          case ("initial_condition")
+            if (timestep /= 0) then
+             cycle
+            end if 
+            field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+            if (has_scalar_field(states(state_id), field_deriv_name)) then
+              if (trim(field_name) == "Velocity") then
+                adj_sfield => extract_scalar_field(states(state_id), "Adjoint" // trim(field_name))
+                sfield => extract_scalar_field(states(state_id), field_deriv_name) ! Output field
+                positions => extract_vector_field(matrices, "Coordinate")
+                call set(sfield, adj_sfield)
+
+                ! Populate boundary conditions
+                call allocate(dummy_u, adj_sfield%mesh, "DummyVelocity")
+                dummy_u%option_path =  "/material_phase::Fluid/scalar_field::AdjointVelocity"
+                call insert(dummy_state(1), positions, "Coordinate")
+                call insert(dummy_state(1), dummy_u, "AdjointVelocity")
+                call populate_boundary_conditions(dummy_state)
+
+                ! Now we need to zero the derivative field everywhere we have Dirichlet conditions
+                call mangle_dirichlet_rows(sfield, dummy_u)
+
+                call deallocate(dummy_u)
+                call deallocate(dummy_state(1))
+              else
+                FLAbort("Sorry, I do not know how to compute the intial condition control for " // trim(field_name) // ".")
+              end if
+            else
+              FLAbort("The control derivative field " // trim(field_deriv_name) // " specified for " // trim(control_deriv_name) // " is not a field in the state.")
+            end if
+
+          !!!!!!!!!!!!! Boundary condition !!!!!!!!!!!!
+          case("boundary_condition")
+            FLAbort("Boundary condition control not implemented yet.")
+
+          !!!!!!!!!!!!!       Source       !!!!!!!!!!!!
+          case("source_term")
+            field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+            sfield => extract_scalar_field(states(state_id), field_deriv_name, stat=stat) ! Output field
+            if (stat /= 0) then
+              FLAbort("The control derivative field " // trim(field_deriv_name) // " specified for " // trim(control_deriv_name) // " is not a field in the state.")
+            end if
+
+            if (timestep == 0) then
+             call zero(sfield)
+            else
+              if (trim(field_name) == "VelocitySource") then
+                adj_sfield => extract_scalar_field(states(state_id), "AdjointVelocity")
+                h_mass_mat => extract_csr_matrix(matrices, "MassMatrix")
+                positions => extract_vector_field(matrices, "Coordinate")
+                
+                ! Populate boundary conditions
+                call allocate(dummy_u, adj_sfield%mesh, "DummyVelocity")
+                dummy_u%option_path =  "/material_phase::Fluid/scalar_field::AdjointVelocity"
+                call insert(dummy_state(1), positions, "Coordinate")
+                call insert(dummy_state(1), dummy_u, "AdjointVelocity")
+                call populate_boundary_conditions(dummy_state)
+
+                ! Mangle the mass matrix
+                call allocate(h_mangled_mass_mat, h_mass_mat%sparsity, name="MangledMassMatrix")
+                call set(h_mangled_mass_mat, h_mass_mat)
+                call mangle_dirichlet_rows(h_mangled_mass_mat, dummy_u, keep_diag=.false.)
+                call mult_T(sfield, h_mangled_mass_mat, adj_sfield)
+
+                ! Deallocate temporary variables
+                call deallocate(dummy_state(1))
+                call deallocate(dummy_u)
+                call deallocate(h_mangled_mass_mat)
+
+              else
+                FLAbort("Sorry, I do not know how to compute the source condition control for " // trim(field_name) // ".")
+              end if
+            end if
+        end select
+      end do
+
+    end subroutine burgers_adjoint_timestep_callback 
+
+end module burgers_adjoint_controls
 
=== modified file 'adjoint/Forward_Main_Loop.F90'
--- adjoint/Forward_Main_Loop.F90	2011-05-13 17:52:57 +0000
+++ adjoint/Forward_Main_Loop.F90	2011-07-12 15:56:32 +0000
@@ -40,6 +40,7 @@
     use adjoint_global_variables
     use adjoint_functional_evaluation
     use populate_state_module
+    use boundary_conditions_from_options
     use signal_vars
     use mangle_options_tree
     use mangle_dirichlet_rows_module
@@ -53,12 +54,13 @@
     implicit none
 
     private
+    public :: forward_main_loop_register_diagnostic
 #ifdef HAVE_ADJOINT
-    public :: compute_forward
+    public :: compute_forward, calculate_functional_values, register_functional_callbacks
 #endif
 
-#ifdef HAVE_ADJOINT
     contains
+#ifdef HAVE_ADJOINT
 
     subroutine compute_forward(state)
       type(state_type), dimension(:), intent(inout) :: state
@@ -153,6 +155,11 @@
 
               if (has_path) then
                 sfield_soln%option_path = trim(path)
+                ! We need to populate the BC values:
+                call insert(state(1), sfield_soln, trim(sfield_soln%name))
+                call populate_boundary_conditions(state)
+                call set_boundary_conditions_values(state, shift_time=.false.)
+                sfield_soln = extract_scalar_field(state(1), trim(sfield_soln%name))
               end if
 
               select case(lhs%klass)
@@ -168,6 +175,10 @@
                       call adj_chkierr(ierr)
                     end if
 
+                    sfield_rhs%bc => sfield_soln%bc
+                    call set_dirichlet_consistent(sfield_rhs)
+                    sfield_rhs%bc => null()
+
                     call petsc_solve(sfield_soln, csr_mat, sfield_rhs, option_path=path)
                     call compute_inactive_rows(sfield_soln, csr_mat, sfield_rhs)
                   endif
@@ -178,6 +189,7 @@
               end select
 
               call insert(state(1), sfield_soln, trim(sfield_soln%name))
+
               soln = field_to_adj_vector(sfield_soln)
               ierr = adj_storage_memory_incref(soln, storage)
               call adj_chkierr(ierr)
@@ -197,6 +209,11 @@
 
               if (has_path) then
                 vfield_soln%option_path = trim(path)
+                ! We need to populate the BC values:
+                call insert(state(1), vfield_soln, trim(vfield_soln%name))
+                call populate_boundary_conditions(state)
+                call set_boundary_conditions_values(state, shift_time=.false.)
+                vfield_soln = extract_vector_field(state(1), trim(vfield_soln%name))
               end if
 
               select case(lhs%klass)
@@ -212,6 +229,10 @@
                       call adj_chkierr(ierr)
                     end if
 
+                    vfield_rhs%bc => vfield_soln%bc
+                    call set_dirichlet_consistent(vfield_rhs)
+                    vfield_rhs%bc => null()
+
                     call petsc_solve(vfield_soln, csr_mat, vfield_rhs, option_path=path)
                     !call compute_inactive_rows(vfield_soln, csr_mat, vfield_rhs)
                   endif
@@ -233,6 +254,7 @@
               end select
 
               call insert(state(1), vfield_soln, trim(vfield_soln%name))
+
               soln = field_to_adj_vector(vfield_soln)
               ierr = adj_storage_memory_incref(soln, storage)
               call adj_chkierr(ierr)
@@ -264,8 +286,17 @@
         call set_prescribed_field_values(state, exclude_interpolated=.true., exclude_nonreprescribed=.true., time=current_time)
         call calculate_diagnostic_variables(state, exclude_nonrecalculated = .true.)
         call calculate_diagnostic_variables_new(state, exclude_nonrecalculated = .true.)
+        ! The first timestep is the initialisation of the model.
+        ! We skip the evaluation of the functional at timestep zero to get the correct value.
+        if (timestep > 0) then
+          call calculate_functional_values(timestep-1)
+          ! The last timestep is a the dummy timestep added at the end to act
+          ! as a container for the last equation
+          if (start_time == end_time) then
+            assert(timestep == no_timesteps-1)
+          end if
+        end if
         call write_diagnostics(state, current_time, dt, equation+1)
-
         if (do_write_state(current_time, timestep)) then
           call write_state(dump_no, state)
         endif
@@ -286,5 +317,69 @@
       call get_option("/timestepping/finish_time", finish_time)
       assert(current_time == finish_time)
     end subroutine compute_forward
-#endif
+
+    subroutine register_functional_callbacks()                                                                                                                                                                 
+      integer :: no_functionals, functional
+      character(len=ADJ_NAME_LEN) :: functional_name
+      integer :: ierr
+
+      no_functionals = option_count("/adjoint/functional")
+      do functional=0,no_functionals-1
+        if (have_option("/adjoint/functional[" // int2str(functional) // "]/functional_value")) then
+          call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+          ! Register the callback to compute J
+          ierr = adj_register_functional_callback(adjointer, trim(functional_name), c_funloc(libadjoint_evaluate_functional))
+          call adj_chkierr(ierr)
+        end if
+      end do
+    end subroutine register_functional_callbacks
+
+    subroutine calculate_functional_values(timestep)
+      integer, intent(in) :: timestep
+      integer :: functional, no_functionals
+      character(len=OPTION_PATH_LEN) :: functional_name
+      real, dimension(:), pointer :: fn_value
+      real :: J
+      integer :: ierr
+
+      no_functionals = option_count("/adjoint/functional")
+      do functional=0,no_functionals-1
+        if (have_option("/adjoint/functional[" // int2str(functional) // "]/functional_value")) then
+          call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+          ierr = adj_evaluate_functional(adjointer, timestep, functional_name, J)
+          call adj_chkierr(ierr)
+          ! So we've computed the component of the functional associated with this timestep.
+          ! We also want to sum them all up ...
+          call set_diagnostic(name=trim(functional_name) // "_component", statistic="value", value=(/J/))
+          fn_value => get_diagnostic(name=trim(functional_name), statistic="value")
+          J = J + fn_value(1)
+          call set_diagnostic(name=trim(functional_name), statistic="value", value=(/J/))
+        end if
+      end do
+    end subroutine calculate_functional_values
+  
+#endif
+
+    ! Register a diagnostic variable for each functional.
+    subroutine forward_main_loop_register_diagnostic
+      integer :: functional, no_functionals
+      character(len=OPTION_PATH_LEN) :: functional_name
+
+#ifdef HAVE_ADJOINT
+      no_functionals = option_count("/adjoint/functional")
+      
+      do functional=0,no_functionals-1
+        ! Register a diagnostic for each functional 
+        if (have_option("/adjoint/functional[" // int2str(functional) // "]/functional_value")) then
+          call get_option("/adjoint/functional[" // int2str(functional) // "]/name", functional_name)
+          call register_diagnostic(dim=1, name=trim(functional_name) // "_component", statistic="value")
+          call register_diagnostic(dim=1, name=trim(functional_name), statistic="value")
+          ! The functional value will be accumulated, so initialise it with zero.
+          call set_diagnostic(name=trim(functional_name) // "_component", statistic="value", value=(/0.0/))
+          call set_diagnostic(name=trim(functional_name), statistic="value", value=(/0.0/))
+        end if
+      end do
+#endif
+   end subroutine forward_main_loop_register_diagnostic
+   
 end module forward_main_loop
 
=== modified file 'adjoint/Libadjoint_Data_Callbacks.F90'
--- adjoint/Libadjoint_Data_Callbacks.F90	2011-05-24 10:42:33 +0000
+++ adjoint/Libadjoint_Data_Callbacks.F90	2011-07-12 15:56:32 +0000
@@ -88,6 +88,8 @@
     call adj_chkierr(ierr)
     ierr = adj_register_data_callback(adjointer, ADJ_VEC_GET_NORM_CB, c_funloc(femtools_vec_getnorm_proc))
     call adj_chkierr(ierr)
+    ierr = adj_register_data_callback(adjointer, ADJ_VEC_GET_SIZE_CB, c_funloc(femtools_vec_get_size_proc))
+    call adj_chkierr(ierr)
     ierr = adj_register_data_callback(adjointer, ADJ_VEC_DOT_PRODUCT_CB, c_funloc(femtools_vec_dot_product_proc))
     call adj_chkierr(ierr)
     ierr = adj_register_data_callback(adjointer, ADJ_VEC_SET_RANDOM_CB, c_funloc(femtools_vec_set_random_proc))
@@ -317,6 +319,27 @@
     end if
   end subroutine femtools_vec_set_random_proc
 
+  subroutine femtools_vec_get_size_proc(vec, sz) bind(c)
+    type(adj_vector), intent(in), value :: vec
+    integer(kind=c_int), intent(out) :: sz
+    type(scalar_field), pointer :: scalar
+    type(vector_field), pointer :: vector
+    type(tensor_field), pointer :: tensor
+
+    if (vec%klass == ADJ_SCALAR_FIELD) then
+      call c_f_pointer(vec%ptr, scalar)
+      sz = node_count(scalar)
+    else if (vec%klass == ADJ_VECTOR_FIELD) then
+      call c_f_pointer(vec%ptr, vector)
+      sz = node_count(vector) * vector%dim
+    else if (vec%klass == ADJ_TENSOR_FIELD) then
+      call c_f_pointer(vec%ptr, tensor)
+      sz = node_count(tensor) * tensor%dim(1) * tensor%dim(2)
+    else
+      FLAbort("adj_vector class not supported.")
+    end if
+  end subroutine femtools_vec_get_size_proc
+
   subroutine femtools_mat_axpy_proc(Y, alpha, X) bind(c)  
     ! Computes Y = alpha*X + Y.
     use iso_c_binding
 
=== modified file 'adjoint/Makefile.dependencies'
--- adjoint/Makefile.dependencies	2011-05-23 14:44:04 +0000
+++ adjoint/Makefile.dependencies	2011-07-12 15:56:32 +0000
@@ -1,4 +1,12 @@
 # Dependencies generated by create_makefile.py. DO NOT EDIT
+../include/adjoint_controls.mod: Adjoint_Controls.o
+	@true
+
+Adjoint_Controls.o ../include/adjoint_controls.mod: Adjoint_Controls.F90 \
+   ../include/boundary_conditions.mod ../include/global_parameters.mod \
+   ../include/fields.mod ../include/spud.mod ../include/python_state.mod \
+   ../include/fdebug.h ../include/state_module.mod
+
 ../include/adjoint_functional_evaluation.mod: Adjoint_Functional_Evaluation.o
 	@true
 
@@ -22,7 +30,7 @@
 
 Adjoint_Main_Loop.o ../include/adjoint_main_loop.mod: Adjoint_Main_Loop.F90 \
    ../include/populate_state_module.mod ../include/sparse_matrices_fields.mod \
-   ../include/diagnostic_fields_new.mod \
+   ../include/diagnostic_fields_new.mod ../include/adjoint_controls.mod \
    ../include/adjoint_functional_evaluation.mod ../include/sparse_tools.mod \
    ../include/fields.mod ../include/diagnostic_variables.mod \
    ../include/mangle_options_tree.mod ../include/write_state_module.mod \
@@ -44,12 +52,25 @@
 
 Burgers_Adjoint_Callbacks.o ../include/burgers_adjoint_callbacks.mod: \
    Burgers_Adjoint_Callbacks.F90 ../include/libadjoint_data_callbacks.mod \
-   ../include/adjoint_global_variables.mod ../include/fields.mod \
-   ../include/spud.mod ../include/mangle_options_tree.mod \
-   ../include/sparse_matrices_fields.mod \
-   ../include/mangle_dirichlet_rows_module.mod ../include/fdebug.h \
+   ../include/simple_advection.mod ../include/burgers_assembly.mod \
+   ../include/populate_state_module.mod \
+   ../include/adjoint_global_variables.mod ../include/global_parameters.mod \
+   ../include/fields.mod ../include/simple_advection_b.mod ../include/spud.mod \
+   ../include/mangle_options_tree.mod ../include/sparse_matrices_fields.mod \
+   ../include/mangle_dirichlet_rows_module.mod \
+   ../include/simple_advection_d.mod ../include/fdebug.h \
    ../include/state_module.mod
 
+../include/burgers_adjoint_controls.mod: Burgers_Control_Callbacks.o
+	@true
+
+Burgers_Control_Callbacks.o ../include/burgers_adjoint_controls.mod: \
+   Burgers_Control_Callbacks.F90 ../include/mangle_dirichlet_rows_module.mod \
+   ../include/manifold_projections.mod ../include/global_parameters.mod \
+   ../include/fields.mod ../include/spud.mod \
+   ../include/sparse_matrices_fields.mod ../include/state_module.mod \
+   ../include/python_state.mod ../include/fdebug.h
+
 ../include/forward_main_loop.mod: Forward_Main_Loop.o
 	@true
 
@@ -58,6 +79,7 @@
    ../include/diagnostic_fields_new.mod \
    ../include/adjoint_functional_evaluation.mod ../include/sparse_tools.mod \
    ../include/fields.mod ../include/diagnostic_variables.mod \
+   ../include/boundary_conditions_from_options.mod \
    ../include/mangle_options_tree.mod ../include/write_state_module.mod \
    ../include/spud.mod ../include/mangle_dirichlet_rows_module.mod \
    ../include/state_module.mod ../include/libadjoint_data_callbacks.mod \
@@ -105,3 +127,14 @@
    ../include/manifold_projections.mod ../include/fdebug.h \
    ../include/state_module.mod
 
+../include/shallow_water_adjoint_controls.mod: \
+   Shallow_Water_Control_Callbacks.o
+	@true
+
+Shallow_Water_Control_Callbacks.o \
+   ../include/shallow_water_adjoint_controls.mod: \
+   Shallow_Water_Control_Callbacks.F90 ../include/manifold_projections.mod \
+   ../include/global_parameters.mod ../include/fields.mod ../include/spud.mod \
+   ../include/sparse_matrices_fields.mod ../include/python_state.mod \
+   ../include/fdebug.h ../include/state_module.mod
+
 
=== modified file 'adjoint/Makefile.in'
--- adjoint/Makefile.in	2011-05-19 14:29:56 +0000
+++ adjoint/Makefile.in	2011-07-12 15:56:32 +0000
@@ -43,10 +43,16 @@
 
 LIB = ../lib/lib@FLUIDITY@.a
 
-OBJS    = $(patsubst %.F90,%.o,$(wildcard *.F90)) $(patsubst %.c,%.o,$(wildcard *.c))
-
-.SUFFIXES: .F90 .c .cpp .o .a
-
+OBJS    = $(patsubst %.F90,%.o,$(wildcard *.F90)) $(patsubst %.f90,%.o,$(wildcard *.f90)) $(patsubst %.c,%.o,$(wildcard *.c)) $(patsubst %.f,%.o,$(wildcard *.f)) 
+
+.SUFFIXES: .f .f90 .F90 .c .cpp .o .a
+
+.f90.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) $(GENFLAGS) -c $<
+.f.o:
+	@echo "        FC $<"
+	$(FC) $(FCFLAGS) $(GENFLAGS) -c $<
 .F90.o:
 	@echo "        FC $<"
 	$(FC) $(FCFLAGS) $(GENFLAGS) -c $< 
@@ -71,4 +77,16 @@
 clean:
 	rm -f *.o *.d *.mod
 
+# Instructions for doing AD with tapenade
+reverse:
+	tapenade -b -head "advection_action" -outvars "ac" -vars "u" simple_advection.f90 -ext PUSHPOPGeneralLib -extAD PUSHPOPADLib-fixinterface
+	echo "Don't forget to remove the zeroing of acb."
+
+forward:
+	tapenade -d -head "advection_action" -outvars "ac" -vars "u" simple_advection.f90
+
+../include/simple_advection.mod: simple_advection.o
+../include/simple_advection_d.mod: simple_advection_d.o
+../include/simple_advection_b.mod: simple_advection_b.o
+
 include Makefile.dependencies
 
=== modified file 'adjoint/Mangle_Dirichlet_Rows_Module.F90'
--- adjoint/Mangle_Dirichlet_Rows_Module.F90	2011-03-31 17:44:57 +0000
+++ adjoint/Mangle_Dirichlet_Rows_Module.F90	2011-07-12 15:56:32 +0000
@@ -35,7 +35,7 @@
   implicit none
 
   interface mangle_dirichlet_rows
-    module procedure mangle_dirichlet_rows_csr_scalar
+    module procedure mangle_dirichlet_rows_csr_scalar, mangle_dirichlet_rows_scalar
   end interface
 
   interface set_inactive_rows
@@ -90,6 +90,28 @@
         end do
     end subroutine mangle_dirichlet_rows_csr_scalar
 
+    subroutine mangle_dirichlet_rows_scalar(scalar, bc_field)
+        type(scalar_field), intent(inout) :: scalar
+        type(scalar_field), intent(in) :: bc_field
+
+        integer :: i, j, node
+        character(len=FIELD_NAME_LEN) :: bctype
+        integer, dimension(:), pointer :: node_list
+
+        assert(node_count(scalar) == node_count(bc_field))
+
+        do i=1,get_boundary_condition_count(bc_field)
+          call get_boundary_condition(bc_field, i, type=bctype, surface_node_list=node_list)
+
+          if (bctype /= "dirichlet") cycle
+
+          do j=1,size(node_list)
+            node = node_list(j)
+            call set(scalar, node, 0.0)
+          end do
+        end do
+    end subroutine mangle_dirichlet_rows_scalar
+
     subroutine set_inactive_rows_csr_scalar(matrix, field)
       ! Any nodes associated with Dirichlet BCs, we make them inactive
       type(csr_matrix), intent(inout) :: matrix
 
=== modified file 'adjoint/Mangle_Options_Tree.F90'
--- adjoint/Mangle_Options_Tree.F90	2011-03-10 10:07:19 +0000
+++ adjoint/Mangle_Options_Tree.F90	2011-07-12 15:56:32 +0000
@@ -139,12 +139,6 @@
       deallocate(fields_to_delete)
     end do
 
-    ! We delete this because if it exists in the forward run,
-    ! write_diagnostics wants to print it out in the stat file.
-    do functional=0,option_count("/adjoint/functional")-1
-      call delete_option("/adjoint/functional[" // int2str(functional) // "]/functional_value", stat=stat)
-    end do
-
   end subroutine mangle_options_tree_forward
 
   subroutine mangle_options_tree_adjoint
@@ -262,8 +256,13 @@
         path = "/material_phase[" // int2str(state) // "]/scalar_field::"
         call move_option(trim(path) // trim(fields_to_move(field)), trim(path) // "Adjoint" // trim(fields_to_move(field)), stat=stat)
         assert(stat == SPUD_NO_ERROR)
-        ! Delete any sources specified for the forward model -- we deal with these ourselves
-        call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/scalar_field::Source", stat=stat)
+        ! Delete the sources specified for the forward model -- if desired
+        if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/scalar_field::Source/")) then
+          if (have_option(trim(complete_field_path(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/scalar_field::Source/")) &
+                        & // "adjoint_storage/exists_in_forward")) then
+            call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/scalar_field::Source", stat=stat)
+          end if
+        end if
         ! And delete any initial conditions -- we will also specify these
         if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition")) then
           call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition", stat=stat)
@@ -294,7 +293,12 @@
         call move_option(trim(path) // trim(fields_to_move(field)), trim(path) // "Adjoint" // trim(fields_to_move(field)), stat=stat)
         assert(stat == SPUD_NO_ERROR)
         ! Delete any sources specified for the forward model -- we deal with these ourselves
-        call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/vector_field::Source", stat=stat)
+        if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/vector_field::Source/")) then
+          if (have_option(trim(complete_field_path(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/vector_field::Source/")) &
+                        & // "adjoint_storage/exists_in_forward")) then
+            call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/vector_field::Source", stat=stat)
+          end if
+        end if
         ! And delete any initial conditions -- we will also specify these
         if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition")) then
           call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition", stat=stat)
@@ -325,7 +329,12 @@
         call move_option(trim(path) // trim(fields_to_move(field)), trim(path) // "Adjoint" // trim(fields_to_move(field)), stat=stat)
         assert(stat == SPUD_NO_ERROR)
         ! Delete any sources specified for the forward model -- we deal with these ourselves
-        call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/tensor_field::Source", stat=stat)
+        if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/tensor_field::Source/")) then
+          if (have_option(trim(complete_field_path(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/tensor_field::Source/")) &
+                        & // "adjoint_storage/exists_in_forward")) then
+            call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/tensor_field::Source", stat=stat)
+          end if
+        end if
         ! And delete any initial conditions -- we will also specify these
         if (have_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition")) then
           call delete_option(trim(path) // "Adjoint" // trim(fields_to_move(field)) // "/prognostic/initial_condition", stat=stat)
 
=== modified file 'adjoint/Shallow_Water_Adjoint_Callbacks.F90'
--- adjoint/Shallow_Water_Adjoint_Callbacks.F90	2011-05-13 14:06:12 +0000
+++ adjoint/Shallow_Water_Adjoint_Callbacks.F90	2011-07-12 15:56:32 +0000
@@ -982,5 +982,6 @@
         end select
       end if
     end subroutine shallow_water_forward_source
+    
 #endif
 end module shallow_water_adjoint_callbacks
 
=== added file 'adjoint/Shallow_Water_Control_Callbacks.F90'
--- adjoint/Shallow_Water_Control_Callbacks.F90	1970-01-01 00:00:00 +0000
+++ adjoint/Shallow_Water_Control_Callbacks.F90	2011-07-12 15:56:32 +0000
@@ -0,0 +1,203 @@
+!    Copyright (C) 2006 Imperial College London and others.
+!
+!    Please see the AUTHORS file in the main source directory for a full list
+!    of copyright holders.
+!
+!    Prof. C Pain
+!    Applied Modelling and Computation Group
+!    Department of Earth Science and Engineering
+!    Imperial College London
+!
+!    [e-mail]
+!
+!    This library is free software; you can redistribute it and/or
+!    modify it under the terms of the GNU Lesser General Public
+!    License as published by the Free Software Foundation,
+!    version 2.1 of the License.
+!
+!    This library is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!    Lesser General Public License for more details.
+!
+!    You should have received a copy of the GNU Lesser General Public
+!    License along with this library; if not, write to the Free Software
+!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+!    USA
+
+#include "fdebug.h"
+
+module shallow_water_adjoint_controls
+  use fields
+  use global_parameters, only : PYTHON_FUNC_LEN, OPTION_PATH_LEN
+  use spud
+  use state_module
+  use python_state
+  use sparse_matrices_fields
+  use manifold_projections
+
+    implicit none
+
+    private
+
+    public :: shallow_water_adjoint_timestep_callback 
+
+    contains 
+    
+    ! Data is a void pointer used to pass variables into the callback
+    subroutine shallow_water_adjoint_timestep_callback(states, timestep, functional_name, data)
+      use iso_c_binding, only: c_ptr
+      use global_parameters, only: OPTION_PATH_LEN
+      use state_module
+      type(state_type), dimension(:), intent(inout) :: states
+      integer, intent(in) :: timestep
+      character(len=*), intent(in) :: functional_name
+      type(c_ptr), intent(in) :: data
+      
+      type(state_type), pointer :: matrices
+      character(len=OPTION_PATH_LEN) :: field_name, control_type, control_deriv_name, field_deriv_name, name, material_phase_name
+      integer :: nb_controls
+      integer :: i, state_id, s_idx
+      type(scalar_field), pointer :: sfield, adj_sfield, adj_sfield2
+      type(vector_field), pointer :: vfield, adj_vfield, positions
+      type(tensor_field), pointer :: tfield, adj_tfield
+      type(block_csr_matrix), pointer :: big_mat, div_mat, u_mass_mat
+      type(csr_matrix), pointer :: h_mass_mat
+      
+      type(vector_field) :: local_src_tmp, local_src_tmp2, src_tmp
+      real :: theta, d0, dt 
+     
+      ! Cast the data to the matrices state 
+      call c_f_pointer(data, matrices)
+      
+
+      nb_controls = option_count("/adjoint/controls/control")
+      do i = 0, nb_controls-1
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/name", control_deriv_name)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type/name", control_type)
+        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type::" // trim(control_type) // "/field_name", name)
+        s_idx = scan(trim(name), "::")
+        if (s_idx == 0) then 
+          FLAbort("The control " // trim(control_deriv_name) // " uses an invalid field_name. It should be of the form Materialphase::Fieldname")
+        end if
+        material_phase_name = name(1:s_idx - 1)
+        field_name = name(s_idx + 2:len_trim(name))
+        ! Find state associated with the material phase
+        do state_id = 1, size(states) 
+          if (trim(states(state_id)%name) == trim(material_phase_name)) then
+            exit
+          end if
+        end do
+        ! Make sure we found the state
+        if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
+          FLAbort("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_deriv_name) // ".")
+        end if
+        
+        control_deriv_name = trim(control_deriv_name) // "_TotalDerivative"
+        select case(trim(control_type))
+          !!!!!!!!!!!!! Initial condition !!!!!!!!!!!!
+          case ("initial_condition")
+            if (timestep /= 0) then
+             cycle
+            end if 
+            field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+            if (has_scalar_field(states(state_id), field_deriv_name)) then
+              if (trim(field_name) == "LayerThickness") then
+                adj_sfield => extract_scalar_field(states(state_id), "Adjoint" // trim(field_name))
+                sfield => extract_scalar_field(states(state_id), field_deriv_name) ! Output field
+                h_mass_mat => extract_csr_matrix(matrices, "LayerThicknessMassMatrix")
+                call mult(sfield, h_mass_mat, adj_sfield)
+              else
+                FLAbort("Sorry, I do not know how to compute the intial condition control for " // trim(field_name) // ".")
+              end if
+            elseif (has_vector_field(states(state_id), field_deriv_name)) then
+              if (trim(field_name) == "Velocity") then
+                adj_vfield => extract_vector_field(states(state_id), "Adjoint" // trim(field_name))
+                vfield => extract_vector_field(states(state_id), field_deriv_name) ! Output field
+                u_mass_mat => extract_block_csr_matrix(matrices, "CartesianVelocityMassMatrix")
+                call mult(vfield, u_mass_mat, adj_vfield)
+              else
+                FLAbort("Sorry, I do not know how to compute the intial condition control for " // trim(field_name) // ".")
+              end if
+            elseif (has_tensor_field(states(state_id), field_deriv_name)) then
+              tfield => extract_tensor_field(states(state_id), field_deriv_name) ! Output field
+              FLAbort("Sorry, I do not know how to compute the intial condition control for " // trim(field_name) // ".")
+            else
+              FLAbort("The control derivative field " // trim(field_deriv_name) // " specified for " // trim(control_deriv_name) // " is not a field in the state.")
+            end if
+          !!!!!!!!!!!!! Boundary condition !!!!!!!!!!!!
+          case("boundary_condition")
+            FLAbort("Boundary condition control not implemented yet.")
+          !!!!!!!!!!!!! Source !!!!!!!!!!!!
+          case("source_term")
+            if (timestep == 0) then
+             cycle
+            end if 
+            field_deriv_name = trim(functional_name) // "_" // control_deriv_name 
+            call get_option("/timestepping/timestep", dt)
+            if (has_scalar_field(states(state_id), field_deriv_name)) then
+              if (trim(field_name) == "LayerThicknessSource") then
+                adj_sfield => extract_scalar_field(states(state_id), "AdjointLayerThicknessDelta")
+                sfield => extract_scalar_field(states(state_id), field_deriv_name) ! Output field
+                h_mass_mat => extract_csr_matrix(matrices, "LayerThicknessMassMatrix")
+                call mult_T(sfield, h_mass_mat, adj_sfield)
+                call scale(sfield, abs(dt))
+              else
+                FLAbort("Sorry, I do not know how to compute the source condition control for " // trim(field_name) // ".")
+              end if
+            elseif (has_vector_field(states(state_id), field_deriv_name)) then
+              if (trim(field_name) == "VelocitySource") then
+                ! We want to compute:
+                ! \lambda_{delta \eta} \Delta t^2 d_0 \theta C^T (M+\theta \Delta t L)^{-1} M P_{Cart to local}
+                vfield => extract_vector_field(states(state_id), field_deriv_name) ! Output field
+                adj_sfield => extract_scalar_field(states(state_id), "AdjointLayerThicknessDelta")
+                adj_vfield => extract_vector_field(states(state_id), "AdjointLocalVelocityDelta") 
+                ! We need the AdjointVelocity for the option path only
+                adj_sfield2 => extract_scalar_field(states(state_id), "AdjointLayerThickness")
+                call get_option(trim(adj_sfield2%option_path) // "/prognostic/temporal_discretisation/theta", theta)
+                call get_option(trim(adj_sfield2%option_path) // "/prognostic/mean_layer_thickness", d0)
+                
+                big_mat => extract_block_csr_matrix(matrices, "InverseBigMatrix")
+                div_mat => extract_block_csr_matrix(matrices, "DivergenceMatrix")
+                h_mass_mat => extract_csr_matrix(matrices, "LayerThicknessMassMatrix")
+                u_mass_mat => extract_block_csr_matrix(matrices, "LocalVelocityMassMatrix")
+                positions => extract_vector_field(matrices, "Coordinate")
+
+                call allocate(local_src_tmp, adj_vfield%dim, adj_vfield%mesh, "TemporaryLocalVelocity")
+                call allocate(local_src_tmp2, adj_vfield%dim, adj_vfield%mesh, "TemporaryLocalVelocity2")
+                call allocate(src_tmp, vfield%dim, vfield%mesh, "TemporaryVelocity")
+                call mult_T(local_src_tmp, div_mat, adj_sfield)
+                call mult_T(local_src_tmp2, big_mat, local_src_tmp)
+                call mult_T(local_src_tmp, u_mass_mat, local_src_tmp2)
+                call project_cartesian_to_local(positions, src_tmp, local_src_tmp, transpose=.true.)
+                call scale(src_tmp, factor = dt*dt * d0 * theta)
+                call zero(vfield)
+                ! TODO: Where does this minus come from?
+                call addto(vfield, src_tmp, scale=-1.0)
+               
+                ! Now add the part from \lambda_{\Delta u}
+                ! \lambda_{\Delta u} M  M_{big} M P_{cart to local}
+                call mult_T(local_src_tmp, u_mass_mat, adj_vfield)
+                call mult_T(local_src_tmp2, big_mat, local_src_tmp)
+                call mult_T(local_src_tmp, u_mass_mat, local_src_tmp2)
+                call project_cartesian_to_local(positions, src_tmp, local_src_tmp, transpose = .true.)
+                call addto(vfield, src_tmp, scale = abs(dt))
+
+                call deallocate(local_src_tmp)
+                call deallocate(local_src_tmp2)
+                call deallocate(src_tmp)
+              else
+                FLAbort("Sorry, I do not know how to compute the source condition control for " // trim(field_name) // ".")
+              end if
+            elseif (has_tensor_field(states(state_id), field_deriv_name)) then
+              tfield => extract_tensor_field(states(state_id), field_deriv_name) ! Output field
+              FLAbort("Sorry, I do not know how to compute the source condition control for " // trim(field_name) // ".")
+            else
+              FLAbort("The control derivative field " // trim(field_deriv_name) // " specified for " // trim(control_deriv_name) // " is not a field in the state.")
+            end if
+        end select
+      end do
+
+    end subroutine shallow_water_adjoint_timestep_callback 
+
+end module shallow_water_adjoint_controls
 
=== added file 'adjoint/adBuffer.c'
--- adjoint/adBuffer.c	1970-01-01 00:00:00 +0000
+++ adjoint/adBuffer.c	2011-07-12 15:56:32 +0000
@@ -0,0 +1,476 @@
+static char adSid[]="$Id: adBuffer.c $";
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "adStack.h"
+
+/************ MEASUREMENT OF PUSH/POP TRAFFIC *************/
+
+static int mmftraffic = 0 ;
+static int mmftrafficM = 0 ;
+
+void addftraffic(int n) {
+  mmftraffic = mmftraffic+n ;
+  while (mmftraffic >= 1000000) {
+    mmftraffic = mmftraffic-1000000 ;
+    ++mmftrafficM ;
+  }
+  while (mmftraffic < 0) {
+    mmftraffic = mmftraffic+1000000 ;
+    --mmftraffic ;
+  }
+}
+
+void printtraffic() {
+  printctraffic_() ;
+  printf(" F Traffic: ") ;
+  printbigbytes(mmftrafficM, 1000000, mmftraffic) ;
+  printf(" bytes\n") ;
+}
+
+/************************** integer*4 ************************/
+static int adi4buf[512] ;
+static int adi4ibuf = 0 ;
+static int adi4lbuf[512] ;
+static int adi4ilbuf = -1 ;
+static int adi4inlbuf = 0 ;
+
+void pushinteger4(int x) {
+  addftraffic(4) ;
+  if (adi4ilbuf != -1) {
+    adi4ilbuf = -1 ;
+    adi4inlbuf = 0 ;
+  }
+  if (adi4ibuf >= 511) {
+    adi4buf[511] = x ;
+    pushNarray((char*)adi4buf, 512*4) ;
+    addftraffic(-512*4) ;
+    adi4ibuf = 0 ;
+  } else {
+    adi4buf[adi4ibuf] = x ;
+    ++adi4ibuf ;
+  }
+}
+
+void lookinteger4(int *x) {
+  if (adi4ilbuf == -1) {
+    adi4ilbuf = adi4ibuf ;
+    resetadlookstack_() ;
+  }
+  if (adi4ilbuf <= 0) {
+    lookNarray((char*)adi4lbuf, 512*4) ;
+    adi4inlbuf = 1 ;
+    adi4ilbuf = 511 ;
+    *x = adi4lbuf[511] ;
+  } else {
+    --adi4ilbuf ;
+    if (adi4inlbuf)
+      *x = adi4lbuf[adi4ilbuf] ;
+    else
+      *x = adi4buf[adi4ilbuf] ;
+  }
+}
+
+void popinteger4(int *x) {
+  if (adi4ilbuf != -1) {
+    adi4ilbuf = -1 ;
+    adi4inlbuf = 0 ;
+  }
+  if (adi4ibuf <= 0) {
+    popNarray((char*)adi4buf, 512*4) ;
+    adi4ibuf = 511 ;
+    *x = adi4buf[511] ;
+  } else {
+    --adi4ibuf ;
+    *x = adi4buf[adi4ibuf] ;
+  }
+}
+
+/*************************** bits *************************/
+static unsigned int adbitbuf = 0 ;
+static int adbitibuf = 1 ;
+static unsigned int adbitlbuf = 0 ;
+static int adbitilbuf = -1 ;
+static int adbitinlbuf = 0 ;
+
+void pushbit(int bit) {
+  if (adbitilbuf != -1) {
+    adbitilbuf = -1 ;
+    adbitinlbuf = 0 ;
+  }
+  adbitbuf<<=1 ;
+  if (bit) adbitbuf++ ;
+  if (adbitibuf>=32) {
+    pushinteger4(adbitbuf) ;
+    adbitbuf = 0 ;
+    adbitibuf = 1 ;
+  } else
+    adbitibuf++ ;
+}
+
+int lookbit() {
+  if (adbitilbuf==-1) {
+    adbitilbuf=adbitibuf ;
+    adbitlbuf = adbitbuf ;
+  }
+  if (adbitilbuf<=1) {
+    lookinteger4(&adbitlbuf) ;
+    adbitilbuf = 32 ;
+  } else
+    adbitilbuf-- ;
+  int bit = adbitlbuf%2 ;
+  adbitlbuf>>=1 ;
+  return bit ;
+}
+
+int popbit() {
+  if (adbitilbuf != -1) {
+    adbitilbuf = -1 ;
+    adbitinlbuf = 0 ;
+  }
+  if (adbitibuf<=1) {
+    popinteger4(&adbitbuf) ;
+    adbitibuf = 32 ;
+  } else
+    adbitibuf-- ;
+  int bit = adbitbuf%2 ;
+  adbitbuf>>=1 ;
+  return bit ;
+}
+
+/************************* controls ***********************/
+
+void pushcontrol1b(int cc) {
+  pushbit(cc) ;
+}
+
+void popcontrol1b(int *cc) {
+  *cc = (popbit()?1:0) ;
+}
+
+void lookcontrol1b(int *cc) {
+  *cc = (lookbit()?1:0) ;
+}
+
+void pushcontrol2b(int cc) {
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc) ;
+}
+
+void popcontrol2b(int *cc) {
+  *cc = (popbit()?2:0) ;
+  if (popbit()) (*cc)++ ;
+}
+
+void lookcontrol2b(int *cc) {
+  *cc = (lookbit()?2:0) ;
+  if (lookbit()) (*cc)++ ;
+}
+
+void pushcontrol3b(int cc) {
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc) ;
+}
+
+void popcontrol3b(int *cc) {
+  *cc = (popbit()?2:0) ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+}
+
+void lookcontrol3b(int *cc) {
+  *cc = (lookbit()?2:0) ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+}
+
+void pushcontrol4b(int cc) {
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc) ;
+}
+
+void popcontrol4b(int *cc) {
+  *cc = (popbit()?2:0) ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+}
+
+void lookcontrol4b(int *cc) {
+  *cc = (lookbit()?2:0) ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+}
+
+void pushcontrol5b(int cc) {
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc%2) ;
+  cc >>= 1 ;
+  pushbit(cc) ;
+}
+
+void popcontrol5b(int *cc) {
+  *cc = (popbit()?2:0) ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (popbit()) (*cc)++ ;
+}
+
+void lookcontrol5b(int *cc) {
+  *cc = (lookbit()?2:0) ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+  (*cc) <<= 1 ;
+  if (lookbit()) (*cc)++ ;
+}
+
+/************************** real*4 ************************/
+static float adr4buf[512] ;
+static int adr4ibuf = 0 ;
+static float adr4lbuf[512] ;
+static int adr4ilbuf = -1 ;
+static int adr4inlbuf = 0 ;
+
+void pushreal4(float x) {
+  addftraffic(4) ;
+  if (adr4ilbuf != -1) {
+    adr4ilbuf = -1 ;
+    adr4inlbuf = 0 ;
+  }
+  if (adr4ibuf >= 511) {
+    adr4buf[511] = x ;
+    pushNarray((char*)adr4buf, 512*4) ;
+    addftraffic(-512*4) ;
+    adr4ibuf = 0 ;
+  } else {
+    adr4buf[adr4ibuf] = x ;
+    ++adr4ibuf ;
+  }
+}
+
+void lookreal4(float *x) {
+  if (adr4ilbuf == -1) {
+    adr4ilbuf = adr4ibuf ;
+    resetadlookstack_() ;
+  }
+  if (adr4ilbuf <= 0) {
+    lookNarray((char*)adr4lbuf, 512*4) ;
+    adr4inlbuf = 1 ;
+    adr4ilbuf = 511 ;
+    *x = adr4lbuf[511] ;
+  } else {
+    --adr4ilbuf ;
+    if (adr4inlbuf)
+      *x = adr4lbuf[adr4ilbuf] ;
+    else
+      *x = adr4buf[adr4ilbuf] ;
+  }
+}
+
+void popreal4(float *x) {
+  if (adr4ilbuf != -1) {
+    adr4ilbuf = -1 ;
+    adr4inlbuf = 0 ;
+  }
+  if (adr4ibuf <= 0) {
+    popNarray((char*)adr4buf, 512*4) ;
+    adr4ibuf = 511 ;
+    *x = adr4buf[511] ;
+  } else {
+    --adr4ibuf ;
+    *x = adr4buf[adr4ibuf] ;
+  }
+}
+
+/************************** real*8 ************************/
+static double adr8buf[512] ;
+static int adr8ibuf = 0 ;
+static double adr8lbuf[512] ;
+static int adr8ilbuf = -1 ;
+static int adr8inlbuf = 0 ;
+
+void pushreal8(double x) {
+  addftraffic(8) ;
+  if (adr8ilbuf != -1) {
+    adr8ilbuf = -1 ;
+    adr8inlbuf = 0 ;
+  }
+  if (adr8ibuf >= 511) {
+    adr8buf[511] = x ;
+    pushNarray((char*)adr8buf, 512*8) ;
+    addftraffic(-4096) ;
+    adr8ibuf = 0 ;
+  } else {
+    adr8buf[adr8ibuf] = x ;
+    ++adr8ibuf ;
+  }
+}
+
+void lookreal8(double *x) {
+  if (adr8ilbuf == -1) {
+    adr8ilbuf = adr8ibuf ;
+    resetadlookstack_() ;
+  }
+  if (adr8ilbuf <= 0) {
+    lookNarray((char*)adr8lbuf, 512*8) ;
+    adr8inlbuf = 1 ;
+    adr8ilbuf = 511 ;
+    *x = adr8lbuf[511] ;
+  } else {
+    --adr8ilbuf ;
+    if (adr8inlbuf)
+      *x = adr8lbuf[adr8ilbuf] ;
+    else
+      *x = adr8buf[adr8ilbuf] ;
+  }
+}
+
+void popreal8(double *x) {
+  if (adr8ilbuf != -1) {
+    adr8ilbuf = -1 ;
+    adr8inlbuf = 0 ;
+  }
+  if (adr8ibuf <= 0) {
+    popNarray((char*)adr8buf, 512*8) ;
+    adr8ibuf = 511 ;
+    *x = adr8buf[511] ;
+  } else {
+    --adr8ibuf ;
+    *x = adr8buf[adr8ibuf] ;
+  }
+}
+
+/********* PRINTING THE SIZE OF STACKS AND BUFFERS ********/
+
+void printbuffertop() {
+  int size = 0 ;
+  size += adi4ibuf*4 ;
+  size += adr4ibuf*4 ;
+  size += adr8ibuf*8 ;
+  printf("Buffer size:%i bytes i.e. %i Kbytes\n",
+         size, (int) (size/1024.0)) ;
+}
+
+void showallstacks() {
+  int i ;
+  printf("BIT STACK      : %x == %i\n",adbitbuf,adbitbuf) ; 
+  printf("INTEGER*4 BUFFER[%i]:",adi4ibuf) ;
+  for (i=0 ; i<adi4ibuf ; ++i) printf(" %i",adi4buf[i]) ;
+  printf("\n") ;
+  printf("REAL*8 BUFFER[%i]:",adr8ibuf) ;
+  for (i=0 ; i<adr8ibuf ; ++i) printf(" %d",(int) adr8buf[i]) ;
+  printf("\n") ;
+  printf("REAL*4 BUFFER[%i]:",adr4ibuf) ;
+  for (i=0 ; i<adr4ibuf ; ++i) printf(" %f",adr4buf[i]) ;
+  printf("\n") ;
+  showrecentcstack_() ;
+}
+
+/**********************************************************
+ *        HOW TO CREATE PUSH* LOOK* POP* SUBROUTINES
+ *              YET FOR OTHER DATA TYPES
+ * Duplicate and uncomment the commented code below.
+ * In the duplicated code, replace:
+ *   ctct -> C type name (e.g. float double, int...)
+ *   TTTT -> BASIC TAPENADE TYPE NAME
+ *     (in character, boolean, integer, real, complex, pointer,...)
+ *   z7   -> LETTER-SIZE FOR TYPE
+ *     (in s,         b,       i,       r,    c,       p,      ...)
+ *   7    -> TYPE SIZE IN BYTES
+ * Don't forget to insert the corresponding lines in
+ * procedure printbuffertop(), otherwise the contribution of
+ * this new type to buffer occupation will not be seen.
+ * (not very important anyway...)
+ **********************************************************/
+
+/************************** TTTT*7 ************************/
+/*
+static ctct adz7buf[512] ;
+static int adz7ibuf = 0 ;
+static ctct adz7lbuf[512] ;
+static int adz7ilbuf = -1 ;
+static int adz7inlbuf = 0 ;
+
+void pushTTTT7(ctct x) {
+  addftraffic(7) ;
+  if (adz7ilbuf != -1) {
+    adz7ilbuf = -1 ;
+    adz7inlbuf = 0 ;
+  }
+  if (adz7ibuf >= 511) {
+    adz7buf[511] = x ;
+    pushNarray((char*)adz7buf, 512*7) ;
+    addftraffic(-512*7) ;
+    adz7ibuf = 0 ;
+  } else {
+    adz7buf[adz7ibuf] = x ;
+    ++adz7ibuf ;
+  }
+}
+
+void lookTTTT7(ctct *x) {
+  if (adz7ilbuf == -1) {
+    adz7ilbuf = adz7ibuf ;
+    resetadlookstack_() ;
+  }
+  if (adz7ilbuf <= 0) {
+    lookNarray((char*)adz7lbuf, 512*7) ;
+    adz7inlbuf = 1 ;
+    adz7ilbuf = 511 ;
+    *x = adz7lbuf[511] ;
+  } else {
+    --adz7ilbuf ;
+    if (adz7inlbuf)
+      *x = adz7lbuf[adz7ilbuf] ;
+    else
+      *x = adz7buf[adz7ilbuf] ;
+  }
+}
+
+void popTTTT7(ctct *x) {
+  if (adz7ilbuf != -1) {
+    adz7ilbuf = -1 ;
+    adz7inlbuf = 0 ;
+  }
+  if (adz7ibuf <= 0) {
+    popNarray((char*)adz7buf, 512*7) ;
+    adz7ibuf = 511 ;
+    *x = adz7buf[511] ;
+  } else {
+    --adz7ibuf ;
+    *x = adz7buf[adz7ibuf] ;
+  }
+}
+*/
+
+/**********************************************************/
 
=== added file 'adjoint/adBufferFortran.f'
--- adjoint/adBufferFortran.f	1970-01-01 00:00:00 +0000
+++ adjoint/adBufferFortran.f	2011-07-12 15:56:32 +0000
@@ -0,0 +1,1848 @@
+C$Id: adBuffer.f 3922 2011-05-19 08:54:39Z llh $
+
+c PISTES D'AMELIORATIONS:
+c  Attention aux IF qui peuvent couter cher.
+c  On pourrait aussi bufferiser les bits avec N entiers,
+c   (1 bit par entier), passer tout le paquet a C et laisser
+c   C faire les jongleries de bitsets.
+c  On pourrait aussi optimiser en -O3 les primitives de ADFirstAidKit
+c  Regarder l'assembleur (option -S (et -o toto.s))
+c  Pourchasser les divisions!
+
+      BLOCK DATA LOOKINGORNOT
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+      DATA looking/.FALSE./
+      END
+
+c======================== BITS ==========================:
+      BLOCK DATA BITS
+      INTEGER*4 adbitbuf, adbitlbuf
+      INTEGER adbitibuf, adbitilbuf
+      LOGICAL adbitinlbuf
+      COMMON /adbitfbuf/adbitbuf,adbitlbuf,
+     +       adbitibuf,adbitilbuf,adbitinlbuf
+      DATA adbitbuf/0/
+      DATA adbitlbuf/0/
+      DATA adbitibuf/0/
+      DATA adbitilbuf/-1/
+      DATA adbitinlbuf/.FALSE./
+      END
+
+c [0,31] are the bit indices we can use in an INTEGER
+
+      SUBROUTINE PUSHBIT(bit)
+      LOGICAL bit
+      INTEGER*4 adbitbuf, adbitlbuf
+      INTEGER adbitibuf, adbitilbuf
+      LOGICAL adbitinlbuf
+      COMMON /adbitfbuf/adbitbuf,adbitlbuf,
+     +       adbitibuf,adbitilbuf,adbitinlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adbitilbuf.ne.-1) THEN
+         adbitilbuf = -1
+         adbitinlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (bit) THEN
+         adbitbuf = IBSET(adbitbuf, adbitibuf)
+      ELSE
+         adbitbuf = IBCLR(adbitbuf, adbitibuf)
+      ENDIF
+      IF (adbitibuf.ge.31) THEN
+         CALL PUSHINTEGER4(adbitbuf)
+         adbitbuf = 0
+         adbitibuf = 0
+      ELSE
+         adbitibuf = adbitibuf+1
+      ENDIF
+      END
+
+      LOGICAL FUNCTION LOOKBIT()
+      INTEGER*4 adbitbuf, adbitlbuf
+      INTEGER adbitibuf, adbitilbuf
+      LOGICAL adbitinlbuf
+      COMMON /adbitfbuf/adbitbuf,adbitlbuf,
+     +       adbitibuf,adbitilbuf,adbitinlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adbitilbuf.eq.-1) THEN
+         adbitilbuf=adbitibuf
+         adbitlbuf = adbitbuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adbitilbuf.le.0) THEN
+         CALL LOOKINTEGER4(adbitlbuf)
+         adbitilbuf = 31
+      ELSE
+         adbitilbuf = adbitilbuf-1
+      ENDIF
+      LOOKBIT = BTEST(adbitlbuf, adbitilbuf)
+      END
+
+      LOGICAL FUNCTION POPBIT()
+      INTEGER*4 adbitbuf, adbitlbuf
+      INTEGER adbitibuf, adbitilbuf
+      LOGICAL adbitinlbuf
+      COMMON /adbitfbuf/adbitbuf,adbitlbuf,
+     +       adbitibuf,adbitilbuf,adbitinlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adbitilbuf.ne.-1) THEN
+         adbitilbuf = -1
+         adbitinlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adbitibuf.le.0) THEN
+         CALL POPINTEGER4(adbitbuf)
+         adbitibuf = 31
+      ELSE
+         adbitibuf = adbitibuf-1
+      ENDIF
+      POPBIT = BTEST(adbitbuf, adbitibuf)
+      END
+
+c====================== CONTROL =========================:
+
+      SUBROUTINE PUSHCONTROL1B(cc)
+      INTEGER cc
+      CALL PUSHBIT(cc.ne.0)
+      END
+
+      SUBROUTINE POPCONTROL1B(cc)
+      INTEGER cc
+      LOGICAL POPBIT
+      IF (POPBIT()) THEN
+         cc = 1
+      ELSE
+         cc = 0
+      ENDIF
+      END
+
+      SUBROUTINE LOOKCONTROL1B(cc)
+      INTEGER cc
+      LOGICAL LOOKBIT
+      IF (LOOKBIT()) THEN
+         cc = 1
+      ELSE
+         cc = 0
+      ENDIF
+      END
+
+      SUBROUTINE PUSHCONTROL2B(cc)
+      INTEGER cc
+      CALL PUSHBIT(BTEST(cc,0))
+      CALL PUSHBIT(BTEST(cc,1))
+      END
+
+      SUBROUTINE POPCONTROL2B(cc)
+      INTEGER cc
+      LOGICAL POPBIT
+      IF (POPBIT()) THEN
+         cc = 2
+      ELSE
+         cc = 0
+      ENDIF
+      IF (POPBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE LOOKCONTROL2B(cc)
+      INTEGER cc
+      LOGICAL LOOKBIT
+      IF (LOOKBIT()) THEN
+         cc = 2
+      ELSE
+         cc = 0
+      ENDIF
+      IF (LOOKBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE PUSHCONTROL3B(cc)
+      INTEGER cc
+      CALL PUSHBIT(BTEST(cc,0))
+      CALL PUSHBIT(BTEST(cc,1))
+      CALL PUSHBIT(BTEST(cc,2))
+      END
+
+      SUBROUTINE POPCONTROL3B(cc)
+      INTEGER cc
+      LOGICAL POPBIT
+      IF (POPBIT()) THEN
+         cc = 4
+      ELSE
+         cc = 0
+      ENDIF
+      IF (POPBIT()) cc = IBSET(cc,1)
+      IF (POPBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE LOOKCONTROL3B(cc)
+      INTEGER cc
+      LOGICAL LOOKBIT
+      IF (LOOKBIT()) THEN
+         cc = 4
+      ELSE
+         cc = 0
+      ENDIF
+      IF (LOOKBIT()) cc = IBSET(cc,1)
+      IF (LOOKBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE PUSHCONTROL4B(cc)
+      INTEGER cc
+      CALL PUSHBIT(BTEST(cc,0))
+      CALL PUSHBIT(BTEST(cc,1))
+      CALL PUSHBIT(BTEST(cc,2))
+      CALL PUSHBIT(BTEST(cc,3))
+      END
+
+      SUBROUTINE POPCONTROL4B(cc)
+      INTEGER cc
+      LOGICAL POPBIT
+      IF (POPBIT()) THEN
+         cc = 8
+      ELSE
+         cc = 0
+      ENDIF
+      IF (POPBIT()) cc = IBSET(cc,2)
+      IF (POPBIT()) cc = IBSET(cc,1)
+      IF (POPBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE LOOKCONTROL4B(cc)
+      INTEGER cc
+      LOGICAL LOOKBIT
+      IF (LOOKBIT()) THEN
+         cc = 8
+      ELSE
+         cc = 0
+      ENDIF
+      IF (LOOKBIT()) cc = IBSET(cc,2)
+      IF (LOOKBIT()) cc = IBSET(cc,1)
+      IF (LOOKBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE PUSHCONTROL5B(cc)
+      INTEGER cc
+      CALL PUSHBIT(BTEST(cc,0))
+      CALL PUSHBIT(BTEST(cc,1))
+      CALL PUSHBIT(BTEST(cc,2))
+      CALL PUSHBIT(BTEST(cc,3))
+      CALL PUSHBIT(BTEST(cc,4))
+      END
+
+      SUBROUTINE POPCONTROL5B(cc)
+      INTEGER cc
+      LOGICAL POPBIT
+      IF (POPBIT()) THEN
+         cc = 16
+      ELSE
+         cc = 0
+      ENDIF
+      IF (POPBIT()) cc = IBSET(cc,3)
+      IF (POPBIT()) cc = IBSET(cc,2)
+      IF (POPBIT()) cc = IBSET(cc,1)
+      IF (POPBIT()) cc = IBSET(cc,0)
+      END
+
+      SUBROUTINE LOOKCONTROL5B(cc)
+      INTEGER cc
+      LOGICAL LOOKBIT
+      IF (LOOKBIT()) THEN
+         cc = 16
+      ELSE
+         cc = 0
+      ENDIF
+      IF (LOOKBIT()) cc = IBSET(cc,3)
+      IF (LOOKBIT()) cc = IBSET(cc,2)
+      IF (LOOKBIT()) cc = IBSET(cc,1)
+      IF (LOOKBIT()) cc = IBSET(cc,0)
+      END
+
+c======================= BOOLEANS =========================
+
+      SUBROUTINE PUSHBOOLEAN(x)
+      LOGICAL x
+      CALL PUSHBIT(x)
+      END
+
+      SUBROUTINE LOOKBOOLEAN(x)
+      LOGICAL x, LOOKBIT
+      x = LOOKBIT()
+      END
+
+      SUBROUTINE POPBOOLEAN(x)
+      LOGICAL x, POPBIT
+      x = POPBIT()
+      END
+
+c===================== CHARACTERS =======================:
+      BLOCK DATA CHARACTERS
+      CHARACTER ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+      DATA ads1ibuf/1/
+      DATA ads1ilbuf/-1/
+      DATA ads1inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHCHARACTER(x)
+      CHARACTER x, ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(1)
+      IF (ads1ilbuf.ne.-1) THEN
+         ads1ilbuf = -1
+         ads1inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (ads1ibuf.ge.512) THEN
+         ads1buf(512) = x
+         CALL PUSHCHARACTERARRAY(ads1buf, 512)
+      CALL addftraffic(-512)
+         ads1ibuf = 1
+      ELSE
+         ads1buf(ads1ibuf) = x
+         ads1ibuf = ads1ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKCHARACTER(x)
+      CHARACTER x, ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (ads1ilbuf.eq.-1) THEN
+         ads1ilbuf=ads1ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (ads1ilbuf.le.1) THEN
+         CALL LOOKCHARACTERARRAY(ads1lbuf, 512)
+         ads1inlbuf = .TRUE.
+         ads1ilbuf = 512
+         x = ads1lbuf(512)
+      ELSE
+         ads1ilbuf = ads1ilbuf-1
+         if (ads1inlbuf) THEN
+            x = ads1lbuf(ads1ilbuf)
+         ELSE
+            x = ads1buf(ads1ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPCHARACTER(x)
+      CHARACTER x, ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (ads1ilbuf.ne.-1) THEN
+         ads1ilbuf = -1
+         ads1inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (ads1ibuf.le.1) THEN
+         CALL POPCHARACTERARRAY(ads1buf, 512)
+         ads1ibuf = 512
+         x = ads1buf(512)
+      ELSE
+         ads1ibuf = ads1ibuf-1
+         x = ads1buf(ads1ibuf)
+      ENDIF
+      END
+
+c======================= INTEGER*4 =========================:
+      BLOCK DATA INTEGERS4
+      INTEGER*4 adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      DATA adi4ibuf/1/
+      DATA adi4ilbuf/-1/
+      DATA adi4inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHINTEGER4(x)
+      INTEGER*4 x, adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(4)
+      IF (adi4ilbuf.ne.-1) THEN
+         adi4ilbuf = -1
+         adi4inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adi4ibuf.ge.512) THEN
+         adi4buf(512) = x
+         CALL PUSHINTEGER4ARRAY(adi4buf, 512)
+         CALL addftraffic(-2048)
+         adi4ibuf = 1
+      ELSE
+         adi4buf(adi4ibuf) = x
+         adi4ibuf = adi4ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKINTEGER4(x)
+      INTEGER*4 x, adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adi4ilbuf.eq.-1) THEN
+         adi4ilbuf=adi4ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adi4ilbuf.le.1) THEN
+         CALL LOOKINTEGER4ARRAY(adi4lbuf, 512)
+         adi4inlbuf = .TRUE.
+         adi4ilbuf = 512
+         x = adi4lbuf(512)
+      ELSE
+         adi4ilbuf = adi4ilbuf-1
+         if (adi4inlbuf) THEN
+            x = adi4lbuf(adi4ilbuf)
+         ELSE
+            x = adi4buf(adi4ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPINTEGER4(x)
+      INTEGER*4 x, adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adi4ilbuf.ne.-1) THEN
+         adi4ilbuf = -1
+         adi4inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adi4ibuf.le.1) THEN
+         CALL POPINTEGER4ARRAY(adi4buf, 512)
+         adi4ibuf = 512
+         x = adi4buf(512)
+      ELSE
+         adi4ibuf = adi4ibuf-1
+         x = adi4buf(adi4ibuf)
+      ENDIF
+      END
+
+c======================= INTEGER*8 =========================
+      BLOCK DATA INTEGERS8
+      INTEGER*8 adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+      DATA adi8ibuf/1/
+      DATA adi8ilbuf/-1/
+      DATA adi8inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHINTEGER8(x)
+      INTEGER*8 x, adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(8)
+      IF (adi8ilbuf.ne.-1) THEN
+         adi8ilbuf = -1
+         adi8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adi8ibuf.ge.512) THEN
+         adi8buf(512) = x
+         CALL PUSHINTEGER8ARRAY(adi8buf, 512)
+      CALL addftraffic(-4096)
+         adi8ibuf = 1
+      ELSE
+         adi8buf(adi8ibuf) = x
+         adi8ibuf = adi8ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKINTEGER8(x)
+      INTEGER*8 x, adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adi8ilbuf.eq.-1) THEN
+         adi8ilbuf=adi8ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adi8ilbuf.le.1) THEN
+         CALL LOOKINTEGER8ARRAY(adi8lbuf, 512)
+         adi8inlbuf = .TRUE.
+         adi8ilbuf = 512
+         x = adi8lbuf(512)
+      ELSE
+         adi8ilbuf = adi8ilbuf-1
+         if (adi8inlbuf) THEN
+            x = adi8lbuf(adi8ilbuf)
+         ELSE
+            x = adi8buf(adi8ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPINTEGER8(x)
+      INTEGER*8 x, adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adi8ilbuf.ne.-1) THEN
+         adi8ilbuf = -1
+         adi8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adi8ibuf.le.1) THEN
+         CALL POPINTEGER8ARRAY(adi8buf, 512)
+         adi8ibuf = 512
+         x = adi8buf(512)
+      ELSE
+         adi8ibuf = adi8ibuf-1
+         x = adi8buf(adi8ibuf)
+      ENDIF
+      END
+
+c======================= REAL*4 =========================
+      BLOCK DATA REALS4
+      REAL*4 adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      DATA adr4ibuf/1/
+      DATA adr4ilbuf/-1/
+      DATA adr4inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHREAL4(x)
+      REAL*4 x, adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(4)
+      IF (adr4ilbuf.ne.-1) THEN
+         adr4ilbuf = -1
+         adr4inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr4ibuf.ge.512) THEN
+         adr4buf(512) = x
+         CALL PUSHREAL4ARRAY(adr4buf, 512)
+      CALL addftraffic(-2048)
+         adr4ibuf = 1
+      ELSE
+         adr4buf(adr4ibuf) = x
+         adr4ibuf = adr4ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKREAL4(x)
+      REAL*4 x, adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr4ilbuf.eq.-1) THEN
+         adr4ilbuf=adr4ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adr4ilbuf.le.1) THEN
+         CALL LOOKREAL4ARRAY(adr4lbuf, 512)
+         adr4inlbuf = .TRUE.
+         adr4ilbuf = 512
+         x = adr4lbuf(512)
+      ELSE
+         adr4ilbuf = adr4ilbuf-1
+         if (adr4inlbuf) THEN
+            x = adr4lbuf(adr4ilbuf)
+         ELSE
+            x = adr4buf(adr4ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPREAL4(x)
+      REAL*4 x, adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr4ilbuf.ne.-1) THEN
+         adr4ilbuf = -1
+         adr4inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr4ibuf.le.1) THEN
+         CALL POPREAL4ARRAY(adr4buf, 512)
+         adr4ibuf = 512
+         x = adr4buf(512)
+      ELSE
+         adr4ibuf = adr4ibuf-1
+         x = adr4buf(adr4ibuf)
+      ENDIF
+      END
+
+c======================= REAL*8 =========================
+      BLOCK DATA REALS8
+      REAL*8 adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      DATA adr8ibuf/1/
+      DATA adr8ilbuf/-1/
+      DATA adr8inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHREAL8(x)
+      REAL*8 x, adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(8)
+      IF (adr8ilbuf.ne.-1) THEN
+         adr8ilbuf = -1
+         adr8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr8ibuf.ge.512) THEN
+         adr8buf(512) = x
+         CALL PUSHREAL8ARRAY(adr8buf, 512)
+      CALL addftraffic(-4096)
+         adr8ibuf = 1
+      ELSE
+         adr8buf(adr8ibuf) = x
+         adr8ibuf = adr8ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKREAL8(x)
+      REAL*8 x, adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr8ilbuf.eq.-1) THEN
+         adr8ilbuf=adr8ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adr8ilbuf.le.1) THEN
+         CALL LOOKREAL8ARRAY(adr8lbuf, 512)
+         adr8inlbuf = .TRUE.
+         adr8ilbuf = 512
+         x = adr8lbuf(512)
+      ELSE
+         adr8ilbuf = adr8ilbuf-1
+         if (adr8inlbuf) THEN
+            x = adr8lbuf(adr8ilbuf)
+         ELSE
+            x = adr8buf(adr8ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPREAL8(x)
+      REAL*8 x, adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr8ilbuf.ne.-1) THEN
+         adr8ilbuf = -1
+         adr8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr8ibuf.le.1) THEN
+         CALL POPREAL8ARRAY(adr8buf, 512)
+         adr8ibuf = 512
+         x = adr8buf(512)
+      ELSE
+         adr8ibuf = adr8ibuf-1
+         x = adr8buf(adr8ibuf)
+      ENDIF
+      END
+
+c======================= REAL*16 =========================
+      BLOCK DATA REALS16
+      DOUBLE PRECISION adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+      DATA adr16ibuf/1/
+      DATA adr16ilbuf/-1/
+      DATA adr16inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHREAL16(x)
+      DOUBLE PRECISION x, adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(16)
+      IF (adr16ilbuf.ne.-1) THEN
+         adr16ilbuf = -1
+         adr16inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr16ibuf.ge.512) THEN
+         adr16buf(512) = x
+         CALL PUSHREAL16ARRAY(adr16buf, 512)
+      CALL addftraffic(-8192)
+         adr16ibuf = 1
+      ELSE
+         adr16buf(adr16ibuf) = x
+         adr16ibuf = adr16ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKREAL16(x)
+      DOUBLE PRECISION x, adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr16ilbuf.eq.-1) THEN
+         adr16ilbuf=adr16ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adr16ilbuf.le.1) THEN
+         CALL LOOKREAL16ARRAY(adr16lbuf, 512)
+         adr16inlbuf = .TRUE.
+         adr16ilbuf = 512
+         x = adr16lbuf(512)
+      ELSE
+         adr16ilbuf = adr16ilbuf-1
+         if (adr16inlbuf) THEN
+            x = adr16lbuf(adr16ilbuf)
+         ELSE
+            x = adr16buf(adr16ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPREAL16(x)
+      DOUBLE PRECISION x, adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adr16ilbuf.ne.-1) THEN
+         adr16ilbuf = -1
+         adr16inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adr16ibuf.le.1) THEN
+         CALL POPREAL16ARRAY(adr16buf, 512)
+         adr16ibuf = 512
+         x = adr16buf(512)
+      ELSE
+         adr16ibuf = adr16ibuf-1
+         x = adr16buf(adr16ibuf)
+      ENDIF
+      END
+
+c======================= COMPLEX*8 =========================
+      BLOCK DATA COMPLEXS8
+      COMPLEX*8 adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      DATA adc8ibuf/1/
+      DATA adc8ilbuf/-1/
+      DATA adc8inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHCOMPLEX8(x)
+      COMPLEX*8 x, adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(8)
+      IF (adc8ilbuf.ne.-1) THEN
+         adc8ilbuf = -1
+         adc8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adc8ibuf.ge.512) THEN
+         adc8buf(512) = x
+         CALL PUSHCOMPLEX8ARRAY(adc8buf, 512)
+      CALL addftraffic(-4096)
+         adc8ibuf = 1
+      ELSE
+         adc8buf(adc8ibuf) = x
+         adc8ibuf = adc8ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKCOMPLEX8(x)
+      COMPLEX*8 x, adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adc8ilbuf.eq.-1) THEN
+         adc8ilbuf=adc8ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adc8ilbuf.le.1) THEN
+         CALL LOOKCOMPLEX8ARRAY(adc8lbuf, 512)
+         adc8inlbuf = .TRUE.
+         adc8ilbuf = 512
+         x = adc8lbuf(512)
+      ELSE
+         adc8ilbuf = adc8ilbuf-1
+         if (adc8inlbuf) THEN
+            x = adc8lbuf(adc8ilbuf)
+         ELSE
+            x = adc8buf(adc8ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPCOMPLEX8(x)
+      COMPLEX*8 x, adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adc8ilbuf.ne.-1) THEN
+         adc8ilbuf = -1
+         adc8inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adc8ibuf.le.1) THEN
+         CALL POPCOMPLEX8ARRAY(adc8buf, 512)
+         adc8ibuf = 512
+         x = adc8buf(512)
+      ELSE
+         adc8ibuf = adc8ibuf-1
+         x = adc8buf(adc8ibuf)
+      ENDIF
+      END
+
+c======================= COMPLEX*16 =========================
+      BLOCK DATA COMPLEXS16
+      COMPLEX*16 adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+      DATA adc16ibuf/1/
+      DATA adc16ilbuf/-1/
+      DATA adc16inlbuf/.FALSE./
+      END
+
+      SUBROUTINE PUSHCOMPLEX16(x)
+      COMPLEX*16 x, adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      CALL addftraffic(16)
+      IF (adc16ilbuf.ne.-1) THEN
+         adc16ilbuf = -1
+         adc16inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adc16ibuf.ge.512) THEN
+         adc16buf(512) = x
+         CALL PUSHCOMPLEX16ARRAY(adc16buf, 512)
+      CALL addftraffic(-8192)
+         adc16ibuf = 1
+      ELSE
+         adc16buf(adc16ibuf) = x
+         adc16ibuf = adc16ibuf+1
+      ENDIF
+      END
+
+      SUBROUTINE LOOKCOMPLEX16(x)
+      COMPLEX*16 x, adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adc16ilbuf.eq.-1) THEN
+         adc16ilbuf=adc16ibuf
+         IF (.not.looking) THEN
+            CALL RESETADLOOKSTACK()
+            looking = .TRUE.
+         ENDIF
+      ENDIF
+      IF (adc16ilbuf.le.1) THEN
+         CALL LOOKCOMPLEX16ARRAY(adc16lbuf, 512)
+         adc16inlbuf = .TRUE.
+         adc16ilbuf = 512
+         x = adc16lbuf(512)
+      ELSE
+         adc16ilbuf = adc16ilbuf-1
+         if (adc16inlbuf) THEN
+            x = adc16lbuf(adc16ilbuf)
+         ELSE
+            x = adc16buf(adc16ilbuf)
+         ENDIF
+      ENDIF
+      END
+
+      SUBROUTINE POPCOMPLEX16(x)
+      COMPLEX*16 x, adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+      LOGICAL looking
+      COMMON /lookingfbuf/looking
+c
+      IF (adc16ilbuf.ne.-1) THEN
+         adc16ilbuf = -1
+         adc16inlbuf = .FALSE.
+         looking = .FALSE.
+      ENDIF
+      IF (adc16ibuf.le.1) THEN
+         CALL POPCOMPLEX16ARRAY(adc16buf, 512)
+         adc16ibuf = 512
+         x = adc16buf(512)
+      ELSE
+         adc16ibuf = adc16ibuf-1
+         x = adc16buf(adc16ibuf)
+      ENDIF
+      END
+
+C=========== MEASUREMENT OF PUSH/POP TRAFFIC ==========
+
+      BLOCK DATA MEMTRAFFIC
+      INTEGER*8 mmftraffic,mmftrafficM
+      COMMON /mmcomtraffic/mmftraffic,mmftrafficM
+      DATA mmftraffic/0/
+      DATA mmftrafficM/0/
+      END
+
+      subroutine addftraffic(n)
+      INTEGER n
+      INTEGER*8 mmftraffic,mmftrafficM
+      COMMON /mmcomtraffic/mmftraffic,mmftrafficM
+c
+      mmftraffic = mmftraffic+n
+      if (mmftraffic.ge.1000000) then
+ 100     mmftraffic = mmftraffic-1000000
+         mmftrafficM = mmftrafficM+1
+         if (mmftraffic.ge.1000000) then
+            goto 100
+         else
+            goto 300
+         endif
+      else if (mmftraffic.lt.0) then
+ 200     mmftraffic = mmftraffic+1000000
+         mmftrafficM = mmftrafficM-1
+         if (mmftraffic.lt.0) then
+            goto 200
+         else
+            goto 300
+         endif
+      endif
+ 300  continue
+      END
+
+      SUBROUTINE PRINTTRAFFIC()
+      INTEGER*8 mmftraffic,mmftrafficM
+      COMMON /mmcomtraffic/mmftraffic,mmftrafficM
+      CALL printctraffic()
+      CALL printftrafficinc(mmftrafficM, 1000000, mmftraffic)
+c       write (6,1001) ' F Traffic: ',mmftrafficM,' Mb and ',
+c      +     (((mmftraffic*1000)/1024)*1000)/1024, ' millionths'
+c 1001  format(a,i6,a,i6,a)
+      END
+
+C ============ PRINTING THE SIZE OF STACKS AND BUFFERS ==========
+
+      SUBROUTINE PRINTBUFFERTOP()
+      integer*4 SMALLSTACKSIZE
+      integer*4 size
+
+      size = SMALLSTACKSIZE()
+      print *,'Buffer size:',size,' bytes i.e. ',size/1024.0,' Kbytes'
+      END
+
+      FUNCTION SMALLSTACKSIZE()
+      CHARACTER ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+c      LOGICAL adl4buf(512), adl4lbuf(512)
+c      INTEGER adl4ibuf,adl4ilbuf
+c      LOGICAL adl4inlbuf
+c      COMMON /adl4fbuf/adl4buf,adl4lbuf,
+c     +       adl4ibuf,adl4ilbuf,adl4inlbuf
+      INTEGER*4 adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      INTEGER*8 adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+c           INTEGER*16 adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+      REAL*4 adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      REAL*8 adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      DOUBLE PRECISION adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+c           REAL*32 x, adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           COMPLEX*4 adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+      COMPLEX*8 adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      COMPLEX*16 adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+c           COMPLEX*32 adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+      integer*4 smallstacksize
+c     
+      smallstacksize = 0
+      smallstacksize = smallstacksize + (ads1ibuf-1)*1
+c           smallstacksize = smallstacksize + (adl4ibuf-1)*4
+      smallstacksize = smallstacksize + (adi4ibuf-1)*4
+      smallstacksize = smallstacksize + (adi8ibuf-1)*8
+c           smallstacksize = smallstacksize + (adi16ibuf-1)*16
+      smallstacksize = smallstacksize + (adr4ibuf-1)*4
+      smallstacksize = smallstacksize + (adr8ibuf-1)*8
+      smallstacksize = smallstacksize + (adr16ibuf-1)*16
+c           smallstacksize = smallstacksize + (adr32ibuf-1)*32
+c           smallstacksize = smallstacksize + (adc4ibuf-1)*4
+      smallstacksize = smallstacksize + (adc8ibuf-1)*8
+      smallstacksize = smallstacksize + (adc16ibuf-1)*16
+c           smallstacksize = smallstacksize + (adc32ibuf-1)*32
+c
+      end
+
+c Very complete display of the current size of the
+c  push/look/pop local Fortran stacks and global C stack.
+      SUBROUTINE PRINTALLBUFFERS()
+      CHARACTER ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+c      LOGICAL adl4buf(512), adl4lbuf(512)
+c      INTEGER adl4ibuf,adl4ilbuf
+c      LOGICAL adl4inlbuf
+c      COMMON /adl4fbuf/adl4buf,adl4lbuf,
+c     +       adl4ibuf,adl4ilbuf,adl4inlbuf
+      INTEGER*4 adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      INTEGER*8 adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+c           INTEGER*16 adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+      REAL*4 adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      REAL*8 adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      DOUBLE PRECISION adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+c           REAL*32 x, adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           COMPLEX*4 adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+      COMPLEX*8 adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      COMPLEX*16 adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+c           COMPLEX*32 adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+      integer*4 bsize,lookbsize
+      integer*4 cblocks, csize, lookcblocks, lookcsize
+c     
+      call getbigcsizes(cblocks,csize,lookcblocks,lookcsize)
+      write (6,'(a,i8,a,i5,a,i8,a,i5,a)')
+     +     'MAIN stack size (in C) :',cblocks,'B +',csize,
+     +     '  (looking:',lookcblocks,'B +',lookcsize,')'
+      bsize = (ads1ibuf-1)*1
+      lookbsize = -999
+      if (ads1inlbuf.or.ads1ilbuf.gt.-1) lookbsize=(ads1ilbuf-1)*1
+      write (6,'(a,i4,a,i4,a)') 'CHARs  :',bsize,
+     +     ' (looking:',lookbsize,')'
+c           bsize = (adl4ibuf-1)*4
+      bsize = (adi4ibuf-1)*4
+      lookbsize = -999
+      if (adi4inlbuf.or.adi4ilbuf.gt.-1) lookbsize=(adi4ilbuf-1)*4
+      write (6,'(a,i4,a,i4,a)') 'INTs4  :',bsize,
+     +     ' (looking:',lookbsize,')'
+      bsize = (adi8ibuf-1)*8
+      lookbsize = -999
+      if (adi8inlbuf.or.adi8ilbuf.gt.-1) lookbsize=(adi8ilbuf-1)*8
+      write (6,'(a,i4,a,i4,a)') 'INTs8  :',bsize,
+     +     ' (looking:',lookbsize,')'
+c           bsize = (adi16ibuf-1)*16
+      bsize = (adr4ibuf-1)*4
+      lookbsize = -999
+      if (adr4inlbuf.or.adr4ilbuf.gt.-1) lookbsize=(adr4ilbuf-1)*4
+      write (6,'(a,i4,a,i4,a)') 'REALs4 :',bsize,
+     +     ' (looking:',lookbsize,')'
+      bsize = (adr8ibuf-1)*8
+      lookbsize = -999
+      if (adr8inlbuf.or.adr8ilbuf.gt.-1) lookbsize=(adr8ilbuf-1)*8
+      write (6,'(a,i4,a,i4,a)') 'REALs8 :',bsize,
+     +     ' (looking:',lookbsize,')'
+      bsize = (adr16ibuf-1)*16
+      lookbsize = -999
+      if (adr16inlbuf.or.adr16ilbuf.gt.-1) lookbsize=(adr16ilbuf-1)*16
+      write (6,'(a,i4,a,i4,a)') 'REALs16:',bsize,
+     +     ' (looking:',lookbsize,')'
+c           bsize = (adr32ibuf-1)*32
+c           bsize = (adc4ibuf-1)*4
+      bsize = (adc8ibuf-1)*8
+      lookbsize = -999
+      if (adc8inlbuf.or.adc8ilbuf.gt.-1) lookbsize=(adc8ilbuf-1)*8
+      write (6,'(a,i4,a,i4,a)') 'CPLXs8 :',bsize,
+     +     ' (looking:',lookbsize,')'
+      bsize = (adc16ibuf-1)*16
+      lookbsize = -999
+      if (adc16inlbuf.or.adc16ilbuf.gt.-1) lookbsize=(adc16ilbuf-1)*16
+      write (6,'(a,i4,a,i4,a)') 'CPLXs16:',bsize,
+     +     ' (looking:',lookbsize,')'
+c           bsize = (adc32ibuf-1)*32
+c
+      end
+
+C FOR INTERNAL DEBUGS ONLY:
+      SUBROUTINE SHOWALLSTACKS()
+      INTEGER*4 adbitbuf, adbitlbuf
+      INTEGER adbitibuf, adbitilbuf
+      LOGICAL adbitinlbuf
+      COMMON /adbitfbuf/adbitbuf,adbitlbuf,
+     +       adbitibuf,adbitilbuf,adbitinlbuf
+      CHARACTER ads1buf(512), ads1lbuf(512)
+      INTEGER ads1ibuf,ads1ilbuf
+      LOGICAL ads1inlbuf
+      COMMON /ads1fbuf/ads1buf,ads1lbuf,
+     +       ads1ibuf,ads1ilbuf,ads1inlbuf
+      INTEGER*4 adi4buf(512), adi4lbuf(512)
+      INTEGER adi4ibuf,adi4ilbuf
+      LOGICAL adi4inlbuf
+      COMMON /adi4fbuf/adi4buf,adi4lbuf,
+     +       adi4ibuf,adi4ilbuf,adi4inlbuf
+      INTEGER*8 adi8buf(512), adi8lbuf(512)
+      INTEGER adi8ibuf,adi8ilbuf
+      LOGICAL adi8inlbuf
+      COMMON /adi8fbuf/adi8buf,adi8lbuf,
+     +       adi8ibuf,adi8ilbuf,adi8inlbuf
+      REAL*4 adr4buf(512), adr4lbuf(512)
+      INTEGER adr4ibuf,adr4ilbuf
+      LOGICAL adr4inlbuf
+      COMMON /adr4fbuf/adr4buf,adr4lbuf,
+     +       adr4ibuf,adr4ilbuf,adr4inlbuf
+      REAL*8 adr8buf(512), adr8lbuf(512)
+      INTEGER adr8ibuf,adr8ilbuf
+      LOGICAL adr8inlbuf
+      COMMON /adr8fbuf/adr8buf,adr8lbuf,
+     +       adr8ibuf,adr8ilbuf,adr8inlbuf
+      DOUBLE PRECISION adr16buf(512), adr16lbuf(512)
+      INTEGER adr16ibuf,adr16ilbuf
+      LOGICAL adr16inlbuf
+      COMMON /adr16fbuf/adr16buf,adr16lbuf,
+     +       adr16ibuf,adr16ilbuf,adr16inlbuf
+      COMPLEX*8 adc8buf(512), adc8lbuf(512)
+      INTEGER adc8ibuf,adc8ilbuf
+      LOGICAL adc8inlbuf
+      COMMON /adc8fbuf/adc8buf,adc8lbuf,
+     +       adc8ibuf,adc8ilbuf,adc8inlbuf
+      COMPLEX*16 adc16buf(512), adc16lbuf(512)
+      INTEGER adc16ibuf,adc16ilbuf
+      LOGICAL adc16inlbuf
+      COMMON /adc16fbuf/adc16buf,adc16lbuf,
+     +       adc16ibuf,adc16ilbuf,adc16inlbuf
+      INTEGER i
+c
+      write (6,1010) 'BIT STACK       : ',adbitbuf,'==',adbitbuf,
+     +     ' (',adbitibuf,')'
+1010  format(a,i20,a,z16,a,i2,a)
+      write (6,1011) 'INTEGER*8 BUFFER[',adi8ibuf-1,']: ',
+     +     (adi8buf(i),i=1,adi8ibuf-1)
+      write (6,1011) 'INTEGER*4 BUFFER[',adi4ibuf-1,']: ',
+     +     (adi4buf(i),i=1,adi4ibuf-1)
+1011  format(a,i3,a,512(i40))
+      write (6,1012) 'REAL*16   BUFFER:[',adr16ibuf-1,']: ',
+     +     (adr16buf(i),i=1,adr16ibuf-1)
+      write (6,1012) 'REAL*8    BUFFER:[',adr8ibuf-1, ']: ',
+     +     (adr8buf(i),i=1,adr8ibuf-1)
+      write (6,1012) 'REAL*4    BUFFER:[',adr4ibuf-1, ']: ',
+     +     (adr4buf(i),i=1,adr4ibuf-1)
+1012  format(a,i3,a,512(e8.2))
+      call showrecentcstack()
+c
+      END
+
+C========================================================
+C PUSH* POP* SUBROUTINES FOR OTHER DATA TYPES
+C Uncomment if these types are available on your compiler
+C  and they are needed by the reverse differentiated code
+C Don't forget to uncomment the corresponding lines in
+C  subroutine PRINTBUFFERTOP, otherwise these types'
+C  contribution to buffer occupation will not be seen.
+C  (not very important anyway...)
+
+c======================= INTEGER*16 =========================
+c           BLOCK DATA INTEGERS16
+c           INTEGER*16 adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+c           DATA adi16ibuf/1/
+c           DATA adi16ilbuf/-1/
+c           DATA adi16inlbuf/.FALSE./
+c           END
+c     c
+c           SUBROUTINE PUSHINTEGER16(x)
+c           INTEGER*16 x, adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           CALL addftraffic(16)
+c           IF (adi16ilbuf.ne.-1) THEN
+c              adi16ilbuf = -1
+c              adi16inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adi16ibuf.ge.512) THEN
+c              adi16buf(512) = x
+c              CALL PUSHINTEGER16ARRAY(adi16buf, 512)
+c           CALL addftraffic(-8192)
+c              adi16ibuf = 1
+c           ELSE
+c              adi16buf(adi16ibuf) = x
+c              adi16ibuf = adi16ibuf+1
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE LOOKINTEGER16(x)
+c           INTEGER*16 x, adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adi16ilbuf.eq.-1) THEN
+c              adi16ilbuf=adi16ibuf
+c              IF (.not.looking) THEN
+c                 CALL RESETADLOOKSTACK()
+c                looking = .TRUE.
+c              ENDIF
+c           ENDIF
+c           IF (adi16ilbuf.le.1) THEN
+c              CALL LOOKINTEGER16ARRAY(adi16lbuf, 512)
+c              adi16inlbuf = .TRUE.
+c              adi16ilbuf = 512
+c              x = adi16lbuf(512)
+c           ELSE
+c              adi16ilbuf = adi16ilbuf-1
+c              if (adi16inlbuf) THEN
+c                 x = adi16lbuf(adi16ilbuf)
+c              ELSE
+c                 x = adi16buf(adi16ilbuf)
+c              ENDIF
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE POPINTEGER16(x)
+c           INTEGER*16 x, adi16buf(512), adi16lbuf(512)
+c           INTEGER adi16ibuf,adi16ilbuf
+c           LOGICAL adi16inlbuf
+c           COMMON /adi16fbuf/adi16buf,adi16lbuf,
+c          +       adi16ibuf,adi16ilbuf,adi16inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adi16ilbuf.ne.-1) THEN
+c              adi16ilbuf = -1
+c              adi16inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adi16ibuf.le.1) THEN
+c              CALL POPINTEGER16ARRAY(adi16buf, 512)
+c              adi16ibuf = 512
+c              x = adi16buf(512)
+c           ELSE
+c              adi16ibuf = adi16ibuf-1
+c              x = adi16buf(adi16ibuf)
+c           ENDIF
+c           END
+
+c======================= REAL*32 =========================
+c           BLOCK DATA REALS32
+c           REAL*32 adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           DATA adr32ibuf/1/
+c           DATA adr32ilbuf/-1/
+c           DATA adr32inlbuf/.FALSE./
+c           END
+c     c
+c           SUBROUTINE PUSHREAL32(x)
+c           REAL*32 x, adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           CALL addftraffic(32)
+c           IF (adr32ilbuf.ne.-1) THEN
+c              adr32ilbuf = -1
+c              adr32inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adr32ibuf.ge.512) THEN
+c              adr32buf(512) = x
+c              CALL PUSHREAL32ARRAY(adr32buf, 512)
+c           CALL addftraffic(-16384)
+c              adr32ibuf = 1
+c           ELSE
+c              adr32buf(adr32ibuf) = x
+c              adr32ibuf = adr32ibuf+1
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE LOOKREAL32(x)
+c           REAL*32 x, adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adr32ilbuf.eq.-1) THEN
+c              adr32ilbuf=adr32ibuf
+c              IF (.not.looking) THEN
+c                 CALL RESETADLOOKSTACK()
+c                looking = .TRUE.
+c              ENDIF
+c           ENDIF
+c           IF (adr32ilbuf.le.1) THEN
+c              CALL LOOKREAL32ARRAY(adr32lbuf, 512)
+c              adr32inlbuf = .TRUE.
+c              adr32ilbuf = 512
+c              x = adr32lbuf(512)
+c           ELSE
+c              adr32ilbuf = adr32ilbuf-1
+c              if (adr32inlbuf) THEN
+c                 x = adr32lbuf(adr32ilbuf)
+c              ELSE
+c                 x = adr32buf(adr32ilbuf)
+c              ENDIF
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE POPREAL32(x)
+c           REAL*32 x, adr32buf(512), adr32lbuf(512)
+c           INTEGER adr32ibuf,adr32ilbuf
+c           LOGICAL adr32inlbuf
+c           COMMON /adr32fbuf/adr32buf,adr32lbuf,
+c          +       adr32ibuf,adr32ilbuf,adr32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adr32ilbuf.ne.-1) THEN
+c              adr32ilbuf = -1
+c              adr32inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adr32ibuf.le.1) THEN
+c              CALL POPREAL32ARRAY(adr32buf, 512)
+c              adr32ibuf = 512
+c              x = adr32buf(512)
+c           ELSE
+c              adr32ibuf = adr32ibuf-1
+c              x = adr32buf(adr32ibuf)
+c           ENDIF
+c           END
+
+c======================= COMPLEX*4 =========================
+c           BLOCK DATA COMPLEXS4
+c           COMPLEX*4 adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+c           DATA adc4ibuf/1/
+c           DATA adc4ilbuf/-1/
+c           DATA adc4inlbuf/.FALSE./
+c           END
+c     c
+c           SUBROUTINE PUSHCOMPLEX4(x)
+c           COMPLEX*4 x, adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           CALL addftraffic(4)
+c           IF (adc4ilbuf.ne.-1) THEN
+c              adc4ilbuf = -1
+c              adc4inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adc4ibuf.ge.512) THEN
+c              adc4buf(512) = x
+c              CALL PUSHCOMPLEX4ARRAY(adc4buf, 512)
+c           CALL addftraffic(-2048)
+c              adc4ibuf = 1
+c           ELSE
+c              adc4buf(adc4ibuf) = x
+c              adc4ibuf = adc4ibuf+1
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE LOOKCOMPLEX4(x)
+c           COMPLEX*4 x, adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adc4ilbuf.eq.-1) THEN
+c              adc4ilbuf=adc4ibuf
+c              IF (.not.looking) THEN
+c                 CALL RESETADLOOKSTACK()
+c                looking = .TRUE.
+c              ENDIF
+c           ENDIF
+c           IF (adc4ilbuf.le.1) THEN
+c              CALL LOOKCOMPLEX4ARRAY(adc4lbuf, 512)
+c              adc4inlbuf = .TRUE.
+c              adc4ilbuf = 512
+c              x = adc4lbuf(512)
+c           ELSE
+c              adc4ilbuf = adc4ilbuf-1
+c              if (adc4inlbuf) THEN
+c                 x = adc4lbuf(adc4ilbuf)
+c              ELSE
+c                 x = adc4buf(adc4ilbuf)
+c              ENDIF
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE POPCOMPLEX4(x)
+c           COMPLEX*4 x, adc4buf(512), adc4lbuf(512)
+c           INTEGER adc4ibuf,adc4ilbuf
+c           LOGICAL adc4inlbuf
+c           COMMON /adc4fbuf/adc4buf,adc4lbuf,
+c          +       adc4ibuf,adc4ilbuf,adc4inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adc4ilbuf.ne.-1) THEN
+c              adc4ilbuf = -1
+c              adc4inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adc4ibuf.le.1) THEN
+c              CALL POPCOMPLEX4ARRAY(adc4buf, 512)
+c              adc4ibuf = 512
+c              x = adc4buf(512)
+c           ELSE
+c              adc4ibuf = adc4ibuf-1
+c              x = adc4buf(adc4ibuf)
+c           ENDIF
+c           END
+
+c======================= COMPLEX*32 =========================
+c           BLOCK DATA COMPLEXS32
+c           COMPLEX*32 adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+c           DATA adc32ibuf/1/
+c           DATA adc32ilbuf/-1/
+c           DATA adc32inlbuf/.FALSE./
+c           END
+c     c
+c           SUBROUTINE PUSHCOMPLEX32(x)
+c           COMPLEX*32 x, adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           CALL addftraffic(32)
+c           IF (adc32ilbuf.ne.-1) THEN
+c              adc32ilbuf = -1
+c              adc32inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adc32ibuf.ge.512) THEN
+c              adc32buf(512) = x
+c              CALL PUSHCOMPLEX32ARRAY(adc32buf, 512)
+c           CALL addftraffic(-16384)
+c              adc32ibuf = 1
+c           ELSE
+c              adc32buf(adc32ibuf) = x
+c              adc32ibuf = adc32ibuf+1
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE LOOKCOMPLEX32(x)
+c           COMPLEX*32 x, adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adc32ilbuf.eq.-1) THEN
+c              adc32ilbuf=adc32ibuf
+c              IF (.not.looking) THEN
+c                 CALL RESETADLOOKSTACK()
+c                looking = .TRUE.
+c              ENDIF
+c           ENDIF
+c           IF (adc32ilbuf.le.1) THEN
+c              CALL LOOKCOMPLEX32ARRAY(adc32lbuf, 512)
+c              adc32inlbuf = .TRUE.
+c              adc32ilbuf = 512
+c              x = adc32lbuf(512)
+c           ELSE
+c              adc32ilbuf = adc32ilbuf-1
+c              if (adc32inlbuf) THEN
+c                 x = adc32lbuf(adc32ilbuf)
+c              ELSE
+c                 x = adc32buf(adc32ilbuf)
+c              ENDIF
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE POPCOMPLEX32(x)
+c           COMPLEX*32 x, adc32buf(512), adc32lbuf(512)
+c           INTEGER adc32ibuf,adc32ilbuf
+c           LOGICAL adc32inlbuf
+c           COMMON /adc32fbuf/adc32buf,adc32lbuf,
+c          +       adc32ibuf,adc32ilbuf,adc32inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adc32ilbuf.ne.-1) THEN
+c              adc32ilbuf = -1
+c              adc32inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adc32ibuf.le.1) THEN
+c              CALL POPCOMPLEX32ARRAY(adc32buf, 512)
+c              adc32ibuf = 512
+c              x = adc32buf(512)
+c           ELSE
+c              adc32ibuf = adc32ibuf-1
+c              x = adc32buf(adc32ibuf)
+c           ENDIF
+c           END
+
+C========================================================
+C        HOW TO CREATE PUSH* POP* SUBROUTINES
+C              YET FOR OTHER DATA TYPES
+C ** Duplicate the commented program lines below
+c ** In the duplicated subroutines, replace:
+c            TTTT by the basic name of the type
+c            z9   by the initial and size of the type
+c                (integer:i real:r complex:c boolean:b character:s)
+c            9    by the size of the type
+c ** Uncomment the duplicated subroutines
+C ** Don't forget to insert the corresponding lines in
+C  subroutine PRINTBUFFERTOP, otherwise these types'
+C  contribution to buffer occupation will not be seen.
+C  (not very important anyway...)
+
+c======================= TTTT*9 =========================
+c           BLOCK DATA TTTTS9
+c           TTTT*9 adz9buf(512), adz9lbuf(512)
+c           INTEGER adz9ibuf,adz9ilbuf
+c           LOGICAL adz9inlbuf
+c           COMMON /adz9fbuf/adz9buf,adz9lbuf,
+c          +       adz9ibuf,adz9ilbuf,adz9inlbuf
+c           DATA adz9ibuf/1/
+c           DATA adz9ilbuf/-1/
+c           DATA adz9inlbuf/.FALSE./
+c           END
+c     c
+c           SUBROUTINE PUSHTTTT9(x)
+c           TTTT*9 x, adz9buf(512), adz9lbuf(512)
+c           INTEGER adz9ibuf,adz9ilbuf
+c           LOGICAL adz9inlbuf
+c           COMMON /adz9fbuf/adz9buf,adz9lbuf,
+c          +       adz9ibuf,adz9ilbuf,adz9inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           CALL addftraffic(9)
+c           IF (adz9ilbuf.ne.-1) THEN
+c              adz9ilbuf = -1
+c              adz9inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adz9ibuf.ge.512) THEN
+c              adz9buf(512) = x
+c              CALL PUSHTTTT9ARRAY(adz9buf, 512)
+c           CALL addftraffic(-9*512)
+c              adz9ibuf = 1
+c           ELSE
+c              adz9buf(adz9ibuf) = x
+c              adz9ibuf = adz9ibuf+1
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE LOOKTTTT9(x)
+c           TTTT*9 x, adz9buf(512), adz9lbuf(512)
+c           INTEGER adz9ibuf,adz9ilbuf
+c           LOGICAL adz9inlbuf
+c           COMMON /adz9fbuf/adz9buf,adz9lbuf,
+c          +       adz9ibuf,adz9ilbuf,adz9inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adz9ilbuf.eq.-1) THEN
+c              adz9ilbuf=adz9ibuf
+c              IF (.not.looking) THEN
+c                 CALL RESETADLOOKSTACK()
+c                looking = .TRUE.
+c              ENDIF
+c           ENDIF
+c           IF (adz9ilbuf.le.1) THEN
+c              CALL LOOKTTTT9ARRAY(adz9lbuf, 512)
+c              adz9inlbuf = .TRUE.
+c              adz9ilbuf = 512
+c              x = adz9lbuf(512)
+c           ELSE
+c              adz9ilbuf = adz9ilbuf-1
+c              if (adz9inlbuf) THEN
+c                 x = adz9lbuf(adz9ilbuf)
+c              ELSE
+c                 x = adz9buf(adz9ilbuf)
+c              ENDIF
+c           ENDIF
+c           END
+c     
+c           SUBROUTINE POPTTTT9(x)
+c           TTTT*9 x, adz9buf(512), adz9lbuf(512)
+c           INTEGER adz9ibuf,adz9ilbuf
+c           LOGICAL adz9inlbuf
+c           COMMON /adz9fbuf/adz9buf,adz9lbuf,
+c          +       adz9ibuf,adz9ilbuf,adz9inlbuf
+c           LOGICAL looking
+c           COMMON /lookingfbuf/looking
+c     c
+c           IF (adz9ilbuf.ne.-1) THEN
+c              adz9ilbuf = -1
+c              adz9inlbuf = .FALSE.
+c              looking = .FALSE.
+c           ENDIF
+c           IF (adz9ibuf.le.1) THEN
+c              CALL POPTTTT9ARRAY(adz9buf, 512)
+c              adz9ibuf = 512
+c              x = adz9buf(512)
+c           ELSE
+c              adz9ibuf = adz9ibuf-1
+c              x = adz9buf(adz9ibuf)
+c           ENDIF
+c           END
 
=== added file 'adjoint/adStack.c'
--- adjoint/adStack.c	1970-01-01 00:00:00 +0000
+++ adjoint/adStack.c	2011-07-12 15:56:32 +0000
@@ -0,0 +1,683 @@
+static char adSid[]="$Id: adStack.c 3919 2011-05-18 15:51:18Z llh $";
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#define ONE_BLOCK_SIZE 16384
+#ifndef STACK_SIZE_TRACING
+#define STACK_SIZE_TRACING 1
+#endif
+/* The main stack is a double-chain of DoubleChainedBlock objects.
+ * Each DoubleChainedBlock holds an array[ONE_BLOCK_SIZE] of char. */
+typedef struct _doubleChainedBlock{
+  struct _doubleChainedBlock *prev ;
+  char                       *contents ;
+  struct _doubleChainedBlock *next ;
+} DoubleChainedBlock ;
+
+/* Globals that define the current position in the stack: */
+static DoubleChainedBlock *curStack = NULL ;
+static char               *curStackTop    = NULL ;
+/* Globals that define the current LOOKing position in the stack: */
+static DoubleChainedBlock *lookStack = NULL ;
+static char               *lookStackTop    = NULL ;
+
+static long int mmctraffic = 0 ;
+static long int mmctrafficM = 0 ;
+#ifdef STACK_SIZE_TRACING
+long int bigStackSize = 0;
+#endif
+
+/* PUSHes "nbChars" consecutive chars from a location starting at address "x".
+ * Resets the LOOKing position if it was active.
+ * Checks that there is enough space left to hold "nbChars" chars.
+ * Otherwise, allocates the necessary space. */
+void pushNarray(char *x, unsigned int nbChars) {
+  unsigned int nbmax = (curStack)?ONE_BLOCK_SIZE-(curStackTop-(curStack->contents)):0 ;
+#ifdef STACK_SIZE_TRACING
+  bigStackSize += nbChars;
+#endif
+
+  mmctraffic += nbChars ;
+  while (mmctraffic >= 1000000) {
+     mmctraffic -= 1000000 ;
+     mmctrafficM++ ;
+  }
+
+  lookStack = NULL ;
+  if (nbChars <= nbmax) {
+    memcpy(curStackTop,x,nbChars) ;
+    curStackTop+=nbChars ;
+  } else {
+    char *inx = x+(nbChars-nbmax) ;
+    if (nbmax>0) memcpy(curStackTop,inx,nbmax) ;
+    while (inx>x) {
+      if ((curStack == NULL) || (curStack->next == NULL)) {
+        /* Create new block: */
+	DoubleChainedBlock *newStack ;
+	char *contents = (char*)malloc(ONE_BLOCK_SIZE*sizeof(char)) ;
+	newStack = (DoubleChainedBlock*)malloc(sizeof(DoubleChainedBlock)) ;
+	if ((contents == NULL) || (newStack == NULL)) {
+	  DoubleChainedBlock *stack = curStack ;
+	  int nbBlocks = (stack?-1:0) ;
+	  while(stack) {
+	      stack = stack->prev ;
+	      nbBlocks++ ;
+	  }
+	  printf("Out of memory (allocated %i blocks of %i bytes)\n",
+		 nbBlocks, ONE_BLOCK_SIZE) ;
+          exit(0);
+	}
+	if (curStack != NULL) curStack->next = newStack ;
+	newStack->prev = curStack ;
+	newStack->next = NULL ;
+	newStack->contents = contents ;
+	curStack = newStack ;
+        /* new block created! */
+      } else
+	curStack = curStack->next ;
+      inx -= ONE_BLOCK_SIZE ;
+      if(inx>x)
+	memcpy(curStack->contents,inx,ONE_BLOCK_SIZE) ;
+      else {
+	unsigned int nbhead = (inx-x)+ONE_BLOCK_SIZE ;
+	curStackTop = curStack->contents ;
+	memcpy(curStackTop,x,nbhead) ;
+	curStackTop += nbhead ;
+      }
+    }
+  }
+}
+
+/* POPs "nbChars" consecutive chars to a location starting at address "x".
+ * Resets the LOOKing position if it was active.
+ * Checks that there is enough data to fill "nbChars" chars.
+ * Otherwise, pops as many blocks as necessary. */
+void popNarray(char *x, unsigned int nbChars) {
+  unsigned int nbmax = curStackTop-(curStack->contents) ;
+#ifdef STACK_SIZE_TRACING
+  bigStackSize -= nbChars;
+#endif
+  lookStack = NULL ;
+  if (nbChars <= nbmax) {
+    curStackTop-=nbChars ;
+    memcpy(x,curStackTop,nbChars);
+  } else {
+    char *tlx = x+nbChars ;
+    if (nbmax>0) memcpy(x,curStack->contents,nbmax) ;
+    x+=nbmax ;
+    while (x<tlx) {
+      curStack = curStack->prev ;
+      if (curStack==NULL) printf("Popping from an empty stack!!!") ;
+      if (x+ONE_BLOCK_SIZE<tlx) {
+	memcpy(x,curStack->contents,ONE_BLOCK_SIZE) ;
+	x += ONE_BLOCK_SIZE ;
+      } else {
+	unsigned int nbtail = tlx-x ;
+	curStackTop=(curStack->contents)+ONE_BLOCK_SIZE-nbtail ;
+	memcpy(x,curStackTop,nbtail) ;
+	x = tlx ;
+      }
+    }
+  }
+}
+
+/* LOOKs "nbChars" consecutive chars to a location starting at address "x".
+ * Activates the LOOKing position if it was reset.
+ * LOOKing is just like POPping, except that the main pointer
+ * remains in place, so that the value is not POPped.
+ * Further PUSHs or POPs will start from the same place as if
+ * no LOOK had been made. */
+void lookNarray(char *x, unsigned int nbChars) {
+  unsigned int nbmax ;
+  if (lookStack == NULL) {
+    lookStack = curStack ;
+    lookStackTop = curStackTop ;
+  }
+  nbmax = lookStackTop-(lookStack->contents) ;
+  if (nbChars <= nbmax) {
+    lookStackTop-=nbChars ;
+    memcpy(x,lookStackTop,nbChars);
+  } else {
+    char *tlx = x+nbChars ;
+    if (nbmax>0) memcpy(x,lookStack->contents,nbmax) ;
+    x+=nbmax ;
+    while (x<tlx) {
+      lookStack = lookStack->prev ;
+      if (lookStack==NULL) printf("Looking into an empty stack!!!") ;
+      if (x+ONE_BLOCK_SIZE<tlx) {
+	memcpy(x,lookStack->contents,ONE_BLOCK_SIZE) ;
+	x += ONE_BLOCK_SIZE ;
+      } else {
+	unsigned int nbtail = tlx-x ;
+	lookStackTop=(lookStack->contents)+ONE_BLOCK_SIZE-nbtail ;
+	memcpy(x,lookStackTop,nbtail) ;
+	x = tlx ;
+      }
+    }
+  }
+}
+
+void resetadlookstack_() {
+  lookStack=NULL ;
+}
+
+/****** Exported PUSH/POP/LOOK functions for ARRAYS: ******/
+/*   --> Called from FORTRAN:                             */
+
+void pushcharacterarray_(char *x, unsigned int *n) {
+  pushNarray(x,*n) ;
+}
+void popcharacterarray_(char *x, unsigned int *n) {
+  popNarray(x,*n) ;
+}
+void lookcharacterarray_(char *x, unsigned int *n) {
+  lookNarray(x,*n) ;
+}
+
+void pushbooleanarray_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*4)) ;
+}
+void popbooleanarray_(char *x, unsigned int *n) {
+  popNarray(x,(*n*4)) ;
+}
+void lookbooleanarray_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*4)) ;
+}
+
+void pushinteger4array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*4)) ;
+}
+void popinteger4array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*4)) ;
+}
+void lookinteger4array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*4)) ;
+}
+
+void pushinteger8array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*8)) ;
+}
+void popinteger8array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*8)) ;
+}
+void lookinteger8array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*8)) ;
+}
+
+void pushinteger16array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*16)) ;
+}
+void popinteger16array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*16)) ;
+}
+void lookinteger16array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*16)) ;
+}
+
+void pushreal4array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*4)) ;
+}
+void popreal4array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*4)) ;
+}
+void lookreal4array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*4)) ;
+}
+
+void pushreal8array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*8)) ;
+}
+void popreal8array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*8)) ;
+}
+void lookreal8array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*8)) ;
+}
+
+void pushreal16array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*16)) ;
+}
+void popreal16array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*16)) ;
+}
+void lookreal16array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*16)) ;
+}
+
+void pushreal32array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*32)) ;
+}
+void popreal32array_(char *x, unsigned int *n) {
+  popNarray(x,(*n*32)) ;
+}
+void lookreal32array_(char *x, unsigned int *n) {
+  lookNarray(x,(*n*32)) ;
+}
+
+void pushcomplex4array_(char *x, unsigned int *n) {
+  pushNarray(x,(*n*4)) ;