anchor.diff:         ENDDO
   
      ENDDO
  -
  -
  -
  -   END SUBROUTINE BD_StaticElementMatrix
  +   CALL Cleanup()
  +CONTAINS
  +   SUBROUTINE Cleanup()
  +      IF(ALLOCATED(gp       ))  DEALLOCATE(gp       )
  +      IF(ALLOCATED(gw       ))  DEALLOCATE(gw       )
  +      IF(ALLOCATED(hhx      ))  DEALLOCATE(hhx      )
  +      IF(ALLOCATED(hpx      ))  DEALLOCATE(hpx      )
  +      IF(ALLOCATED(GLL_temp ))  DEALLOCATE(GLL_temp )
  +      IF(ALLOCATED(w_temp   ))  DEALLOCATE(w_temp   )
  +      IF(ALLOCATED(temp_Nvvv))  DEALLOCATE(temp_Nvvv)
  +      IF(ALLOCATED(temp_Naaa))  DEALLOCATE(temp_Naaa)
  +   END SUBROUTINE Cleanup
  +END SUBROUTINE BD_StaticElementMatrix
   
      SUBROUTINE BD_StaticUpdateConfiguration(uinc,uf,node_total,dof_node)
   
  
positive.code:    CHARACTER(ErrMsgLen)    :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER :: RoutineName = 'BD_StaticElementMatrix'
  
     vvv(:)   = 0.0D0
     Nvvv(:)  = 0.0D0
     elk(:,:) = 0.0D0
     elf(:)   = 0.0D0
     CALL BD_GenerateGLL(node_elem-1,GLL_temp,w_temp,ErrStat2,ErrMsg2)
     CALL BD_GaussPointWeight(ngp,gp,gw,ErrStat2,ErrMsg2)
     DO igp=1,ngp
         gpr = gp(igp)
         CALL BD_ComputeJacobian(gpr,Nuu0,node_elem,dof_node,gp,GLL_temp,ngp,igp,hhx,hpx,Jacobian,ErrStat2,ErrMsg2)
         CALL BD_GaussPointDataAt0(hhx,hpx,Nuu0,Nrr0,node_elem,dof_node,uu0,E10,ErrStat2,ErrMsg2)
         Stif(1:6,1:6) = EStif0_GL(1:6,1:6,igp)
         CALL BD_GaussPointData(hhx,hpx,Nuuu,Nrrr,uu0,E10,node_elem,dof_node,&
                                uuu,uup,E1,RR0,kapa,Stif,cet,ErrStat2,ErrMsg2)
         CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,Fc,Fd,Oe,Pe,Qe,ErrStat2,ErrMsg2)
         mmm          = EMass0_GL(1,1,igp)
         mEta(2)      =-EMass0_GL(1,6,igp)
         mEta(3)      = EMass0_GL(1,5,igp)
         rho(1:3,1:3) = EMass0_GL(4:6,4:6,igp)
         CALL BD_GaussPointDataMass(hhx,hpx,Nvvv,temp_Naaa,RR0,node_elem,dof_node,&
                                    vvv,temp_aaa,vvp,mmm,mEta,rho,ErrStat2,ErrMsg2)
         CALL BD_GravityForce(mmm,mEta,gravity,Fg,ErrStat2,ErrMsg2)
         Fd(:) = Fd(:) - Fg(:) - Distr_GL(:,igp)
  
         DO i=1,node_elem
  
negative.code:        ENDDO
         CALL BD_StaticElementMatrix(uuN0(:,nelem),Nuuu,Nrr0,Nrrr,DistrLoad_GL,gravity,&
                                     Mass0(:,:,temp_id+1:temp_id+ngp),Stif0(:,:,temp_id+1:temp_id+ngp),&
                                     ngp,node_elem,dof_node,elk,elf,ErrStat2,ErrMsg2)
  
         CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,elk,StifK,ErrStat2,ErrMsg2)
         CALL BD_AssembleRHS(nelem,dof_elem,node_elem,dof_node,elf,RHS,ErrStat2,ErrMsg2)
  