Fortran How does Fortran know the Inputs and Outputs of Subroutines

[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.

 

[jedishrfu]

Another note on the common block, it was used in cases where there were large numbers of variables where it became painful to pass them all in a subroutine call. Doing it the functional way i.e. Not the common block way would often take a lot more time to run the program. In other words common blocks could speed up computationally intense programs due to the register indirection used in subroutine calls.

When common block became hell was when you had the same common block referenced in two or more subroutines use different names for the same variables. Invariably someone would add a new variable in the middle of the block and forget to add it to all the other same named blocks they would use fortrans implicit datatyping convention meaning one block thought it was an integer and another thought it was a float.

Common blocks are a primitive version of the union+structs used in C code for low level data manipulation.

 

[BvU]

The output is strange.

The one I posted in #34 or the one you didn't ?

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.

Value for what ? The integrated apparent reflectance 0.0330946 ?

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.

This is hard to follow. There's a lot of loops in that routine. Concrete code ? Is this in fortran or in MATLAB ?

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

The arguments for fortran subroutine calls all go on the stack, so if you are a dozen levels deep and the arguments are huge double precision arrays you can run out of stack space.

I am using the MatLab software, by the way.

Are you saying you ported some 2 MB of fortran to MATLAB 'by the way' ?

[edit]PS In the output I posted in #34 I found the 11.44 for Chi (surface Polarization) rather suspect. It's due to a dangerous bug: main.f reads in inhomo at label 30 and jumps to label 31 if it's > 0. idirec is not initialized in that case, but it is checked later on: after label 37 its uninitialized value (!) determines if ropq and ropu (surface Polarization Q and U) are initialized to 0 or not. And in the not case it continues to do things based on the also uninitialized value of ibrdf .
Basically you are in the other branch of the tree and don't care much (unless some xx variable gets a value you later break your neck over ) but it looks soooo ugly ! (Someone should alert the authors about this ?).

 

[jbriggs444]

The arguments for fortran subroutine calls all go on the stack, so if you are a dozen levels deep and the arguments are huge double precision arrays you can run out of stack space.

In the compilers and machine architecture I am familiar with, Fortran uses call by reference and static allocation for local storage. There is nothing to cause excessive stack utilization. However, my experience is decades old. Possibly there are re-entrant Fortran implementations today with stack-allocated local storage.

 

[BvU]

Hi JB, I'm a VAX Fortran 77 native speaker too, but nowadays only have access to an Intel compiler on Windows. It appears to use the stack for actual arguments and the default size is only 1 MB. SInce upping it solved my problem I'm inlclined to believe the documentation.

 

[ecastro]

The one I posted in #34 or the one you didn't ?

It is the 34th post.

Value for what ? The integrated apparent reflectance 0.0330946 ?

Yes, and also the apparent radiance. I have values reaching magnitudes of 10, if I remember correctly.

This is hard to follow. There's a lot of loops in that routine. Concrete code ? Is this in fortran or in MATLAB ?

Yes, I re-checked the code and there seems to loops in my MatLab code that are wrong. I will try to correct this first before sharing the code.

PS In the output I posted in #34 I found the 11.44 for Chi (surface Polarization) rather suspect. It's due to a dangerous bug: main.f reads in inhomo at label 30 and jumps to label 31 if it's > 0. idirec is not initialized in that case, but it is checked later on: after label 37 its uninitialized value (!) determines if ropq and ropu (surface Polarization Q and U) are initialized to 0 or not. And in the not case it continues to do things based on the also uninitialized value of ibrdf

I also encountered this 'check', the 'idirec' variable has no value. However, I have set to zero, and the values I had seemed correct.

 

[ecastro]

I have another question regarding this. The subroutine 'atmref' calls the 'opsol' subroutine with three of its outputs have the same name (the first call of the subroutine, for example), 'rolutd'. In the code of the 'opsol' subroutine, these three outputs are three different variables, so which one of these will have be considered 'rolutd'?

 

[jbriggs444]

I have another question regarding this. The subroutine 'atmref' calls the 'opsol' subroutine with three of its outputs have the same name (the first call of the subroutine, for example), 'rolutd'. In the code of the 'opsol' subroutine, these three outputs are three different variables, so which one of these will have be considered 'rolutd'?

All three.

The preferred terminology for the variables passed in a procedure call is "actual parameters" (referring to what is placed into the parameter list by the caller) and "formal parameters" or "dummy parameters" (referring to the parameter names used by the subroutine).

If you use 'rolutd' three times as an actual parameter and if the corresponding formal parameters in the subroutine are, for instance, 'a', 'b' and 'c' then a, b and c will refer to the same block of storage. [Here I assumed call by reference. If rolutd is a scalar variable rather than an array and if the compiler uses copy-in/copy-back semantics for scalars then you would have three different [notional] memory locations populated with the same initial value].

Note that if the actual parameters are variables or array elements (as opposed to being a whole array) and if the subroutine modifies anyone of the corresponding formal parameters then the standard forbids you to use the same actual parameter more than once. Without having read the relevant section of the standard or the minds of its creators, the idea is to give maximum flexibility to the compiler writers when choosing a parameter passing mechanism. Call by reference, and call by value with copy-in/copy-back would give different results if the programmer did that. So the standard forbids the programmer to do that -- or declares that the results of doing so are "undefined" or that a program that depends on the results is "erroneous".

See the section on "aliasing of dummy arguments" in https://www.ibiblio.org/pub/languages/fortran/ch1-8.html

 

[ecastro]

I am a bit lost. So, for example I have a subroutine ABC with dummy parameters as (x, y, z) and more or less looks like this:

subroutine ABC (x, y, z)
...
x = 1
y = 2
z = 3
...
end

Then this subroutine is called in the main program,

call ABC (d, d, d)

Then, what is the value of 'd'?

 

[jbriggs444]

I am a bit lost. So, for example I have a subroutine ABC with dummy parameters as (x, y, z) and more or less looks like this:

subroutine ABC (x, y, z)
...
x = 1
y = 2
z = 3
...
end

Then this subroutine is called in the main program,

call ABC (d, d, d)

Then, what is the value of 'd'?

Undefined. You have written an illegal program.

 

[ecastro]

Even if 'x', 'y', and 'z' are arrays in the subroutine?

 

[jbriggs444]

Even if 'x', 'y', and 'z' are arrays in the subroutine?

How are you going to make sense of the statement: "x = 1" in such a case?

 

[DrClaude]

I am a bit lost. So, for example I have a subroutine ABC with dummy parameters as (x, y, z) and more or less looks like this:

subroutine ABC (x, y, z)
...
x = 1
y = 2
z = 3
...
end

Then this subroutine is called in the main program,

call ABC (d, d, d)

Then, what is the value of 'd'?

Undefined. You have written an illegal program.

I beg to differ. This is a perfectly valid call in Fortran, and the value of d after the call will be 3. All three variables, x, y and z, point to the same memory location and the instructions are expected to be executed in order

This no different from

equivalence (x, y), (x, z)
call ABC (x,y,z)

That said, if ABC does something else to/with x, y, and z, the outcome can become undefined.

 

[jbriggs444]

I beg to differ. This is a perfectly valid call in Fortran, and the value of d affecter the call will be 3. All three variables, x, y and z, point to the same memory location and the instructions are expected to be executed in order

I am unable to locate any section in the standard which explicitly declares this usage improper. However, I would not want to trust that a compiler would produce that result without having tested it first.

Optimizing compilers can make unexpected choices. With separate compilation, the compiler for the subroutine may not be aware that its dummy parameters are aliases for the same actual parameter.

Edit: Found it.

https://www.fortran.com/F77_std/rjcnf0001-sh-15.html#sh-15.9.3.6

15.9.3.6 Restrictions on Association of Entities.
If a subprogram reference causes a dummy argument in the referenced subprogram to become associated with another dummy argument in the referenced subprogram, neither dummy argument may become defined during execution of that subprogram. For example, if a subroutine is headed by
SUBROUTINE XYZ (A,B)
and is referenced by
CALL XYZ (C,C)
then the dummy arguments A and B each become associated with the same actual argument C and therefore with each other. Neither A nor B may become defined during this execution of subroutine XYZ or by any procedures referenced by XYZ.

 

[Mark44]

I agree with @DrClaude. Inside the ABC routine, x is set to 1, causing d in the main program to be set to 1. Then y is set to 2, causing d to be reset to 2. Then z = set to 3, causing d to be reset again to 3.

 

[jbriggs444]

I agree with @DrClaude. Inside the ABC routine, x is set to 1, causing d in the main program to be set to 1. Then y is set to 2, causing d to be reset to 2. Then z = set to 3, causing d to be reset again to 3.

The standard disagrees.

 

[DrClaude]

The standard disagrees.

Where?

2.3.4 Execution sequence

Execution of an executable program begins with the first executable construct of the main program. The execution of a main program or subprogram involves execution of the executable constructs within its scoping unit. When a procedure is invoked, execution begins with the first executable construct appearing after the invoked entry point. With the following exceptions, the effect of execution is as if the executable constructs are executed in the order in which they appear in the main program or subprogram until a STOP, RETURN, or END statement is executed. The exceptions are:

1. (1) Execution of a branching statement (8.2) changes the execution sequence. These statements explicitly specify a new starting place for the execution sequence.
2. (2) IF constructs, CASE constructs, and DO constructs contain an internal statement structure and execution of these constructs involves implicit (i.e., automatic) internal branching. See Section 8 for the detailed semantics of each of these constructs.
3. (3) Alternate return and END=, ERR=, and EOR= specifiers may result in a branch.
   

Internal subprograms may precede the END statement of a main program or a subprogram. The execution sequence excludes all such definitions.

(emphasis added)

 

[Mark44]

I stand corrected... Thanks for digging this up, @jbriggs444.

 

[jbriggs444]

Where?

https://www.fortran.com/F77_std/rjcnf0001-sh-15.html#sh-15.9.3.6

[I'd also editted that into response #48 earlier]

 

[DrClaude]

I also stand corrected. Thanks!

 

[Mark44]

Here's how C deals with a very similar situation.

// Test.c
#include <stdio.h>
void func(int *, int *, int *);

int main(void)
{
   int y = 1;
   func(&y, &y, &y);
   printf("y = %d", y);
   return 0;
}

void func(int * a, int * b, int * c)
{
   *a = 2;
   *b = 3;
   *c = 4;
}

Output: y = 4

 

[BvU]

I have another question regarding this. The subroutine 'atmref' calls the 'opsol' subroutine with three of its outputs have the same name (the first call of the subroutine, for example), 'rolutd'. In the code of the 'opsol' subroutine, these three outputs are three different variables, so which one of these will have be considered 'rolutd'?

I agree with @jbriggs444 that the standard doesn't allow it. But my fortran compiler (Intel visual fortran 2011) let's it pass and indeed, OSPOL (not opsol) is called twice as you describe. The effect -- IF the code to write to the rolut(m,l) is executed -- is as @DrClaude describes and the lookup table would get completely wrong values. a disastrous bug. However, in the example we use,

C Look up table generation
      do m=1,mu
      do l=1,nfilut(m)
      phimul=filut(m,l)*pi/180.
      rolut(m,l)=rolut(m,l)+delta0s*i3(m)*cos(is*(phimul+pi))
      rolutq(m,l)=rolutq(m,l)+delta0s*q3(m)*cos(is*(phimul+pi))
      rolutu(m,l)=rolutu(m,l)+delta0s*u3(m)*sin(is*(phimul+pi))
      enddo
      enddo
C end of look up table generation

nfilut(m) is all zeroes in all calls from atmref: it looks as if the reference table is not used for this set of inputs. Pretty inefficient (or a bug!)

I can't oversee what the first two calls to ospol from atmref are supposed to do and if they should be the same as the third one, or instead need dummy arrays (of the proper size) as the last two arguments.

 

[BvU]

@ecastro : my advice is to get into contact with the authors. If they published the stuff in 2015 that means there are still users and perhaps also developers working with/on this stuff. If I were one of them, I would be extremely grateful for this kind of detailed feedback. And perhaps they will also be interested in what you are trying to make !

 

[ecastro]

Thank you all for your help. I also wanted to contact the authors, but I wanted to confirm if it is a programming error or if I understand how Fortran works. I just found out that 'rolutd' is not used in 'atmref' (like a dummy output), that is maybe why the compiler let it pass?

 

[jbriggs444]

Thank you all for your help. I also wanted to contact the authors, but I wanted to confirm if it is a programming error or if I understand how Fortran works. I just found out that 'rolutd' is not used in 'atmref' (like a dummy output), that is maybe why the compiler let it pass?

Compilers do not detect all standards violations. Their job is simpler than that -- to take compliant programs and generate object code that will produce compliant results.

In the development environment that I am familiar with, the Fortran compiler would (or could) do separate compilation of each subroutine. When compiling 'atmref', the compiler would have knowledge of where and how 'atmref' will be called. So the compiler would not generate a compile time error message for 'atmref'. When compiling the caller of 'atmref', the compiler would have no knowledge of the inner workings of 'atmref'. So the compiler would not generate a compile time error message when compiling the caller.

It would be possible for the compiler to build in run time checks to detect and flag attempts to modify aliased dummy parameters. But the standard does not require that this be done. The standard does not specify the required behavior when running non-compliant code.

However, it is also plausible that your re-examination of the code is correct, that there is no attempt to modify an aliased dummy variable and that this particular standards violation does not exist in the code.

 

[aphirst]

I would just like to chip into contest the frequently asserted claim that "Fortran is a language where calls are intrinsically by-reference", by "referring" to a blog post by Steve Lionel a.k.a. Doctor Fortran (formerly?) of Intel. https://software.intel.com/en-us/blogs/2009/03/31/doctor-fortran-in-ive-come-here-for-an-argument

On a fundamental level, argument association is pretty simple. Although most people would tell you that Fortran uses "pass by reference", meaning that the address of the actual argument is used for the dummy argument, that's not quite true. The standard says that argument association is "usually similar to call by reference" and then adds "The actual mechanism by which this happens is determined by the processor." [Fortran 2003 Note 12.22]

There are some cases where what actually gets passed is not the original argument. The most common case for this is when a non-contiguous array pointer or a array slice is passed to an argument for which the compiler does not see an explicit interface specifying that the corresponding dummy argument is an assumed-shape array. Because the called routine is expecting a contiguous array, the compiler needs to make a copy of the actual argument, pass the copy, and then copy back any changes. The compiler can warn you about this at run-time if you use the option /check:arg_temp_created (Windows) or -check arg_temp_created (Linux/Mac OS). The compiler will generate a run-time test to see if the argument is actually contiguous and skip the copy if it is.

He goes on to give some more examples then emphasises the importance of explicitly declaring INTENT, but I think the point is clear.