text:
module fishpack_precision

    use, intrinsic :: ISO_C_binding, only: &
        C_INT, &
        C_DOUBLE, &
        C_DOUBLE_COMPLEX

    ! Explicit typing only
    implicit none

    ! Everything is private unless stated otherwise
    private
    public :: wp ! Working precision
    public :: ip ! Integer precision
    public :: PI, TWO_PI, HALF_PI, EPS, MACHINE_EPSILON
    public :: get_pi

    ! Floating point precision constants
    integer, parameter :: FLOAT128 = selected_real_kind(p=33, r=4931) ! 33 digits, range [10^(-4931), 10^(+4931) - 1], 128 bits
    integer, parameter :: FLOAT64 = selected_real_kind(p=15, r=307) ! 15 digits, range [10^(-307) , 10^(+307)  - 1], 64 bits
    integer, parameter :: FLOAT32 = selected_real_kind(p=6, r=37) ! 6 digits, range [10^(-37)  , 10^(+37)   - 1],  32 bits
    integer, parameter :: wp = C_DOUBLE ! Default floating point precision
    integer, parameter :: cwp = C_DOUBLE_COMPLEX

    ! Integer precision constants
    integer, parameter :: INT64 = selected_int_kind(r=18) ! 19 digits plus sign, range [-2^(63),+2^(63) - 1], 64 bits
    integer, parameter :: INT32 = selected_int_kind(r=9) ! 10 digits plus sign, range [-2^(31),+2^(31) - 1], 32 bits
    integer, parameter :: INT16 = selected_int_kind(r=4) ! 5 digits plus sign, range [-2^(15),+2^(15) - 1], 16 bits
    integer, parameter :: INT8 = selected_int_kind(r=2) ! 3 digits plus sign, range [-2^(7) ,+2^(7)  - 1], 8 bits
    integer, parameter :: ip = C_INT ! Default integer precision

    ! Calculation constants
    real(wp), parameter :: ONE = 1.0_wp
    real(wp), parameter :: TWO = 2.0_wp
    real(wp), parameter :: PI = acos(-ONE)
    real(wp), parameter :: HALF_PI = PI/TWO
    real(wp), parameter :: TWO_PI = TWO * PI
    real(wp), parameter :: MACHINE_EPSILON = epsilon(ONE)
    real(wp), parameter :: EPS = MACHINE_EPSILON

contains

    pure function get_pi() &
        result (return_value)

        ! Dummy arguments
        real(wp) :: return_value

        return_value = 3.141592653589793238462643383279502884197169399375105820974_wp

    end function get_pi

end module fishpack_precision

max_stars_repo_name: jlokimlin/fishpack

max_issues_repo_name: jlokimlin/fishpack


text:
SUBROUTINE qloop(nelbnd,cdbnd,nlbnd)
! delete nodes in the loop of contour lines
IMPLICIT NONE

  !--- Dummy variables
  INTEGER(kind=4),INTENT(INOUT):: nelbnd  !number of segments defining contour
  REAL(kind=8),POINTER:: cdbnd(:,:)       !coordinates of segments defining contour
  INTEGER(kind=4),POINTER:: nlbnd(:)      !segments in each contour lines
  !--- Constants
  REAL(kind=8),PARAMETER:: tol=1d-6
  !--- Local variables
  INTEGER(kind=4):: i, j, l, ii, jj, nl, segbnd, frsts, frstn !elo, elf,
  REAL(kind=8):: aux, long
  REAL(kind=8),ALLOCATABLE:: tt(:,:), nn(:,:)
  LOGICAL,ALLOCATABLE:: ds(:)

  INTERFACE
    INCLUDE 'del_e_mrea.h'
  END INTERFACE

  !=============================
  nl = SIZE(nlbnd,1)
  ALLOCATE(tt(2,nelbnd),nn(2,nelbnd),ds(nelbnd))
  tt = 0d0
  nn = 0d0
  ds = .FALSE.

  !for each line in boundary apply the qloop
  frsts = 0
  frstn = 1
  DO l = 1,nl
  segbnd = frsts + nlbnd(l)  !upper limit for segments in line
  !compute tangent vector and normal vector
  long = 0d0
  DO i=frstn,segbnd
    j = MOD(i,segbnd)+1
    tt(1:2,i) = cdbnd(1:2,j) - cdbnd(1:2,i)
    aux = SQRT(DOT_PRODUCT(tt(1:2,i),tt(1:2,i)))
    long = long + aux
    IF (aux /= 0) tt(1:2,i)=tt(1:2,i)/aux
    nn(1,i) =-tt(2,i)
    nn(2,i) = tt(1,i)
  END DO

  ! compute intersections
  DO i=frstn,segbnd-1
    ii = MOD(i,segbnd)+1  !segment ii
    DO j=i+2,segbnd
      jj = MOD(j,segbnd)+1 !segment jj
      IF (ds(j)) CYCLE
      IF (jj == 1 .AND. i == frstn) CYCLE
      IF (jj == 1) jj = frstn   ! close each line in multiline boundary (WC)
      IF (jj == i) CYCLE  !this ensure that the last segment not compares with first
      IF (DMAX1(cdbnd(1,j),cdbnd(1,jj)) < DMIN1(cdbnd(1,i),cdbnd(1,ii))) CYCLE
      IF (DMAX1(cdbnd(1,i),cdbnd(1,ii)) < DMIN1(cdbnd(1,j),cdbnd(1,jj))) CYCLE
      IF (DMAX1(cdbnd(2,j),cdbnd(2,jj)) < DMIN1(cdbnd(2,i),cdbnd(2,ii))) CYCLE
      IF (DMAX1(cdbnd(2,i),cdbnd(2,ii)) < DMIN1(cdbnd(2,j),cdbnd(2,jj))) CYCLE
      !Calculo de la interseccion
      aux = DOT_PRODUCT(cdbnd(1:2,jj)-cdbnd(1:2,j),nn(1:2,i))
      IF (aux == 0) CYCLE
      aux = DOT_PRODUCT(cdbnd(1:2,i) - cdbnd(1:2,j),nn(1:2,i))/aux
      IF ((aux < -tol) .OR. (aux > 1d0+tol)) CYCLE
!      cd(1:2) = (cdbnd(1:2,jj)-cdbnd(1:2,j))*aux + cdbnd(1:2,j)

      aux = DOT_PRODUCT(cdbnd(1:2,ii)-cdbnd(1:2,i),nn(1:2,j))
      IF (aux == 0) CYCLE
      aux = DOT_PRODUCT(cdbnd(1:2,j) - cdbnd(1:2,i),nn(1:2,j))/aux
      IF ((aux < -tol) .OR. (aux > 1d0+tol)) CYCLE
!      cd(1:2) = ((cdbnd(1:2,ii)-cdbnd(1:2,i))*aux + cdbnd(1:2,i) + cd(1:2))*0.5d0

      ds(i) = .TRUE. ! mark the segments that exhibits
      ds(j) = .TRUE. ! loop intersection
      EXIT
    END DO
  END DO

  frstn = frstn + nlbnd(l)
  frsts = frsts + nlbnd(l)
  END DO

  !if EXIST a loop in boundary STOP the program
  IF (ANY(ds(1:nelbnd))) &
    CALL runen3('QLOOP: BOUNDARY LOOP DETECTED IN REMESHING TASK.')

  !memory release
  DEALLOCATE(tt,nn,ds)

RETURN
END SUBROUTINE  qloop

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
C
C $Id: frame.f,v 1.9 2008-07-27 00:17:24 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE FRAME
      COMMON /GFLASH/MODEF,IOPWKS(100),IOACT(100),NUMOP,IWISSI
C
C  FRAME is designed to effect a break in the picture drawing
C  sequence depending upon whether the workstation type is
C  MO, or whether it is an OUTPUT or OUTIN workstation.
C
C  An UPDATE WORKSTATION and CLEAR WORKSTATION is done on all
C  metafiles and all workstations of type OUTPUT.  For metafiles
C  this inserts an END PICTURE into the metafile.
C
C  If there are any OUTIN workstations, all of them are updated
C  with an UPDATE WORKSTATION and a pause is done on the OUTIN
C  workstation of most recent creation.  After return from the
C  pause, a CLEAR WORKSTATION is done on all OUTIN workstations.
C
      INTEGER WKID
      CHARACTER*80 DATREC,STR,ISTR
C
C Check for an uncleared prior error.
C
      IF (ICFELL('FRAME - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C  First, flush the pen-move buffer.
C
      CALL PLOTIF (0.,0.,2)
      IF (ICFELL('FRAME',2).NE.0) RETURN
C
C  If no workstations are open, return.
C
      CALL GQOPWK (1,IE,NO,ID)
      IF (IE.NE.0) THEN
        CALL SETER ('FRAME - ERROR EXIT FROM GQOPWK',3,1)
        RETURN
      END IF
      IF (NO .EQ. 0) RETURN
C
C  Update all workstations.
C
      DO 200 I=1,NO
C
C  Get the workstation ID.
C
        CALL GQOPWK (I,IE,NO,WKID)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQOPWK',4,1)
          RETURN
        END IF
C
C  Get workstation type.
C
        CALL GQWKC (WKID,IE,ICON,ITYPE)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKC',5,1)
          RETURN
        END IF
C
C  Get workstation category (0=output; 2=out/in; 4=metafile).
C
        CALL GQWKCA (ITYPE,IE,ICAT)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKCA',6,1)
          RETURN
        END IF
C
        IF (ICAT .EQ. 4) THEN
C
C  Illegal to call FRAME while a FLASH buffer is open.
C
          IF (MODEF .EQ. 1) THEN
            CALL SETER
     -    ('FRAME - ILLEGAL TO CALL FRAME WHILE A FLASH BUFFER IS OPEN',
     -      7,1)
          ENDIF
          CALL GCLRWK(WKID,1)
        ELSE IF (ICAT.EQ.0 .OR. ICAT.EQ.2) THEN
          CALL GUWK(WKID,0)
          IF (ICAT .EQ. 0) THEN
            CALL GCLRWK(WKID,1)
          ENDIF
        ENDIF
  200 CONTINUE
C
C  Pause on the OUTIN workstation of most recent creation.
C
      DO 100 I=NO,1,-1
        CALL GQOPWK (I,IE,NO,WKID)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQOPWK',8,1)
          RETURN
        END IF
        CALL GQWKC (WKID,IE,ICON,ITYPE)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKC',9,1)
          RETURN
        END IF
        CALL GQWKCA (ITYPE,IE,ICAT)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKCA',10,1)
          RETURN
        END IF
        IF (ICAT.EQ.2) THEN
          ISTR(1:1) = CHAR(0)
          CALL GINST(WKID,1,0,ISTR,1,0.,1279.,0.,1023.,1,1,1,DATREC)
          CALL GSSTM(WKID,1,0,0)
          CALL GRQST(WKID,1,ISTAT,LOSTR,STR)
          GO TO 110
        ENDIF
  100 CONTINUE
  110 CONTINUE
C
C  Clear all OUTIN workstations.
C
      DO 300 I=1,NO
        CALL GQOPWK (I,IE,NO,WKID)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQOPWK',11,1)
          RETURN
        END IF
        CALL GQWKC (WKID,IE,ICON,ITYPE)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKC',12,1)
          RETURN
        END IF
        CALL GQWKCA (ITYPE,IE,ICAT)
        IF (IE.NE.0) THEN
          CALL SETER ('FRAME - ERROR EXIT FROM GQWKCA',13,1)
          RETURN
        END IF
        IF (ICAT.EQ.2) THEN
          CALL GCLRWK(WKID,1)
        ENDIF
  300 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      double complex function TWbbpZab(p1,p2,p3,p4,p5,p6,p7,p8)
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      include 'sprods_com.f'
C     Author: R.K. Ellis Feb, 2013
C     written by program WbbZdiags.frm
C     These are the Abelian diagrams with both the W and Z coming
C     off light line, only three diagrams with Z emitted before W
C     Calculation is performed for LH light-line (perforce because of W)
C     Calculation is performed for RH bbbar-line
C     Calculation is performed for LH Z-dcay line
      integer p1,p2,p3,p4,p5,p6,p7,p8
      double precision s3,s4,s34,s56,s134,s234,s567,s568,s5678
      double complex d9p,d10p
C     statement functions
      s3(p1,p2,p3)=s(p1,p2)+s(p1,p3)+s(p2,p3)
      s4(p1,p2,p3,p4)=s(p1,p2)+s(p1,p3)+s(p1,p4)
     &               +s(p2,p3)+s(p2,p4)+s(p3,p4)
C     end statement functions
      s34=s(p3,p4)
      s56=s(p5,p6)
      s134=s3(p1,p3,p4)
      s234=s3(p2,p3,p4)
      s567=s3(p5,p6,p7)
      s568=s3(p5,p6,p8)
      s5678=s4(p5,p6,p7,p8)
      d9p= + s34**(-1)*s56**(-1)*s134**(-1)*s567**(-1) * ( za(p1,p3)*
     &    za(p2,p8)*za(p5,p6)*zb(p1,p4)*zb(p1,p6)*zb(p6,p7)*s5678**(-1)
     &     + za(p1,p3)*za(p2,p8)*za(p5,p7)*zb(p1,p4)*zb(p1,p7)*zb(p6,p7
     &    )*s5678**(-1) - za(p2,p8)*za(p3,p4)*za(p5,p6)*zb(p1,p4)*zb(p4
     &    ,p6)*zb(p6,p7)*s5678**(-1) - za(p2,p8)*za(p3,p4)*za(p5,p7)*
     &    zb(p1,p4)*zb(p4,p7)*zb(p6,p7)*s5678**(-1) )
      d9p = d9p + s34**(-1)*s56**(-1)*s134**(-1)*s568**(-1) * ( za(p1,
     &    p3)*za(p2,p5)*za(p5,p8)*zb(p1,p4)*zb(p1,p7)*zb(p5,p6)*
     &    s5678**(-1) - za(p1,p3)*za(p2,p8)*za(p5,p8)*zb(p1,p4)*zb(p1,
     &    p7)*zb(p6,p8)*s5678**(-1) - za(p2,p5)*za(p3,p4)*za(p5,p8)*zb(
     &    p1,p4)*zb(p4,p7)*zb(p5,p6)*s5678**(-1) + za(p2,p8)*za(p3,p4)*
     &    za(p5,p8)*zb(p1,p4)*zb(p4,p7)*zb(p6,p8)*s5678**(-1) )

      d10p= + s34**(-1)*s56**(-1)*s234**(-1)*s567**(-1) * (  - za(p2,p3
     &    )*za(p2,p8)*za(p5,p6)*zb(p1,p6)*zb(p2,p4)*zb(p6,p7)*
     &    s5678**(-1) - za(p2,p3)*za(p2,p8)*za(p5,p7)*zb(p1,p7)*zb(p2,
     &    p4)*zb(p6,p7)*s5678**(-1) - za(p2,p3)*za(p3,p8)*za(p5,p6)*zb(
     &    p1,p6)*zb(p3,p4)*zb(p6,p7)*s5678**(-1) - za(p2,p3)*za(p3,p8)*
     &    za(p5,p7)*zb(p1,p7)*zb(p3,p4)*zb(p6,p7)*s5678**(-1) )
      d10p = d10p + s34**(-1)*s56**(-1)*s234**(-1)*s568**(-1) * (  -
     &    za(p2,p3)*za(p2,p5)*za(p5,p8)*zb(p1,p7)*zb(p2,p4)*zb(p5,p6)*
     &    s5678**(-1) + za(p2,p3)*za(p2,p8)*za(p5,p8)*zb(p1,p7)*zb(p2,
     &    p4)*zb(p6,p8)*s5678**(-1) - za(p2,p3)*za(p3,p5)*za(p5,p8)*zb(
     &    p1,p7)*zb(p3,p4)*zb(p5,p6)*s5678**(-1) + za(p2,p3)*za(p3,p8)*
     &    za(p5,p8)*zb(p1,p7)*zb(p3,p4)*zb(p6,p8)*s5678**(-1) )

      TWbbpZab=d9p+d10p
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


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 IPModel_ASAP
!X
!% Interface to ASAP potential.
!% P. Tangney and S. Scandolo,
!% An ab initio parametrized interatomic force field for silica
!% J. Chem. Phys, 117, 8898 (2002).
!%
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"

module IPModel_ASAP_module

use libAtoms_module
use QUIP_Common_module

implicit none
private

include 'IPModel_interface.h'

logical, private :: asap_initialised = .false.

public :: IPModel_ASAP
type IPModel_ASAP
  integer :: n_types = 0, n_atoms = 0
  real(dp) :: betapol, tolpol, yukalpha, yuksmoothlength
  integer :: maxipol, pred_order
  integer, allocatable :: atomic_num(:), type_of_atomic_num(:)
  real(dp), allocatable, dimension(:) :: pol, z
  real(dp), allocatable, dimension(:,:) :: D_ms, gamma_ms, R_ms, B_pol, C_pol
  logical :: tewald, tdip_sr
  real :: raggio, a_ew, gcut
  integer :: iesr(3)

  real(dp) :: cutoff_coulomb, cutoff_ms

  character(len=STRING_LENGTH) :: label
  logical :: initialised

end type IPModel_ASAP

logical, private :: parse_in_ip, parse_matched_label
type(IPModel_ASAP), private, pointer :: parse_ip

interface Initialise
  module procedure IPModel_ASAP_Initialise_str
end interface Initialise

interface Finalise
  module procedure IPModel_ASAP_Finalise
end interface Finalise

interface Print
  module procedure IPModel_ASAP_Print
end interface Print

interface Calc
  module procedure IPModel_ASAP_Calc
end interface Calc

contains

#ifdef HAVE_ASAP
  subroutine asap_singlepoint_finalise()

    use atoms
    use verlet
    use neighbour
    use pot_parameters
    use stuff
    use energy_i
    use print65
    use electric_field
    use nndim
    use neighbour3
    use fixpar
    use distortion, only : taimsp,rmin_aim
    use thermal_cond
    use shakes
    use compress
    use quench
    use forcetest
    use iteration
    use changepress
    use metric
    use netquant
    use universal
    use initstuff
    use plot_efield
    use polar
    use structopt
    use parameters
    use presstemp
    use testf
    use minimiser
    implicit none

    deallocate(mass,spind,objind,numobj,numobjsp)
    deallocate(nnlist,nnlist2,nnlist3)
    deallocate(nnat,nnat2,nnat3,nn_imag)
    deallocate(nn_imag2,nn_imag3)
    deallocate(s,sm)
    deallocate(r,rm)
    deallocate(Z,bij,Cij)
    deallocate(Dij,alphaij)
    deallocate(Eij,Nij)
    deallocate(pol,bpol,cpol)
    deallocate(theta0jik,bjik)
    deallocate(b_tt,r_ms)
    deallocate(gamma_ms,d_ms)
    deallocate(c_harm,rmin,rmin_aim)
    deallocate(adist,bdist,cdist,ddist)
    deallocate(bu1,alphau1)
    deallocate(sigma1,sigma2,sigma3,sigma4)
    deallocate(bu2,alphau2)
    deallocate(taimsp,minr,mindistl)
    deallocate(elements)
    deallocate(aelements)
    deallocate(ielements)
    deallocate(dt2bym, aumass)
    deallocate(parf, tparf)
    deallocate(posobj,velobj,dipobj)
    deallocate(timposepos,timposevel,timposedip)
    deallocate(zerovec,falseimp)
    if (allocated(dip_perm)) deallocate(dip_perm)
    if (t_therm) deallocate(nzave,ztempave)

    if (allocated(efield_old)) deallocate(efield_old)
    if (allocated(efield)) deallocate(efield)
    if (allocated(fqdip)) deallocate(fqdip)
    if (allocated(fdipdip)) deallocate(fdipdip)
    if (allocated(force_ew)) deallocate(force_ew)
    if (allocated(efield_ew)) deallocate(efield_ew)
    if (allocated(dip)) deallocate(dip)
    if (allocated(stress_ew)) deallocate(stress_ew)
    if (allocated(dip_sr)) deallocate(dip_sr)
    if (allocated(efield_old)) deallocate(efield_old)
    if (allocated(force_pol)) deallocate(force_pol)
    if (allocated(stress_pol)) deallocate(stress_pol)
    if (allocated(testt)) deallocate(testt)

  end subroutine asap_singlepoint_finalise

#endif

subroutine IPModel_ASAP_Initialise_str(this, args_str, param_str)
  use system_module, only : dp
  type(IPModel_ASAP), intent(inout) :: this
  character(len=*), intent(in) :: args_str, param_str

  type(Dictionary) :: params

  this%initialised = .false.

  call Finalise(this)

  call initialise(params)
  this%label=''
  call param_register(params, 'label', '', this%label, help_string="No help yet.  This source file was $LastChangedBy$")
  if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='IPModel_ASAP_Initialise_str args_str')) then
    call system_abort("IPModel_ASAP_Initialise_str failed to parse label from args_str="//trim(args_str))
  endif
  call finalise(params)

  call IPModel_ASAP_read_params_xml(this, param_str)
  this%n_atoms = 0
  this%initialised = .true.

end subroutine IPModel_ASAP_Initialise_str

subroutine IPModel_ASAP_Finalise(this)
  use system_module, only : dp
  type(IPModel_ASAP), intent(inout) :: this

#ifdef HAVE_ASAP
  if (asap_initialised) then
     call asap_singlepoint_finalise()
     asap_initialised = .false.
  end if
  this%initialised = .false.
#else
  call system_abort('ASAP potential is not compiled in. Recompile with HAVE_ASAP=1')
#endif

  if (allocated(this%atomic_num)) deallocate(this%atomic_num)
  if (allocated(this%type_of_atomic_num)) deallocate(this%type_of_atomic_num)
  if (allocated(this%pol)) deallocate(this%pol)
  if (allocated(this%z)) deallocate(this%z)
  if (allocated(this%D_ms)) deallocate(this%D_ms)
  if (allocated(this%gamma_ms)) deallocate(this%gamma_ms)
  if (allocated(this%R_ms)) deallocate(this%R_ms)
  if (allocated(this%B_pol)) deallocate(this%B_pol)
  if (allocated(this%C_pol)) deallocate(this%C_pol)
  this%n_types = 0
  this%label = ''
end subroutine IPModel_ASAP_Finalise

subroutine IPModel_ASAP_Calc(this, at, e, local_e, f, virial, local_virial, args_str, mpi, error)
#ifdef HAVE_ASAP
  use neighbour!, only : nsp,nat,iesr,spind,rcut
  use neighbour3
  use pot_parameters
  use logical_stuff
  use polar
  use stuff
  use distortion, only :  xgmin,xgmax,taimsp
  use parameters
  use testf
  use initstuff
  use atoms
  use verlet
  use iteration
  use nndim
  use metric
  use fixpar, only: ntype, tpbc, treadpar
  use energy_i
  use print65
  use electric_field
  use nndim
  use distortion, only : taimsp,rmin_aim
  use thermal_cond
  use shakes
  use compress
  use quench
  use forcetest
  use iteration
  use changepress
  use metric
  use netquant
  use universal, only: au_to_kbar, small
  use initstuff
  use plot_efield
  use polar
  use structopt
  use parameters
  use presstemp
  use testf
  use minimiser

#endif
   use libAtoms_module, only : myAtoms => Atoms, Dictionary, BOHR, HARTREE, &
        ElementMass, massconvert, initialise, param_read_line, finalise, &
        ElementName, print, operator(//), cell_volume, has_property, &
        assign_pointer, add_property, system_abort, fit_box_in_cell, &
        param_register, check_size, operator(.fne.), assign_property_pointer

   type(IPModel_ASAP), intent(inout):: this
   type(myAtoms), intent(inout)      :: at
   real(dp), intent(out), optional :: e, local_e(:)
   real(dp), intent(out), optional :: f(:,:), local_virial(:,:)   !% Forces, dimensioned as \texttt{f(3,at%N)}, local virials, dimensioned as \texttt{local_virial(9,at%N)}
   real(dp), intent(out), optional :: virial(3,3)
   character(len=*), optional, intent(in) :: args_str
   type(MPI_Context), intent(in), optional :: mpi

   integer, intent(out), optional :: error
   type(Dictionary) :: params
   real(dp) :: asap_e, asap_stress(3,3), a_ew, f_ew
   logical smlra,smlrc,smlgc,smlrc2,smlrc3,smlrc4
   logical nsmlra,nsmlrc,nsmlgc
   real*8 rarec,gcrec,rcrec,root2
   real*8 raggio_in,rcut_in,gcut_in
   real(dp), allocatable :: asap_f(:,:)
   real(dp), pointer :: dipoles_ptr(:,:), ext_efield_ptr(:,:)
   integer :: i, ti, tj
   logical :: do_restart, calc_dipoles
   integer at0, atf
   integer idebug
   real(dp) dtold,dtnew
   logical  tzeroc
   logical sumewald
   real(dp) mass_cel
   logical tscaled
   logical readnat
   logical texist
   logical tpow,tgmin
   integer nesr

   logical :: has_atom_mask_name, applied_efield
   character(STRING_LENGTH) :: atom_mask_name
   real(dp) :: r_scale, E_scale
   logical :: do_rescale_r, do_rescale_E

#ifdef HAVE_ASAP

   INIT_ERROR(error)

   if (present(e)) e = 0.0_dp
   if (present(local_e)) then
      call check_size('Local_E',local_e,(/at%N/),'IPModel_ASAP_Calc', error)
      local_e = 0.0_dp
   endif
   if (present(f)) then
      call check_size('Force',f,(/3,at%Nbuffer/),'IPModel_ASAP_Calc', error)
      f = 0.0_dp
   end if
   if (present(virial)) virial = 0.0_dp
   if (present(local_virial)) then
      call check_size('Local_virial',local_virial,(/9,at%Nbuffer/),'IPModel_ASAP_Calc', error)
      local_virial = 0.0_dp
      RAISE_ERROR("IPModel_ASAP_Calc: local_virial calculation requested but not supported yet.", error)
   endif

   allocate(asap_f(3,at%N))

   call initialise(params)
   this%label=''
   call param_register(params, 'restart', 'F', do_restart, help_string="No help yet.  This source file was $LastChangedBy$")
   call param_register(params, 'calc_dipoles', 'T', calc_dipoles, help_string="No help yet.  This source file was $LastChangedBy$")
   call param_register(params, 'atom_mask_name', 'NONE', atom_mask_name, has_value_target=has_atom_mask_name, help_string="No help yet.  This source file was $LastChangedBy$")
   call param_register(params, 'r_scale', '1.0',r_scale, has_value_target=do_rescale_r, help_string="Recaling factor for distances. Default 1.0.")
   call param_register(params, 'E_scale', '1.0',E_scale, has_value_target=do_rescale_E, help_string="Recaling factor for energy. Default 1.0.")
   call param_register(params, 'applied_efield', 'F', applied_efield, help_string="No help yet.  This source file was $LastChangedBy$")

   if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='IPModel_ASAP_Calc args_str')) then
      call system_abort("IPModel_ASAP_Initialise_str failed to parse args_str="//trim(args_str))
   endif
   call finalise(params)

   if(has_atom_mask_name) then
      RAISE_ERROR('IPModel_ASAP_Calc: atom_mask_name found, but not supported', error)
   endif
   if (do_rescale_r .or. do_rescale_E) then
      RAISE_ERROR("IPModel_ASAP_Calc: rescaling of potential with r_scale and E_scale not yet implemented!", error)
   end if

   if (.not. asap_initialised .or. this%n_atoms /= at%n) then

      if (asap_initialised) then
         call asap_singlepoint_finalise()
         asap_initialised = .false.
      end if

      this%n_atoms = at%n

      idebug = 0
      nat = this%n_atoms
      nsp = this%n_types
      at0 = 1
      atf = nat
      nobj = 1

      nthermo = min(atf-at0+1,nat-1)
      allocate(mass(this%n_types))

      tsinglepoint = .true.
      tangstrom = .false.
      irestart = -1
      itsave = 0
      readunit = 90
      saveunit = 90
      readnat = .false.
      dtold = 0.0d0
      dtnew = 0.0d0
      ntstep = 0
      iprint = 0
      iprintxyz = 0
      deltat = dtnew
      tcenter = .false.
      timposepbc = .false.
      tremovetrans = .false.
      treadxyz = .false.
      tfirst_xyzprint = .true.
      treadpar = .false.
      tappendfiles = .false.
      tstructopt = .false.
      tsdp = .false.
      tpbc = .true.
      tpow = .false.
      tgmin = .false.
      testewald = .false.
      ttime = .true.
      tforcetest = .false.
      tscaled = .false.
      tangstrom = .false.
      texist = .false.
      tprint65 = .false.
      tenergy_i = .false.
      tbin65 = .false.
      iprint65 = 1
      tzeroc = .false.
      tcel = .false.
      tsdc = .false.
      cell_dir = 3
      press = 0.0d0
      mass_cel = 2.0d6
      press = press / au_to_kbar
      wc = mass_cel
      trescale = .false.
      tnosep = .false.
      tboltz = .false.
      tshake = .false.
      tcompress = .false.
      tquench = .false.
      tchangeP= .false.
      tpolwrite = .false.
      t_therm = .false.
      tshowforce = .false.
      tplot_efield = .false.
      tfirst_efield_plot = .true.
      tewald = .false.
      write_sr = .false.
      read_sr= .false.
      calc_sr=.false.
      calc_pol = .false.
      write_pol=.false.
      read_pol=.false.
      calc_gvec=.false.
      write_gvec=.false.
      read_gvec=.false.
      write_aim = .false.
      read_aim= .false.
      calc_aim=.false.
      tzvar=.false.
      tpolvar=.false.
      tsrvar=.false.
      tpol=.false.
      tsr=.false.
      taim = .false.
      tz=.false.
      hafta=.false.
      tbegin = .true.
      tyukawa = .true.

      sumewald = .true.
      allocate(spind(nat),objind(nat),numobj(nobj),numobjsp(nobj,nsp))
      allocate(z(nsp),bij(nsp,nsp),Cij(nsp,nsp))
      allocate(Dij(nsp,nsp),alphaij(nsp,nsp))
      allocate(Eij(nsp,nsp),Nij(nsp,nsp))
      allocate(pol(nsp),bpol(nsp,nsp),cpol(nsp,nsp))
      allocate(theta0jik(nsp,nsp,nsp),bjik(nsp,nsp,nsp))
      allocate(b_tt(3,nsp,nsp),r_ms(nsp,nsp))
      allocate(gamma_ms(nsp,nsp),d_ms(nsp,nsp))
      allocate(c_harm(nsp,nsp),rmin(nsp,nsp))
      allocate(rmin_aim(nsp,nsp))
      rmin_aim = 0.d0
      allocate(adist(nsp,nsp),bdist(nsp,nsp))
      allocate(cdist(nsp),ddist(nsp,nsp))
      allocate(bu1(nsp,nsp),alphau1(nsp,nsp))
      allocate(sigma1(nsp),sigma2(nsp),sigma3(nsp),sigma4(nsp))
      allocate(bu2(nsp,nsp),alphau2(nsp,nsp))
      allocate(taimsp(nsp))
      allocate(minr(nsp,nsp),mindistl(nsp,nsp))
      mindistl = 0.d0

      allocate(s(3,nat),sm(3,nat))
      allocate(r(3,nat),rm(3,nat))
      allocate(dt2bym(nat),aumass(nat))
      allocate(elements(nsp))
      allocate(aelements(nat))
      allocate(ielements(nat))
      nparm = (nsp*nsp*nsp+21*nsp*nsp+28*nsp)/2 + 4*nsp
      allocate (parf(nparm),tparf(nparm))
      tgen = .false.
      testforce = .false.
      npar = 1

      ntnlist = (/1,1/)
      nnatmax = 100000
      nnatmax_sr = 30000
      allocate(posobj(3,nobj),velobj(3,nobj),dipobj(3,nobj))
      allocate(timposepos(nobj),timposevel(nobj),timposedip(nobj))
      allocate(zerovec(3,nobj),falseimp(nobj))
      zerovec = 0.0d0!

      falseimp = .false.
      XNOS2M = 0.D0
      XNOSM  = 0.D0
      XNOS0  = 0.D0
      VRNOS  = 0.D0

      it = 1

      c_harm = 0.0_dp
      rmin = 0.0_dp

      ! Per type parameters
      do ti=1,this%n_types
         mass(ti)     = ElementMass(this%atomic_num(ti))/MASSCONVERT
         elements(ti) = ElementName(this%atomic_num(ti))
         z(ti) = this%z(ti)
         pol(ti) = this%pol(ti)
         ntype(ti) = count(at%Z == this%atomic_num(ti))
      end do

      write(6,'(/," Nos of each species:",10i6)') (ntype(i),i=1,nsp)
      write(6,'(/," Masses :",2(1x,f10.5))') mass

      ! Ewald parameters
      tewald = this%tewald
      raggio = this%raggio
      a_ew = this%a_ew
      gcut = this%gcut
      rcut(1) = this%cutoff_coulomb
      rcut(3) = this%cutoff_ms
      iesr = this%iesr

      call print('tewald :'//tewald)
      write(6,'(/," Raggio :",1x,f10.5)') raggio
      write(6,'(" A_ew   :",1x,e10.5)') a_ew
      write(6,'(" iesr   :",3i3)') iesr
      write(6,'(" Gcut   :",1x,f10.5,/)') gcut

      smlgc = gcut.lt.small
      smlrc = rcut(1).lt.small
      smlra = raggio.lt.small
      smlrc2 = rcut(2).lt.small
      smlrc3 = rcut(3).lt.small
      smlrc4 = rcut(4).lt.small

      nsmlgc = .not.smlgc      ! iesr will be computed later from rcut and cell

      nsmlrc = .not.smlrc
      nsmlra = .not.smlra

      f_ew = dsqrt(-1.0d0*log(a_ew))
      root2 = dsqrt(2.0d0)

      rarec = 0.0d0
      rcrec = 0.0d0
      gcrec = 0.0d0

      raggio_in = raggio
      gcut_in = gcut
      rcut_in = rcut(1)

      if (smlrc.and.smlgc.and.smlra) then
       rarec = 3.0d0
       gcrec = f_ew*root2/rarec
       rcrec = gcut*rarec*rarec
      else if (nsmlrc.and.smlgc.and.smlra) then
       rarec = rcut(1)/f_ew/root2
       gcrec = f_ew*root2/rarec
       raggio = rarec
       gcut = gcrec
      else if (smlrc.and.nsmlgc.and.smlra) then
       rarec  = f_ew*root2/gcut
       rcrec  = gcut*rarec*rarec
       raggio = rarec
       rcut(1) = rcrec
      else if (smlrc.and.smlgc.and.nsmlra) then
       gcrec = f_ew*dsqrt(2.0d0)/raggio
       rcrec = gcrec*raggio*raggio
       gcut = gcrec
       rcut(1) = rcrec
      else if (nsmlrc.and.nsmlgc.and.smlra) then
       rarec = rcut(1)/f_ew/root2
       gcrec = f_ew*root2/raggio
       raggio = rarec
      else if (nsmlrc.and.smlgc.and.nsmlra) then
       rarec = rcut(1)/f_ew/root2
       gcrec = f_ew*root2/raggio
       gcut = gcrec
      else if (smlrc.and.nsmlgc.and.nsmlra) then
       gcrec = f_ew*root2/raggio
       rcrec = gcrec*raggio*raggio
       rcut(1) = rcrec
      else if (nsmlrc.and.nsmlgc.and.nsmlra) then
       write(6,*) 'raggio, rcut, gcut, are all defined!!'
      endif

      if (smlrc2) rcut(2) = rcut(1)
      if (smlrc3) rcut(3) = rcut(1)
      if (smlrc4) rcut(4) = rcut(1)

      write(6,'(/," Input Raggio                : ",f10.5)')raggio_in
      write(6,'(" Input r-space cutoff        : ",f10.5)')rcut_in
      write(6,'(" Input g-space cutoff        : ",f10.5)')gcut_in
      write(6,'(/," Recommended Raggio          : ",f10.5)')rarec
      write(6,'(" Recommended r-space cutoff  : ",f10.5)')rcrec
      write(6,'(" Recommended g-space cutoff  : ",f10.5)')gcrec
      write(6,'(/," Using Raggio                : ",f10.5)')raggio
      write(6,'(" Using r-space cutoff        : ",f10.5)')rcut(1)
      write(6,'(" Using g-space cutoff        : ",f10.5)')gcut
      write(6,'(/," Secondary r-space cutoff    : ",f10.5)')rcut(2)
      write(6,'(" Tertiary  r-space cutoff    : ",f10.5)')rcut(3)
      write(6,'(" Quaternary r-space cutoff   : ",f10.5)')rcut(4)

      netcharge = 0.0d0
      do i=1,this%n_types
         netcharge = netcharge + z(i)*dfloat(ntype(i))
         if (dabs(z(i)).gt.1.d-10) tz = .true.
      end do

      write(6,'(/," Charges : ",100(e13.5))') z
      write(6,'(/," Net charge ",e13.5)') netcharge
      write(6,*)

      ! Short range parameters
      tsr = .true.
      alphaij = 0.0_dp
      bij = 0.0_dp
      cij = 0.0_dp
      dij = 0.0_dp
      eij = 0.0_dp
      nij = 0.0_dp
      b_tt = 0.0_dp
      gamma_ms = 0.d0
      r_ms = 0.0d0
      d_ms= 0.0d0
      do ti=1,this%n_types
         do tj=1,this%n_types
            d_ms(ti,tj) = this%d_ms(ti,tj)
            gamma_ms(ti,tj) = this%gamma_ms(ti,tj)
            R_ms(ti,tj) = this%R_ms(ti,tj)
         end do
      end do

      ! Polarisation
      betapol = this%betapol
      maxipol = this%maxipol
      tolpol = this%tolpol
      pred_order = this%pred_order

      tpol = any(dabs(pol) > 1.0e-6_dp)
      tpolvar = tpol
      call Print('tpol = '//tpol)
      call Print('tpolvar = '//tpolvar)
      call Print('Polarisation: '//pol)

      bpol = 0.0_dp
      cpol = 0.0_dp
      do ti=1,this%n_types
         do tj=1,this%n_types
            bpol(ti,tj) = this%B_pol(ti,tj)
            cpol(ti,tj) = this%C_pol(ti, tj)
         end do
      end do

      ! Smoothing
      smooth = .false.

      ! Yukawa parameters
      yukalpha = this%yukalpha
      yuksmoothlength = this%yuksmoothlength
      tdip_sr = this%tdip_sr

      asap_initialised = .true.
   end if

   spind = 0
   do i=1,this%n_types
      where(at%Z == this%atomic_num(i)) spind = i
   end do

   ! ASAP uses atomic units - lengths are in Bohr, energies in Hartree,
   ! forces in Hartree/Bohr and stress in Hartree/Bohr**3

   r(:,:) = at%pos(:,:)/BOHR         ! positions
   rm(:,:) = at%pos(:,:)/BOHR
   htm = transpose(at%lattice/BOHR)  ! lattice
   ht = transpose(at%lattice/BOHR)
   restart = do_restart   ! reset electric field

   if (applied_efield) then
      call assign_property_pointer(at, 'ext_efield', ext_efield_ptr, error)
      PASS_ERROR(error)

      ! check electric field is constant
      if (maxval(abs(sum(ext_efield_ptr,2)/real(size(ext_efield_ptr,2),dp) - ext_efield_ptr(:,1))) > 1.0e-6_dp) then
         RAISE_ERROR('IPModel_ASAP_calc: external efield must be constant (same for all atoms)', error)
      end if

      telectric = .true.
      tfieldion = .true.
      e_applied = ext_efield_ptr(:,1)/(HARTREE/BOHR) ! external applied electric field
   else
      telectric = .false.
      tfieldion = .false.
      e_applied(:) = 0.0_dp
   end if

   objind = 1
   numobj(1) = nat
   ntype = 1
   numobjsp = 0
   do i=1,nat
      ntype(spind(i)) = ntype(spind(i)) + 1
      numobjsp(1,spind(i)) = numobjsp(1,spind(i)) + 1
   end do

   call inv3(ht,htm1,omega)

   if (all(iesr == -1)) then
      ! given rcut(1) and lattice, compute iesr
      call fit_box_in_cell(rcut(1),rcut(1),rcut(1), at%lattice, iesr(1), iesr(2), iesr(3))
      iesr = iesr/2
   end if

   nesr = (2*iesr(1)+1)*(2*iesr(2)+1)*(2*iesr(3)+1)
   nnatmax = min(nnatmax,nat*nesr)
   nnatmax = nnatmax + 10

   if (allocated(nnat)) deallocate(nnat)
   if (allocated(nnat2)) deallocate(nnat2)
   if (allocated(nnat3)) deallocate(nnat3)
   if (allocated(nnlist)) deallocate(nnlist)
   if (allocated(nnlist2)) deallocate(nnlist2)
   if (allocated(nnlist3)) deallocate(nnlist3)
   if (allocated(nn_imag)) deallocate(nn_imag)
   if (allocated(nn_imag2)) deallocate(nn_imag2)
   if (allocated(nn_imag3)) deallocate(nn_imag3)

   call label_elements
   call init
   call prepare_traj

   if((mod(it,ntnlist(1)).eq.1).or.(ntnlist(1).eq.1)) then
      call nbrlist(r,nnatmax)
      call nbrlist3(r)
   endif

   if(mod(it,ntnlist(2)).eq.1) then
      call nbrlist3(r)
   endif

   call force_ft(asap_e, asap_f, asap_stress, 0)
   it = it + 1

   ! Convert to {eV, A, fs} units
   if (present(e)) e = asap_e*HARTREE
   if (present(f)) f = asap_f*(HARTREE/BOHR)
   if (present(virial)) virial = asap_stress*(HARTREE/(BOHR**3))*cell_volume(at)
   if (present(local_e)) local_e = 0.0_dp

   if (calc_dipoles .and. tpol) then
      if (.not. has_property(at, 'dipoles')) call add_property(at, 'dipoles', 0.0_dp, n_cols=3)
      call assign_property_pointer(at, 'dipoles', dipoles_ptr, error)
      PASS_ERROR(error)
      dipoles_ptr = dip*BOHR
   end if

   deallocate(asap_f)
#else
  RAISE_ERROR('ASAP potential is not compiled in. Recompile with HAVE_ASAP=1', error)
#endif

end subroutine IPModel_ASAP_Calc

subroutine IPModel_ASAP_Print(this, file)
  use system_module, only : dp
  type(IPModel_ASAP), intent(in) :: this
  type(Inoutput), intent(inout),optional :: file

  integer :: ti, tj

  call Print("IPModel_ASAP : ASAP Potential", file=file)
  call Print("IPModel_ASAP : n_types = " // this%n_types //" n_atoms = "//this%n_atoms, file=file)
  call Print("IPModel_ASAP : betapol = "//this%betapol//" maxipol = "//this%maxipol//" tolpol = "//this%tolpol//" pred_order = "//this%pred_order, file=file)
  call Print("IPModel_ASAP : yukalpha = "//this%yukalpha//" yuksmoothlength = "//this%yuksmoothlength, file=file)
  call Print("IPModel_ASAP : cutoff_coulomb = "//this%cutoff_coulomb//" cutoff_ms = "//this%cutoff_ms, file=file)

  do ti=1, this%n_types
    call Print ("IPModel_ASAP : type " // ti // " atomic_num " // this%atomic_num(ti), file=file)
    call Print ("IPModel_ASAP : pol = "//this%pol(ti), file=file)
    call Print ("IPModel_ASAP : z   = "//this%z(ti), file=file)
   call verbosity_push_decrement()
    do tj =1,this%n_types
       call Print ("IPModel_ASAP : pair interaction ti tj " // ti // " " // tj // " Zi Zj " // this%atomic_num(ti) //&
            " " // this%atomic_num(tj), file=file)
       call Print ("IPModel_ASAP : pair " // this%D_ms(ti,tj) // " " // this%gamma_ms(ti,tj) // " " &
            // this%R_ms(ti,tj) // " " // this%B_pol(ti,tj) // " " // this%C_pol(ti, tj), file=file)
    end do
   call verbosity_pop()
  end do

end subroutine IPModel_ASAP_Print

subroutine IPModel_ASAP_read_params_xml(this, param_str)
  use system_module, only : dp
  type(IPModel_ASAP), intent(inout), target :: this
  character(len=*), intent(in) :: param_str

  type(xml_t) :: fxml

  if (len(trim(param_str)) <= 0) return

  parse_in_ip = .false.
  parse_matched_label = .false.
  parse_ip => this

  call open_xml_string(fxml, param_str)
  call parse(fxml,  &
    startElement_handler = IPModel_startElement_handler, &
    endElement_handler = IPModel_endElement_handler)
  call close_xml_t(fxml)

  if (this%n_types == 0) then
    call system_abort("IPModel_ASAP_read_params_xml parsed file, but n_types = 0")
  endif

end subroutine IPModel_ASAP_read_params_xml

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!% XML param reader functions
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
subroutine IPModel_startElement_handler(URI, localname, name, attributes)
  use system_module, only : dp
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name
  type(dictionary_t), intent(in) :: attributes

  integer :: status
  character(len=STRING_LENGTH) :: val

  integer ti, tj, Zi, Zj

  if (name == 'TS_params') then ! new ASAP stanza

    if (parse_matched_label) return ! we already found an exact match for this label

    call Quip_fox_get_value(attributes, 'label', val, status)
    if (status /= 0) val = ''

    if (len(trim(parse_ip%label)) > 0) then ! we were passed in a label
      if (val == parse_ip%label) then ! exact match
        parse_matched_label = .true.
        parse_in_ip = .true.
      else ! no match
        parse_in_ip = .false.
      endif
    else ! no label passed in
      parse_in_ip = .true.
    endif

    if (parse_in_ip) then
      if (parse_ip%n_types /= 0) then
        call finalise(parse_ip)
      endif

      call Quip_fox_get_value(attributes, 'n_types', val, status)
      if (status == 0) then
        read (val, *), parse_ip%n_types
      else
        call system_abort("Can't find n_types in TS_params")
      endif

      allocate(parse_ip%atomic_num(parse_ip%n_types))
      parse_ip%atomic_num = 0

      allocate(parse_ip%pol(parse_ip%n_types))
      parse_ip%pol = 0.0_dp
      allocate(parse_ip%z(parse_ip%n_types))

      allocate(parse_ip%D_ms(parse_ip%n_types,parse_ip%n_types))
      parse_ip%D_ms = 0.0_dp
      allocate(parse_ip%gamma_ms(parse_ip%n_types,parse_ip%n_types))
      parse_ip%gamma_ms = 0.0_dp
      allocate(parse_ip%R_ms(parse_ip%n_types,parse_ip%n_types))
      parse_ip%R_ms = 0.0_dp
      allocate(parse_ip%B_pol(parse_ip%n_types,parse_ip%n_types))
      parse_ip%B_pol = 0.0_dp
      allocate(parse_ip%C_pol(parse_ip%n_types,parse_ip%n_types))
      parse_ip%C_pol = 0.0_dp

      call QUIP_FoX_get_value(attributes, "cutoff_coulomb", val, status)
      if (status /= 0) call system_abort ("IPModel_TS_read_params_xml cannot find cutoff_coulomb")
      read (val, *) parse_ip%cutoff_coulomb

      call QUIP_FoX_get_value(attributes, "cutoff_ms", val, status)
      if (status /= 0) call system_abort ("IPModel_TS_read_params_xml cannot find cutoff_ms")
      read (val, *) parse_ip%cutoff_ms

      call Quip_fox_get_value(attributes, "betapol", val, status)
      if (status == 0) read (val, *) parse_ip%betapol

      call Quip_fox_get_value(attributes, "maxipol", val, status)
      if (status == 0) read (val, *) parse_ip%maxipol

      call Quip_fox_get_value(attributes, "tolpol", val, status)
      if (status == 0) read (val, *) parse_ip%tolpol

      call Quip_fox_get_value(attributes, "pred_order", val, status)
      if (status == 0) read (val, *) parse_ip%pred_order

      call Quip_fox_get_value(attributes, "yukalpha", val, status)
      if (status == 0) read (val, *) parse_ip%yukalpha

      call Quip_fox_get_value(attributes, "yuksmoothlength", val, status)
      if (status == 0) read (val, *) parse_ip%yuksmoothlength

      parse_ip%tewald = .false.
      call Quip_fox_get_value(attributes, "tewald", val, status)
      if (status == 0) read (val, *), parse_ip%tewald

      parse_ip%tdip_sr = .true.
      call Quip_fox_get_value(attributes, "tdip_sr", val, status)
      if (status == 0) read (val, *), parse_ip%tdip_sr

      parse_ip%raggio = 0.0_dp
      call Quip_fox_get_value(attributes, "raggio", val, status)
      if (status == 0) read (val, *), parse_ip%raggio

      parse_ip%a_ew = 0.0_dp
      call Quip_fox_get_value(attributes, "a_ew", val, status)
      if (status == 0) read (val, *), parse_ip%a_ew

      parse_ip%gcut = 0.0_dp
      call Quip_fox_get_value(attributes, "gcut", val, status)
      if (status == 0) read (val, *), parse_ip%gcut

      parse_ip%iesr = 0
      call Quip_fox_get_value(attributes, "iesr", val, status)
      if (status == 0) read (val, *), parse_ip%iesr

    endif

  elseif (parse_in_ip .and. name == 'per_type_data') then

    call Quip_fox_get_value(attributes, "type", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find type")
    read (val, *) ti

    call Quip_fox_get_value(attributes, "atomic_num", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find atomic_num")
    read (val, *) parse_ip%atomic_num(ti)

    call Quip_fox_get_value(attributes, "pol", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find pol")
    read (val, *) parse_ip%pol(ti)

    call Quip_fox_get_value(attributes, "z", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find z")
    read (val, *) parse_ip%z(ti)

    if (allocated(parse_ip%type_of_atomic_num)) deallocate(parse_ip%type_of_atomic_num)
    allocate(parse_ip%type_of_atomic_num(maxval(parse_ip%atomic_num)))
    parse_ip%type_of_atomic_num = 0
    do ti=1, parse_ip%n_types
      if (parse_ip%atomic_num(ti) > 0) &
        parse_ip%type_of_atomic_num(parse_ip%atomic_num(ti)) = ti
    end do

  elseif (parse_in_ip .and. name == 'per_pair_data') then

    call Quip_fox_get_value(attributes, "atnum_i", val, status)
    if (status /= 0) call system_abort ("IPModel_SW_read_params_xml cannot find atnum_i")
    read (val, *) Zi
    call Quip_fox_get_value(attributes, "atnum_j", val, status)
    if (status /= 0) call system_abort ("IPModel_SW_read_params_xml cannot find atnum_j")
    read (val, *) Zj

    ti = get_type(parse_ip%type_of_atomic_num,Zi)
    tj = get_type(parse_ip%type_of_atomic_num,Zj)

    call Quip_fox_get_value(attributes, "D_ms", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find D_ms")
    read (val, *) parse_ip%D_ms(ti,tj)
    call Quip_fox_get_value(attributes, "gamma_ms", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find gamma_ms")
    read (val, *) parse_ip%gamma_ms(ti,tj)
    call Quip_fox_get_value(attributes, "R_ms", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find R_ms")
    read (val, *) parse_ip%R_ms(ti,tj)
    call Quip_fox_get_value(attributes, "B_pol", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find B_pol")
    read (val, *) parse_ip%B_pol(ti,tj)
    call Quip_fox_get_value(attributes, "C_pol", val, status)
    if (status /= 0) call system_abort ("IPModel_ASAP_read_params_xml cannot find C_pol")
    read (val, *) parse_ip%C_pol(ti,tj)

    if (ti /= tj) then
      parse_ip%D_ms(tj,ti) = parse_ip%D_ms(ti,tj)
      parse_ip%gamma_ms(tj,ti) = parse_ip%gamma_ms(ti,tj)
      parse_ip%R_ms(tj,ti) = parse_ip%R_ms(ti,tj)
      parse_ip%B_pol(tj,ti) = parse_ip%B_pol(ti,tj)
      parse_ip%C_pol(tj,ti) = parse_ip%C_pol(ti,tj)
    endif

  endif

end subroutine IPModel_startElement_handler

subroutine IPModel_endElement_handler(URI, localname, name)
  use system_module, only : dp
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name

  if (parse_in_ip) then
    if (name == 'TS_params') then
      parse_in_ip = .false.
    end if
  endif

end subroutine IPModel_endElement_handler

end module IPModel_ASAP_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
      subroutine MATL03(eps,trace,td,d,ud,hn,h1,nh,sig,dd,isw)

c      * * F E A P * * A Finite Element Analysis Program

c     Michael Scott
c     CE 233

c-----[--.----+----.----+----.-----------------------------------------]
c     Purpose: Linear Isotropic Hardening Material Model

c     Input:
c          eps(*)  -  Current strains at point
c          trace   -  Trace of strain at point
c          td      -  Temperature change
c          d(*)    -  Program material parameters (ndd)
c          ud(*)   -  User material parameters (nud)
c          hn(nh)  -  History terms at point: t_n
c          h1(nh)  -  History terms at point: t_n+1
c          nh      -  Number of history terms
c          isw     -  Solution option from element

c     Output: (N.B. Arrays are set to zero before call)
c          sig(6)  -  Stresses at point.
c          dd(6,6) -  Current material tangent moduli
c-----[--.----+----.----+----.-----------------------------------------]

      implicit none

      include  'iofile.h'
      include  'tdata.h'
      include  'counts.h'

      integer  nh, isw
      real*8   trace, td
      real*8   eps(6),d(*),ud(4),hn(nh),h1(nh),sig(6),dd(6,6)

      real*8   K,G,sigY,H,twoG
      real*8   epn(9),en1(9),syn,kpn,kpn1
      real*8   pn1,sn1(9)
      real*8   one3,two3,root23

      real*8   evn1,nsn1,Fn1,nn1(9),dlam

      integer  i,j
      real*8   tmp1,tmp2,a(6)

c     // Compute useful constants
      one3 = 1.d0/3.d0
      two3 = 2.d0*one3
      root23 = sqrt(two3)

c     // Read in material parameters
      K = ud(1)
      G = ud(2)
      sigY = ud(3)
      H = ud(4)
      twoG = 2.d0*G

c     // Add pressure contribution to material stress
c     // \sigma += m p
      pn1 = K*trace
      do i = 1,3
         sig(i) = pn1
      end do

c     // Add bulk term to material tangent
c     // D_t += K * mm^T
      do i = 1,3
         do j = 1,i
            dd(i,j) = K
         end do
      end do

c     // Compute trial deviatoric strains in 9-component form
c     // e_{n+1} = I_{dev} \epsilon = (I - 1/3 mm^T) \epsilon
      evn1 = trace*one3
      en1(1) = eps(1) - evn1
      en1(2) = eps(2) - evn1
      en1(3) = eps(3) - evn1
      en1(4) = 0.5d0*eps(4)
      en1(5) = en1(4)
      en1(6) = 0.5d0*eps(5)
      en1(7) = en1(6)
      en1(8) = 0.5d0*eps(6)
      en1(9) = en1(8)

c     // Get committed deviatoric plastic strains from history
c     // and map from 6- to 9-component form
      do i = 1,3
         epn(i) = hn(i)
      end do
      epn(4) = 0.5d0*hn(4)
      epn(5) = epn(4)
      epn(6) = 0.5d0*hn(5)
      epn(7) = epn(6)
      epn(8) = 0.5d0*hn(6)
      epn(9) = epn(8)

c     // Get committed effective plastic strain and compute yield stress
      kpn = hn(7)
      syn = sigY + H*kpn

c     // Compute trial deviatoric stress in 9-component form
c     // s_{n+1} = 2 G (e_{n+1} - e^p_n)
c     // Do the bulk terms ...
      do i = 1,9
         sn1(i) = twoG*(en1(i)-epn(i))
      end do

c     // Compute norm of trial deviatoric stress
      nsn1 = 0.d0
      do i = 1,9
         nsn1 = nsn1 + sn1(i)*sn1(i)
      end do
      nsn1 = sqrt(nsn1)

c     // Check yield condition
c     // F = \|s\| - sqrt(2/3)*\sigma_y
      Fn1 = nsn1 - root23*syn

c     // First iteration or an elastic step
      if (niter.eq.0 .or. Fn1.lt.0.d0) then

c       // No change in history variables n+1 <-- n
         do i = 1,7
	    h1(i) = hn(i)
         end do

c       // Add deviatoric part to material stress after mapping
c       // to 6-component form
         do i = 1,3
	    sig(i) = sig(i) + sn1(i)
         end do
         sig(4) = sig(4) + 0.5d0*(sn1(4)+sn1(5))
         sig(5) = sig(5) + 0.5d0*(sn1(6)+sn1(7))
         sig(6) = sig(6) + 0.5d0*(sn1(8)+sn1(9))

c       // Add deviatoric part to material tangent after mapping to
c       // 6-component form
c       // D_t += 2 G I_o I_{dev}
c       // Do the bulk terms ...
         tmp1 = twoG*two3
         do i = 1,3
            dd(i,i) = dd(i,i) + tmp1
         end do
c       // ... and do the deviatoric terms ...
         do i = 4,6
            dd(i,i) = dd(i,i) + G
         end do
c       // ... and finally do the off-diagonal terms
         tmp1 = G*two3          ! G*two3 == 2.d0*G*one3
         dd(2,1) = dd(2,1) - tmp1
         dd(3,1) = dd(3,1) - tmp1
         dd(3,2) = dd(3,2) - tmp1

c     // Plastic step
      else

c       // Compute normal to yield surface
c       // n_{n+1} = \frac{s_{n+1}}{\|s_{n+1}\|}
         tmp1 = 1.d0/nsn1       ! Save some flops
         do i = 1,9
	    nn1(i) = sn1(i)*tmp1
         end do

c       // Compute consistency parameter
         dlam = Fn1 / (twoG+two3*H)

c       // Update plastic deviatoric strains in 9-component form
         do i = 1,9
	    epn(i) = epn(i) + dlam*nn1(i)
         end do

c       // Update effective plastic strain
c       // \kappa{n+1} = \kappa_n + sqrt(2/3) (\Delta \lambda)
         kpn1 = kpn + root23*dlam

c       // Put trial history variables into vector after mapping back
c       // to 6-component form
         do i = 1,3
	    h1(i) = epn(i)
         end do
         h1(4) = 2.d0*epn(4)
         h1(5) = 2.d0*epn(6)
         h1(6) = 2.d0*epn(8)
         h1(7) = kpn1

c       // Return deviatoric stress, in 9-component form, to
c       // yield surface
         tmp1 = twoG*dlam
         do i = 1,9
            sn1(i) = sn1(i) - tmp1*nn1(i)
         end do

c       // Add deviatoric part to material stress after mapping to
c       // 6-component form
         do i = 1,3
	    sig(i) = sig(i) + sn1(i)
         end do
         sig(4) = sig(4) + 0.5d0*(sn1(4)+sn1(5))
         sig(5) = sig(5) + 0.5d0*(sn1(6)+sn1(7))
         sig(6) = sig(6) + 0.5d0*(sn1(8)+sn1(9))

c       // Add deviatoric part to material tangent
c       // D_t += tmp_1 I_o I_{dev}
         tmp1 = twoG*(1.d0-twoG*dlam/nsn1)
c       // Do the bulk terms ...
         tmp2 = two3*tmp1
         do i = 1,3
            dd(i,i) = dd(i,i) + tmp2
         end do
c       // ... and do the deviatoric terms ...
         tmp2 = 0.5d0*tmp1
         do i = 4,6
            dd(i,i) = dd(i,i) + tmp2
         end do
c       // ... and finally do the off-diagonal terms
         tmp2 = one3*tmp1
         dd(2,1) = dd(2,1) - tmp2
         dd(3,1) = dd(3,1) - tmp2
         dd(3,2) = dd(3,2) - tmp2

c       // Map normal vector to 6-component form
c       // n_6 = P^T n_9
         nn1(4) = 0.5d0*(nn1(4)+nn1(5))
         nn1(5) = 0.5d0*(nn1(6)+nn1(7))
         nn1(6) = 0.5d0*(nn1(8)+nn1(9))

c       // Compute temporary 6-component vector
c       // a = n^T I_{dev} = n - 1/3 (n^T m) m
         tmp1 = one3*(nn1(1)+nn1(2)+nn1(3))
         do i = 1,3
	    a(i) = nn1(i)-tmp1
         end do
         do i = 4,6
	    a(i) = nn1(i)
         end do

c       // Add deviatoric part to material tangent
c       // D_t += tmp_2 n n^T I_{dev} = tmp2 n a^T
         tmp2 = twoG*G*(1.d0/(G+one3*H)-2.d0*dlam/nsn1)
         do i = 1,6
	    do j = 1,i
               dd(i,j) = dd(i,j) - tmp2*nn1(i)*a(j)
	    end do
         end do

      end if

c     // Complete symmetric tangent matrix
      do i = 1,6
         do j = 1,i-1
	    dd(j,i) = dd(i,j)
         end do
      end do

      end

max_stars_repo_name: kuanshi/ductile-fracture


text: The Diversity Coalition is Blacks for Effective Community Action (BECA), http://www.yolanocdr.org/ Californians for Disability Rights, Council on American/Islamic Relations, Davis Asians for Racial Equality, FORGE (For Gay Equality), Mexican American Concilio of Yolo County, Yolo County Rainbow Coalition.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
# include "define.inc"

! RN 2008/05/26:
!  give fpp -DVAMPIR to use Vampire analyzer
!  Vampire analyzer is now replaced by Intel Trace Collector and Analyzer,
!  but this may work with them.

module prof
  implicit none

  public :: init_prof
  public :: prof_entering
  public :: prof_leaving

  private

# ifdef VAMPIR
  integer, parameter :: maxregions = 20

  character (len=80), dimension (maxregions,2) :: regions
  integer :: nregions
# endif

contains

  subroutine init_prof
# ifdef VAMPIR
    implicit none
    nregions = 0
# endif
  end subroutine init_prof

# ifdef VAMPIR
  subroutine find_region (name, group, iregion)
    implicit none
    character(*), intent (in) :: name, group
    integer, intent (out) :: iregion
    integer :: i, ierr

    do i = 1, nregions
       if (name == trim(regions(i,1)) .and. group == trim(regions(i,2))) then
          iregion = i
          return
       end if
    end do
    if (nregions >= maxregions) then
       iregion = 0
       return
    end if
    nregions = nregions + 1
    regions(nregions,1) = name
    regions(nregions,2) = group
    iregion = nregions
    call vtsymdef (iregion, name, group, ierr)
    return
  end subroutine find_region
#endif

  subroutine prof_entering (name, group)
    implicit none
    character(*), intent (in) :: name
    character(*), intent (in), optional :: group
# ifdef VAMPIR
    integer :: i, ierr

    if (present(group)) then
       call find_region (name, group, i)
    else
       call find_region (name, "user", i)
    end if
    if (i == 0) return
    call vtbegin (i, ierr)
#endif
  end subroutine prof_entering

  subroutine prof_leaving (name, group)
    implicit none
    character(*), intent (in) :: name
    character(*), intent (in), optional :: group
# ifdef VAMPIR
    integer :: i, ierr

    if (present(group)) then
       call find_region (name, group, i)
    else
       call find_region (name, "user", i)
    end if
    if (i == 0) return
    call vtend (i, ierr)
#endif
  end subroutine prof_leaving

end module prof

ext: fpp


text:
program bugSize
  integer, pointer :: x(:)
  allocate(x(4))
  call f(x)

contains

subroutine f(x)
integer, intent(in) :: x(:)
real :: E(size(x))
E = x
E = [ 1, 2, 3, 4 ]
print*, E
end subroutine

end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


!-------------------------------------------------------------------------------
!>
!! Calendar module
!!
!! @par Description
!!         This module provides the subroutines for calendar.
!!
!! @author NICAM developers
!<
!-------------------------------------------------------------------------------
module mod_calendar
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use mod_precision
!ESC!  use mod_stdio
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: calendar_setup
  public :: calendar_ointvl
  public :: calendar_ss2cc
  public :: calendar_ssaft
  public :: calendar_secdy
  public :: calendar_ss2ds
  public :: calendar_ss2yd
  public :: calendar_ym2dd
  public :: calendar_ds2ss
  public :: calendar_dayyr
  public :: calendar_ss2yh
  public :: calendar_ss2ym
  public :: calendar_xx2ss
  public :: calendar_cc2yh
  public :: calendar_yh2ss
  public :: calendar_PERPR
  public :: calendar_dd2ym
  public :: calendar_daymo

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, private :: isecdy, isechr, jy, jy4, jcent, jcent4
  integer, private :: idays0, idy, ileap, id, m, idayyr, jyear, jmonth
  !
  !--- flag of automatic or not
  logical, private, save ::  oauto = .true.
  !<---                      yr=0-999     : 360day
  !<---                      yr=1000-1899 : 365day
  !<---                      yr=1900-     : gregorian
  !
  !--- flag of the gregorian calendar or not
  logical, private, save :: ogrego = .true.
  !
  !--- number of days in a month
  integer, private, save :: monday ( 12,2 )
  data   monday /                            &
       31,28,31,30,31,30,31,31,30,31,30,31,  &
       31,29,31,30,31,30,31,31,30,31,30,31 /
  !
  !--- flag of ideal calender (n day per month)
  logical, private, save :: oideal = .false.
  !------ 1 month = x days in the ideal case
  integer, private, save :: idaymo = 30
  !------ 1 year = x months in the ideal case
  integer, private, save :: imonyr = 12
  !
  !--- flag of perpetual or not
  logical, private, save :: operpt = .false.
  !------ perpetual date(year)
  integer, private, save :: iyrpp  = 0
  !------ perpetual date(month)
  integer, private, save :: imonpp = 3
  !------ perpetual date(day)
  integer, private, save :: idaypp = 21
  !
  !--- 1 minute = x sec.
  integer, private, save :: isecmn = 60
  !--- 1 hour = x minutes
  integer, private, save :: iminhr = 60
  !--- 1 day = x hours
  integer, private, save :: ihrday = 24
  !
  logical, private, save :: ooperz = .false.
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> setup
  subroutine calendar_setup
    implicit none

    namelist  /nm_calendar/  &
         oauto ,             &
         ogrego,             &
         oideal,             &
         operpt,             &
         idaymo,             &
         imonyr,             &
         iyrpp,              &
         imonpp,             &
         idaypp,             &
         isecmn,             &
         iminhr,             &
         ihrday

    integer :: ierr

!ESC!    rewind(IO_FID_CONF)
!ESC!    read(IO_FID_CONF, nm_calendar, iostat=ierr)
!ESC!    if ( ierr < 0 ) then
!ESC!       write(IO_FID_LOG,*) 'Msg : Sub[calendar_setup]/Mod[calendar]'
!ESC!       write(IO_FID_LOG,*) ' *** Not found namelist.'
!ESC!       write(IO_FID_LOG,*) ' *** Use default values.'
!ESC!    elseif( ierr > 0 ) then
!ESC!       write(*,*) 'Msg : Sub[calendar_setup]/Mod[calendar]'
!ESC!       write(*,*) ' *** WARNING : Not appropriate names in namelist!! CHECK!!'
!ESC!    endif

    return
  end subroutine calendar_setup

  !-----------------------------------------------------------------------------
  subroutine calendar_PERPR( &    !! calendar, refer to fixed date
       IYEAR ,               & !--- OUT
       IMONTH,               & !--- OUT
       IDAY,                 & !--- OUT
       OOPERP                & !--- OUT
       )
    implicit none

    integer, intent(out) :: IYEAR
    integer, intent(out) :: IMONTH
    integer, intent(out) :: IDAY
    logical, intent(out) :: OOPERP

    OOPERP = OPERPT
    IYEAR  = IYRPP
    IMONTH = IMONPP
    IDAY   = IDAYPP

    RETURN
  end subroutine calendar_PERPR

  !-----------------------------------------------------------------------------
  subroutine calendar_daymo(&
       ndaymo,              & !--- OUT : day
       iyear,               & !--- IN : year
       imonth               & !--- IN : month
       )
    !--- calendar, no.of day in a month
    implicit none

    integer, intent(out) :: ndaymo
    integer, intent(in) :: iyear
    integer, intent(in) :: imonth

    if ( oauto  ) then
       if ( iyear >= 1900 ) then
          ogrego = .true.
       else
          ogrego = .false.
          if ( iyear >= 1000 ) then
             oideal = .false.
          else
             oideal = .true.
             idaymo = 30
             imonyr = 12
          endif
       endif
    endif

    if ( ogrego ) then
       if ( ocleap( iyear ) ) then
          ndaymo = monday( imonth,2 )
       else
          ndaymo = monday( imonth,1 )
       endif
    elseif( .not. oideal ) then
       ndaymo = monday( imonth,1 )
    else
       ndaymo = idaymo
    endif

    return
  end subroutine calendar_daymo

  !-----------------------------------------------------------------------------
  subroutine calendar_dayyr(&
       ndayyr,              & !-- OUT : day
       iyear                & !---IN : year
       )
    !--- calendar, no.of day in an year
    implicit none

    integer, intent(out) :: ndayyr
    integer, intent(in) :: iyear

    if ( oauto  ) then
       if ( iyear >= 1900 ) then
          ogrego = .true.
       else
          ogrego = .false.
          if ( iyear >= 1000 ) then
             oideal = .false.
          else
             oideal = .true.
             idaymo = 30
             imonyr = 12
          endif
       endif
    endif

    if ( ogrego ) then
       if ( ocleap( iyear ) ) then
          ndayyr = 366
       else
          ndayyr = 365
       endif
    elseif( .not. oideal ) then
       ndayyr = 365
    else
       ndayyr = idaymo*imonyr
    endif

    return
  end subroutine calendar_dayyr

  !-----------------------------------------------------------------------------
  subroutine calendar_monyr( &
       nmonyr,               & !--- OUT : month
       iyear                 & !--- IN : year
       )
    !--- calendar, no.of month in an year
    implicit none

    integer, intent(out) :: nmonyr
    integer, intent(in) :: iyear

    if ( oauto  ) then
       nmonyr = 12
    else
       nmonyr = imonyr
    endif

    return
  end subroutine calendar_monyr

  !-----------------------------------------------------------------------------
  subroutine calendar_secdy( &
       nsecdy                & !--- OUT : sec
       )
    !--- calendar, no.of sec. in a day
    implicit none

    integer, intent(out) :: nsecdy
    nsecdy = isecmn*iminhr*ihrday

    return
  end subroutine calendar_secdy

  !-----------------------------------------------------------------------------
  subroutine calendar_secmi( &
       nsecmi                & !--- OUT
       )
    !--- calendar, no of sec. in a minute
    implicit none

    integer, intent(out) :: nsecmi
    nsecmi = isecmn

    return
  end subroutine calendar_secmi

  !-----------------------------------------------------------------------------
  subroutine calendar_sechr( &
       nsechr                &
       )
    !--- calendar, no.of sec. in an hour
    implicit none

    integer, intent(out) :: nsechr
    nsechr = isecmn*iminhr

    return
  end subroutine calendar_sechr

  !-----------------------------------------------------------------------------
  subroutine calendar_ss2ds( &
       idays, &
       rsec,  &
       dsec   )
    implicit none

    integer,  intent(out) :: idays
    real(DP), intent(out) :: rsec
    real(DP), intent(in)  :: dsec
    !---------------------------------------------------------------------------

    isecdy = isecmn * iminhr * ihrday
    idays  = int( dsec / real(isecdy,kind=DP) ) + 1
    rsec   = dsec - real((idays-1),kind=DP) * real(isecdy,kind=DP)

    if ( nint(rsec) >= isecdy ) then
       idays = idays + 1
       rsec  = rsec - real(isecdy,kind=DP)
    endif

    return
  end subroutine calendar_ss2ds

  !-----------------------------------------------------------------------------
  subroutine calendar_ds2ss( &
       dsec,  &
       idays, &
       rsec   )
    implicit none

    real(DP), intent(out) :: dsec
    integer,  intent(in)  :: idays
    real(DP), intent(in)  :: rsec
    !---------------------------------------------------------------------------

    isecdy = isecmn * iminhr * ihrday
    dsec   = real(idays-1,kind=DP) * real(isecdy,kind=DP) + real(rsec,kind=DP)

    return
  end subroutine calendar_ds2ss

  !-----------------------------------------------------------------------------
  subroutine calendar_rs2hm( &
       ihour, &
       imin,  &
       isec,  &
       rsec   )
    implicit none

    integer,  intent(out) :: ihour
    integer,  intent(out) :: imin
    integer,  intent(out) :: isec
    real(DP), intent(in)  :: rsec
    !---------------------------------------------------------------------------

    isechr = isecmn * iminhr
    ihour  = int ( rsec / real(isechr,kind=DP) )
    imin   = int ( ( rsec - real(ihour*isechr,kind=DP) ) / real(isecmn,kind=DP) )
    isec   = nint( rsec - real(ihour*isechr,kind=DP) - real(imin*isecmn,kind=DP) )

    if ( isec >= isecmn ) then
       imin  = imin + 1
       isec  = isec - isecmn
    endif

    if ( imin == iminhr ) then
       ihour = ihour + 1
       imin  = imin  - iminhr
    endif

    return
  end subroutine calendar_rs2hm

  !-----------------------------------------------------------------------------
  subroutine calendar_hm2rs( &
       rsec  ,               & !--- OUT
       ihour ,               & !--- IN
       imin  ,               & !--- IN
       isec                  & !--- IN
       )
    !--- calendar, hhmmss -> sec.
    implicit none

    real(DP), intent(out) :: rsec
    integer, intent(in) :: ihour
    integer, intent(in) :: imin
    integer, intent(in) :: isec

    rsec = real( ihour*isecmn*iminhr + imin*isecmn + isec, kind=DP )

    return
  end subroutine calendar_hm2rs

  !-----------------------------------------------------------------------------
  subroutine calendar_dd2ym( &
       iyear,  &
       imonth, &
       iday,   &
       idays   )
    implicit none

    integer, intent(out) :: iyear
    integer, intent(out) :: imonth
    integer, intent(out) :: iday
    integer, intent(in)  :: idays
    !---------------------------------------------------------------------------

    if ( oauto ) then
       if ( idays >= 693961 ) then       !" 1900*365+1900/4-19+5
          ogrego = .true.
       else
          ogrego = .false.
          if ( idays >= 1000*365 ) then
             oideal = .false.
          else
             oideal = .true.
             idaymo = 30
             imonyr = 12
          endif
       endif
    endif

    if ( operpt .AND. ooperz ) then
       iyear  = iyrpp
       imonth = imonpp
       iday   = idaypp
       return
    endif

    if ( ogrego ) then
       jy     = int(real(idays,kind=DP)/365.24_DP)
1100   continue
       jy4    = (jy+3)/4
       jcent  = (jy+99)/100
       jcent4 = (jy+399)/400
       idays0 = jy*365+jy4-jcent+jcent4
       if ( idays <= idays0 ) then
          jy = jy -1
          if ( jy >= 0 ) goto 1100
       endif
       iyear = jy
       idy   = idays - idays0
       if ( ocleap( iyear ) ) then
          ileap  = 2
       else
          ileap  = 1
       endif
    elseif( .not. oideal ) then
       iyear = idays/365
       idy   = idays - iyear*365
       ileap = 1
    endif
    if ( ogrego .OR. .not. oideal ) then
       id = 0
       do m = 1, 12
          id = id + monday(m,ileap)
          if ( idy <= id ) then
             imonth = m
             iday   = idy-id+monday(m,ileap)
             exit
          endif
       enddo
    else
       idayyr = idaymo*imonyr
       iyear  = ( idays-1 ) / idayyr
       imonth = ( idays-1 - iyear*idayyr )/idaymo+1
       iday   = idays - iyear*idayyr - (imonth-1)*idaymo
    endif

    return
  end subroutine calendar_dd2ym

  !-----------------------------------------------------------------------------
  subroutine calendar_ym2dd( &
       idays ,               & !--- OUT
       iyear ,               & !--- IN
       imonth,               & !--- IN
       iday                  & !--- IN
       )
    !--- calendar, yymmdd -> day
    implicit none
    integer, intent(out) :: idays
    integer, intent(in) :: iyear
    integer, intent(in) :: imonth
    integer, intent(in) :: iday
    if ( oauto  ) then
       if ( iyear >= 1900 ) then
          ogrego = .true.
       else
          ogrego = .false.
          if ( iyear >= 1000 ) then
             oideal = .false.
          else
             oideal = .true.
             idaymo = 30
             imonyr = 12
          endif
       endif
    endif
    if ( ogrego .OR. .not. oideal ) then
       if ( imonth > 0 ) then
          jyear  = iyear + (imonth-1)/12
          jmonth = mod(imonth-1,12)+1
       else
          jyear  = iyear - (-imonth)/12 - 1
          jmonth = 12-mod(-imonth,12)
       endif
    endif

    if ( ogrego ) then
       jy4    = (jyear+3)/4
       jcent  = (jyear+99)/100
       jcent4 = (jyear+399)/400
       idays0 = jyear*365+jy4-jcent+jcent4
       if ( ocleap( jyear ) ) then
          ileap = 2
       else
          ileap = 1
       endif
    elseif( .not. oideal ) then
       idays0 = jyear*365
       ileap  = 1
    endif

    if ( ogrego .OR. .not. oideal ) then
       id = 0
       do m = 1, jmonth-1
          id = id + monday(m,ileap)
       enddo
    else
       idays0 = iyear*idaymo*imonyr
       id     = (imonth-1)*idaymo
    endif
    idays = idays0 + id + iday

    return
  end subroutine calendar_ym2dd
  !-----------------------------------------------------------------------------
  subroutine calendar_ym2yd( &
       idaysy,               & !--- OUT
       iyear ,               & !--- IN
       imonth,               & !--- IN
       iday                  & !--- IN
       )
    !--- calendar, yymmdd -> yydd
    implicit none
    integer, intent(out) :: idaysy
    integer, intent(in) :: iyear
    integer, intent(in) :: imonth
    integer, intent(in) :: iday
    if ( oauto  ) then
       if ( iyear >= 1900 ) then
          ogrego = .true.
       else
          ogrego = .false.
          if ( iyear >= 1000 ) then
             oideal = .false.
          else
             oideal = .true.
             idaymo = 30
             imonyr = 12
          endif
       endif
    endif
    if ( ogrego .OR. .not. oideal ) then
       if ( imonth > 0 ) then
          jyear  = iyear + (imonth-1) / 12
          jmonth = mod(imonth-1,12) + 1
       else
          jyear  = iyear - (-imonth) / 12 - 1
          jmonth = 12 - mod(-imonth,12)
       endif
    endif

    if ( ogrego ) then
       if ( ocleap( jyear ) ) then
          ileap = 2
       else
          ileap = 1
       endif
    elseif( .not. oideal ) then
       ileap  = 1
    endif

    if ( ogrego .OR. .not. oideal ) then
       id = 0
       do m = 1, jmonth-1
          id = id + monday(m,ileap)
       enddo
    else
       id     = (imonth-1)*idaymo
    endif
    idaysy = id + iday

    return
  end subroutine calendar_ym2yd

  !-----------------------------------------------------------------------------
  subroutine calendar_ss2yh( &
       idate, &
       dsec   )
    implicit none

    integer,  intent(out) :: idate(6) ! yymmddhhmmss
    real(DP), intent(in)  :: dsec     ! time(second)

    integer  :: idays ! serial number of day
    real(DP) :: rsec  ! seconds in a day
    !---------------------------------------------------------------------------

    call calendar_ss2ds( idays, rsec, dsec )
    call calendar_dd2ym( idate(1), idate(2), idate(3), idays )
    call calendar_rs2hm( idate(4), idate(5), idate(6), rsec  )

    return
  end subroutine calendar_ss2yh

  !-----------------------------------------------------------------------------
  subroutine calendar_yh2ss( &
       dsec  ,               & !--- OUT
       idate                 & !--- IN
       )
    !--- calendar, date -> sec.
    implicit none

    real(DP), intent(out) :: dsec    !" time
    integer, intent(in) :: idate(6) !" yymmddhhmmss

    integer    idays             !" serial no.of day
    real(DP)    rsec              !" no. of sec. in a day

    call calendar_ym2dd &
         ( idays   , &
         idate(1), idate(2), idate(3) )

    call calendar_hm2rs &
         ( rsec    , &
         idate(4), idate(5), idate(6) )

    call calendar_ds2ss &
         ( dsec  , &
         idays , rsec   )

    return
  end subroutine calendar_yh2ss

  !-----------------------------------------------------------------------------
  subroutine calendar_ss2yd( &
       iyear ,               & !--- OUT
       idaysy,               & !--- OUT
       dsec                  & !--- IN
       )
    !--- calendar, sec. -> yydd
    implicit none

    integer, intent(out) :: iyear
    integer, intent(out) :: idaysy
    real(DP), intent(in) :: dsec

    integer :: imonth, iday
    call calendar_ss2ym &
         ( iyear , imonth, iday  , &
         dsec                    )
    call calendar_ym2yd &
         ( idaysy, &
         iyear , imonth, iday    )

    return
  end subroutine calendar_ss2yd

  !-----------------------------------------------------------------------------
  subroutine calendar_ss2ym( &
       iyear ,               & !--- OUT
       imonth,               & !--- OUT
       iday  ,               & !--- OUT
       dsec                  & !--- IN
       )
    !--- calendar, sec. -> yymmdd
    implicit none

    integer, intent(out) :: iyear
    integer, intent(out) :: imonth
    integer, intent(out) :: iday
    real(DP), intent(in) :: dsec

    integer :: idays
    real(DP) :: rsec

    call calendar_ss2ds &
         ( idays , rsec  , &
         dsec            )
    call calendar_dd2ym &
         ( iyear , imonth, iday  , &
         idays                  )

    return
  end subroutine calendar_ss2ym

  !-----------------------------------------------------------------------------
  subroutine calendar_xx2ss( &
       ddsec, &
       rtdur, &
       hunit, &
       dsec   )
    implicit none

    real(DP),         intent(out) :: ddsec
    real(DP),         intent(in)  :: rtdur
    character(len=*), intent(in)  :: hunit
    real(DP),         intent(in)  :: dsec

    character(len=10) :: hunitx
    integer :: isecmi, isechr, isecdy
    integer :: iyear, imonth, iday, ndaymo, ndayyr
    !---------------------------------------------------------------------------

    hunitx = hunit

    if    ( hunitx(1:1) == 's' .OR. hunitx(1:1) == 'S' ) then

       ddsec = real(rtdur,kind=DP)

    elseif( hunitx(1:2) == 'mi' .OR. hunitx(1:2) == 'MI' ) then

       call calendar_secmi( isecmi )
       ddsec = real(rtdur,kind=DP)*real(isecmi,kind=DP)

    elseif( hunitx(1:1) == 'h' .OR. hunitx(1:1) == 'H' ) then

       call calendar_sechr( isechr )
       ddsec = real(rtdur,kind=DP)*real(isechr,kind=DP)

    elseif( hunitx(1:1) == 'd' .OR. hunitx(1:1) == 'D' ) then

       call calendar_secdy( isecdy )
       ddsec = real(rtdur,kind=DP)*real(isecdy,kind=DP)

    elseif( hunitx(1:2) == 'mo' .OR. hunitx(1:2) == 'MO' ) then

       call calendar_ss2ym( iyear, imonth, iday, dsec )
       call calendar_daymo( ndaymo, iyear, imonth )
       call calendar_secdy( isecdy )
       ddsec = real(rtdur,kind=DP) * real(ndaymo,kind=DP) * real(isecdy,kind=DP)

    elseif( hunitx(1:1) == 'y' .OR. hunitx(1:1) == 'Y' ) then

       call calendar_ss2ym( iyear, imonth, iday, dsec )
       call calendar_dayyr( ndayyr, iyear )
       call calendar_secdy( isecdy )
       ddsec = real(rtdur,kind=DP)*real(ndayyr,kind=DP)*real(isecdy,kind=DP)

    else
       write(*,*) 'xxx cxx2ss: invalid unit : ', hunit, ' [sec] assumed'
       ddsec = rtdur
    endif

    return
  end subroutine calendar_xx2ss

  !-----------------------------------------------------------------------------
  subroutine calendar_cc2yh( &
       itime ,               & !--- OUT
       htime                 & !--- IN
       )
    !--- calendar, character -> date
    implicit none

    integer, intent(out) :: itime ( 6 )
    character(len=*), intent(in) :: htime

    integer :: i
    read ( htime, 2600 ) (itime(i),i=1,6)
2600 format( i4.4,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2 )

    return
  end subroutine calendar_cc2yh

  !-----------------------------------------------------------------------------
  !> second to character
  subroutine calendar_ss2cc( &
       htime,      &
       dsec,       &
       number_only )
    implicit none

    character(len=*), intent(out) :: htime
    real(DP),         intent(in)  :: dsec
    logical,          intent(in), optional :: number_only

    integer :: itime(6)
    !---------------------------------------------------------------------------

    call calendar_ss2yh( itime, dsec )

    if ( present(number_only) ) then
       if ( number_only ) then
          write(htime,'(I4.4,5(I2.2))') itime(1), itime(2), itime(3), &
                                        itime(4), itime(5), itime(6)
          return
       endif
    endif

    write(htime,'(I4.4,5(A,I2.2))') itime(1), '/', itime(2), '/', itime(3), &
                               '-', itime(4), ':', itime(5), ':', itime(6)

    return
  end subroutine calendar_ss2cc

  !-----------------------------------------------------------------------------
  function ocleap( &
       iyear       & !--- IN
       )
    !--- calendar :bissextile or not
    implicit none

    logical :: ocleap
    integer, intent(in) :: iyear

    integer :: iy, iycen, icent

    iy     = mod(iyear,4)
    iycen  = mod(iyear,100)
    icent  = mod(iyear/100,4)
    if ( iy == 0 .AND. ( iycen /= 0 .OR. icent == 0 ) ) then
       ocleap = .true.
    else
       ocleap = .false.
    endif

    return
  end function ocleap

  !-----------------------------------------------------------------------------
  subroutine calendar_ssaft( &
       dseca ,               & !--- OUT
       dsec  ,               & !--- IN
       raftr ,               & !--- IN
       hunit                 & !--- IN
       )
    !--- calendar, time advancing
    implicit none

    real(DP), intent(out) :: dseca
    real(DP), intent(in) :: dsec
    real(DP), intent(in) :: raftr
    character(len=*), intent(in) :: hunit

    integer :: idays, iyear, imonth, iday
    real(DP) :: rsec
    real(DP) :: ddtime

    if ( hunit(1:1) == 'y' .OR. hunit(1:1) == 'Y' &
         .OR. hunit(1:1) == 'mo' .OR. hunit(1:2) == 'MO' ) then
       call calendar_ss2ds &
            ( idays , rsec  , &
            dsec            )
       call calendar_dd2ym &
            ( iyear , imonth, iday  , &
            idays                  )
       if ( hunit(1:1) == 'y' .OR. hunit(1:1) == 'Y' ) then
          iyear  = iyear  + int(raftr)
       elseif( hunit(1:2) == 'mo' .OR. hunit(1:2) == 'MO' ) then
          imonth = imonth + int(raftr)
       endif
       call calendar_ym2dd &
            ( idays , &
            iyear , imonth, iday   )
       call calendar_ds2ss &
            ( dseca , &
            idays , rsec    )
    else
       call calendar_xx2ss( ddtime, raftr, hunit, dsec )
       dseca = dsec + ddtime
    endif

    return
  end subroutine calendar_ssaft

  !-----------------------------------------------------------------------------
  function   calendar_ointvl( &
       dtime ,                & !--- IN
       dtprev ,               & !--- IN
       dtorgn,                & !--- IN
       rintv ,                & !--- IN
       htunit                 & !--- IN
       )
    !--- time step passed ?
    implicit none

    logical :: calendar_ointvl
    real(DP), intent(in) :: dtime
    real(DP), intent(in) :: dtprev
    real(DP), intent(in) :: dtorgn
    real(DP), intent(in) :: rintv
    character(len=*), intent(in) :: htunit

    real(DP) :: ddtime
    character(len=5) :: hunit
    integer :: iyear, imon, iday, iyearp, imonp, idayp
    integer :: iy, imo
    integer :: nmonyr, ndayyr, ndaymo
    real(DP) :: ry, rmo

    hunit = htunit

    if ( dtime == dtprev ) then
       calendar_ointvl = .true.
       return
    endif

    calendar_ointvl = .false.

    call calendar_ss2ym ( iyear, imon, iday, dtime  )
    call calendar_ss2ym ( iyearp, imonp, idayp, dtprev )
    call calendar_xx2ss ( ddtime, rintv, hunit, dtime  )

    if ( dtime >= dtorgn ) then
       if      ( hunit(1:1) == 'y' .OR. hunit(1:1) == 'Y' ) then
          call calendar_monyr( nmonyr, iyear )
          call calendar_dayyr( ndayyr, iyear )
          ry = real( iyear-iyearp, kind=DP )&
             + real( imon-imonp, kind=DP ) / real( nmonyr, kind=DP ) &
             + real( iday-idayp, kind=DP ) / real( ndayyr, kind=DP )
          if ( ry >= rintv ) then
             calendar_ointvl = .true.
          endif
       elseif( hunit(1:2) == 'mo' .OR. hunit(1:2) == 'MO' ) then
          imo = 0
          do iy = iyearp, iyear-1
             call calendar_monyr( nmonyr, iy )
             imo = imo + nmonyr
          enddo
          call calendar_daymo( ndaymo, iyear, imon )
          rmo = real( imon-imonp+imo, kind=DP ) &
              + real( iday-idayp, kind=DP ) / real( ndaymo, kind=DP )
          if ( rmo >= rintv ) then
             calendar_ointvl = .true.
          endif
       elseif(      calendar_dgaus((dtime -dtorgn)/ddtime)  &
            > calendar_dgaus((dtprev-dtorgn)/ddtime) ) then
          calendar_ointvl = .true.
       endif
    endif

    return
  end function calendar_ointvl

  !-----------------------------------------------------------------------------
  function calendar_dgaus( &
        dx                 & !--- IN
        )
    !--- dicard gaussian
    implicit none
    real(DP) :: calendar_dgaus
    real(DP), intent(in) :: dx

    calendar_dgaus = aint(dx) + aint(dx - aint(dx) + 1.0_DP) - 1.0_DP

  end function calendar_dgaus

end module mod_calendar