!****h* ROBODoc/H5O (F03)
!
! NAME
!  H5O_PROVISIONAL
!
! PURPOSE
!  This file contains Fortran 90 and Fortran 2003 interfaces for H5O functions.
!  It contains the same functions as H5Off_F90.f90 but includes the
!  Fortran 2003 functions and the interface listings. This file will be compiled
!  instead of H5Off_F90.f90 if Fortran 2003 functions are enabled.
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the files COPYING and Copyright.html.  COPYING can be found at the root   *
!   of the source code distribution tree; Copyright.html can be found at the  *
!   root level of an installed copy of the electronic HDF5 document set and   *
!   is linked from the top-level documents page.  It can also be found at     *
!   http://hdfgroup.org/HDF5/doc/Copyright.html.  If you do not have          *
!   access to either file, you may request a copy from help@hdfgroup.org.     *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! NOTES
!                         *** IMPORTANT ***
!  If you add a new H5P function you must add the function name to the
!  Windows dll file 'hdf5_fortrandll.def.in' in the fortran/src directory.
!  This is needed for Windows based operating systems.
!
!*****

MODULE H5O_PROVISIONAL

  USE H5GLOBAL
  USE, INTRINSIC :: ISO_C_BINDING

  IMPLICIT NONE

  enum, bind(c)
     enumerator :: H5O_TYPE_UNKNOWN_F = -1
     enumerator :: H5O_TYPE_GROUP_F, H5O_TYPE_DATASET_F, H5O_TYPE_NAMED_DATATYPE_F, H5O_TYPE_NTYPES_F
  end enum

!****t* H5T (F03)/h5o_info_t
!
! Fortran2003 Derived Type:
!
  TYPE, BIND(C) :: space_t
     INTEGER(hsize_t) :: total ! Total space for storing object header in file
     INTEGER(hsize_t) :: meta  ! Space within header for object header metadata information
     INTEGER(hsize_t) :: mesg  ! Space within header for actual message information
     INTEGER(hsize_t) :: free  ! Free space within object header
  END TYPE space_t

  TYPE, BIND(C) :: mesg_t
     INTEGER(c_int64_t) :: present ! Flags to indicate presence of message type in header
     INTEGER(c_int64_t) :: shared  ! Flags to indicate message type is shared in header
  END TYPE mesg_t

  TYPE, BIND(C) :: hdr_t
     INTEGER :: version ! Version number of header format in file
     INTEGER :: nmesgs  ! Number of object header messages
     INTEGER :: nchunks ! Number of object header chunks
     INTEGER :: flags   ! Object header status flags
     TYPE(space_t)  :: space
     TYPE(mesg_t)   :: mesg
  END TYPE hdr_t

  ! Extra metadata storage for obj & attributes
  TYPE, BIND(C) :: H5_ih_info_t
     INTEGER(hsize_t) :: index_size ! btree and/or list
     INTEGER(hsize_t) :: heap_size
  END TYPE H5_ih_info_t

  TYPE, BIND(C) :: meta_size_t
     TYPE(H5_ih_info_t) :: obj  ! v1/v2 B-tree & local/fractal heap for groups, B-tree for chunked datasets
     TYPE(H5_ih_info_t) :: attr ! v2 B-tree & heap for attributes
  ENDTYPE meta_size_t

  TYPE, BIND(C) :: h5o_info_t
     INTEGER(C_LONG)  :: fileno     ! File number that object is located in
     INTEGER(haddr_t) :: addr       ! Object address in file
     INTEGER(C_INT)   :: type       ! Basic object type (group, dataset, etc.)
     INTEGER          :: rc         ! Reference count of object

     INTEGER, DIMENSION(8) :: atime ! Access time         !    -- NOTE --
     INTEGER, DIMENSION(8) :: mtime ! Modification time   ! Returns an integer array
     INTEGER, DIMENSION(8) :: ctime ! Change time         ! as specified in the Fortran
     INTEGER, DIMENSION(8) :: btime ! Birth time          ! intrinsic DATE_AND_TIME(VALUES)

     INTEGER(hsize_t) :: num_attrs  ! # of attributes attached to object

     TYPE(hdr_t) :: hdr

     TYPE(meta_size_t) :: meta_size
  END TYPE h5o_info_t

!*****

CONTAINS

!****s* H5O (F03)/h5ovisit_f_F03
!
! NAME
!  h5ovisit_f
!
! PURPOSE
!  Recursively visits all objects starting from a specified object.
!
! Inputs:
!  object_id  - Identifier of the object at which the recursive iteration begins.
!  index_type - Type of index; valid values include:
!                H5_INDEX_NAME_F
!                H5_INDEX_CRT_ORDER_F
!  order      - Order in which index is traversed; valid values include:
!                H5_ITER_DEC_F
!                H5_ITER_INC_F
!                H5_ITER_NATIVE_F
!  op 	      - Callback function passing data regarding the group to the calling application
!  op_data    - User-defined pointer to data required by the application for its processing of the group
!
! Outputs:
!  return_value - returns the return value of the first operator that returns a positive value, or
!                 zero if all members were processed with no operator returning non-zero.
!  hdferr       - Returns 0 if successful and -1 if fails
! AUTHOR
!  M. Scot Breitenfeld
!  November 19, 2008
!
! Fortran2003 Interface:
  SUBROUTINE h5ovisit_f(object_id, index_type, order, op, op_data, return_value, hdferr)
    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER(HID_T), INTENT(IN) :: object_id
    INTEGER, INTENT(IN) :: index_type
    INTEGER, INTENT(IN) :: order

    TYPE(C_FUNPTR):: op
    TYPE(C_PTR)   :: op_data
    INTEGER, INTENT(OUT) :: return_value
    INTEGER, INTENT(OUT) :: hdferr
!*****

    INTERFACE
       INTEGER FUNCTION h5ovisit_c(object_id, index_type, order, op, op_data)
         USE, INTRINSIC :: ISO_C_BINDING
         USE H5GLOBAL
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5OVISIT_C'::h5ovisit_c
         !DEC$ENDIF
         INTEGER(HID_T), INTENT(IN) :: object_id
         INTEGER, INTENT(IN) :: index_type
         INTEGER, INTENT(IN) :: order
         TYPE(C_FUNPTR), VALUE :: op
         TYPE(C_PTR), VALUE :: op_data
       END FUNCTION h5ovisit_c
    END INTERFACE

    return_value = h5ovisit_c(object_id, index_type, order, op, op_data)

    IF(return_value.GE.0)THEN
       hdferr = 0
    ELSE
       hdferr = -1
    END IF

  END SUBROUTINE h5ovisit_f

!****s* H5O (F03)/h5oget_info_by_name_f_F03
!
! NAME
!  h5oget_info_by_name_f
!
! PURPOSE
!  Retrieves the metadata for an object, identifying the object by location and relative name.
!
! Inputs:
!  loc_id      - File or group identifier specifying location of group
!                in which object is located.
!  name        - Name of group, relative to loc_id.
!
! Outputs:
!  object_info - Buffer in which to return object information.
!  hdferr      - Returns 0 if successful and -1 if fails.
!
! Optional parameters:
!  lapl_id     - Link access property list.
!
! AUTHOR
!  M. Scot Breitenfeld
!  December 1, 2008
!
! Fortran2003 Interface:
  SUBROUTINE h5oget_info_by_name_f(loc_id, name, object_info, hdferr, lapl_id)

    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER(HID_T)  , INTENT(IN)            :: loc_id
    CHARACTER(LEN=*), INTENT(IN)            :: name
    TYPE(h5o_info_t), INTENT(OUT), TARGET   :: object_info
    INTEGER         , INTENT(OUT)           :: hdferr
    INTEGER(HID_T)  , INTENT(IN) , OPTIONAL :: lapl_id
!*****
    INTEGER         :: corder_valid
    INTEGER(SIZE_T) :: namelen
    INTEGER(HID_T)  :: lapl_id_default
    TYPE(C_PTR)     :: ptr

    INTERFACE
       INTEGER FUNCTION h5oget_info_by_name_c(loc_id, name, namelen, lapl_id_default, &
           object_info)
         USE H5GLOBAL
         USE, INTRINSIC :: ISO_C_BINDING
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5OGET_INFO_BY_NAME_C'::h5oget_info_by_name_c
         !DEC$ENDIF
         !DEC$ATTRIBUTES reference :: name
         INTEGER(HID_T)  , INTENT(IN)  :: loc_id
         CHARACTER(LEN=*), INTENT(IN)  :: name
         INTEGER(SIZE_T) , INTENT(IN)  :: namelen
         INTEGER(HID_T)  , INTENT(IN)  :: lapl_id_default
         TYPE(C_PTR),VALUE             :: object_info

       END FUNCTION h5oget_info_by_name_c
    END INTERFACE

    namelen = LEN(name)

    lapl_id_default = H5P_DEFAULT_F
    IF(PRESENT(lapl_id)) lapl_id_default = lapl_id

    ptr = C_LOC(object_info)

    hdferr = H5Oget_info_by_name_c(loc_id, name, namelen, lapl_id_default, ptr)

  END SUBROUTINE H5Oget_info_by_name_f

!****s* H5O (F03)/h5oget_info_f_F03
!
! NAME
!  h5oget_info_f
!
! PURPOSE
!  Retrieves the metadata for an object specified by an identifier.
!
! Inputs:
!  object_id   - Identifier for target object.
!
! Outputs:
!  object_info - Buffer in which to return object information.
!  hdferr      - Returns 0 if successful and -1 if fails.
!
! AUTHOR
!  M. Scot Breitenfeld
!  May 11, 2012
!
! Fortran2003 Interface:
  SUBROUTINE h5oget_info_f(object_id, object_info, hdferr)

    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER(HID_T)  , INTENT(IN)            :: object_id
    TYPE(h5o_info_t), INTENT(OUT), TARGET   :: object_info
    INTEGER         , INTENT(OUT)           :: hdferr
!*****
    TYPE(C_PTR) :: ptr

    INTERFACE
       INTEGER FUNCTION h5oget_info_c(object_id, object_info)
         USE H5GLOBAL
         USE, INTRINSIC :: ISO_C_BINDING
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5OGET_INFO_C'::h5oget_info_c
         !DEC$ENDIF
         INTEGER(HID_T), INTENT(IN)  :: object_id
         TYPE(C_PTR), VALUE          :: object_info

       END FUNCTION h5oget_info_c
    END INTERFACE

    ptr = C_LOC(object_info)
    hdferr = H5Oget_info_c(object_id, ptr)

  END SUBROUTINE H5Oget_info_f

!****s* H5O (F03)/h5oget_info_by_idx_f_F03
!
! NAME
!  h5oget_info_by_idx_f
!
! PURPOSE
!  Retrieves the metadata for an object, identifying the object by an index position.
!
! Inputs:
!  loc_id      - File or group identifier specifying location of group
!                in which object is located.
!  group_name  - Name of group in which object is located.
!  index_field - Index or field that determines the order.
!  order       - Order within field or index.
!  n           - Object for which information is to be returned
!
! Outputs:
!  object_info - Buffer in which to return object information.
!  hdferr      - Returns 0 if successful and -1 if fails.
!
! Optional parameters:
!  lapl_id     - Link access property list. (Not currently used.)
!
! AUTHOR
!  M. Scot Breitenfeld
!  May 11, 2012
!
! Fortran2003 Interface:
  SUBROUTINE h5oget_info_by_idx_f(loc_id, group_name, index_field, order, n, &
       object_info, hdferr, lapl_id)

    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER(HID_T)  , INTENT(IN)            :: loc_id
    CHARACTER(LEN=*), INTENT(IN)            :: group_name
    INTEGER         , INTENT(IN)            :: index_field
    INTEGER         , INTENT(IN)            :: order
    INTEGER(HSIZE_T), INTENT(IN)            :: n
    TYPE(h5o_info_t), INTENT(OUT), TARGET   :: object_info
    INTEGER         , INTENT(OUT)           :: hdferr
    INTEGER(HID_T)  , INTENT(IN) , OPTIONAL :: lapl_id
!*****
    INTEGER         :: corder_valid
    INTEGER(SIZE_T) :: namelen
    INTEGER(HID_T)  :: lapl_id_default
    TYPE(C_PTR)     :: ptr

    INTERFACE
       INTEGER FUNCTION h5oget_info_by_idx_c(loc_id, group_name, namelen, &
            index_field, order, n, lapl_id_default, object_info)
         USE H5GLOBAL
         USE, INTRINSIC :: ISO_C_BINDING
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5OGET_INFO_BY_IDX_C'::h5oget_info_by_idx_c
         !DEC$ENDIF
         !DEC$ATTRIBUTES reference :: group_name
         INTEGER(HID_T)  , INTENT(IN)  :: loc_id
         CHARACTER(LEN=*), INTENT(IN)  :: group_name
         INTEGER(SIZE_T) , INTENT(IN)  :: namelen
         INTEGER         , INTENT(IN)  :: index_field
         INTEGER         , INTENT(IN)  :: order
         INTEGER(HSIZE_T), INTENT(IN)  :: n
         INTEGER(HID_T)  , INTENT(IN)  :: lapl_id_default
         TYPE(C_PTR), VALUE            :: object_info

       END FUNCTION h5oget_info_by_idx_c
    END INTERFACE

    namelen = LEN(group_name)

    lapl_id_default = H5P_DEFAULT_F
    IF(PRESENT(lapl_id)) lapl_id_default = lapl_id

    ptr = C_LOC(object_info)
    hdferr = H5Oget_info_by_idx_c(loc_id, group_name, namelen, index_field, order, n, lapl_id_default, ptr)

  END SUBROUTINE H5Oget_info_by_idx_f

!****s* H5O (F03)/h5ovisit_by_name_f_F03
!
! NAME
!  h5ovisit_by_name_f
!
! PURPOSE
!  Recursively visits all objects starting from a specified object.
!
! Inputs:
!  loc_id      - Identifier of a file or group.
!  object_name - Name of the object, generally relative to loc_id, that will serve as root of the iteration
!  index_type  - Type of index; valid values include:
!                 H5_INDEX_NAME_F
!                 H5_INDEX_CRT_ORDER_F
!  order       - Order in which index is traversed; valid values include:
!                 H5_ITER_DEC_F
!                 H5_ITER_INC_F
!                 H5_ITER_NATIVE_F
!  op 	       - Callback function passing data regarding the group to the calling application
!  op_data     - User-defined pointer to data required by the application for its processing of the group
!
! Outputs:
!  return_value - Returns the return value of the first operator that returns a positive value, or
!                 zero if all members were processed with no operator returning non-zero.
!  hdferr       - Returns 0 if successful and -1 if fails
!
! Optional parameters:
!  lapl_id      - Link access property list identifier.
!
! AUTHOR
!  M. Scot Breitenfeld
!  November 19, 2008
!
! Fortran2003 Interface:
  SUBROUTINE h5ovisit_by_name_f(loc_id, object_name, index_type, order, op, op_data, &
       return_value, hdferr, lapl_id)
    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER(HID_T)  , INTENT(IN)             :: loc_id
    CHARACTER(LEN=*), INTENT(IN)             :: object_name
    INTEGER         , INTENT(IN)             :: index_type
    INTEGER         , INTENT(IN)             :: order

    TYPE(C_FUNPTR)                           :: op
    TYPE(C_PTR)                              :: op_data
    INTEGER         , INTENT(OUT)            :: return_value
    INTEGER         , INTENT(OUT)            :: hdferr
    INTEGER(HID_T)  , INTENT(IN) , OPTIONAL  :: lapl_id
!*****

    INTEGER(SIZE_T) :: namelen
    INTEGER(HID_T)  :: lapl_id_default
    TYPE(C_PTR)     :: ptr

    INTERFACE
       INTEGER FUNCTION h5ovisit_by_name_c(loc_id, object_name, namelen, index_type, order, &
            op, op_data, lapl_id)
         USE, INTRINSIC :: ISO_C_BINDING
         USE H5GLOBAL
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5OVISIT_BY_NAME_C'::h5ovisit_by_name_c
         !DEC$ENDIF
         !DEC$ATTRIBUTES reference :: object_name
         INTEGER(HID_T)  , INTENT(IN) :: loc_id
         CHARACTER(LEN=*), INTENT(IN) :: object_name
         INTEGER(SIZE_T)              :: namelen
         INTEGER         , INTENT(IN) :: index_type
         INTEGER         , INTENT(IN) :: order
         TYPE(C_FUNPTR)  , VALUE      :: op
         TYPE(C_PTR)     , VALUE      :: op_data
         INTEGER(HID_T)  , INTENT(IN) :: lapl_id
       END FUNCTION h5ovisit_by_name_c
    END INTERFACE

    namelen = LEN(object_name)

    lapl_id_default = H5P_DEFAULT_F
    IF(PRESENT(lapl_id)) lapl_id_default = lapl_id

    return_value = h5ovisit_by_name_c(loc_id, object_name, namelen, index_type, order, &
         op, op_data, lapl_id_default)

    IF(return_value.GE.0)THEN
       hdferr = 0
    ELSE
       hdferr = -1
    END IF

  END SUBROUTINE h5ovisit_by_name_f

END MODULE H5O_PROVISIONAL


text:
      PROGRAM CalIntegrator
	  REAL a, b, s, pi
	  INTEGER n
	  pi = 3.141592654
	  a = 0
	  b = 2
	  CALL qtrap(a, b, s)
	  WRITE (*,*) s
	  END PROGRAM CalIntegrator

	  FUNCTION xFunc(x)
	  		REAL :: xFunc
			REAL :: x
			xFunc = (x**4)*log(x + sqrt(x**2+1))
			RETURN
	  END FUNCTION xFunc

	  SUBROUTINE trapzd(a, b, s, n)
	  INTEGER n
	  REAL a, b, s
	  INTEGER it, j
	  REAL del, sum, tnm, x
	  if (n.eq.1) then
C
	  			    s = 0.5*(b-a)*(xFunc(a) +xFunc(b))
	  else
	  				it = 2**(n-2)
					tnm=it
					del=(b-a)/tnm
					x=a+0.5*del
					sum = 0.
					do j = 1, it
C
									    sum = sum + xFunc(x)
										x = x + del
					end do
					s = 0.5*(s+(b-a)*sum/tnm)
	  endif
	  return
	  END SUBROUTINE trapzd

	  SUBROUTINE qtrap(a, b, s)
	  INTEGER JMAX
	  REAL a, b, s, EPS
	  PARAMETER (EPS=1.e-6, JMAX=20)

	  INTEGER k
	  REAL olds
	  olds=-1.e30
	  do j = 1, JMAX
	  			call trapzd(a, b, s, j)
				jj: if (j.gt.5) then
					kk: if (abs(s-olds).lt.EPS*abs(olds)) then
					         return
					end if kk
					mm: if (s.eq.0 .and. olds.eq.0) then
							return
					end if mm
				end if jj
				olds = s
	 end do
	 END SUBROUTINE qtrap

ext: f

max_stars_repo_name: tqbdev/Calculate-Integration-in-Fortran

max_issues_repo_name: tqbdev/Calculate-Integration-in-Fortran

max_forks_repo_name: tqbdev/Calculate-Integration-in-Fortran


  ! Copyright 2020 University of Warwick

  ! Licensed under the Apache License, Version 2.0 (the "License");
  ! you may not use this file except in compliance with the License.
  ! You may obtain a copy of the License at

  !    http://www.apache.org/licenses/LICENSE-2.0

  ! Unless required by applicable law or agreed to in writing, software
  ! distributed under the License is distributed on an "AS IS" BASIS,
  ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  ! See the License for the specific language governing permissions and
  ! limitations under the License.

MODULE openboundary

  USE shared_data

  IMPLICIT NONE

  REAL(num) :: pfar, rhofar, efar, uxfar, uyfar, uzfar
  REAL(num) :: bxfar, byfar, bzfar, pbc, vnorm

  REAL(num), DIMENSION(0:1) :: vxbc, vybc, vzbc
  REAL(num), DIMENSION(0:1) :: bxbc, bybc, bzbc
  REAL(num), DIMENSION(0:1) :: rbc, ebc

CONTAINS

  SUBROUTINE open_bcs

    REAL(num) :: bperp

    ! Update ghost cells based on Riemann problem with farfield.
    ! Only expected to work perfectly for problems with straight B field
    ! through boundaries which do not drastically change shape during the
    ! simulation.

    ! x_min boundary
    IF (xbc_min == BC_OPEN .AND. proc_x_min == MPI_PROC_NULL) THEN
      DO iy = 0, ny+1
        ! Variables carried out of domain by Riemann invariants
        vxbc(1) = -vx(0,iy)
        vybc(1) =  vy(0,iy)
        vzbc(1) =  vz(0,iy)
        bxbc(1) = -bx(0,iy)
        bybc(1) =  by(1,iy)
        bzbc(1) =  bz(1,iy)
        rbc (1) = rho(1,iy)
        ebc (1) = energy(1,iy)

        pbc  = (gamma - 1.0_num) * energy( 1,iy) * rho( 1,iy)

        ! Farfield values carried into domain
        pfar = (gamma - 1.0_num) * energy(-1,iy) * rho(-1,iy)

        uxfar  = -vx(-2,iy)
        uyfar  =  vy(-2,iy)
        uzfar  =  vz(-2,iy)
        bxfar  =  -bx(-2,iy)
        byfar  =  by(-1,iy)
        bzfar  =  bz(-1,iy)
        rhofar = rho(-1,iy)
        efar = energy(-1,iy)

        vnorm = -vx(0,iy)

        ! Select correct open bc solver
        bperp = SQRT(0.25_num * (by(-1,iy) + by(-1,iy-1))**2 + bzfar**2)

        IF(bperp <= none_zero) THEN
          CALL open_bcs_alfven
        ELSE
          CALL open_bcs_fast
        END IF

        bx (-1,iy) = -bxbc(0)
        by ( 0,iy) =  bybc(0)
        bz ( 0,iy) =  bzbc(0)
        rho( 0,iy) =   rbc(0)
        energy(0,iy) = ebc(0)
      END DO
    END IF

    ! x_max boundary
    IF (xbc_max == BC_OPEN .AND. proc_x_max == MPI_PROC_NULL) THEN
      DO iy = 0, ny+1
        ! Variables carried out of domain by Riemann invariants
        vxbc(1) =  vx(nx,iy)
        vybc(1) =  vy(nx,iy)
        vzbc(1) =  vz(nx,iy)
        bxbc(1) =  bx(nx,iy)
        bybc(1) =  by(nx,iy)
        bzbc(1) =  bz(nx,iy)
        rbc (1) = rho(nx,iy)
        ebc (1) = energy(nx,iy)

        pbc  = (gamma - 1.0_num) * energy(nx  ,iy) * rho(nx  ,iy)

        ! Farfield values carried into domain
        pfar = (gamma - 1.0_num) * energy(nx+2,iy) * rho(nx+2,iy)

        uxfar  =  vx(nx+2,iy)
        uyfar  =  vy(nx+2,iy)
        uzfar  =  vz(nx+2,iy)
        bxfar  =  bx(nx+2,iy)
        byfar  =  by(nx+2,iy)
        bzfar  =  bz(nx+2,iy)
        rhofar = rho(nx+2,iy)
        efar = energy(nx+2,iy)

        vnorm = vx(nx,iy)

        ! Select correct open bc solver
        bperp = SQRT(0.25_num * (by(nx+2,iy) + by(nx+2,iy-1))**2 + bzfar**2)

        IF(bperp <= none_zero) THEN
          CALL open_bcs_alfven
        ELSE
          CALL open_bcs_fast
        END IF

        bx (nx+1,iy) = bxbc(0)
        by (nx+1,iy) = bybc(0)
        bz (nx+1,iy) = bzbc(0)
        rho(nx+1,iy) =  rbc(0)
        energy(nx+1,iy) = ebc(0)
      END DO
    END IF

    ! y_min boundary
    IF (ybc_min == BC_OPEN .AND. proc_y_min == MPI_PROC_NULL) THEN
      DO ix = 0, nx+1
        ! Variables carried out of domain by Riemann invariants
        vxbc(1) = -vy(ix,0)
        vybc(1) =  vx(ix,0)
        vzbc(1) =  vz(ix,0)
        bxbc(1) = -by(ix,0)
        bybc(1) =  bx(ix,1)
        bzbc(1) =  bz(ix,1)
        rbc (1) = rho(ix,1)
        ebc (1) = energy(ix,1)

        pbc  = (gamma - 1.0_num) * energy(ix, 1) * rho(ix, 1)

        ! Farfield values carried into domain
        pfar = (gamma - 1.0_num) * energy(ix,-1) * rho(ix,-1)

        uxfar  = -vy(ix,-2)
        uyfar  =  vx(ix,-2)
        uzfar  =  vz(ix,-2)
        bxfar  = -by(ix,-2)
        byfar  =  bx(ix,-1)
        bzfar  =  bz(ix,-1)
        rhofar = rho(ix,-1)
        efar = energy(ix,-1)

        vnorm = -vy(ix,0)

        ! Select correct open bc solver
        bperp = SQRT(0.25_num * (bx(ix,-1) + bx(ix-1,-1))**2 + bzfar**2)

        IF(bperp <= none_zero) THEN
          CALL open_bcs_alfven
        ELSE
          CALL open_bcs_fast
        END IF

        bx (ix, 0) =  bybc(0)
        by (ix,-1) = -bxbc(0)
        bz (ix, 0) =  bzbc(0)
        rho(ix, 0) =   rbc(0)
        energy(ix,0) = ebc(0)
      END DO
    END IF

    ! y_max boundary
    IF (ybc_max == BC_OPEN .AND. proc_y_max == MPI_PROC_NULL) THEN
      DO ix = 0, nx+1
        ! Variables carried out of domain by Riemann invariants
        vxbc(1) =  vy(ix,ny)
        vybc(1) =  vx(ix,ny)
        vzbc(1) =  vz(ix,ny)
        bxbc(1) =  by(ix,ny)
        bybc(1) =  bx(ix,ny)
        bzbc(1) =  bz(ix,ny)
        rbc (1) = rho(ix,ny)
        ebc (1) = energy(ix,ny)

        pbc  = (gamma - 1.0_num) * energy(ix,ny  ) * rho(ix,ny  )

        ! Farfield values carried into domain
        pfar = (gamma - 1.0_num) * energy(ix,ny+2) * rho(ix,ny+2)

        uxfar  =  vy(ix,ny+2)
        uyfar  =  vx(ix,ny+2)
        uzfar  =  vz(ix,ny+2)
        bxfar  =  by(ix,ny+2)
        byfar  =  bx(ix,ny+2)
        bzfar  =  bz(ix,ny+2)
        rhofar = rho(ix,ny+2)
        efar = energy(ix,ny+2)

        vnorm = vy(ix,ny)

        ! Select correct open bc solver
        bperp = SQRT(0.25_num * (bx(ix,ny+2) + bx(ix-1,ny+2))**2 + bzfar**2)

        IF(bperp <= none_zero) THEN
          CALL open_bcs_alfven
        ELSE
          CALL open_bcs_fast
        END IF

        bx (ix,ny+1) = bybc(0)
        by (ix,ny+1) = bxbc(0)
        bz (ix,ny+1) = bzbc(0)
        rho(ix,ny+1) =  rbc(0)
        energy(ix,ny+1) = ebc(0)
      END DO
    END IF

  END SUBROUTINE open_bcs

  SUBROUTINE open_bcs_fast

    ! Open bc when bx = 0

    REAL(num) :: c0, ct, cf
    REAL(num) :: pg, rhog, cffar, c0far, ctfar
    REAL(num), DIMENSION(3) :: vtest, pstar, vstar, rhostar, pmagstar
    REAL(num), DIMENSION(3) :: bystar, bzstar
    INTEGER :: i

    c0far = SQRT(gamma * pfar / rhofar)
    ctfar = SQRT((byfar**2 + bzfar**2) / rhofar)
    cffar = SQRT(c0far**2 + ctfar**2)

    c0 = SQRT(gamma * (gamma - 1.0_num) * ebc(1))
    ct = SQRT((bybc(1)**2 + bzbc(1)**2) / rbc(1))
    cf = SQRT(c0**2 + ct**2)

    ! Define the speeds of the characteristics to be checked along
    vtest(1) = vnorm + cf
    vtest(2) = vnorm - cf
    vtest(3) = vnorm

    DO i = 1, 3
      IF (vtest(i) >= 0.0_num) THEN
        pstar(i) =  pbc + 0.5_num * (bybc(1)**2 + bzbc(1)**2)
        pmagstar(i) =  0.5_num * (bybc(1)**2 + bzbc(1)**2)
        vstar(i) = vxbc(1)
        rhostar(i) = rbc(1)
        bystar(i) = bybc(1)
        bzstar(i) = bzbc(1)
      ELSE
        pstar(i) = pfar + 0.5_num * (byfar**2 + bzfar**2)
        pmagstar(i) = 0.5_num * (byfar**2 + bzfar**2)
        vstar(i) = uxfar
        rhostar(i) = rhofar
        bystar(i) = byfar
        bzstar(i) = bzfar
      END IF
    END DO

    bxbc(0) = 0.5_num * (bxbc(1) + bxfar)
    bybc(0) = 0.5_num * (bystar(1) + bystar(2))
    bzbc(0) = 0.5_num * (bzstar(1) + bzstar(2))

    pg = 0.5_num &
        * (pstar(1) + pstar(2) + rhofar * cffar * (vstar(1) - vstar(2))) &
        - 0.5_num * (bybc(0)**2 + bzbc(0)**2)
    pg = MAX(pg, none_zero)

    rhog = rhostar(3) + (pg - (pstar(3) - pmagstar(3))) / c0far**2
    rbc(0) = MAX(rhog, none_zero)

    ebc(0) = pg / (gamma - 1.0_num) / rbc(0)

  END SUBROUTINE open_bcs_fast

  SUBROUTINE open_bcs_alfven

    ! Open bc when bperp = 0

    REAL(num) :: lambdayfar, lambdazfar
    REAL(num) :: c0, cx
    REAL(num) :: pg, rhog, c0far, cxfar
    REAL(num) :: lambdag
    REAL(num), DIMENSION(5) :: vtest, pstar, uxstar, rhostar, pmagstar
    REAL(num), DIMENSION(5) :: uystar, lambdaystar, lambdazstar, uzstar
    INTEGER :: i

    lambdayfar = -bxfar * byfar
    lambdazfar = -bxfar * bzfar
    c0far = SQRT(gamma * pfar / rhofar)
    cxfar = SQRT((bxfar**2)/ rhofar)

    c0 = SQRT(gamma * (gamma - 1.0_num) * ebc(1))
    cx = SQRT((bxbc(1)**2)/ rbc(1))

    ! Define the speeds of the characteristics to be checked along
    vtest(1) = vnorm + c0
    vtest(2) = vnorm - c0
    vtest(3) = vnorm + cx
    vtest(4) = vnorm - cx
    vtest(5) = vnorm

    DO i = 1, 5
      IF (vtest(i) >= 0.0_num) THEN
        pstar(i) = pbc + 0.5_num * (bybc(1)**2 + bzbc(1)**2 - bxbc(1)**2)
        pmagstar(i) = 0.5_num * (bybc(1)**2 + bzbc(1)**2 - bxbc(1)**2)
        lambdaystar(i) = -bybc(1) * bxbc(1)
        lambdazstar(i) = -bzbc(1) * bxbc(1)
        uxstar(i) = vxbc(1)
        uystar(i) = vybc(1)
        uzstar(i) = vzbc(1)
        rhostar(i) = rbc(1)
      ELSE
        pstar(i) = pfar + 0.5_num * (byfar**2 + bzfar**2 - bxfar**2)
        pmagstar(i) = 0.5_num * (byfar**2 + bzfar**2 - bxfar**2)
        lambdaystar(i) = lambdayfar
        lambdazstar(i) = lambdazfar
        uxstar(i) = uxfar
        uystar(i) = uyfar
        uzstar(i) = uzfar
        rhostar(i) = rhofar
      END IF
    END DO

    bxbc(0) = 0.5_num * (bxbc(1) + bxfar)

    lambdag = 0.5_num * (lambdaystar(3) + lambdaystar(4) &
        + rhofar * cxfar * (uystar(3) - uystar(4)))
    IF (ABS(bxbc(0)) <= none_zero) bxbc(0) = none_zero
    bybc(0) = -lambdag / bxbc(0)

    lambdag = 0.5_num * (lambdazstar(3) + lambdazstar(4) &
        + rhofar * cxfar * (uzstar(3) - uzstar(4)))
    IF (ABS(bxbc(0)) <= none_zero) bxbc(0) = none_zero
    bzbc(0) = -lambdag  / bxbc(0)

    pg = 0.5_num &
        * (pstar(1) + pstar(2) + rhofar * c0far * (uxstar(1) - uxstar(2))) &
        - 0.5_num * (bybc(0)**2 + bzbc(0)**2 - bxbc(0)**2)
    pg = MAX(pg, none_zero)

    rhog = rhostar(5) + (pg - (pstar(5) - pmagstar(5))) / c0far**2
    rbc(0) = MAX(rhog, none_zero)
    ebc(0) = pg / (gamma - 1.0_num) / rbc(0)

  END SUBROUTINE open_bcs_alfven

END MODULE openboundary


text:
      SUBROUTINE HODOTB( SHYDRA, DOTDEN, LDOT, LLINE, LCEN, LSPIKE,
     &                   STRMAX, STRNOP, STRNON, STRCEN, STRRAD,
     &                   STRLVC, STRLBO, STRSBO, STRBRD, CVECT,
     &                   ATMAX, ATNO, ATXYZ, ATVDW, NOUT, LERR)
      IMPLICIT NONE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Birkbeck College. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Birkbeck College are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 1994 Oliver Smart & Birkbeck College, All rights reserved    *
C * (c) 1996 Oliver Smart & Birkbeck College, All rights reserved    *
C *                                                                  *
C ********************************************************************
C
C Modification history:
C
C Date	Author		Modification
C 12/93	O.S. Smart	Original version
C
C This s/r contains code to output molecular graphics output
C (centre line and dot surface) of HOLE's result.
C This routine does job for sphere-box's (s/r hodotu is the one for sphere).

C passed vbles ****************

C stream number for output hydra/quanta binary plot file
C return unchanged!
      INTEGER			SHYDRA

C a dot density for dot surface
C same definition as in hydra - value between 5 and 30
C default to 10
      INTEGER                   DOTDEN

C
C if are drawing graphical representation of results
C can turn off the dot surface by setting ldot false
C turn off lines joining each sphere edge to
C   the two closest atoms by setting ldot false
C can turn off line joining the sphere centres
      LOGICAL                   LDOT, LLINE, LCEN

C New option April 1993 - instead of drawing dots
C draw a line from dot to surface of atom.
C Turn on by setting lspike .true..
      LOGICAL                   LSPIKE

C store for sphere centres and radii
C as this routine only deals with spheres do not need spher-box
C store the initial point as no 0 and subseqently points in the +ve
C bit. Points below the initial point are stored in -ve bit.
C maximum number of entries
      INTEGER                   STRMAX
C the no. of entries, +Ve entries -Ve entries
      INTEGER                   STRNOP, STRNON
C the centres and radii
      DOUBLE PRECISION          STRCEN(3,-STRMAX:STRMAX),
     &                          STRRAD(-STRMAX:STRMAX)

C vbles to store the dimensions found in the spherebox
C options.
      DOUBLE PRECISION          STRLVC(3,-STRMAX:STRMAX),
     &                          STRLBO(-STRMAX:STRMAX),
     &                          STRSBO(-STRMAX:STRMAX),
     &                          STRBRD(-STRMAX:STRMAX)

C the channel vector
      DOUBLE PRECISION		CVECT(3)

C need atom numbers, co-ords & radius
C maximum no. of atoms
      INTEGER                   ATMAX
C number of atoms we have
      INTEGER                   ATNO
C co-ordinates
      DOUBLE PRECISION          ATXYZ( 3, ATMAX)
C vdw radii of each atom
      DOUBLE PRECISION          ATVDW(ATMAX)

C output stream number
      INTEGER			NOUT
C error indicator
      LOGICAL			LERR

C internal vbles **************

C program uses a sphere of isotropically distributed points
C on a sphere of unit radius. PTMAX is the array bound,
C PTNO is the number of dots (dependent on dotden)
      INTEGER                   PTMAX, PTNO
      PARAMETER(                PTMAX = 4000)
      DOUBLE PRECISION          PTXYZ( 3, PTMAX)

C for drawing move/draws/dots etc. in binary plot file
C n.b. need real rather than d.p's
      REAL                      RVEC4(0:3)

C loop counts
      INTEGER			ICOUNT, JCOUNT, CCOUNT

C centre of sphere box
      DOUBLE PRECISION                  CENTRE(3)

C unit vector giving directions of axes which define sphere box
C first axis is vector (axis(1,1),axis(2,1),axis(3,1)
C N.B. no check is made in this s/r that axes are unit
C and orthogonal - This must be made certain in calling routine.
      DOUBLE PRECISION                  AXIS(3,3)

C radius of sphere box, its additional dimensions along the three axis
      DOUBLE PRECISION                  BOXRAD, BOXDIM(3)

C cloest point to sphere box, another one, distance
      DOUBLE PRECISION			XCLOS(3), DCLOS(3), DIST

C point on sphere
      DOUBLE PRECISION			SPOINT(3)

C maximum radius of sphere-box
      DOUBLE PRECISION			MAXRAD

C end of decs *****************

C line joining centres (just like uniform)
      IF (LCEN) THEN
C change to colour 2
C Change 26/10/93 do in colour 7 - quanta default bright green
        WRITE(SHYDRA) 1.0, 7.0, 0.0, 0.0
C move to first centre
	RVEC4(0) = 2.0
	RVEC4(1) = STRCEN(1,-STRNON)
	RVEC4(2) = STRCEN(2,-STRNON)
	RVEC4(3) = STRCEN(3,-STRNON)
	WRITE(SHYDRA) RVEC4
C draw between the next
	DO 50 ICOUNT = -STRNON+1, STRNOP
	  RVEC4(0) = 3.0
	  RVEC4(1) = STRCEN(1,ICOUNT)
	  RVEC4(2) = STRCEN(2,ICOUNT)
	  RVEC4(3) = STRCEN(3,ICOUNT)
	  WRITE(SHYDRA) RVEC4
50	CONTINUE
      ENDIF

      IF (LDOT) THEN
C load up the set of points =lly distributed in space
C density determined by dotden
        CALL PTGEN( DOTDEN, PTMAX, PTNO, PTXYZ, NOUT, LERR)
        IF (LERR) THEN
          WRITE( NOUT,'(A)')
     &' Have found error in s/r PTGEN',
     &'   Presumably dotden set too high - will try to recover',
     &'   By setting it to 10'
          LERR = .FALSE.
          DOTDEN = 10
          CALL PTGEN( DOTDEN, PTMAX, PTNO, PTXYZ, NOUT, LERR)
          IF (LERR) THEN
            WRITE( NOUT, '(A)')
     &' HAS NOT WORKED! ABORTING'
            GOTO 55555
          ENDIF
        ENDIF

C change to colour 3
C Change 26/10/93 do in colour 8 - quanta default purple
        WRITE(SHYDRA) 1.0, 8.0, 0.0, 0.0
C do not draw dots for the end points
C to avoid bone shaped object
	DO 60 ICOUNT = -STRNON+1, STRNOP-1
C do not plot if radius -ve
	  IF (STRRAD(ICOUNT).GT.0.) THEN
C go through all ptno equally distributed points
	    DO 30 JCOUNT = 1, PTNO
C load up axes for sphere-box (do for each dot
C as  use same vbles to check for overlap).
C axis #1 is cvect
              AXIS(1,1) = CVECT(1)
              AXIS(2,1) = CVECT(2)
              AXIS(3,1) = CVECT(3)
C axis#2 is the long axis
              AXIS(1,2) = STRLVC(1,ICOUNT)
              AXIS(2,2) = STRLVC(2,ICOUNT)
              AXIS(3,2) = STRLVC(3,ICOUNT)
C third axis is the cross product of the other two
              CALL DCROSS( AXIS(1,1), AXIS(1,2), AXIS(1,3))
C just check that this has length 1 - if not then stop
              IF ( ABS(1.-(AXIS(1,3)**2+AXIS(2,3)**2+AXIS(3,3)**2))
     &     .GT. 1E-06) THEN
                WRITE( NOUT, *) 'ERROR in s/r hodotb'
                WRITE( NOUT, *) 'CVECT and LVECT not orthogonal!'
                WRITE( NOUT,*) ' CVECT  = ', CVECT
                WRITE( NOUT,*) ' STRLVC  = ', STRLVC(1,ICOUNT),
     &                       STRLVC(2,ICOUNT), STRLVC(3,ICOUNT)
                LERR = .TRUE.
                GOTO 55555
              ENDIF
C load up box lengths - length along cvect (axis#1) is the long length
              BOXDIM(1) = ABS(STRLBO(ICOUNT))
              BOXDIM(2) = ABS(STRLBO(ICOUNT))
              BOXDIM(3) = ABS(STRSBO(ICOUNT))
              BOXRAD    = ABS(STRBRD(ICOUNT))
C centre of sphere-box
              CENTRE(1) =  STRCEN(1,ICOUNT)
              CENTRE(2) =  STRCEN(2,ICOUNT)
              CENTRE(3) =  STRCEN(3,ICOUNT)

C work out 'maximum' radius of sphere box
C which has box-length strsbo & strlbo twice
              MAXRAD = SQRT( 2.*STRLBO(ICOUNT)**2 +
     &                          STRSBO(ICOUNT)**2   ) +
     &                 STRBRD(ICOUNT)

C We are doing sphere-boxes so consider point closest to
C a point on sphere radius maxrad
	      SPOINT(1) = STRCEN(1,ICOUNT) + MAXRAD* PTXYZ(1,JCOUNT)
	      SPOINT(2) = STRCEN(2,ICOUNT) + MAXRAD* PTXYZ(2,JCOUNT)
	      SPOINT(3) = STRCEN(3,ICOUNT) + MAXRAD* PTXYZ(3,JCOUNT)

C find closest point on sphere-box (XCLOS)
              CALL SBDCLO( DIST, XCLOS, SPOINT,
     &                    AXIS, BOXRAD, BOXDIM, CENTRE)

C o.k. have prospective dot XCLOS - shall we plot it?
C only plot if point is not within neighbouring sphere-box
C have to replace old simple distance calculation with
C calls to sbdclo
C preceeding 20 centres
              DO 35 CCOUNT = ICOUNT-1, MAX(ICOUNT-20,-STRNON), -1
C load up axes for this sphere box
C axis #1 is cvect
                AXIS(1,1) = CVECT(1)
                AXIS(2,1) = CVECT(2)
                AXIS(3,1) = CVECT(3)
C axis#2 is the long axis
                AXIS(1,2) = STRLVC(1,CCOUNT)
                AXIS(2,2) = STRLVC(2,CCOUNT)
                AXIS(3,2) = STRLVC(3,CCOUNT)
C third axis is the cross product of the other two
                CALL DCROSS( AXIS(1,1), AXIS(1,2), AXIS(1,3))
C just check that this has length 1 - if not then stop
                IF ( ABS(1.-(AXIS(1,3)**2+AXIS(2,3)**2+AXIS(3,3)**2))
     &     .GT. 1E-06) THEN
                  WRITE( NOUT, *) 'ERROR in s/r hodotb'
                  WRITE( NOUT, *) 'CVECT and LVECT not orthogonal! (2)'
                  WRITE( NOUT,*) ' CVECT  = ', CVECT
                  WRITE( NOUT,*) ' STRLVC  = ', STRLVC(1,CCOUNT),
     &                       STRLVC(2,CCOUNT), STRLVC(3,CCOUNT)
                  LERR = .TRUE.
                  GOTO 55555
                ENDIF
C load up box lengths - length along cvect (axis#1) is the long length
                BOXDIM(1) = ABS(STRLBO(CCOUNT))
                BOXDIM(2) = ABS(STRLBO(CCOUNT))
                BOXDIM(3) = ABS(STRSBO(CCOUNT))
                BOXRAD =    ABS(STRBRD(CCOUNT))
C centre of sphere-box
                CENTRE(1) =  STRCEN(1,CCOUNT)
                CENTRE(2) =  STRCEN(2,CCOUNT)
                CENTRE(3) =  STRCEN(3,CCOUNT)
C use s/r sbdclo to find the distance from our prospective dot
C to the surface of sphere-box# ccount (dclos is the
C closest point - not used)
                CALL SBDCLO( DIST, DCLOS, XCLOS,
     &                      AXIS, BOXRAD, BOXDIM, CENTRE)

C if the distance is less than zero then should not plot
                IF (DIST.LT.0.) GOTO 30
35            CONTINUE

C following 20 centres
              DO 36 CCOUNT = ICOUNT+1, MIN(ICOUNT+20,STRNOP)
C load up axes for this sphere box
C axis #1 is cvect
                AXIS(1,1) = CVECT(1)
                AXIS(2,1) = CVECT(2)
                AXIS(3,1) = CVECT(3)
C axis#2 is the long axis
                AXIS(1,2) = STRLVC(1,CCOUNT)
                AXIS(2,2) = STRLVC(2,CCOUNT)
                AXIS(3,2) = STRLVC(3,CCOUNT)
C third axis is the cross product of the other two
                CALL DCROSS( AXIS(1,1), AXIS(1,2), AXIS(1,3))
C just check that this has length 1 - if not then stop
                IF ( ABS(1.-(AXIS(1,3)**2+AXIS(2,3)**2+AXIS(3,3)**2))
     &     .GT. 1E-06) THEN
                  WRITE( NOUT, *) 'ERROR in s/r hodotb'
                  WRITE( NOUT, *) 'CVECT and LVECT not orthogonal! (3)'
                  WRITE( NOUT,*) ' CVECT  = ', CVECT
                  WRITE( NOUT,*) ' STRLVC  = ', STRLVC(1,CCOUNT),
     &                       STRLVC(2,CCOUNT), STRLVC(3,CCOUNT)
                  LERR = .TRUE.
                  GOTO 55555
                ENDIF
C load up box lengths - length along cvect (axis#1) is the long length
                BOXDIM(1) = ABS(STRLBO(CCOUNT))
                BOXDIM(2) = ABS(STRLBO(CCOUNT))
                BOXDIM(3) = ABS(STRSBO(CCOUNT))
                BOXRAD =    ABS(STRBRD(CCOUNT))
C centre of sphere-box
                CENTRE(1) =  STRCEN(1,CCOUNT)
                CENTRE(2) =  STRCEN(2,CCOUNT)
                CENTRE(3) =  STRCEN(3,CCOUNT)
C use s/r sbdclo to find the distance from our prospective dot
C to the surface of sphere-box# ccount (dclos is the
C closest point - not used)
                CALL SBDCLO( DIST, DCLOS, XCLOS,
     &                      AXIS, BOXRAD, BOXDIM, CENTRE)

C if the distance is less than zero then should not plot
                IF (DIST.LT.0.) GOTO 30
36            CONTINUE

C write the dot - n.b. spike option not considered for
C sphere-box at present
C indicate dot by making 1st record 4.0
              RVEC4(0) = 4.0
              RVEC4(1) = XCLOS(1)
              RVEC4(2) = XCLOS(2)
              RVEC4(3) = XCLOS(3)
              WRITE(SHYDRA) RVEC4

30          CONTINUE
	  ENDIF
60	CONTINUE
      ENDIF

C sphere-box option draw the long and short axes for each sphere-box
      WRITE(SHYDRA) 1.0, 9.0, 0.0, 0.0
C go thru centre's
      DO 77 ICOUNT = -STRNON, STRNOP
C draw line for long axis
          RVEC4(0) = 2.0
          RVEC4(1) = STRCEN(1,ICOUNT) - STRLBO(ICOUNT)*STRLVC(1,ICOUNT)
          RVEC4(2) = STRCEN(2,ICOUNT) - STRLBO(ICOUNT)*STRLVC(2,ICOUNT)
          RVEC4(3) = STRCEN(3,ICOUNT) - STRLBO(ICOUNT)*STRLVC(3,ICOUNT)
          WRITE(SHYDRA) RVEC4
          RVEC4(0) = 3.0
          RVEC4(1) = STRCEN(1,ICOUNT) + STRLBO(ICOUNT)*STRLVC(1,ICOUNT)
          RVEC4(2) = STRCEN(2,ICOUNT) + STRLBO(ICOUNT)*STRLVC(2,ICOUNT)
          RVEC4(3) = STRCEN(3,ICOUNT) + STRLBO(ICOUNT)*STRLVC(3,ICOUNT)
          WRITE(SHYDRA) RVEC4
C work out short box length axis
C axis #1 is cvect
          AXIS(1,1) = CVECT(1)
          AXIS(2,1) = CVECT(2)
          AXIS(3,1) = CVECT(3)
C axis#2 is the long axis
          AXIS(1,2) = STRLVC(1,ICOUNT)
          AXIS(2,2) = STRLVC(2,ICOUNT)
          AXIS(3,2) = STRLVC(3,ICOUNT)
C third axis is the cross product of the other two
          CALL DCROSS( AXIS(1,1), AXIS(1,2), AXIS(1,3))

C draw line for this
          RVEC4(0) = 2.0
          RVEC4(1) = STRCEN(1,ICOUNT) - STRSBO(ICOUNT)*AXIS(1,3)
          RVEC4(2) = STRCEN(2,ICOUNT) - STRSBO(ICOUNT)*AXIS(2,3)
          RVEC4(3) = STRCEN(3,ICOUNT) - STRSBO(ICOUNT)*AXIS(3,3)
          WRITE(SHYDRA) RVEC4
          RVEC4(0) = 3.0
          RVEC4(1) = STRCEN(1,ICOUNT) + STRSBO(ICOUNT)*AXIS(1,3)
          RVEC4(2) = STRCEN(2,ICOUNT) + STRSBO(ICOUNT)*AXIS(2,3)
          RVEC4(3) = STRCEN(3,ICOUNT) + STRSBO(ICOUNT)*AXIS(3,3)
          WRITE(SHYDRA) RVEC4
77      CONTINUE

55555 RETURN
      END


  MODULE soundcb_module

  implicit none

  private
  public :: soundcb

  CONTAINS

      subroutine soundcb(dt,xh,arh1,arh2,uh,ruh,xf,uf,yh,vh,rvh,yf,vf,    &
                        rds,sigma,rdsf,sigmaf,zh,mh,rmh,c1,c2,mf,zf,      &
                        pi0,rho0,rr0,rf0,rrf0,th0,rth0,zs,                &
                        gz,rgz,gzu,rgzu,gzv,rgzv,                         &
                        dzdx,dzdy,gx,gxu,gy,gyv,                          &
                        radbcw,radbce,radbcs,radbcn,                      &
                        dum1,dum2,dum3,dum4,dum5,dum6,                    &
                        ppd ,fpk ,qk ,pk1,pk2,ftk,                        &
                        u0,rru,ua,u3d,uten,                               &
                        v0,rrv,va,v3d,vten,                               &
                        rrw,wa,w3d,wten,                                  &
                        ppi,pp3d,piadv,ppten,ppx,phi1,phi2,               &
                        thv,ppterm,nrk,dttmp,rtime,mtime,get_time_avg,    &
                        pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2,reqs_p)
      use input
      use constants
      use misclibs , only : convinitu,convinitv,get_wnudge,getdiv
      use bc_module
      use comm_module
      implicit none

      real, intent(in) :: dt
      real, intent(in), dimension(ib:ie) :: xh,arh1,arh2,uh,ruh
      real, intent(in), dimension(ib:ie+1) :: xf,uf
      real, intent(in), dimension(jb:je) :: yh,vh,rvh
      real, intent(in), dimension(jb:je+1) :: yf,vf
      real, intent(in), dimension(kb:ke) :: rds,sigma
      real, intent(in), dimension(kb:ke+1) :: rdsf,sigmaf
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: zh,mh,rmh,c1,c2
      real, intent(in), dimension(ib:ie,jb:je,kb:ke+1) :: mf,zf
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: pi0,rho0,rr0,rf0,rrf0,th0,rth0
      real, intent(in), dimension(ib:ie,jb:je) :: zs
      real, intent(in), dimension(itb:ite,jtb:jte) :: gz,rgz,gzu,rgzu,gzv,rgzv,dzdx,dzdy
      real, intent(in), dimension(itb:ite,jtb:jte,ktb:kte) :: gx,gxu,gy,gyv
      real, intent(inout), dimension(jb:je,kb:ke) :: radbcw,radbce
      real, intent(inout), dimension(ib:ie,kb:ke) :: radbcs,radbcn
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke) :: dum1,dum2,dum3,dum4,dum5,dum6
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke) :: ppd,fpk,qk,pk1,pk2,ftk
      real, intent(in), dimension(ib:ie+1,jb:je,kb:ke) :: u0
      real, intent(inout), dimension(ib:ie+1,jb:je,kb:ke) :: rru,ua,u3d,uten
      real, intent(in), dimension(ib:ie,jb:je+1,kb:ke) :: v0
      real, intent(inout), dimension(ib:ie,jb:je+1,kb:ke) :: rrv,va,v3d,vten
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke+1) :: rrw,wa,w3d,wten
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke) :: ppi,pp3d,piadv,ppten,ppx
      real, intent(inout), dimension(ibph:ieph,jbph:jeph,kbph:keph) :: phi1,phi2
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: thv,ppterm
      integer, intent(in) :: nrk
      real, intent(in)  :: dttmp,rtime
      double precision, intent(in) :: mtime
      logical, intent(in) :: get_time_avg
      real, intent(inout), dimension(jmp,kmp) :: pw1,pw2,pe1,pe2
      real, intent(inout), dimension(imp,kmp) :: ps1,ps2,pn1,pn2
      integer, intent(inout), dimension(rmp) :: reqs_p

!-----

      integer :: i,j,k,n,nloop
      real :: tem,tem1,tem2,tem3,tem4,r1,r2,dts

      real :: temx,temy,tavg

!---------------------------------------------------------------------

    IF( nrkmax.eq.3 )THEN
      if(nrk.eq.1)then
!!!        nloop=1
!!!        dts=dt/3.
        nloop=nint(float(nsound)/3.0)
        dts=dt/(nloop*3.0)
        if( dts.gt.(dt/nsound) )then
          nloop=nloop+1
          dts=dt/(nloop*3.0)
        endif
      elseif(nrk.eq.2)then
        nloop=0.5*nsound
        dts=dt/nsound
      elseif(nrk.eq.3)then
        nloop=nsound
        dts=dt/nsound
      endif
    ELSE
      stop 97394
    ENDIF

!!!      print *,'  nloop,dts,dttmp = ',nloop,dts,nloop*dts

!---------------------------------------------------------------------
!  Arrays for vadv:

      IF(terrain_flag)THEN

!$omp parallel do default(shared)   &
!$omp private(i,j,k,tem,tem1,r1,r2)
            do j=1,nj
            do i=1,ni
              dum3(i,j,1)=0.0
              dum3(i,j,nk+1)=0.0
            enddo
            enddo

      ENDIF

!---------------------------------------------------------------------
!  Prepare for acoustic steps

      if( nrk.eq.1 )then

        !$omp parallel do default(shared)   &
        !$omp private(i,j,k)
        do k=1,nk
        do j=0,nj+1
        do i=0,ni+1
          ppd(i,j,k)=ppx(i,j,k)
        enddo
        enddo
        enddo

      else

        !$omp parallel do default(shared)   &
        !$omp private(i,j,k)
        do k=1,nk
          do j=1,nj
          do i=1,ni+1
            u3d(i,j,k)=ua(i,j,k)
          enddo
          enddo
          IF(axisymm.eq.0)THEN
            ! Cartesian grid:
            do j=1,nj+1
            do i=1,ni
              v3d(i,j,k)=va(i,j,k)
            enddo
            enddo
          ENDIF
          IF(k.ge.2)THEN
            do j=1,nj
            do i=1,ni
              w3d(i,j,k)=wa(i,j,k)
            enddo
            enddo
          ENDIF
          do j=1,nj
          do i=1,ni
            phi2(i,j,k)=phi1(i,j,k)
          enddo
          enddo
          do j=0,nj+1
          do i=0,ni+1
            ppd(i,j,k)=ppx(i,j,k)
          enddo
          enddo
        enddo

      endif

!---------------------------------------------------------------------
!  time-averaged velocities:

      IF( get_time_avg )THEN
!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        DO k=1,nk
          do j=1,nj
          do i=1,ni+1
            rru(i,j,k)=0.0
          enddo
          enddo
          IF(axisymm.eq.0)THEN
            do j=1,nj+1
            do i=1,ni
              rrv(i,j,k)=0.0
            enddo
            enddo
          ENDIF
          IF(k.ge.2)THEN
            do j=1,nj
            do i=1,ni
              rrw(i,j,k)=0.0
            enddo
            enddo
          ENDIF
        ENDDO
      ENDIF
      if(timestats.ge.1) time_sound=time_sound+mytime()

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  Begin small steps:

      small_step_loop:  DO N=1,NLOOP

!-----

        if(irbc.eq.2)then

          if(ibw.eq.1 .or. ibe.eq.1) call radbcew(radbcw,radbce,u3d)

          if(ibs.eq.1 .or. ibn.eq.1) call radbcns(radbcs,radbcn,v3d)

        endif

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  Open boundary conditions:

        IF(wbc.eq.2.and.ibw.eq.1)THEN
          ! west open bc tendency:
          call   ssopenbcw(uh,rds,sigma,rdsf,sigmaf,gz,rgzu,gx,radbcw,dum1,u3d,uten,dts)
        ENDIF
        IF(ebc.eq.2.and.ibe.eq.1)THEN
          ! east open bc tendency:
          call   ssopenbce(uh,rds,sigma,rdsf,sigmaf,gz,rgzu,gx,radbce,dum1,u3d,uten,dts)
        ENDIF

        IF(roflux.eq.1)THEN
          call restrict_openbc_we(rvh,rmh,rho0,u3d)
        ENDIF

!-----

      IF(axisymm.eq.0)THEN
        IF(sbc.eq.2.and.ibs.eq.1)THEN
          ! south open bc tendency:
          call   ssopenbcs(vh,rds,sigma,rdsf,sigmaf,gz,rgzv,gy,radbcs,dum1,v3d,vten,dts)
        ENDIF
        IF(nbc.eq.2.and.ibn.eq.1)THEN
          ! north open bc tendency:
          call   ssopenbcn(vh,rds,sigma,rdsf,sigmaf,gz,rgzv,gy,radbcn,dum1,v3d,vten,dts)
        ENDIF

        IF(roflux.eq.1)THEN
          call restrict_openbc_sn(ruh,rmh,rho0,v3d)
        ENDIF
      ENDIF

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  integrate u,v forward in time:

        if( n.ne.1 )then
          if(timestats.ge.1) time_sound=time_sound+mytime()
          call comm_1p_end(ppd,pw1,pw2,pe1,pe2,   &
                               ps1,ps2,pn1,pn2,reqs_p)
        endif

!-----

    IF(.not.terrain_flag)THEN

      IF(axisymm.eq.0)THEN
        ! Cartesian grid without terrain:

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        do k=1,nk
          do j=1,nj
          do i=1+ibw,ni+1-ibe
            u3d(i,j,k)=u3d(i,j,k)+dts*( uten(i,j,k)      &
                    -rdx*(ppd(i,j,k)-ppd(i-1,j,k))*uf(i) )
          enddo
          enddo
          do j=1+ibs,nj+1-ibn
          do i=1,ni
            v3d(i,j,k)=v3d(i,j,k)+dts*( vten(i,j,k)      &
                    -rdy*(ppd(i,j,k)-ppd(i,j-1,k))*vf(j) )
          enddo
          enddo
        enddo

      ELSE
        ! axisymmetric grid:

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        do k=1,nk
          do j=1,nj
          do i=1+ibw,ni+1-ibe
            u3d(i,j,k)=u3d(i,j,k)+dts*( uten(i,j,k)      &
                    -rdx*(ppd(i,j,k)-ppd(i-1,j,k))*uf(i) )
          enddo
          enddo
        enddo

      ENDIF

    ELSE

        ! Cartesian grid with terrain:

!$omp parallel do default(shared)   &
!$omp private(i,j,k,r1,r2)
        do j=0,nj+1
          do k=2,nk
          do i=0,ni+1
            dum1(i,j,k) = (ppd(i,j,k)-ppd(i,j,k-1))*rds(k)
          enddo
          enddo
          do i=0,ni+1
            dum1(i,j,1) = 0.0
            dum1(i,j,nk+1) = 0.0
          enddo
        enddo

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        do k=1,nk
          ! x-dir
          do j=1,nj
          do i=1+ibw,ni+1-ibe
            u3d(i,j,k)=u3d(i,j,k)+dts*( uten(i,j,k)         &
                  -( (ppd(i,j,k)-ppd(i-1,j,k))*rdx*uf(i)    &
              +0.125*( (dum1(i,j,k+1)+dum1(i-1,j,k+1))      &
                      +(dum1(i,j,k  )+dum1(i-1,j,k  )) )    &
                    *(gxu(i,j,k)+gxu(i,j,k+1))    ) )
          enddo
          enddo
          do j=1+ibs,nj+1-ibn
          do i=1,ni
            v3d(i,j,k)=v3d(i,j,k)+dts*( vten(i,j,k)         &
                  -( (ppd(i,j,k)-ppd(i,j-1,k))*rdy*vf(j)    &
              +0.125*( (dum1(i,j,k+1)+dum1(i,j-1,k+1))      &
                      +(dum1(i,j,k  )+dum1(i,j-1,k  )) )    &
                    *(gyv(i,j,k)+gyv(i,j,k+1))    ) )
          enddo
          enddo
        enddo

    ENDIF

!----------------------------------------------
!  convergence forcing:

        IF( convinit.eq.1 )THEN
          IF( rtime.le.convtime .and. nx.gt.1 )THEN
            call convinitu(myid,ib,ie,jb,je,kb,ke,ni,nj,nk,ibw,ibe,   &
                           zdeep,lamx,lamy,xcent,ycent,aconv,    &
                           xf,yh,zh,u0,u3d)
          ENDIF
        ENDIF

!----------------------------------------------
!  convergence forcing:

      IF(axisymm.eq.0)THEN
        ! Cartesian grid:
        IF( convinit.eq.1 )THEN
          IF( rtime.le.convtime .and. ny.gt.1 )THEN
            call convinitv(myid,ib,ie,jb,je,kb,ke,ni,nj,nk,ibs,ibn,   &
                           zdeep,lamx,lamy,xcent,ycent,aconv,    &
                           xh,yf,zh,v0,v3d)
          ENDIF
        ENDIF

!----------------------------------------------

      ENDIF

      if(timestats.ge.1) time_sound=time_sound+mytime()

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  integrate w forward in time:

      IF( wnudge.eq.1 )THEN
        !  updraft nudging tendency:
        IF( (mtime+dt).le.t2_wnudge )THEN
          !$omp parallel do default(shared)   &
          !$omp private(i,j,k)
          do k=2,nk
          do j=1,nj
          do i=1,ni
            dum1(i,j,k)=0.0
          enddo
          enddo
          enddo
          call get_wnudge(mtime,dts,xh,yh,zf,w3d,dum1)
        ENDIF
      ENDIF

      IF(.not.terrain_flag)THEN
        ! without terrain:

!$omp parallel do default(shared)   &
!$omp private(i,j,k,tem1)
        do k=2,nk
        tem1 = rdz*mf(1,1,k)
        do j=1,nj
        do i=1,ni
          w3d(i,j,k)=w3d(i,j,k)+dts*( wten(i,j,k)            &
                  -tem1*(ppd(i,j,k)-ppd(i,j,k-1)) )
        enddo
        enddo
        enddo

      ELSE
        ! with terrain:

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        do k=2,nk
        do j=1,nj
        do i=1,ni
          w3d(i,j,k)=w3d(i,j,k)+dts*( wten(i,j,k)                     &
                -rds(k)*(ppd(i,j,k)-ppd(i,j,k-1))*gz(i,j) )
        enddo
        enddo
        enddo
        if(timestats.ge.1) time_sound=time_sound+mytime()

        call bcwsfc(gz,dzdx,dzdy,u3d,v3d,w3d)

      ENDIF

      IF( wnudge.eq.1 )THEN
        !  apply updraft nudging:
        IF( (mtime+dt).le.t2_wnudge )THEN
          !$omp parallel do default(shared)   &
          !$omp private(i,j,k)
          do k=2,nk
          do j=1,nj
          do i=1,ni
            w3d(i,j,k)=w3d(i,j,k)+dum1(i,j,k)
          enddo
          enddo
          enddo
        ENDIF
      ENDIF

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  get divergence:

        call     getdiv(arh1,arh2,uh,vh,mh,u3d,v3d,w3d,dum1,dum2,dum3,dum4,  &
                        rds,rdsf,sigma,sigmaf,gz,rgzu,rgzv,dzdx,dzdy)

!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  get new pp,th

      tem = dts*csound*csound

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
      DO k=1,nk
        do j=1,nj
        do i=1,ni
          ppd(i,j,k)=phi2(i,j,k)
          phi2(i,j,k)=phi2(i,j,k)-tem*dum4(i,j,k)
          if(abs(phi2(i,j,k)).lt.smeps) phi2(i,j,k)=0.0
          dum1(i,j,k)=kdiv*( phi2(i,j,k)-ppd(i,j,k) )
          ppd(i,j,k)=phi2(i,j,k)+dum1(i,j,k)
        enddo
        enddo
      ENDDO
      if(timestats.ge.1) time_sound=time_sound+mytime()

        IF( n.lt.nloop )THEN
          call bcs(ppd)
          call comm_1s_start(ppd,pw1,pw2,pe1,pe2,   &
                                 ps1,ps2,pn1,pn2,reqs_p)
        ENDIF

!--------------------------------------------------------------------
!  time-averaged velocities:

      IF( get_time_avg )THEN

      if( n.lt.nloop )then
!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        DO k=1,nk
          do j=1,nj
          do i=1,ni+1
            rru(i,j,k)=rru(i,j,k)+u3d(i,j,k)
          enddo
          enddo
          IF( axisymm.eq.0 )THEN
            do j=1,nj+1
            do i=1,ni
              rrv(i,j,k)=rrv(i,j,k)+v3d(i,j,k)
            enddo
            enddo
          ENDIF
          IF( k.ge.2 )THEN
            do j=1,nj
            do i=1,ni
              rrw(i,j,k)=rrw(i,j,k)+w3d(i,j,k)
            enddo
            enddo
          ENDIF
        ENDDO
      else
        tavg = 1.0/float(nloop)
!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        DO k=1,nk
          do j=1,nj
          do i=1,ni+1
            rru(i,j,k)=(rru(i,j,k)+u3d(i,j,k))*tavg
          enddo
          enddo
          IF( axisymm.eq.0 )THEN
            do j=1,nj+1
            do i=1,ni
              rrv(i,j,k)=(rrv(i,j,k)+v3d(i,j,k))*tavg
            enddo
            enddo
          ENDIF
          IF( k.ge.2 )THEN
            do j=1,nj
            do i=1,ni
              rrw(i,j,k)=(rrw(i,j,k)+w3d(i,j,k))*tavg
            enddo
            enddo
          ENDIF
        ENDDO
      endif

      ENDIF

!--------------------------------------------------------------------

      ENDDO  small_step_loop

!  end of small steps
!--------------------------------------------------------------------

      if(timestats.ge.1) time_sound=time_sound+mytime()
      call bcs(phi2)
      call comm_1s_start(phi2,pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2,reqs_p)
      call comm_1s_end(  phi2,pw1,pw2,pe1,pe2,ps1,ps2,pn1,pn2,reqs_p)

      IF( nrk.eq.nrkmax )THEN

        ! pressure tendency term: save for next timestep:
        ! 170823:  also, step forward phi

!$omp parallel do default(shared)   &
!$omp private(i,j,k)
        do k=1,nk
          do j=1,nj
          do i=1,ni
            ppx(i,j,k)=dum1(i,j,k)
          enddo
          enddo
          do j=0,nj+1
          do i=0,ni+1
            phi1(i,j,k)=phi2(i,j,k)
          enddo
          enddo
        enddo

      ENDIF
      if(timestats.ge.1) time_sound=time_sound+mytime()

      end subroutine soundcb

  END MODULE soundcb_module


text:
program main
  implicit none
  intrinsic:: sum, product, maxval, minval
  integer, dimension(0:3) :: vec
  integer, dimension(0:2, 0:2) :: arr
  integer :: vecSum, vecMax, vecProd, vecMin, i

  vec(0) = 2
  vec(1) = 3
  vec(2) = 4
  vec(3) = 6

  vecSum = 0
  do i = 0, size(vec) - 1
    vecSum = vecSum + vec(i)
  end do

  vecProd = 1 ! do not start with 0
  do i = 0, size(vec) - 1
    vecProd = vecProd * vec(i)
  end do

  vecMax = vec(0)
  do i = 1, size(vec) - 1
    if (vecMax < vec(i)) then
      vecMax = vec(i)
    end if
  end do

  vecMin = vec(0)
  do i = 1, size(vec) - 1
    if (vecMin > vec(i)) then
      vecMin = vec(i)
    end if
  end do

  write (*, *) "Sum of all elements (intrinsic)", sum(vec)
  write (*, *) "Sum of all elements (own)", vecSum
  write (*, *) "Product of all elements (intrinsic)", product(vec)
  write (*, *) "Product of all elements (own)", vecProd
  write (*, *) "Minimal: ", minval(vec), "Maximal value: ", maxval(vec)
  write (*, *) "Minimal: ", vecMin, "Maximal value: ", vecMax
  write (*, *) ""

  ! 1 2 3
  ! 4 5 6
  ! 7 8 9
  arr(0, 0) = 1
  arr(0, 1) = 2
  arr(0, 2) = 3
  arr(1, 0) = 4
  arr(1, 1) = 5
  arr(1, 2) = 6
  arr(2, 0) = 7
  arr(2, 1) = 8
  arr(2, 2) = 9
  write (*, *) 'Full matrix ', arr ! goes down cols
  write (*, *) "Sum of first col", sum(arr(0:2, 0))
end program main

max_stars_repo_name: marvinfriede/projects

max_issues_repo_name: marvinfriede/projects

max_forks_repo_name: marvinfriede/projects