Fix: uninitialised variables in cable_canopy.F90 (#352)

# CABLE

Thank you for submitting a pull request to the CABLE Project.

## Description

Fixes #351. 

Move the definition of sum_rad_gradis and sum_rad_rniso out of the NITER
loop in `define_canopy()` since these variables do not change in the
loop. Having a quick comparison between the outputs before and after
this commit [the effect is
negligeable](#352 (comment)).

Use `sum_rad_gradis` and `sum_rad_rniso` as calculated in
define_canopy() everywhere in define_canopy(), dryLeaf() and wetLeaf(),
instead of recalculating them sometimes.

Move the initialisation of `canopy%DvLitt` and `canopy%kthLitt` to the
top of `define_canopy()` since these variables are constant throughout
the run. We also want these variables to be defined in all cases to
avoid errors in function calls, but only used with `cable_user%litter`
turned on.

## Type of change

Please delete options that are not relevant.

- [X] Bug fix

## Testing

benchcab/me.org analysis:
https://modelevaluation.org/analyses/anywhere/BPLvifw3DyARokYdg/s6k22L3WajmiS9uGv/CJGXP5GQWhGf3nH28/all

Tested in benchcab the changes for DvLitt and kthLitt on their own.
benchcab showed bitwise identical results when moving the initialisation
compared to the previous implementation with the litter option on.

Please add a reviewer when ready for review.


<!-- readthedocs-preview cable start -->
----
📚 Documentation preview 📚:
https://cable--352.org.readthedocs.build/en/352/

<!-- readthedocs-preview cable end -->

src/science/canopy/cable_canopy.F90

@@ -186,6 +186,16 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
     canopy%tv = met%tvair
     canopy%fwsoil = 1.0
 
+    ! Initialise canopy%DvLitt and canopy%kthLitt. This value is only used if
+    ! cable_user%litter is .TRUE.
+    ! Reference:
+    ! Matthews (2006), A process-based model of fine fuel moisture,
+    !                 International Journal of Wildland Fire 15,155-168
+    !                 https://doi.org/10.1071/WF05063
+    ! assuming here u=1.0 ms-1, bulk litter density 63.5 kgm-3
+    canopy%kthLitt = 0.3_r_2 ! ~ 0.2992125984251969 = 0.2+0.14*0.045*1000.0/63.5
+    canopy%DvLitt = 3.1415841138194147e-05_r_2 ! = 2.17e-5*exp(1.0*2.6)*exp(-0.5*(2.08+(1.0*2.38)))
+
     CALL define_air (met, air)
 
     CALL qsatfjh(mp, qstvair, CRMH2o, Crmair, CTETENA, CTETENB, CTETENC, met%tvair-CTfrz,met%pmb)
@@ -234,6 +244,9 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
     canopy%zetash(:,1) = CZETA0 ! stability correction terms
     canopy%zetash(:,2) = CZETPOS + 1
 
+    sum_rad_rniso = SUM(rad%rniso,2)
+    sum_rad_gradis = SUM(rad%gradis,2)
+
 
     DO iter = 1, NITER
 
@@ -312,14 +325,6 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
        ELSE ! NOT sli
           rt0 = MAX(rt_min,rough%rt0us / canopy%us)
 
-          IF (cable_user%litter) THEN
-             ! Mathews (2006), A process-based model of offine fuel moisture,
-             !                 International Journal of Wildland Fire 15,155-168
-             ! assuming here u=1.0 ms-1, bulk litter density 63.5 kgm-3
-             canopy%kthLitt = 0.3_r_2 ! ~ 0.2992125984251969 = 0.2+0.14*0.045*1000.0/63.5
-             canopy%DvLitt = 3.1415841138194147e-05_r_2 ! = 2.17e-5*exp(1.0*2.6)*exp(-0.5*(2.08+(1.0*2.38)))
-          ENDIF
-
        ENDIF
 
 !       ! Aerodynamic resistance (sum 3 height integrals)/us
@@ -385,17 +390,15 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
        ENDDO
 
 
-       rny = SUM(rad%rniso,2) ! init current estimate net rad
+       rny = sum_rad_rniso ! init current estimate net rad
        hcy = 0.0              ! init current estimate lat heat
        ecy = rny - hcy        ! init current estimate lat heat
 
-       sum_rad_rniso = SUM(rad%rniso,2)
-
        CALL dryLeaf( dels, rad, rough, air, met,                             &
                   veg, canopy, soil, ssnow, dsx,                             &
                   fwsoil, tlfx, tlfy, ecy, hcy,                              &
                   rny, gbhu, gbhf, csx, cansat,                              &
-            ghwet,  iter,climate )
+            ghwet,  iter,climate, sum_rad_gradis, sum_rad_rniso )
 
 
       CALL wetLeaf( dels,                                 &
@@ -418,8 +421,6 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
        canopy%fnv = REAL(ftemp)
 
        ! canopy rad. temperature calc from long-wave rad. balance
-       sum_rad_gradis = SUM(rad%gradis,2)
-
        DO j=1,mp
 
           IF ( canopy%vlaiw(j) > CLAI_THRESH .AND.                             &
@@ -649,7 +650,7 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
 
 
 
-       canopy%rniso = SUM(rad%rniso,2) + rad%qssabs + rad%transd*met%fld + &
+       canopy%rniso = sum_rad_rniso + rad%qssabs + rad%transd*met%fld + &
             (1.0-rad%transd)*CEMLEAF* &
             CSBOLTZ*met%tvrad**4 - CEMSOIL*CSBOLTZ*met%tvrad**4
 
@@ -1008,11 +1009,11 @@ SUBROUTINE define_canopy(bal,rad,rough,air,met,dels,ssnow,soil,veg, canopy,clima
 ssnow%dfe_dtg = ssnow%dfe_ddq * ssnow%ddq_dtg
 canopy%dgdtg = ssnow%dfn_dtg - ssnow%dfh_dtg - ssnow%dfe_dtg
 
-bal%drybal = REAL(ecy+hcy) - SUM(rad%rniso,2)                               &
-     + CCAPP*Crmair*(tlfy-met%tk)*SUM(rad%gradis,2)  ! YP nov2009
+bal%drybal = REAL(ecy+hcy) - sum_rad_rniso                               &
+     + CCAPP*Crmair*(tlfy-met%tk)*sum_rad_gradis  ! YP nov2009
 
-bal%wetbal = canopy%fevw + canopy%fhvw - SUM(rad%rniso,2) * canopy%fwet      &
-     + CCAPP*Crmair * (tlfy-met%tk) * SUM(rad%gradis,2) *          &
+bal%wetbal = canopy%fevw + canopy%fhvw - sum_rad_rniso * canopy%fwet      &
+     + CCAPP*Crmair * (tlfy-met%tk) * sum_rad_gradis *          &
      canopy%fwet  ! YP nov2009
 
 DEALLOCATE(cansat,gbhu)

src/science/canopy/cbl_dryLeaf.F90

@@ -11,7 +11,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
        veg, canopy, soil, ssnow, dsx,                             &
        fwsoil, tlfx,  tlfy,  ecy, hcy,                            &
        rny, gbhu, gbhf, csx,                                      &
-       cansat, ghwet, iter,climate )
+       cansat, ghwet, iter,climate, sum_rad_gradis, sum_rad_rniso )
 
     USE cable_def_types_mod
     USE cable_common_module
@@ -71,6 +71,10 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
 
     REAL, INTENT(IN)     :: dels ! integration time step (s)
 
+    ! Sums of radiation quantities over rad bands
+    REAL, INTENT(IN) :: sum_rad_rniso(mp)
+    REAL, INTENT(IN) :: sum_rad_gradis(mp)
+
     !local variables
     REAL, PARAMETER  ::                                                         &
          co2cp3 = 0.0,  & ! CO2 compensation pt C3
@@ -88,8 +92,6 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
          deltlfy,       & ! del temp successive iter.
          gras,          & ! Grashof coeff
          evapfb,        & !
-         sum_rad_rniso, & !
-         sum_rad_gradis,& !
          gwwet,         & ! cond for water for a wet canopy
          ghrwet,        & ! wet canopy cond: heat & thermal rad
          sum_gbh,       & !
@@ -170,13 +172,12 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
     lower_limit2 = rad%scalex * gsw_term
     gswmin = MAX(1.e-6,lower_limit2)
 
-
     gw = 1.0e-3 ! default values of conductance
     gh = 1.0e-3
     ghr= 1.0e-3
     rdx = 0.0
     anx = 0.0
-    rnx = SUM(rad%rniso,2)
+    rnx = sum_rad_rniso
     abs_deltlf = 999.0
 
 
@@ -185,7 +186,7 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
     hcx = 0.0              ! init sens heat iteration memory variable
     hcy = 0.0
     rdy = 0.0
-    ecx = SUM(rad%rniso,2) ! init lat heat iteration memory variable
+    ecx = sum_rad_rniso ! init lat heat iteration memory variable
     tlfxx = tlfx
     psycst(:,:) = SPREAD(air%psyc,2,mf)
     canopy%fevc = 0.0
@@ -197,8 +198,6 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
     canopy%fevw = 0.0
     canopy%fhvw = 0.0
     sum_gbh = SUM((gbhu+gbhf),2)
-    sum_rad_rniso = SUM(rad%rniso,2)
-    sum_rad_gradis = SUM(rad%gradis,2)
 
     DO kk=1,mp
 
@@ -528,17 +527,17 @@ SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
 
              ENDIF
              ! Update canopy sensible heat flux:
-             hcx(i) = (SUM(rad%rniso(i,:))-ecx(i)                               &
+             hcx(i) = (sum_rad_rniso(i)-ecx(i)                               &
                   - Ccapp*Crmair*(met%tvair(i)-met%tk(i))                       &
-                  * SUM(rad%gradis(i,:)))                                         &
+                  * sum_rad_gradis(i))                                         &
                   * SUM(gh(i,:))/ SUM(ghr(i,:))
              ! Update leaf temperature:
              tlfx(i)=met%tvair(i)+REAL(hcx(i))/(Ccapp*Crmair*SUM(gh(i,:)))
 
              ! Update net radiation for canopy:
-             rnx(i) = SUM( rad%rniso(i,:)) -                                    &
+             rnx(i) = sum_rad_rniso(i) -                                    &
                   CCAPP * Crmair *( tlfx(i)-met%tk(i) ) *                 &
-                  SUM( rad%gradis(i,:) )
+                  sum_rad_gradis(i)
 
              ! Update leaf surface vapour pressure deficit:
              dsx(i) = met%dva(i) + air%dsatdk(i) * (tlfx(i)-met%tvair(i))