text:
David Martin is the owner of Cobalt Salon & Gallery.  David is a full service hairstylist, specializing in color, color correction, color education,  wet and dry cutting techniques, and photo styling. Everything David does is %100 guaranteed!  (If you dont love what I do Ill fix it or refund it, period.David)

David has worked in Davis for over 20 years.  He started as an apprentice with Bunny Jean Cunningham, of J. Cunningham Hair Salon at her downtown salon, J. Cunningham & Co., on the second floor of the Mansion Square Center, in 1987.  David worked with Jean for 4 years before moving to the salon that he would work in for the next fifteen years, All That Jaz.  He owned All That Jaz under the name All That Hair Artistry for the last five years he was there. After selling All That/Hair Artistry, David moved back downtown and worked briefly with Roxie Battaglia in her downtown salon, Roxies, before leaving to work with his friend Tenaya Freitas, at Davids Haircutting, in the Orange Court Center.  David and Tenaya moved to Salon At Maple Court the following year.  Here he would work with his old friends and mentors from J. Cunningham & Co., Sharon Jackson, and Cecilia Scott.  When circumstances forced them out of their H Street location David opened his own salon, Cobalt Salon & Gallery and invited his friends Tenaya, Sharon, and Cecilia to join him.

He is planning on staying at his current location for many, many, years.

David specializes in:
  http://davidmartinsstyle.com/techniques.htm haircolor (both permanent and semi/demi permanent techniques depending on the clients style choice)
  color correction
  hilighting, lolighting, paneling, baliage, ombre, and blockcolor techniques
  (very rarely) styling for special occasions such as updos and/or wedding hairstyles (full or half day packages if within a few hours drive) or for a night out.
  both wet and dry cutting and razor techniques.  I love Curly hair, too!!
David is one of the best haircolor artists in town, and his haircuts are always perfect, guaranteed!

Please visit Davids website at http://www.davidmartinsstyle.com for more information.

Links
  http://www.davislifemagazine.com  Cary Arnold is a long time client/friend, and David wrote for Davis Life Magazine for a time.
  http://music.ucdavis.edu/concertband/  Davids husband of 15 years Pete Nowlen
  http://www.groundcontrolstore.com/index.php  sister and brotherinlaws business
  http://www.capradio.org/programs/bios/biodetail.aspx?bioid7   Michael Frost is a client and friend of many years.
  http://ichthyologistbright.blogspot.com/   Davids exwifes blog, interesting, funny, infuriating, you decide...

Artwork
David works mainly in oil on canvas; abstract and realistic portraits are his specialty. David has shown at All That/Hair Artistry in July of 2006, and his work can be viewed by appointment in his studio in Sacramento, or at Cobalt Salon & Gallery.

To learn about other businesses in Davis that provide professional aesthetic services, visit our Beauty Salons and Barber Shops page.

20070515 18:44:11 nbsp Ive been going to David for longer than I care to remember.  Hes taken great care of my hair, and he makes it look nice even when it (the hair) suddenly adopts a new growth pattern!  Last summer I got some strange new waves to my hair when we cut it short for a trip.  It looked good even when I didnt blow it dry.  He also has good advice for keeping it healthy. Users/Christal

20070623 13:19:46 nbsp Ive been going to David for a haircut since I moved to Davis.  I have always been happy with the results! Users/DanWillenbring

20071103 21:31:03 nbsp I highly recommend David to all my friends. I LOVE my new haircut. I always had long hair as a child so I was very nervous about having short hair. I brought in pictures and told him roughly what I wanted. I ended up having a hair cut that was EXACTLY what I wanted. He will listen to how you want your hair and tell you all you need to know about what would look best for you. I am very excited to show everyone my new hair cut.  Users/tessa

20071210 12:53:49 nbsp I went to David for a cut last week.  I have curly hair and usually wear it short.  Most barbers/stylists dont know what to do with it.  David did a great job!  In addition to the great cut, the overall experience was wonderful.  I felt very comfortable and relaxed while there. Users/CharlesMcLaughlin

20080105 07:30:49 nbsp David always gets it right! I leave his salon loving my hair. Thanks David! Users/CaryArnold

20080319 16:20:23 nbsp I LOVE him! He does such a great job. I went to David to have may hair cut short and it turned out amazing. I would definitely recommend him.  Users/atnuech

20081104 09:36:22 nbsp All I can say is WOW!  David is amazing!  He performed magic on my hair for my wedding.  He made me feel like a total princess.  He not only did my hair but did 5 bridesmaids, 3 moms and two flower girls (with lots of help from my friend Sharon Gamble!!David).  He is totally professional and I cant believe how much he was able to do in a few hours!  I would highly recommend David to everyone I know!  Thank you David  you helped turn my day into a dream! Users/CaraJones

20081114 19:02:01 nbsp I went in for anything different.  David talked with me and was quite funny and friendly.  But most importantly, he gave me the most amazing haircut.  Its different from what I had before and totally me.  I DEFINITELY recommend him. Users/Tanguera

20090123 20:54:30 nbsp I have difficult hair to work with and I want some pretty challenging cuts. I asked him to mimic a photo of a previous cut with some alterations. He was able to handle the cut and do in 1.5 hours what took 2.5 at a highend salon in Austin (at a lower price too). David Martin is an artist. Users/Stefphanni

20090413 09:56:27 nbsp I had my first appointment with David just over two weeks ago  the cut and highlights were exactly what I wanted, the colour is particularly gorgeous! He was great to chat with, would highly recommend. Users/LizA

20090807 06:56:43 nbsp I was referred to David by a friend.  He did an outstanding job on her hair, especially on the color.  Now it is two weeks after my own appt. with him and I am very happy with both the cut and color.  I also appreciate how he gets his clients.  He listens.  Really listens.  My husband was, and still is, thrilled with Davids work.  He doesnt typically comment on my hair, but he calls me glamorous.  There have been subsequent date nights and a romantic weekend on the Mendocino coast.  With a three yearold in tow, but what the heck.  I recommend David without reservation, especially if you have fine hair and want to work with a true professional. Users/MrsC

20090811 17:28:18 nbsp Ive worked with David for years. Hes great. Very conscientious. He listens and finds out exactly what it is you want.   Users/TenayaFreitas

20100109 19:24:04 nbsp David is fantastic. He is an extremely talented hair stylist and colorist.  And a super nice guy too!  I highly recommend David.   Users/EBP

20100416 21:14:24   I have been devoted to Davids hair styling since I moved to Davis and feel that I was so lucky to have chosen him. I trust him completely and he never fails to put a smile on my face when I get out of his salon. Not only he is in total command of his craft, he completely devotes himself to make his clients happy. The quality of his professional work is exceptional. —ascaglione@ucdavis.edu

20100429 22:05:59 nbsp Had my first appointment with David today and must say that I love the work that he did with the color and cut. Super professional and knows what he is doing. I highly recommend him. Users/Yuriko

20100929 16:01:11   David is WONDERFUL! I got a hair cut with him as a gift from a friend and I couldnt be happier. He is professional, funny, kind, and really knows about hair! My boyfriend decided he finally needed to cut off his pony tail but was nervous as he hadnt had short hair in nearly a decade. He went to David for a consult and cut and came away with a great cut that works perfectly with his type of hair. I would highly recommend David! —SaraRoscoe

20110323 08:26:19 nbsp I try to search for a good salon in Davis. First of all, I chose Luci salon because my friends said they know how to style your hair. However, as I went online, I saw lots of good reviews about the Cobalt salon, so I decided to change my mind and go to this salon instead. It turns out that I got the worst haircut ever from this place. Although I brought the picture with me to the place, my haircut turns out so short and not even looks like the picture. It is completely a crazy short hair without any styling or any creative,or not even look like the pictures that I brought ( I brought 2 pictures, but if its not possible to cut as the pix, just tell me before the process, and dont say that its possible, but then it turns out to be worst at the end)The fullness of the hair in the back is also weird; the half top is flatted, and the half bottom is having fullness. My mom and aunts also agree so. I was crying whenever I look at the mirror.

I have to go to my regular stylist to ask if he can fix it for me, but it is impossible because if he fixes the hair, it will even shorter. Afterward, I had to go to Sally to ask for a shampoo to induce hair grow so it can grow faster and i can fix it. This is the first time that I have to pay for an expensive haircut and I dont like it at all. Im so disappointed, upset, and sorry for my self for not going to the suggested salon or stay with my stylist in my hometown.... So, Be careful with all the good reviews online.
PS. the owner, Davis is nice and friendly, hes just not a good hair stylist Users/ChauDoan

Ms. Doan, I am so very sorry for your bad experience with your haircut.  As I stated in the email I sent you, I would be happy to  refund your money if you do not want me to fix any problems with your cut.
 David Martin, owner

20110324 12:10:14 nbsp Id like to take a moment to come to Davids defense regarding the above comment.  David is wonderful and kind and has obviously tried to make things right with this unsatisfied customer.  Ive had bad hair cuts before (never from David) so I know how upsetting it can be, especially if you get a really short dramatic hair cut.  But David is a very talented stylist and I think her rant is unwarranted.  Also, hes only received one bad review in years of amazing styling, dont take it too hard David.   Users/SaraRoscoe
  How do you know her rant is unwarranted?  Everyone has a bad day, even the best of us, at what we do.  Perhaps she really wasnt satisfied with her hair and ya know what, thats alright in the grand scheme of things.  No one business is going to please everyone.users/OliviaY

20110324 13:12:31 nbsp Im in agreement with Sara on this one also. Ive been having David cut my hair for over a year now and have bringing my mother into him for nearly as long. Both of us have loved every single hair cut. Hes been a great listener, very knowledgeable, professional, and is funny and kind as well. Before David I was terrified of every little snip of my hair and now I look forward to working with him every time I get my hair cut. David, please dont let one bad review get you down because the rest of us think youre amazing! Users/ElysiaFrink

20110325 03:07:55 nbsp I call shannigans on Doans post. Youre so upset that youre going to create 5 new accounts on different sites and copy and paste your review in 5 different spots? And then refuse to take a refund? People lose creditability imo when they specifically reference other stylists/salons in a review. Nice plug.

While I do appreciate Bretts defense of me, I dont believe the salon referenced in Ms. Doans edit had anything to do with her review.David Martin

David is a good hair stylist. Maybe hes not good for you (hmmm). But you cant generalize for everyone else. Or try to discount all the good reviews for one bad (laughable).

Not only is David one of best cutter/stylists both me and my gf have ever been to. Hes pretty badass at color too. I should know cause ive had twotone mohawks and weve had all sorts of styles and colors.

Bottom line. David rocks. Users/BrettBretterson


text:
!==============================================================================!
  subroutine Comm_Mod_Open_File_Read(fh, file_name)
!------------------------------------------------------------------------------!
!   Open file for reading for parallel runs.                                   !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  integer   :: fh             ! file handle
  character :: file_name*(*)  ! file name
!-----------------------------------[Locals]-----------------------------------!
  integer :: error
!==============================================================================!

  ! Open file with MPI
  call Mpi_File_Open(MPI_COMM_WORLD,   &
                     file_name,        &
                     MPI_MODE_RDONLY,  &
                     MPI_INFO_NULL,    &
                     fh,               &
                     error)

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
      SUBROUTINE CHKMCV (STPATH, EXPTHS)

C  CHKMCV checks that the correct sequence of found paths is returned
C  by ISS as a result of the starting path.  It issues pass or fail
C  depending on whether the actual found paths match those expected.
C  It always uses a search reference point of 0,0,0 and distance of 2.

C  Input parameters:
C    STPATH : starting path for ISS
C    EXPTHS : expected found paths

C clipping indicator
C                noclip      clip
      INTEGER    PNCLIP,     PCLIP
      PARAMETER (PNCLIP = 0, PCLIP = 1)

      CHARACTER  STPATH*(*), EXPTHS*(*)

      CALL ISPTHS (0.,0.,0., 2., 1, PCLIP, STPATH, EXPTHS)

      END

ext: f


text:
PROGRAM kinds
IMPLICIT NONE
WRITE(*, '("The KIND for single precision is", I2)') KIND(0.0)
WRITE(*, '("The KIND for double precision is", I2)') KIND(0.0D0)
END PROGRAM kinds

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice


!!  Exponential integrator sweeper  module

! MODULE: pf_mod_exp
! !> @author
! Tommaso Buvoli
!
! Last Modified: Dec 28, 2018
!
!> Exponential integrator module
!!
!!  This module extends pf_sweeper_t and is used for creating an exponential sweeper
!!  that solves equations of the form
!!         $$   y' = L y + f(t,y)  $$
!!  When extending this class, you must supply the functions phib, swpPhib, and resPhib
!!  that each compute matrix-vector products of the form
!!         $$ \sum_{i=0}^n t^i \varphi_i(tL)b_i $$
!!  in addition to the function f_eval for compluting the nonlinear term N(t,y).
!!  The complete description of these three functions is contained below.
module pf_mod_fexp_sweeper

  use pf_mod_dtype
  use pf_mod_utils

  implicit none

  !> Exponential SDC sweeper type, extends abstract pf_sweeper_t
  type, extends(pf_sweeper_t), abstract :: pf_fexp_sweeper_t

    real(pfdp),        allocatable :: nodes(:)   ! nodes
    real(pfdp),        allocatable :: eta(:)     ! normalized substeps (on interval [0, 1])
    class(pf_encap_t), allocatable :: F_old(:)   ! scratch space for storing nonlinear terms
    class(pf_encap_t), allocatable :: newF       ! scratch space for storing new function evaluations

  contains

    ! Specialized procedures for exponential integrator
    procedure(pf_f_eval_p),           deferred :: f_eval    ! computes nonlinear term in equation
!    procedure(pf_comp_dt_p),          deferred :: comp_dt      !!  computes the time step
    procedure(pf_expSweepSubstep),    deferred :: expSweepSubstep
    procedure(pf_expResidualSubstep), deferred :: expResidualSubstep

    ! Generic Functions
    procedure :: initialize   => exp_initialize
    procedure :: sweep        => exp_sweep
    procedure :: evaluate     => exp_evaluate
    procedure :: integrate    => exp_integrate
    procedure :: residual     => exp_residual
    procedure :: spreadq0     => exp_spreadq0
    procedure :: compute_dt   => exp_compute_dt
    procedure :: evaluate_all => exp_evaluate_all
    procedure :: destroy      => exp_destroy

    ! functions that can be accessed directly by types that inherit pf_fexp_sweeper_t
    procedure :: exp_destroy
    procedure :: exp_initialize

  end type pf_fexp_sweeper_t

  interface ! DESCRIPTION OF REQUIRED FUNCTIONS

    ! =================================================================================
    ! EXPSWEEPSUBSTEP: Computes the jth substep of the exponential sweeper
    !
    !     		exp(x L) y_{n,j}^{k+1} + h \varphi_1(x L) [N_{n,j}^{[k+1]} - N_{n,j}^{[k]}] + I_{n,j}
    !
    !	where x = dt * (t_{n,j+1} - t_{n,j}), I_{j} are the exponential integrals
    !
    !		I_{n,j} = x \int_{0}^{1} exp(x L(1-s)) Q_j(s) ds
    !
    !	and Q_j(s) is a Lagrange interpolating polynomial that satisfies
    !
    !		Q_j(\tau_{j,i}) = N(y(t_{n,i})) for i = 1 ... m
    !
    !	for \tau_{j,i} = (t_{n,i} - t_{n,j}) / (t{n,j+1} - t_{n,j})
    !
    !	NOTE: The operator L is not passed in as a parameter, and must be
    !       implemented appropriately within the class.
    !
    ! Arguments
    !
    !   j   	(input) Integer
    !       	substep index for determining t: t = t_{n,j+1} - t_{n,j}
    !
    !   y_jp1 (inout) pf_encap_t
    !		      stores the solution y_{n,j+1}^{k+1}
    !
    !   dt  	(input) real
    !       	stepsize
    !
    !   y_j  (input) pf_encap_t
    !		     stores the solution y_{n,j}^{k+1}
    !
    !   F    (input) pf_encap_t(:,:)
    !		     should be the lev%F matrix of dimension (nnodes-1 x 2). First component
    !        stores the nonlinear term, and the second component is used by this
    !        function to store the exponential integral terms I_{n,j}
    !
    !   Nk  	(input) pf_encap_t(:)
    !		stores N(y^{k+1}_j) for j = 1 ... m
    ! =================================================================================

    subroutine pf_expSweepSubstep(this, y_jp1, j, dt, y_j, F, Nk)
      import pf_fexp_sweeper_t, pf_encap_t, pfdp
      class(pf_fexp_sweeper_t),   intent(inout) :: this
      class(pf_encap_t),  intent(inout) :: y_jp1
      integer,            intent(in)    :: j
      real(pfdp),         intent(in)    :: dt
      class(pf_encap_t),  intent(in)    :: y_j
      class(pf_encap_t),  intent(inout) :: F(:,:)
      class(pf_encap_t),  intent(inout) :: Nk(:)
    end subroutine pf_expSweepSubstep

    ! =================================================================================
    ! EXPRESSUBSTEP: Computes the jth residual
    !
    ! Arguments
    !
    !   j   	(input) Integer
    !       	substep index for determining t: t = t_{n,j+1} - t_{n,j}
    !
    !   dt  	(input) real
    !       	stepsize
    !
    !   y_n   	(input) pf_encap_t
    !		stores the solution y_n^{k+1}
    !
    !   F    (input) pf_encap_t(:,:)
    !		     should be the lev%F matrix of dimension (nnodes-1 x 2). First component
    !        stores the nonlinear term, and the second component is used by this
    !        function to store the exponential integral terms I_{n,j}
    !
    ! =================================================================================

    subroutine pf_expResidualSubstep(this, y_jp1, j, dt, y_j, F)
      import pf_fexp_sweeper_t, pf_encap_t, pfdp
      class(pf_fexp_sweeper_t),   intent(inout) :: this
      class(pf_encap_t),  intent(inout) :: y_jp1
      integer,            intent(in)    :: j
      real(pfdp),         intent(in)    :: dt
      class(pf_encap_t),  intent(in)    :: y_j
      class(pf_encap_t),  intent(in)    :: F(:,:)
    end subroutine pf_expResidualSubstep

    ! =================================================================================
    ! f_eval: computes the equations nonlinear term N(t,y)
    !
    ! Arguments
    !
    !   y     (input) pf_encap_t
    !         solution y(t)
    !
    !   t     (input) DOUBLE
    !         time t
    !
    !   level_index (input) INTEGER
    !         current level index
    !
    !   f     (output) pf_encap_t
    !         N(t,y)
    ! =================================================================================

    subroutine pf_f_eval_p(this, y, t, level_index, f)

      import pf_fexp_sweeper_t, pf_encap_t, pfdp
      class(pf_fexp_sweeper_t),   intent(inout) :: this

      class(pf_encap_t), intent(in)    :: y
      real(pfdp),        intent(in)    :: t
      integer,           intent(in)    :: level_index
      class(pf_encap_t), intent(inout) :: f
    end subroutine pf_f_eval_p
    subroutine pf_comp_dt_p(this,y, t, level_index, dt)
      !>  Set the time step
      import pf_fexp_sweeper_t, pf_encap_t, pfdp
      class(pf_fexp_sweeper_t),  intent(inout) :: this
      class(pf_encap_t), intent(in   ) :: y        !!  Argument for evaluation
      real(pfdp),        intent(in   ) :: t        !!  Time at evaluation
      integer,    intent(in   ) :: level_index     !!  Level index
      real(pfdp),        intent(inout) :: dt       !!  time step chosen
    end subroutine pf_comp_dt_p

  end interface

contains

  ! =================================================================================
  ! INITIALIZE: initializes the following internal parameters
  !      w        DOUBLE(:,:,:)     contains FD weights for computing local derivatives at t_{n,j}
  !      nodes    DOUBLE(:)         sdc nodes
  !      eta      DOUBLE(:)         normalized substeps (t_{n,j+1} - t_{n,j})/h
  !      npieces  INTEGER           number of RHS peices (always will be one)
  !      newF     pf_encap_t        stores new function evaluations
  !      b        pf_encap_t(:)     stores vectors b for computing phi products
  ! =================================================================================

  subroutine exp_initialize(this, pf,level_index)

    ! arguments
    class(pf_fexp_sweeper_t),   intent(inout) :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index

    ! local variables
    integer :: i, nnodes
    real(pfdp), allocatable :: q(:)

    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)
    nnodes = lev%nnodes
    allocate(this%eta(nnodes - 1))
    allocate(this%nodes(nnodes))
    allocate(q(nnodes))
    ! set nodes and substeps
    this%nodes = lev%sdcmats%qnodes
    this%eta = this%nodes(2 : nnodes) - this%nodes(1 : nnodes - 1) ! substeps
    ! set number of rhs components
    this%npieces = 2
    ! initialize temporary storage objects
    call lev%ulevel%factory%create_single(this%newF, lev%index, lev%lev_shape)
    call lev%ulevel%factory%create_array(this%F_old, nnodes, lev%index, lev%lev_shape)

  end subroutine exp_initialize

  ! SWEEP: exponential sweep subroutine ===============================================
  subroutine exp_sweep(this, pf, level_index, t0, dt,nsweeps, flags)

    use pf_mod_timer
    use pf_mod_hooks

    ! arguments
    class(pf_fexp_sweeper_t),   intent(inout) :: this
    type(pf_pfasst_t), intent(inout), target :: pf   !!  PFASST structure
    integer,           intent(in)    :: level_index  !!  which level to sweep on
    real(pfdp),        intent(in)    :: t0           !!  Time at beginning of time step
    real(pfdp),        intent(in)    :: dt           !!  time step size
    integer,           intent(in)    :: nsweeps      !!  number of sweeps to do
    integer, optional, intent(in)    :: flags

    ! local variables
    type(pf_level_t), pointer :: lev
    integer                    :: m, nnodes, j, k
    real(pfdp)                 :: t

    lev => pf%levels(level_index)
    nnodes = lev%nnodes
    ! error sweeps

    do k = 1, nsweeps
       call call_hooks(pf, level_index, PF_PRE_SWEEP)
       if (pf%save_timings > 1) call pf_start_timer(pf, T_SWEEP,level_index)
       pf%state%sweep=k

      ! NOTE: ensure that lev%F has been properly initialized here
      do j = 1, nnodes
        call this%F_old(j)%copy(lev%F(j,1))  ! Save old f
      end do
      if (k .eq. 1) then ! Is this necessary? it seems that lev%F(j,1) = F(lev%(Q,q)) or F_old definition above line would be incorrect or is it zero initially?
         call lev%Q(1)%copy(lev%q0)
         if (pf%save_timings > 1) call pf_start_timer(pf, T_FEVAL,level_index)
         call this%f_eval(lev%Q(1), t0, level_index, lev%F(1,1))      ! compute F_j^{[k+1]}
         if (pf%save_timings > 1) call pf_stop_timer(pf, T_FEVAL,level_index)
      end if
      t = t0
      do j = 1, nnodes - 1
        t = t0 + dt * this%nodes(j+1)

        call this%expSweepSubstep(lev%Q(j+1), j, dt, lev%Q(j), lev%F, this%F_old)	! compute exp(h_j L) y_{j} + \varphi_1(h_j L) * (F^{k+1}_j - F^{k}_j) + I_j
        !  Now we have to add in the tauQ
        if (level_index < pf%state%finest_level) then
          call lev%Q(j+1)%axpy(1.0_pfdp, lev%tauQ(j))
          if (j > 1) then     ! The tau is not node to node, so subtract out
            call lev%Q(j+1)%axpy(-1.0_pfdp, lev%tauQ(j-1))
          end if
       end if
       if (pf%save_timings > 1) call pf_start_timer(pf, T_FEVAL,level_index)
       call this%f_eval(lev%Q(j+1), t, level_index, lev%F(j+1,1))      			! compute F_j^{[k+1]}
       if (pf%save_timings > 1) call pf_stop_timer(pf, T_FEVAL,level_index)
      end do  !  Substepping over nodes

      call pf_residual(pf, level_index, dt)
      call lev%qend%copy(lev%Q(lev%nnodes))
      if (pf%save_timings > 1) call pf_stop_timer(pf, T_SWEEP,level_index)
      call call_hooks(pf, level_index, PF_POST_SWEEP)

   end do  !  Sweeps

  end subroutine exp_sweep

  ! =================================================================================
  ! INTEGRATE: computes the integrals for the exponential Picard residual
  !
  !         r_j = A_j - B_j
  !
  !     where
  !
  !         A_j = \left[ \exp(t h L) y(t_n) - \int_{t_n}^{t_{n,j}} \exp(L(t - t_n)) P(t) dt \right]
  !         B_j = y(t_{n,j})
  !
  !     NOTE: This procedure computes the expression \hat{r}_j = A_j - y(t_n). The
  !     term y(t_n) is subtracted from result since the generic calling function
  !     compute the residual as
  !
  !        r_j = y(t_n) + exp_integrate() - y(t_{n,j})
  !
  !     thus incorrectly adding the term y_n
  ! =================================================================================

  subroutine exp_integrate(this, pf,level_index, qSDC, fSDC, dt, fintsdc, flags)
    class(pf_fexp_sweeper_t),   intent(inout) :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    class(pf_encap_t), intent(in   ) :: qSDC(:)      !!  Solution values
    class(pf_encap_t), intent(in   ) :: fSDC(:, :)   !!  RHS Function values
    real(pfdp),        intent(in   ) :: dt           !!  Time step
    class(pf_encap_t), intent(inout) :: fintsdc(:)   !!  Integral from t_n to t_m
    integer, optional, intent(in   ) :: flags

    ! local variables
    integer :: i, nnodes
    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer

    nnodes = lev%nnodes
    do i = 1, nnodes - 1 ! loop over integrals : compute \int_{t_{n,i}}^{t_{n, i + 1}}
       call this%expResidualSubstep(fintsdc(i), i, dt, qSDC(i), fSDC)
       if (i .eq. 1) then
          call fintsdc(i)%axpy(-1.0_pfdp,qSDC(i))
       else
!      if (i > 1) then
          call fintsdc(i)%axpy(-1.0_pfdp,qSDC(i))
          call fintsdc(i)%axpy(1.0_pfdp,fintsdc(i-1))
      end if
    end do

  end subroutine exp_integrate

  ! RESIDUAL: compute  residual (generic) ====================================
  subroutine exp_residual(this, pf, level_index, dt, flags)

    class(pf_fexp_sweeper_t),  intent(inout)  :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    real(pfdp),             intent(in)    :: dt
    integer, intent(in), optional    :: flags

    integer :: m
    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer

    !>  Compute the integral of F from t_n to t_m at each node
    call lev%ulevel%sweeper%integrate(pf,level_index, lev%Q, lev%F, dt, lev%I, flags)

    !> add tau if it exists
    if (level_index < pf%state%finest_level) then
      do m = 1, lev%nnodes-1
        call lev%I(m)%axpy(1.0_pfdp, lev%tauQ(m), flags)
      end do
    end if

    !> subtract out the solution value
    do m = 1, lev%nnodes-1
      call lev%R(m)%copy(lev%I(m))
      call lev%R(m)%axpy(-1.0_pfdp, lev%Q(m+1))
       if (present(flags)) then
          if (flags .eq. 0) then
             call lev%R(m)%axpy(1.0_pfdp, lev%q0)
          else
             call lev%R(m)%axpy(1.0_pfdp, lev%Q(1))
          end if
       else
          call lev%R(m)%axpy(1.0_pfdp, lev%Q(1))
       end if
    end do

  end subroutine exp_residual

  ! SPREADQ: spread solution (generic) ======================================
  subroutine exp_spreadq0(this, pf,level_index, t0, flags, step)
    class(pf_fexp_sweeper_t),  intent(inout)  :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    real(pfdp),        intent(in   ) :: t0
    integer, optional,   intent(in)  :: flags, step

    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer

    call pf_generic_spreadq0(this, pf,level_index, t0)

  end subroutine exp_spreadq0
  !> Set the time step
  subroutine exp_compute_dt(this,pf,level_index,  t0, dt,flags)
    class(pf_fexp_sweeper_t),  intent(inout) :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    real(pfdp),        intent(in   ) :: t0
    real(pfdp),        intent(inout) :: dt
    integer, optional,   intent(in)    :: flags

    type(pf_level_t),    pointer :: lev
    lev => pf%levels(level_index)   !!  Assign level pointer
    !  Do nothing now
    return
  end subroutine exp_compute_dt

  ! EVALUATE: evaluate the nonlinear term at node m ========================
  subroutine exp_evaluate(this, pf,level_index, t, m, flags, step)
    ! arguments
    class(pf_fexp_sweeper_t),   intent(inout) :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    real(pfdp),        intent(in   ) :: t    !!  Time at which to evaluate
    integer,           intent(in   ) :: m    !!  Node at which to evaluate
    integer, intent(in), optional    :: flags, step

    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer
    if (pf%save_timings > 1) call pf_start_timer(pf, T_FEVAL,level_index)
    call this%f_eval(lev%Q(m), t, level_index, lev%F(m,1))
    if (pf%save_timings > 1) call pf_stop_timer(pf, T_FEVAL,level_index)
  end subroutine exp_evaluate

  ! EVALUATE_ALL: evaluate the nonlinear term at all nodes =================
  subroutine exp_evaluate_all(this, pf,level_index, t, flags, step)
    ! arguments
    class(pf_fexp_sweeper_t),  intent(inout)  :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index
    real(pfdp),        intent(in   ) :: t(:)     !!  Array of times at each node
    integer, intent(in), optional    :: flags, step

    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer

    call pf_generic_evaluate_all(this,pf, level_index, t)

  end subroutine exp_evaluate_all

  ! DEALLOCATE: deallocate sweeper variables
  subroutine exp_destroy(this, pf,level_index)
    ! arguments
    class(pf_fexp_sweeper_t),   intent(inout) :: this
    type(pf_pfasst_t), target, intent(inout) :: pf
    integer,              intent(in)    :: level_index

    type(pf_level_t), pointer :: lev
    lev => pf%levels(level_index)   !  Assign level pointer

    if(allocated(this%eta)) deallocate(this%eta)
    if(allocated(this%eta)) deallocate(this%newF)
    call lev%ulevel%factory%destroy_array(this%F_old)
  end subroutine exp_destroy

end module pf_mod_fexp_sweeper


text:
      subroutine wrihyd2(lundia,lun   ,bndval,typbnd,nambnd,nobnd ,
     *                   notims,kmax  ,nolay ,tstart,dtmin ,itdate)
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: wrihyd2.f 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/tools_gpl/nesthd2/packages/nesthd2/src/wrihyd2.f $
!***********************************************************************
! Deltares                         marine and coastal management
!
! subroutine         : wrihyd
! version            : v1.0
! date               : July 1997
! programmer         : Theo van der Kaaij
!
! function           : Writes hydrodynamic bc. to TRISULA attribute file
!                      (new format)
! limitations        :
! subroutines called :
!***********************************************************************
!
      character(len=37), dimension(13) :: fmtbct
      character(len=20), dimension(16) :: keywrd

      real         bndval (nobnd ,notims,kmax  ,1     ,2)

      character*1  eol
      character*1  typbnd (nobnd )
      character*1  quote
      character*3  cdum
      character*17 layer
      character*20 nambnd (nobnd )
      character*40 cntain

      data         fmtbct/'(a1,i5,                       t83,a1)',
     *                    '(a20,a1,a42,a3,a1,            t83,a1)',
     *                    '(2(a20,a1),a40,               t83,a1)',
     *                    '(2(a20,a1),                   t83,a1)',
     *                    '(a20,                         t83,a1)',
     *                    '(a20,                         t83,a1)',
     *                    '(a20,a1,a15,a1,               t83,a1)',
     *                    '(a20,i9,                      t83,a1)',
     *                    '(a20,a1,a7,a1,                t83,a1)',
     *                    '(a20,a1,a6,a1,                t83,a1)',
     *                    '(a20,a1,a37,a1,1x,2(a20,a1),  t83,a1)',
     *                    '(a20,i6,                      t83,a1)',
     *                    '(f12.4,    2ES14.6,           t83,a1)'/

      data         keywrd/'table-name          ',
     *                    'contents            ',
     *                    'location            ',
     *                    'geo-coordinates     ',
     *                    'metric-coordinates  ',
     *                    'frequency           ',
     *                    'frequency-unit      ',
     *                    'time-function       ',
     *                    'reference-time      ',
     *                    'time-unit           ',
     *                    'time-step           ',
     *                    'interpolation       ',
     *                    'extrapolation       ',
     *                    'parameter           ',
     *                    'unit                ',
     *                    'records-in-table    '/

!----------------------------------------------------------------
!     Initialisation
!----------------------------------------------------------------
      write (*     ,'('' >>> Writing hydrodynamic boundary '',
     *                ''conditions <<<'')')
      write (lundia,'('' >>> Writing hydrodynamic boundary '',
     *                ''conditions <<<'')')
      cdum  = '   '
      quote = CHAR (39)
      eol   = ' '
!
      cntain = ' # at ends A&B of open boundary segment '

      do 10 ifmt = 1, 13
         write (fmtbct(ifmt)(29:29),'(a1)') 't'
         write (fmtbct(ifmt)(30:33),'(i4.4)') max(12 + 2*kmax*14 + 1,83)
   10 continue

!----------------------------------------------------------------------
!     Write to time-series file
!----------------------------------------------------------------------

      do 100 ibnd = 1, nobnd
         do 100 itim = 1, notims
            if (itim .eq. 1) then
                write(lun   ,fmtbct( 2)) keywrd( 1),
     *                quote,'Nested flow BC '//nambnd(ibnd)//
     *                '      ',
     *                cdum,quote,eol
                if (typbnd(ibnd) .eq. 'z' .or.
     *              typbnd(ibnd) .eq. 'u' ) then
                   write(lun   ,fmtbct( 3)) keywrd( 2),
     *                   quote,'uniform             ',quote,cntain,eol
                else
                   write(lun   ,fmtbct( 3)) keywrd( 2),
     *                   quote,'3d-profile          ',quote,cntain,eol
                endif
                write(lun   ,fmtbct( 4)) keywrd( 3),
     *                quote,nambnd(ibnd),quote,eol
                write(lun   ,fmtbct( 7)) keywrd( 8),
     *                quote,'non-equidistant',quote,eol
                write(lun   ,fmtbct( 8)) keywrd( 9),
     *                itdate,eol
                write(lun   ,fmtbct( 9)) keywrd(10),
     *                quote,'minutes',quote,eol
                write(lun   ,fmtbct(10)) keywrd(12),
     *                quote,'linear',quote,eol
                write(lun   ,fmtbct(11)) keywrd(14),
     *                quote,'time                '//
     *                      '                 ',
     *                quote,keywrd(15),quote,' ',quote,eol
                if (typbnd(ibnd) .eq. 'z') then
                   layer = ' End A           '
                   write(lun   ,fmtbct(11)) keywrd(14),
     *                   quote,'Water elevation (Z) '//
     *                   layer,
     *                   quote,keywrd(15),quote,' ',quote,eol
                   layer = ' End B           '
                   write(lun   ,fmtbct(11)) keywrd(14),
     *                   quote,'Water elevation (Z) '//
     *                   layer,
     *                   quote,keywrd(15),quote,' ',quote,eol
                else
                   layer = ' End A layer:    '
                   do 110 k = 1, nolay
                      write (layer(15:17),'(a3)') cdum
                      write (layer(15:17),'(i3)') k
                      write(lun   ,fmtbct(11)) keywrd(14),
     *                   quote,'Current         (C) '//
     *                   layer,
     *                   quote,keywrd(15),quote,' ',quote,eol
  110              continue
                   layer = ' End B layer:    '
                   do 120 k = 1, nolay
                      write (layer(15:17),'(a3)') cdum
                      write (layer(15:17),'(i3)') k
                      write(lun   ,fmtbct(11)) keywrd(14),
     *                   quote,'Current         (C) '//
     *                   layer,
     *                   quote,keywrd(15),quote,' ',quote,eol
  120              continue
                endif
                write(lun   ,fmtbct(12)) keywrd(16),
     *                notims,eol

            endif

            if (typbnd(ibnd) .eq. 'z') then
              write (fmtbct(13) (10:12),'(i3.3)') 2
              write(lun   ,fmtbct(13))
     *              tstart +  (itim -1)*dtmin,
     *              (bndval (ibnd,itim,1,1,isize), isize=1,2),eol
            else
               write (fmtbct(13) (10:12),'(i3.3)') 2*nolay
               write(lun   ,fmtbct(13))
     *              tstart +  (itim -1)*dtmin,
     *              ((bndval (ibnd,itim,k,1,isize), k=    1,nolay),
     *                                              isize=1,2    ),eol
            endif
  100 continue

      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                                !!
!!  This file is part of SciFT project                                               !!
!!  Copyright (c) 2011-2016 Nestor F. Aguirre (nfaguirrec@gmail.com)                 !!
!!                                                                                !!
!!  Redistribution and use in source and binary forms, with or without               !!
!!  modification, are permitted provided that the following conditions are met:      !!
!!                                                                                !!
!!  1. Redistributions of source code must retain the above copyright notice, this   !!
!!     list of conditions and the following disclaimer.                              !!
!!  2. Redistributions in binary form must reproduce the above copyright notice,     !!
!!     this list of conditions and the following disclaimer in the documentation     !!
!!     and/or other materials provided with the distribution.                        !!
!!  3. Neither the name of the copyright holders nor the names of its contributors   !!
!!     may be used to endorse or promote products derived from this software         !!
!!     without specific prior written permission.                                    !!
!!                                                                                !!
!!  The copyright holders provide no reassurances that the source code provided      !!
!!  does not infringe any patent, copyright, or any other intellectual property      !!
!!  rights of third parties.  The copyright holders disclaim any liability to any    !!
!!  recipient for claims brought against recipient by any third party for            !!
!!  infringement of that parties intellectual property rights.                       !!
!!                                                                                !!
!!  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND  !!
!!  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED    !!
!!  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE           !!
!!  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR  !!
!!  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES   !!
!!  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;     !!
!!  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND      !!
!!  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT       !!
!!  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS    !!
!!  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                     !!
!!                                                                                !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>
!! @brief Test program
!!
program main
	use GOptions_
	use RNFunction_
	use CNFunction_
	implicit none

	integer :: argc
	character(1000) :: fileNameA, fileNameB
	type(RNFunction) :: rNFuncA, rNFuncB
	type(CNFunction) :: cNFuncA, cNFuncB
	integer :: fileTypeA, fileTypeB
	real(8) :: diff

    argc = command_argument_count()

    if( argc < 2 ) then
            write(*,*) "usage:"
            write(*,*) "   n1df.compare fileA fileB"
            write(*,*) ""
            write(*,*) "   diff² = 1.0-|<fA|fB>|²/( |fA|² * |fB|² )"
            stop
    end if

    call get_command_argument( 1, fileNameA )
    call get_command_argument( 2, fileNameB )

    write(0,*) "### ERROR ### Command not implemented yet"//trim(fileNameB)

! 	fileTypeA = cNFunc.checkTypeN1DF( fileNameA )
! 	fileTypeB = cNFunc.checkTypeN1DF( fileNameB )
!
! 	if( fileTypeA == 0 ) then
! 		call rNFunc.init( fileNameA )
! 	else if( fileTypeA == 1 ) then
		call cNFuncA.init( fileNameA )
! 	else
! 		write(0,*) "### ERROR ### unknown format for "//trim(fileNameA)
! 		stop
! 	end if
!
! 	if( fileTypeB == 0 ) then
! 		call rNFunc.init( fileNameB )
! 	else if( fileTypeB == 1 ) then
		call cNFuncB.init( fileNameB )
! 	else
! 		write(0,*) "### ERROR ### unknown format for "//trim(fileNameB)
! 		stop
! 	end if

	diff = sqrt( 1.0_8-abs( cNFuncA.innerProduct( cNFuncB ) )**2/cNFuncA.norm()**2/cNFuncB.norm()**2 )
	write(*,*) diff
end program main

max_stars_repo_name: nfaguirrec/scift


subroutine assert_true(x)
    implicit none
    logical :: x

    if (.not. x) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_true

subroutine assert_false(x)
    implicit none
    logical :: x

    if (x) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_false

subroutine assert_eq_byte(x, y)
    implicit none
    integer(kind=selected_int_kind(2)) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_byte

subroutine assert_eq_int16(x, y)
    implicit none
    integer(kind=selected_int_kind(4)) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_int16

subroutine assert_eq_integer(x, y)
    implicit none
    integer :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_integer

subroutine assert_eq_long_int(x, y)
    implicit none
    integer (selected_int_kind(18)) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_long_int

subroutine assert_eq_logical(x, y)
    implicit none
    logical :: x, y

    if (x .neqv. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_logical

subroutine assert_eq_character(x, y)
    implicit none
    character(len=*) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_character

subroutine assert_eq_single(x, y)
    implicit none
    real (selected_real_kind(6)) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_single

subroutine assert_eq_double(x, y)
    implicit none
    real (selected_real_kind(15)) :: x, y

    if (x .ne. y) then
        print *, 'Assertion failed'
        stop
    end if
end subroutine assert_eq_double

subroutine test_data_basic_types
    use libmuscle

    implicit none

    integer, parameter :: byte = selected_int_kind(2)
    integer, parameter :: short_int = selected_int_kind(4)
    integer, parameter :: long_int = selected_int_kind(18)
    type(LIBMUSCLE_Data) :: d1, d2
    integer :: i1, err_code
    logical :: l1
    character(len=:), allocatable :: s1, err_msg
    integer (kind=byte) :: bi1
    integer (kind=short_int) :: si1
    integer (kind=long_int) :: li1
    real (selected_real_kind(6)) :: sr1
    real (selected_real_kind(15)) :: dr1

    print *, '[  RUN     ] data.nil'
    d1 = LIBMUSCLE_Data_create()
    call assert_true(LIBMUSCLE_Data_is_nil(d1))
    call assert_false(LIBMUSCLE_Data_is_a_int(d1))
    call LIBMUSCLE_Data_set(d1, 10)
    call assert_false(LIBMUSCLE_Data_is_nil(d1))
    call LIBMUSCLE_Data_set_nil(d1)
    call assert_true(LIBMUSCLE_Data_is_nil(d1))
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.nil'

    print *, '[  RUN     ] data.logical'
    d1 = LIBMUSCLE_Data_create(.true.)
    call assert_true(LIBMUSCLE_Data_is_a_bool(d1))
    call assert_false(LIBMUSCLE_Data_is_a_int(d1))
    l1 = LIBMUSCLE_Data_as_bool(d1)
    call assert_eq_logical(l1, .true.)
    s1 = LIBMUSCLE_Data_as_string(d1, err_code, err_msg)
    call assert_eq_integer(err_code, LIBMUSCLE_runtime_error)
    deallocate(err_msg)
    call LIBMUSCLE_Data_set(d1, .false.)
    call assert_eq_logical(LIBMUSCLE_Data_as_bool(d1), .false.)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.logical'

    print *, '[  RUN     ] data.string'
    d1 = LIBMUSCLE_Data_create('Testing')
    call assert_true(LIBMUSCLE_Data_is_a_string(d1))
    call assert_false(LIBMUSCLE_Data_is_a_bool(d1))
    s1 = LIBMUSCLE_Data_as_string(d1)
    call assert_eq_character(s1, 'Testing')
    deallocate(s1)
    i1 = LIBMUSCLE_Data_as_int(d1, err_code)
    call assert_eq_integer(err_code, LIBMUSCLE_runtime_error)
    call LIBMUSCLE_Data_set(d1, 'Something else')
    s1 = LIBMUSCLE_Data_as_string(d1)
    call assert_eq_character(LIBMUSCLE_Data_as_string(d1), 'Something else')
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.string'

    print *, '[  RUN     ] data.int8'
    d1 = LIBMUSCLE_Data_create(121_byte)
    call assert_true(LIBMUSCLE_Data_is_a_char(d1))
    call assert_true(LIBMUSCLE_Data_is_a_int(d1))
    call assert_false(LIBMUSCLE_Data_is_a_string(d1))
    bi1 = LIBMUSCLE_Data_as_char(d1)
    call assert_eq_byte(bi1, 121_byte)
    call LIBMUSCLE_Data_set(d1, 43_byte)
    call assert_eq_byte(LIBMUSCLE_Data_as_char(d1), 43_byte)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.int8'

    print *, '[  RUN     ] data.int16'
    d1 = LIBMUSCLE_Data_create(1313_short_int)
    call assert_true(LIBMUSCLE_Data_is_a_int16(d1))
    si1 = LIBMUSCLE_Data_as_int16(d1)
    call assert_eq_int16(si1, 1313_short_int)
    call LIBMUSCLE_Data_set(d1, 24242_short_int)
    call assert_eq_int16(LIBMUSCLE_Data_as_int16(d1), 24242_short_int)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.int16'

    print *, '[  RUN     ] data.integer'
    d1 = LIBMUSCLE_Data_create(13)
    call assert_true(LIBMUSCLE_Data_is_a_int(d1))
    call assert_false(LIBMUSCLE_Data_is_a_bool(d1))
    i1 = LIBMUSCLE_Data_as_int(d1)
    call assert_eq_integer(i1, 13)
    call LIBMUSCLE_Data_set(d1, 242424)
    call assert_eq_integer(LIBMUSCLE_Data_as_int(d1), 242424)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.integer'

    print *, '[  RUN     ] data.integer64'
    d1 = LIBMUSCLE_Data_create(131313131313131313_long_int)
    call assert_true(LIBMUSCLE_Data_is_a_int64(d1))
    call assert_false(LIBMUSCLE_Data_is_a_string(d1))
    li1 = LIBMUSCLE_Data_as_int64(d1)
    call assert_eq_long_int(li1, 131313131313131313_long_int)
    call LIBMUSCLE_Data_set(d1, 242424242424242424_long_int)
    call assert_eq_long_int(LIBMUSCLE_Data_as_int64(d1), 242424242424242424_long_int)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.integer64'

    print *, '[  RUN     ] data.single'
    d1 = LIBMUSCLE_Data_create(42.0)
    call assert_true(LIBMUSCLE_Data_is_a_float(d1))
    call assert_false(LIBMUSCLE_Data_is_a_int(d1))
    sr1 = LIBMUSCLE_Data_as_float(d1)
    call assert_eq_single(sr1, 42.0)
    call LIBMUSCLE_Data_set(d1, 3.1415926536)
    call assert_eq_single(LIBMUSCLE_Data_as_float(d1), 3.1415926536)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.single'

    print *, '[  RUN     ] data.double'
    d1 = LIBMUSCLE_Data_create(42.0d0)
    call assert_true(LIBMUSCLE_Data_is_a_double(d1))
    call assert_false(LIBMUSCLE_Data_is_a_float(d1))
    dr1 = LIBMUSCLE_Data_as_double(d1)
    call assert_eq_double(dr1, 42.0d0)
    call LIBMUSCLE_Data_set(d1, 3.1415926536d0)
    call assert_eq_double(LIBMUSCLE_Data_as_double(d1), 3.1415926536d0)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.double'

end subroutine test_data_basic_types

subroutine test_data_copy_assign
    use libmuscle

    implicit none
    type(LIBMUSCLE_Data) :: d1, d2

    print *, '[  RUN     ] data.copy_constructor'
    d1 = LIBMUSCLE_Data_create(42.0d0)
    d2 = LIBMUSCLE_Data_create(d1)
    call LIBMUSCLE_Data_free(d1)
    call assert_true(LIBMUSCLE_Data_is_a_double(d2))
    call assert_false(LIBMUSCLE_Data_is_a_float(d2))
    call LIBMUSCLE_Data_free(d2)
    print *, '[       OK ] data.copy_constructor'

    print *, '[  RUN     ] data.assign'
    d1 = LIBMUSCLE_Data_create(42.0d0)
    d2 = LIBMUSCLE_Data_create()
    call LIBMUSCLE_Data_set(d2, d1)
    call LIBMUSCLE_Data_free(d1)
    call assert_true(LIBMUSCLE_Data_is_a_double(d2))
    call assert_false(LIBMUSCLE_Data_is_a_float(d2))
    call LIBMUSCLE_Data_free(d2)
    print *, '[       OK ] data.assign'
end subroutine test_data_copy_assign

subroutine test_data_dict
    use libmuscle

    implicit none
    integer, parameter :: byte = selected_int_kind(2)
    integer, parameter :: short_int = selected_int_kind(4)
    integer, parameter :: long_int = selected_int_kind(18)
    type(LIBMUSCLE_Data) :: d1, d2
    integer (kind=long_int) :: i

    print *, '[  RUN     ] data.dict'
    d1 = LIBMUSCLE_Data_create_dict()
    call assert_true(LIBMUSCLE_Data_is_a_dict(d1))
    call LIBMUSCLE_Data_set_item(d1, 'key1', 1)
    call assert_eq_long_int(LIBMUSCLE_Data_size(d1), 1_long_int)

    call LIBMUSCLE_Data_set_item(d1, 'key1', .true.)
    d2 = LIBMUSCLE_Data_get_item(d1, 'key1')
    call assert_true(LIBMUSCLE_Data_as_bool(d2))
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_set_item(d1, 'key2', 'test')
    d2 = LIBMUSCLE_Data_get_item(d1, 'key2')
    call assert_eq_character(LIBMUSCLE_Data_as_string(d2), 'test')
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_set_item(d1, 'key1', 63_byte)
    d2 = LIBMUSCLE_Data_get_item(d1, 'key1')
    call assert_eq_byte(LIBMUSCLE_Data_as_char(d2), 63_byte)
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_set_item(d1, 'key1', 30000_short_int)
    d2 = LIBMUSCLE_Data_get_item(d1, 'key1')
    call assert_eq_int16(LIBMUSCLE_Data_as_int16(d2), 30000_short_int)
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_set_item(d1, 'key2', 1000030000)
    d2 = LIBMUSCLE_Data_get_item(d1, 'key2')
    call assert_eq_integer(LIBMUSCLE_Data_as_int(d2), 1000030000)
    call LIBMUSCLE_Data_free(d2)

    call assert_eq_character(LIBMUSCLE_Data_key(d1, 1_long_int), 'key1')
    d2 = LIBMUSCLE_Data_value(d1, 1_long_int)
    call assert_eq_int16(LIBMUSCLE_Data_as_int16(d2), 30000_short_int)
    call LIBMUSCLE_Data_free(d2)

    call assert_eq_character(LIBMUSCLE_Data_key(d1, 2_long_int), 'key2')
    d2 = LIBMUSCLE_Data_value(d1, 2_long_int)
    call assert_eq_integer(LIBMUSCLE_Data_as_int(d2), 1000030000)
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.dict'
end subroutine test_data_dict

subroutine test_data_list
    use libmuscle

    implicit none
    integer, parameter :: byte = selected_int_kind(2)
    integer, parameter :: short_int = selected_int_kind(4)
    integer, parameter :: long_int = selected_int_kind(18)
    integer, parameter :: real4 = selected_real_kind(6)
    integer, parameter :: real8 = selected_real_kind(15)
    type(LIBMUSCLE_Data) :: d1, d2

    print *, '[  RUN     ] data.list'
    d1 = LIBMUSCLE_Data_create_list()
    call assert_true(LIBMUSCLE_Data_is_a_list(d1))
    call assert_false(LIBMUSCLE_Data_is_a_float(d1))
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.list'

    print *, '[  RUN     ] data.nils'
    d1 = LIBMUSCLE_Data_create_nils(100_long_int)
    call assert_true(LIBMUSCLE_Data_is_a_list(d1))
    call assert_false(LIBMUSCLE_Data_is_a_float(d1))
    call assert_eq_long_int(LIBMUSCLE_Data_size(d1), 100_long_int)
    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.nils'

    print *, '[  RUN     ] data.list_get_item'
    d1 = LIBMUSCLE_Data_create_nils(10_long_int)
    d2 = LIBMUSCLE_Data_get_item(d1, 1_long_int)
    call assert_true(LIBMUSCLE_Data_is_nil(d2))
    call LIBMUSCLE_Data_free(d1)
    call LIBMUSCLE_Data_free(d2)
    print *, '[       OK ] data.list_get_item'

    print *, '[  RUN     ] data.list_set_item'
    d1 = LIBMUSCLE_Data_create_nils(9_long_int)
    call LIBMUSCLE_Data_set_item(d1, 1_long_int, .true.)
    call LIBMUSCLE_Data_set_item(d1, 2_long_int, 'Testing')
    call LIBMUSCLE_Data_set_item(d1, 3_long_int, 42_byte)
    call LIBMUSCLE_Data_set_item(d1, 4_long_int, 6565_short_int)
    call LIBMUSCLE_Data_set_item(d1, 5_long_int, 1313131313)
    call LIBMUSCLE_Data_set_item(d1, 6_long_int, 100100100100100_long_int)
    call LIBMUSCLE_Data_set_item(d1, 7_long_int, 12.34_real4)
    call LIBMUSCLE_Data_set_item(d1, 8_long_int, 56.78901234_real8)
    d2 = LIBMUSCLE_Data_create_dict()
    call LIBMUSCLE_Data_set_item(d1, 9_long_int, d2)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 1_long_int)
    call assert_true(LIBMUSCLE_Data_as_bool(d2))
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 2_long_int)
    call assert_eq_character(LIBMUSCLE_Data_as_string(d2), 'Testing')
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 3_long_int)
    call assert_eq_byte(LIBMUSCLE_Data_as_char(d2), 42_byte)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 4_long_int)
    call assert_eq_int16(LIBMUSCLE_Data_as_int16(d2), 6565_short_int)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 5_long_int)
    call assert_eq_integer(LIBMUSCLE_Data_as_int(d2), 1313131313)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 6_long_int)
    call assert_eq_long_int(LIBMUSCLE_Data_as_int64(d2), 100100100100100_long_int)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 7_long_int)
    call assert_eq_single(LIBMUSCLE_Data_as_float(d2), 12.34_real4)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 8_long_int)
    call assert_eq_double(LIBMUSCLE_Data_as_double(d2), 56.78901234_real8)
    call LIBMUSCLE_Data_free(d2)

    d2 = LIBMUSCLE_Data_get_item(d1, 9_long_int)
    call assert_true(LIBMUSCLE_Data_is_a_dict(d2))
    call assert_eq_long_int(LIBMUSCLE_Data_size(d2), 0_long_int)
    call LIBMUSCLE_Data_free(d2)

    call LIBMUSCLE_Data_free(d1)
    print *, '[       OK ] data.list_set_item'
end subroutine test_data_list

subroutine test_data_byte_array
    use libmuscle

    implicit none
    integer, parameter :: long_int = selected_int_kind(18)
    character(len=1), dimension(:), allocatable :: buf
    integer :: err_code
    type(LIBMUSCLE_Data) :: d1

    print *, '[  RUN     ] data.byte_array'
    d1 = LIBMUSCLE_Data_create_byte_array(1024_long_int)
    call assert_true(LIBMUSCLE_Data_is_a_byte_array(d1))
    call assert_false(LIBMUSCLE_Data_is_a_int(d1))
    call assert_eq_long_int(LIBMUSCLE_Data_size(d1), 1024_long_int)

    allocate(buf(LIBMUSCLE_Data_size(d1)))
    call LIBMUSCLE_Data_as_byte_array(d1, buf, err_code)
    call assert_eq_integer(err_code, LIBMUSCLE_success)
    call LIBMUSCLE_Data_free(d1)
    deallocate(buf)
    print *, '[       OK ] data.byte_array'
end subroutine test_data_byte_array

subroutine test_data
    call test_data_basic_types
    call test_data_copy_assign
    call test_data_dict
    call test_data_list
    call test_data_byte_array
end subroutine

program test_fortran_api
    implicit none

    print *, ''
    print *, '[==========] Fortran API test'

    call test_data

    print *, '[==========] Fortran API test'
    print *, '[  PASSED  ] Fortran API test'
    print *, ''
end program test_fortran_api


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine domain

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------

      nx = nx0
      ny = ny0
      nz = nz0

c---------------------------------------------------------------------
c   check the sub-domain size
c---------------------------------------------------------------------
      if ( ( nx .lt. 4 ) .or.
     >     ( ny .lt. 4 ) .or.
     >     ( nz .lt. 4 ) ) then
         write (*,2001) nx, ny, nz
 2001    format (5x,'SUBDOMAIN SIZE IS TOO SMALL - ',
     >        /5x,'ADJUST PROBLEM SIZE OR NUMBER OF PROCESSORS',
     >        /5x,'SO THAT NX, NY AND NZ ARE GREATER THAN OR EQUAL',
     >        /5x,'TO 4 THEY ARE CURRENTLY', 3I3)
         stop
      end if

      if ( ( nx .gt. isiz1 ) .or.
     >     ( ny .gt. isiz2 ) .or.
     >     ( nz .gt. isiz3 ) ) then
         write (*,2002) nx, ny, nz
 2002    format (5x,'SUBDOMAIN SIZE IS TOO LARGE - ',
     >        /5x,'ADJUST PROBLEM SIZE OR NUMBER OF PROCESSORS',
     >        /5x,'SO THAT NX, NY AND NZ ARE LESS THAN OR EQUAL TO ',
     >        /5x,'ISIZ1, ISIZ2 AND ISIZ3 RESPECTIVELY.  THEY ARE',
     >        /5x,'CURRENTLY', 3I4)
         stop
      end if

c---------------------------------------------------------------------
c   set up the start and end in i and j extents for all processors
c---------------------------------------------------------------------
      ist = 2
      iend = nx - 1

      jst = 2
      jend = ny - 1

      ii1 = 2
      ii2 = nx0 - 1
      ji1 = 2
      ji2 = ny0 - 2
      ki1 = 3
      ki2 = nz0 - 1

      return
      end

max_stars_repo_name: jonarbo/KUBE

max_issues_repo_name: jonarbo/KUBE

max_forks_repo_name: jonarbo/KUBE


text:
      module inversion_module

c***********************************************************************
c
c     dl_poly module for defining inversion potential arrays
c     copyright - daresbury laboratory
c
c     author  - w. smith  sep 2003
c     adapted - w. smith  jun 2008 : solvation, free energy, excitation
c     adapted - w. smith  jan 2011 : metadynamics
c
c***********************************************************************

      use config_module
      use error_module
      use metafreeze_module
      use parse_module
      use property_module
      use setup_module
      use site_module
      use solvation_module

      implicit none

      real(8), allocatable :: prminv(:,:)
      integer, allocatable :: listinv(:,:)
      integer, allocatable :: numinv(:),keyinv(:),lstinv(:,:)

      save prminv,listinv,numinv,keyinv,lstinv

      contains

      subroutine alloc_inv_arrays(idnode,mxnode)

      implicit none

      logical safe
      integer i,fail,idnode,mxnode
      dimension fail(5)

      safe=.true.

c     allocate arrays

      fail(:)=0

      allocate (prminv(mxtinv,mxpinv),stat=fail(1))
      allocate (numinv(mxtmls),stat=fail(2))
      allocate (keyinv(mxtinv),stat=fail(3))
      allocate (lstinv(mxtinv,4),stat=fail(4))
      allocate (listinv(mxinv,5),stat=fail(5))

      if(any(fail.gt.0))safe=.false.
      if(mxnode.gt.1)call gstate(safe)
      if(.not.safe)call error(idnode,1120)

c     initialise numin array

      do i=1,mxtmls
         numinv(i)=0
      enddo

      end subroutine alloc_inv_arrays

      subroutine define_inversions
     x  (safe,idnode,itmols,ninver,nsite,engunit)

c***********************************************************************
c
c     dl_poly subroutine for defining inversion angle potentials
c     note: includes calcite planar potential which is not strictly
c     an inversion potential and is different in implementation
c
c     copyright - daresbury laboratory
c     author    - w. smith    aug 2003
c
c***********************************************************************

      implicit none

      logical safe
      character*8 keyword
      character*1 message(80)
      integer idnode,itmols,ninver,nsite,ntmp,inv,j,keytmp
      integer iatm1,iatm2,iatm3,iatm4,isite1,isite2,isite3,isite4
      integer ia,ja,idum
      real(8) engunit,parpot(mxpinv)

      ntmp=intstr(record,lenrec,idum)
      numinv(itmols)=numinv(itmols)+ntmp

      if(idnode.eq.0)then
        write(nrite,"(/,1x,'number of inversion terms',
     x    6x,i10)")numinv(itmols)
        Write(nrite,"(/,1x,'inversion angle details:',
     x       /,/,12x,'unit',5x,'key',5x,'index',5x,'index',5x,'index',
     x       5x,'index',7x,'f-const',8x,'angle',8x,'factor',/)")
      endif

      do inv=1,ntmp

        ninver=ninver+1
        if(ninver.gt.mxtinv)call error(idnode,73)

c     read inversion potential parameters

        call getrec(safe,idnode,nfield)
        if(.not.safe)return

        call copystring(record,message,80)
        call lowcase(record,4)
        call getword(keyword,record,4,lenrec)

        if(keyword(1:4).eq.'harm')then
          keytmp=1
        elseif(keyword(1:4).eq.'hcos')then
          keytmp=2
        elseif(keyword(1:4).eq.'plan')then
          keytmp=3
        elseif(keyword(1:4).eq.'calc')then
          keytmp=4
        else
          if(idnode.eq.0)write(nrite,*)record
          call error(idnode,449)
        endif

        iatm1=intstr(record,lenrec,idum)
        iatm2=intstr(record,lenrec,idum)
        iatm3=intstr(record,lenrec,idum)
        iatm4=intstr(record,lenrec,idum)
        parpot(1)=dblstr(record,lenrec,idum)
        parpot(2)=dblstr(record,lenrec,idum)

        isite1=nsite-numsit(itmols)+iatm1
        isite2=nsite-numsit(itmols)+iatm2
        isite3=nsite-numsit(itmols)+iatm3
        isite4=nsite-numsit(itmols)+iatm4

c     test for frozen atom pairs

        if(lfzsit(isite1)*lfzsit(isite2)*
     x    lfzsit(isite3)*lfzsit(isite4).ne.0)then

          numinv(itmols)=numinv(itmols)-1
          if(idnode.eq.0)
     x      write(nrite,"(4x,a8,i10,4x,a4,4i10,10f15.6)")
     x      '*frozen*',inv,keyword(1:4),iatm1,iatm2,iatm3,iatm4,
     x      (parpot(j),j=1,mxpinv)

        else

          if(idnode.eq.0)
     x      write(nrite,"(12x,i10,4x,a4,4i10,10f15.6)")
     x      inv,keyword(:4),iatm1,iatm2,iatm3,iatm4,
     x      (parpot(j),j=1,mxpinv)

        endif

        keyinv(ninver)=keytmp
        lstinv(ninver,1)=iatm1
        lstinv(ninver,2)=iatm2
        lstinv(ninver,3)=iatm3
        lstinv(ninver,4)=iatm4
        prminv(ninver,1)=parpot(1)
        prminv(ninver,2)=parpot(2)

c     convert energies to internal units and angles to radians

        prminv(ninver,1)=prminv(ninver,1)*engunit

        if(keytmp.eq.2)then

          prminv(ninver,2)=cos(prminv(ninver,2)*(pi/180.d0))

        elseif(keytmp.eq.4)then

          prminv(ninver,2)=prminv(ninver,2)*engunit

        endif

      enddo

      return
      end subroutine define_inversions

      subroutine invfrc
     x  (lsolva,lfree,lexcite,idnode,imcon,mxnode,ntinv,enginv,virinv)

c***********************************************************************
c
c     dl_poly subroutine for calculating inversion energy and force
c     terms in molecular dynamics.
c
c     copyright - daresbury laboratory 1996
c     author    - w. smith       may   1996
c     modified  - w.smith        jan   2011 : metadynamics
c
c***********************************************************************

      implicit none

      logical safe,lsolva,lfree,lexcite,lselect,lcalcite
      logical idrive,jdrive,kdrive,ldrive
      integer idnode,imcon,mxnode,ntinv,fail1,fail2
      integer fail3,inv1,inv2,i,ii,ia,id,kk,ib,ic
      real(8) strs(6),strs_loc(6)
      real(8) xab,yab,zab,rab2,rrab,xac,yac,zac,rac2,rrac,xad,yad
      real(8) zad,rad2,rrad,rbc,rcd,rdb,ubx,uby,ubz,ubn,rub,vbx
      real(8) vby,vbz,rvb,wwb,ucx,ucy,ucz,ucn,vcx,vcy,vcz,rvc,wwc
      real(8) udx,udy,udz,udn,vdx,vdy,vdz,vdn,rvd,wwd,cosb,cosc
      real(8) cosd,thb,thc,thd,gamb,gamc,gamd,rubc,rubd,rucd,rucb
      real(8) rudb,rudc,rvbc,rvbd,rvcd,rvcb,rvdb,rvdc,fax,fay,faz
      real(8) fbx,fby,fbz,fcx,fcy,fcz,fdx,fdy,fdz,vbn,vcn,ruc,rud
      real(8) uuu,uu2,uun,uux,uuy,uuz,enginv,virinv,omega,gamma
      real(8) gamvir

      real(8), allocatable :: xdab(:),ydab(:),zdab(:)
      real(8), allocatable :: xdac(:),ydac(:),zdac(:)
      real(8), allocatable :: xdad(:),ydad(:),zdad(:)

      data fail1,fail2,fail3/0,0,0/

      allocate (xdab(msbad),ydab(msbad),zdab(msbad),stat=fail1)
      allocate (xdac(msbad),ydac(msbad),zdac(msbad),stat=fail2)
      allocate (xdad(msbad),ydad(msbad),zdad(msbad),stat=fail3)
      if(fail1.ne.0.or.fail2.ne.0.or.fail3.ne.0)
     x     call error(idnode,1130)

c     check size of work arrays

      if((ntinv-mxnode+1)/mxnode.gt.msbad)call error(idnode,427)

c     block indices

      inv1=(idnode*ntinv)/mxnode+1
      inv2=((idnode+1)*ntinv)/mxnode

      safe=.true.

c     initialise accumulators

      enginv=0.d0
      virinv=0.d0
      inv_fre=0.d0
      inv_vir=0.d0
      strs(:)=0.d0
      strs_loc(:)=0.d0

      if(lsolva)then

        lcomp(4)=.true.
        inv_sol(:)=0.d0
        if(lexcite)inv_exc(:)=0.d0

      endif

c     calculate bond vectors

      ii=0
      do i=inv1,inv2

        ii=ii+1

c     potential energy function type

        kk=listinv(ii,1)

c     indices of bonded atoms

        ia=listinv(ii,2)
        ib=listinv(ii,3)
        ic=listinv(ii,4)
        id=listinv(ii,5)

c     define components of bond vectors

        xdab(ii)=xxx(ib)-xxx(ia)
        ydab(ii)=yyy(ib)-yyy(ia)
        zdab(ii)=zzz(ib)-zzz(ia)

        if(keyinv(kk).eq.4)then

          xdac(ii)=xxx(ic)-xxx(ib)
          ydac(ii)=yyy(ic)-yyy(ib)
          zdac(ii)=zzz(ic)-zzz(ib)

          xdad(ii)=xxx(id)-xxx(ib)
          ydad(ii)=yyy(id)-yyy(ib)
          zdad(ii)=zzz(id)-zzz(ib)

        else

          xdac(ii)=xxx(ic)-xxx(ia)
          ydac(ii)=yyy(ic)-yyy(ia)
          zdac(ii)=zzz(ic)-zzz(ia)

          xdad(ii)=xxx(id)-xxx(ia)
          ydad(ii)=yyy(id)-yyy(ia)
          zdad(ii)=zzz(id)-zzz(ia)

        endif

      enddo

c     periodic boundary condition

      call images(imcon,0,1,ii,cell,xdab,ydab,zdab)
      call images(imcon,0,1,ii,cell,xdac,ydac,zdac)
      call images(imcon,0,1,ii,cell,xdad,ydad,zdad)

c     loop over all specified inversions

      ii=0
      do i=inv1,inv2

        ii=ii+1

c     select potential energy function type

        kk=listinv(ii,1)
        lcalcite=(keyinv(kk).eq.4)

c     define components of bond vectors

        xab=xdab(ii)
        yab=ydab(ii)
        zab=zdab(ii)
        rab2=xab*xab+yab*yab+zab*zab
        rrab=1.d0/sqrt(rab2)

        xac=xdac(ii)
        yac=ydac(ii)
        zac=zdac(ii)
        rac2=xac*xac+yac*yac+zac*zac
        rrac=1.d0/sqrt(rac2)

        xad=xdad(ii)
        yad=ydad(ii)
        zad=zdad(ii)
        rad2=xad*xad+yad*yad+zad*zad
        rrad=1.d0/sqrt(rad2)

        if(lcalcite)then

c     calculate vector normal to plane

          uux=yac*zad-zac*yad
          uuy=zac*xad-xac*zad
          uuz=xac*yad-yac*xad
          uun=1.d0/sqrt(uux**2+uuy**2+uuz**2)
          uux=uun*uux
          uuy=uun*uuy
          uuz=uun*uuz
          uuu=xab*uux+yab*uuy+zab*uuz

        else

c     scalar products of bond vectors

          rbc=xab*xac+yab*yac+zab*zac
          rcd=xac*xad+yac*yad+zac*zad
          rdb=xad*xab+yad*yab+zad*zab

c     calculate bond-angle-plane vectors

          ubx=xac*rrac+xad*rrad
          uby=yac*rrac+yad*rrad
          ubz=zac*rrac+zad*rrad
          ubn=1.d0/sqrt(ubx**2+uby**2+ubz**2)
          ubx=ubn*ubx
          uby=ubn*uby
          ubz=ubn*ubz
          rub=xab*ubx+yab*uby+zab*ubz

          vbx=xac*rrac-xad*rrad
          vby=yac*rrac-yad*rrad
          vbz=zac*rrac-zad*rrad
          vbn=1.d0/sqrt(vbx**2+vby**2+vbz**2)
          vbx=vbn*vbx
          vby=vbn*vby
          vbz=vbn*vbz
          rvb=xab*vbx+yab*vby+zab*vbz
          wwb=sqrt(rub**2+rvb**2)

          ucx=xad*rrad+xab*rrab
          ucy=yad*rrad+yab*rrab
          ucz=zad*rrad+zab*rrab
          ucn=1.d0/sqrt(ucx**2+ucy**2+ucz**2)
          ucx=ucn*ucx
          ucy=ucn*ucy
          ucz=ucn*ucz
          ruc=xac*ucx+yac*ucy+zac*ucz

          vcx=xad*rrad-xab*rrab
          vcy=yad*rrad-yab*rrab
          vcz=zad*rrad-zab*rrab
          vcn=1.d0/sqrt(vcx**2+vcy**2+vcz**2)
          vcx=vcn*vcx
          vcy=vcn*vcy
          vcz=vcn*vcz
          rvc=xac*vcx+yac*vcy+zac*vcz
          wwc=sqrt(ruc**2+rvc**2)

          udx=xab*rrab+xac*rrac
          udy=yab*rrab+yac*rrac
          udz=zab*rrab+zac*rrac
          udn=1.d0/sqrt(udx**2+udy**2+udz**2)
          udx=udn*udx
          udy=udn*udy
          udz=udn*udz
          rud=xad*udx+yad*udy+zad*udz

          vdx=xab*rrab-xac*rrac
          vdy=yab*rrab-yac*rrac
          vdz=zab*rrab-zac*rrac
          vdn=1.d0/sqrt(vdx**2+vdy**2+vdz**2)
          vdx=vdn*vdx
          vdy=vdn*vdy
          vdz=vdn*vdz
          rvd=xad*vdx+yad*vdy+zad*vdz
          wwd=sqrt(rud**2+rvd**2)

c     calculate inversion angle cosines

          cosb=wwb*rrab
          cosc=wwc*rrac
          cosd=wwd*rrad

        endif

c     calculate potential energy and scalar force term

        if(keyinv(kk).eq.1)then

c     key=1 for harmonic inversion potential

          thb=acos(cosb)
          thc=acos(cosc)
          thd=acos(cosd)
          omega=0.5d0*prminv(kk,1)*((thb-prminv(kk,2))**2+
     x      (thc-prminv(kk,2))**2+(thd-prminv(kk,2))**2)/3.d0
          gamvir=0.d0

          gamb=0.d0
          if(abs(thb).gt.1.d-12)
     x      gamb=prminv(kk,1)*(thb-prminv(kk,2))/(3.d0*sin(thb))
          gamc=0.d0
          if(abs(thc).gt.1.d-12)
     x      gamc=prminv(kk,1)*(thc-prminv(kk,2))/(3.d0*sin(thc))
          gamd=0.d0
          if(abs(thd).gt.1.d-12)
     x      gamd=prminv(kk,1)*(thd-prminv(kk,2))/(3.d0*sin(thd))

        else if(keyinv(kk).eq.2)then

c     key=2 for harmonic cosine inversion potential

          omega=0.5d0*prminv(kk,1)*((cosb-prminv(kk,2))**2+
     x      (cosc-prminv(kk,2))**2+(cosb-prminv(kk,2))**2)/3.d0
          gamvir=0.d0

          gamb=-prminv(kk,1)*(cosb-prminv(kk,2))/3.d0
          gamc=-prminv(kk,1)*(cosc-prminv(kk,2))/3.d0
          gamd=-prminv(kk,1)*(cosd-prminv(kk,2))/3.d0

        else if(keyinv(kk).eq.3)then

c     key=3 for planar inversion potentials

          omega=prminv(kk,1)*((1.d0-cosb)+(1.d0-cosc)+(1.d0-cosd))/3.d0
          gamvir=0.d0

          gamb=prminv(kk,1)/3.d0
          gamc=prminv(kk,1)/3.d0
          gamd=prminv(kk,1)/3.d0

        else if(keyinv(kk).eq.4)then

c     key=4 for planar calcite potential

          uu2=uuu*uuu
          omega=uu2*(prminv(kk,1)+prminv(kk,2)*uu2)
          gamvir=uu2*(2.d0*prminv(kk,1)+4.d0*prminv(kk,2)*uu2)
          gamma=-uuu*(2.d0*prminv(kk,1)+4.d0*prminv(kk,2)*uu2)

        else

c     undefined potential

          safe=.false.
          gamvir=0.d0
          gamb=0.d0
          gamc=0.d0
          gamd=0.d0

        endif

c     indices of bonded atoms

        ia=listinv(ii,2)
        ib=listinv(ii,3)
        ic=listinv(ii,4)
        id=listinv(ii,5)

c     metadynamics local definitions

        if(lmetadyn)then

          idrive=driven(ltype(ia))
          jdrive=driven(ltype(ib))
          kdrive=driven(ltype(ic))
          ldrive=driven(ltype(id))

        endif

c     set selection control

        lselect=.true.

        if(lexcite)then

c     selected excitation option

          if((atm_fre(ia).ne.1).and.(atm_fre(ib).ne.1).and.
     x      (atm_fre(ic).ne.1).and.(atm_fre(id).ne.1))then

c     reset selection control

            lselect=(atm_fre(ia)+atm_fre(ib)+atm_fre(ic)+
     x        atm_fre(id).eq.0)

            if(lsolva)then
              inv_exc(atmolt(ia))=inv_exc(atmolt(ia))+omega
            endif

          endif

        elseif(lfree)then

c     selected free energy option

          if((atm_fre(ia).eq.1).or.(atm_fre(ib).eq.1).or.
     x      (atm_fre(ic).eq.1).or.(atm_fre(id).eq.1))then

c     set hamiltonian mixing parameter

            inv_fre=inv_fre-omega
            omega=lambda1*omega
            gamb=lambda1*gamb
            gamc=lambda1*gamc
            gamd=lambda1*gamd

          elseif((atm_fre(ia).eq.2).or.(atm_fre(ib).eq.2).or.
     x        (atm_fre(ic).eq.2).or.(atm_fre(id).eq.2))then

c     set hamiltonian mixing parameter

            inv_fre=inv_fre+omega
            omega=lambda2*omega
            gamb=lambda2*gamb
            gamc=lambda2*gamc
            gamd=lambda2*gamd

          endif

        endif

c     calculate bond and u,v scalar products

        if(.not.lcalcite)then

          rubc=xab*ucx+yab*ucy+zab*ucz
          rubd=xab*udx+yab*udy+zab*udz
          rucd=xac*udx+yac*udy+zac*udz
          rucb=xac*ubx+yac*uby+zac*ubz
          rudb=xad*ubx+yad*uby+zad*ubz
          rudc=xad*ucx+yad*ucy+zad*ucz

          rvbc=xab*vcx+yab*vcy+zab*vcz
          rvbd=xab*vdx+yab*vdy+zab*vdz
          rvcd=xac*vdx+yac*vdy+zac*vdz
          rvcb=xac*vbx+yac*vby+zac*vbz
          rvdb=xad*vbx+yad*vby+zad*vbz
          rvdc=xad*vcx+yad*vcy+zad*vcz

        endif

        if(lselect)then

c     calculate potential energy and virial

          enginv=enginv+omega
          virinv=virinv+gamvir

c     calculate solvation energy

          if(lsolva)then
            inv_sol(atmolt(ia))=inv_sol(atmolt(ia))+omega
          endif

c     calculate atomic forces

          if(lcalcite)then

            fax=-gamma*uux
            fay=-gamma*uuy
            faz=-gamma*uuz

            fcx=gamma*uun*((yad*zab-zad*yab)-uuu*(yad*uuz-zad*uuy))
            fcy=gamma*uun*((zad*xab-xad*zab)-uuu*(zad*uux-xad*uuz))
            fcz=gamma*uun*((xad*yab-yad*xab)-uuu*(xad*uuy-yad*uux))

            fdx=gamma*uun*((yab*zac-zab*yac)-uuu*(zac*uuy-yac*uuz))
            fdy=gamma*uun*((zab*xac-xab*zac)-uuu*(xac*uuz-zac*uux))
            fdz=gamma*uun*((xab*yac-yab*xac)-uuu*(yac*uux-xac*uuy))

            fbx=-(fax+fcx+fdx)
            fby=-(fay+fcy+fdy)
            fbz=-(faz+fcz+fdz)

c     stress tensor calculation for calcite terms

            strs(1)=strs(1)+uuu*gamma*uux*uux
            strs(2)=strs(2)+uuu*gamma*uux*uuy
            strs(3)=strs(3)+uuu*gamma*uux*uuz
            strs(4)=strs(4)+uuu*gamma*uuy*uuy
            strs(5)=strs(5)+uuu*gamma*uuy*uuz
            strs(6)=strs(6)+uuu*gamma*uuz*uuz

          else

            fbx=gamb*(-cosb*xab*rrab**2+rrab*(rub*ubx+rvb*vbx)/wwb)
     x      +(ruc*ucn*rrab*(xac-ruc*ucx-(rbc-ruc*rubc)*xab*rrab**2)
     x      -rvc*vcn*rrab*(xac-rvc*vcx-(rbc-rvc*rvbc)*xab*rrab**2))
     x      *gamc*rrac/wwc
     x      +(rud*udn*rrab*(xad-rud*udx-(rdb-rud*rubd)*xab*rrab**2)
     x      +rvd*vdn*rrab*(xad-rvd*vdx-(rdb-rvd*rvbd)*xab*rrab**2))
     x      *gamd*rrad/wwd

            fby=gamb*(-cosb*yab*rrab**2+rrab*(rub*uby+rvb*vby)/wwb)
     x      +(ruc*ucn*rrab*(yac-ruc*ucy-(rbc-ruc*rubc)*yab*rrab**2)
     x      -rvc*vcn*rrab*(yac-rvc*vcy-(rbc-rvc*rvbc)*yab*rrab**2))
     x      *gamc*rrac/wwc
     x      +(rud*udn*rrab*(yad-rud*udy-(rdb-rud*rubd)*yab*rrab**2)
     x      +rvd*vdn*rrab*(yad-rvd*vdy-(rdb-rvd*rvbd)*yab*rrab**2))
     x      *gamd*rrad/wwd

            fbz=gamb*(-cosb*zab*rrab**2+rrab*(rub*ubz+rvb*vbz)/wwb)
     x      +(ruc*ucn*rrab*(zac-ruc*ucz-(rbc-ruc*rubc)*zab*rrab**2)
     x      -rvc*vcn*rrab*(zac-rvc*vcz-(rbc-rvc*rvbc)*zab*rrab**2))
     x      *gamc*rrac/wwc
     x      +(rud*udn*rrab*(zad-rud*udz-(rdb-rud*rubd)*zab*rrab**2)
     x      +rvd*vdn*rrab*(zad-rvd*vdz-(rdb-rvd*rvbd)*zab*rrab**2))
     x      *gamd*rrad/wwd

            fcx=gamc*(-cosc*xac*rrac**2+rrac*(ruc*ucx+rvc*vcx)/wwc)
     x      +(rud*udn*rrac*(xad-rud*udx-(rcd-rud*rucd)*xac*rrac**2)
     x      -rvd*vdn*rrac*(xad-rvd*vdx-(rcd-rvd*rvcd)*xac*rrac**2))
     x      *gamd*rrad/wwd
     x      +(rub*ubn*rrac*(xab-rub*ubx-(rbc-rub*rucb)*xac*rrac**2)
     x      +rvb*vbn*rrac*(xab-rvb*vbx-(rbc-rvb*rvcb)*xac*rrac**2))
     x      *gamb*rrab/wwb

            fcy=gamc*(-cosc*yac*rrac**2+rrac*(ruc*ucy+rvc*vcy)/wwc)
     x      +(rud*udn*rrac*(yad-rud*udy-(rcd-rud*rucd)*yac*rrac**2)
     x      -rvd*vdn*rrac*(yad-rvd*vdy-(rcd-rvd*rvcd)*yac*rrac**2))
     x      *gamd*rrad/wwd
     x      +(rub*ubn*rrac*(yab-rub*uby-(rbc-rub*rucb)*yac*rrac**2)
     x      +rvb*vbn*rrac*(yab-rvb*vby-(rbc-rvb*rvcb)*yac*rrac**2))
     x      *gamb*rrab/wwb

            fcz=gamc*(-cosc*zac*rrac**2+rrac*(ruc*ucz+rvc*vcz)/wwc)
     x      +(rud*udn*rrac*(zad-rud*udz-(rcd-rud*rucd)*zac*rrac**2)
     x      -rvd*vdn*rrac*(zad-rvd*vdz-(rcd-rvd*rvcd)*zac*rrac**2))
     x      *gamd*rrad/wwd
     x      +(rub*ubn*rrac*(zab-rub*ubz-(rbc-rub*rucb)*zac*rrac**2)
     x      +rvb*vbn*rrac*(zab-rvb*vbz-(rbc-rvb*rvcb)*zac*rrac**2))
     x      *gamb*rrab/wwb

            fdx=gamd*(-cosd*xad*rrad**2+rrad*(rud*udx+rvd*vdx)/wwd)
     x      +(rub*ubn*rrad*(xab-rub*ubx-(rdb-rub*rudb)*xad*rrad**2)
     x      -rvb*vbn*rrad*(xab-rvb*vbx-(rdb-rvb*rvdb)*xad*rrad**2))
     x      *gamb*rrab/wwb
     x      +(ruc*ucn*rrad*(xac-ruc*ucx-(rcd-ruc*rudc)*xad*rrad**2)
     x      +rvc*vcn*rrad*(xac-rvc*vcx-(rcd-rvc*rvdc)*xad*rrad**2))
     x      *gamc*rrac/wwc

            fdy=gamd*(-cosd*yad*rrad**2+rrad*(rud*udy+rvd*vdy)/wwd)
     x      +(rub*ubn*rrad*(yab-rub*uby-(rdb-rub*rudb)*yad*rrad**2)
     x      -rvb*vbn*rrad*(yab-rvb*vby-(rdb-rvb*rvdb)*yad*rrad**2))
     x      *gamb*rrab/wwb
     x      +(ruc*ucn*rrad*(yac-ruc*ucy-(rcd-ruc*rudc)*yad*rrad**2)
     x      +rvc*vcn*rrad*(yac-rvc*vcy-(rcd-rvc*rvdc)*yad*rrad**2))
     x      *gamc*rrac/wwc

            fdz=gamd*(-cosd*zad*rrad**2+rrad*(rud*udz+rvd*vdz)/wwd)
     x      +(rub*ubn*rrad*(zab-rub*ubz-(rdb-rub*rudb)*zad*rrad**2)
     x      -rvb*vbn*rrad*(zab-rvb*vbz-(rdb-rvb*rvdb)*zad*rrad**2))
     x      *gamb*rrab/wwb
     x      +(ruc*ucn*rrad*(zac-ruc*ucz-(rcd-ruc*rudc)*zad*rrad**2)
     x      +rvc*vcn*rrad*(zac-rvc*vcz-(rcd-rvc*rvdc)*zad*rrad**2))
     x      *gamc*rrac/wwc

            fax=-(fbx+fcx+fdx)
            fay=-(fby+fcy+fdy)
            faz=-(fbz+fcz+fdz)

c     stress tensor calculation for inversion terms

            strs(1)=strs(1)+xab*fbx+xac*fcx+xad*fdx
            strs(2)=strs(2)+yab*fbx+yac*fcx+yad*fdx
            strs(3)=strs(3)+zab*fbx+zac*fcx+zad*fdx
            strs(4)=strs(4)+yab*fby+yac*fcy+yad*fdy
            strs(5)=strs(5)+yab*fbz+yac*fcz+yad*fdz
            strs(6)=strs(6)+zab*fbz+zac*fcz+zad*fdz

          endif

          fxx(ia)=fxx(ia)+fax
          fyy(ia)=fyy(ia)+fay
          fzz(ia)=fzz(ia)+faz

          fxx(ib)=fxx(ib)+fbx
          fyy(ib)=fyy(ib)+fby
          fzz(ib)=fzz(ib)+fbz

          fxx(ic)=fxx(ic)+fcx
          fyy(ic)=fyy(ic)+fcy
          fzz(ic)=fzz(ic)+fcz

          fxx(id)=fxx(id)+fdx
          fyy(id)=fyy(id)+fdy
          fzz(id)=fzz(id)+fdz

        endif

c     metadynamics local parameters

        if(lmetadyn.and.(idrive.or.jdrive.or.kdrive.or.ldrive))then

c     local energy and virial

          eng_loc=eng_loc+omega
          vir_loc=vir_loc+gamvir

c     local forces

          fxx_loc(ia)=fxx_loc(ia)+fax
          fyy_loc(ia)=fyy_loc(ia)+fay
          fzz_loc(ia)=fzz_loc(ia)+faz

          fxx_loc(ib)=fxx_loc(ib)+fbx
          fyy_loc(ib)=fyy_loc(ib)+fby
          fzz_loc(ib)=fzz_loc(ib)+fbz

          fxx_loc(ic)=fxx_loc(ic)+fcx
          fyy_loc(ic)=fyy_loc(ic)+fcy
          fzz_loc(ic)=fzz_loc(ic)+fcz

          fxx_loc(id)=fxx_loc(id)+fdx
          fyy_loc(id)=fyy_loc(id)+fdy
          fzz_loc(id)=fzz_loc(id)+fdz

c     local stress tensor

          strs_loc(1)=strs_loc(1)+xab*fbx+xac*fcx+xad*fdx
          strs_loc(2)=strs_loc(2)+yab*fbx+yac*fcx+yad*fdx
          strs_loc(3)=strs_loc(3)+zab*fbx+zac*fcx+zad*fdx
          strs_loc(4)=strs_loc(4)+yab*fby+yac*fcy+yad*fdy
          strs_loc(5)=strs_loc(5)+yab*fbz+yac*fcz+yad*fdz
          strs_loc(6)=strs_loc(6)+zab*fbz+zac*fcz+zad*fdz

        endif

      enddo

c     complete stress tensor

      stress(1)=stress(1)+strs(1)
      stress(2)=stress(2)+strs(2)
      stress(3)=stress(3)+strs(3)
      stress(4)=stress(4)+strs(2)
      stress(5)=stress(5)+strs(4)
      stress(6)=stress(6)+strs(5)
      stress(7)=stress(7)+strs(3)
      stress(8)=stress(8)+strs(5)
      stress(9)=stress(9)+strs(6)

      if(lmetadyn)then

        stress_loc(1)=stress_loc(1)+strs_loc(1)
        stress_loc(2)=stress_loc(2)+strs_loc(2)
        stress_loc(3)=stress_loc(3)+strs_loc(3)
        stress_loc(4)=stress_loc(4)+strs_loc(2)
        stress_loc(5)=stress_loc(5)+strs_loc(4)
        stress_loc(6)=stress_loc(6)+strs_loc(5)
        stress_loc(7)=stress_loc(7)+strs_loc(3)
        stress_loc(8)=stress_loc(8)+strs_loc(5)
        stress_loc(9)=stress_loc(9)+strs_loc(6)

      endif

c     check for undefined potentials

      if(mxnode.gt.1)call gstate(safe)
      if(.not.safe)call error(idnode,449)

c     sum contributions over all nodes

      if(mxnode.gt.1)then

        buffer(1)=enginv
        buffer(2)=virinv
        buffer(3)=inv_fre
        buffer(4)=inv_vir
        call gdsum(buffer(1),4,buffer(5))
        enginv=buffer(1)
        virinv=buffer(2)
        inv_fre=buffer(3)
        inv_vir=buffer(4)

        if(lsolva)then

          call gdsum(inv_sol(1),mxtmls,buffer(1))
          if(lexcite)call gdsum(inv_exc(1),mxtmls,buffer(1))

        endif

      endif

      deallocate (xdab,ydab,zdab,stat=fail1)
      deallocate (xdac,ydac,zdac,stat=fail2)
      deallocate (xdad,ydad,zdad,stat=fail3)

      return
      end subroutine invfrc

      end module inversion_module