0

VERY new to fortran and feeling lost. What I am trying to do is calculate pressure at water depth based on initial atmospheric pressure. Then, each pressure will be based on the one previous to it. Depth will increase by (-5) at each step. The compiler does not like the notation of z(i)=(0:5:-50), which I intend to mean "jumps by -5 from 0 to -50.

This process needs to be repeated for 5 different initial atmospheric pressures. What is the best way to then organize this info into columns to be printed?

Does this look like a reasonable way to go about accomplishing this goal?

Many thanks!

       program pressure
       real pk(11), zp(11), g, rz,z,dz, p_atm
       g= 9.81
       z(i)=(0:5:-50)
       dz= -5
       p_atm= (/0.97, 1.02, 1.04, 1.03, 1.01/)
       
       do 
         i= 1,11
         d= z(i)
         
         if(d.GT.-10.0) then
           rz= 1020.0
         elseif(d.LT.-10.0 .and. d.GT.-50.0) then 
           rz= 1020.0+(0.25*(ABS(d+10.0)))
     
         if d == 0
           p(i)= -(((-rz*g)*dz))+p_atm)
         else
           p(i)= -(((-rz*g)*dz))+p(i-1)
         
         endif

       open(20,file='pressure.txt', status='unknown')
       rewind(20)
       write(20,25) 'Model Day', '1', '2', '3', '4', '5'
       write(20,26) 'Patm (10^5 Pa)', '0.97', '1.02', '1.04', '1.03', '1.01'
  26   format(a14,12F9.2)
  25   format(a10,12F9.0)
  • Calculus has been around for a while now. Rather than computationally approximate integration, why not just express the integral as code in a single calculation? – Bohemian Feb 11 '22 at 20:37
  • 1
    You should not just leave your old question open and create a new one, you should have used [edit] to make the old one better. I vote to close the old one as a duplicate of this one, but it is better to delete it. – Vladimir F Героям слава Feb 11 '22 at 20:40
  • 2
    You should at least make the code you paste here syntactically correct. I assume the `i= 1,11` belongs to the `do` on the previous line. But it must be on the some line. You are missing the `end do`. – Vladimir F Героям слава Feb 11 '22 at 20:42
  • As for the logic itself, in your if condition, you should rather test the value of i, not the value of d. Comparison of floating point numbers is often inexact. Even better, just declare `p` to start at 0 and not at 1 and declare the value of `p(0)` somehow. You might be interested how I do a very similar kind of pressure integration at https://bitbucket.org/LadaF/elmm/src/313324325cde1cc7a6f4ea52642018abd29ea557/src/pressure.f90#lines-867 and later https://bitbucket.org/LadaF/elmm/src/313324325cde1cc7a6f4ea52642018abd29ea557/src/pressure.f90#lines-903 – Vladimir F Героям слава Feb 11 '22 at 20:51

1 Answers1

0

Okay, let's go from the top.

  1. Always use implicit none. If you don't use implicit none, you will get bugs.
  2. Your indenting looks a bit like fixed-form fortran, but you're actually using free-form fortran (and we are in 2022, everyone should be using free-form fortran), so your indentation should really be consistent.
  3. You don't have an end program line or an enddo line, presumably because this is just a snippet of code.
  4. As VladimirF points out, the line do i=1,11 must be one line, not two.

Lets fix these first, to give

program pressure
  implicit none
  
  real pk(11), zp(11), g, rz,z,dz, p_atm
  
  g= 9.81
  z(i)=(0:5:-50)
  dz= -5
  p_atm= (/0.97, 1.02, 1.04, 1.03, 1.01/)
  
  do i= 1,11
    d= z(i)
    
    if(d.GT.-10.0) then
      rz= 1020.0
    elseif(d.LT.-10.0 .and. d.GT.-50.0) then 
      rz= 1020.0+(0.25*(ABS(d+10.0)))
     
    if d == 0
      p(i)= -(((-rz*g)*dz))+p_atm)
    else
      p(i)= -(((-rz*g)*dz))+p(i-1)
    endif

    open(20,file='pressure.txt', status='unknown')
    rewind(20)
    write(20,25) 'Model Day', '1', '2', '3', '4', '5'
    write(20,26) 'Patm (10^5 Pa)', '0.97', '1.02', '1.04', '1.03', '1.01'
    26   format(a14,12F9.2)
    25   format(a10,12F9.0)
  enddo
end program

Okay, now we can try compiling and see what happens. The compiler is going to shout at us for a number of reasons:

  1. The variables i, d and p have not been declared.
  2. z(i) either means a function call or an array element access (from context, you want the latter), but z is defined as a scalar variable, not a function or an array.
  3. p_atm is being initialised as an array, but was declared as a scalar. You're also trying to add p_atm as a scalar, so let's add a loop over these pressures.
  4. The if and endif statements don't match; you need to endif your first if block.
  5. Your format strings don't match: you're printing character strings, but formatting them as floats ('0.97' is a character string, 0.97 is a float).
  6. Your second if statement needs to have appropriate brackets and a then.
  7. You close more brackets than you open in your first p(i)=... line. I'm going to have to guess here.
  8. The comparison d==0 is a really bad idea. Never compare floats for equality. Let's replace this with i==1.

Okay, let's fix these problems:

program pressure
  implicit none
  
  real pk(11), zp(11), g, rz, z(11), dz, p_atm(5), p(11), d
  integer i,j
  
  g= 9.81
  z(i)=(0:5:-50)
  dz= -5
  p_atm= (/0.97, 1.02, 1.04, 1.03, 1.01/)
  
  do j=1,5
    do i= 1,11
      d= z(i)
      
      if(d.GT.-10.0) then
        rz= 1020.0
      elseif(d.LT.-10.0 .and. d.GT.-50.0) then 
        rz= 1020.0+(0.25*(ABS(d+10.0)))
      endif
     
      if (i == 1) then
        p(i)= -(((-rz*g)*dz))+p_atm(j)
      else
        p(i)= -(((-rz*g)*dz))+p(i-1)
      endif
      
      open(20,file='pressure.txt', status='unknown')
      rewind(20)
      write(20,'(a14,5a9)') 'Model Day', '1', '2', '3', '4', '5'
      write(20,'(a14,5a9)') 'Patm (10^5 Pa)', '0.97', '1.02', '1.04', '1.03', '1.01'
    enddo
  enddo
end program

Okay, now we get to initialising z. The most straightforward way of doing this is with a simple loop:

do i=1,11
  z(i) = 5-5*i
enddo

but if you want to get fancy with it, you can use an implicit do loop:

z = [(5-5*i, i=1, 11)]

Also, given your output file header, it probably makes sense to store p for all values of i and j, so let's make it a 2D array. Let's also move the file-writing to after we've calculated all values of p. Then, making a few modernisations along the way, your program looks like:

program pressure
  implicit none
  
  real :: pk(11), zp(11), g, rz, z(11), dz, p_atm(5), p(11,5), d
  integer :: i,j
  
  g = 9.81
  z = [(5-5*i, i=1, 11)]
  dz = -5
  p_atm = [0.97, 1.02, 1.04, 1.03, 1.01]
  
  do j=1,5
    do i=1,11
      d = z(i)
      
      if (d > -10.0) then
        rz = 1020.0
      elseif (d < -10.0 .and. d > -50.0) then 
        rz = 1020.0 + 0.25*abs(d+10.0)
      endif
     
      if (i == 1) then
        p(i,j) = -(-rz*g*dz)+p_atm(j)
      else
        p(i,j) = -(-rz*g*dz)+p(i-1,j)
      endif
    enddo
  enddo
  
  open(20, file='pressure.txt', status='unknown')
  write(20,'(a14,5a9)') 'Model Day', '1', '2', '3', '4', '5'
  write(20,'(a14,5a9)') 'Patm (10^5 Pa)', '0.97', '1.02', '1.04', '1.03', '1.01'
  do i=1,11
    write(20,'(5f9.2)') p(i,:)
  enddo
  close(20)
end program
veryreverie
  • 2,871
  • 2
  • 13
  • 26
  • All I might suggest adding to this comprehensive (and very generous) answer would be some clarification for OP along the lines *Your question title suggests a confusion between the concepts of recursion and iteration (the fancy word for 'looping'). In this case iteration is the more natural approach. Lessons in recursion can wait.* – High Performance Mark Feb 12 '22 at 11:17
  • @HighPerformanceMark I think anything I write on that would broadly be "see HPM's comment". It also feels like a bit of a tangent from the rest of my own answer. I wonder if you might like to add your own answer covering the point? – veryreverie Feb 12 '22 at 12:53