!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS 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 General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! NAME
!     SUBROUTINE create_global_domain_structure
!
! DESCRIPTION
!   The processes implemented in this subroutine calculate the matrices which are required to 
!   characterise the multi-conductor cable propagation in both a local, domain based system for the Spice model development
!   and in a global basis for analysis of the bundle as a whole for the validation test circuit.
!   The matrices characterising the propagation are :
!       the domain based 'high frequency' L and C matrices,
!       the domain based frequency dependent Z and Y matrices,
!       the global  'high frequency' L and C matrices,
!       the global frequency dependent Z and Y matrices,
!       the voltage and current domain decomposition matrices MI and MV
!
!   For shielded domains the L, C, Z, and Y matrices are known from the cable data structures
!   For any overshield domains and the external domain these matrices are calculated here via
!   calls to the appopriate PUL calculation subroutines
!
!   The processes implemented in this subroutine are summarised as follows:
!
! 1. Count the total number of conductors in the bundle and Allocate and reset the local data structures 
! 2. Loop over the cables filling the referencing arrays with local information
! 3. For each cable, work out which domain it is in and hence its reference conductor number.
! 4. Count the number of conductors in the external domain to check whether it is viable i.e. at least 2 conductors
! 5. build the global_domain numbering, ignoring any unused domains
! 6. build the domain based L and C matrices plus the global MI and MV matrices and the global L and C matrices
! 7. Calculate the inductance and capacitance matrices for the external domains (domains within overshields and the external domain)
! 8. Copy the cable based conductor impedance (loss) models to the bundle structure
! 8b. Copy the cable based conductor labels to the bundle structure
! 9  save numbering information required for the transfer impedance calculation
! 10. deallocate local memory
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!     27/4/2016 CJS: Include a conductor based impedance (loss) model
!     24/5/2016 CJS: Fix errors in MI and MV matrices for multi-conductor overshield domains
!     7/7/2016 CJS: Allow the use of the Laplace solver for L,C,G matrix calculation in external and overshield domains.
!     9/9/2016 CJS: Fix problem with MV and MI matrix for the twisted pair model. Elements were being overwritten incorrectly.
!                   when adding the external conductor contribution to the domain decomposition matrices   
!     4/10/2016 CJS: Save the local domain conductor numbering
!     17/5/2017 CJS: Work out which domains are Shielded Twisted Pair differential mode domains and flag these
!                    so that we can work out the is_shield flag properly in all circumstances.
!     18/10/2017 CJS: include 8b. Copy the cable based conductor labels to the bundle structure
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!     16/3/2018 CJS  add y offset for ML_flex_cable
!     22/10/2019 CJS Fix bug with numbering here: overshields were numbered using terminal conductor number not global domain number
!     24/10/2019 CJS Allow mutliple layers of overshields (previously overshields could not be 'nested' 
!                    (though they could cover shielded cables).
!

SUBROUTINE create_global_domain_structure(bundle)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE PUL_parameter_module
USE filter_module
USE maths

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(INOUT)    :: bundle
  
! local variables

! global cable, conductor and domain numbers and loop variables
  integer :: cable,cable_t
  integer :: tot_n_cables  
  integer :: tot_n_cables_without_ground_plane
  
  integer :: tot_n_conductors,conductor,conductor_count,terminal_conductor
  integer :: tot_n_domains,tot_n_viable_domains,domain,domain_count
  integer :: tot_n_internal_domains

! local conductor and domain numbers and loop variables

  integer :: local_n_conductors,local_conductor,n_cable_conductors
  integer :: local_cable_conductor,local_domain_conductor
  integer :: local_n_domains,local_domain
  integer :: first_external_domain
  integer :: reference_conductor,global_reference_conductor
  integer :: first_external_conductor
  
! variables for looking at overshields and determining which domain cables are in
  integer :: tot_n_overshields,n_overshield_domains,overshield,reference_overshield
  real(dp)::radius_reference_overshield
  
  integer,allocatable  :: overshield_shape(:)
  real(dp),allocatable :: overshield_x(:)
  real(dp),allocatable :: overshield_y(:)
  real(dp),allocatable :: overshield_r(:)
  real(dp),allocatable :: overshield_w(:)
  real(dp),allocatable :: overshield_w2(:)
  real(dp),allocatable :: overshield_h(:)
  integer,allocatable  :: overshield_domain(:)
  integer,allocatable  :: overshield_reference_terminal_conductor(:)
  integer,allocatable  :: overshield_to_cable_number(:)

  integer,allocatable  :: overshield_n_conductors(:)
  integer,allocatable  :: overshield_terminal_conductor(:,:) 
  
  real(dp) :: cable_x,cable_y,dist_cable_to_overshield_centre,cable_r
  logical :: is_overshield_domain

! referencing arrays
  integer,allocatable :: terminal_conductor_to_cable(:)  
  integer,allocatable :: terminal_conductor_to_cable_local_domain(:)  
  integer,allocatable :: terminal_conductor_to_global_domain(:)  
  integer,allocatable :: terminal_conductor_to_global_domain_conductor(:)  
  integer,allocatable :: terminal_conductor_to_reference_global_domain_conductor(:)  
  integer,allocatable :: terminal_conductor_to_reference_terminal_conductor(:)  
  integer,allocatable :: global_domain_conductor_to_terminal_conductor(:)  
  
  logical,allocatable :: terminal_conductor_is_reference_conductor(:)  
  logical,allocatable :: domain_is_TP_differential_mode(:)  

  integer,allocatable :: cable_reference_domain(:) 
  integer,allocatable :: cable_reference_conductor(:) 
  
  integer,allocatable :: external_terminal_conductor(:)
  logical :: is_external_domain
  integer :: external_conductor_count
  
  integer,allocatable :: global_domain_reference_conductor(:)

! structure used by the Per-Unit-Length parameter calculation for overshield and external domains
  type(PUL_type)    :: PUL

! loop variables for matrix operations
  integer :: dim
  integer :: row,col
  integer :: row_l,col_l
  integer :: row_g,col_g,row_g_ref
  
  character(LEN=3) :: conductor_string  ! used to generate the conductor labels in the bundle
  
  real(dp),allocatable :: M1(:,:),M1I(:,:),M2(:,:)
  integer :: Mdim
  integer :: ierr
  real(dp) :: Mdiff

! START

  if(verbose) write(*,*)'CALLED: create_global_domain_structure'
  
! Copy some data to local variables for clarity of notation

  tot_n_cables=bundle%n_cables
  
  if (bundle%ground_plane_present) then
    tot_n_cables_without_ground_plane=tot_n_cables-1
  else
    tot_n_cables_without_ground_plane=tot_n_cables
  end if
  
  if(verbose) write(*,*)' Total number of cables=',tot_n_cables
  if(verbose) write(*,*)' Total number of cables not including ground plane=',tot_n_cables_without_ground_plane
  
! 1a. Count the total number of conductors in the bundle

  tot_n_conductors=0
  do cable=1,tot_n_cables  
    tot_n_conductors = tot_n_conductors + bundle%cable(cable)%tot_n_conductors                   
  end do
  
  if(verbose) write(*,*)' Total number of conductors=',tot_n_conductors
  
  bundle%tot_n_conductors=tot_n_conductors
  bundle%system_dimension=tot_n_conductors-1
  
! 1b Allocate and reset the local data structures 

  ALLOCATE( terminal_conductor_to_cable(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_to_cable_local_domain(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_to_global_domain(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_to_global_domain_conductor(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_to_reference_global_domain_conductor(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_to_reference_terminal_conductor(1:tot_n_conductors) )
  ALLOCATE( terminal_conductor_is_reference_conductor(1:tot_n_conductors) )
  ALLOCATE( global_domain_conductor_to_terminal_conductor(1:tot_n_conductors) )
  terminal_conductor_to_cable(1:tot_n_conductors)=0
  terminal_conductor_to_cable_local_domain(1:tot_n_conductors)=0
  terminal_conductor_to_global_domain(1:tot_n_conductors)=0
  terminal_conductor_to_global_domain_conductor(1:tot_n_conductors)=0
  terminal_conductor_to_reference_global_domain_conductor(1:tot_n_conductors)=0
  terminal_conductor_to_reference_terminal_conductor(1:tot_n_conductors)=0
  terminal_conductor_is_reference_conductor(1:tot_n_conductors)=.FALSE.
  global_domain_conductor_to_terminal_conductor(1:tot_n_conductors)=0
   
! numbering information required for the transfer impedance calculation
  ALLOCATE( bundle%terminal_conductor_is_shield_flag(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_inner_domain(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_outer_domain(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_global_domain_conductor(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_local_domain_conductor(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_reference_terminal_conductor(1:bundle%tot_n_conductors) )
  bundle%terminal_conductor_is_shield_flag(1:bundle%tot_n_conductors)=.FALSE.
  bundle%terminal_conductor_to_inner_domain(1:bundle%tot_n_conductors)=0
  bundle%terminal_conductor_to_outer_domain(1:bundle%tot_n_conductors)=0
  bundle%terminal_conductor_to_global_domain_conductor(1:bundle%tot_n_conductors)=0
  bundle%terminal_conductor_to_local_domain_conductor(1:bundle%tot_n_conductors)=0
  bundle%terminal_conductor_to_reference_terminal_conductor(1:bundle%tot_n_conductors)=0
  
! 2. Loop over the cables filling the referencing arrays with local information
! Also count the total number of internal domains in the bundle
  conductor=0
  tot_n_domains=0
  
  do cable=1,tot_n_cables
  
    local_n_domains=bundle%cable(cable)%tot_n_domains
    tot_n_domains=tot_n_domains+bundle%cable(cable)%n_internal_domains
    
    do local_domain=1,local_n_domains
    
      local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included
      
      do local_conductor=1,local_n_conductors
    
        conductor=conductor+1
        terminal_conductor_to_cable(conductor)=cable
        terminal_conductor_to_cable_local_domain(conductor)=local_domain
        
      end do
      
    end do ! next conductor in this cable
    
  end do
  
  tot_n_internal_domains=tot_n_domains
  first_external_domain=tot_n_domains+1
  
  if(verbose) then
    write(*,*)'Total number of internal domains=',tot_n_domains
    write(*,*)'First external domain=',first_external_domain
  end if
  
! check counting...
  if (conductor.NE.tot_n_conductors) then
    run_status='ERROR in create_global_domain_structure: last conductor.NE.tot_n_conductors'
    CALL write_program_status()
    STOP 1
  end if

! 3. For each cable, work out which domain it is in and hence its reference conductor number.

! 3a. first count the overshields

  tot_n_overshields=0
  do cable=1,tot_n_cables
    if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_overshield) then 
      tot_n_overshields=tot_n_overshields+1
    end if
  end do
  
  if(verbose) write(*,*)'Total number of over-shields=',tot_n_overshields  

! 3b. Get the position of the overshields
  if (tot_n_overshields.GT.0) then
  
    ALLOCATE( overshield_shape(1:tot_n_overshields)  )
    ALLOCATE( overshield_x(1:tot_n_overshields)  )
    ALLOCATE( overshield_y(1:tot_n_overshields)  )
    ALLOCATE( overshield_r(1:tot_n_overshields)  )
    ALLOCATE( overshield_w(1:tot_n_overshields)  )
    ALLOCATE( overshield_w2(1:tot_n_overshields)  )
    ALLOCATE( overshield_h(1:tot_n_overshields)  )
    ALLOCATE( overshield_reference_terminal_conductor(1:tot_n_overshields) )
    ALLOCATE( overshield_to_cable_number(1:tot_n_overshields) )
    
    overshield=0
    conductor_count=0
    do cable=1,tot_n_cables
      if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_overshield) then 
        overshield=overshield+1
        overshield_to_cable_number(overshield)=cable
        overshield_shape(overshield)=circle
        overshield_x(overshield)=bundle%cable_x_offset(cable)
        overshield_y(overshield)=bundle%cable_y_offset(cable)
        overshield_r(overshield)=bundle%cable(cable)%parameters(1)  ! overshield radius
        overshield_w(overshield)=0d0
        overshield_w2(overshield)=0d0
        overshield_h(overshield)=0d0
        overshield_reference_terminal_conductor(overshield)=conductor_count+1
        
      else

! This is not an overshield

        if ( (bundle%cable(cable)%cable_type.NE.cable_geometry_type_ground_plane) .AND. (overshield.NE.0) ) then

! check that oveshields are the last conductors defined in the bundle_spec file
          run_status='ERROR, Overshields should be the last cables defined in the bundle specification'
          CALL write_program_status()
          STOP 1
        
        end if      
      
      end if
      conductor_count=conductor_count+bundle%cable(cable)%tot_n_conductors
    end do
    
    ALLOCATE( overshield_domain(1:tot_n_overshields)  )
    overshield_domain(1:tot_n_overshields)=0
    ALLOCATE( overshield_n_conductors(1:tot_n_overshields)  )
    overshield_n_conductors(1:tot_n_overshields)=0
    ALLOCATE( overshield_terminal_conductor(1:tot_n_overshields,1:tot_n_conductors)  ) ! over the top memory allocation here but simpler this way
    overshield_terminal_conductor(1:tot_n_overshields,1:tot_n_conductors)=0
    
  end if  ! n_overshields.GT.0

! allocate array for domain reference conductors assuming worst case that all overshields and external domain are viable  
  ALLOCATE( global_domain_reference_conductor(1:tot_n_domains+tot_n_overshields+1) )
  global_domain_reference_conductor(1:tot_n_domains+tot_n_overshields+1)=0
  
! 3c. Loop over cables and see whether it sits within an overshield or not  
  
  ALLOCATE( cable_reference_domain(1:bundle%n_cables)  )
  cable_reference_domain(1:bundle%n_cables)=0
  ALLOCATE( cable_reference_conductor(1:bundle%n_cables)  )
  cable_reference_conductor(1:bundle%n_cables)=0
  n_overshield_domains=0     ! counter for overshield domains
  
  write(*,*)'**********************************************************************'
  
  do cable=1,tot_n_cables
  
    write(*,*)
    write(*,*)'Checking cable number',cable
  
    if (bundle%cable(cable)%cable_type.NE.cable_geometry_type_overshield) then 
    
! IF THIS CABLE IS NOT AN OVERSHIELD THEN SEE WHETHER IT IS WITHIN AN OVERSHIELD
  
      write(*,*)'THIS CABLE IS NOT AN OVERSHIELD'

      cable_x=bundle%cable_x_offset(cable)
      cable_y=bundle%cable_y_offset(cable)

! find the smallest overshield enclosing this cable      

      reference_overshield=0
      radius_reference_overshield=1d30
      
      do overshield=1,tot_n_overshields
            
        dist_cable_to_overshield_centre=sqrt( (cable_x-overshield_x(overshield))**2+  &
                                                (cable_y-overshield_y(overshield))**2  )
                                              
        if ( dist_cable_to_overshield_centre.LT.overshield_r(overshield) ) then
! we conclude that this cable is within the overshield

          write(*,*)'This cable is inside overshield number',overshield

          if (radius_reference_overshield.GT.overshield_r(overshield)) then
! this overshield is smaller than the current reference overshield so this becomes the new reference overshield
            
            reference_overshield=overshield
            radius_reference_overshield=overshield_r(overshield)
            
            write(*,*)'This cable reference overshield has been updated to',reference_overshield

          end if  ! is this the smallest overshield?
          
        end if  ! is cable within this overshield?
                 
      end do
      
      if (reference_overshield.NE.0) then
! this cable has been identified as being within an overshield

        if (overshield_domain(reference_overshield).EQ.0) then 
          
! this overshield domain is unallocated so create a new domain
  
          tot_n_domains=tot_n_domains+1
          n_overshield_domains=n_overshield_domains+1
          if(verbose) write(*,*)'Creating overshield domain=',tot_n_domains
          overshield_domain(reference_overshield)=tot_n_domains
            
          global_domain_reference_conductor(tot_n_domains)=overshield_reference_terminal_conductor(reference_overshield)
            
          write(*,*)'This overshield domain is unallocated so create a new domain, domain=',tot_n_domains
                        
        end if ! this overshield domain is unallocated
          
        cable_reference_domain(cable)=overshield_domain(reference_overshield)
          
        write(*,*)'Setting the cable reference domain to',overshield_domain(reference_overshield)
          
      end if   ! this cable is within the overshield
             
      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count
        
    else  ! THIS CABLE IS AN OVERSHIELD SO WE NEED TO CHECK WHETHER IT IS INSIDE ANY OTHER OVERSHIELDS   
   
      write(*,*)'THIS CABLE IS AN OVERSHIELD'
       
      cable_x=bundle%cable_x_offset(cable)
      cable_y=bundle%cable_y_offset(cable)
      cable_r=bundle%cable(cable)%parameters(1)  ! overshield radius
      
! find the smallest overshield enclosing this overshield if there is one   

      reference_overshield=0
      radius_reference_overshield=1d30
      
      do overshield=1,tot_n_overshields
      
        cable_t=overshield_to_cable_number(overshield)
        if (cable_t.NE.cable) then   ! don't check an overshield against itself
      
          dist_cable_to_overshield_centre=sqrt( (cable_x-overshield_x(overshield))**2+  &
                                                (cable_y-overshield_y(overshield))**2  )
                                              
          if ( dist_cable_to_overshield_centre.LT.overshield_r(overshield) ) then
! we conclude that this cable is within the overshield

            write(*,*)'This overshield centre is inside overshield number',overshield

!! Check that can be removed when we are happy with nexted overshields code            
!            run_status='ERROR, Nested overshields cannot be used in SACAMOS at the moment'
!            CALL write_program_status()
!            STOP 1

            if (radius_reference_overshield.GT.overshield_r(overshield)) then
! this overshield is smaller than the current reference overshield so this becomes the new reference overshield
            
              if (overshield_r(overshield).GT.cable_r) then
! the cable is smaller than the overshield and hence is inside it

! Check that the lowest number overshield is the smaller. Overshields should be numbered inside to out

                if (cable_t.LT.cable) then
              
                  run_status='ERROR, Nested overshields should be numbered in the order inside to out'
                  CALL write_program_status()
                  STOP 1
              
                end if

                write(*,*)'This overshield is smaller than the reference overshield and is therefore inside it'
           
                reference_overshield=overshield
                radius_reference_overshield=overshield_r(overshield)
            
                write(*,*)'This cable reference overshield has been updated to',reference_overshield
              
              end if  ! This overshield is smaller than the reference overshield and is therefore inside it
              
            end if  ! is this the smallest overshield?
          
          end if  ! is cable centre within this overshield?
          
        end if  ! don't check an overshield against itself
       
      end do
      
      if (reference_overshield.NE.0) then
! this cable has been identified as being within an overshield

        if (overshield_domain(reference_overshield).EQ.0) then 
          
! this overshield domain is unallocated so create a new domain
  
          tot_n_domains=tot_n_domains+1
          n_overshield_domains=n_overshield_domains+1
          if(verbose) write(*,*)'Creating overshield domain=',tot_n_domains
          overshield_domain(reference_overshield)=tot_n_domains
            
          global_domain_reference_conductor(tot_n_domains)=overshield_reference_terminal_conductor(reference_overshield)
            
          write(*,*)'This overshield domain is unallocated so create a new domain, domain=',tot_n_domains
                        
        end if ! this overshield domain is unallocated
          
        cable_reference_domain(cable)=overshield_domain(reference_overshield)
          
        write(*,*)'Setting the cable reference domain to',overshield_domain(reference_overshield)
          
      end if   ! this cable is within the overshield
             
      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count
     
       
    end if ! not an overshield
    
  end do ! next cable

  if (verbose) write(*,*) 'n_overshield_domains=',n_overshield_domains
  
  write(*,*)'**********************************************************************'

! 3c part 2. loop over the overshield domains and count the number of external conductors in each one  

  do overshield=1,tot_n_overshields
 
    terminal_conductor=0   ! count the terminal conductor number
    conductor_count=0      ! count the number of condcutors referenced to this overshield
    
    do cable=1,tot_n_cables_without_ground_plane
      
      if (cable_reference_domain(cable).EQ.overshield_domain(overshield)) then
! this cable contributes its external conductors to the overshield

        first_external_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors-bundle%cable(cable)%n_external_conductors+1
        
        do conductor=1,bundle%cable(cable)%n_external_conductors
          conductor_count=conductor_count+1
          overshield_terminal_conductor(overshield,conductor_count)=first_external_conductor+(conductor-1)
        end do
        
      end if
      
      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count

    end do ! next cable  
    
    overshield_n_conductors(overshield)=conductor_count
    
  end do ! next overshield
  
  if (verbose) then
    write(*,*)
    write(*,*)'Overshield conductor information'
    do overshield=1,tot_n_overshields
      write(*,*)'Overshield:',overshield,' Number of condcutors=',overshield_n_conductors(overshield)
      write(*,*)' Overshield_conductor    terminal_conductor '
      do conductor=1,overshield_n_conductors(overshield)
        write(*,*)conductor,overshield_terminal_conductor(overshield,conductor)
      end do
    end do ! next overshield
    write(*,*)
  end if
  
! 3d. We now have unallocated reference domains, these are not in overshields so must be in the external domain

! 4. Count the number of conductors in the external domain to check whether it is viable i.e. at least 2 conductors

  bundle%tot_n_external_conductors=0
  do cable=1,tot_n_cables    
    if (cable_reference_domain(cable).EQ.0) then    
! this cable reference domain must be the external domain
      bundle%tot_n_external_conductors=bundle%tot_n_external_conductors+bundle%cable(cable)%n_external_conductors     
    end if ! reference domain not yet set    
  end do ! next cable
  
  if(verbose) write(*,*)'Number of conductors in external domain=',bundle%tot_n_external_conductors

! 4b.     Create an external domain and allocate unallocated reference domains now...
  ALLOCATE( external_terminal_conductor(1:bundle%tot_n_external_conductors) )
  external_terminal_conductor(1:bundle%tot_n_external_conductors)=0
  
  tot_n_viable_domains=tot_n_domains
  tot_n_domains=tot_n_domains+1
  
  if (bundle%tot_n_external_conductors.GT.1) then
    tot_n_viable_domains=tot_n_viable_domains+1
  end if
  
  if(verbose) write(*,*)'Creating external domain=',tot_n_domains
    
  terminal_conductor=0
  external_conductor_count=0
    
  do cable=1,tot_n_cables
    
    is_external_domain=.FALSE.
    if (cable_reference_domain(cable).EQ.0) then ! this cable reference domain must be the external domain
      cable_reference_domain(cable)=tot_n_domains
      is_external_domain=.TRUE.
    end if
      
    do local_conductor=1,bundle%cable(cable)%tot_n_conductors
    
      terminal_conductor=terminal_conductor+1
        
      if ( (is_external_domain).AND.(bundle%cable(cable)%local_reference_conductor(local_conductor).EQ.0) ) then        
        external_conductor_count=external_conductor_count+1
        external_terminal_conductor(external_conductor_count)=terminal_conductor      
      end if
        
    end do
        
    if (is_external_domain) then ! this cable reference domain must be the external domain
      global_domain_reference_conductor(tot_n_domains)=terminal_conductor
    end if
      
  end do ! next cable
    
! check the conductor count
  if (external_conductor_count.NE.bundle%tot_n_external_conductors) then
    run_status='ERROR in create_global_domain_structure: external_conductor_count.NE.bundle%tot_n_external_conductors'
    CALL write_program_status()
    STOP 1
  end if
  
  if(verbose) then  
    write(*,*)'External domain terminal conductors:'
    do conductor=1,bundle%tot_n_external_conductors
      write(*,*)'conductor=',external_terminal_conductor(conductor)
    end do ! next cable
  end if
  
  bundle%tot_n_domains=tot_n_viable_domains
  
  if(verbose) then  
    write(*,*)'Total number of domains=',tot_n_domains
    write(*,*)'Total number of viable domains=',tot_n_viable_domains
    write(*,*)'Cable reference domain numbers:'
    do cable=1,tot_n_cables  
      write(*,*)'cable=',cable,' reference domain=',cable_reference_domain(cable) 
    end do ! next cable
  end if
  
  ALLOCATE( domain_is_TP_differential_mode(1:tot_n_domains) )
  domain_is_TP_differential_mode(1:tot_n_domains)=.FALSE.

! 5. build the global_domain numbering, ignoring any unused domains

  conductor=0       ! this is the external conductor number
  domain_count=0
  conductor_count=0 ! this is the conductor count for the domain based numbering
  
! 5a. Count conductors from internal domains initially    
  do cable=1,tot_n_cables
  
    local_cable_conductor=0
    local_n_domains=bundle%cable(cable)%tot_n_domains
        
! 17/5/2017 CJS Add a check here to flag twisted pair differential mode domains
! The check is hrd wired based on cable types
 
    if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_twisted_pair) then 
    
      domain_is_TP_differential_mode(domain_count+1)=.TRUE.
      
    else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_shielded_twisted_pair) then 
    
      domain_is_TP_differential_mode(domain_count+1)=.TRUE.

    else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_spacewire) then 
    
      domain_is_TP_differential_mode(domain_count+1)=.TRUE.
      domain_is_TP_differential_mode(domain_count+3)=.TRUE.
      domain_is_TP_differential_mode(domain_count+5)=.TRUE.
      domain_is_TP_differential_mode(domain_count+7)=.TRUE.
    
    end if

! exclude the external domain  
    do local_domain=1,local_n_domains
    
      if (local_domain.NE.local_n_domains) then

! This is an internal domain so add to the domain count
        domain_count=domain_count+1
        domain=domain_count
      
        local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included
      
        do local_conductor=1,local_n_conductors
    
          conductor=conductor+1
          local_cable_conductor=local_cable_conductor+1
          conductor_count=conductor_count+1
          terminal_conductor_to_global_domain(conductor)=domain
          terminal_conductor_to_cable_local_domain(conductor)=local_conductor
          terminal_conductor_to_global_domain_conductor(conductor)=conductor_count                        
          bundle%cable(cable)%global_domain_conductor(local_cable_conductor)=conductor_count
          bundle%cable(cable)%terminal_conductor(local_cable_conductor)=conductor
       
        end do
        
      else
! This is an external domain so leave until later but count the global conductors.

        local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included     
        do local_conductor=1,local_n_conductors    
          conductor=conductor+1        
          local_cable_conductor=local_cable_conductor+1
        end do
      
      end if ! internal or external domain
      
    end do ! next conductor in this cable
    
  end do ! next cable

! 5b. Include the 'external' domain conductors now  
  do domain=first_external_domain,tot_n_domains
  
    conductor=0       ! this is the external conductor number
    do cable=1,tot_n_cables
  
      local_cable_conductor=0
      local_n_domains=bundle%cable(cable)%tot_n_domains

      do local_domain=1,local_n_domains
    
        if (local_domain.NE.local_n_domains) then

! This is an internal domain so this is already included, only count the conductors
      
          local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included
      
          do local_conductor=1,local_n_conductors
    
            conductor=conductor+1
            local_cable_conductor=local_cable_conductor+1
        
          end do
        
        else
! This is an external domain so add the conductors if it is the currrent external domain

          local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included     
          do local_conductor=1,local_n_conductors    
            conductor=conductor+1       
            local_cable_conductor=local_cable_conductor+1
            
            if (cable_reference_domain(cable).EQ.domain) then
            
              conductor_count=conductor_count+1       
              terminal_conductor_to_global_domain(conductor)=domain
              terminal_conductor_to_cable_local_domain(conductor)=local_conductor
              terminal_conductor_to_global_domain_conductor(conductor)=conductor_count                       
              bundle%cable(cable)%global_domain_conductor(local_cable_conductor)=conductor_count
              bundle%cable(cable)%terminal_conductor(local_cable_conductor)=conductor
              
            end if ! this external domain is the current one to be included
             
          end do ! next local conductor
      
        end if ! internal or external domain
        
      end do ! next local domain
           
    end do ! next cable
    
  end do ! next external domain
  
! 5c. Set the  global_domain_conductor_to_terminal_conductor from terminal_conductor_to_global_domain_conductor
  
  if (verbose) then
  
    write(*,*)'conductor,terminal_conductor_to_global_domain_conductor(conductor)'
    
    do conductor=1,tot_n_conductors
      write(*,*)conductor,terminal_conductor_to_global_domain_conductor(conductor)
    end do 
        
  end if 
  
! loop over terminal conductors
  do conductor=1,tot_n_conductors
    global_domain_conductor_to_terminal_conductor(terminal_conductor_to_global_domain_conductor(conductor))=conductor
  end do 
  
  if (verbose) then
  
    write(*,*)'conductor,global_domain_conductor_to_terminal_conductor(conductor)'
    
    do conductor=1,tot_n_conductors
      write(*,*)conductor,global_domain_conductor_to_terminal_conductor(conductor)
    end do 
        
  end if 

! 22/10/2019 CJS Fix bug with numbering here: overshields were numbered using terminal conductor number
!                not global domain conductor as it should be at this stage
! 5d. Renumber the reference conductors for overshield domains. At the moment these are terminal conductors
!     not global domain conductors

  do overshield=1,tot_n_overshields
  
    domain=overshield_domain(overshield)

    reference_conductor=global_domain_reference_conductor(domain)
      
    global_reference_conductor=terminal_conductor_to_global_domain_conductor(reference_conductor)
    global_domain_reference_conductor(domain)=global_reference_conductor

  end do ! next overshield 
   
! 5e. Set the global_domain_reference_conductor information, also flag all the reference conductors in internal domains
  
  conductor=0       ! this is the external conductor number
  
  write(*,*)
  write(*,*)
  write(*,*)'    terminal     cable      local    local_domain   global      cable     local_ref   ref terminal   ref global     '
  write(*,*)'    conductor              domain     conductor     domain     conductor  conductor    conductor    domain conductor'
  
  do cable=1,tot_n_cables
  
    local_n_domains=bundle%cable(cable)%tot_n_domains
    
    local_cable_conductor=0

    do local_domain=1,local_n_domains
          
      local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included
                  
      do local_domain_conductor=1,local_n_conductors
      
        conductor=conductor+1
        local_cable_conductor=local_cable_conductor+1
        
        reference_conductor=bundle%cable(cable)%local_reference_conductor(local_cable_conductor)
        
        if (reference_conductor.NE.0) then
! the reference conductor is defined locally        
          terminal_conductor_to_reference_global_domain_conductor(conductor)=   &
                      bundle%cable(cable)%global_domain_conductor(reference_conductor)

! 27/4/2016 CJS                      
          domain=terminal_conductor_to_global_domain(conductor)
! 16/5/2016 this should be the global not the local number...          global_domain_reference_conductor(domain)=reference_conductor
          global_domain_reference_conductor(domain)=bundle%cable(cable)%global_domain_conductor(reference_conductor)
                               
        else
        
! the reference conductor is defined in terms of the reference domain of the cable         
   
          domain=cable_reference_domain(cable) 
          
          global_reference_conductor=global_domain_reference_conductor(domain)     
          
          reference_conductor=global_domain_conductor_to_terminal_conductor(global_reference_conductor)
          
          terminal_conductor_to_reference_global_domain_conductor(conductor)=  &
                    terminal_conductor_to_global_domain_conductor(reference_conductor)
                    
          reference_conductor=0
                    
        end  if
        
        local_conductor=terminal_conductor_to_reference_global_domain_conductor(conductor)
        terminal_conductor_to_reference_terminal_conductor(conductor)=  &
          global_domain_conductor_to_terminal_conductor(local_conductor)
          
        terminal_conductor_is_reference_conductor(terminal_conductor_to_reference_terminal_conductor(conductor))=.TRUE.
        
        write(*,8000)bundle%cable(cable)%terminal_conductor(local_cable_conductor),&
                     cable,local_domain,local_domain_conductor, &
                     terminal_conductor_to_global_domain(conductor),    &
                     local_cable_conductor,reference_conductor, &
                     terminal_conductor_to_reference_terminal_conductor(conductor), &
                     terminal_conductor_to_reference_global_domain_conductor(conductor)
8000 format(5I11,4I13)
                            
      end do ! next local domain conductor
        
    end do ! next local domain
           
  end do ! next cable
  
   
! 6.  Now that we have counted the number of domains we can allocate the 
!     domain based number of conductors and L and C matrices plus
!     the global MI and MV matrices 

! Note: this allocation includes non-viable domains
  ALLOCATE( bundle%n_conductors(1:tot_n_domains) ) 
  bundle%n_conductors(1:tot_n_domains)=0
  
  ALLOCATE( bundle%L(1:tot_n_domains) ) 
  ALLOCATE( bundle%C(1:tot_n_domains) )
  ALLOCATE( bundle%Z(1:tot_n_domains) ) 
  ALLOCATE( bundle%Y(1:tot_n_domains) )

! 6a. Count the number of conductors in each domain
  if (verbose) write(*,*)'terminal_conductor   domain   local_domain_conductor'
  do conductor=1,tot_n_conductors 
    domain=terminal_conductor_to_global_domain(conductor)
    bundle%n_conductors(domain)=bundle%n_conductors(domain)+1

! 4/10/2016 Save the local domain coonductor numbering
    bundle%terminal_conductor_to_local_domain_conductor(conductor)=bundle%n_conductors(domain)
    
    if (verbose) write(*,*)conductor,domain,bundle%terminal_conductor_to_local_domain_conductor(conductor)

  end do ! next conductor
  
! Note there may be unviable overshield domains which need to be excluded in the following loop...  
  
! internal  domains and overshield domains don't include the reference conductor so add this here
  do domain=1,tot_n_internal_domains+n_overshield_domains
    bundle%n_conductors(domain)=bundle%n_conductors(domain)+1
  end do
  
  if (verbose) then
    write(*,*)
    write(*,*)'Number of conductors in each domain'
    write(*,*)'       domain     n_conductors'
    do domain=1,tot_n_domains
      write(*,*)domain,bundle%n_conductors(domain)
    end do
  end if
  
! allocate the global voltage and current domain transformation matrices
! strictly the dimension should be bundle%tot_n_conductors-1
! this would require many checks later on to prevent array bounds problems with no ground plane specified
! so we allocate the extra space then ignore the last row and col when doing the domain decomposition

  dim=bundle%tot_n_conductors
  
  if(verbose) then
    write(*,*)'ALLOCATING global_MI and global_MV, dimension:',dim
  end if

  bundle%global_MI%dim=dim
  ALLOCATE( bundle%global_MI%mat(1:dim,1:dim) ) 
  bundle%global_MI%mat(1:dim,1:dim)=0d0

  bundle%global_MV%dim=dim
  ALLOCATE( bundle%global_MV%mat(1:dim,1:dim) ) 
  bundle%global_MV%mat(1:dim,1:dim)=0d0

! Allocate the conductor based impedance (loss) model and reset all the parameters
  ALLOCATE( bundle%conductor_impedance(1:dim) )
  do conductor=1,bundle%tot_n_conductors
    bundle%conductor_impedance(conductor)%radius=0d0
    bundle%conductor_impedance(conductor)%width=0d0
    bundle%conductor_impedance(conductor)%height=0d0
    bundle%conductor_impedance(conductor)%conductivity=0d0
    bundle%conductor_impedance(conductor)%thickness=0d0
    bundle%conductor_impedance(conductor)%Resistance_multiplication_factor=1d0
  end do
  
! allocate the global inductance and capacitance matrices - these are used in the analytic solution
! which is used for validation

  dim=bundle%tot_n_conductors-1 
  
  if(verbose) then
    write(*,*)'ALLOCATING global_L and global_C, dimension:',dim
  end if

  bundle%global_L%dim=dim
  ALLOCATE( bundle%global_L%mat(1:dim,1:dim) ) 
  bundle%global_L%mat(1:dim,1:dim)=0d0

  bundle%global_C%dim=dim
  ALLOCATE( bundle%global_C%mat(1:dim,1:dim) ) 
  bundle%global_C%mat(1:dim,1:dim)=0d0

  if(verbose) then
    write(*,*)'ALLOCATING global_Z and global_Y, dimension:',dim
  end if
  bundle%global_Z%dim=dim
  ALLOCATE( bundle%global_Z%sfilter_mat(1:dim,1:dim) ) 
  do row=1,dim
    do col=1,dim
      bundle%global_Z%sfilter_mat(row,col)=0d0         ! set all filters to zero filter for now
    end do
  end do
  
  bundle%global_Y%dim=dim
  ALLOCATE( bundle%global_Y%sfilter_mat(1:dim,1:dim) ) 
  do row=1,dim
    do col=1,dim
      bundle%global_Y%sfilter_mat(row,col)=0d0        ! set all filters to zero filter for now
    end do
  end do
  
! Allocate the domain based terminal conductor list  CJS 27/4/2016
  ALLOCATE( bundle%terminal_conductor_list(1:bundle%tot_n_domains))
  
! 6b. Copy the cable based internal domain L and C matrices to bundle structure
!     and construct the global MI and MV matrices from the cable based structures 
 
  terminal_conductor=0       ! this is the external conductor number
  domain_count=0
  
  if (verbose) then
    write(*,*)''
    write(*,*)'SET MI, MV matrices'
    write(*,*)''
  end if
  
  do cable=1,tot_n_cables_without_ground_plane
  
    local_n_domains=bundle%cable(cable)%tot_n_domains
    
    local_cable_conductor=0
    n_cable_conductors=bundle%cable(cable)%tot_n_conductors
        
! copy the cable based MI and MV matrices to the global structure. The cable based
! MI and MV matrices are on the basis of terminal_conductor numbering
    if (verbose) then
      write(*,*)''
      write(*,*)'Copy MI, MV, cable=',cable,' of', tot_n_cables_without_ground_plane
      write(*,*)'Matrix dimension',n_cable_conductors+1
      write(*,*)'intial terminal_conductor count=',terminal_conductor
    end if
    
    do row_l=1,n_cable_conductors+1 ! add 1 as the reference conductor must be included
                
      do col_l=1,n_cable_conductors+1 ! add 1 as the reference conductor must be included

! In the global MI, MV structure rows correspond to the global domain conductor number, cols correspond to the terminal conductor number 
        
        if (row_l.NE.n_cable_conductors+1) then
          row_g=terminal_conductor_to_global_domain_conductor(terminal_conductor+row_l)
        else
          row_g=terminal_conductor_to_reference_global_domain_conductor(terminal_conductor+row_l-1) ! last conductor for this cable and the reference    
        end if
     
        if (col_l.NE.n_cable_conductors+1) then
          col_g=terminal_conductor+col_l
        else
          col_g=terminal_conductor_to_reference_terminal_conductor(terminal_conductor+col_l-1) ! last conductor for this cable and the reference
        end if
        
        bundle%global_MI%mat(row_g,col_g)=bundle%cable(cable)%MI%mat(row_l,col_l)
        bundle%global_MV%mat(row_g,col_g)=bundle%cable(cable)%MV%mat(row_l,col_l)
        
!      if (verbose) write(*,*)'Copy MV element',row_g,col_g,bundle%global_MV%mat(row_g,col_g)

          if (verbose) then
          write(*,'(A,A,I3,A,I3,A,I3,A,I3,F5.2,F5.2)')'Copy MI element',' row_g=',row_g,' col_g=',col_g,      &
                  ' row_l=',row_l,' col_l=',col_l,bundle%cable(cable)%MI%mat(row_l,col_l),bundle%cable(cable)%MV%mat(row_l,col_l)
        end if
        
      end do ! next col
      
    end do ! next row
    
    if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_overshield) then 
! This is an overshield so invoke the special process for overshields

! Get the global domain number of the shield conductor
      row_g=terminal_conductor_to_global_domain_conductor(terminal_conductor+1)
      
! Get the global domain number of the shield reference conductor
      row_g_ref=terminal_conductor_to_reference_global_domain_conductor(terminal_conductor+1)
      
! Add all the non-zero elements of the MI matrix shield row to the shield reference row

      write(*,*)'Overshield cable',cable,' Copying row',row_g,' to row',row_g_ref
                      
      do col_g=1,bundle%global_MI%dim

! In the global MI, MV structure rows correspond to the global domain conductor number, cols correspond to the terminal conductor number 
                
        if (bundle%global_MI%mat(row_g,col_g).NE.0d0) then
        
          bundle%global_MI%mat(row_g_ref,col_g)=bundle%global_MI%mat(row_g,col_g)
        
        end if
                
      end do ! next col
      
    end if ! overshield
    
    do local_domain=1,local_n_domains  ! exclude the external domain for now
    
      if (local_domain.NE.local_n_domains) then

! This is an internal domain so add to the domain count
        
        domain_count=domain_count+1 
        domain=domain_count
        
! 6c copy the domain based L and C matrices to the domain based L and C in the bundle structure                       #

        if (verbose) write(*,*)'Copy the domain based L and C matrices to the domain based L and C in the bundle structure'

        if (verbose) write(*,*)'domain=',domain   
 
        dim=bundle%cable(cable)%L_domain(local_domain)%dim
        ALLOCATE( bundle%L(domain)%mat(dim,dim) ) 
        bundle%L(domain)%dim=dim
        if (verbose) write(*,*)'Allocating bundle%L(domain)        ,dim=',dim
        
        dim=bundle%cable(cable)%C_domain(local_domain)%dim
        ALLOCATE( bundle%C(domain)%mat(dim,dim) ) 
        bundle%C(domain)%dim=dim
        if (verbose) write(*,*)'Allocating bundle%C(domain)        ,dim=',dim
 
        dim=bundle%cable(cable)%Z_domain(local_domain)%dim
        ALLOCATE( bundle%Z(domain)%sfilter_mat(dim,dim) ) 
        bundle%Z(domain)%dim=dim
        if (verbose) write(*,*)'Allocating bundle%Z(domain)        ,dim=',dim
        
        dim=bundle%cable(cable)%Y_domain(local_domain)%dim
        ALLOCATE( bundle%Y(domain)%sfilter_mat(dim,dim) ) 
        bundle%Y(domain)%dim=dim
        if (verbose) write(*,*)'Allocating bundle%Y(domain)        ,dim=',dim
      
        local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)
      
        do row_l=1,local_n_conductors-1 ! subtract 1 as the reference conductor is included
                    
          do col_l=1,local_n_conductors-1 ! subtract 1 as the reference conductor is included
          
            bundle%L(domain)%mat(row_l,col_l)=bundle%cable(cable)%L_domain(local_domain)%mat(row_l,col_l)      
            bundle%C(domain)%mat(row_l,col_l)=bundle%cable(cable)%C_domain(local_domain)%mat(row_l,col_l)
            
            bundle%Z(domain)%sfilter_mat(row_l,col_l)=bundle%cable(cable)%Z_domain(local_domain)%sfilter_mat(row_l,col_l)      
            bundle%Y(domain)%sfilter_mat(row_l,col_l)=bundle%cable(cable)%Y_domain(local_domain)%sfilter_mat(row_l,col_l)

! copy the domain based L and C matrices to the global L and C in the bundle structure                        

! in the global L,C structure rows and columns correspond to the global_domain conductor number
            row_g=terminal_conductor_to_global_domain_conductor(terminal_conductor+row_l)
            col_g=terminal_conductor_to_global_domain_conductor(terminal_conductor+col_l)

            bundle%global_L%mat(row_g,col_g)=bundle%L(domain)%mat(row_l,col_l)
            bundle%global_C%mat(row_g,col_g)=bundle%C(domain)%mat(row_l,col_l)

            bundle%global_Z%sfilter_mat(row_g,col_g)=bundle%Z(domain)%sfilter_mat(row_l,col_l)
            bundle%global_Y%sfilter_mat(row_g,col_g)=bundle%Y(domain)%sfilter_mat(row_l,col_l)
         
          end do ! next column of matrix
          
        end do ! next row of matrix

! Added CJS 27/4/2016
        bundle%terminal_conductor_list(domain)%n_elements=local_n_conductors
        ALLOCATE( bundle%terminal_conductor_list(domain)%element(1:local_n_conductors) )
        
        do conductor=1,local_n_conductors-1   ! do reference conductor separately    
        
          bundle%terminal_conductor_list(domain)%element(conductor)=terminal_conductor+conductor    
            
        end do ! next conductor in this domain
        
! reference conductor        
        reference_conductor=global_domain_reference_conductor(domain)       ! note domain based numbering...
        write(*,*)'Domain=',domain
        write(*,*)'Global domain reference conductor=',reference_conductor
        bundle%terminal_conductor_list(domain)%element(local_n_conductors)=  &
                                               global_domain_conductor_to_terminal_conductor(reference_conductor)
          
        terminal_conductor=terminal_conductor+local_n_conductors-1  ! subtract 1 as the count includes the reference conductor

      else
! This is an external domain so leave until later but count the global conductors.

        local_n_conductors=bundle%cable(cable)%local_domain_n_conductors(local_domain)-1 ! subtract 1 as the reference conductor is included     
        do local_conductor=1,local_n_conductors    
          terminal_conductor=terminal_conductor+1        
        end do
      
      end if ! internal or external domain

    end do ! next local domain
    
  end do ! next cable
      
!  if (verbose) then

    write(*,*)''
    write(*,*)'CHECK CONSISTENCY OF MI, MV matrices, (MV)T=(MI^-1)'
    write(*,*)''
    
    Mdim=bundle%global_MI%dim
    
! This is the matrix dimensions however if there is no ground plane then the last conductor is
! the reference and hence the matrix dimension is reduced by one

    if (.NOT.bundle%ground_plane_present) then
      Mdim=Mdim-1
    end if
    
    ALLOCATE( M1(1:Mdim,1:Mdim) )
    ALLOCATE( M2(1:Mdim,1:Mdim) )
    ALLOCATE( M1I(1:Mdim,1:Mdim) )    
    
    M1(1:Mdim,1:Mdim)=bundle%global_MI%mat(1:Mdim,1:Mdim)

    do row_l=1,Mdim
      do col_l=1,Mdim
        M2(row_l,col_l)=bundle%global_MV%mat(col_l,row_l)
      end do
    end do
    
    CALL dinvert_Gauss_Jordan(M1,Mdim,M1I,Mdim,ierr) 
    
    Mdiff=0d0
    do row_l=1,Mdim
      do col_l=1,Mdim
        Mdiff=Mdiff+abs( M1I(row_l,col_l)-M2(row_l,col_l) )
      end do
    end do
    
    write(*,*)'MI='
    CALL dwrite_matrix(M1,Mdim,Mdim,Mdim,0)
    write(*,*)''
    write(*,*)'(MI)^-1^='
    CALL dwrite_matrix(M1I,Mdim,Mdim,Mdim,0)
    write(*,*)''
    write(*,*)'(MV)T='
    CALL dwrite_matrix(M2,Mdim,Mdim,Mdim,0)
    
    write(*,*)''
    write(*,*)'Matrix difference=',Mdiff
    write(*,*)''
        
    DEALLOCATE( M1 )
    DEALLOCATE( M2 )
    DEALLOCATE( M1I )
    
    if (Mdiff.GT.0.001) then
      write(*,*)'ERROR in create_global_domain_structure'
      write(*,*)'(MV)T .NE. (MI)^-1'
      STOP
    end if  
    
!  end if

! 7. Calculate the inductance and capacitance matrices for the external domains (domains within overshields and the external domain)

  if(verbose) write(*,*)'first external domain:',first_external_domain
  if(verbose) write(*,*)'total number of domains:',tot_n_domains
  
  overshield=0    ! counter for overshield domains

  do domain=first_external_domain,tot_n_viable_domains   ! exclude the external domain if there is only one conductor
  
! allocate memory for the PUL parameter solver interface

    if(verbose) write(*,*)'Domain:',domain
    if(verbose) write(*,*)'Allocating PUL data structure for ',bundle%n_conductors(domain),' conductors'
    
    CALL allocate_and_reset_PUL_data(PUL,bundle%n_conductors(domain))
     
! copy external conductor information to the PUL structure

    conductor=0        ! this provides a count of the external domain conductor numbering for PUL parameter calculation
  
! loop over cables  
    do cable=1,bundle%n_cables
    
      if (cable_reference_domain(cable).EQ.domain) then
      
! The reference domain of this cable is the domain we are currently working on so add the external 
! conductor information for this cable to the PUL structure

        do local_conductor=1,bundle%cable(cable)%n_external_conductors
      
          if (bundle%cable(cable)%cable_type.NE.cable_geometry_type_ground_plane) then
             
            conductor=conductor+1  ! add another conductor to the external domain
         
            if(verbose) write(*,*)'Adding conductor number',local_conductor,' from cable',cable,'. Domain conductor',conductor
      
            PUL%shape(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_type
            PUL%x(conductor)=bundle%cable_x_offset(cable)    
            PUL%y(conductor)=bundle%cable_y_offset(cable)   
            PUL%rtheta(conductor)=bundle%cable_angle(cable)
            
            PUL%ox(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox
            PUL%oy(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_oy
            PUL%r(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_radius
            PUL%rw(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width
            PUL%rw2(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width2
            PUL%rh(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_height
            PUL%rd(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_radius
            PUL%rdw(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_width
            PUL%rdh(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_height
            PUL%rdox(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_ox
            PUL%rdoy(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_oy
            PUL%epsr(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_epsr 
               
          end if ! not a ground plane    
              
        end do ! next cable external conductor
              
      end if  ! cable contributes to this domain
          
    end do ! next cable
   
! Calculate the per-unit-length parameters (inductance and capacitance matrices)
    
! if this is an overshield domain then
    if (verbose) write(*,*)'domain=',domain,' last domain(check)=',first_external_domain+tot_n_overshields
    if(domain.LT.(first_external_domain+n_overshield_domains)) then
    
      overshield=overshield+1

      if (verbose) write(*,*)'PUL parameter calculation for oversheild domain'
! no ground plane
      PUL%ground_plane_present=.FALSE.
      PUL%ground_plane_angle  =0d0
      PUL%ground_plane_offset =0d0
     
! add overshield information   
      if (verbose) write(*,*)'Add overshield information'
      is_overshield_domain=.TRUE.
      PUL%overshield_present=.TRUE.
      overshield=domain-first_external_domain+1
      PUL%overshield_shape = overshield_shape(overshield)
      PUL%overshield_x = overshield_x(overshield)
      PUL%overshield_y = overshield_y(overshield)
      PUL%overshield_r = overshield_r(overshield)
      PUL%overshield_w = overshield_w(overshield)
      PUL%overshield_w2 = overshield_w2(overshield)
      PUL%overshield_h = overshield_h(overshield)
      
      PUL%epsr_background = 1d0 ! background permittivity =1.0 within an overshield i.e. cables are in air
      
      if (use_Laplace) then
      
         CALL PUL_LC_Laplace(PUL,bundle%bundle_name,bundle%Y_fit_model_order,bundle%Y_fit_freq_spec,domain)
      
      else
      
        if (verbose) write(*,*)'CALL PUL_LC_calc_overshield_wide_separation_approximation'
        CALL PUL_LC_calc_overshield_wide_separation_approximation(PUL)
        
      end if ! use_Laplace or not...

    else
! this is an external domain    
    
! copy ground plane information
      PUL%ground_plane_present=bundle%ground_plane_present
      PUL%ground_plane_angle  =bundle%ground_plane_angle
      PUL%ground_plane_offset =bundle%ground_plane_offset
     
! No overshield information   
      is_overshield_domain=.FALSE.
      PUL%overshield_present=.FALSE.
      PUL%overshield_shape=0
      PUL%overshield_r= 0d0
      PUL%overshield_w= 0d0
      PUL%overshield_w2= 0d0
      PUL%overshield_h= 0d0
      
      PUL%epsr_background = 1d0 ! background permittivity =1.0 i.e. cables are in air

      if (use_Laplace) then
      
         CALL PUL_LC_Laplace(PUL,bundle%bundle_name,bundle%Y_fit_model_order,bundle%Y_fit_freq_spec,domain)
      
      else

        CALL PUL_LC_calc_wide_separation_approximation(PUL)

      end if ! use_Laplace or not...

    end if  ! last external domain (i.e. not an overshield domain)

! We now have the L and C matrices for this external domain        
! Copy the domain based L and C matrices to the domain based L and C in the bundle structure        
                
    if (verbose) write(*,*)'domain=',domain
    
    dim=PUL%L%dim
    if (verbose)  write(*,*)'Allocating bundle%L(domain)        ,dim=',dim
    ALLOCATE( bundle%L(domain)%mat(dim,dim) ) 
    bundle%L(domain)%dim=dim
    
    dim=PUL%C%dim
    if (verbose) write(*,*)'Allocating bundle%C(domain)        ,dim=',dim
    ALLOCATE( bundle%C(domain)%mat(dim,dim) ) 
    bundle%C(domain)%dim=dim
    
    dim=PUL%Zfilter%dim
    if (verbose) write(*,*)'Allocating bundle%Z(domain)        ,dim=',dim
    Bundle%Z(domain)%dim=dim
    ALLOCATE(Bundle%Z(domain)%sfilter_mat(dim,dim))
    
    dim=PUL%Yfilter%dim
    if (verbose) write(*,*)'Allocating bundle%Y(domain)        ,dim=',dim
    Bundle%Y(domain)%dim=dim
    ALLOCATE(Bundle%Y(domain)%sfilter_mat(dim,dim))
     
    local_n_conductors=bundle%n_conductors(domain) 
    
    if (is_overshield_domain) then
    
      if (verbose) write(*,*)'Copying L, C, Z, Y matrices in overshield domain, n_conductors=',local_n_conductors, &
                             overshield_n_conductors(overshield)+1
      
      do row_l=1,local_n_conductors-1 ! Note that the reference conductor is included in the conductor count here
               
        do col_l=1,local_n_conductors-1 ! Note that the reference conductor is included in the conductor count here
     
! In the global L,C structure rows and columns correspond to the global_domain conductor number
          row_g=terminal_conductor_to_global_domain_conductor(overshield_terminal_conductor(overshield,row_l))
          col_g=terminal_conductor_to_global_domain_conductor(overshield_terminal_conductor(overshield,col_l))
          
          bundle%L(domain)%mat(row_l,col_l)=PUL%L%mat(row_l,col_l)      
          bundle%C(domain)%mat(row_l,col_l)=PUL%C%mat(row_l,col_l)
          
          bundle%Z(domain)%sfilter_mat(row_l,col_l)=PUL%Zfilter%sfilter_mat(row_l,col_l)
          bundle%Y(domain)%sfilter_mat(row_l,col_l)=PUL%Yfilter%sfilter_mat(row_l,col_l)
          
! copy the domain based L and C matrices to the global L and C in the bundle structure                        
          bundle%global_L%mat(row_g,col_g)=bundle%L(domain)%mat(row_l,col_l)
          bundle%global_C%mat(row_g,col_g)=bundle%C(domain)%mat(row_l,col_l)
          
! copy the domain based Z and Y filter matrices to the global Z and Y in the bundle structure                        
          bundle%global_Z%sfilter_mat(row_g,col_g)=bundle%Z(domain)%sfilter_mat(row_l,col_l)
          bundle%global_Y%sfilter_mat(row_g,col_g)=bundle%Y(domain)%sfilter_mat(row_l,col_l)
         
        end do ! next column of matrix
          
      end do ! next row of matrix

! Added CJS 27/4/2016
      bundle%terminal_conductor_list(domain)%n_elements=local_n_conductors
      ALLOCATE( bundle%terminal_conductor_list(domain)%element(1:local_n_conductors) )
        
      do conductor=1,local_n_conductors-1   ! do reference conductor separately    
        
        bundle%terminal_conductor_list(domain)%element(conductor)=overshield_terminal_conductor(overshield,conductor)
            
      end do ! next conductor in this domain
        
! reference conductor        
      reference_conductor=global_domain_reference_conductor(domain)       ! note domain based numbering...
      bundle%terminal_conductor_list(domain)%element(local_n_conductors)=  &
                                               global_domain_conductor_to_terminal_conductor(reference_conductor)

    else
! external domain
     
      write(*,*)'External domain, n_conductors=',local_n_conductors,' dim=',local_n_conductors-1 
     
      do row_l=1,local_n_conductors-1 ! subtract 1 as the reference conductor is included
               
        do col_l=1,local_n_conductors-1 ! subtract 1 as the reference conductor is included
        
! In the global L,C structure rows and columns correspond to the bundle_domain conductor number
          row_g=terminal_conductor_to_global_domain_conductor(external_terminal_conductor(row_l))
          col_g=terminal_conductor_to_global_domain_conductor(external_terminal_conductor(col_l))
          
          if (verbose) write(*,*)'Copy L,C from PUL structure'
          bundle%L(domain)%mat(row_l,col_l)=PUL%L%mat(row_l,col_l)      
          bundle%C(domain)%mat(row_l,col_l)=PUL%C%mat(row_l,col_l)
          
          if (verbose) write(*,*)'Copy Z,Y from PUL structure'
          bundle%Z(domain)%sfilter_mat(row_l,col_l)=PUL%Zfilter%sfilter_mat(row_l,col_l)
          bundle%Y(domain)%sfilter_mat(row_l,col_l)=PUL%Yfilter%sfilter_mat(row_l,col_l)

! copy the domain based L and C matrices to the global L and C in the bundle structure                        
          if (verbose) write(*,*)'Copy L,C from bundle structure'
          bundle%global_L%mat(row_g,col_g)=bundle%L(domain)%mat(row_l,col_l)
          bundle%global_C%mat(row_g,col_g)=bundle%C(domain)%mat(row_l,col_l)
          
          if (verbose) write(*,*)'Copy Z,Y from bundle structure'
 ! copy the domain based Z and Y filter matrices to the global Z and Y in the bundle structure                        
          bundle%global_Z%sfilter_mat(row_g,col_g)=bundle%Z(domain)%sfilter_mat(row_l,col_l)
          bundle%global_Y%sfilter_mat(row_g,col_g)=bundle%Y(domain)%sfilter_mat(row_l,col_l)
         
        end do ! next column of matrix
          
      end do ! next row of matrix
      
! Added CJS 27/4/2016
      bundle%terminal_conductor_list(domain)%n_elements=local_n_conductors
      ALLOCATE( bundle%terminal_conductor_list(domain)%element(1:local_n_conductors) )
        
      do conductor=1,local_n_conductors-1   ! do reference conductor separately    
        
        bundle%terminal_conductor_list(domain)%element(conductor)=external_terminal_conductor(conductor)  
            
      end do ! next conductor in this domain
        
! reference conductor        
      reference_conductor=global_domain_reference_conductor(domain)       ! note domain based numbering...
      bundle%terminal_conductor_list(domain)%element(local_n_conductors)=  &
                                               global_domain_conductor_to_terminal_conductor(reference_conductor)
      
    end if ! external domain
    
    CALL deallocate_PUL_data(PUL)

  end do ! next external domain

! 8. Copy the cable based conductor impedance (loss) models to the bundle structure

  conductor=0  ! count of terminal conductors
  
  do cable=1,tot_n_cables
                    
      do local_conductor=1,bundle%cable(cable)%tot_n_conductors
      
        conductor=conductor+1
! copy the conductor impedance model for this cable conductor to the bundle structure        
        
        bundle%conductor_impedance(conductor)%impedance_model_type=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%impedance_model_type
        

        bundle%conductor_impedance(conductor)%radius=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%radius

        bundle%conductor_impedance(conductor)%width=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%width

        bundle%conductor_impedance(conductor)%height=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%height
                    
        bundle%conductor_impedance(conductor)%conductivity=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%conductivity
                    
        bundle%conductor_impedance(conductor)%thickness=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%thickness
              
        bundle%conductor_impedance(conductor)%Resistance_multiplication_factor=  &
                    bundle%cable(cable)%conductor_impedance(local_conductor)%Resistance_multiplication_factor
                    
        if((bundle%conductor_impedance(conductor)%impedance_model_type.EQ.impedance_model_type_filter).OR. &  
           (bundle%conductor_impedance(conductor)%impedance_model_type.EQ.impedance_model_type_cylindrical_shield)) then
                                
          bundle%conductor_impedance(conductor)%ZT_filter=   &
                 bundle%cable(cable)%conductor_impedance(local_conductor)%ZT_filter
                               
        end if
    
      end do ! next conductor on this cable
      
  end do ! next cable

! 8b. Copy the cable based conductor labels to the bundle structure

  ALLOCATE( bundle%conductor_label(1:bundle%tot_n_conductors) )
  conductor=0  ! count of terminal conductors
  
  do cable=1,tot_n_cables
                    
      do local_conductor=1,bundle%cable(cable)%tot_n_conductors
      
        conductor=conductor+1
! copy the conductor impedance model for this cable conductor to the bundle structure        
        write(conductor_string,'(I3)')conductor
        bundle%conductor_label(conductor)='Conductor number '//conductor_string//'. '//  &
                                          trim(bundle%cable(cable)%conductor_label(local_conductor))
    
      end do ! next conductor on this cable
      
  end do ! next cable
  
! 9 save numbering information required for the transfer impedance calculation

! first loop over shielded domains. The reference conductor for each domain is a shield (by definition) 
! apart from the shielded twisted pair differential mode exception
! so record this hence set out the inner domain number for shields

write(*,*)'****************************************************************'

  do domain=1,bundle%tot_n_domains-1   ! note we exclude the exterior domain here
  
    reference_conductor=global_domain_reference_conductor(domain)
    terminal_conductor=global_domain_conductor_to_terminal_conductor(reference_conductor)  
    
    write(*,*)'domain=',domain,' reference_conductor=',reference_conductor,' terminal_conductor=',terminal_conductor
    
    if (.NOT.domain_is_TP_differential_mode(domain)) then
      bundle%terminal_conductor_is_shield_flag(terminal_conductor)=.TRUE.
      bundle%terminal_conductor_to_inner_domain(terminal_conductor)=domain
    end if
    
  end do ! next conductor
  
write(*,*)'****************************************************************'
 
! work out the outer domain for each conductor and the associated reference conductor
  do conductor=1,bundle%tot_n_conductors
  
    bundle%terminal_conductor_to_outer_domain(conductor)=terminal_conductor_to_global_domain(conductor)   
    bundle%terminal_conductor_to_global_domain_conductor(conductor)=terminal_conductor_to_global_domain_conductor(conductor)   
    bundle%terminal_conductor_to_reference_terminal_conductor(conductor)=   &
           terminal_conductor_to_reference_terminal_conductor(conductor)    
           
  end do ! next conductor
    
! 10. finish up

  if (allocated( overshield_shape )) DEALLOCATE( overshield_shape)
  if (allocated( overshield_x )) DEALLOCATE( overshield_x)
  if (allocated( overshield_y )) DEALLOCATE( overshield_y)
  if (allocated( overshield_r )) DEALLOCATE( overshield_r)
  if (allocated( overshield_w )) DEALLOCATE( overshield_w)
  if (allocated( overshield_w2 )) DEALLOCATE( overshield_w2)
  if (allocated( overshield_h )) DEALLOCATE( overshield_h)
  if (allocated( overshield_domain )) DEALLOCATE( overshield_domain)
  if (allocated( overshield_reference_terminal_conductor)) DEALLOCATE( overshield_reference_terminal_conductor )
  if (allocated( overshield_to_cable_number)) DEALLOCATE( overshield_to_cable_number )
  if (allocated( global_domain_reference_conductor )) DEALLOCATE( global_domain_reference_conductor )
  DEALLOCATE( cable_reference_conductor  )
  DEALLOCATE( cable_reference_domain  )
 
! Deallocate the referencing arrays
 
  DEALLOCATE( terminal_conductor_to_cable )  
  DEALLOCATE( terminal_conductor_to_cable_local_domain )  
  DEALLOCATE( terminal_conductor_to_global_domain )  
  DEALLOCATE( terminal_conductor_to_global_domain_conductor )  
  DEALLOCATE( terminal_conductor_to_reference_global_domain_conductor )  
  DEALLOCATE( terminal_conductor_to_reference_terminal_conductor )  
  DEALLOCATE( terminal_conductor_is_reference_conductor )  
  DEALLOCATE( global_domain_conductor_to_terminal_conductor )  
  DEALLOCATE( external_terminal_conductor )

  if (allocated( domain_is_TP_differential_mode )) DEALLOCATE( domain_is_TP_differential_mode )

END SUBROUTINE create_global_domain_structure