text:
! Test program for using customized calling subroutine, use
! the Rosenbrock function from driver 1, added constants

!
! CHANGES
! 21.12.2018 - Miha Polajnar
! Initial version
!
program custom1
  use iso_fortran_env, only: rk => real64
  use l_bfgs_b_ag, only: l_bfgs_b, optim_prob_t, lbfgsb_opt_t
  implicit none
  integer, parameter :: n = 2
  integer :: i
  type :: rconst_t
    real(rk) :: a, b
  end type
  type(optim_prob_t) :: optim_prob
  type(lbfgsb_opt_t) :: lbfgsb_opt
  !
  lbfgsb_opt = lbfgsb_opt_t(iprint=-200,factr=1.e7_rk,pgtol=1.e-5_rk,&
    maxeval=500)
  !
  optim_prob%name = 'NOCEDAL EX. 17.1'
  optim_prob%n = n
  optim_prob%x = [(0._rk,i=1,n)]
  optim_prob%nbd=[(2,i=1,n)]
  allocate(optim_prob%lb(n),optim_prob%ub(n))
  optim_prob%lb = [(-5.12_rk,i=1,n)]
  optim_prob%ub = [(5.12_rk,i=1,n)]
  optim_prob%func => nocedal_ex_17_1
  optim_prob%grad => grad_nocedal_ex_17_1
  optim_prob%neqc = 1
  optim_prob%eq_const => nocedal_ex_17_1_eq_const
  optim_prob%grad_eq_const => grad_nocedal_ex_17_1_eq_const
  allocate(optim_prob%lagr(optim_prob%neqc),source=10._rk)
  allocate(optim_prob%dat,source=2._rk)
  ! Run the algorithm
  call l_bfgs_b(optim_prob,lbfgsb_opt)
contains

  real(rk) pure function nocedal_ex_17_1(x,dat)
    real(rk), intent(in) :: x(:)
    class(*), intent(in) :: dat
    nocedal_ex_17_1 = x(1) + x(2)
  end function

  pure subroutine grad_nocedal_ex_17_1(x,dat,grad)
    real(rk), intent(in) :: x(:)
    class(*), intent(in) :: dat
    real(rk), intent(out) :: grad(:)
    grad(1) = 1._rk
    grad(2) = 1._rk
  end subroutine

  pure subroutine nocedal_ex_17_1_eq_const(x,dat,ce)
    real(rk), intent(in) :: x(:)
    class(*), intent(in) :: dat
    real(rk), intent(out) :: ce(:)
    select type(dat)
    type is (real(rk))
      ce(1) = x(1)**2 + x(2)**2 - dat
    class default
      error stop 'Nocedal EX. 17.1 wrong type.'
    end select
  end subroutine

  pure subroutine grad_nocedal_ex_17_1_eq_const(x,dat,grad)
    real(rk), intent(in) :: x(:)
    class(*), intent(in) :: dat
    real(rk), intent(out) :: grad(:,:)
    grad(1,1) = 2._rk * x(1)
    grad(1,2) = 2._rk * x(2)
  end subroutine
end program

max_stars_repo_name: spule/l-fbgs-bc

max_issues_repo_name: spule/l-fbgs-bc

max_forks_repo_name: spule/l-fbgs-bc


text:
Hibbert Lumber Co. is a familyowned business offering a selection of lumber and hardware for both contractors and doityourselfers. The staff is both knowledgeable and helpful. The business has been around since 1947, established by Jack and Carol Hibbert. It originally existed in the form of a small metal shed and barn located at Third Street and I Street. This first incarnation of the business burned down in 1955. This prompted a move to its current location, where the lumber barn was destroyed by a fire in 1979, this time set by an arson arsonist.

They also own a the The Door Shop Door Shop which is located across the railroad tracks from the main location.

In the 19th century, this site was the location of one of two very large scales that weighed wagons bearing grain.  After the weighing process, the product was then shipped out via the nearby railroad.  The other scale was in the location of where The Secretariat now stands.

20051016 22:28:34 nbsp I have been continually astonished by how crappy Hibberts lumber is.  Their wood is splintered, warped, overpriced, falsely marketed.  the staff inside the store is very helpful though. Users/JosephAaronElijahStewart

20051211 02:30:29 nbsp The wood is stored outside, so you have to let it acclimate to your workshop.  Hibberts and Ace are the only hardware stores that I know of in Davis.  Hibbert has a much wider selection of hardwoods and you can pick out your own pieces. Users/OwenScharlotte

20061205 19:51:29 nbsp These guys are way more helpful than those over at Ace. Highly recommended. Users/WilliamLewis

20070324 07:21:39 nbsp Hibbert is a great place to buy building supplies and lumber. I’ve been buying from them for eight years – small projects to a complete house remodel. The staff is very knowledgeable and has great recommendations for solving problems. Their lumber is excellent quality compared to the big box stores. Returns are easy. Also, they offer free delivery in town  a real plus if you dont have a truck. Generally, their prices are about 510% more than the big box stores, but it is worth it for the excellent service, better quality lumber, and shorter drive. Users/MkLindq

20081123 08:36:21 nbsp Corrugated roofing material at Hibbert lumber is the same price as it is in Home Depot.  This nugget wont be relevant if you want to buy some on a Sunday though. Users/RobiPochapin

20090517 20:48:10 nbsp I am a big fan of the hardware store here. The selection isnt the biggest, but they always have what youre looking for. Staff is knowledgeable and friendly. After having just a few too many negative experiences at Ace, I would rather put a Sunday project on hold because Hibberts is closed than head over to an alternative store. Prices may be slightly higher, but youll save with staff insight and reliable products in the end. Users/MichaelAWoods

20091103 21:20:28 nbsp This is how lumber yards and hardware shops used to be: The place to find people who know the building trades. When I ask for a powertool part, it can be ordered in minutes. Something that cant be done at Ace or Home Depot. When I ask for options for lumber, I get proper answers. Hibberts has treated me well, and their prices typically are not higher than anywhere else for the same quality product. Users/DarellDickey

20100223 23:52:25 nbsp Davis is extremely fortunate to have two excellent hardware stores.  Hibbert Lumber has always had both the items Ive needed and people at the front counter who can take me right to them.  On my last trip I purchased four cap nuts.  It was great to buy only the four I needed (and not two pastic bubble packs of 3 each).  It was also nice to be able to hold the various sizes to my part so that I knew I was buying the right size.  Users/RMS

20100718 13:02:04 nbsp Hibberts is always on point and super friendly! A charming shop, we got a 4x8 3/4{{{}}} ACX ply sheet and some 4x4s for a little frame work bed install for the back of truck and everyone was friendly and awesome. You cant go wrong if your a contractor headed to hibberts, but if your the little guy who is clueless they are a solid bet. Stay local! Users/StevenDaubert

20110715 21:01:06 nbsp An excellent family business that really cares about its customers and its community.  They are always friendly and ready to help.  The guys in the yard are extremely helpful.  I have had several projects that couldnt have gotten finished without them.  Great way to support the local economy. Users/SilvaJ

20120721 13:56:46 nbsp I went to Ace around this time last summer asking them how to remove a caster that I had in my Ikea desk (I had to completely take the desk apart in order to get it out of my apartment). The dude at Ace told me I needed to buy a power drillwhich I thought was a load of crap. Then I went to Hibbert, showed them a picture of the caster in my desk (like I did at Ace), and the owner physically finds me a caster and shows me how to remove it with a screwdriver. It was so easy to remove and saved my sanity. So definitely not going to Ace anymore for help with stuff like that.  Users/smcosta

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE DRWPLB (NUMBUN, BUNDIS, BUNDIF)

C  DRWPLB draws a set of actual polylines using the specified bundle
C  values, and a set of expected polylines, using the attribute
C  values returned by <inquire predefined polyline representation>.
C  It draws the expected polyline incorrectly at the specified
C  position.
C
C  Input parameters:
C    NUMBUN : number of entries in BUNDIS
C    BUNDIS : list of polyline bundle indices to be displayed
C    BUNDIF : position of polyline to be drawn incorrectly

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

      INTEGER  NUMBUN, BUNDIS(NUMBUN), ERRIND, PDLT,PDCI, BUNDIF, IX

      REAL     XA(2),YA(2), YINCR,YLOC, PDLW

C  draw and label actual; use individual attributes
      CALL SETASF (PBUNDL)
      XA(1) = 0.25
      XA(2) = 0.55
      YINCR = 0.8/NUMBUN
      YLOC  = 0.9
      CALL NUMLAB (NUMBUN, 0.2, YLOC, YINCR)
      DO 100 IX = 1,NUMBUN
         YA(1) = YLOC
         YA(2) = YLOC
         CALL PSPLI (BUNDIS(IX))
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
100   CONTINUE

C  now draw expected results - with individual attributes from WDT
C  use individual attributes
      CALL SETASF (PINDIV)
      XA(1) = .6
      XA(2) = .9
      YLOC  = .9

      DO 200 IX = 1, NUMBUN
         CALL PQPPLR (SPECWT, BUNDIS(IX), ERRIND, PDLT, PDLW, PDCI)
         CALL CHKINQ ('pqpplr', ERRIND)
         IF (IX .EQ. BUNDIF) THEN
            PDLT = MOD(PDLT,4) + 1
            PDLW = PDLW * 1.5
            PDCI = MOD(PDCI,5) + 1
         ENDIF
         CALL PSLN   (PDLT)
         CALL PSLWSC (PDLW)
         CALL PSPLCI (PDCI)
         YA(1) = YLOC
         YA(2) = YLOC
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
200   CONTINUE

      END


module CFDElementM

  use UtilitiesM
  use SparseKit

  use Triangle2D3NodeM
  use Triangle2D6NodeM
  use Quadrilateral2D4NodeM
  use Quadrilateral2D8NodeM

  use IntegratorPtrM

  use LeftHandSideM
  use ProcessInfoM

  use PointM
  use NodeM
  use NodePtrM

  use SourceM
  use SourcePtrM

  use ElementM

  use CFDMaterialM

  implicit none

  private
  public :: CFDElementDT, cfdElement, initGeometries, allocateMass, allocateStabMat, zeroStabMat &
       , lumpedMass

  type, extends(ElementDT) :: CFDElementDT
     class(CFDMaterialDT), pointer                :: material
     real(rkind)    , dimension(:,:), allocatable :: dNdx, dNdy
     real(rkind)    , dimension(:)  , allocatable :: jacobianDet
     real(rkind)    , dimension(:,:), allocatable :: Tau
     real(rkind)    , dimension(:)  , allocatable :: nu
   contains
     procedure, public  :: init
     procedure, public  :: calculateLHS
     procedure, public  :: calculateRHS
     procedure, public  :: calculateLocalSystem
     procedure, public  :: calculateResults
     procedure, public  :: calculateDT
     procedure, public  :: calculateMass
     procedure, public  :: calculateStabMat
     procedure, private :: setupIntegration
     procedure, private :: getValuedSource
     procedure, private :: vectorCalculus
  end type CFDElementDT

  interface cfdElement
     procedure :: constructor
  end interface cfdElement

  real(rkind)                 , dimension(:)  , allocatable :: lumpedMass
  real(rkind)                 , dimension(:,:), allocatable :: stabMat
  type(Triangle2D3NodeDT)     , target        , save        :: myTriangle2D3Node
  type(Triangle2D6NodeDT)     , target        , save        :: myTriangle2D6Node
  type(Quadrilateral2D4NodeDT), target        , save        :: myQuadrilateral2D4Node
  type(Quadrilateral2D8NodeDT), target        , save        :: myQuadrilateral2D8Node

contains

  type(CFDElementDT) function constructor(id, node, material)
    implicit none
    integer(ikind)                          , intent(in) :: id
    type(NodePtrDT)           , dimension(:), intent(in) :: node
    class(CFDMaterialDT), target            , intent(in) :: material
    call constructor%init(id, node, material)
  end function constructor

  subroutine init(this, id, node, material)
    implicit none
    class(CFDElementDT)                     , intent(inout) :: this
    integer(ikind)                          , intent(in)    :: id
    type(NodePtrDT)           , dimension(:), intent(in)    :: node
    class(CFDMaterialDT), target            , intent(in)    :: material
    this%id = id
    this%node = node
    this%material => material
    if(size(node) == 3) then
       this%geometry => myTriangle2D3Node
    else if(size(node) == 4) then
       this%geometry => myQuadrilateral2D4Node
    else if(size(node) == 6) then
       this%geometry => myTriangle2D6Node
    else if(size(node) == 8) then
       this%geometry => myQuadrilateral2D8Node
    end if
    allocate(this%source(1))
    call this%vectorCalculus()
  end subroutine init

  subroutine initGeometries(nGauss)
    implicit none
    integer(ikind), intent(in) :: nGauss
    myTriangle2D3Node = triangle2D3Node(nGauss)
    myTriangle2D6Node = triangle2D6Node(nGauss)
    myQuadrilateral2D4Node = quadrilateral2D4Node(nGauss)
    myQuadrilateral2D8Node = quadrilateral2D8Node(nGauss)
  end subroutine initGeometries

  subroutine calculateLocalSystem(this, processInfo, lhs, rhs)
    implicit none
    class(CFDElementDT)                       , intent(inout) :: this
    type(ProcessInfoDT)                       , intent(inout) :: processInfo
    type(LeftHandSideDT)                      , intent(inout) :: lhs
    real(rkind)    , dimension(:), allocatable, intent(inout) :: rhs
    type(IntegratorPtrDT)                                     :: integrator
    integer(ikind)                                            :: i, j, k, elemID
    integer(ikind)                                            :: nNode, nDof, iNodeID
    real(rkind)    , dimension(:,:)  , allocatable            :: U, theta
    real(rkind)                                               :: gamma, V_sq
    real(rkind)                                               :: bi, ci, v1, v2, e, rho
    real(rkind)                                               :: thetaK(4), AiUi(4), AiUi_theta(4)
    real(rkind)                                               :: U_k(4), theta_k(4), Ux(4), Uy(4)
    real(rkind)                                               :: A1AiUi_theta(4), A2AiUi_theta(4)
    real(rkind)                                               :: K1jUj(4), K2jUj(4), lambda, mu, temp
    real(rkind)                                               :: Cv, T_inf, r
    !Esto de la línea siguiente no hace nada creo..
    !OMP DEFAULT(PRIVATE) SHARED(lambda, mu, temp, Cf, t_inf, gamma, lumpedMass, stabMat)
    nNode = this%getnNode()
    nDof  = 4
    allocate(U(nDof,nNode))
    allocate(rhs(nNode*nDof))
    integrator  = this%getIntegrator()

    lambda  = this%material%k
    mu      = this%material%mu
    temp    = 0._rkind
    Cv      = this%material%Cv
    T_inf   = this%material%T_inf
    r       = this%material%R

    gamma   = this%material%gamma
    elemID  = this%getID()

    rhs = 0._rkind
    do i = 1, nNode
       do j = 1, nDof
          !U(j,i) = this%node(i)%ptr%dof(j)%getVal()
          U(j,i) = this%node(i)%ptr%dof(j)%val
!!$          write(11,*) 'elem -> ', this%getID()
!!$          write(11,*) 'dof -> ', this%node(i)%ptr%dof(j)%getVal()
       end do
       temp = temp + (((U(4,i)/U(1,i))-0.5*((U(2,i)/U(1,i))**2+(U(3,i)/U(1,i))**2))/Cv)/nNode
       !temp = temp + (((U(4,i)/U(1,i))-0.5*((U(2,i)/U(1,i))**2 &
       !     +(U(3,i)/U(1,i))**2))*(gamma-1._rkind)/r)/nNode

!!$       write(11,*) 'elem -> ', elemID
!!$       write(11,*) 'node -> ', i
!!$       write(11,*) 'U(2,i) -> ', U(2,i)
!!$       write(11,*) 'U(3,i) -> ', U(3,i)
!!$       write(11,*) 'temp -> ', temp
    end do

    lambda = lambda*(temp/T_inf)**1.5*(T_inf + 194)/(temp + 194)
!!$    write(11,*) 'temp -> ', temp
!!$    write(11,*) 'lambda -> ', lambda

    do k = 1, integrator%getIntegTerms()

       Ux      = 0._rkind
       Uy      = 0._rkind
       U_k     = 0._rkind
       theta_k = 0._rkind
       do i = 1, nNode
          iNodeID = this%getNodeID(i)
          U_k     = U_k     + integrator%ptr%shapeFunc(k,i)*U(:,i)
          theta_k = theta_k + integrator%ptr%shapeFunc(k,i)*stabMat(:,iNodeID)
          Ux = Ux + U(:,i)*this%dNdx(i,k)
          Uy = Uy + U(:,i)*this%dNdy(i,k)
!!$          write(11,*) 'iNodeID -> ', iNodeID
!!$          write(11,*) 'lumpedMass -> ', lumpedMass(iNodeID)
!!$          write(11,*) 'stabMat ->', stabMat(1:nNode,iNodeID)
!!$          write(11,*) 'theta_k ->', theta_k
       end do

       rho  = U_k(1)
       v1   = U_k(2)/rho
       v2   = U_k(3)/rho
       e    = U_k(4)/rho
       V_sq = v1**2 + v2**2
       AiUi(1) = Ux(2) + Uy(3)
       AiUi(2) = (1._rkind/2._rkind)*Ux(1)*(V_sq*(gamma - 1) - 2*v1**2) - Ux(2)*v1*( &
            gamma - 3) - Ux(3)*v2*(gamma - 1) + Ux(4)*(gamma - 1) - Uy(1)*v1*v2 + &
            Uy(2)*v2 + Uy(3)*v1
       AiUi(3) = -Ux(1)*v1*v2 + Ux(2)*v2 + Ux(3)*v1 + (1._rkind/2._rkind)*Uy(1)*(V_sq*( &
            gamma - 1) - 2*v2**2) - Uy(2)*v1*(gamma - 1) - Uy(3)*v2*(gamma - 3) + &
            Uy(4)*(gamma - 1)
       AiUi(4) = Ux(1)*v1*(V_sq*(gamma - 1) - e*gamma) - 1._rkind/2._rkind*Ux(2)*(V_sq &
            *(gamma - 1) - 2*e*gamma + 2*v1**2*(gamma - 1)) - Ux(3)*v1*v2*( &
            gamma - 1) + Ux(4)*gamma*v1 + Uy(1)*v2*(V_sq*(gamma - 1) - e*gamma) - &
            Uy(2)*v1*v2*(gamma - 1) - 1._rkind/2._rkind*Uy(3)*(V_sq*(gamma - 1) - 2*e* &
            gamma + 2*v2**2*(gamma - 1)) + Uy(4)*gamma*v2

       AiUi_theta = theta_k + AiUi

       A1AiUi_theta(1) = AiUi_theta(2)
       A1AiUi_theta(2) = v1*(-gamma + 3)*AiUi_theta(2) - v2*(gamma - 1)* &
            AiUi_theta(3) + (gamma - 1)*AiUi_theta(4) + ((1._rkind/2._rkind)* &
            V_sq*(gamma - 1) - v1**2)*AiUi_theta(1)
       A1AiUi_theta(3) = -v1*v2*AiUi_theta(1) + v1*AiUi_theta(3) + v2* &
            AiUi_theta(2)
       A1AiUi_theta(4) = gamma*v1*AiUi_theta(4) - v1*v2*(gamma - 1)* &
            AiUi_theta(3) + v1*(V_sq*(gamma - 1) - e*gamma)*AiUi_theta(1) &
            + (-1._rkind/2._rkind*V_sq*(gamma - 1) + e*gamma - v1**2*(gamma - 1 &
            ))*AiUi_theta(2)

       A2AiUi_theta(1) = AiUi_theta(3)
       A2AiUi_theta(2) = -v1*v2*AiUi_theta(1) + v1*AiUi_theta(3) + v2* &
            AiUi_theta(2)
       A2AiUi_theta(3) = -v1*(gamma - 1)*AiUi_theta(2) + v2*(-gamma + 3)* &
            AiUi_theta(3) + (gamma - 1)*AiUi_theta(4) + ((1._rkind/2._rkind)* &
            V_sq*(gamma - 1) - v2**2)*AiUi_theta(1)
       A2AiUi_theta(4) = gamma*v2*AiUi_theta(4) - v1*v2*(gamma - 1)* &
            AiUi_theta(2) + v2*(V_sq*(gamma - 1) - e*gamma)*AiUi_theta(1) &
            + (-1._rkind/2._rkind*V_sq*(gamma - 1) + e*gamma - v2**2*(gamma - 1 &
            ))*AiUi_theta(3)

       A1AiUi_theta(1) = A1AiUi_theta(1)*this%Tau(1,k)
       A1AiUi_theta(2:3) = A1AiUi_theta(2:3)*this%Tau(2,k)
       A1AiUi_theta(4) = A1AiUi_theta(4)*this%Tau(3,k)

       A2AiUi_theta(1) = A2AiUi_theta(1)*this%Tau(1,k)
       A2AiUi_theta(2:3) = A2AiUi_theta(2:3)*this%Tau(2,k)
       A2AiUi_theta(4) = A2AiUi_theta(4)*this%Tau(3,k)

       K1jUj(1) = 0._rkind
       K1jUj(2) = (2._rkind/3._rkind)*mu*(-2*Ux(1)*v1 + 2*Ux(2) + Uy(1)*v2 - Uy(3))/rho
       K1jUj(3) = mu*(-Ux(1)*v2 + Ux(3) - Uy(1)*v1 + Uy(2))/rho
       K1jUj(4) = (1._rkind/3._rkind)*(Cv*mu*(-Uy(1)*v1*v2 + 3*Uy(2)*v2 - 2*Uy(3)*v1) &
            - Ux(1)*(Cv*mu*(3*V_sq + v1**2) - 3*lambda*(V_sq - e)) + Ux(2)*v1*(4*Cv*mu &
            - 3*lambda) + 3*Ux(3)*v2*(Cv*mu - lambda) + 3*Ux(4)*lambda)/(Cv*rho)
       K2jUj(1) = 0._rkind
       K2jUj(2) = mu*(-Ux(1)*v2 + Ux(3) - Uy(1)*v1 + Uy(2))/rho
       K2jUj(3) = (2._rkind/3._rkind)*mu*(Ux(1)*v1 - Ux(2) - 2*Uy(1)*v2 + 2*Uy(3))/rho
       K2jUj(4) = (1._rkind/3._rkind)*(Cv*mu*(-Ux(1)*v1*v2 - 2*Ux(2)*v2 + 3*Ux(3)*v1) &
            - Uy(1)*(Cv*mu*(3*V_sq + v2**2) - 3*lambda*(V_sq - e)) + 3*Uy(2)*v1*(Cv*mu &
            - lambda) + Uy(3)*v2*(4*Cv*mu - 3*lambda) + 3*Uy(4)*lambda)/(Cv*rho)

       do i = 1, nNode
          iNodeID = this%getNodeID(i)
          rhs(i*nDof-3:i*nDof) = rhs(i*nDof-3:i*nDof) - (integrator%getShapeFunc(k,i)*AiUi     &
               + this%dNdx(i,k)*A1AiUi_theta + this%dNdy(i,k)*A2AiUi_theta                     &
               + this%nu(k)*(this%dNdx(i,k)*Ux + this%dNdy(i,k)*Uy)                            &
               + (this%dNdx(i,k)*K1jUj + this%dNdy(i,k)*K2jUj))                                &
               * this%jacobianDet(k)*integrator%getWeight(k)/lumpedMass(iNodeID)
       end do
    end do
    deallocate(U)
  end subroutine calculateLocalSystem

  subroutine calculateDT(this, processInfo)
    implicit none
    class(CFDElementDT) , intent(inout)        :: this
    class(ProcessInfoDT), intent(inout)        :: processInfo
    integer(ikind)                             :: i, nNode
    real(rkind)                                :: alpha, deltaTU
    real(rkind)                                :: V, deltaTC, dt, fmu, et, pe, hh
    real(rkind)                                :: Vx, Vy, VxMAX, VyMAX, fSafe
    type(IntegratorPtrDT)                      :: integrator
    integrator  = this%getIntegrator()
    nNode       = this%getnNode()
    fSafe       = this%material%fSafe
    VxMAX       = 0._rkind
    VyMAX       = 0._rkind
    do i = 1, nNode
       Vx = abs(this%node(i)%ptr%dof(2)%getVal()/this%node(i)%ptr%dof(1)%getVal())
       Vy = abs(this%node(i)%ptr%dof(3)%getVal()/this%node(i)%ptr%dof(1)%getVal())
       if (Vx.gt.VxMAX) VxMAX = Vx
       if (Vy.gt.VyMAX) VyMAX = Vy
    end do
    V  = (VxMAX**2+VyMAX**2)**0.5_rkind

    hh = sqrt(2._rkind*this%geometry%getLenght(this%node)) !Esto es solo así para triang lin..
    fmu = 0.001_rkind
    et  = fmu
    pe  = (v*hh)/(2._rkind*et)
    alpha = min(pe/3._rkind,1._rkind)
    deltaTU = 1._rkind/(4._rkind*et/hh**2._rkind+alpha*v/hh)
    deltaTC = 1._rkind/(4._rkind*et/hh**2._rkind)
    dt = fSafe/(1._rkind/deltaTC+1._rkind/deltaTU)
    call processInfo%setMinimumDT(dt)

!!$    do i = 1, integrator%getIntegTerms()
!!$       alpha   = min((V*sqrt(this%jacobianDet(i)))/(2._rkind*0.001_rkind)/3._rkind,1._rkind)
!!$       deltaTU = 1._rkind/(4._rkind*0.001_rkind/sqrt(this%jacobianDet(i))**2._rkind &
!!$            +ALPHA*V/sqrt(this%jacobianDet(i)))
!!$       deltaTC = 1._rkind/(4._rkind*0.001_rkind/sqrt(this%jacobianDet(i))**2._rkind)
!!$       dt      = fSafe/(1._rkind/deltaTC+1._rkind/deltaTU)
!!$       call processInfo%setMinimumDt(dt)
!!$    end do
  end subroutine calculateDT

  subroutine allocateMass(n)
    implicit none
    integer(ikind), intent(in) :: n
    if(.not.allocated(lumpedMass)) allocate(lumpedMass(n))
    lumpedMass = 0._rkind
  end subroutine allocateMass

  subroutine calculateMass(this)
    implicit none
    class(CFDElementDT)                    , intent(inout) :: this
    type(IntegratorPtrDT)                                  :: integrator
    integer(ikind)                                         :: i, j, k, iNodeID
    integer(ikind)                                         :: nNode, nDof
    real(rkind)                                            :: adder, area
    nNode = this%getnNode()
!!$    integrator = this%getIntegrator()
!!$    do i = 1, nNode
!!$       adder = 0._rkind
!!$       do j = 1, nNode
!!$          do k = 1, integrator%getIntegTerms()
!!$             adder = adder + (integrator%getWeight(k)*this%jacobianDet(k)      &
!!$                  *integrator%getShapeFunc(k,i)*integrator%getShapeFunc(k,j))
!!$             write(11,*) 'k -> ', k
!!$             write(11,*) 'getWeight -> ', integrator%getWeight(k)
!!$             write(11,*) 'jacobianDet -> ', this%jacobianDet(k)
!!$             write(11,*) 'sfi -> ', integrator%getShapeFunc(k,i)
!!$             write(11,*) 'sfj -> ', integrator%getShapeFunc(k,j)
!!$             write(11,*) 'term -> ', integrator%getWeight(k)*this%jacobianDet(k)       &
!!$                  *integrator%getShapeFunc(k,i)*integrator%getShapeFunc(k,j)
!!$          end do
!!$       end do
!!$       iNodeID = this%getNodeID(i)
!!$       lumpedMass(iNodeID) = lumpedMass(iNodeID) + adder
!!$       write(11,*) 'CONTRIBUTION TO LUMPED MASS ON'
!!$       write(11,'(A,I0)') 'iNodeID = ', iNodeID
!!$       write(11,*) 'is -------> ', adder
!!$    end do
!!$
    area = this%geometry%getLenght(this%node)
    do i = 1, nNode
       iNodeID = this%getNodeID(i)
       lumpedMass(iNodeID) = lumpedMass(iNodeID) + area/nNode
    end do

  end subroutine calculateMass

  subroutine allocateStabMat(n,m)
    implicit none
    integer(ikind), intent(in) :: n
    integer(ikind), intent(in) :: m
    if(.not.allocated(stabMat)) allocate(stabMat(n,m))
    stabMat = 0._rkind
  end subroutine allocateStabMat

  subroutine zeroStabMat()
    implicit none
    stabMat = 0._rkind
  end subroutine zeroStabMat

!!$  subroutine calculateStabMat(this, processInfo)
!!$    implicit none
!!$    class(CFDElementDT)                  , intent(inout) :: this
!!$    class(processinfoDT)                 , intent(inout) :: processInfo
!!$    integer(ikind)                                       :: i, j, k, iNodeID
!!$    integer(ikind)                                       :: nNode, nDof
!!$    real(rkind)          , dimension(:,:)  , allocatable :: U
!!$    real(rkind)          , dimension(:)    , allocatable :: T
!!$    real(rkind)                                          :: U_loc(4)
!!$    real(rkind)                                          :: Ux(4), Uy(4), AiUi(4)
!!$    real(rkind)                                          :: v1, v2, V_sq, e
!!$    real(rkind)                                          :: rho, h_jgn, h_rgn, h_rgne
!!$    real(rkind)                                          :: deltaTC, dtmin, Tinf
!!$    real(rkind)                                          :: Vx, Vy, gamma, Cv, Vc, dt2, cte
!!$    real(rkind)                                          :: fmu, smu, bi, ci, rrr, zzz, shoc
!!$    real(rkind)                                          :: rhoinf, Tau, t_sugn1, t_sugn2, t_sugn3
!!$    real(rkind)                                          :: term1, term2, tr1, tau_sung3, tau_sung3_e
!!$    real(rkind)                                          :: h_rgn1, h_rgn2, resumen, term_1, term_2
!!$    real(rkind)                                          :: R, temp, drx, dry, dr2, rjy, dtx, dty
!!$    real(rkind)                                          :: rtx, rty, dux, duy, du2, rux, ruy, rjx
!!$    type(IntegratorPtrDT)                                :: integrator
!!$    nNode = this%getnNode()
!!$    nDof  = 4
!!$    allocate(U(4,nNode), T(nNode))
!!$    integrator  = this%getIntegrator()
!!$    gamma       = this%material%gamma
!!$    do i = 1, nNode
!!$       do j = 1, nDof
!!$          U(j,i) = this%node(i)%ptr%dof(j)%getVal()
!!$       end do
!!$    end do
!!$    do k = 1, integrator%getIntegTerms()
!!$       Ux          = 0._rkind
!!$       Uy          = 0._rkind
!!$       U_loc       = 0._rkind
!!$       do i = 1, nNode
!!$          U_loc = U_loc + integrator%ptr%shapeFunc(k,i)*U(:,i)
!!$          Ux = Ux + U(:,i)*this%dNdx(i,k)
!!$          Uy = Uy + U(:,i)*this%dNdy(i,k)
!!$       end do
!!$       rho  = U_loc(1)
!!$       v1   = U_loc(2)/rho
!!$       v2   = U_loc(3)/rho
!!$       e    = U_loc(4)/rho
!!$       V_sq = v1**2 + v2**2
!!$       AiUi(1) = Ux(2) + Uy(3)
!!$       AiUi(2) = (1._rkind/2._rkind)*Ux(1)*(V_sq*(gamma - 1) - 2*v1**2) - Ux(2)*v1*( &
!!$            gamma - 3) - Ux(3)*v2*(gamma - 1) + Ux(4)*(gamma - 1) - Uy(1)*v1*v2 + &
!!$            Uy(2)*v2 + Uy(3)*v1
!!$       AiUi(3) = -Ux(1)*v1*v2 + Ux(2)*v2 + Ux(3)*v1 + (1._rkind/2._rkind)*Uy(1)*(V_sq*( &
!!$            gamma - 1) - 2*v2**2) - Uy(2)*v1*(gamma - 1) - Uy(3)*v2*(gamma - 3) + &
!!$            Uy(4)*(gamma - 1)
!!$       AiUi(4) = Ux(1)*v1*(V_sq*(gamma - 1) - e*gamma) - 1._rkind/2._rkind*Ux(2)*(V_sq &
!!$            *(gamma - 1) - 2*e*gamma + 2*v1**2*(gamma - 1)) - Ux(3)*v1*v2*( &
!!$            gamma - 1) + Ux(4)*gamma*v1 + Uy(1)*v2*(V_sq*(gamma - 1) - e*gamma) - &
!!$            Uy(2)*v1*v2*(gamma - 1) - 1._rkind/2._rkind*Uy(3)*(V_sq*(gamma - 1) - 2*e* &
!!$            gamma + 2*v2**2*(gamma - 1)) + Uy(4)*gamma*v2
!!$       do i = 1, nNode
!!$          iNodeID = this%getNodeID(i)
!!$          stabMat(:,iNodeID) = stabMat(:,iNodeID)                       &
!!$               - integrator%getShapeFunc(k,i)*AiUi                       &
!!$               * this%jacobianDet(k)*integrator%getWeight(k)/lumpedMass(iNodeID)
!!$       end do
!!$    end do
!!$    !calculo tau y nu
!!$    dtmin       = processInfo%getDt()
!!$    cte         = this%material%constant
!!$    R           = this%material%R
!!$    Cv          = this%material%Cv
!!$    Tinf        = this%material%T_inf
!!$    rhoinf      = this%material%rho
!!$
!!$    rho         = 0._rkind
!!$    Vx          = 0._rkind
!!$    Vy          = 0._rkind
!!$    temp        = 0._rkind
!!$    do i = 1, nNode
!!$       rho = rho +  U(1,i)
!!$       Vx  = Vx  + (U(2,i)/U(1,i))
!!$       Vy  = Vy  + (U(3,i)/U(1,i))
!!$    end do
!!$    rho = rho/nNode
!!$    Vx  = Vx/nNode
!!$    Vy  = Vy/nNode
!!$    V2  = sqrt(Vx*Vx+Vy*Vy)
!!$
!!$    do i = 1, nNode
!!$       T(i) = ((U(4,i)/U(1,i))-0.5_rkind*((U(2,i)/U(1,i))**2+(U(3,i)/U(1,i))**2))/Cv
!!$    end do
!!$    do i = 1, nNode
!!$       temp = temp + T(i)
!!$    end do
!!$    temp = temp/nNode
!!$    Vc   = sqrt(gamma*R*temp)
!!$
!!$    do k = 1, integrator%getIntegTerms()
!!$       Tau         = 0._rkind
!!$       drx         = 0._rkind
!!$       dry         = 0._rkind
!!$       dTx         = 0._rkind
!!$       dTy         = 0._rkind
!!$       dUx         = 0._rkind
!!$       dUy         = 0._rkind
!!$       h_rgne      = 1.d-10
!!$       h_rgn       = 1.d-10
!!$       h_jgn       = 1.d-10
!!$
!!$
!!$       do i = 1, nNode
!!$          drx  = drx + U(1,i)*this%dNdx(i,k)
!!$          dry  = dry + U(1,i)*this%dNdy(i,k)
!!$       end do
!!$       dr2  = sqrt(drx*drx+dry*dry)+1.d-20
!!$       rjx  = drx/dr2
!!$       rjy  = dry/dr2
!!$       do i = 1, nNode
!!$          dTx  = dTx + T(i)*this%dNdx(i,k)
!!$          dTy  = dTy + T(i)*this%dNdy(i,k)
!!$       end do
!!$       dT2  = sqrt(dTx*dTx+dTy*dTy)+1.d-20
!!$       rTx  = dTx/dT2
!!$       rTy  = dTy/dT2
!!$       do i = 1, nNode
!!$          dUx  = dUx + V2*this%dNdx(i,k)
!!$          dUy  = dUy + V2*this%dNdy(i,k)
!!$       end do
!!$       dU2  = sqrt(dUx*dUx+dUy*dUy)+1.d-20
!!$       rUx  = dUx/dU2
!!$       rUy  = dUy/dU2
!!$       smu  = 110._rkind
!!$       fmu  = 0.41685_rkind*(temp/Tinf)**1.5d0*(Tinf+smu)/(temp+smu)
!!$       do i = 1, nNode
!!$          term_1 = abs(Vx *this%dNdx(i,k)+Vy *this%dNdy(i,k))
!!$          term_2 = abs(rjx*this%dNdx(i,k)+rjy*this%dNdy(i,k))
!!$          h_rgn1 = abs(rtx*this%dNdx(i,k)+rty*this%dNdy(i,k))
!!$          h_rgn2 = abs(rux*this%dNdx(i,k)+ruy*this%dNdy(i,k))
!!$          Tau    = Tau+term_1+term_2*Vc
!!$          h_rgne = h_rgne+h_rgn1
!!$          h_rgn  = h_rgn+h_rgn2
!!$          h_jgn  = h_jgn+term_2
!!$       end do
!!$       Tau     = 1._rkind/Tau
!!$       h_rgne  = 2._rkind/h_rgne
!!$       h_rgn   = 2._rkind/h_rgn
!!$       if (h_rgn .gt. 1.d3) h_rgn = 0._rkind
!!$       h_jgn   = 2._rkind/h_jgn
!!$       if (h_jgn .gt. 1.d10) h_jgn = 0._rkind
!!$       tr1     = dr2*h_jgn/rho
!!$       zzz     = h_jgn/2._rkind
!!$       shoc    = (sqrt(tr1)+tr1**2)*Vc*2*zzz
!!$       resumen = (1._rkind/Tau)**2._rkind +(2._rkind/dtmin)**2._rkind
!!$       rrr     = resumen**(-0.5_rkind)
!!$       t_sugn1 = rrr
!!$       t_sugn2 = rrr
!!$       t_sugn3 = rrr
!!$       if(fmu.ne.0._rkind)then
!!$          tau_sung3   = h_rgn**2._rkind/(4._rkind*fmu/rhoinf)
!!$          tau_sung3_e = h_rgne**2._rkind/(4._rkind*fmu/rhoinf)
!!$          t_sugn2     = (resumen+1._rkind/tau_sung3**2._rkind)**(-0.5_rkind)
!!$          t_sugn3     = (resumen+1._rkind/tau_sung3_e**2._rkind)**(-0.5_rkind)
!!$       end if
!!$       this%Tau(1,k) = t_sugn1
!!$       this%Tau(2,k) = t_sugn2
!!$       this%Tau(3,k) = t_sugn3
!!$       this%nu(k)    = shoc*cte
!!$    end do
!!$    deallocate(U, T)
!!$  end subroutine calculateStabMat

  subroutine calculateStabMat(this, processInfo)
    implicit none
    class(CFDElementDT)                  , intent(inout) :: this
    class(processinfoDT)                 , intent(inout) :: processInfo
    integer(ikind)                                       :: i, j, k, iNodeID
    integer(ikind)                                       :: nNode, nDof
    real(rkind)          , dimension(:,:)  , allocatable :: U
    real(rkind)          , dimension(:)    , allocatable :: T
    real(rkind)                                          :: U_loc(4)
    real(rkind)                                          :: Ux(4), Uy(4), AiUi(4)
    real(rkind)                                          :: v1, v2, V_sq, e
    real(rkind)                                          :: rho, h_jgn, h_rgn, h_rgne
    real(rkind)                                          :: deltaTC, dtmin, Tinf
    real(rkind)                                          :: Vx, Vy, gamma, Cv, Vc, dt2, cte
    real(rkind)                                          :: fmu, smu, bi, ci, rrr, zzz, shoc
    real(rkind)                                          :: rhoinf, Tau, t_sugn1, t_sugn2, t_sugn3
    real(rkind)                                          :: term1, term2, tr1, tau_sung3, tau_sung3_e
    real(rkind)                                          :: h_rgn1, h_rgn2, resumen, term_1, term_2
    real(rkind)                                          :: R, temp, drx, dry, dr2, rjy, dtx, dty
    real(rkind)                                          :: rtx, rty, dux, duy, du2, rux, ruy, rjx
    type(IntegratorPtrDT)                                :: integrator
    nNode = this%getnNode()
    nDof  = 4
    allocate(U(4,nNode), T(nNode))
    integrator  = this%getIntegrator()
    gamma       = this%material%gamma
    dtmin       = processInfo%getDt()
    cte         = this%material%constant
    R           = this%material%R
    Cv          = this%material%Cv
    Tinf        = this%material%T_inf
    rhoinf      = this%material%rho
    temp        = 0._rkind
    do i = 1, nNode
       do j = 1, nDof
          U(j,i) = this%node(i)%ptr%dof(j)%getVal()
       end do
       T(i) = ((U(4,i)/U(1,i))-0.5d0*((U(2,i)/U(1,i))**2+(U(3,i)/U(1,i))**2))/Cv
       temp = temp + T(i)/nNode
    end do
    Vc   = sqrt(gamma*R*temp)
    do k = 1, integrator%getIntegTerms()
       Ux          = 0._rkind
       Uy          = 0._rkind
       U_loc       = 0._rkind
       do i = 1, nNode
          U_loc = U_loc + integrator%ptr%shapeFunc(k,i)*U(:,i)
          Ux = Ux + U(:,i)*this%dNdx(i,k)
          Uy = Uy + U(:,i)*this%dNdy(i,k)
       end do
       rho  = U_loc(1)
       v1   = U_loc(2)/rho
       v2   = U_loc(3)/rho
       e    = U_loc(4)/rho
       V_sq = v1**2 + v2**2
       !StabMat:
       AiUi(1) = Ux(2) + Uy(3)
       AiUi(2) = (1._rkind/2._rkind)*Ux(1)*(V_sq*(gamma - 1) - 2*v1**2) - Ux(2)*v1*( &
            gamma - 3) - Ux(3)*v2*(gamma - 1) + Ux(4)*(gamma - 1) - Uy(1)*v1*v2 + &
            Uy(2)*v2 + Uy(3)*v1
       AiUi(3) = -Ux(1)*v1*v2 + Ux(2)*v2 + Ux(3)*v1 + (1._rkind/2._rkind)*Uy(1)*(V_sq*( &
            gamma - 1) - 2*v2**2) - Uy(2)*v1*(gamma - 1) - Uy(3)*v2*(gamma - 3) + &
            Uy(4)*(gamma - 1)
       AiUi(4) = Ux(1)*v1*(V_sq*(gamma - 1) - e*gamma) - 1._rkind/2._rkind*Ux(2)*(V_sq &
            *(gamma - 1) - 2*e*gamma + 2*v1**2*(gamma - 1)) - Ux(3)*v1*v2*( &
            gamma - 1) + Ux(4)*gamma*v1 + Uy(1)*v2*(V_sq*(gamma - 1) - e*gamma) - &
            Uy(2)*v1*v2*(gamma - 1) - 1._rkind/2._rkind*Uy(3)*(V_sq*(gamma - 1) - 2*e* &
            gamma + 2*v2**2*(gamma - 1)) + Uy(4)*gamma*v2
       do i = 1, nNode
          iNodeID = this%getNodeID(i)
          stabMat(:,iNodeID) = stabMat(:,iNodeID)                       &
               - integrator%getShapeFunc(k,i)*AiUi                       &
               * this%jacobianDet(k)*integrator%getWeight(k)/lumpedMass(iNodeID)
       end do
       !Tau & Nu:
       Tau    = 0._rkind
       drx    = Ux(1)
       dry    = Uy(1)
       dTx    = 0._rkind
       dTy    = 0._rkind
       dUx    = 0._rkind
       dUy    = 0._rkind
       h_rgne = 1.d-10
       h_rgn  = 1.d-10
       h_jgn  = 1.d-10
       dr2    = sqrt(drx*drx+dry*dry)+1.d-20
       rjx    = drx/dr2
       rjy    = dry/dr2
       do i = 1, nNode
          dTx  = dTx + T(i)*this%dNdx(i,k)
          dTy  = dTy + T(i)*this%dNdy(i,k)
       end do
       dT2  = sqrt(dTx*dTx+dTy*dTy)+1.d-20
       rTx  = dTx/dT2
       rTy  = dTy/dT2
       do i = 1, nNode
          dUx  = dUx + V_sq*this%dNdx(i,k)
          dUy  = dUy + V_sq*this%dNdy(i,k)
       end do
       dU2  = sqrt(dUx*dUx+dUy*dUy)+1.d-20
       rUx  = dUx/dU2
       rUy  = dUy/dU2
       smu  = 110._rkind
       fmu  = 0.41685_rkind*(temp/Tinf)**1.5_rkind*(Tinf+smu)/(temp+smu)
       do i = 1, nNode
          term_1 = abs(v1 *this%dNdx(i,k)+v1 *this%dNdy(i,k))
          term_2 = abs(rjx*this%dNdx(i,k)+rjy*this%dNdy(i,k))
          h_rgn1 = abs(rtx*this%dNdx(i,k)+rty*this%dNdy(i,k))
          h_rgn2 = abs(rux*this%dNdx(i,k)+ruy*this%dNdy(i,k))
          Tau    = Tau+term_1+term_2*Vc
          h_rgne = h_rgne+h_rgn1
          h_rgn  = h_rgn+h_rgn2
          h_jgn  = h_jgn+term_2
       end do
       Tau     = 1._rkind/Tau
       h_rgne  = 2._rkind/h_rgne
       h_rgn   = 2._rkind/h_rgn
       if (h_rgn .gt. 1.d3) h_rgn = 0._rkind
       h_jgn   = 2._rkind/h_jgn
       if (h_jgn .gt. 1.d10) h_jgn = 0._rkind
       tr1     = dr2*h_jgn/rho
       zzz     = h_jgn/2._rkind
       shoc    = (sqrt(tr1)+tr1**2)*Vc*2*zzz
       resumen = (1._rkind/Tau)**2._rkind +(2._rkind/dtmin)**2._rkind
       rrr     = resumen**(-0.5_rkind)
       t_sugn1 = rrr
       t_sugn2 = rrr
       t_sugn3 = rrr
       if(fmu.ne.0._rkind)then
          tau_sung3   = h_rgn**2._rkind/(4._rkind*fmu/rhoinf)
          tau_sung3_e = h_rgne**2._rkind/(4._rkind*fmu/rhoinf)
          t_sugn2     = (resumen+1._rkind/tau_sung3**2._rkind)**(-0.5_rkind)
          t_sugn3     = (resumen+1._rkind/tau_sung3_e**2._rkind)**(-0.5_rkind)
       end if
       this%Tau(1,k) = t_sugn1
       this%Tau(2,k) = t_sugn2
       this%Tau(3,k) = t_sugn3
       this%nu(k)    = shoc*cte
    end do
    deallocate(U, T)
  end subroutine calculateStabMat

  subroutine calculateResults(this, processInfo, resultMat)
    implicit none
    class(CFDElementDT)                       , intent(inout) :: this
    type(ProcessInfoDT)                       , intent(inout) :: processInfo
    real(rkind), dimension(:,:,:), allocatable, intent(inout) :: resultMat
    real(rkind), dimension(:,:)  , allocatable                :: U
    real(rkind)                                               :: rho, Vx, Vy, T
    real(rkind)                                               :: P, E, M, v2
    real(rkind)                                               :: R, Cv, Vc, gamma
    integer(ikind)                                            :: nNode, nDof, iNode
    integer(ikind)                                            :: i, j
    nDof = 4
    nNode = this%getnNode()
    allocate(resultMat(7,nNode,1), U(nDof,nNode))
    do i = 1, nNode
       do j = 1, nDof
          U(j,i) = this%node(i)%ptr%dof(j)%val
       end do
    end do
    R     = this%material%R
    Cv    = this%material%Cv
    Vc    = this%material%Vc
    gamma = this%material%gamma
    resultMat = 0._rkind
    do iNode = 1, nNode
       rho = U(1,iNode)
       Vx  = U(2,iNode)/rho
       Vy  = U(3,iNode)/rho
       E   = U(4,iNode)/rho
       v2  = Vx**2+Vy**2
       P   = rho*(gamma-1)*(E-.5_rkind*v2)
       T   = P/(rho*R)
       M   = sqrt(v2/(T*gamma*R))
       resultMat(1,iNode,1) = rho
       resultMat(2,iNode,1) = Vx
       resultMat(3,iNode,1) = Vy
       resultMat(4,iNode,1) = E
       resultMat(5,iNode,1) = T
       resultMat(6,iNode,1) = M
       resultMat(7,iNode,1) = P
    end do
  end subroutine calculateResults

  subroutine calculateLHS(this, processInfo, lhs)
    implicit none
    class(CFDElementDT) , intent(inout) :: this
    type(ProcessInfoDT) , intent(inout) :: processInfo
    type(LeftHandSideDT), intent(inout) :: lhs
  end subroutine calculateLHS

  subroutine calculateRHS(this, processInfo, rhs)
    implicit none
    class(CFDElementDT)                             , intent(inout) :: this
    type(ProcessInfoDT)                             , intent(inout) :: processInfo
    real(rkind)        , dimension(:)  , allocatable, intent(inout) :: rhs
  end subroutine calculateRHS

  subroutine setupIntegration(this, integrator, valuedSource)
    implicit none
    class(CFDElementDT)                                 , intent(inout) :: this
    type(IntegratorPtrDT)                               , intent(in)    :: integrator
    real(rkind), dimension(4,integrator%getIntegTerms()), intent(out)   :: valuedSource
    integer(ikind)                                                      :: i, nNode
    type(NodePtrDT), dimension(:), allocatable                          :: nodalPoints
    nNode = this%getnNode()
    valuedSource = this%getValuedSource(integrator)
  end subroutine setupIntegration

  function getValuedSource(this, integrator)
    implicit none
    class(CFDElementDT), intent(inout) :: this
    type(IntegratorPtrDT) , intent(in) :: integrator
    real(rkind), dimension(4,integrator%getIntegTerms()) :: getValuedSource
    integer(ikind) :: i, j, nNode
    real(rkind) :: x, y
    type(NodePtrDT), dimension(:), allocatable :: node
    nNode = this%getnNode()
    do i = 1, integrator%getIntegTerms()
       node = this%node
       x = 0
       y = 0
       do j = 1, nNode
          x = x + integrator%getShapeFunc(i,j)*node(j)%ptr%getx()
          y = y + integrator%getShapeFunc(i,j)*node(j)%ptr%gety()
       end do
       getValuedSource(1,i) = this%source(1)%evaluate(1, (/x,y/))
       getValuedSource(2,i) = this%source(1)%evaluate(2, (/x,y/))
       getValuedSource(3,i) = this%source(1)%evaluate(3, (/x,y/))
       getValuedSource(4,i) = this%source(1)%evaluate(4, (/x,y/))
    end do
  end function getValuedSource

  subroutine vectorCalculus(this)
    implicit none
    class(CFDElementDT), intent(inout) :: this
    type(IntegratorPtrDT)                          :: integrator
    type(NodePtrDT), dimension(:)    , allocatable :: nodalPoints
    real(rkind)    , dimension(:,:,:), allocatable :: jacobian
    real(rkind)    , dimension(:)    , allocatable :: jacobianDet
    real(rkind)                                    :: bi, ci
    integer(ikind)                                 :: i, j, nNode, nDof
    nNode = this%getnNode()
    nDof  = this%node(1)%getnDof()
    integrator    = this%getIntegrator()
    allocate(this%dNdx(nNode, integrator%getIntegTerms()))
    allocate(this%dNdy(nNode, integrator%getIntegTerms()))
    allocate(this%jacobianDet(integrator%getIntegTerms()))
    allocate(this%tau(3, integrator%getIntegTerms()))
    allocate(this%nu(integrator%getIntegTerms()))
    jacobian      = this%geometry%jacobianAtGPoints(this%node)
    jacobianDet   = this%geometry%jacobianDetAtGPoints(jacobian)
    do i = 1, nNode
       do j = 1, integrator%getIntegTerms()
          bi = jacobian(j,2,2)*integrator%getDShapeFunc(j,1,i) &
               - jacobian(j,1,2)*integrator%getDShapeFunc(j,2,i)
          ci = jacobian(j,1,1)*integrator%getDShapeFunc(j,2,i) &
               - jacobian(j,2,1)*integrator%getDShapeFunc(j,1,i)
          this%dNdx(i,j) = bi/jacobianDet(j)
          this%dNdy(i,j) = ci/jacobianDet(j)
       end do
    end do
    do j = 1, integrator%getIntegTerms()
       this%jacobianDet(j) = jacobianDet(j)
    end do

  end subroutine vectorCalculus

end module CFDElementM


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X     Learn-on-the-fly (LOTF) hybrid molecular dynamics code
!X
!X
!X     Authors: Gabor Csanyi, Alessio Commisso, Steven Winfield
!X     James Kermode, Gianpietro Moras, Michael Payne, Alessandro De Vita
!X
!X     Copyright 2005, All Rights Reserved
!X
!X     This source code is confidential, all distribution is
!X     prohibited. Making unauthorized copies is also prohibited
!X
!X     When using this software, the following should be referenced:
!X
!X     Gabor Csanyi, Tristan Albaret, Mike C. Payne and Alessandro De Vita
!X     "Learn on the fly": a hybrid classical and quantum-mechanical
!X         molecular dynamics simulation
!X     Physical Review Letters 93 p. 175503 (2004) >>PDF [626 KB]
!X
!X     Gabor Csanyi, T. Albaret, G. Moras, M. C. Payne, A. De Vita
!X     Multiscale hybrid simulation methods for material systems
!X     J. Phys. Cond. Mat. 17 R691-R703 Topical Review (2005)
!X
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X  AdjustablePotential_linear module
!X
!X  a set of linear potentials that can be adapted a reproduce
!X  given forces on atoms
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

module AdjustablePotential_module
  use system_module, only : dp, print, inoutput, PRINT_VERBOSE, PRINT_NERD, line, optional_default, mainlog, value, system_abort, reallocate, print_warning, operator(//), mpi_id, abort_on_mpi_error
  use linearalgebra_module
  use table_module
  use sparse_module
  use minimization_module

  use atoms_types_module
  use atoms_module
  use clusters_module

#ifdef _MPI
#ifndef _OLDMPI
  use mpi
#endif
#endif
  implicit none
#ifdef _MPI
#ifdef _OLDMPI
include 'mpif.h'
#endif
#endif
  private
  SAVE

  public :: adjustable_potential_init, adjustable_potential_optimise, adjustable_potential_force, adjustable_potential_finalise

  ! ratio of smallest to largest eigenvalue of spring space spanning matrix
  ! for an atom to be considered not well described
  real(dp), parameter :: AP_bad_atom_threshold = 0.1_dp

  ! Threshold for spring constants above which a warning will be issued
  real(dp), parameter :: adjustable_potential_max_spring_constant = 1.0_dp

  integer, parameter :: adjustable_potential_nparams = 1

  ! are we initialised or not
  logical::adjustable_potential_initialised = .false.

  ! have we got some parameters saved
  logical::adjustable_potential_parameters_saved = .false.

  ! this table contains the pairs of atoms for which we have a two-body
  ! adjustable spring potential, together with their distance
  ! and the potential parameters. the parameters are just the y values of the spring
  ! !!!the integer indices refer to the order of atoms in the fitlist
  ! argument of the init() routine!!!
  !
  !                    int               real
  !
  !fitlist indexes ->  fi,fj,shift       rij, fij <- fij: variable parameter
  type(table)::twobody

  ! this table is used as a store for old data which can be used to initialise
  ! parameters to values from the last interpolation step. It also stores
  ! parameters from one end point during interpolation
  type(table)::twobody_old

  real(dp), dimension(:,:), allocatable:: fitforce_old

  !The exclusion list table is used to prevent springs being created between
  !certain pairs of atoms, e.g. if there is a bond length constraint between
  !them.
  type(table) :: exclusion_list
  logical     :: exclusions = .false.

  ! Operation                                            twobody   twobody_old
  !
  ! Extrapolate from t1 -> t2 using parameters P1           P1  --,     --
  ! Compute parameters P2' at t2                            P2'   `-->  P1
  ! Interpolate from t1 -> t2 using P1 and P2'              P2' --,     P1
  ! Call CalcConnect                                        --    `-->  P2'
  ! Create parameters P2 using P2' as initial values        P2          P2'
  ! Extrapolate from t2 -> t3 using parameters P2           P2          --
  ! ...

  ! local copy of target force to fit to: F1x,F1yF1z,F2x,F2y,F2z...
  real(dp), allocatable, dimension(:)   ::target_force

  ! local copy of fitlist indices
  ! atomlist%int(1,fi)   = i,    fi: index in fitlist, i: index in Atoms structure
  type(Table) :: atomlist, atomlist_old
!!$  type(Table) :: fitlist_old

  ! forcematrix multiplied onto the vector of
  ! parameters (fij...)
  ! gives the force vector on the atoms. the number of columns
  ! is the number of parameters, the number of rows is 3*the number of atoms
  !
  ! force = forcematrix * params
  !
  ! ...0..y1_x y2_x...0......
  ! ...0..y1_y y2_y...0......
  ! ...0..y1_z y2_z...0......
  type(sparse)                          ::forcematrix ! linear part

contains

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X Initialize adjustable potential
  !X
  !X inputs are an atoms structure and a table with a single list integers
  !X of atoms for which target forces are to be fitted.
  !X
  !X The adjustable potential will be fitted to the difference between
  !X a quantum mechanical force and a classical force.
  !X
  !X An adjustable potential will be created between every pair of neighbouring
  !X atoms that are given in the fitlist.
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !

  ! this is called with a table
  subroutine adjustable_potential_init(atoms_in, fitlist, directionN, spring_hops, &
       map, method, nnonly)
    type(Atoms),       intent(in) :: atoms_in
    type(Table),       intent(in) :: fitlist
    integer,           intent(in) :: directionN
    integer, optional, intent(in) :: spring_hops
    logical, optional, intent(in) :: map, nnonly
    character(len=*),optional,intent(in):: method

    real(dp):: ratio, evals(3), evecs(3,3), small_evec(3), r_ij, u_ij(3), cos_theta, oldratio
    real(dp), allocatable, dimension(:)::default_row
    real(dp), allocatable, dimension(:,:) :: df
    integer::Nfit_components, i, j, fi, fj, nj, n, hops, atomcount, oldn, &
         shift(3), k, best_j, best_fj, best_shift(3), sloc(1), sx, sy, sz
    type(Table) :: springs, newsprings, tmpsprings,cand_springs
    logical :: do_map, do_nnonly
    logical :: do_refit

    do_map = .false.
    if (present(map)) do_map = map

    do_nnonly = .true.
    if (present(nnonly)) do_nnonly = nnonly

    hops = 2
    if(present(spring_hops)) then
       hops = spring_hops
    end if

    ! If QM region hasn't changed since we last made out spring list, then
    ! there's nothing to be done now.

    if(do_map) then
       ! this save makes the twobody_old equal to twobody
       call adjustable_potential_save_parameters
       call print("Saving old parameters ....", PRINT_VERBOSE)
    end if

!!$    if (first_time) then
!!$       first_time = .false.
!!$    else if (fitlist%N == fitlist_old%N) then
!!$
!!$       allocate(s1(fitlist%N),s2(fitlist%N))
!!$       s1 = fitlist%int(1,1:fitlist%N)
!!$       s2 = atomlist%int(1,1:atomlist%N)
!!$       call sort_array(s1)
!!$       call sort_array(s2)
!!$
!!$       if (all(int_part(fitlist) == int_part(fitlist_old))) then
!!$          call print('fitlist unchanged, nothing to be done in adjustable_potential_init', PRINT_VERBOSE)
!!$          deallocate(s1)
!!$          deallocate(s2)

          ! Atom set matches, but it's possible some of the shifts have changed
          ! if any fit atoms have crossed a periodic boundary since last call to
          ! adjustable_potential_init. Let's overwrite with new shifts.
!!$          do n=1,twobody%N
!!$             i = atomlist%int(1,twobody%int(1,n))
!!$             j = atomlist%int(1,twobody%int(2,n))
!!$             ! overwrite old shift in twobody%int(3:5,n) with correct shift
!!$             r = distance_min_image(atoms_in, i, j, shift=shift)
!!$
!!$             if (abs(r-twobody%real(1,i)) < 0.1_dp) then
!!$                call print('shift update: '//shift//' '//twobody%int(3:5,i)//' '//r//' '//twobody%real(1,i))
!!$                twobody%int(3:5,i) = shift
!!$             end if
!!$
!!$             !! for the thin crack system, not all springs are min
!!$             !! image springs; some of them span periodic boundary.
!!$             !! how can we get round this?
!!$
!!$          end do
!!$
!!$          return
!!$       end if
!!$
!!$       deallocate(s1)
!!$       deallocate(s2)
!!$    end if
!!$    call wipe(fitlist_old)
!!$    call append(fitlist_old, fitlist)

    call print("Setting up adjustable potential....", PRINT_VERBOSE)

    ! test table size
    if(fitlist%intsize /= 4 .and. fitlist%intsize /= 1)  &
         call system_abort("adjustable_potential_init: fitlist%intsize must be 1 or 4")
    if(fitlist%realsize /= 0) &
         call system_abort("adjustable_potential_init: fitlist%realsize /= 0")

    if (multiple_images(fitlist)) &
         call system_abort('adjustable_potential_init: fitlist contains repeated atom indices')

    call print("Fitting on atoms:", PRINT_VERBOSE)
    call print(fitlist, PRINT_VERBOSE)
    call print("", PRINT_VERBOSE)

    ! allocate globals
    call print("Allocating target_force", PRINT_NERD)
    Nfit_components = fitlist%N*3

    call reallocate(target_force,Nfit_components)
    target_force = 0.0_dp

    ! copy atom indices from argument to global

    call wipe(atomlist)
    call allocate(atomlist,1,0,0,0,fitlist%N)
    call append(atomlist, fitlist%int(1,1:fitlist%N))

    ! set up twobody table that holds pairs of atom indices and shiftts
    ! and the associated potential parameters (initialized to 0)

    ! wipe previous table
    call print("Wiping twobody table", PRINT_NERD)
    call wipe(twobody)
    call allocate(twobody, 5, 1+2*adjustable_potential_nparams,0,0)

    call allocate(springs, 4, 0, 0, 0)
    call allocate(newsprings, 4, 0, 0, 0)
    call allocate(tmpsprings, 4, 0, 0, 0)

    ! default parameters
    allocate(default_row(1+2*adjustable_potential_nparams))
    default_row = 0.0_dp

    ! First we try to add springs between pairs of atoms in fitlist within 'hops'
    ! bond hops of one another
    call print('Adding default springs')
    do fi=1,atomlist%N ! loop through fit atoms
       i = atomlist%int(1,fi)

       call wipe(newsprings)
       call append(newsprings, (/i,0,0,0/))
       call bfs_grow(atoms_in, newsprings, hops, nneighb_only = do_nnonly)

       do nj = 1, newsprings%N

          j = newsprings%int(1,nj)
          shift = newsprings%int(2:4,nj)

          fj = find(atomlist, j)
          if (fj == 0) cycle ! not in fitlist
          if (i == j) cycle  ! don't join to self or images of self

          ! Have we already got this spring?
          ! always store in order fi, fj with fi < fj
          if (find(twobody, (/min(fi,fj), max(fi,fj), sign(1,fj-fi)*shift/)) /= 0) cycle

          ! Is this pair excluded?
          if (exclusions) then
             if (is_excluded(i,j,shift)) cycle
          end if

          default_row(1) = distance(atoms_in, i, j, shift)
          call append(twobody, (/min(fi,fj), max(fi,fj), sign(1,fj-fi)*shift/), default_row)
       end do

    end do
    write(line, '(a,i0,a)') 'Got ', twobody%N, ' springs' ; call print(line)
    write(line, '(a,f0.2)') 'Searching for difficult atoms, anisotropy threshold: ', AP_bad_atom_threshold
    call print(line, PRINT_VERBOSE)

    oldn = twobody%N
    atomcount = 0
    ! For those fit atoms that we require directionality for, check how well the
    ! springs connected to it span 3D space. If necessary add more springs for those atoms.
    do fi=1,directionN
       i = atomlist%int(1, fi)

       ! Find who atom i is connected to with springs, and get associated shifts
       call wipe(springs)
       do j = 1, twobody%N
          if (twobody%int(1,j) == fi) &
               call append(springs, (/atomlist%int(1,twobody%int(2,j)), twobody%int(3:5,j)/))
          if (twobody%int(2,j) == fi) &
               call append(springs, (/atomlist%int(1,twobody%int(1,j)),-twobody%int(3:5,j)/))
       end do

       ! How well do these springs span 3D space?
       if (springs%N > 2) then
          call directionality(atoms_in, i, springs, evals, evecs, method=2)
          ratio = minval(evals)/maxval(evals)
       else
          ratio = 0.0_dp ! less than three springs - really bad!
       end if

       sloc = minloc(evals)
       small_evec = evecs(:,sloc(1))

       if (ratio < AP_bad_atom_threshold) then

          write (line,'(a,i6,a,i6,a,f10.3)') 'Atom ', i, ' with ', springs%N, &
               ' springs has ratio ', ratio
          call print(line, PRINT_VERBOSE)
          call print('Smallest evect ='//small_evec, PRINT_VERBOSE)

          ! Not well enough
          atomcount = atomcount + 1

          ! Tabulate candidate springs
          call allocate(cand_springs, 5, 2, 0, 0)
          do fj = 1,atomlist%N
             j = atomlist%int(1,fj)

             if (i == j) cycle ! Don't join to self

             do sx=-1,1
                do sy=-1,1
                   do sz=-1,1

                      shift = (/sx,sy,sz/)

                      if (find(springs, (/j,shift/)) /= 0) cycle ! We've tried this spring already

                      ! Don't add if this pair is excluded
                      if (exclusions) then
                         if (is_excluded(i,j,shift)) cycle
                      end if

                      r_ij = distance(atoms_in,i,j,shift)
                      u_ij = diff(atoms_in,i,j,shift)
                      cos_theta = abs(small_evec .dot. u_ij) / (norm(small_evec)*norm(u_ij))

                      call append(cand_springs, (/j,fj,shift/), (/r_ij,cos_theta/))

                   end do
                end do
             end do
          end do

          ! Try to add springs until ratio exceeds thresholds
          n = 0
          do while (ratio < AP_bad_atom_threshold .and. n < 8)

             ! Find candidate spring with largest cos_theta
             sloc = maxloc(cand_springs%real(2,1:cand_springs%N))
             k = sloc(1)

             best_j     = cand_springs%int(1,k)
             best_fj    = cand_springs%int(2,k)
             best_shift = cand_springs%int(3:5,k)
             r_ij       = cand_springs%real(1,k)

             call print('Trying spring to atom '//best_j//' (index '//best_fj//') shift '//best_shift//' length '//r_ij)

             call delete(cand_springs, k)

             call append(springs, (/best_j,best_shift/))

             oldratio = ratio

             call directionality(atoms_in, i, springs, evals, evecs, 2)
             ratio = minval(evals)/maxval(evals)

             ! Add this spring to twobody only if it improves ratio
             if (ratio > oldratio) then
                default_row(1) = r_ij
                call append(twobody,(/min(fi,best_fj),max(fi,best_fj),sign(1,best_fj-fi)*best_shift/), default_row)

                write (line, '(a,i0,a,f0.3)') '  Now got ', springs%N,' springs. Ratio improved to ', ratio
                call print(line, PRINT_VERBOSE)
             end if

             n = n + 1

!!$     newsprings = springs
!!$
!!$          call bfs_step(atoms_in, newsprings, tmpsprings, nneighb_only = .true.)
!!$          call append(newsprings, tmpsprings)
!!$
!!$          if (tmpsprings%N == 0) then
!!$             ! Didn't manage to add any atoms, so try nneigb_only = false, since
!!$             ! these are the 'bad' atoms and their neighbours are
!!$             ! probably some distance away
!!$
!!$             call bfs_step(atoms_in, newsprings, tmpsprings, nneighb_only = .false.)
!!$             call append(newsprings, tmpsprings)
!!$          end if
!!$
!!$          do nj = 1, newsprings%N
!!$             j = newsprings%int(1,nj)
!!$             shift = newsprings%int(2:4,nj)
!!$             fj = find(atomlist, j)
!!$             if (fj == 0) cycle ! not in fitlist
!!$             if (i == j) cycle ! don't join to self or images of self
!!$
!!$             ! already got this spring
!!$             if (find(twobody, (/min(fi,fj), max(fi,fj), sign(1,fj-fi)*shift/)) /= 0) cycle
!!$
!!$             ! Don't add if this pair is excluded
!!$             if (exclusions) then
!!$                if (is_excluded(i,j,shift)) cycle
!!$             end if
!!$
!!$             default_row(1) = distance(atoms_in, i, j, shift)
!!$             call append(twobody,(/min(fi,fj),max(fi,fj),sign(1,fj-fi)*shift/), default_row)
!!$          end do
!!$
!!$          ! Work out the springs again
!!$          call wipe(springs)
!!$          do j = 1, twobody%N
!!$             if (twobody%int(1,j) == fi) &
!!$                  call append(springs, (/atomlist%int(1,twobody%int(2,j)), twobody%int(3:5,j)/))
!!$             if (twobody%int(2,j) == fi) &
!!$                  call append(springs, (/atomlist%int(1,twobody%int(1,j)),-twobody%int(3:5,j)/))
!!$          end do
!!$
!!$          ! Did we fail to add any springs at all? If so system has probably exploded
!!$          if (springs%N == 0) then
!!$             write(line, '(a,i0)') 'adjustable_potential_init: can''t add any springs for atom ', i
!!$             call system_abort(line)
!!$          end if
!!$
!!$          call directionality(atoms_in, i, springs, evals, evecs, 2)
!!$          ratio = minval(evals)/maxval(evals)
!!$
!!$          write (line, '(a,i0,a,f0.3)') '  Now got ', springs%N,' springs. Ratio improved to ', ratio
!!$          call print(line, PRINT_VERBOSE)
!!$
!!$          n = n + 1
          end do
          if (n == 8 .and. ratio < AP_bad_atom_threshold) then
             call system_abort('adjustable_potential_init: failed to sufficiently improve bad spring space spanning')
          end if

       end if

    end do

    if(atomcount > 0) then
       write(line, '(a,i0,a,i0,a)') 'Added ', twobody%N-oldn, ' springs to fix ', atomcount, ' atoms'
       call print(line)
    end if

    call finalise(springs)
    call finalise(newsprings)
    call finalise(tmpsprings)
    deallocate(default_row)

    call print("Allocating twobody table to proper size", PRINT_NERD)
    call allocate(twobody) ! reduce table to proper length

    write(line, '(a,i0)') "Number of force components: ", size(target_force) ; call print(line)
    write(line, '(a,i0)') "Number of parameters:       ", twobody%N*adjustable_potential_nparams ; call print(line)

    adjustable_potential_initialised = .true.

    if(do_map) then
       ! this maps from the twobody_old with the old structure into twobody with the new structure
       call adjustable_potential_map_parameters

       ! Only need to refit if spline set has changed since last time
       do_refit = .false.
       ! if sizes don't match definitely need to refit
       if (twobody%N /= twobody_old%N) then
          do_refit = .true.
       else
          do i=1,twobody%N
             ! As soon as we can't find a new spring in twobody_old, know we have to refit
             if (find(twobody_old,twobody%int(:,i)) == 0) then
                do_refit = .true.
                exit
             end if
          end do
       end if

       ! reoptimise with new springs to reproduce old forces at new positions
       ! if fit region has changed, just copy the target force components that
       ! are in both the old and new regions
       if (do_refit) then
          call print('Springs changed, refitting old forces with new spring set')
          allocate(df(3,atomlist%N))
          df = 0.0_dp
          do i=1,atomlist_old%N
             fi = find(atomlist, atomlist_old%int(:,i)) ! find old fit atom in new fitlist
             if (fi /= 0) df(:,fi) = fitforce_old(:,i)         ! if found, copy target force
          end do
          call adjustable_potential_optimise(atoms_in, df, method=method)
          deallocate(df)
       end if

    endif

    ! this save makes the twobody_old equal to twobody (both new structure)
    call adjustable_potential_save_parameters

  end subroutine adjustable_potential_init

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_finalise()
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !
  subroutine adjustable_potential_finalise()
    if(allocated(target_force)) deallocate(target_force)
    call finalise(atomlist)
    call finalise(atomlist_old)
    call finalise(twobody)
    call finalise(twobody_old)
    call adjustable_potential_delete_saved_parameters
    call sparse_finalise(forcematrix)
    if (exclusions) call finalise(exclusion_list)
  end subroutine adjustable_potential_finalise

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_delete_saved_parameters
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_delete_saved_parameters
    call finalise(twobody_old)
    call finalise(atomlist_old)
    adjustable_potential_parameters_saved = .false.
  end subroutine adjustable_potential_delete_saved_parameters

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_save_parameters
  !X
  !X Makes a copy of the twobody and atomlist data, avoiding a
  !X reallocation if possible
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_save_parameters
    if (.not.adjustable_potential_initialised) &
         call system_abort('Adjustable_Potential_Save_Parameters: No parameters exist to save')
    call wipe(twobody_old)
    call append(twobody_old,twobody)

    call wipe(atomlist_old)
    call append(atomlist_old, atomlist)

    adjustable_potential_parameters_saved = .true.

  end subroutine adjustable_potential_save_parameters

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_map_parameters
  !X
  !X Given an old set of data (twobody_old, atomlist_old), fill in the
  !X y1 and y2 values which were used previously
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_map_parameters
    integer :: i, j, fi, fj, n, m, k, shift(3)

    if(.not.adjustable_potential_initialised) &
         call system_abort('Adjustable_Potential_Map_Parameters: No parameters to map!')
    if (twobody_old%N == 0) &
         call system_abort('Adjustable_Potential_Map_Parameters: No parameters are saved')

    if (twobody%N >= twobody_old%N) then  ! loop over the smallest table

       !Loop over the entries in the old table
       do n = 1, twobody_old%N

          !look up the atom indices
          i     = atomlist_old%int(1,twobody_old%int(1,n))
          j     = atomlist_old%int(1,twobody_old%int(2,n))
          shift = twobody_old%int(3:5,n)

          !try to find the corresponding atoms in the new list
          fi = find(atomlist,i)
          if (fi == 0) cycle !atom i is not present in the new list
          fj = find(atomlist,j)
          if (fj == 0) cycle !atom j is not present in the new list

          !find the entry in the new twobody table
          m = find(twobody,(/min(fi,fj),max(fi,fj),sign(1,fj-fi)*shift/))
          if (m == 0) cycle ! i and j are no longer connected

          ! copy the parameter values across
          do k=1,2*adjustable_potential_nparams
             twobody%real(1+k,m) = twobody_old%real(1+k,n)
          end do

       end do

    else

       !Loop over the entries in the new table
       do n = 1, twobody%N

          !look up the atom indices
          i     = atomlist%int(1,twobody%int(1,n))
          j     = atomlist%int(1,twobody%int(2,n))
          shift = twobody%int(3:5,n)

          !try to find the corresponding atoms in the old list
          fi = find(atomlist_old, i)
          if (fi == 0) cycle !atom i is not present in the old list
          fj = find(atomlist_old, j)
          if (fj == 0) cycle !atom j is not present in the old list

          !find the entry in the old twobody table
          m = find(twobody_old,(/min(fi,fj),max(fi,fj),sign(1,fj-fi)*shift/))
          if (m == 0) cycle ! i and j were not neighbours

          !copy the param values across
          do k=1,2*adjustable_potential_nparams
             twobody%real(1+k,n) = twobody_old%real(1+k,m)
          end do

       end do

    end if

  end subroutine adjustable_potential_map_parameters

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_force(atoms, force)
  !X
  !X force and energy from the adjustable potential, using the interatomic
  !X distance information from atoms. If interp is given, then interpolate
  !X the parameters between twobody_old and twobody
  !X
  !X the force(:) array is the same size as the stored atomlist.
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_force(atoms_in, force, interp, interp_space, interp_order, energy, power)
    type(Atoms),                            intent(in)  :: atoms_in
    real(dp), dimension(3,atomlist%N), intent(out) :: force
    real(dp), optional,                     intent(in)  :: interp
    logical, optional,                      intent(in)  :: interp_space
    character(len=*),optional,              intent(in)  :: interp_order
    real(dp), optional, intent(out) :: energy, power

    ! local
    integer:: i, a1, a2, fi, fj, shift(3), my_interp_order
    real(dp), dimension(3)::f
    real(dp)::fij, r, interp_point, dalpha
    logical::interpolate, do_interp_space

    do_interp_space = optional_default(.false., interp_space)

    interpolate = present(interp) .or. do_interp_space
    if(interpolate .and. .not.adjustable_potential_parameters_saved) &
         call System_Abort('Adjustable_Potential_Force: No parameters saved to interpolate with')

    my_interp_order = 1
    if(present(interp_order)) then
       if(trim(interp_order) .eq. 'linear') then
          my_interp_order = 1
       else if(trim(interp_order) .eq. 'quadratic') then
          my_interp_order = 2
       else if(trim(interp_order) .eq. 'cubic') then
          my_interp_order = 3
       else
          call System_Abort("Invalid interpolation order `"//interp_order//"'in Adjustable_Potential_Force")
       end if
    end if

    force = 0.0_dp

    if(present(energy)) energy = 0.0_dp
    if(present(power))   power = 0.0_dp

    do i=1,twobody%N

       fi = twobody%int(1,i)
       fj = twobody%int(2,i)
       shift = twobody%int(3:5,i)

       a1 = atomlist%int(1,fi)
       a2 = atomlist%int(1,fj)
       r = distance(atoms_in, a1, a2, shift)

       if(interpolate) then
          if(do_interp_space) then
             interp_point = (r - twobody_old%real(1,i))/(twobody%real(1,i)-twobody_old%real(1,i))
             fij = linear_interpolate(0.0_dp,twobody_old%real(2,i),&   ! x0,y0
               1.0_dp,twobody%real(2,i),interp_point) !x1, y1, x
          else
             select case(my_interp_order)
                case(1) ! linear interpolation
                   fij = linear_interpolate(0.0_dp,twobody_old%real(2,i),&   ! x0,y0
                        1.0_dp,twobody%real(2,i),interp) !x1, y1, x

                   dalpha = (twobody%real(2,i) - twobody_old%real(2,i))
                case(2) ! quadratic using the derivative info at start point
                   fij = twobody_old%real(2,i)+interp*twobody_old%real(3,i)&
                        +interp*interp*(twobody%real(2,i)-twobody_old%real(2,i)-twobody_old%real(3,i))

                   dalpha = 0.0_dp
                case(3) ! cubic with zero derivatives at endpoints
                   fij = cubic_interpolate(0.0_dp,twobody_old%real(2,i),&   ! x0,y0
                        1.0_dp,twobody%real(2,i),interp) ! x1,y1, x

                   dalpha = 0.0_dp
                case default
                   call System_Abort("Invalid interpolation order in Adjustable_Potential_Force")
             end select
          end if
       else
          ! Extrapolate
          fij = twobody%real(2,i)
          dalpha = 0.0_dp ! Parameters don't change during extrap
       end if

       ! optionally calculate energy and power
       if(present(energy)) energy = energy + r*fij
       if(present(power))  power  = power  + r*dalpha ! (since dV_dalpha = r, and dt = 1 here)

       if (value(mainlog%verbosity_stack) >= PRINT_NERD) &
            call print(i//' '//a1//' '//a2//' '//r//' '//fij//' <-- SPRING', PRINT_NERD)

       ! get force
       f = direction_cosines(atoms_in,a1,a2,shift)*fij
       force(:,fi) = force(:,fi) + f
       force(:,fj) = force(:,fj) - f
    end do

  end subroutine adjustable_potential_force

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_force_sparse(params)
  !X
  !X force from the adjustable potential, using the precomputed
  !X sparse matrix
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !
  function adjustable_potential_force_sparse(params) result(force)
    real(dp)::params(:)
    real(dp)::force(size(target_force))

    if(.not.adjustable_potential_initialised) &
         call system_abort("adjustable_potential_force_sparse: not initialised!")

    force = (forcematrix .mult. params)
  end function adjustable_potential_force_sparse

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_force_error(params)
  !X
  !X normsq() of difference between target and adjustable force
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !
  function adjustable_potential_force_error(params,data) result(error)
    real(dp)::params(:)
    character,optional :: data(:)
    real(dp)::error

    if(.not.adjustable_potential_initialised) &
         call system_abort("adjustable_potential_force_error: not initialised!")

    if(value(mainlog%verbosity_stack) .ge. PRINT_NERD) then
       write(line, *) "adjustable_potential_force_error: params(", size(params), ")"
       call print(line, PRINT_NERD)
       call print(params, PRINT_NERD)
    end if

    error = normsq(target_force - (forcematrix .mult. params))

  end function adjustable_potential_force_error

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_force_error_deriv(params)
  !X
  !X derivative of the force_error wrt variable parameters
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !
  function adjustable_potential_force_error_deriv(params,data) result(deriv)
    real(dp)::params(:)
    character,optional :: data(:)
    real(dp)::force(size(target_force))
    real(dp)::deriv(size(params))

    if(.not.adjustable_potential_initialised) &
         call system_abort("adjustable_potential_force_error_deriv: not initialised!")

    if(value(mainlog%verbosity_stack) .ge. PRINT_NERD) then
       write(line, *) "adjustable_potential_force_error_deriv: params(", size(params), ")"
       call print(line, PRINT_NERD)
       call print(params, PRINT_NERD)
    end if

    force = (forcematrix .mult. params)
    deriv = (-2.0_dp * (target_force - force)) .mult. forcematrix
  end function adjustable_potential_force_error_deriv

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_create_forcematrix()
  !X
  !X creates the sparse matrix forcematrix that we
  !X use to compute the adjustable forces quickly for different
  !X parameters, but unchanging atomic positions
  !X
  !X forcematrix has as many rows as the 3*number of atoms to be fitted
  !X on, so same as size(target_force), and as many columns as the
  !X number of adjustable parameters, so twobody%N*adjustable_potential_nparams
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !

  subroutine adjustable_potential_create_forcematrix(atoms_in, precondition)
    type(Atoms), intent(IN)::atoms_in
    real(dp), intent(in) :: precondition(:)

    ! temporary hack to go via a full matrix. to be redone
    real(dp), allocatable, dimension(:,:)::tmpmatrix
    integer, parameter::nparams = adjustable_potential_nparams
    integer::i, j, k, a1, a2, fi, fj, shift(3)
    real(dp)::cosi(3)
    integer, allocatable::colcounter(:), tmpcolumns(:,:)

    write(line, *) "Setting up adjustable potential forcematrix.." ; call print(line, PRINT_VERBOSE)
    write(line, *) ; call print(line, PRINT_VERBOSE)

    allocate(tmpmatrix(size(target_force), twobody%N*nparams))
    allocate(colcounter(size(target_force)))
    allocate(tmpcolumns(size(target_force), twobody%N*nparams))
    tmpmatrix = 0.0_dp
    colcounter = 0
    tmpcolumns = 0

    ! loop through the springs
    do i=1,twobody%N
       ! we take and cache the distance from atoms_in, as it might have changed since time
       ! we initialised the twobody table
       fi = twobody%int(1,i)
       fj = twobody%int(2,i)
       shift = twobody%int(3:5,i)
       a1 = atomlist%int(1,fi)
       a2 = atomlist%int(1,fj)

       twobody%real(1,i) = distance(atoms_in, a1, a2, shift)
       cosi = direction_cosines(atoms_in, a1, a2, shift)

       fi = (twobody%int(1,i)-1)  ! subtract 1 so we can use it as index
       do k=1,3
          colcounter(fi*3+k) = colcounter(fi*3+k) + 1
          ! Multiply by preconditioning factor for this spring
          tmpmatrix (fi*3+k, colcounter(fi*3+k)) = cosi(k)*precondition(i)
          tmpcolumns(fi*3+k, colcounter(fi*3+k)) = (i-1)*nparams+1
       end do
       fj = (twobody%int(2,i)-1) ! subtract 1 so we can use it as index
       do k=1,3
          colcounter(fj*3+k) = colcounter(fj*3+k) + 1
          ! Multiply by preconditioning factor for this spring
          tmpmatrix (fj*3+k, colcounter(fj*3+k)) = -cosi(k)*precondition(i)
          tmpcolumns(fj*3+k, colcounter(fj*3+k)) = (i-1)*nparams+1
       end do
    end do ! i=1,twobody%N

    ! fill in sparse matrix from values in tmpmatrix and tmpcolumns
    call sparse_init(forcematrix, size(target_force), twobody%N*nparams)
    do i=1,size(target_force)
       forcematrix%rows(i) = forcematrix%table%N+1
       do j=1,colcounter(i)
          call append(forcematrix%table, tmpcolumns(i,j), tmpmatrix(i,j))
       end do
    end do
    ! put in last row value
    forcematrix%rows(i) = forcematrix%table%N+1

    deallocate(tmpmatrix)
    deallocate(tmpcolumns)
    deallocate(colcounter)

    ! print it
    if(value(mainlog%verbosity_stack) .ge. PRINT_NERD) then
       call print('Forcematrix:',PRINT_NERD)
       call print(forcematrix, PRINT_NERD)
       call print('',PRINT_NERD)
    end if
  end subroutine adjustable_potential_create_forcematrix

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_optimise()
  !X
  !X the whole point of this module: vary the adjustable parameters
  !X to minimize the force_error()
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_optimise(atoms_in, fitforce, fromscratch, method, prec_func, parallel_svd)
    type(Atoms),              intent(in):: atoms_in
    real(dp), dimension(:,:), intent(in):: fitforce
    logical, optional                   :: fromscratch

    logical                             :: do_fromscratch
    character(len=*),optional,intent(in):: method
    logical, optional, intent(in)       :: parallel_svd
    logical::do_svd, do_parallel_svd
    integer, parameter::nparam = adjustable_potential_nparams
    integer::nparamtot, bad_spring
    real(dp), allocatable::params(:)
    integer::steps, i, shift(3), k
    real(dp):: err, ferr(3, atomlist%N)
    real(dp), parameter::eigenvalue_threshold=1e-10_dp
    real(dp), allocatable::full_forcematrix(:,:), X(:,:), S(:), WORK(:), precondition(:)
    integer::LWORK, INFO, maxmn, minmn, M,N, RANK, a1, a2, error_code
    optional :: prec_func

    interface
       function prec_func(r)
         use System_module
         real(dp), intent(in) :: r
         real(dp) :: prec_func
       end function prec_func
    end interface

    if(twobody%N == 0) then
       call print("There are no parameters to optimise")
       return
    end if

    do_fromscratch = .true.
    if(present(fromscratch)) do_fromscratch = fromscratch

    do_svd = .true.
    if(present(method)) then
       if(trim(method) .eq. "SVD") then
          do_svd = .true.
       else if(trim(method) .eq. "minim") then
          do_svd = .false.
       else
          call system_abort("Adjustable_potential_optimise: invalid method '" // trim(method) // "'")
       end if
    end if

    !If the chosen method is SVD then the default is to run serially to
    !avoid memory bottlenecks when all processors do SVD.
    do_parallel_svd = optional_default(.false.,parallel_svd)

    if(size(target_force) /= size(fitforce)) &
         call system_abort("adjustable_potential_optimise: target force is wrong size")

    ! Calculate preconditioning factors for each spring
    allocate(precondition(twobody%N))
    if (present(prec_func)) then
       do i=1,twobody%N
          a1 = atomlist%int(1,twobody%int(1,i))
          a2 = atomlist%int(1,twobody%int(2,i))
          shift = twobody%int(3:5,i)
          precondition(i) = prec_func(distance(atoms_in, a1, a2, shift))
       end do
    else
       precondition = 1.0_dp
    end if

    ! create matrices for fast evaluation
    call adjustable_potential_create_forcematrix(atoms_in, precondition)

    target_force = reshape(fitforce, (/size(target_force)/))

    nparamtot = twobody%N*nparam
    allocate(params(nparamtot))

    if(.not.adjustable_potential_initialised) &
         call system_abort("adjustable_potential_optimise: not initialised!")

    call print("Target force:", PRINT_NERD)
    call print(target_force, PRINT_NERD)

    call print('Optimising '//nparamtot//' adjustable parameters')
    call print("RMS force component error before optimisation : "// &
         sqrt(normsq(target_force)/(real(atomlist%N,dp)*3.0_dp)))
    call print("Max force component error before optimisation : "//&
         maxval(abs(target_force)))

    if(do_fromscratch) then
       params = 0.0_dp
    else
       params = reshape(twobody%real(2:1+nparam,1:twobody%N), (/nparamtot/))
    endif

    if(do_svd) then

#ifdef _MPI

       if (do_parallel_svd .or. mpi_id()==0) then

#endif

       ! Solution by SVD
       call print('Using SVD for least squares fit, eigenvalue threshold = '//eigenvalue_threshold)
       allocate(full_forcematrix(size(target_force), size(params)))
       full_forcematrix = forcematrix
       M = size(full_forcematrix,1)
       N = size(full_forcematrix,2)
       minmn = minval((/M,N/))
       maxmn = maxval((/M,N/))
       allocate(X(maxmn,1))
       allocate(S(minmn))

       allocate(WORK(1))

       X(1:M,1) = target_force

       ! Do a workspace query first
       call dgelss(M, N, 1, full_forcematrix, M, X,maxmn,S, eigenvalue_threshold, RANK, WORK, -1, INFO )

       ! Allocate optimal size workspace
       LWORK = WORK(1)
       deallocate(WORK)
       allocate(WORK(LWORK))

       ! Do the SVD
       call dgelss(M, N, 1, full_forcematrix, M, X,maxmn,S, eigenvalue_threshold, RANK, WORK, LWORK, INFO )

       if (INFO /= 0) then
          call system_abort('adjustable_potential_optimise: DGELSS failed with error code '//INFO)
       end if

       params = X(1:N,1)
       deallocate(WORK)
       deallocate(S)
       deallocate(X)
       deallocate(full_forcematrix)

#ifdef _MPI

       endif

       !Broadcast the params to each process if SVD only done on proc 0
       if (.not.do_parallel_svd) then
          call MPI_BCAST(params,size(params),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,error_code)
          !check for errors
          call abort_on_mpi_error(error_code, "Adjustable_Potentia_Optimise: Params MPI_Bcast()")
       end if

#endif

    else
       !! Solution by minimisation
       steps = minim(params,&
            adjustable_potential_force_error, &
            adjustable_potential_force_error_deriv,&
            "cg",1e-10_dp,1000,"FAST_LINMIN")
    end if

    ! compute RMS after optimisation
    err = sqrt(adjustable_potential_force_error(params)/(real(atomlist%N,dp)*3.0_dp))

    ! Print all magnitudes of force errors if we're a NERD
    call print('Force errors', PRINT_NERD)
    ferr = reshape(target_force - (forcematrix .mult. params), &
         (/3, atomlist%N/))

!!$    call print('sparse force=')
!!$    call Print(forcematrix .mult. params)
!!$    call print('')

    do i=1,atomlist%N
       write (line, '(i4,f12.6)') i, sqrt(normsq(ferr(:,i)))
       call print(line, PRINT_NERD)
    end do
    call Print('', PRINT_NERD)

    write(line, '(a,e10.2)') "RMS force component error after  optimisation : ", err
    call print(line)
    write(line, '(a,e10.2)') "Max force component error after  optimisation : ", &
         maxval(abs(target_force - (forcematrix .mult. params)))
    call print(line)

    if (value(mainlog%verbosity_stack) >= PRINT_NERD) then
       call Print('Final Springs:', PRINT_NERD)
       do i=1,twobody%N
          write (line, '(3i8,4f16.8)') i, atomlist%int(1,twobody%int(1,i)), atomlist%int(1,twobody%int(2,i)), twobody%real(1,i), precondition(i), params(i), params(i)*precondition(i)
          call Print(line, PRINT_NERD)
       end do
       call Print('', PRINT_NERD)
    end if

    ! Get physical spring constants by multiplying by precondition factor
    do i=1,twobody%N
       params(i) = params(i)*precondition(i)
    end do

    twobody%real(2:1+nparam,1:twobody%N) = reshape(params, (/nparam,twobody%N/))

    write(line, '(a,e10.2)') "Max abs spring constant   after  optimisation : ", maxval(abs(params))
    call print(line)

    if (maxval(abs(params)) > adjustable_potential_max_spring_constant) then
       write(line, '(a,f10.2)') "WARNING: Max spring constant value after optimisation is large : ", &
            maxval(abs(params))
       call print(line)
       bad_spring = maxloc(abs(params),1)
       write(line,'(2(a,i0))')'         Max spring constant belongs to spring ', &
                               atomlist%int(1,twobody%int(1,bad_spring)),'--',atomlist%int(1,twobody%int(2,bad_spring))
       call print(line)
    end if

    call print('')

    deallocate(params)
    deallocate(precondition)

    ! calculate new parameter derivatives
    do i=1,min(twobody%N,twobody_old%N)
       do k=1,nparam
          twobody%real(1+nparam+k,i) = twobody_old%real(1+nparam+k, i)+2.0_dp*(twobody%real(1+k,i)-twobody_old%real(1+k,i)-twobody_old%real(1+nparam+k,i))
       end do
    end do

    ! Save the old fit force
    if (allocated(fitforce_old)) deallocate(fitforce_old)
    allocate(fitforce_old(size(fitforce,1),size(fitforce,2)))
    fitforce_old = fitforce

  end subroutine adjustable_potential_optimise

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X adjustable_potential_print_params()
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_print_params()
    real(dp)::f(size(target_force))
    integer::np

    if(.not.adjustable_potential_initialised) &
         call system_abort("adjustable_potential_print_params: not initialised!")

    write(line,*) "Atoms              Shift                   Parameters" ; call print(line)
    write(line,*) "-----------------------------------------------------"
    call print(line)
    call print(twobody)
    write(line, *) ; call print(line)

    write(line, *) "Force matrix: " ; call print(line)
    write(line, *) ; call print(line)
    call print_full(forcematrix)
    write(line, *) ; call print(line)

    write(line,*) "Forces from force matrix:" ; call print( line)
    np = twobody%N*adjustable_potential_nparams
    f = adjustable_potential_force_sparse(reshape(twobody%real(2:1+adjustable_potential_nparams,1:twobody%N), (/np/)))
    call print(transpose(reshape(f, (/3,atomlist%N/))))
    write(line, *) ; call print(line)

  end subroutine adjustable_potential_print_params

  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X Set up, alter and query the exclusion list
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

  subroutine adjustable_potential_add_exclusion(atom1,atom2,shift)

    integer, intent(in) :: atom1,atom2,shift(3)

    if (atom1 == atom2) then
       call print_warning('Adjustable_Potential_Add_Exclusion: Atom1 = Atom2')
       return
    end if

    if (.not.exclusions) then !Set up the exclusion table
       call allocate(exclusion_list,5,0,0,0)
       exclusions = .true.
    end if

    !Append (/a, b, shift/) where a < b
    if (atom1 < atom2) then
       call append(exclusion_list,(/atom1,atom2,shift/))
    else
       call append(exclusion_list,(/atom2,atom1,-shift/))
    end if

  end subroutine adjustable_potential_add_exclusion

  subroutine adjustable_potential_delete_exclusion(atom1,atom2,shift)

    integer, intent(in) :: atom1, atom2, shift(3)

    if (atom1 < atom2) then
       call delete(exclusion_list,(/atom1,atom2,shift/))
    else
       call delete(exclusion_list,(/atom2,atom1,-shift/))
    end if

    if (exclusion_list%N == 0) then
       call finalise(exclusion_list)
       exclusions = .false.
    end if

  end subroutine adjustable_potential_delete_exclusion

  function is_excluded(atom1,atom2, shift)

    integer, intent(in) :: atom1, atom2, shift(3)
    logical             :: Is_Excluded

    if (atom1 < atom2) then
       if (find(exclusion_list,(/atom1,atom2,shift/))==0) then
          Is_Excluded = .false.
       else
          Is_Excluded = .true.
       end if
    else
       if (find(exclusion_list,(/atom2,atom1,-shift/))==0) then
          Is_Excluded = .false.
       else
          Is_Excluded = .true.
       end if

    end if

  end function is_excluded

end module AdjustablePotential_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text: #redirect Latin American Restaurant Guide


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
        subroutine get_uv_3d(i4time,imax,jmax,kmax,u_3d,v_3d,ext,istatus
     1)

!       Returns U and V wind components from a 3D gridded LAPS dataset

!       Steve Albers            1990

        logical ltest_vertical_grid

        character*150 DIRECTORY
        character*31 EXT      ! Input extension of file (normally 3 characters)

        character*125 comment_3d(kmax)
        character*10 units_3d(kmax)
        character*3 var_u(kmax)
        character*3 var_v(kmax)
        integer LVL_3d(kmax)
        character*4 LVL_COORD_3d(kmax)

        real u_3d(imax,jmax,kmax)
        real v_3d(imax,jmax,kmax)

        call get_directory(ext,directory,len_dir)

        do k = 1,kmax
            units_3d(k)   = 'M/S'
            comment_3d(k) = '3DWIND'
            if(ltest_vertical_grid('HEIGHT'))then
                lvl_3d(k) = zcoord_of_level(k)/10
                lvl_coord_3d(k) = 'MSL'
            elseif(ltest_vertical_grid('PRESSURE'))then
                lvl_3d(k) = nint(zcoord_of_level(k))/100
                lvl_coord_3d(k) = 'HPA'
            else
                write(6,*)' Error, vertical grid not supported,'
     1                   ,' this routine supports PRESSURE or HEIGHT'
                istatus = 0
                return
            endif

            var_u(k) = 'U3' ! newvar = 'U3', oldvar = 'U'
            var_v(k) = 'V3' ! newvar = 'V3', oldvar = 'V'

        enddo ! k

!       Read in 3d U array
        CALL READ_LAPS_DATA(I4TIME,DIRECTORY,
     1          EXT,imax,jmax,kmax,kmax,
     1          VAR_u,LVL_3d,LVL_COORD_3d,UNITS_3d,COMMENT_3d,
     1          u_3d,ISTATUS)
        IF(ISTATUS .ne. 1)THEN
            write(6,*)
     1  ' Sorry, file has not yet been generated this hour'
        else
            write(6,*)
     1  ' 3D - U analysis successfully read in ',EXT(1:3)
        endif

!       Read in 3d V array
        CALL READ_LAPS_DATA(I4TIME,DIRECTORY,
     1          EXT,imax,jmax,kmax,kmax,
     1          VAR_v,LVL_3d,LVL_COORD_3d,UNITS_3d,COMMENT_3d,
     1          v_3d,ISTATUS)
        IF(ISTATUS .ne. 1)THEN
            write(6,*)
     1  ' Sorry, file has not yet been generated this hour'
!           stop
        else
            write(6,*)
     1  ' 3D - V analysis successfully read in ',EXT(1:3)
        endif

        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      DOUBLE PRECISION FUNCTION PHO_BESSI1(X)
C**********************************************************************
C
C      Bessel Function I1
C
C**********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION ax , besli1 , X , y
      SAVE

      ax = ABS(X)

      IF ( ax.LT.3.75D0 ) THEN
         y = (X/3.75D0)**2
         besli1 = ax*
     &            (0.5D0+y*(0.87890594D0+y*(0.51498869D0+y*(0.15084934D0
     &            +y*(0.2658733D-1+y*(0.301532D-2+y*0.32411D-3))))))
      ELSE
         y = 3.75D0/ax
         besli1 = 0.2282967D-1 +
     &            y*(-0.2895312D-1+y*(0.1787654D-1-y*0.420059D-2))
         besli1 = 0.39894228D0 +
     &            y*(-0.3988024D-1+y*(-0.362018D-2+y*(0.163801D-2+
     &            y*(-0.1031555D-1+y*besli1))))
         besli1 = besli1*EXP(ax)/SQRT(ax)
      END IF
      IF ( X.LT.0.D0 ) besli1 = -besli1

      PHO_BESSI1 = besli1

      END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET