Recently I had to implement a scientific paper The influence of vegetation, fire spread and fire behaviour on biomass burning and trace gas emissions: results from a process-based model into Fortran code. However, several challenges came up. Specially, edge cases not dealt within several equations. As an example:
As seen, the equation shown above is not complicated, as it is made of well known operations.
In order to illustrate what’s being expressed:
Where:
- (Input) CL: crown length in meters
- (Input) Height: tree height in meters
- (Input) Scorch Height: maximum height of the flames in meters
- (Output) CK: crown scorch as fraction from 0 to 1
Implementing such function, as seen:
FUNCTION calc_prop_ck(nc, scorch_height, height, crown_length) RESULT (prop_ck)
INTEGER, INTENT(in) :: nc !< total number of gridcells
REAL(wp), INTENT(in) :: scorch_height(nc) !< [m]
REAL(wp), INTENT(in) :: height(nc) !< [m]
REAL(wp), INTENT(in) :: crown_length(nc) !< [m]
REAL(wp) :: prop_ck(nc) !< output
prop_ck = (scorch_height - height + crown_length) / crown_length ! no (Eq. 17 from Thonicke et al. (2010))
END FUNCTION calc_prop_ck
As simple it might, there are some edge cases not accounted for:
- What if SH is taller than H? If that’s the case, there might be CK > 1
- What if SH is shorter than CL? If nothing’s changed, then CK < 0!
- What if CL is 0? An error would be triggered due to a division by 0
For that reason, such cases must be dealt. However, someone might be tempted to filter out those values above 1 and below 0. Nonetheless, such situation might not be ideal at all. Other inconsistencies could be masked. As a consequence, no error might be raised at all, preventing the developers to be aware of it.
For that reason, provide a proper solution for each case. Overall, this is how it could be done:
FUNCTION calc_prop_ck(nc, scorch_height, height, crown_length) RESULT (prop_ck)
INTEGER, INTENT(in) :: nc !< total number of gridcells
REAL(wp), INTENT(in) :: scorch_height(nc) !< [m]
REAL(wp), INTENT(in) :: height(nc) !< [m]
REAL(wp), INTENT(in) :: crown_length(nc) !< [m]
REAL(wp) :: prop_ck(nc) !< output
LOGICAL :: are_flames_taller_than_tree(nc)
LOGICAL :: are_flames_shorter_than_crown_length(nc)
are_flames_taller_than_tree(:) = (scorch_height(:) > height(:))
are_flames_shorter_than_crown_length(:) = (height > (scorch_height + crown_length))
WHERE(are_flames_taller_than_tree)
prop_ck = 1._wp ! yes
ELSEWHERE(are_flames_shorter_than_crown_length)
prop_ck = 0._wp ! yes
ELSEWHERE
prop_ck = (scorch_height - height + crown_length) / crown_length ! no (Eq. 17 from Thonicke et al. (2010))
ENDWHERE
END FUNCTION calc_prop_ck
Finally, what happens when CL is 0? As already mentioned, that would trigger an exception as it would be divided by 0. Nonetheless, I would leave outside the scope of this function. The reason for that, preventing technical errors is not the main point of this function. On top of that, It is actually interesting an error is raised when it happens, as it would let us know there is a consistency problem within our program. In consequence, somewhere else there would be a proper check on this.