Bill Sacks < [email protected] > - 2013-05-01 15:24:28 -0600
Bugzilla Id: 1679
Bugzilla CC: [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected],
This started with a search through the CNFireMod code for potentially-incorrect uses of weights, now that I have introduced 'active' flags. But in looking for this, I noticed a number of questionable things in this code, so decided to review the entire module. Here are some problems I found.
Some of these are cosmetic (e.g., variable names); others are not bugs per se, but make the code fragile to future changes; and others look like bugs (and if they aren't bugs, then they need some explanation about why they're correct even though they look incorrect).
(1) Many variables are declared as "! local pointers to implicit in scalars", when in fact they are implicit out or in/out
(2) There are a number of cryptic variable names, some with no description. For example:
real(r8) :: lh !
real(r8) :: fs !
real(r8) :: ig !
(3) There are a number of places where the code checks whether the current pft is crop or natural veg. This is done like this:
! For crop veg types
if( ivt(p) > nc4_grass )then
! For natural vegetation (non-crop)
if( ivt(p) >= ndllf_evr_tmp_tree .and. ivt(p) <= nc4_grass )then
This is very fragile code, which will break if anyone adds new natural PFTs to the beginning or end of the list. A possible replacement for crop is - I think - to use npcropmin & npcropmax (defined in pftvarcon.F90) - though I'm not sure whether that does the right thing when running without the specific crop types?
(4) Similarly, checks of the type of pft:
if( ivt(p) .ge. nbrdlf_evr_shrub )then !for shurb and grass
...
else ! for trees
This assumes that people adding PFTs in the future will add them in the "right" place, where "right" is determined by the logic embedded here. This should probably be handled by adding constants to the pft physiology file.
The idea of someone wanting to add PFTs is not purely theoretical: Andy Jones just told me today that he is adding some PFTs for some of his work. So I am leery of code that makes this difficult or error-prone.
(5) There are some seemingly ad-hoc conditionals, some of which seem incorrect. In particular, there seem to be terms added to conditionals that don't need to (or shouldn't) be there. For example:
(a)
! For crop PFT's
if( ivt(p) > nc4_grass .and. wtcol(p) > 0._r8 .and. leafc_col(c) > 0._r8 )then
fuelc_crop(c)=fuelc_crop(c) + (leafc(p) + leafc_storage(p) + &
leafc_xfer(p))*wtcol(p)/cropf_col(c) + &
totlitc(c)*leafc(p)/leafc_col(c)*wtcol(p)/cropf_col(c)
end if
Why should this only be done if leafc_col(c) > 0? Note that the first parenthesized term can be non-zero even if leafc_col(c) == 0, because leafc_storage or leafc_xfer could be > 0.
Also:
(b)
if( .not. shr_infnan_isnan(btran2(p)) .and. btran2(p) .le. 1._r8 )then
Why is this check needed? Since this code is in a loop over filter_soilc points, it shouldn't include points that have NaN values; and I can't see how btran2 could be greater than 1 -- and if it ever were greater than 1, I would question whether it's really correct to exclude points that have values greater than 1: is it really correct to treat points with values of 1.0000000001 (maybe due to rounding error??) fundamentally differently from points with values of 1.0?
(6) More intermediate variables (and/or comments) are needed to make it more clear what this code is doing, so people other than the author can modify it in the future. For example, in the example in #5a, I have no idea what this term is trying to accomplish:
totlitc(c)*leafc(p)/leafc_col(c)*wtcol(p)/cropf_col(c)
As another example, this equation is very hard to read:
baf_peatf(c) = boreal_peatfire_c*exp(-SHR_CONST_PI*(max(wf2(c),0._r8)/0.3_r8))* &
max(0._r8,min(1._r8,(tsoi17(c)-SHR_CONST_TKFRZ)/10._r8))*peatf_lf(c)* &
(1._r8-fsat(c))
(7) There are many places where long expressions do two things at once: (1) compute some pft-level variable with a complex expression, and (2) average from pft to column -- all in a single line of code. The example in #5a, above, illustrates this.
The code would be clearer and easier to modify if these expressions were split into two pieces: (1) first compute pft-level variable for all quantities, and then (2) average to column, preferably via a call to p2c.
(6) The use of cwtgcell in this code is suspicious to me:
if( ivt(p) == nbrdlf_evr_trp_tree .and. wtcol(p) .gt. 0._r8 )then
trotr1_col(c)=trotr1_col(c)+wtcol(p)*cwtgcell(c)
end if
if( ivt(p) == nbrdlf_dcd_trp_tree .and. wtcol(p) .gt. 0._r8 )then
trotr2_col(c)=trotr2_col(c)+wtcol(p)*cwtgcell(c)
end if
if( ivt(p) == nbrdlf_evr_trp_tree .or. ivt(p) == nbrdlf_dcd_trp_tree )then
if(lfpftd(p).gt.0._r8)then
dtrotr_col(c)=dtrotr_col(c)+lfpftd(p)*cwtgcell(c)
end if
end if
There appears to be no other science code (i.e., non-infrastructure code) that depends on weights on the grid cell. I believe that all code is intended to compute per-area quantities, so a column's weight on the grid cell should be irrelevant in computing a column-level quantity.
If this code is truly correct, it needs more explanation.
(7) Why does the normalization differ for these two averages from pft to col:
rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + &
frootc_xfer(p) + deadcrootc(p) + &
deadcrootc_storage(p) + deadcrootc_xfer(p) + &
livecrootc(p)+livecrootc_storage(p) + &
livecrootc_xfer(p))*wtcol(p)
fsr_col(c) = fsr_col(c) + fsr_pft(ivt(p))*wtcol(p)/(1.0_r8-cropf_col(c))
i.e., why does the second normalize by (1-cropf_col), but the first does not?
(8) There are lots of magic numbers. For example:
fire_m = exp(-SHR_CONST_PI *(m/0.69_r8)**2)*(1.0_r8 - max(0._r8, &
min(1._r8,(forc_rh(g)-30._r8)/(70._r8-30._r8))))* &
min(1._r8,exp(SHR_CONST_PI*(forc_t(g)-SHR_CONST_TKFRZ)/10._r8))
lh = 0.0035_r8*6.8_r8*hdmlf**(0.43_r8)/30._r8/24._r8
fs = 1._r8-(0.01_r8+0.98_r8*exp(-0.025_r8*hdmlf))
ig = (lh+forc_lnfm(g)*0.25_r8)*(1._r8-fs)*(1._r8-cropf_col(c))
(9) In this code and some following code, normalizations are sometimes done by (1 - cropf_col), and sometimes by lfwt:
if( ivt(p) .ge. nbrdlf_evr_shrub )then !for shurb and grass
lgdp_col(c) = lgdp_col(c) + (0.1_r8 + 0.9_r8* &
exp(-1._r8*SHR_CONST_PI* &
(gdp_lf(c)/8._r8)**0.5_r8))*wtcol(p) &
/(1.0_r8 - cropf_col(c))
lgdp1_col(c) = lgdp1_col(c) + (0.2_r8 + 0.8_r8* &
exp(-1._r8*SHR_CONST_PI* &
(gdp_lf(c)/7._r8)))*wtcol(p)/lfwt(c)
lpop_col(c) = lpop_col(c) + (0.2_r8 + 0.8_r8* &
exp(-1._r8*SHR_CONST_PI* &
(hdmlf/450._r8)**0.5_r8))*wtcol(p)/lfwt(c)
Is lfwt = (1 - cropf_col)?? If so, these normalizations should be made consistent. If not, why do they differ, and why are some normalizations done by one factor and others by another?
(10) The computations of lgdp_col, lgdp1_col and lpop_col are hard to decipher, and could use some explanatory comments and/or code cleanup (e.g., use of intermediate variables) to make it more clear what's going on here.
(11) The following code needs some explanation:
if (fpftdyn /= ' ') then !true when landuse data is used
do fc = 1,num_soilc
c = filter_soilc(fc)
if( dtrotr_col(c) .gt. 0._r8 )then
if( kmo == 1 .and. kda == 1 .and. mcsec == 0)then
lfc(c) = 0._r8
end if
if( kmo == 1 .and. kda == 1 .and. mcsec == dt)then
lfc(c) = dtrotr_col(c)dayspyrsecspday/dt
end if
else
lfc(c)=0._r8
end if
end do
end if
It looks like this is trying to re-initialize lfc at the beginning of each year. However, I'm concerned about whether this will behave correctly if drotr_col(c) == 0 on the first and/or second timesteps of the year. In this case the variables won't be initialized for that year. Maybe that's okay - it's just not clear to me if that's what's intended.
(12) It is odd to me that deforestation fires cannot happen on Jan 1 at midnight:
!
! if landuse change data is used, calculate deforestation fires and
! add it in the total of burned area fraction
!
if (fpftdyn /= ' ') then !true when landuse change data is used
if( trotr1_col(c)+trotr2_col(c) > 0.6_r8 )then
if(( kmo == 1 .and. kda == 1 .and. mcsec == 0) .or. &
dtrotr_col(c) <=0._r8 )then
fbac1(c) = 0._r8
farea_burned(c) = baf_crop(c)+baf_peatf(c)
but perhaps there is some good reason for this
(13) Subroutine CNFireFluxes hard-codes a lot of information about the structure of the CN pools. I believe the same is done elsewhere in the CN code. This makes it very hard for anyone to add or remove a carbon or nitrogen pool, because they have to understand and modify code spread all over the model. I'm not sure what the solution is, but I feel like the current implementation is not sustainable, so this is worth some thought.
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent