Merge branch 'main' into 368-merge-cablegithub-with-esm15-cable-for-e…

…sm16

documentation/docs/user_guide/inputs/cable_nml.md

@@ -134,7 +134,7 @@ applications. The following are annotated examples of cable.nml:
 | cable_user%access13roots         | logical            | .TRUE. .FALSE.                                             | .FALSE.                          | Switch to use ACCESS1.3 %froot.                                                                         |
 | cable_user%l_limit_labile        | logical            | .TRUE. .FALSE.                                             | .FALSE.                          | Limit labile in spinup.                                                                                 |
 | cable_user%NtilesThruMetFile     | logical            | .TRUE. .FALSE.                                             | .FALSE.                          | Specify Ntiles through met file.                                                                        |
-
+| cable_user%l_ice_consistency     | logical            | .TRUE. .FALSE.                                             | .FALSE.                          | If true, ensures consistency between soil and vegetation tiles with permanent ice. All tiles with permanent ice for soil will have ice for vegetation and vice-versa. All the parameters for these new ice tiles are updated to the ice parameters.                       |
 
 ## For offline applications ##
 

src/offline/cable.nml

@@ -68,4 +68,5 @@
    cable_user%CASA_DUMP_READ = .FALSE.      ! TRUE reads CASA forcing from netcdf format
    cable_user%CASA_DUMP_WRITE = .FALSE.      ! TRUE outputs CASA forcing in netcdf format
    cable_user%SSNOW_POTEV= 'HDM'      ! Humidity Deficit Method
+   cable_user%l_ice_consistency = .FALSE.  ! TRUE add QA step for ice tiles input.
 &end

src/offline/cable_input.F90

@@ -2818,15 +2818,20 @@ SUBROUTINE load_parameters(met,air,ssnow,veg,climate,bgc,soil,canopy,rough,rad,
     END IF
     WRITE(logn,*)
 
-    !| 9. Construct derived parameters and zero initialisations for the
+    !> 9. Ensure the consistency of ice points between soil and vegetation 
+    IF (cable_user%l_ice_consistency) THEN
+      CALL consistency_ice_veg_soil(soil, veg)
+    END IF
+
+    !| 10. Construct derived parameters and zero initialisations for the
     ! groundwater routine, regardless of where parameters and other
     ! initialisations have loaded from
     CALL derived_parameters(soil,sum_flux,bal,ssnow,veg,rough)
 
-    !> 10. Check for basic inconsistencies in parameter values
+    !> 11. Check for basic inconsistencies in parameter values
     CALL check_parameter_values(soil,veg,ssnow)
 
-    !> 11. Write per-site parameter values to log file if requested
+    !> 12. Write per-site parameter values to log file if requested
     CALL report_parameters(logn,soil,veg,bgc,rough,ssnow,canopy, &
          casamet,casapool,casaflux,phen,vegparmnew,verbose)
 

src/offline/cable_parameters.F90

@@ -67,7 +67,7 @@ MODULE cable_param_module
   PRIVATE
   PUBLIC get_default_params, write_default_params, derived_parameters,         &
        check_parameter_values, report_parameters, parID_type,                &
-       write_cnp_params
+       write_cnp_params, consistency_ice_veg_soil 
   INTEGER :: patches_in_parfile=4 ! # patches in default global parameter
   ! file
 
@@ -2064,7 +2064,195 @@ SUBROUTINE check_parameter_values(soil, veg, ssnow)
     END WHERE
 
   END SUBROUTINE check_parameter_values
+
+  !===============================================================================
+
+  SUBROUTINE consistency_ice_veg_soil(soil, veg)  
+    ! Ensure that when an active patch has a veg type of ice then its soil type is also ice and vice versa
+    ! Any change effected to enforce this consistency includes correcting the appropriate paramter values 
+
+    USE grid_constants_mod_cbl, ONLY : ICE_SoilType, ICE_VegType
+    USE cable_phys_constants_mod, ONLY : csice, density_ice
+
+    TYPE (soil_parameter_type), INTENT(INOUT) :: soil  ! soil parameter data
+    TYPE (veg_parameter_type),  INTENT(INOUT) :: veg   ! vegetation parameter
+
+    INTEGER :: i, j, k, kIVP
+    LOGICAL :: LIceVegPatch
+
+
+    LIceVegPatch = .FALSE.
+
+    look_for_ice_veg_patch:     DO i = 1, mland                 ! loop over all gridcells
+       DO j = 1, landpt(i)%nap      ! loop over all active patches in every gridcell
+
+          k = landpt(i)%cstart + j - 1    ! absolute position of active veg patch 
+
+          IF(veg%iveg(k) == ICE_VegType) THEN  ! check to see if there is at least one ice veg patch already
+             LIceVegPatch = .TRUE.
+             kIVP = k    ! remember which patch it is 
+             EXIT look_for_ice_veg_patch ! exit as we only need to find one ice veg patch 
+          END IF
+
+       END DO
+    END DO look_for_ice_veg_patch
+
+
+    DO i = 1, mland                 ! loop over all gridcells 
+       DO j = 1, landpt(i)%nap      ! loop over all active patches in every gridcell 
+
+          k = landpt(i)%cstart + j - 1    ! absolute position of active patch 
+
+          IF(soil%isoilm(k) == ICE_SoilType) THEN  ! check to see if ice soil patch
+
+             WRITE(logn,*) 'SUBROUTINE load_parameters:'
+             WRITE(logn,*) 'At land point number ', i
+             WRITE(logn,*) 'And patch number ', k
+             WRITE(logn,*) 'isoilm is ICE_SoilType'
+             WRITE(logn,*) 'Set rhosoil = density_ice from CABLE physical constants'
+             WRITE(logn,*) 'Set css     = csice from CABLE physical constants'
+
+             soil%rhosoil(k) = density_ice   
+             soil%css(k) =     csice         
+
+          END IF
+
+       END DO
+    END DO
+
+
+    DO i = 1, mland                 ! loop over all gridcells 
+       DO j = 1, landpt(i)%nap      ! loop over all active patches in every gridcell 
+
+          k = landpt(i)%cstart + j - 1    ! absolute position of active patch 
+
+          IF((veg%iveg(k) == ICE_VegType) .AND. (soil%isoilm(k) /= ICE_SoilType)) THEN
+
+
+             WRITE(logn,*) 'SUBROUTINE load_parameters:'
+             WRITE(logn,*) 'At land point number ', i
+             WRITE(logn,*) 'And patch number ', k
+             WRITE(logn,*) 'iveg is ICE_VegType but isoilm is not ICE_SoilType'
+             WRITE(logn,*) 'Changed isoilm to ICE_SoilType with appropriate parameter corrections'
+
+             soil%isoilm(k) = ICE_SoilType
+
+      ! correct appropriately parameters and derived parameters 
+
+             soil%rhosoil(k) = density_ice   
+             soil%css(k) =     csice         
+
+
+          ELSE IF ((veg%iveg(k) /= ICE_VegType) .AND. (soil%isoilm(k) == ICE_SoilType)) THEN
+
+             WRITE(logn,*) 'SUBROUTINE load_parameters:'
+             WRITE(logn,*) 'At land point number ', i
+             WRITE(logn,*) 'And patch number ', k
+             WRITE(logn,*) 'iveg is not ICE_VegType but isoilm is ICE_SoilType'
+             WRITE(logn,*) 'Changed iveg to ICE_VegType with appropriate parameter corrections'
+
+
+             veg%iveg(k) = ICE_VegType
+
+             ! correct appropriately parameters and derived parameters 
+
+             IF (LIceVegPatch) THEN   ! if there is at least one ice veg patch already
+                                      ! use the parameter values from this ice veg patch for the new ice veg patch
+
+                veg%frac4(k)    =        veg%frac4(kIVP)   
+                veg%taul(k,1)   =        veg%taul(kIVP,1)  
+                veg%taul(k,2)   =        veg%taul(kIVP,2)  
+                veg%refl(k,1)   =        veg%refl(kIVP,1)  
+                veg%refl(k,2)   =        veg%refl(kIVP,2)  
+                veg%canst1(k)   =        veg%canst1(kIVP)  
+                veg%dleaf(k)    =        veg%dleaf(kIVP)   
+                veg%vcmax(k)    =        veg%vcmax(kIVP)   
+                veg%ejmax(k)    =        veg%ejmax(kIVP)   
+                veg%hc(k)       =        veg%hc(kIVP)      
+                veg%xfang(k)    =        veg%xfang(kIVP)   
+                veg%vbeta(k)    =        veg%vbeta(kIVP)   
+                veg%xalbnir(k)  =        veg%xalbnir(kIVP) 
+                veg%rp20(k)     =        veg%rp20(kIVP)    
+                veg%rpcoef(k)   =        veg%rpcoef(kIVP)  
+                veg%rs20(k)     =        veg%rs20(kIVP)    
+                veg%shelrb(k)   =        veg%shelrb(kIVP)  
+                veg%wai(k)      =        veg%wai(kIVP)     
+                veg%a1gs(k)     =        veg%a1gs(kIVP)    
+                veg%d0gs(k)     =        veg%d0gs(kIVP)    
+                veg%vegcf(k)    =        veg%vegcf(kIVP)   
+                veg%extkn(k)    =        veg%extkn(kIVP)   
+                veg%tminvj(k)   =        veg%tminvj(kIVP)  
+                veg%tmaxvj(k)   =        veg%tmaxvj(kIVP)  
+                veg%g0(k)       =        veg%g0(kIVP)      
+                veg%g1(k)       =        veg%g1(kIVP)      
+                veg%a1gs(k)     =        veg%a1gs(kIVP)    
+                veg%d0gs(k)     =        veg%d0gs(kIVP)    
+                veg%alpha(k)    =        veg%alpha(kIVP)   
+                veg%convex(k)   =        veg%convex(kIVP)  
+                veg%cfrd(k)     =        veg%cfrd(kIVP)    
+                veg%gswmin(k)   =        veg%gswmin(kIVP)  
+                veg%conkc0(k)   =        veg%conkc0(kIVP)  
+                veg%conko0(k)   =        veg%conko0(kIVP)  
+                veg%ekc(k)      =        veg%ekc(kIVP)     
+                veg%eko(k)      =        veg%eko(kIVP)     
+                veg%rootbeta(k) =        veg%rootbeta(kIVP)
+                veg%zr(k)       =        veg%zr(kIVP)      
+                veg%clitt(k)    =        veg%clitt(kIVP)   
+
+             ELSE           ! otherwise use default parameter values for ice 
+
+                veg%frac4(k)    = vegin%frac4(ICE_VegType)
+                veg%taul(k,1)   = vegin%taul1(ICE_VegType)
+                veg%taul(k,2)   = vegin%taul2(ICE_VegType)
+                veg%refl(k,1)   = vegin%refl1(ICE_VegType)
+                veg%refl(k,2)   = vegin%refl2(ICE_VegType)
+                veg%canst1(k)   = vegin%canst1(ICE_VegType)
+                veg%dleaf(k)    = vegin%dleaf(ICE_VegType)
+                veg%vcmax(k)    = vegin%vcmax(ICE_VegType)
+                veg%ejmax(k)    = vegin%ejmax(ICE_VegType)
+                veg%hc(k)       = vegin%hc(ICE_VegType)
+                veg%xfang(k)    = vegin%xfang(ICE_VegType)
+                veg%vbeta(k)    = vegin%vbeta(ICE_VegType)
+                veg%xalbnir(k)  = vegin%xalbnir(ICE_VegType)
+                veg%rp20(k)     = vegin%rp20(ICE_VegType)
+                veg%rpcoef(k)   = vegin%rpcoef(ICE_VegType)
+                veg%rs20(k)     = vegin%rs20(ICE_VegType)
+                veg%shelrb(k)   = vegin%shelrb(ICE_VegType)
+                veg%wai(k)      = vegin%wai(ICE_VegType)
+                veg%a1gs(k)     = vegin%a1gs(ICE_VegType)
+                veg%d0gs(k)     = vegin%d0gs(ICE_VegType)
+                veg%vegcf(k)    = vegin%vegcf(ICE_VegType)
+                veg%extkn(k)    = vegin%extkn(ICE_VegType)
+                veg%tminvj(k)   = vegin%tminvj(ICE_VegType)
+                veg%tmaxvj(k)   = vegin%tmaxvj(ICE_VegType)
+                veg%g0(k)       = vegin%g0(ICE_VegType) 
+                veg%g1(k)       = vegin%g1(ICE_VegType) 
+                veg%a1gs(k)     = vegin%a1gs(ICE_VegType)
+                veg%d0gs(k)     = vegin%d0gs(ICE_VegType)
+                veg%alpha(k)    = vegin%alpha(ICE_VegType)
+                veg%convex(k)   = vegin%convex(ICE_VegType)
+                veg%cfrd(k)     = vegin%cfrd(ICE_VegType)
+                veg%gswmin(k)   = vegin%gswmin(ICE_VegType)
+                veg%conkc0(k)   = vegin%conkc0(ICE_VegType)
+                veg%conko0(k)   = vegin%conko0(ICE_VegType)
+                veg%ekc(k)      = vegin%ekc(ICE_VegType)
+                veg%eko(k)      = vegin%eko(ICE_VegType)
+                veg%rootbeta(k) = vegin%rootbeta(ICE_VegType)
+                veg%zr(k)       = vegin%zr(ICE_VegType)
+                veg%clitt(k)    = vegin%clitt(ICE_VegType)
+
+             END IF
+
+          END IF
+
+       END DO
+    END DO
+
+
+  END SUBROUTINE consistency_ice_veg_soil
+
   !===============================================================================
+
   SUBROUTINE report_parameters(logn, soil, veg, bgc, rough,                    &
        ssnow, canopy, casamet, casapool, casaflux,     &
        phen, vegparmnew, verbose )

src/params/grid_constants_cbl.F90

@@ -47,6 +47,8 @@ MODULE grid_constants_mod_cbl
 INTEGER, PARAMETER :: ICE_SoilType = 9 ! SoilType Index (soilparm_cable.nml JAC)
 INTEGER, PARAMETER :: lakes_cable  = 16! SoilType Index (soilparm_cable.nml JAC)
 
+INTEGER, PARAMETER :: ICE_VegType = 17  ! permanent ice index for veg
+
 INTEGER, PARAMETER :: mf = 2          ! # leaves (sunlit, shaded)
 INTEGER, PARAMETER :: niter = 4       ! number of iterations for za/L
 

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)