Dear Friends,
I can’t find the error in the first subroutine MaxMin.
The output for maxmom and minmom is OK.
The output for maxposx, maxposa, minposx, minposa stays the same for whatever value of P1, P2, P3, P4.
I’m no software guru but a retired structural engineer doing this as a hobby so don’t be to hard on me. 
! DLL Max. moments and shear forces in 2-span beam with equal lengths
! under up to 4 moving loads
! Version Author Date
! A R.M. 14/06/2024
! B R.M. 02/05/2025
module mom_m
implicit none
save
contains
! Subroutine to calculate moments maxmin
Subroutine MaxMin(L,P1,P2,P3,P4,a1,a2,a3,t,ch1,ar)
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: MaxMin
!DIR$ ATTRIBUTES ALIAS: 'MaxMin' :: Maxmin
!DIR$ ATTRIBUTES REFERENCE :: ch1, ar
real, intent(in) :: L, P1, P2, P3, P4, a1, a2, a3, t
real, intent(out) ::ar(6), ch1(2)
integer :: n
real, dimension(:,:), allocatable :: moments
real :: maxmom
real :: minmom
integer :: maxpos(2)
integer :: minpos(2)
real :: maxposx
real :: maxposa
real :: minposx
real :: minposa
integer :: ia, ix
n = int(2*L/t + 1)
! initialize
ar = 0.0
maxpos = 0
minpos = 0
maxposx = 0.0
maxposa = 0.0
minposx = 0.0
minposa = 0.0
allocate (moments(n,n))
do ix = 1, n
do ia = 1, n
moments(ia, ix) = Mx((ix-1)*t, (ia-1)*t, L, P1) + &
Mx((ix-1)*t, (ia-1)*t - a1, L, P2) + &
Mx((ix-1)*t, (ia-1)*t - (a1+a2), L, P3) + &
Mx((ix-1)*t, (ia-1)*t - (a1+a2+a3), L, P4)
end do
end do
maxmom=maxval(moments)
minmom=minval(moments)
! maxpos=maxloc(moments,mask=moments.gt.0)
maxpos = findloc(moments, maxmom)
! minpos=minloc(moments,mask=moments.lt.0)
minpos = findloc(moments, minmom)
maxposx=(maxpos(2)-1) *t
maxposa=(maxpos(1)-1) *t
minposx=(minpos(2)-1) *t
minposa=(minpos(1)-1) *t
ar= (/maxmom, minmom, maxposx, maxposa, minposx, minposa/)
ch1 = minpos
return
end subroutine MaxMin
function DkEnd(L, P1, P2, P3, P4,a1, a2, a3, t)
! Function to calculate max. shear force at end support
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT:: DkEnd
!DIR$ ATTRIBUTES ALIAS: 'DkEnd' :: DkEnd
real :: DkEnd
integer :: ia
integer :: n
real, intent(in) :: L,P1,P2,P3,P4,a1,a2,a3,t
real:: MaxVa
real, dimension(:), allocatable :: dka
n = int(2*L/t + 1)
allocate (dka(n))
do ia = 1, n
dka(ia) = Va((ia-1)*t, L, P1) + &
Va((ia-1)*t - a1, L, P2) + &
Va((ia-1)*t - (a1+a2), L, P3) + &
Va((ia-1)*t - (a1+a2+a3), L, P4)
end do
MaxVa = minval(dka)
DkEnd=MaxVa
end function DkEnd
function DkMid(L, P1, P2, P3, P4,a1, a2, a3, t)
! Function to calculate max. shear force at mid support
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: DkMid
!DIR$ ATTRIBUTES ALIAS: 'DkMid' :: DkMid
real :: DkMid
integer :: ia
integer :: n
real, intent(in) :: L,P1,P2,P3,P4,a1,a2,a3,t
real MaxVb
real, dimension(:), allocatable :: dkb
n = int(2*L/t + 1)
allocate (dkb(n))
do ia = 1, n
dkb(ia) = Vb((ia-1)*t, L, P1) + &
Vb((ia-1)*t - a1, L, P2) + &
Vb((ia-1)*t - (a1+a2), L, P3) + &
Vb((ia-1)*t - (a1+a2+a3), L, P4)
end do
MaxVb = minval(dkb)
DkMid=MaxVb
end function DkMid
real function ReacB(L, P1, P2, P3, P4,a1, a2, a3, t)
! Function to calculate max. reaction at mid support
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: ReacB
!DIR$ ATTRIBUTES ALIAS: 'ReacB' :: ReacB
integer :: ia
integer :: n
real, intent(in) :: L,P1,P2,P3,P4,a1,a2,a3,t
real MaxRb
real, dimension(:), allocatable :: recb
n = int(2*L/t + 1)
allocate (recb(n))
do ia = 1, n
recb(ia) = Rb((ia-1)*t, L, P1) + &
Rb((ia-1)*t - a1, L, P2) + &
Rb((ia-1)*t - (a1+a2), L, P3) + &
Rb((ia-1)*t - (a1+a2+a3), L, P4)
end do
MaxRb = maxval(recb)
ReacB=MaxRb
end function ReacB
real function Mx(x,a,L,P)
! Function to calcute moment under one moving point load
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: Mx
!DIR$ ATTRIBUTES ALIAS: 'Mx' :: Mx
!Calling parameters
real :: x ! Position where moment is calculated
real :: a ! Position point load
real :: L ! Span
real :: P ! Value load
real :: b
! Evaluate Mx
if (a < L) then
b = L-a
else
a = 2.*L-a
b = L-a
x = 2.*L-x
endif
if (a < 0 .or. a > 2.*L) then
Mx = 0.0
else if (x <= a) then
Mx = ((P * b) / (4. * (L ** 3.))) * ((4. * (L ** 2.)) - (a * (L + a))) * x
else if ((x > a) .and. (x <= L)) then
Mx = ((P * b) / (4. * (L ** 3.))) * ((4. * (L ** 2.)) - (a * (L + a))) * x - P * (x - a)
else if (x > L) then
Mx = -(((P * b * a) / (4 * L ** 3.)) * (L + a)) * (2. * L - x)
end if
end function Mx
real function Va(a,L,P)
! Function to calcute max. shear at end support under one moving point load
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: Va
!DIR$ ATTRIBUTES ALIAS: 'Va' :: Va
!Calling parameters
real :: a, L, P, b ! Position point load, Span, load
Va = 0.0
! Evaluate Va
if (a <= L) then
b = L-a
else
a = 2*L-a
b = a-L
endif
if (a < 0 .or. a > 2*L) then
Va = 0.0
else if (a <= L) then
Va = -((P*b)/(4*L**3))*(4*L**2-a*(L+a))
else if (a > L) then
Va = -((P*a*b)/(4*L**3))*(L+a)
end if
end function Va
real function Vb(a,L,P)
! Function to calcute max. shear at mid support under one moving point load
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: Vb
!DIR$ ATTRIBUTES ALIAS: 'Vb' :: Vb
!Calling parameters
real :: a, L, P, b
Vb = 0.0
! Evaluate Vb
if (a <= L) then
b = L-a
else
a = 2*L-a
b = a-L
endif
if (a < 0 .or. a > 2*L) then
Vb = 0.0
else if (a <= L) then
Vb = -((P*a)/(4*L**3))*(4*L**2+b*(L+a))
else if (a > L) then
Vb = -((P*a*b)/(4*L**3))*(L+a)
end if
end function Vb
real function Rb(a,L,P)
!Function to calcute max. reaction at mid support under one moving point load
!DIR$ ATTRIBUTES STDCALL, DLLEXPORT :: Rb
!DIR$ ATTRIBUTES ALIAS: 'Rb' :: Rb
!Calling parameters
real :: a, L, P, b
Rb = 0.0
! Evaluate rb
if (a <= L) then
b = L-a
else
a = 2*L-a
b = a-L
endif
if (a < 0 .or. a > 2*L) then
Rb = 0.0
else if (a <= L) then
Rb = ((P*a)/(2*L**3))*(2*L**2+b*(L+a))
else if (a > L) then
Rb = ((P*a)/(2*L**3))*(2*L**2+b*(L+a))
end if
end function Rb
subroutine poutrel(name, arp)
!DIR$ ATTRIBUTES REFERENCE, STDCALL, DLLEXPORT :: poutrel
!DIR$ ATTRIBUTES ALIAS: 'poutrel' :: poutrel
!DIR$ ATTRIBUTES MIXED_STR_LEN_ARG :: poutrel
!DIR$ ATTRIBUTES NOINLINE :: procedure
use, intrinsic :: iso_fortran_env, only: error_unit
implicit none
! Calling parameter
character(len=*), intent(in) :: name ! Keep as INTENT(IN) since we won't modify it
real, intent(out) :: arp(26)
type :: profiel
character(10) :: name
real :: A, h, b, tw, tf, r
end type profiel
integer, parameter :: NPROF = 90
type(profiel) :: profielen(NPROF)
! The table is generated automatically from poutr.csv and
! included using a data statement
include "profielen.fi"
real :: A = 0.0, h = 0.0, b = 0.0, tw = 0.0, tf = 0.0, r = 0.0
real :: Iy = 0.0, Avz = 0.0, Wy = 0.0, Wply = 0.0, iry = 0.0
real :: Iz = 0.0, Wz = 0.0, Wplz = 0.0, irz = 0.0
real :: Ss = 0.0, IT = 0.0, Iw = 0.0, G = 0.0
real :: Auf = 0.0, Izuf = 0.0, irzuf = 0.0, Wzuf = 0.0, Wzpluf = 0.0, Sy = 0.0, Syr = 0.0
real :: Mplyd = 0.0, Mplzd = 0.0, Mplwd = 0.0, Vplzd = 0.0, Vplyd = 0.0
real, parameter :: pi = 3.1415926535
character(len=10) :: name_to_find
integer :: idx
! Initialize output array to zeros
arp = 0.0
! Initialize variables
name_to_find = trim(name) ! Use a local copy of the input parameter
! Debug output
print *, "Looking for profile: '", trim(name_to_find), "'"
! Search for the profile
idx = findloc(profielen%name, name_to_find, dim=1)
if (idx == 0) then
print *, "Invalid input: Profile not found."
! Initialize with zeros in case profile not found
arp = 0.0
return
end if
! Profile was found
A = profielen(idx)%A
h = profielen(idx)%h
b = profielen(idx)%b
tw = profielen(idx)%tw
tf = profielen(idx)%tf
r = profielen(idx)%r
! Shear area
Avz = A-2*b*tf+(tw+2*r)*tf
! Weight per meter
G = A*8000*1.0e-6 ! Fixed precision issue
! Second moment of inertia
Iy =(1.0/12.0)*(b*h**3.0-(b-tw)*(h-2.0*tf)**3.0)+0.03*r**4.0+0.2146*r**2.0*(h-2.0*tf-0.4468*r)**2.0
! Second moment of inertia about z-axis
Iz = (1.0/12.0)*(2*tf*b**3.0+(h-2*tf)*tw**3.0) + &
0.03*r**4.0 + 0.2146*r**2.0*(tw-0.4468*r)**2.0
! Radius of gyration
iry = sqrt(Iy/A)
irz = sqrt(Iz/A)
! Tortional constant
IT = (2.0/3.0)*(b-0.63*tf)*tf**3.0+1.0/3.0*(h-2*tf)*tw**3.0+2*tw/tf*(0.145+0.1*r/tf) &
*(((r+tw/2.0)**2.0+(r+tf)**2.0-r**2.0)/(2*r+tf))**4.0
! Warping constant
Iw = ((tf*b**3.0)/24.0)*(h-tf)**2.0
! Length of stiff bearing
Ss = tw+2*tf+(4.0-2.0*sqrt(2.0))*r
! Elastic section modulus
Wy = 2*Iy/h
Wz = 2*Iz/b
! Plastic section modulus
Wply = (tw*h**2.0)/4.0+(b-tw)*(h-tf)*tf+((4.0-pi)/2.0)*r**2.0*(h-2.0*tf)+((3.0*pi-10.0)/3.0)*r**3.0
Wplz = (b**2.0*tf)/2.0 + ((h-2.0*tf)/4.0)*tw**2.0+r**3.0*(10.0/3.0-pi)+(2.0-pi/2.0)*tw*r**2.0
! Static radius
Sy = b*tf*((h-tf)/2.0)+(h/2.0-tf)*tw*((h-2.0*tf)/4.0)+((r**2*(4.0-pi))/4.0)*((h-2.0*tf-r)/2.0)
Syr = Sy - tw*((h/2.0-tf-r)**2)/2.0
! Cross sectional values top flange top flange + 1/5 web)
Auf = b*tf + ((h - 2*tf)/5)*tw
Izuf = Iz/2.0
irzuf = sqrt(Izuf/Auf)
Wzuf = Wz /2.0
Wzpluf = Wplz/2.0
! Output (indexed for easier reference)
arp( 1) = A
arp( 2) = h
arp( 3) = b
arp( 4) = tw
arp( 5) = tf
arp( 6) = r
arp( 7) = Avz
arp( 8) = G
arp( 9) = Iy
arp(10) = Iz
arp(11) = iry
arp(12) = irz
arp(13) = IT
arp(14) = Iw
arp(15) = Ss
arp(16) = Wy
arp(17) = Wz
arp(18) = Wply
arp(19) = Wplz
arp(20) = Auf
arp(21) = Izuf
arp(22) = irzuf
arp(23) = Wzuf
arp(24) = Wzpluf
arp(25) = Sy
arp(26) = Syr
end subroutine poutrel
! Fixed subroutine to properly use profile data
subroutine mplst(name, name_len, Fy, pl)
!DEC$ ATTRIBUTES DLLEXPORT :: mplst
!DIR$ ATTRIBUTES ALIAS: 'mplst' :: mplst
!DEC$ ATTRIBUTES STDCALL :: mplst
!DEC$ ATTRIBUTES REFERENCE :: name
character(len=32), intent(in) :: name ! or 64, or whatever max you need
integer, intent(in) :: name_len
real, intent(in) :: Fy
real, intent(out) :: pl(5)
! Local variables - must declare all at beginning
real :: arp(26)
real :: Wply, Wplz, h, tf, b, Avz, Auf
real :: Mplyd, Mplzd, Mplwd, Vplzd, Vplyd
print *, "mplst called with name: ", trim(name)
print *, "name_len: ", name_len
print *, "Fy: ", Fy
! First get the profile data
call poutrel(name(1:name_len), arp)
if (any(arp == 0.0)) then
print *, "Warning: poutrel returned zero values"
end if
! Get values from the arp array
Wply = arp(18)
Wplz = arp(19)
h = arp(2)
tf = arp(5)
b = arp(3)
Avz = arp(7)
Auf = arp(20)
! Full plastic capacity profile calculations
Mplyd = Wply * Fy
Mplzd = Wplz * Fy
Mplwd = 0.25 * ((h-tf) * b**2 * tf * Fy)
Vplzd = (Avz * Fy) / sqrt(3.0)
Vplyd = (2 * Auf * Fy) / sqrt(3.0)
pl(1) = Mplyd
pl(2) = Mplzd
pl(3) = Mplwd
pl(4) = Vplzd
pl(5) = Vplyd
end subroutine mplst
real function test_dll(x)
!DEC$ ATTRIBUTES DLLEXPORT :: test_dll
!DIR$ ATTRIBUTES ALIAS: 'test_dll' :: test_dll
!DEC$ ATTRIBUTES STDCALL :: test_dll
real, intent(in) :: x
test_dll = x * 2.0
end function test_dll
end module mom_m
! A program to compute moments using the module
program mom
! Need to compile the module first before using it
use mom_m
! Moments in a 2 span continuous beam with equal spans under up to 4 moving point loads.
! 20/05/2024 R.M. Versie A
implicit none
! Declare variables
real :: a1, a2, a3 ! Distance between point loads
real :: P1, P2, P3, P4 ! Point loads
real :: L ! Length span
real :: Fy
real, parameter :: t = 0.05 ! Distance between studied points
real :: ar(6), ch1(2)
real :: arp(26), pl(5) ! Arrays for profile data and plastic capacities
real :: maxposx
character(len=10) :: name = ""
integer :: name_len ! Declare name_len at the beginning
! Retrieve data of European parallel flange beams
! Version Author Date
! A R.M. 21/04/2025
! Test values
name = 'HEB300'
L = 6.0
a1 = 3.6
a2 = 1.0
a3 = 1.0
P1 = 22.5
P2 = 22.5
P3 = 0.0
P4 = 0.0
Fy = 355.0
! Call MaxMin subroutine
call MaxMin(L, P1, P2, P3, P4, a1, a2, a3, t, ch1, ar)
write(*,*) 'Max Min: maxmom, minmom, maxposx, maxposa, minposx, minposa ='
write(*,*) ar
write(*,*) 'Max Va=', DkEnd(L, P1, P2, P3, P4, a1, a2, a3, t)
write(*,*) 'Max Vb=', DkMid(L, P1, P2, P3, P4, a1, a2, a3, t)
write(*,*) 'Max Rb=', ReacB(L, P1, P2, P3, P4, a1, a2, a3, t)
! Get beam profile values
call poutrel(name, arp)
print *, "maxpos= " , maxposx
print *, "Profiel= ", name
print *, "A= ", arp(1)
print *, "h= ", arp(2)
print *, "b= ", arp(3)
print *, "tw= ", arp(4)
print *, "tf= ", arp(5)
print *, "r= ", arp(6)
print *, "Avz= ", arp(7)
print *, "G= ", arp(8)
print *, "Iy= ", arp(9)
print *, "Iz= ", arp(10)
print *, "iy= ", arp(11)
print *, "iz= ", arp(12)
print *, "IT= ", arp(13)
print *, "Iw= ", arp(14)
print *, "Ss= ", arp(15)
print *, "Wy= ", arp(16)
print *, "Wz= ", arp(17)
print *, "Wply= ", arp(18)
print *, "Wplz= ", arp(19)
print *, "Auf= ", arp(20)
print *, "Izuf= ", arp(21)
print *, "irzuf= ", arp(22)
print *, "Wzuf= ", arp(23)
print *, "Wpluf= ", arp(24)
print *, "Sy= ", arp(25)
print *, "Syr= ", arp(26)
! Calculate plastic capacities
name_len = len_trim(name) ! Calculate the length of the name string
call mplst(name(1:name_len), name_len, Fy, pl)
print *, "Mplyd= ", pl(1)
print *, "Mplzd= ", pl(2)
print *, "Mplwd= ", pl(3)
print *, "Vplzd= ", pl(4)
print *, "Vplyd= ", pl(5)
end program mom