379 updated branch with all GW changes

CMakeLists.txt

@@ -146,6 +146,15 @@ if(CABLE_LIBRARY)
         src/util/cable_climate_type_mod.F90
         src/util/masks_cbl.F90
     )
+    target_compile_definitions(cable PRIVATE UM_CBL)
+    target_include_directories(cable PUBLIC "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>")
+
+    install(TARGETS cable
+      EXPORT CABLEFortranTargets
+      LIBRARY DESTINATION lib
+      INCLUDES DESTINATION include)
+
+    install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}/ TYPE INCLUDE)
 
     if(CABLE_LIBRARY_TARGET STREQUAL "access-esm1.6")
         list(APPEND SOURCES

src/offline/cable_checks.F90

@@ -84,7 +84,8 @@ MODULE cable_checks_module
       ESoil = [-0.0015, 0.0015], &
       TVeg = [-0.0003, 0.0003], &
       ECanop = [-0.0003, 0.0003], &
-      PotEvap = [-0.005, 0.005], &  !note should encompass Evap
+      PotEvap = [-0.005, 0.005], &  !note should encompass Evap  ! I have PotEvap = [-0.0006, 0.0006] if it makes a difference - rk4417
+
       ACond = [0.0, 1.0], &
       SoilWet = [-0.4, 1.2], &
       Albedo = [0.0, 1.0], &
@@ -133,12 +134,13 @@ MODULE cable_checks_module
       css = [700.0, 2200.0], &
       rhosoil = [300.0, 3000.0], &
       hyds = [5.0E-7, 8.5E-3], & ! vh_js ! sep14
+! MMY 8.5E-3->8.5 since hyds uses m/s, but hyds_vec uses mm/s
       rs20 = [0.0, 10.0], &
       sand = [0.0, 1.0], &
       sfc = [0.1, 0.5], &
       silt = [0.0, 1.0], &
       ssat = [0.35, 0.5], &
-      sucs = [-0.8, -0.03], &
+      sucs = [-0.8, -0.03], &   ! MMY@23Apr2023 works for CABLE non-GW
       swilt = [0.05, 0.4], &
       froot = [0.0, 1.0], &
       zse = [0.0, 5.0], &
@@ -498,7 +500,10 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
     ! Local variables
     REAL(r_2), DIMENSION(:),POINTER, SAVE     :: owb         ! volumetric soil moisture 
                                                              ! at previous time step
+
     REAL(r_2), DIMENSION(mp)                  :: delwb       ! change in soilmoisture
+   ! To remove qrecharge or not from the mass balance
+    REAL(r_2), DIMENSION(mp)                  :: qrecharge_opt
     ! b/w tsteps
     REAL, DIMENSION(mp)                       :: canopy_wbal !canopy water balance
     INTEGER                                   :: j, k        ! do loop counter
@@ -521,9 +526,14 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
     !      it's included in change in canopy storage calculation))
     ! rml 28/2/11 ! BP changed rnof1+rnof2 to ssnow%runoff which also included rnof5
     ! which is used when nglacier=2 in soilsnow routines (BP feb2011)
+    qrecharge_opt = ssnow%qrecharge
+    IF ( cable_user%GW_MODEL )  THEN
+       qrecharge_opt = 0.0_r_2
+    ENDIF
+
     bal%wbal = REAL(met%precip - canopy%delwc - ssnow%snowd+ssnow%osnowd        &
          - ssnow%runoff-(canopy%fevw+canopy%fevc                                &
-         + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - ssnow%qrecharge)
+         + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - qrecharge_opt)
 
     ! Canopy water balance: precip-change.can.storage-throughfall-evap+dew
     canopy_wbal = REAL(met%precip-canopy%delwc-canopy%through                   &

src/offline/cable_define_types.F90

@@ -217,6 +217,9 @@ MODULE cable_def_types_mod
   TYPE soil_snow_type
 
      INTEGER, DIMENSION(:), POINTER :: isflag ! 0 => no snow 1 => snow
+          !* Flag (0/1) to indicate if snow is present 
+          ! For isflag=1, if there is enough snow, it will use a 3 layer snowpack
+          ! If there isn't enough snow, it will use a single layer snowpack
 
      REAL, DIMENSION(:), POINTER ::                                           &
           iantrct, & ! pointer to Antarctic land points

src/offline/cable_initialise.F90

@@ -141,6 +141,21 @@ SUBROUTINE get_default_inits(met,soil,ssnow,canopy,logn, EMSOIL)
     canopy%fhs     = 0.0   ! sensible heat flux from soil (W/m2)
     canopy%us = 0.1 ! friction velocity (needed in roughness before first call to canopy: should in be in restart?)
 
+! block below added by rk4417 - phase2
+! ccc - These need to be set elsewhere since this subroutine is deprecated
+! (it is never called) and should be removed.
+    canopy%sublayer_dz = 0.001  
+    ssnow%rtevap_sat   = 0.0
+    ssnow%rtevap_unsat = 0.0
+    ssnow%satfrac    = 1.0e-12
+    ssnow%qhz        = 0.0
+    ssnow%wtd        = 1000.0
+    ssnow%wb_hys  = 0.99*soil%ssat_vec
+    ssnow%smp_hys = -0.99*soil%sucs_vec
+!ccc To come back to it    ssnow%ssat_hys = gw_params%ssat_wet_factor*soil%ssat_vec
+    ssnow%watr_hys = soil%watr
+    ssnow%hys_fac = 1.0
+
   END SUBROUTINE get_default_inits
 
   !==============================================================================
@@ -408,20 +423,79 @@ SUBROUTINE get_restart_data(logn,ssnow,canopy,rough,bgc,                       &
     CALL readpar(ncid_rin,'runoff',dummy,ssnow%runoff,filename%restart_in,      &
          max_vegpatches,'def',from_restart,mp)
 
-    !MD
-    ok = NF90_INQ_VARID(ncid_rin,'GWwb',parID)
-    IF(ok == NF90_NOERR) THEN
-       CALL readpar(ncid_rin,'GWwb',dummy,ssnow%GWwb,filename%restart_in,            &
-            max_vegpatches,'def',from_restart,mp)
-    ELSE
-       ssnow%GWwb = 0.95*soil%ssat
+    IF (cable_user%gw_model) THEN   ! inserted by rk4417 - phase2
+      ok = NF90_INQ_VARID(ncid_rin,'GWwb',parID)
+      IF(ok == NF90_NOERR) THEN
+          CALL readpar(ncid_rin,'GWwb',dummy,ssnow%GWwb,filename%restart_in,   &
+               max_vegpatches,'def',from_restart,mp)
+      ELSE
+          ssnow%GWwb = 0.95*soil%ssat
+      END IF
+
+! below part until END IF inserted by rk4417 - phase2
+
+      ok = NF90_INQ_VARID(ncid_rin,'wb_hys',parID)
+      IF(ok == NF90_NOERR) THEN 
+          CALL readpar(ncid_rin,'wb_hys',dummy,ssnow%wb_hys,filename%restart_in,&
+               max_vegpatches,'msd',from_restart,mp)   
+      ELSE
+         ssnow%wb_hys = 0.99*soil%ssat_vec
+      END IF
+
+      ok = NF90_INQ_VARID(ncid_rin,'smp_hys',parID)
+      IF(ok == NF90_NOERR) THEN 
+          CALL readpar(ncid_rin,'smp_hys',dummy,ssnow%smp_hys,filename%restart_in,&
+               max_vegpatches,'msd',from_restart,mp)   
+      ELSE
+         ssnow%smp_hys = -1.0*abs(soil%sucs_vec)*0.99
+      END IF
+
+      ok = NF90_INQ_VARID(ncid_rin,'ssat_hys',parID)
+      IF(ok == NF90_NOERR) THEN 
+          CALL readpar(ncid_rin,'ssat_hys',dummy,ssnow%ssat_hys,filename%restart_in,&
+               max_vegpatches,'msd',from_restart,mp)   
+      ELSE
+          ssnow%ssat_hys = soil%ssat_vec
+      END IF
+
+      ok = NF90_INQ_VARID(ncid_rin,'watr_hys',parID)
+      IF(ok == NF90_NOERR) THEN 
+          CALL readpar(ncid_rin,'watr_hys',dummy,ssnow%watr_hys,filename%restart_in,&
+               max_vegpatches,'msd',from_restart,mp)   
+      ELSE
+          ssnow%watr_hys = soil%watr
+      END IF
+
+
+      ok = NF90_INQ_VARID(ncid_rin,'hys_fac',parID)
+      IF(ok == NF90_NOERR) THEN 
+          CALL readpar(ncid_rin,'hys_fac',dummy,ssnow%hys_fac,filename%restart_in,&
+               max_vegpatches,'msd',from_restart,mp)   
+      ELSE
+          ssnow%hys_fac = 1.0
+      END IF
     END IF
 
 !   IF(cable_user%SOIL_STRUC=='sli'.or.cable_user%FWSOIL_SWITCH=='Haverd2013') THEN
 !      CALL readpar(ncid_rin,'gamma',dummy,veg%gamma,filename%restart_in,           &
 !           max_vegpatches,'def',from_restart,mp)
 !   ENDIF
 
+    IF (cable_user%or_evap) THEN
+      ok = NF90_INQ_VARID(ncid_rin,'sublayer_dz',parID)
+      IF(ok == NF90_NOERR) THEN 
+         CALL readpar(ncid_rin,'sublayer_dz',dummy,canopy%sublayer_dz,filename%restart_in,&
+              max_vegpatches,'def',from_restart,mp)   
+      ELSE
+         canopy%sublayer_dz(:) = 0.01
+      END IF
+
+      IF (ANY(canopy%sublayer_dz < 0.0) .OR. ANY(canopy%sublayer_dz > 0.5)) THEN
+         WRITE(*,*) 'problem with sublayer_dz and restart.  check restart values!'
+      END IF
+    END IF
+
+
     IF(cable_user%SOIL_STRUC=='sli') THEN
        CALL readpar(ncid_rin,'S',dummy,ssnow%S,filename%restart_in, &
             max_vegpatches,'ms',from_restart,mp)

src/offline/cable_input.F90

@@ -2866,6 +2866,11 @@ SUBROUTINE load_parameters(met,air,ssnow,veg,climate,bgc,soil,canopy,rough,rad,
          sum_flux,veg,mp)
     WRITE(logn,*) ' CABLE variables allocated with ', mp, ' patch(es).'
 
+! call below inserted by rk4417 - phase2
+    !CALL for gw_model false and true sets constants when false
+    CALL GWspatialParameters(logn,soil,ssnow)    ! MMY gw_model = True read var from gridinfo
+    ! MMY gw_model = False use default values
+
     IF (icycle > 0 .OR. CABLE_USER%CASA_DUMP_WRITE ) &
          CALL alloc_casavariable(casabiome,casapool,casaflux, &
          casamet,casabal,mp)

src/offline/cable_iovars.F90

@@ -182,7 +182,9 @@ MODULE cable_IO_vars_module
           veg_class,soil_class,mvtype,mstype,patchfrac,                        &
           WatSat,GWWatSat,SoilMatPotSat,GWSoilMatPotSat,                       &
           HkSat,GWHkSat,FrcSand,FrcClay,Clappb,Watr,GWWatr,sfc_vec,forg,swilt_vec, &
-          slope,slope_std,GWdz,SatFracmax,Qhmax,QhmaxEfold,HKefold,HKdepth
+          slope,slope_std,GWdz,SatFracmax,Qhmax,QhmaxEfold,HKefold,HKdepth,   &
+          sand_vec,clay_vec,bch_vec,org_vec,elev,elev_std  ! inserted by rk4417 - phase2
+
      INTEGER :: ishorizon,nhorizons,clitt, &
           zeta,fsatmax, &
           gamma,ZR,F10
@@ -404,14 +406,23 @@ MODULE cable_IO_vars_module
           isoil  = .FALSE.,    & ! soil type from global index
           meth  = .FALSE.,     & ! method for solving turbulence in canopy scheme
           za  = .FALSE.,       & ! something to do with roughness ????
+          elev = .false.,&       !mean subgrid elev    ! inserted 2 lines - rk4417 - phase2
+          elev_std=.false.,&     !stddev of subgrid elev
           slope = .FALSE.,&      !mean subgrid slope
           slope_std=.FALSE.,&    !stddev of subgrid slope
           GWdz=.FALSE.,&         !aquifer thickness
           SatFracmax=.FALSE.,&
           Qhmax=.FALSE.,&
           QhmaxEfold=.FALSE.,&
           HKefold=.FALSE.,&
-          HKdepth
+!          HKdepth ! commented out by rk4417 - phase2
+          HKdepth=.false.,&  ! added this block - rk4417 - phase2
+          SMP=.false.,&
+          SMP_hys=.false.,&
+          WB_hys=.false.,&
+          SSAT_hys=.false.,&
+          WATR_hys=.false.,&
+          hys_fac=.false.
 
   END TYPE output_inclusion_type
 

src/offline/cable_mpicommon.F90

@@ -29,7 +29,10 @@ MODULE cable_mpicommon
 
   ! base number of input fields: must correspond to CALLS to 
   ! MPI_address (field ) in *_mpimaster/ *_mpiworker
-  INTEGER, PARAMETER :: nparam = 347
+! Values from rebase in main. Do we need to update the value from 351 considering this one is higher than expected? ccc
+!  INTEGER, PARAMETER :: nparam = 347
+!  INTEGER, PARAMETER :: nparam = 341 ! replaced by below by ccc - GW
+  INTEGER, PARAMETER :: nparam =351  
 
   ! MPI: extra params sent only if nsoilparmnew is true
   INTEGER, PARAMETER :: nsoilnew = 1
@@ -77,7 +80,8 @@ MODULE cable_mpicommon
   !INTEGER, PARAMETER :: nmat = 29
   ! MPI: CABLE_r491, after following up with Bernard on the new variables
   ! vh sli nmat + 4 36 -> 40
-  INTEGER, PARAMETER :: nmat = 39
+!  INTEGER, PARAMETER :: nmat = 40 ! replaced by below by rk4417 - phase2
+  INTEGER, PARAMETER :: nmat = 45  !hysteresis 41
 
   ! MPI: number of contig vector parts / worker (results)
   !INTEGER, PARAMETER :: nvec = 149
@@ -94,7 +98,11 @@ MODULE cable_mpicommon
   ! vh sli nvec + 6 162 -> 168
   ! INTEGER, PARAMETER :: nvec = 172! 168
   ! INH REV_CORR +3  (SSEB +2 will be needed)
+
+! Keeping the highest value here but might need adjustment. ccc
   INTEGER, PARAMETER :: nvec = 183
+!  INTEGER, PARAMETER :: nvec = 175  ! replaced by below by rk4417 - phase2
+!  INTEGER, PARAMETER :: nvec = 176! 176!175
 
   ! MPI: number of final casa result matrices and vectors to receive
   ! by the master for casa_poolout and casa_fluxout