Uninitialised variables

[dramanica]

I have been playing a little bit with the biome4 subroutine (isolating it from the driver and feeding it some numbers directly), and it seems to be sometimes returning slightly different results when feeding it the same inputs. It is usually the first run that gives a discordant result, and then the output stabilises.

Compilation with -finit-local-zero leads to yet slightly different results, suggesting that the issue might be unitialised variables. Compilation with -Wuninitialized points out a few places in the code that seem problematic.
An easy one to fix (I think) is on line 2022:

c      Set non-co2-dependent parameters
       ts = o2 / (2.*tao)
       z  = cmass*jtoe*dsun*fpar*twigloss*tune

c      Need to change the partial pressure of o2 also:
       mfo2=slo2/1e5
       o2=p*mfo2 

Where the two blocks should be swapped (as they are in another subroutine) so that we first define o2, and then use it to compute ts.

I get a few others, which are a bit more tricky to manage. If @jedokaplan is game, we could have a go at going through each of them and find a fix, but I'd need a bit of wisdom from someone who understands the innards of the model better than I do.

[jedokaplan]

I am sure that there are many bugs in the code like this one. I haven't worked on it since 1999! Sure, pedantic compilation would be a good way to catch the easy ones. Happy to assist @dramanica, let me know when you would like to work on this and what questions you have.

[dramanica]

@jedokaplan It is certainly vintage code, but people keep using it, so a little love won't hurt. I have written an R driver for it that would make it even easier to run it, but that brought up the numerical inconsistency issues. So, let's see whether we can deal with those.
So, I fixed the o2 define discussed above. I then set -fpe0 (in ifort) to start catching problematic operations (and thus eventually track down any undefined variable) . The next bit of code I stumbled upon is in hydrology ,
on line 2371:

       if(demand.gt.supply)then
        a = (1. - supply/(deq(d)*alfam))
        if (a.lt.0.0) a=0.0
        gsurf    =  -gm*log(a)
        aet      =  supply
        gc = gsurf - gmin
c       Constrain gc value to zero!
        if (gc.le.0.0) then
         gc=0.0
         wilt=.true.
        end if
       end if

the problem is on the 3rd line. If a<0, then we set a to 0.0, but then we have log(a), which is -infinity. this propagates to gc = infinity when we get out of the if statement.
Suggestion for a fix that keeps the maths sane?
An Addendum to this. if demand.le.supply, gc is not defined, but it is used later on in the function. What should we initialise it as in hydrology?

[dramanica]

And another little puzzle for you, this time in growth:

       if (meangc(m).ne.0) then
        mcount=mcount+1
        anngasum=anngasum+(mgpp(m)/meangc(m))  !new line for A/g!
       end if

in some instances meangc(m) can be tiny (10^-9). This makes anngasum huge, and things go wrong after that.
Would it make sense to have a cap on how small meangc(m) can be. Something like meangc(m).gt.0.0001 (or a more sensible number; from a quick look, we do go down to 10^-2). Or is that the wrong way to approach it and we just need to cap anngsum?

[dramanica]

Ok, here is an early Christmas present of sorts for @jedokaplan. I have gone through the whole source code of biome4, and tracked down all instances of uninitialised variables and instances of numerical problems. The attached code includes flags for all the locations that lead to problems. There are two flags: @FIXED where I felt I knew how to fix the problem properly (but it would be good to have a check), and @NOTFIXED where I bodged something for the code to work (e.g. initialise a variable to zero) but I felt that my solution might not be ideal (so, they require a better understanding of the code to make sure that the fix is correct).
@jedokaplan I don't think it should take you too long to go through it: I have made notes on the logic of the problem for each location of the code, and suggested possible fixes (or options on what I thought would be required to fix). You can either do that in your own time, or we can plan for a Zoom chat and I can walk you through each of the issues.
Note that I wrapped the subroutines into a module for my convenience, but otherwise kept the old F77 formatting.

biome4.zip

[jedokaplan]

Thanks Andrea, I will see if I can get a bugfix out before the holidays.

[dramanica]

Great, let me know if anything is unclear. I have additional notes for each of the issues, with a lat and long that shows each of the problem in action.