How does Fortran know the Inputs and Outputs of Subroutines

[ecastro]

I have trouble understanding how Fortran knows the inputs and outputs of a subroutine. I am studying a code written in Fortran, but I do not have a compiler, so I cannot check how it works.

Here is my problem; for example, I have a subroutine,

subroutine test(x, y, z, a, b, c)

The variables x, y, and z have values in the main program and are needed in the subroutine, while a, b, and c are the supposed outputs. What troubles me are:
1. Should all six variables be defined in the main program and in the subroutine (defined to be real, integer, complex, etc.)?
2. Is it alright if I define (defined to be real, integer, complex, etc.) a variable in the subroutine and not in the main program? How about vice versa? What are the consequences of this action?
3. If x, y, and z somehow receive a new value in the subroutine, will the previous values be replaced in the main program (x, y, and z become part of the output)?
4. If there's a variable named d in the subroutine, will I be able to call it in the main program even though it's not part of the output?

Thank you in advance.

 

[BvU]

1. Should all six variables be defined in the main program and in the subroutine (defined to be real, integer, complex, etc.)?

Yes, but they can have different names - the names are symbolic and serve as references
The fortran parameter passing mechanism is mainly by reference: the address of a variable is passed to the subroutine.

2. Is it alright if I define (defined to be real, integer, complex, etc.) a variable in the subroutine and not in the main program? How about vice versa? What are the consequences of this action?

yes. no consequences. In some compilers you can not be sure a variable in a subroutine is static (has the same value in a call as it had been given in a preceding call)

3. If x, y, and z somehow receive a new value in the subroutine, will the previous values be replaced in the main program (x, y, and z become part of the output)?

A good environment will protect you against such an access violation, but I'm not 100% sure.
[edit] Ignore. Confusion with newer languages on my part. Vanadium corrects me in post #2 and Anorlunda elucidates in #3.
F90 has some features in this direction

4. If there's a variable named d in the subroutine, will I be able to call it in the main program even though it's not part of the output?

No. Other routines (including the main program) do not know where to find d

 

[Vanadium 50]

3. If x, y, and z somehow receive a new value in the subroutine, will the previous values be replaced in the main program (x, y, and z become part of the output)?

Yes. This is a feature of the language.

 

[anorlunda]

In older FORTRANs, all arguments are by-reference and inherently input and/or output. That is very simple, but it lacks safeguards against errors.

Consider this routine

SUBROUTINE BAZFAX(X,Y,Z)
INTEGER X,Y,Z
Z=0
X=X/Y
RETURN
END

It appears that:
X is input and output
Y is input
Z is output
It is also possible to use IF statements to make it conditional whether an argument is used as input, output or both.

Now think what might happen with these calls

CALL BAZFAZ(1,2,3)
CALL BAZFAZ(X*2,Y*3,SQRT(Z))
CALL BAZFAX(X,X,X)

The first two calls have no arguments suitable for outputs.
The last call will crash with divide by zero because in the SUBROUTINE, Y and Z refer to the same memory address.

The point is that the caller must know how the subroutine is written and which arguments are used as outputs, inputs, or both. You can get into all sorts of trouble if you do it wrong. That is why other languages include safeguards to prohibit some of these kinds of errors.

EDIT: I neglected to mention that you can also declare the types of arguments to be different than the types supplied in the CALL. Even incompatible types. You can also have an inconsistent count of the number of arguments between the caller and the routine.

All these things are consequences of a very simple model. Arguments are merely addresses of a memory region to use for whatever you want. FORTRAN compilers compile programs and subprograms separately, without any prior knowledge of the other, nor any macro-data left behind about the signature. That was the optimum strategy when computer resources were severely limited, and when compiling programs was nearly as difficult as executing them.

 

[Mark44]

Newer versions of Fortran (Fortran 90?) use syntax in the declaration of subroutine and function parameters that indicates whether a parameter should be considered input, output, or both.

subroutine square_cube(i, isquare, icube)
  integer, intent(in)  :: i             ! input
  integer, intent(out) :: isquare, icube ! output
  isquare = i**2
  icube   = i**3
end subroutine square_cube

program xx
  implicit none
  integer :: i, isq, icub
  i = 4
  call square_cube(i, isq, icub)
  print*,"i,i^2,i^3=",i,isq,icub
end program xx

 

[BvU]

I am studying a code written in Fortran

Ha, I missed that in first reading. What's the goal: reverse engineering for technology, porting to a newer environment, ?
[edit] and your source code dates from when, approximately ? Any knowlegde about the environment it was used in (IBM, DEC, CDC... ) ?

 

[ecastro]

Thank you for your comments and precautions on the subroutines. I have another question, can a subroutine access a variable that is not included as an input but defined in the main program?

Ha, I missed that in first reading. What's the goal: reverse engineering for technology, porting to a newer environment, ?
[edit] and your source code dates from when, approximately ? Any knowlegde about the environment it was used in (IBM, DEC, CDC... ) ?

I think more like porting to a newer environment, a checking if the code I am writing is correct. The code was last updated May 2015, but I have no knowledge of its environment.

 

[BvU]

2015 is very recent for a Fortran IV programmer A subroutine does not have access to the adresses (variables) the main program uses (unless they are passed to the subroutine as arguments).

However:
The issue I'm about to adress now requires distinction between Fortran 90 and predecessors like Fortran 77.
One of the gizmos in older Fortran that has both been very useful AND a horrifiying source of errors is the COMMON block. If that doesn't occur in your code, then you're lucky; if it does, then that's a way for the routines (in which such a common block is declared) to all access the same variables.

In fortran 90 that role is taken over (in a more structured way) by modules

porting to a newer environment

I hope you know what you are doing. Porting is error-prone and you may well end up with a program that is less efficient. You sure the investment in a compiler and learning some mixed-language procedure calling isn't much more effective and time-saving ?

 

[Mark44]

I have another question, can a subroutine access a variable that is not included as an input but defined in the main program?

As BvU notes, this is possible to do using a COMMON block. However, it's considered very bad practice, in any programming language, for a function or subroutine to access variables other than its own local variables or those passed as parameters.

 

[BvU]

Dear Mark, OP has indicated he has to deal with whatever is in there already !

Interesting case, isn't it ?

 

[jbriggs444]

Thank you for your comments and precautions on the subroutines. I have another question, can a subroutine access a variable that is not included as an input but defined in the main program?

If one is a tricky sort with either faith or inside knowledge about the storage layout in the main program, one could declare a parameter in a subroutine as an array and proceed to use array references to access variables not appearing in the parameter list.

e.g.

       integer i, j, k
       call sub ( j )
       end
       [...]
       subroutine sub ( m )
       integer m(3)
       m(0) = 1        ! Modifies i (and may require you to have compiled the subroutine with array bounds checking disabled)
       m(1) = 2        ! Modifies j
       m(2) = 3        ! Modifies k
       end

This is, of course, very bad practice. In addition to violating the expectations of all sane code readers, an optimizing compiler could render the approach useless, making memory allocation order unpredictable, moving variables into registers or optimizing them away altogether. One of the virtues (and one of the problems) associated with using common blocks is that memory allocation order is constrained and the compiler is not permitted to make such optimizations.

With DEC Fortran-77 (and presumably others), one could make games like this somewhat more kosher by using directives like %loc(variablename) to obtain the address of a variable. The resulting address could then be passed to a subroutine. The compiler was smart enough to use the appearance of %loc(x) as a hint that x was volatile and subject to modification behind the compiler's back -- so that optimizations for that variable would not be attempted.

 

[ecastro]

Unfortunately, there are multiple common blocks in the subroutine. In my case, there is this:

common /common_name/ x, y, z

This is in both the main program and in the subroutine. If I understand correctly on how these work, you can access the variables if 'common_name' was called inside the subroutine. Does this mean if the variables 'x', 'y', and 'z' changed values inside the subroutine, they will look like outputs? That is, if these variables were called after the subroutine, they will have new values, right?

I hope you know what you are doing. Porting is error-prone and you may well end up with a program that is less efficient. You sure the investment in a compiler and learning some mixed-language procedure calling isn't much more effective and time-saving ?

I am almost done with the code, the problem is my output from my code is not what I expect.

@jbriggs444 I don't understand your example. There is a subroutine named 'sub' which has a single input 'm' (from the main program, the input is 'j'). How does it modify 'i' and 'k'? Isn't the output of your subroutine just modifies the input?

 

[Mark44]

I am almost done with the code, the problem is my output from my code is not what I expect.

I hope you know what you are doing. Porting is error-prone and you may well end up with a program that is less efficient.

Or one that produces output that you don't expect. As BvU said, porting is error-prone...

@jbriggs444 I don't understand your example. There is a subroutine named 'sub' which has a single input 'm' (from the main program, the input is 'j'). How does it modify 'i' and 'k'? Isn't the output of your subroutine just modifies the input?

Fortran is a call-by-reference language. This means that any subroutine or function can modify its parameters. In contrast, C and other languages that stem from C (such as C++ and others) are call-by-value languages, whereby a function cannot modify its parameters.

In call-by-reference what is actually passed in a function or subroutine call is the address of the parameter. "call sub(j)" is passing the address of j. In the sub subroutine, the parameter is declared as a three element array. Fortran arrays have indexes that normally start at 1, so the value for the parameter j would normally be at m(1). By accessing m(0), the subroutine has access to i (in the main program), and m(1) provides access to j (in the main program), and m(2) provides access to k (in the main program). The subroutine doesn't know about the identifiers i, j, and k, but it has access to them via the array m.

 

[ecastro]

In call-by-reference what is actually passed in a function or subroutine call is the address of the parameter. "call sub(j)" is passing the address of j. In the sub subroutine, the parameter is declared as a three element array. Fortran arrays have indexes that normally start at 1, so the value for the parameter j would normally be at m(1). By accessing m(0), the subroutine has access to i (in the main program), and m(1) provides access to j (in the main program), and m(2) provides access to k (in the main program). The subroutine doesn't know about the identifiers i, j, and k, but it has access to them via the array m.

This is problematic... So, it doesn't just change the input 'j' (presumably a single-valued integer value) to an array?

 

[Mark44]

This is problematic... So, it doesn't just change the input 'j' (presumably a single-valued integer value) to an array?

No. The main program (the caller of sub()) knows about i, j, and k. It (main) doesn't know about m. The subroutine, sub(), knows about m, but doesn't know about i, j, or k. By passing addresses, which is how it works in Fortran and other call-by-reference languages, a subroutine doesn't have to know the names of variables, but can still have read or write access to them.

The caller (main program) sends the address of j as the argument to sub() in the call. The subroutine treats the address that is sent as an array with three elements: m(1), m(2), and m(3).

The value stored in the address that m(1) represernts is the value stored in j. Likewise, the value stored in the address that m(2) represents is the value stored in k.m(0) represents the memory address that comes before (i.e., lower in memory) has the value stored in i. This is my understanding of jbriggs888's example, which I'm pretty sure is what he was driving at.

 

[BvU]

Unfortunately, there are multiple common blocks in the subroutine. In my case, there is this:

common /common_name/ x, y, z

This is in both the main program and in the subroutine. If I understand correctly on how these work, you can access the variables if 'common_name' was called inside the subroutine. Does this mean if the variables 'x', 'y', and 'z' changed values inside the subroutine, they will look like outputs? That is, if these variables were called after the subroutine, they will have new values, right?

I see. These variable addresses are known in all routines that have the common declared. Note I use the word address -- the variable names may differ from routine to routine. Not good practice, but it happens.

Your understanding is correct. Such variables 'look like' outputs and also like inputs. They are common (their 'scope' spans the routines that declare the common block). You have to be careful to check which routine is responsible for e,g, initializing them.

I am almost done with the code, the problem is my output from my code is not what I expect.

Painted yourself in a corner. If it's a big program you could spend a long time debugging. If you are unlucky, not just your code but the original code as well. Do you have output from the old program to compare with ? What's the language you are porting to ? Can you restrict your debugging to a limited set of source code and give concrete examples ? Is it all very confidential ? Did you use automatic code converting ?
(there's a lot of very, very experienced fortranners about who'd love to help you by digging into this out of pure sentimentality )

@jbriggs444 I don't understand your example. There is a subroutine named 'sub' which has a single input 'm' (from the main program, the input is 'j'). How does it modify 'i' and 'k'? Isn't the output of your subroutine just modifies the input?

JB example is malignant. Be optimistic (you must be) and give the original programmers the benefit of the doubt.

 

[FactChecker]

Unfortunately, there are multiple common blocks in the subroutine. In my case, there is this:

common /common_name/ x, y, z

This is in both the main program and in the subroutine. If I understand correctly on how these work, you can access the variables if 'common_name' was called inside the subroutine. Does this mean if the variables 'x', 'y', and 'z' changed values inside the subroutine, they will look like outputs? That is, if these variables were called after the subroutine, they will have new values, right?

Yes.

I am almost done with the code, the problem is my output from my code is not what I expect.

@jbriggs444 I don't understand your example. There is a subroutine named 'sub' which has a single input 'm' (from the main program, the input is 'j'). How does it modify 'i' and 'k'? Isn't the output of your subroutine just modifies the input?

It thinks that m is an array and it knows the starting address of the array. So it can modify values from that address and higher when the index of m increases. (The case of m(0) => i is special and only if certain compiler options have been set. Usually the lowest index of a FORTRAN array is 1.)

 

[jbriggs444]

It thinks that m is an array and it knows the starting address of the array. So it can modify values from that address and higher when the index of m increases. (The case of m(0) => i is special and only if certain compiler options have been set. Usually the lowest index of a FORTRAN array is 1.)

Yes.

For the compiler that I used to use in the 1980's (VAX Fortran on VMS) one would have needed to compile with options like:

    $ fortran test /check=nobounds
    $ link test
    $ run test

By default, the compiler would build in run time array bounds checking so that an attempt to access array element zero (one cell before the array started at element one) would throw an error at run time. I no longer remember if it was smart enough to throw the error at compile time when the index was a knowable compile time constant.

While one could declare an array as, for instance,

       subroutine sub ( m )
       int m(0:2)
       m(0) = 1
       m(1) = 2
       m(2) = 3

this would not accomplish the goal of accessing a memory location prior to the start of the array. The parameter passed in the argument list is interpreted as the address of the first array element based on the way the array is declared in the subroutine. As declared above, the first element of m is at index zero. The memory location one cell prior would be m(-1) and would still trigger an array bounds error, by default, if access were attempted.

 

[anorlunda]

Unfortunately, there are multiple common blocks in the subroutine. In my case, there is this:

common /common_name/ x, y, z

It is worse than you imagine.

1. First, FORTRAN also had a feature called EQUIVALENCE. That allowed the programmer to create an indefinite number of aliases for the same memory locations, either COMMON memory or stack memory.
2. Second, it was regular practice to use COMMON in ways to mimic the use of "heap" storage in other languages. In a language like C, one could use x=malloc(1000) to reserve 1000 bytes of memory, and free(x) to release the memory block. In FORTRAN we might say COMON/BAZFAZ/I,J,K in one subprogram, COMMON/BAZFAZ/X(3) in another subprogram, and COMMON/BAZFAZ/A,B(2) in a third subprogram. In all three cases, the same memory block is referenced. That allowed memory to be reused for different purposes in different phases of the program, just like heap storage in other languages.
3. Third, in the days when FORTRAN had no INCLUDE feature, and before the time when text editors had copy/paste, or when we could store source code in files, the COMMON declaration in each subprogram had to be retyped by hand. A single character typo, could be enough to misalign the COMMON declarations leading to bugs that were extraordinarily hard to find. The typos did not cause compilation errors. More than once I helped hapless programmers proof read their COMMON declarations with great tedium, character-by-character subprogram-by-subprogram. With punch cards, you could reliably duplicate cards one-at-a-time but making 200 duplicates of a 100 card long COMMON declaration was its own kind of tedium and still gave lots of opportunities for errors.
4. There was one context where "global COMMON" was not evil. I worked on real time simulators where programs executed 10 to 50 times per second, and could not be slowed down or halted temporarily for debugging purposes. We made 100% of the variables global. Even local temps were made global. That made it possible to make real time debuggers that could capture data for debug purposes from any part of any program without changing the real time performance. It worked remarkably well, streamlining debugging and interprogram communications. Imagine just pointing a D-A converter to a memory address, and then to an oscilloscope or a chart recorder without changing any program code. The evil of global data comes when programmers forget which data is local and which is global, but if 100% is global, there is no confusion.

The OP faced with the task of maintaining a very old program is in a bad spot. Because of call by ref, COMMON and EQUIVALENCE, you can not dependably break down the task into program modules. You may need to examine every line of code in every subprogram for possible hidden interactions with other subprograms.

The SUBROUTINE and FUNCTION features in FORTRAN provided code reuse and separate compilation, but they did not provide modularity.

Thinking back on the constraints we lived with, it sounds impossible that we could achieve what we did, but we did make things work correctly.

 

[FactChecker]

Thinking back on the constraints we lived with, it sounds impossible that we could achieve what we did, but we did make things work correctly.

That being said, there were many times when common blocks were the only realistic way to get things done.
One further complication was that the same position in a common block could be used under different names in different subroutines. And a subroutine that needed only the 58'th and 68'th variables might get them by: common /xxx/ dummy1(57), var1, dummy2(9), var2. (Or something like that. My memory is vague on it.)
I had to write some programs to track down and cross-reference all the usage of common block data elements in large programs.

 

[ecastro]

Painted yourself in a corner. If it's a big program you could spend a long time debugging. If you are unlucky, not just your code but the original code as well. Do you have output from the old program to compare with ? What's the language you are porting to ? Can you restrict your debugging to a limited set of source code and give concrete examples ? Is it all very confidential ? Did you use automatic code converting ?

It is indeed a long program. The code is from NASA, so I trust that it works. The code is not that confidential, it is actually available here: http://6s.ltdri.org/pages/downloads.html, and there is an example output in the user manual. I am learning the physics of the code, and I guess if I were able to replicate the code, I understand it well.

I am actually looking every line of code, and I haven't seen the references of @anorlunda and @FactChecker. As far as I have seen, every time the commons are declared, they have the same set of names and every variable is present, as to how they were declared in the beginning of the main program.

I came asking here because there is a variable, 'roatm', (an array to be specific) that was only defined in the main program, and by following the inputs of the code manual, the array has no values but it is used to calculate a quantity, which returns a not-a-number value.

I am also curious if I understand this method of Fortran correctly,

do 10 i = a, b
if (insert condition) goto 10
[other processes]
10 continue

If the condition is satisfied, the process should jump to the line 10. Then, does this continue to the next iteration, or does it break the iteration?

 

[FactChecker]

I am also curious if I understand this method of Fortran correctly,

do 10 i = a, b
if (insert condition) goto 10
[other processes]
10 continue

If the condition is satisfied, the process should jump to the line 10. Then, does this continue to the next iteration, or does it break the iteration?

Yes, it continues to the next iteration. Line 10 is considered to be still inside the "do loop". To jump out of the loop it has to goto a line other than 10, outside of the loop

 

[BvU]

so I trust that it works

Dangerous ! But it looks decent, is pretty modern and doesn't look devious (no EQUIVALENCE statements).

learning the physics of the code, and I guess if I were able to replicate the code, I understand it well

That's rather roundabout (for the physics especially)

I came asking here because there is a variable, 'roatm', (an array to be specific) that was only defined in the main program, and by following the inputs of the code manual, the array has no values but it is used to calculate a quantity, which returns a not-a-number value.

roatm(3,20) sits is in common /sixs_disc/ which is declared in DISCOM.f and in INTERP.f, as well as in main.f, mainlutareo.f and SPECINTERP.f

Values are written into it in DISCOM

        roatm(1,l)=rorayl
        roatm(2,l)=romix
        roatm(3,l)=roaero

In other words: a more specific question in post #1 would have given you this answer a long time ago ! But you have learned a lot about Fortran in the mean time, so it's not a waste.

 

[ecastro]

roatm(3,20) sits is in common /sixs_disc/ which is declared in DISCOM.f and in INTERP.f, as well as in main.f, mainlutareo.f and SPECINTERP.f

Yes, the values of 'roatm' come from the DISCOM subroutine. However, the whole subroutine is under four conditions, two of the first conditions are for entering the subroutine while the last two skips it. It would seem that I misunderstood the process of the do-go-to loop. I thought that if the code were to have a go-to under a do-loop, the whole do-loop breaks, so I had thought that I missed something on the processes of subroutines.

Thank you for all your help, I hope this works now.

 

[FactChecker]

Did the original NASA code generate NaNs, or was that your code?

 

[Vanadium 50]

The code is from NASA, so I trust that it works.

Three words. Mars. Climate. Orbiter.

 

[BvU]

Story goes: not NASA error, though...

 

[ecastro]

Did the original NASA code generate NaNs, or was that your code?

It was my code that generated NaNs. According to the user manual, the code should generate a numerical value.

Three words. Mars. Climate. Orbiter.

Ah, yes. I was reminded that there was a failed mission due to conversion errors. But, I suppose this code is used for the generation of satellite products and is accompanied by published papers.

 

[BvU]

Can you post your NaN disaster as a more concrete, traceable casus ?

 

[BvU]

Wasted some of my employers time dragging these source files into a DEC compiler (which, yuck, they call an Intel compiler nowadays) and got loads of errors (*), but some of them are pretty serious:
in AKTOOL.f subroutine dakg writes into u(30+) (up to 48) where the calling subroutine alkalbe declares uu(30), a(30)
Doesn't have to be a disaster (if the call isn't executed ever), but sure isn't very robust .

(*) well, 25 of them, and I compiled everything -- and there are clearly versions that probably shouldn't be used, like newMODIS ("subscript out of range: SR(13,...)" -- sr is declared SR (12,1501) ) etc.

I'm not familiar with gfortran, but if you can set compiler switches, you should do a maximum of checking.

 

[ecastro]

Can you post your NaN disaster as a more concrete, traceable casus ?

Do you mean the line of processes for the NaN?

The 'refet' variable must have a quantitative value. If I traced it correctly, the 'refet' variable was divided with the 'seb' variable, which I have checked to be non-zero, so it can't cause the NaN. Before this, it was also calculated with 'romeas2' and 'coef'. Now, 'romeas2' is NaN.

The 'romeas2' variable is a result from 'ratm2', 'rsurf', and 'tgtot', and 'ratm2' is NaN (I haven't checked the other two variables yet).

The 'ratm2' variable is from the result of 'romix', 'rorayl', 'tgtot', and 'tgp1', and 'romix' is also NaN (I haven't also checked the three variables).

Now, 'romix' is an output of the 'interp' subroutine. Near the beginning of the code, the subroutine calculates a coefficient 'alphar', which uses the 'roatm' variable. This part causes the NaN because it is 0/0.

Wasted some of my employers time dragging these source files into a DEC compiler (which, yuck, they call an Intel compiler nowadays) and got loads of errors (*), but some of them are pretty serious:
in AKTOOL.f subroutine dakg writes into u(30+) (up to 48) where the calling subroutine alkalbe declares uu(30), a(30)
Doesn't have to be a disaster (if the call isn't executed ever), but sure isn't very robust .

(*) well, 25 of them, and I compiled everything -- and there are clearly versions that probably shouldn't be used, like newMODIS ("subscript out of range: SR(13,...)" -- sr is declared SR (12,1501) ) etc.

I'm not familiar with gfortran, but if you can set compiler switches, you should do a maximum of checking.

Thank you for the trouble of compiling the codes. As of now, I do not use the subroutines you mentioned; I do not know for the other 24.

 

[BvU]

Do you mean the line of processes for the NaN?

Actually I meant what do you input to run the program (igeom, month, jday, tu and so forth). This hopefully sits in a file, because the reading in from readin=5 goes on and on and on... ?
For now I use the input file as in appendix II of the manual.

In my case it bombs out in DISCOM (I think before Interp is called).

forrtl: severe (408): fort: (2): [B]Subscript #1 of the array WLDIS has value 21 wh
ich is greater than the upper bound of 20[/B]

Image              PC        Routine            Line        Source
Console1.exe       00491B20  Unknown               Unknown  Unknown
Console1.exe       003F99BF  Unknown               Unknown  Unknown
Console1.exe       003D7207  Unknown               Unknown  Unknown
Console1.exe       003D7963  Unknown               Unknown  Unknown
Console1.exe       0030B926  _DISCOM                    65  DISCOM.f
Console1.exe       002D5C71  _MAIN__                  1575  main.f
Console1.exe       00498B73  Unknown               Unknown  Unknown
Console1.exe       00441859  Unknown               Unknown  Unknown
Console1.exe       0044171F  Unknown               Unknown  Unknown
kernel32.dll       768A336A  Unknown               Unknown  Unknown
ntdll.dll          774D9902  Unknown               Unknown  Unknown
ntdll.dll          774D98D5  Unknown               Unknown  Unknown

No wonder:

      common /sixs_disc/ roatm(3,20),dtdir(3,20),dtdif(3,20),
     a utdir(3,20),utdif(3,20),sphal(3,20),wldis(20),trayl(20),
     a traypl(20),rqatm(3,20),ruatm(3,20)
...
c     computation of all scattering parameters at wavelength
c     discrete values,so we can interpolate at any wavelength

      do 50 l=1,20
        wl=wldis(l)
        if ((wlsup.lt.wldis(1)).and.(l.le.2)) goto 30
        if (wlinf.gt.wldis(20).and.(l.ge.19)) goto 30
        if ((l.lt.20).and.(wldis(l).lt.wlinf).and.
     a     (wldis(l+1).lt.wlinf)) goto 50
        if ((l.gt.1).and.(wldis(l).gt.wlsup).and.
     a      (wldis(l-1).gt.wlsup)) goto 50

In the last iteration it picks up wldis(L+1) for value 20 of L ! A crash with my compiler and a random number ( well, whatever is in trayl(1) ) if array bounds aren't checked, which is a dumb thing for NASA not to do !
( mitigating circumstances: some compilers test the first term and don't continue if it gives FALSE
but counting on that invites bad programming: if a further term depends on a subroutine call, such a call is never made )

Should be something like

      do 50 l=1,20
        wl=wldis(l)
        if ((wlsup.lt.wldis(1)).and.(l.le.2)) goto 30
        if (wlinf.gt.wldis(20).and.(l.ge.19)) goto 30

        if ((l.lt.20) then

          if ((wldis(l).lt.wlinf).and.
     a     (wldis(l+1).lt.wlinf)) goto 50

        endif

        if ((l.gt.1).and.(wldis(l).gt.wlsup).and.
     a      (wldis(l-1).gt.wlsup)) goto 50

Now I still crash very hard in DISCOM calling ATMREF calling OSPOL calling KERNELPOL -- not far before the point where you crash. Probably relevant, because in the same loop, after the call to ATMREF the array roatm is filled.

On the other hand: the error I get is stack overflow, the moment the program wants to do the first step in KERNELPOL. I can't find a good reason why this happens.
So your tracing is very useful also. 0/0 makes one suspect some initialization hasn't been performed yet.

        coef=alog(wldis(lsup)/wldis(linf))
...
        alphar=alog(roatm(1,lsup)/roatm(1,linf))/ coef

I would check the values of rorayl, romix and roaero after the return from subroutine atmref

 

[ecastro]

Yes, the inputs I used are the ones presented in the appendix of the manual.

In the last iteration it picks up wldis(L+1) for value 20 of L ! A crash with my compiler and a random number ( well, whatever is in trayl(1) ) if array bounds aren't checked, which is a dumb thing for NASA not to do !
( mitigating circumstances: some compilers test the first term and don't continue if it gives FALSE
but counting on that invites bad programming: if a further term depends on a subroutine call, such a call is never made )

I did not notice this one... I thought in an AND statement when the first condition is false, the other conditions will not be processed. However, I was able to produce an output for this subroutine, and the undefined matrix is fixed. However, I still have some NaN values being generated at the INTERP subroutine, specifically the 'dtott' variable.

 

[BvU]

I thought in an AND statement when the first condition is false, the other conditions will not be processed

Generally true. Not safe to assume always true. But here it was a matter of not even passing the array bounds check ecxecuted by the compiler. "If it works, leave it as it is for now" is probably the motto.

Further 'bad' news:

Now I still crash very hard in DISCOM calling ATMREF calling OSPOL calling KERNELPOL -- not far before the point where you crash. Probably relevant, because in the same loop, after the call to ATMREF the array roatm is filled.

On the other hand: the error I get is stack overflow, the moment the program wants to do the first step in KERNELPOL. I can't find a good reason why this happens.

I got around the stack overflow crash (by raising the Stack Reserve Size from 1 to 4M) and now I can't reproduce your problem. It does thihngs and produces a clean output:

******************************* 6SV version 2.1 *******************************
*                                                                             *
*                       geometrical conditions identity                       *
*                       -------------------------------                       *
*                       user defined conditions                               *
*                                                                             *
*   month:  7 day :  23                                                       *
*   solar zenith angle:   40.00 deg  solar azimuthal angle:      100.00 deg   *
*   view zenith angle:    45.00 deg  view azimuthal angle:        50.00 deg   *
*   scattering angle:    146.49 deg  azimuthal angle difference:  50.00 deg   *
*                                                                             *
*                       atmospheric model description                         *
*                       -----------------------------                         *
*           atmospheric model identity :                                      *
*             user defined water content : uh2o= 3.000 g/cm2                  *
*             user defined ozone content : uo3 = 3.500 cm-atm                 *
*           aerosols type identity :                                          *
*             user-defined aerosol model:                                     *
*                           0.250 % of dust-like                              *
*                           0.250 % of water-soluble                          *
*                           0.250 % of oceanic                                *
*                           0.250 % of soot                                   *
*           optical condition identity :                                      *
*               visibility :  8.49 km  opt. thick. 550 nm :  0.5000           *
*                       spectral condition                                    *
*                       ------------------                                    *
*           avhrr 1 (noaa9)                                                   *
*               value of filter function :                                    *
*                wl inf= 0.530 mic   wl sup= 0.810 mic                        *
*                                                                             *
*                       Surface polarization parameters                       *
*                       ----------------------------------                    *
*                                                                             *
*                                                                             *
* Surface Polarization Q,U,Rop,Chi    0.00000  0.00000  0.00000    11.44      *
*                                                                             *
*                                                                             *
*                       target type                                           *
*                       -----------                                           *
*           homogeneous ground                                                *
*                                                                             *
*                       target elevation description                          *
*                       ----------------------------                          *
*           ground pressure  [mb]  989.01                                     *
*           ground altitude  [km]-0.200                                       *
*                                                                             *
*                       plane simulation description                          *
*                       ----------------------------                          *
*           plane  pressure          [mb]  657.54                             *
*           plane  altitude absolute [km]  3.500                              *
*                atmosphere under plane description:                          *
*                ozone content             0.082                              *
*                h2o   content             2.389                              *
*               aerosol opt. thick. 550nm  0.281                              *
*                                                                             *
*                        atmospheric correction activated                     *
*                        --------------------------------                     *
*                        BRDF coupling correction                             *
*           input apparent reflectance :  0.100                               *
*                                                                             *
*******************************************************************************
*******************************************************************************
*                                                                             *
*                         integrated values of  :                             *
*                         --------------------                                *
*                                                                             *
*       apparent reflectance  0.0330946  appar. rad.(w/m2/sr/mic)   12.751    *
*                   total gaseous transmittance  0.675                        *
*                                                                             *
*******************************************************************************
*                                                                             *
*                         coupling aerosol -wv  :                             *
*                         --------------------                                *
*           wv above aerosol :   0.033     wv mixed with aerosol :   0.033    *
*                       wv under aerosol :   0.033                            *
*******************************************************************************
*                                                                             *
*                         integrated values of  :                             *
*                         --------------------                                *
*                                                                             *
*       app. polarized refl.     Inf    app. pol. rad. (w/m2/sr/mic)    0.065 *
*             direction of the plane of polarization 45.00                    *
*                   total polarization ratio       Inf                        *
*                                                                             *
*******************************************************************************
*                                                                             *
*                         int. normalized  values  of  :                      *
*                         ---------------------------                         *
*                      % of irradiance at ground level                        *
*     % of direct  irr.    % of diffuse irr.    % of enviro. irr              *
*               0.773               0.221               0.005                 *
*                       reflectance at satellite level                        *
*     atm. intrin. ref.   environment ref.  target reflectance                *
*               0.015               0.004               0.014                 *
*                                                                             *
*                         int. absolute values of                             *
*                         -----------------------                             *
*                      irr. at ground level (w/m2/mic)                        *
*     direct solar irr.    atm. diffuse irr.    environment  irr              *
*             453.574             127.136               3.157                 *
*                      rad at satel. level (w/m2/sr/mic)                      *
*     atm. intrin. rad.    environment rad.    target radiance                *
*               5.650               1.633               5.468                 *
*                                                                             *
*                                                                             *
*          int. funct filter (in mic)              int. sol. spect (in w/m2)  *
*             0.1174545                                 185.589               *
*                                                                             *
**************************************************************************************************************************************************************
*                                                                             *
*                          integrated values of  :                            *
*                          --------------------                               *
*                                                                             *
*                             downward        upward          total           *
*      global gas. trans. :     0.68965        0.97248        0.67513         *
*      water   "     "    :     0.98573        0.98623        0.97592         *
*      ozone   "     "    :     0.70609        0.99079        0.70008         *
*      CO2     "     "    :     1.00000        1.00000        1.00000         *
*      oxyg    "     "    :     0.99344        0.99533        0.99179         *
*      no2     "     "    :     1.00000        1.00000        1.00000         *
*      ch4     "     "    :     1.00000        1.00000        1.00000         *
*      co      "     "    :     1.00000        1.00000        1.00000         *
*                                                                             *
*                                                                             *
*      rayl.  sca. trans. :     0.96494        0.93809        0.90520         *
*      aeros. sca.   "    :     0.72090        0.82111        0.59194         *
*      total  sca.   "    :     0.69208        0.81074        0.56110         *
*                                                                             *
*                                                                             *
*                                                                             *
*                             rayleigh       aerosols         total           *
*                                                                             *
*      spherical albedo   :     0.04939        0.04918        0.06820         *
*      optical depth total:     0.05550        0.42021        0.47570         *
*      optical depth plane:     0.01848        0.23606        0.25454         *
*      reflectance I      :     0.01098        0.01327        0.02176         *
*      reflectance Q      :     0.00118        0.00037        0.00122         *
*      reflectance U      :    -0.00156        0.00000       -0.00173         *
*      polarized reflect. :     0.00195        0.00037        0.00212         *
*      degree of polar.   :       17.77           2.76           9.75         *
*      dir. plane polar.  :      -26.48           0.00         -27.43         *
*      phase function I   :     1.26026        0.27565        0.39051         *
*      phase function Q   :    -0.21911       -0.00611       -0.03096         *
*      phase function U   :    -1.19913       -0.15957       -0.28084         *
*      primary deg. of pol:    -0.17386       -0.02215       -0.07927         *
*      sing. scat. albedo :     1.00000        0.52284        0.57850         *
*                                                                             *
*                                                                             *
*******************************************************************************
*******************************************************************************
*******************************************************************************
*                        atmospheric correction result                        *
*                        -----------------------------                        *
*       input apparent reflectance            :    0.100                      *
*       measured radiance [w/m2/sr/mic]       :   38.529                      *
*       atmospherically corrected reflectance                                 *
*       Lambertian case :      0.22186                                        *
*       BRDF       case :      0.22186                                        *
*       coefficients xa xb xc                 :  0.00685  0.03871  0.06820    *
*       y=xa*(measured radiance)-xb;  acr=y/(1.+xc*y)                         *
*       coefficients xap xb xc                :  2.639795  0.038714  0.068196 *
*       y=xap*(measured reflectance)-xb;  acr=y/(1.+xc*y)                     *
*******************************************************************************

Close to the example in the manual. Not exactly, but I suspect that one is a bit obsolete (something to check, ... later).

Question now remains how to get your version to run to a comparable end without crashing.
Increasing the stack size might be worth investigating ?

 

[ecastro]

The output is strange. I got everything correct for the first part of the program (before the apparent reflectance values). I am now able to produce a numerical value, but an incorrect one, way too large compared with the ones shown in the manual. It would seem that the main problem comes from the ISO subroutine; I replaced the 'continue's with 'break's, so the iteration breaks rather than continue to the next loop, which I have corrected.

Increasing the stack size might be worth investigating ?

What do you mean by this? I am using the MatLab software, by the way.