consistency_ice_veg_soil

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_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/util/cable_runtime_opts_mod.F90

@@ -113,7 +113,9 @@ MODULE cable_runtime_opts_mod
   LOGICAL :: access13roots = .FALSE.     !switch to use ACCESS1.3 %froot
 
   LOGICAL :: l_limit_labile = .FALSE.    ! #237: limit Labile in spinup
-  LOGICAL :: NtilesThruMetFile = .FALSE. ! #199: Specify Ntiles thru met file 
+  LOGICAL :: NtilesThruMetFile = .FALSE. ! #199: Specify Ntiles thru met file
+  ! #338 https://github.com/CABLE-LSM/CABLE/issues/338
+  LOGICAL :: l_ice_consistency = .FALSE.
 
 END TYPE kbl_user_switches