Problem compiling code and creating executables using namelist

Hi everyone :smiley:

I currently try to simply a programm calculating low-Prandl number convection using derived types. For the problem I need to take the inputs variables from a external file which I do with a subroutine read_inputs and namelist. I however now unable to create an executable even though the compiler does not indicate an error. I get the following message on the terminal:

/usr/bin/ld:lowPrandt_parameters_FULL.txt: file format not recognized; treating as linker script
/usr/bin/ld:lowPrandt_parameters_FULL.txt:2: syntax error
collect2: error: ld returned 1 exit status
./run: line 4: ./a.out: No such file or directory

Where I was executing the following commands on linux:

gfortran -o ex1 -fbacktrace -g -fcheck=all -fbounds-check -finit-real=snan -ffpe-trap=invalid,zero,overflow -O0 -O -Wall -g  modul_Low_Prandl.f90 exercise_11.f90 lowPrandt_parameters_FULL.txt
./a.out

Do you know how I can resolve this and what the problem is? Also are the used subroutines and functions like second_derivative in the module private? Or if not how can I make sure that they are?

I currently try to run this .f90 file:

program LowPrandl
    use module_LowPrandl
    use iso_fortran_env
    implicit none
    type(grid_variables)    :: all
    real(real32)            :: time
    INTEGER                 :: step

    call all%read_inputs()
    call all%initialise()

    time    = 0.
    step    = 0

    do while (time<all%total_time)
        call all%VelocitySolve()
        call all%Timestep()
        time = time + all%dt
        step = step + 1
        print*, 'Time, step, maxUx, maxUy:', time, step, maxval(abs(all%ux)), maxval(abs(all%uy))  
        if((mod(step,100)==1) .or.( time>=all%total_time)) call all%write(step, time)
    end do
end program LowPrandl

using this module:

module module_LowPrandl
    use iso_fortran_env
    implicit none
    private
    public :: grid_variables
    real(real64)        :: Pi =  atan(1.)*4.

    type :: grid_variables
        !input parameters
        real(real64)        :: total_time = 0.1                     !total integration time
        character(len=10)   :: Tinit='cosine'                       !Temperature initial condition
        integer             :: nx, ny                               !grid spacing
        real(real64)        :: a_dif, a_adv                         !time step constant
        real(real64)        :: err = 1e-6                           !maximum error in Poisson solves = resiude_rms/rhs_rms<err
        real(real64)        :: Ra                                   !Rayleigh number
        real(real64)        :: Pr                                   !Prandt number
        !character(len=100)  :: fname = 'parameter.txt'              !read command-line arguments

        !Variables in type
        real(real64), allocatable   :: T(:,:),dTdx(:,:), d2Tdxdy(:,:),d2Wdxdy(:,:)
        real(real64), allocatable   :: ux(:,:), uy(:,:)
        real(real64), allocatable   :: S(:,:), ugradT(:,:), ugradW(:,:), W(:,:)
        real(real64)                :: dx, dy, h, dt = 0.1, dt_dif = 0.1, dt_adv = 0.1

      contains
        procedure,public,pass(this) :: read_inputs
        procedure,public,pass(this) :: initialise
        procedure,public,pass(this) :: VelocitySolve
        procedure,public,pass(this) :: Timestep
        procedure,public,pass(this) :: write

  end type grid_variables
contains
!the maine subroutines to handle the translation of generic functions for the specific need -> should be public
subroutine read_inputs(this)
    class(grid_variables), INTENT(inout)   :: this
    integer             :: tmp_nx, tmp_ny         
    real(real64)        :: tmp_a_dif, tmp_a_adv, tmp_err,tmp_Ra, tmp_Pr, tmp_total_time
    character(len=100)  :: tmp_fname
    character(len=10)   :: tmp_Tinit

    namelist / inputs / tmp_nx,tmp_ny,tmp_total_time,tmp_Ra,tmp_err,tmp_a_dif,tmp_a_adv,tmp_Tinit,tmp_Pr
    !read command-line arguments
    if ( command_argument_count() > 0 ) then
        call get_command_argument(1,tmp_fname)
    end if

    !inputs is the name of the list
    OPEN(1,file=tmp_fname, status ='old')
    READ(1, inputs) 
    CLOSE(1)

    this%nx           = tmp_nx
    this%ny           = tmp_ny
    this%total_time   = tmp_total_time
    this%Ra           = tmp_Ra
    this%err          = tmp_err
    this%a_dif        = tmp_a_dif
    this%a_adv        = tmp_a_adv
    this%Tinit        = tmp_Tinit
    this%Pr           = tmp_Pr
end subroutine read_inputs

subroutine initialise(this)
    class(grid_variables), INTENT(inout)   :: this
    integer                                :: counter
    allocate(this%T(this%nx,this%ny), this%d2Tdxdy(this%nx, this%ny), this%ugradT(this%nx, this%ny))
    allocate(this%S(this%nx, this%ny), this%ux(this%nx,this%ny), this%uy(this%nx, this%ny),this% W(this%nx,this%ny))
    allocate(this%ugradW(this%nx,this%ny), this%d2Wdxdy(this%nx,this%ny), this%dTdx(this%nx,this%ny))
    !initalize stuff 
    this%S       = 0.
    this%W       = 0.
    this%ux      = 0.
    this%uy      = 0.
    this%ugradT  = 0.
    this%d2Tdxdy = 0.
    this%d2Wdxdy = 0.
    !spatial mesh distance
    this%dy = 1./(this%ny-1)    
    this%dx = this%dy           !h=dy
    this%h =  min(this%dx,this%dy)
    !numerical computation goes here
    if (this%Tinit=='random') then
        call random_number(this%T)
        call random_number(this%W)
        this%T(1,:)         = 0
        this%T(this%nx,:)   = 0
        this%T(:,1)         = 0   
        this%T(:,this%ny)   = 0   
    else if (this%Tinit == 'spike') then
        this%T                      = 0.
        this%T(this%nx/2,this%ny/2) = 1
    else if ( this%Tinit=='cosine' ) then
        do counter = 1, this%nx
            this%T(counter,:) =  0.5*(1+cos(3*Pi*(counter-1)/(this%nx-1)))
        end do
    else
        print*, 'ERROR! Tinit must be either ^spike^, ^cosine^ or ^random^!'
        stop
    end if
    !diffusive time step
    this%dt_dif = abs(this%a_dif*min(this%dx, this%dy)**2/max(this%Pr,1.))
end subroutine initialise

subroutine VelocitySolve(this)
    class(grid_variables), INTENT(inout)   :: this
    this%S = solve_Poisson_Eq_multigrid(this%S, this%W, This%h,this%err)
    call S_to_vel(this%S, this%ux, this%uy, this%nx, this%ny, this%dx, this%dy)
end subroutine VelocitySolve

subroutine Timestep(this)
    class(grid_variables), INTENT(INOUT)   :: this
    !calculate advective timestep and min(adv,diff)
    !take advection and diffusion time step
    if ( (maxval(abs(this%ux(:,:))) <= 0.0001).and.(maxval(abs(this%uy(:,:))) <= 0.0001) ) then !check if velocity = 0
        this%dt = this%dt_dif
    else
        !Calculate dt_adv from vx, vy
        this%dt_adv = this%a_adv*min(this%dx/(maxval(abs(this%ux(:,:)))), this%dy/(maxval(abs(this%uy(:,:)))))
        this%dt = min(this%dt_dif, this%dt_adv)
    end if

    if ( this%dt > 0.6 * this%dx .OR. this%dt > 0.6 * this%dy) then
        this%dt = min(this%dx, this%dy) * 0.6
        print*, 'dt was corrected!'
    end if

    !set boundaries
    this%T(1,:)         = this%T(2,:)
    this%T(this%nx,:)   = this%T(this%nx-1,:)
    this%T(:,1)         = 1.
    this%T(:,this%ny)   = 0.

    this%W(:,1)         = 0. 
    this%W(:,this%ny)   = 0. 
    this%W(1,:)         = 0.
    this%W(this%nx,:)   = 0.
    
    !Calculate del2(T)
    this%d2Tdxdy = second_derivative(this%T,this%dx,this%dy)
    !Calculate del2(W)
    this%d2Wdxdy = second_derivative(this%W,this%dx,this%dy)
    !Calculate u.grad(T)
    this%ugradT = u_grad(this%ux, this%uy, this%T, this%dx, this%dy)
    !Calculate u.grad(W)
    this%ugradW = u_grad(this%ux, this%uy, this%W, this%dx, this%dy)
    !Calculate dTdx
    this%dTdx  = first_derivative_centred(this%T,this%dx, 'x')

    !update variables according to time step
    this%T(2:this%nx-1,2:this%ny-1) = this%T(2:this%nx-1,2:this%ny-1) + this%dt * &
        (this%d2Tdxdy(2:this%nx-1,2:this%ny-1) - this%ugradT(2:this%nx-1,2:this%ny-1))
    this%W(2:this%nx-1,2:this%ny-1) = this%W(2:this%nx-1,2:this%ny-1) + this%dt * &
        (this%Pr*this%d2Wdxdy(2:this%nx-1,2:this%ny-1)-this%ugradW(2:this%nx-1,2:this%ny-1) - &
        this%Ra*this%Pr*this%dTdx(2:this%nx-1,2:this%ny-1))
end subroutine Timestep

subroutine write(time_step, time_iteration, this)
    class(grid_variables), INTENT(inout)   :: this
    INTEGER, INTENT(IN)                    :: time_step
    real(real32),INTENT(IN)                :: time_iteration
    character(len=60)   :: outputfile_timestep_T    = ''
    character(len=60)   :: outputfile_timestep_S    = ''
    character(len=60)   :: outputfile_timestep_W    = ''
    character(len=60)   :: outputfile_uxmax_vs_time = ''
    character(len=60)   :: outputfile_uymax_vs_time = ''
    character(len=10)   :: Pr_number = ''

    !write Prandl number to file
    write(Pr_number,'(F8.3)') this%Pr
    call StripSpaces (Pr_number)

    !write time step to file
    write(Pr_number,'(F8.3)') time_step
    call StripSpaces (Pr_number)

    if ( time_step == 0 ) then
        outputfile_uxmax_vs_time = trim(Pr_number)//'_uxmax_vs_time.data'
        open(11,file=outputfile_uxmax_vs_time,action='write',position='append')
        outputfile_uymax_vs_time = trim(Pr_number)//'_uymax_vs_time.data'
        open(12,file=outputfile_uymax_vs_time,action='write',position='append')
    end if

    WRITE(11,*) maxval(this%ux)
    WRITE(12,*) maxval(this%uy)

    if((mod(time_step,100)==1) .or.( time_iteration>=this%total_time)) then
        outputfile_timestep_T =trim(Pr_number)//'_T_timestep_' //trim(Pr_number)//'.bin'
        outputfile_timestep_S =trim(Pr_number)//'_S_timestep_' //trim(Pr_number)//'.bin'
        outputfile_timestep_W =trim(Pr_number)//'_W_timestep_' //trim(Pr_number)//'.bin'
        print*, 'Calculating.... ', (time_iteration/this%total_time)*100, '%'
        open(1,file=outputfile_timestep_T,form='unformatted',access='stream')
        WRITE(1) this%nx, this%ny, this%T
        close(1)
        open(1,file=outputfile_timestep_S,form='unformatted',access='stream')
        WRITE(1) this%nx, this%ny, this%S
        close(1)
        open(1,file=outputfile_timestep_W,form='unformatted',access='stream')
        WRITE(1) this%nx, this%ny, this%W
        close(1)
    end if
    
    if(time_iteration>=this%total_time) then
        !write results for plot
        WRITE(11,*) maxval(this%ux)
        CLOSE(11)
        WRITE(12,*) maxval(this%uy)
        CLOSE(12)
    end if
    
end subroutine write

!generic functions used -> shuld be private

function second_derivative(func,hx,hy)  
    real(real64), INTENT(IN)                            :: hx, hy
    real(real64)                                        :: func(:,:)
    INTEGER                                             :: np_x, np_y
    real(real64), ALLOCATABLE                           :: func_second_deriv(:,:)
    real(real64)                                        :: second_derivative(size(func,1),size(func,2))
    INTEGER                                             :: i,j !local variable
    ALLOCATE(func_second_deriv(size(func,1),size(func,2)))
    np_x = size(func,1)
    np_y = size(func,2)

    func_second_deriv(1,:)      = 0
    func_second_deriv(np_x,:)   = 0
    func_second_deriv(:,1)      = 0
    func_second_deriv(:,np_y)   = 0

    do i = 2, np_x - 1
        do j = 2, np_y - 1
            func_second_deriv(i,j) =(func(i-1,j)+func(i+1,j) - 2*func(i,j))/hx**2 + &
            (func(i,j-1) + func(i,j+1) - 2*func(i,j))/hy**2
        end do
    end do

    !Assume derivatives at the end points 
    second_derivative = func_second_deriv  
    DEALLOCATE(func_second_deriv)
end function second_derivative

function first_derivative_centred(func,h,variable)
    character(len=1)                            :: variable
    real(real64), INTENT(IN)                    :: h
    real(real64)                                :: func(:,:)
    INTEGER                                     :: np_x, np_y
    real(real64), ALLOCATABLE                   :: func_first_deriv(:,:), first_derivative_centred(:,:) 
    INTEGER                                     :: i,j !local variable

    ALLOCATE(func_first_deriv(size(func,1),size(func,2)),first_derivative_centred(size(func,1),size(func,2)))

    np_x = size(func,1)
    np_y = size(func,2)

    func_first_deriv(1,:)      = 0
    func_first_deriv(np_x,:)   = 0
    func_first_deriv(:,1)      = 0
    func_first_deriv(:,np_y)   = 0
    if ( variable == 'x' ) then
        do i = 2, np_x - 1
            do j = 2, np_y - 1
                func_first_deriv(i,j) = (func(i+1,j) - func(i-1,j)) / (2*h)
            end do
        end do
    else if ( variable == 'y' ) then
        do i = 2, np_x - 1
            do j = 2, np_y - 1
                func_first_deriv(i,j) = (func(i,j+1) - func(i,j-1)) / (2*h)
            end do
        end do
    else
        print*,'Error!!! The variable is not correct!!'
        stop
    end if
    !Assume derivatives at the end points 
    first_derivative_centred = func_first_deriv
    DEALLOCATE(func_first_deriv)
end function first_derivative_centred


function u_grad(velocity_x, velocity_y, T_field, hx, hy)
    implicit none
    INTEGER                     :: i, j
    real(real64), INTENT(IN)    :: velocity_x(:,:), velocity_y(:,:), T_field(:,:), hx, hy
    real(real64), ALLOCATABLE   :: u_grad_cal_x(:,:), u_grad_cal_y(:,:), u_grad(:,:)

    ALLOCATE(u_grad_cal_x(size(T_field,1), size(T_field,2)))
    ALLOCATE(u_grad_cal_y(size(T_field,1), size(T_field,2)))
    ALLOCATE(u_grad(size(T_field,1), size(T_field,2)))

    if (  (maxval(abs(velocity_x(:,:))) <= 0.0001).and.(maxval(abs(velocity_y(:,:))) <= 0.0001) ) then
        u_grad = 0.
    else
        !boundaries
        u_grad_cal_x(:,size(T_field,2)) = 0
        u_grad_cal_x(size(T_field,1),:) = 0
        u_grad_cal_x(:,1) = 0
        u_grad_cal_x(1,:) = 0

        u_grad_cal_y(:,size(T_field,2)) = 0
        u_grad_cal_y(size(T_field,1),:) = 0
        u_grad_cal_y(:,1) = 0
        u_grad_cal_y(1,:) = 0

        do i = 2, size(T_field,1) - 1
            do j = 2, size(T_field,2) - 1
                !ux*grad(T)
                if ( velocity_x(i,j) > 0 ) then
                    u_grad_cal_x(i,j) = velocity_x(i,j) * ( T_field(i,j) - T_field(i-1,j) ) / hx
                else if ( velocity_x(i,j) < 0 ) then
                    u_grad_cal_x(i,j) = velocity_x(i,j) * ( T_field(i+1,j) - T_field(i,j) ) / hx
                end if
                !uy*grad(T)
                if ( velocity_y(i,j) > 0 ) then
                    u_grad_cal_y(i,j) = velocity_y(i,j) * ( T_field(i,j) - T_field(i,j-1) ) / hy
                else if ( velocity_y(i,j) < 0 ) then
                    u_grad_cal_y(i,j) = velocity_y(i,j) * ( T_field(i,j+1) - T_field(i,j) ) / hy
                end if
            end do
        end do
        u_grad = u_grad_cal_x + u_grad_cal_y
    end if
    DEALLOCATE(u_grad_cal_x,u_grad_cal_y)
end function u_grad

recursive function Vcycle_2DPoisson(u_f,rhs,h) result (resV)      
    !resV result to programm
    !u_f is the fine grid approximate solution
    !rhs = f, is the right hand side
    !h grid spacing
    implicit none
    real(real64)                 :: resV
    real(real64),intent(inout)   :: u_f(:,:)          ! arguments
    real(real64),intent(in)      :: rhs(:,:), h
    integer              :: nx,ny,nxc,nyc, i  ! local variables
    real(real64),allocatable     :: res_c(:,:),corr_c(:,:),res_f(:,:),corr_f(:,:)
    real(real64)                 :: alpha=0.7, res_rms

    nx=size(u_f,1); ny=size(u_f,2)  ! must be power of 2 plus 1
    if( nx-1/=2*((nx-1)/2) .or. ny-1/=2*((ny-1)/2) ) stop 'ERROR:not a power of 2'
    nxc=1+(nx-1)/2; nyc=1+(ny-1)/2  ! coarse grid size
    if (min(nx,ny)>5) then  ! not the coarsest level

        allocate(res_f(nx,ny),corr_f(nx,ny), &
        corr_c(nxc,nyc),res_c(nxc,nyc))

        !---------- take 2 iterations on the fine grid--------------
        res_rms = iteration_2DPoisson(u_f,rhs,h,alpha) 
        res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
  
        !---------- restrict the residue to the coarse grid --------
        call residue_2DPoisson(u_f,rhs,h,res_f) 
        call restrict(res_f,res_c)

        !---------- solve for the coarse grid correction -----------
        corr_c = 0.  
        res_rms = Vcycle_2DPoisson(corr_c,res_c,h*2) ! *RECURSIVE CALL*

        !---- prolongate (interpolate) the correction to the fine grid 
        call prolongate(corr_c,corr_f)

        !---------- correct the fine-grid solution -----------------
        u_f = u_f - corr_f  

        !---------- two more smoothing iterations on the fine grid---
        res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
        res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)

        deallocate(res_f,corr_f,res_c,corr_c)

    else  

        !----- coarsest level (ny=5): iterate to get 'exact' solution

        do i = 1,100
            res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
        end do

    end if

    resV = res_rms   ! returns the rms. residue

end function Vcycle_2DPoisson

function iteration_2DPoisson(u,f,h,a)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  Gauss-Seidel iteration                                            !
    !  is a function that returns the root mean square of the residual   !
    !  u is the approximated solution                                    !
    !    f is what is on the right hand side of the equal sine           !
    !  h is the grid spacing                                             !
    !  a is alpha which is the relaxation parameter                      !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    integer                                              :: i, j!local variable
    real(real64)                                         :: h, a, f(:,:)
    real(real64)                                         :: u(:,:)
    real(real64), allocatable                            :: residue (:,:)
    real(real64)                                         :: iteration_2DPoisson
    allocate(residue(size(u,1),size(u,2)))

    residue = 0.

    do concurrent (i = 2:size(u,1)-1, j=2:size(u,2)-1)
        residue(i,j)   = (u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - 4*u(i,j)) / h**2 - f(i,j)
        u(i,j)         = u(i,j) + a * (h**2)/4 * residue(i,j)
    end do
    iteration_2DPoisson = RMS(residue)

    DEALLOCATE(residue)
end function iteration_2DPoisson

subroutine residue_2DPoisson(u,f,h,res)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  calculates the residue in array res                               !
    !  u is the approximated solution                                    !
    !  f is what is on the right hand side of the equal sine             !
    !  h is the grid spacing                                             !
    !  res is residual array                                             !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer                                                    :: i, j           !local variable
    real(real64)                                               :: h, f(:,:)
    real(real64)                                               :: u(:,:)
    real(real64), INTENT(OUT), DIMENSION(size(u,1),size(u,2))  :: res
    !real(real64), INTENT(OUT),ALLOCATABLE                      :: res(:,:)           !output residual

    !ALLOCATE(res(size(u,1),size(u,2)))

    res = 0.
    do concurrent(i = 2:size(u,1)-1, j=2:size(u,2)-1)
        res(i,j) = (u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - 4*u(i,j)) / h**2 - f(i,j)
    end do
end subroutine residue_2DPoisson

subroutine restrict(fine, coarse)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  copies every other point from fine into coarse                    !
    !  fine is the fine grid                                             !
    !  coarse is the coarse grid                                         !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    integer                                                                   :: i_fine, j_fine !local variable
    real(real64), INTENT(IN)                                                  :: fine(:,:)      !input grid
    !real(real64), INTENT(OUT),ALLOCATABLE                                     :: coarse(:,:)    !output grid
    real(real64), INTENT(OUT)        :: coarse(size(fine,1)/2+1,size(fine,2)/2+1)         !output grid

    !ALLOCATE(coarse(size(fine,1)/2+1,size(fine,2)/2+1))

    do concurrent(i_fine=1:size(fine,1), j_fine=1:size(fine,2))
        if ( (mod(i_fine,2) > 0).and.(mod(j_fine,2) > 0)) then
            coarse(i_fine/2+1,j_fine/2+1) = fine(i_fine,j_fine)
        end if
    end do 
end subroutine restrict

subroutine prolongate(coarse, fine)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  copies coarse into every other point in fine                      !
    !  does linear interpolation to fill the other points                !
    !  fine is the fine grid                                             !
    !  coarse is the coarse grid                                         !
    !!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer                                                     :: i, j         !local variable
    real(real64), INTENT(OUT)                                   :: fine(:,:)    !input grid
    real(real64), DIMENSION(size(fine,1)/2+1,size(fine,2)/2+1)  :: coarse       !output grid
    !real(real64), ALLOCATABLE  :: coarse(:,:)       !output grid
    !ALLOCATE(coarse(size(fine,1)/2+1,size(fine,2)/2+1))

    do concurrent(i=1:size(fine,1), j=1:size(fine,2))   
        if ( (mod(i,2) > 0).and.(mod(j,2) > 0)) then
            fine(i,j) = coarse(i/2+1,j/2+1) 
        else if ((mod(i,2) > 0).and.(mod(j,2) == 0)) then 
            fine(i,j) = 0.5*(coarse(i/2+1,j/2) + coarse(i/2+1,j/2+1))
        else if ((mod(i,2) == 0).and.(mod(j,2) > 0)) then
            fine(i,j) = 0.5*(coarse(i/2,j/2+1) + coarse(i/2+1,j/2+1))
        end if
    end do
    do concurrent(i=1:size(fine,1), j=1:size(fine,2))
        if ((mod(i,2) == 0).and.(mod(j,2) == 0)) then
            fine(i,j) = 0.25*(fine(i-1,j)+fine(i+1,j)+fine(i,j-1)+fine(i,j+1))
        end if
    end do
    !DEALLOCATE(coarse)
end subroutine prolongate

real(real64) function RMS(A)
    implicit none
    real(real64), INTENT(IN) :: A(:,:)
    INTEGER                  :: i,j
    real(real64)             :: sum_squared!, RMS
    sum_squared = 0.0

    do concurrent (i = 1:size(A,1), j = 1:size(A,2))
        sum_squared = sum_squared + A(i,j)**2
    end do
    RMS = sqrt(sum_squared/(size(A,1)*size(A,2)))
end function RMS

function solve_Poisson_Eq_multigrid(u,f,h, crit_given)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  Solve Poisson Equation using multigrid                            !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer                                              :: iteration, iteration_max
    real(real64)                                         :: h, crit_given, crit_current, f(:,:)
    real(real64)                                         :: u(:,:), residual_func
    real(real64), allocatable                            :: solve_Poisson_Eq_multigrid(:,:)
    allocate(solve_Poisson_Eq_multigrid(size(u,1),size(u,2)))
    !set values back
    residual_func = 0.   !changed frim residual to residual_func
    iteration_max = 100
    iteration = 0
    crit_current = 1.
    if ( (RMS(f) <= 0.0001).and.(RMS(f) >= - 0.0001) ) then   !check if rms(f)==0
        solve_Poisson_Eq_multigrid = 0.
    else 
        do while ((crit_current >= crit_given).and.(iteration < iteration_max) )
            residual_func = Vcycle_2DPoisson(u,f,h)
            iteration = iteration + 1
            crit_current = residual_func/rms(f)
        end do
        solve_Poisson_Eq_multigrid = u
    end if
end function solve_Poisson_Eq_multigrid

subroutine S_to_vel(S,ux, uy, nx, ny, dx, dy)
    implicit none
    real(real64), allocatable,intent(inout)  :: S(:,:)
    integer, INTENT(IN)                   :: nx, ny
    real(real64),intent(in)               :: dx, dy
    real(real64), allocatable,intent(out) :: ux(:,:), uy(:,:)
    S(1,:)  = 0.
    S(nx,:) = 0.
    S(:,1)  = 0.
    S(:,ny) = 0.
    !Calculate ux, uy from S
    ux =  first_derivative_centred(S(:,:),dy,'y')
    uy =  -1*first_derivative_centred(S(:,:),dx,'x')
end subroutine S_to_vel

subroutine StripSpaces(string)
    character(len=*) :: string
    integer :: stringLen 
    integer :: last, actual
    stringLen = len (string)
    last = 1
    actual = 1
    do while (actual < stringLen)
        if ( string(last:last) == '.' ) then
            string(last:last) = '_'
        end if
        if (string(last:last) == ' ') then
            actual = actual + 1
            string(last:last) = string(actual:actual)
            string(actual:actual) = ' '
        else
            last = last + 1
            if (actual < last) &
                actual = last
        endif
    end do
end subroutine
end module module_LowPrandl

and this input file:

 &inputs
  Pr=0.01   ! <---- CHANGE ONLY THIS 
  nx=257    ny=65
  total_time=0.1
  Ra=1.e6
  err=1.e-3
  a_dif=0.15   a_adv=0.4
  Tinit='cosine'
  /
gfortran -o ex1 <some_options> modul_Low_Prandl.f90 exercise_11.f90 lowPrandt_parameters_FULL.txt
./a.out

You don’t compile a data file such as lowPrandt_parameters_FULL.txt. You open it within your program, as demonstrated here.

Regarding another question you asked, note that procedures can be declared private. What I do is

module m
implicit none
private ! makes private the default
public :: foo,bar ! only declare public the entities that invoked from outside the module
contains
...
end module m
1 Like

thank you very much I can now compile without a problem. Regarding the second question, if I understand it correctly the functions second_derivative and first_derivative_centred are therefore already private? Or do I need to implement it in the following manner?

    use iso_fortran_env
    implicit none
    private
    public :: grid_variables
    real(real64)        :: Pi =  atan(1.)*4.

    type :: grid_variables
        !input parameters
        real(real64)        :: total_time = 0.1                     !total integration time
        character(len=10)   :: Tinit='cosine'                       !Temperature initial condition
        integer             :: nx, ny                               !grid spacing
        real(real64)        :: a_dif, a_adv                         !time step constant
        real(real64)        :: err = 1e-6                           !maximum error in Poisson solves = resiude_rms/rhs_rms<err
        real(real64)        :: Ra                                   !Rayleigh number
        real(real64)        :: Pr                                   !Prandt number
        !character(len=100)  :: fname = 'parameter.txt'              !read command-line arguments

        !Variables in type
        real(real64), allocatable   :: T(:,:),dTdx(:,:), d2Tdxdy(:,:),d2Wdxdy(:,:)
        real(real64), allocatable   :: ux(:,:), uy(:,:)
        real(real64), allocatable   :: S(:,:), ugradT(:,:), ugradW(:,:), W(:,:)
        real(real64)                :: dx, dy, h, dt = 0.1, dt_dif = 0.1, dt_adv = 0.1

      contains
        procedure,public,pass(this) :: read_inputs
        procedure,public,pass(this) :: initialise
        procedure,public,pass(this) :: VelocitySolve
        procedure,public,pass(this) :: Timestep
        procedure,public,pass(this) :: write

        procedure,private,pass(hx) :: second_derivative
        procedure,private,pass(h)  :: first_derivative_centred

  end type grid_variables

if I do ths however I get the following error message:

modul_Low_Prandl.f90:33:17:

   33 |         procedure,private,pass(h)  :: first_derivative_centred
      |                 1
Error: Non-polymorphic passed-object dummy argument of ‘first_derivative_centred’ at (1)
modul_Low_Prandl.f90:32:17:

   32 |         procedure,private,pass(hx) :: second_derivative
      |                 1
Error: Non-polymorphic passed-object dummy argument of ‘second_derivative’ at (1)

How could I avoid that? Or do they appear because the functions do not use the derived type?