Update from master

.github/workflows/build_CoLM_gnu.yml

@@ -13,7 +13,13 @@ on:
   push:
     branches:
       - master
-
+    paths-ignore:
+      - 'postprocess/**'
+      - 'preprocess/**'
+      - 'run/**'
+      - 'README.md'
+      - '.gitignore'
+      - '**/**.sh'
   workflow_dispatch:
 
 jobs:

.gitignore

@@ -7,3 +7,6 @@ __pycache__
 *.mod
 .vscode
 .bld
+lib
+*.a
+*.x

Makefile

@@ -79,6 +79,7 @@ OBJS_MKSRFDATA = \
 				  Aggregation_Topography.o          \
 				  Aggregation_TopographyFactors.o   \
 				  Aggregation_Urban.o               \
+				  Aggregation_SoilTexture.o         \
 				  MOD_MeshFilter.o                  \
 				  MOD_RegionClip.o                  \
 				  MKSRFDATA.o
@@ -143,6 +144,7 @@ OBJS_BASIC =    \
 				 MOD_DBedrockReadin.o           \
 				 MOD_SoilColorRefl.o            \
 				 MOD_SoilParametersReadin.o     \
+				 MOD_SoilTextureReadin.o        \
 				 MOD_HtopReadin.o               \
 				 MOD_UrbanReadin.o              \
 				 MOD_BGC_CNSummary.o            \

include/define.h

@@ -14,6 +14,14 @@
 
 ! 2.1 3D Urban model (put it temporarily here):
 #undef URBAN_MODEL
+!    Dependence:  only LULC_IGBP subgrid type for
+!    single point URBAN_MODEL right now.
+#if (defined URBAN_MODEL && defined SinglePoint)
+#define LULC_IGBP
+#undef LULC_USGS
+#undef LULC_IGBP_PFT
+#undef LULC_IGBP_PC
+#endif
 
 ! 3. If defined, debug information is output.
 #define CoLMDEBUG
@@ -64,10 +72,10 @@
 #undef DataAssimilation
 
 ! 10. Vector write model.
-!     1) "VectorInOneFileP" : write vector data in one file in parallel mode;  
-!     2) "VectorInOneFileS" : write vector data in one file in serial mode;  
-!     3) Neither "VectorInOneFileS" nor "VectorInOneFileP" is defined : 
-!        write vector data in separate files.  
+!     1) "VectorInOneFileP" : write vector data in one file in parallel mode;
+!     2) "VectorInOneFileS" : write vector data in one file in serial mode;
+!     3) Neither "VectorInOneFileS" nor "VectorInOneFileP" is defined :
+!        write vector data in separate files.
 #undef VectorInOneFileP
 !     Conflict
 #ifdef VectorInOneFileP

main/HYDRO/MOD_Catch_SubsurfaceFlow.F90

@@ -20,13 +20,20 @@ MODULE MOD_Catch_SubsurfaceFlow
    IMPLICIT NONE
 
    real(r8), parameter :: e_ice  = 6.0   ! soil ice impedance factor
-   real(r8), parameter :: raniso = 1.    ! anisotropy ratio, unitless
+
+   ! anisotropy ratio of lateral/vertical hydraulic conductivity (unitless)
+   ! for USDA soil texture class:
+   ! 0: undefined
+   ! 1: clay;  2: silty clay;  3: sandy clay;   4: clay loam;   5: silty clay loam;   6: sandy clay loam; &    
+   ! 7: loam;  8: silty loam;  9: sandy loam;  10: silt;       11: loamy sand;       12: sand
+   real(r8), parameter :: raniso(0:12) = (/ 1., &
+                                            48., 40., 28., 24., 20., 14., 12., 10., 4., 2., 3., 2. /)
 
    ! -- neighbour variables --
    type(pointer_real8_1d), allocatable :: agwt_nb    (:)  ! ground water area (for patchtype <= 2) of neighbours [m^2]
    type(pointer_real8_1d), allocatable :: theta_a_nb (:)  ! saturated volume content [-]
    type(pointer_real8_1d), allocatable :: zwt_nb     (:)  ! water table depth [m]
-   type(pointer_real8_1d), allocatable :: Ks_nb      (:)  ! saturated hydraulic conductivity [m/s]
+   type(pointer_real8_1d), allocatable :: Kl_nb      (:)  ! lateral hydraulic conductivity [m/s]
    type(pointer_real8_1d), allocatable :: wdsrf_nb   (:)  ! depth of surface water [m]
    type(pointer_logic_1d), allocatable :: islake_nb  (:)  ! whether a neighbour is water body
 
@@ -59,7 +66,7 @@ SUBROUTINE basin_neighbour_init ()
          CALL allocate_neighbour_data (agwt_nb   )
          CALL allocate_neighbour_data (theta_a_nb)
          CALL allocate_neighbour_data (zwt_nb    )
-         CALL allocate_neighbour_data (Ks_nb     )
+         CALL allocate_neighbour_data (Kl_nb    )
          CALL allocate_neighbour_data (wdsrf_nb  )
          CALL allocate_neighbour_data (islake_nb )
 
@@ -127,25 +134,25 @@ SUBROUTINE subsurface_flow (deltime)
 
    real(r8), allocatable :: theta_a_h (:) 
    real(r8), allocatable :: zwt_h     (:) 
-   real(r8), allocatable :: Ks_h      (:) ! [m/s]
+   real(r8), allocatable :: Kl_h      (:) ! [m/s]
    real(r8), allocatable :: xsubs_h   (:) ! [m/s]
    real(r8), allocatable :: xsubs_fc  (:) ! [m/s]
 
    logical  :: j_is_river
    real(r8) :: theta_s_h, air_h, icefrac, imped, delp
    real(r8) :: sumwt, sumarea, zwt_mean
    real(r8) :: zsubs_h_up, zsubs_h_dn
-   real(r8) :: slope, bdamp, Ks_fc, Ks_in
+   real(r8) :: slope, bdamp, Kl_fc, Kl_in
    real(r8) :: ca, cb
    real(r8) :: alp
 
    real(r8), allocatable :: theta_a_bsn (:) 
    real(r8), allocatable :: zwt_bsn     (:) 
-   real(r8), allocatable :: Ks_bsn      (:) ! [m/s]
+   real(r8), allocatable :: Kl_bsn      (:) ! [m/s]
 
    integer  :: jnb
-   real(r8) :: zsubs_up, zwt_up, Ks_up, theta_a_up, area_up
-   real(r8) :: zsubs_dn, zwt_dn, Ks_dn, theta_a_dn, area_dn
+   real(r8) :: zsubs_up, zwt_up, Kl_up, theta_a_up, area_up
+   real(r8) :: zsubs_dn, zwt_dn, Kl_dn, theta_a_dn, area_dn
    real(r8) :: lenbdr, xsubs_nb
    logical  :: iam_lake, nb_is_lake, has_river
 
@@ -185,7 +192,7 @@ SUBROUTINE subsurface_flow (deltime)
          IF (numbasin > 0) THEN
             allocate (theta_a_bsn (numbasin));  theta_a_bsn = 0.
             allocate (zwt_bsn     (numbasin));  zwt_bsn     = 0.
-            allocate (Ks_bsn      (numbasin));  Ks_bsn      = 0.
+            allocate (Kl_bsn      (numbasin));  Kl_bsn      = 0.
          ENDIF
 
          DO ibasin = 1, numbasin
@@ -199,7 +206,7 @@ SUBROUTINE subsurface_flow (deltime)
 
             allocate (theta_a_h (nhru));  theta_a_h = 0.
             allocate (zwt_h     (nhru));  zwt_h     = 0.
-            allocate (Ks_h      (nhru));  Ks_h      = 0.
+            allocate (Kl_h      (nhru));  Kl_h      = 0.
 
             DO i = 1, nhru
 
@@ -253,23 +260,23 @@ SUBROUTINE subsurface_flow (deltime)
                      ENDIF
                   ENDIF
 
-                  Ks_h(i) = 0.
+                  Kl_h(i) = 0.
                   sumwt   = 0.
                   DO ipatch = ps, pe
                      IF (patchtype(ipatch) <= 2) THEN
                         DO ilev = 1, nl_soil
                            icefrac = min(1., wice_soisno(ilev,ipatch)/denice/dz_soi(ilev)/porsl(ilev,ipatch))
                            imped   = 10.**(-e_ice*icefrac)
-                           Ks_h(i) = Ks_h(i) + hru_patch%subfrc(ipatch) &
+                           Kl_h(i) = Kl_h(i) + hru_patch%subfrc(ipatch) * raniso(soiltext(ipatch)) &
                               * hksati(ilev,ipatch)/1.0e3 * imped * dz_soi(ilev)/zi_soi(nl_soil) 
                         ENDDO
                         sumwt = sumwt + hru_patch%subfrc(ipatch)
                      ENDIF
                   ENDDO
-                  IF (sumwt > 0) Ks_h(i) = Ks_h(i) / sumwt
+                  IF (sumwt > 0) Kl_h(i) = Kl_h(i) / sumwt
                ELSE
                   ! Frozen soil.
-                  Ks_h(i) = 0.
+                  Kl_h(i) = 0.
                ENDIF
 
             ENDDO
@@ -285,11 +292,11 @@ SUBROUTINE subsurface_flow (deltime)
                j = hrus%inext(i)
 
                IF (j <= 0)        CYCLE ! downstream is out of catchment
-               IF (Ks_h(i) == 0.) CYCLE ! this HRU is frozen
+               IF (Kl_h(i) == 0.) CYCLE ! this HRU is frozen
 
                j_is_river = (hrus%indx(j) == 0)
 
-               IF ((.not. j_is_river) .and. (Ks_h(j) == 0.)) CYCLE ! non-river downstream HRU is frozen
+               IF ((.not. j_is_river) .and. (Kl_h(j) == 0.)) CYCLE ! non-river downstream HRU is frozen
 
                zsubs_h_up = hrus%elva(i) - zwt_h(i)
 
@@ -317,27 +324,27 @@ SUBROUTINE subsurface_flow (deltime)
                IF ((zsubs_h_up > zsubs_h_dn) .or. j_is_river) THEN
                   IF (zwt_h(i) > 1.5) THEN
                      ! from Fan et al., JGR 112(D10125)
-                     Ks_fc = raniso * Ks_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp)
+                     Kl_fc = Kl_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp)
                   ELSE
-                     Ks_fc = raniso * Ks_h(i) * ((1.5-zwt_h(i)) + bdamp)
+                     Kl_fc = Kl_h(i) * ((1.5-zwt_h(i)) + bdamp)
                   ENDIF
                ELSE
                   IF (zwt_h(j) > 1.5) THEN
-                     Ks_fc = raniso * Ks_h(j) * bdamp * exp(-(zwt_h(j)-1.5)/bdamp)
+                     Kl_fc = Kl_h(j) * bdamp * exp(-(zwt_h(j)-1.5)/bdamp)
                   ELSE
-                     Ks_fc = raniso * Ks_h(j) * ((1.5-zwt_h(j)) + bdamp)
+                     Kl_fc = Kl_h(j) * ((1.5-zwt_h(j)) + bdamp)
                   ENDIF
                ENDIF
 
-               ca = hrus%flen(i) * Ks_fc / theta_a_h(i) / delp / hrus%agwt(i) * deltime
+               ca = hrus%flen(i) * Kl_fc / theta_a_h(i) / delp / hrus%agwt(i) * deltime
 
                IF (.not. j_is_river) THEN
-                  cb = hrus%flen(i) * Ks_fc / theta_a_h(j) / delp / hrus%agwt(j) * deltime
+                  cb = hrus%flen(i) * Kl_fc / theta_a_h(j) / delp / hrus%agwt(j) * deltime
                ELSE
-                  cb = hrus%flen(i) * Ks_fc / delp / hrus%area(j) * deltime
+                  cb = hrus%flen(i) * Kl_fc / delp / hrus%area(j) * deltime
                ENDIF
 
-               xsubs_fc(i) = (zsubs_h_up - zsubs_h_dn) * hrus%flen(i) * Ks_fc / (1+ca+cb) / delp
+               xsubs_fc(i) = (zsubs_h_up - zsubs_h_dn) * hrus%flen(i) * Kl_fc / (1+ca+cb) / delp
 
                xsubs_h(i) = xsubs_h(i) + xsubs_fc(i) / hrus%agwt(i)
 
@@ -390,9 +397,9 @@ SUBROUTINE subsurface_flow (deltime)
 
                   IF (zwt_h(i) > 1.5) THEN
                      ! from Fan et al., JGR 112(D10125)
-                     Ks_in = raniso * Ks_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp)
+                     Kl_in = Kl_h(i) * bdamp * exp(-(zwt_h(i)-1.5)/bdamp)
                   ELSE
-                     Ks_in = raniso * Ks_h(i) * ((1.5-zwt_h(i)) + bdamp)
+                     Kl_in = Kl_h(i) * ((1.5-zwt_h(i)) + bdamp)
                   ENDIF
 
                   ps = hru_patch%substt(hrus%ihru(i))
@@ -403,7 +410,7 @@ SUBROUTINE subsurface_flow (deltime)
 
                      DO ipatch = ps, pe
                         IF (patchtype(ipatch) <= 2) THEN
-                           xsubs_pch(ipatch) = - Ks_in * (zwt(ipatch) - zwt_mean) *6.0*pi/hrus%agwt(i)
+                           xsubs_pch(ipatch) = - Kl_in * (zwt(ipatch) - zwt_mean) *6.0*pi/hrus%agwt(i)
                            ! Update total subsurface lateral flow (2): Between patches
                            xwsub(ipatch) = xwsub(ipatch) + xsubs_pch(ipatch) * 1.e3 ! m/s to mm/s
                         ENDIF
@@ -416,20 +423,20 @@ SUBROUTINE subsurface_flow (deltime)
             IF (sumarea > 0) THEN
                theta_a_bsn (ibasin) = sum(theta_a_h * hrus%agwt) / sumarea
                zwt_bsn     (ibasin) = sum(zwt_h     * hrus%agwt) / sumarea
-               Ks_bsn      (ibasin) = sum(Ks_h      * hrus%agwt) / sumarea
+               Kl_bsn      (ibasin) = sum(Kl_h      * hrus%agwt) / sumarea
             ENDIF
 
             deallocate (theta_a_h)
             deallocate (zwt_h    )
-            deallocate (Ks_h     )
+            deallocate (Kl_h     )
             deallocate (xsubs_h  )
             deallocate (xsubs_fc )
 
          ENDDO
 
          CALL retrieve_neighbour_data (theta_a_bsn, theta_a_nb)
          CALL retrieve_neighbour_data (zwt_bsn    , zwt_nb    )
-         CALL retrieve_neighbour_data (Ks_bsn     , Ks_nb     )
+         CALL retrieve_neighbour_data (Kl_bsn     , Kl_nb     )
          CALL retrieve_neighbour_data (wdsrf_bsn  , wdsrf_nb  )
 
          DO ibasin = 1, numbasin
@@ -449,7 +456,7 @@ SUBROUTINE subsurface_flow (deltime)
                ENDIF
 
                IF (.not. iam_lake) THEN
-                  Ks_up      = Ks_bsn     (ibasin)
+                  Kl_up      = Kl_bsn    (ibasin)
                   zwt_up     = zwt_bsn    (ibasin)
                   theta_a_up = theta_a_bsn(ibasin)
                   zsubs_up   = elementneighbour(ibasin)%myelva - zwt_up 
@@ -461,7 +468,7 @@ SUBROUTINE subsurface_flow (deltime)
                ENDIF
 
                IF (.not. nb_is_lake) THEN
-                  Ks_dn      = Ks_nb(ibasin)%val(jnb)
+                  Kl_dn      = Kl_nb(ibasin)%val(jnb)
                   zwt_dn     = zwt_nb(ibasin)%val(jnb)
                   theta_a_dn = theta_a_nb(ibasin)%val(jnb)
                   zsubs_dn   = elementneighbour(ibasin)%elva(jnb) - zwt_dn
@@ -474,8 +481,8 @@ SUBROUTINE subsurface_flow (deltime)
 
                IF ((.not. iam_lake)   .and. (area_up <= 0)) CYCLE
                IF ((.not. nb_is_lake) .and. (area_dn <= 0)) CYCLE
-               IF ((.not. iam_lake)   .and. (Ks_up == 0.) ) CYCLE 
-               IF ((.not. nb_is_lake) .and. (Ks_dn == 0.) ) CYCLE
+               IF ((.not. iam_lake)   .and. (Kl_up == 0. )) CYCLE 
+               IF ((.not. nb_is_lake) .and. (Kl_dn == 0. )) CYCLE
 
                ! water body is dry.
                IF (iam_lake .and. (zsubs_up > zsubs_dn) .and. (wdsrf_bsn(ibasin) == 0.)) THEN
@@ -507,22 +514,22 @@ SUBROUTINE subsurface_flow (deltime)
                IF (nb_is_lake .or. ((.not. iam_lake) .and. (zsubs_up > zsubs_dn))) THEN
                   IF (zwt_up > 1.5) THEN
                      ! from Fan et al., JGR 112(D10125)
-                     Ks_fc = raniso * Ks_up * bdamp * exp(-(zwt_up-1.5)/bdamp)
+                     Kl_fc = Kl_up * bdamp * exp(-(zwt_up-1.5)/bdamp)
                   ELSE
-                     Ks_fc = raniso * Ks_up * ((1.5-zwt_up) + bdamp)
+                     Kl_fc = Kl_up * ((1.5-zwt_up) + bdamp)
                   ENDIF
                ELSE
                   IF (zwt_dn > 1.5) THEN
-                     Ks_fc = raniso * Ks_dn * bdamp * exp(-(zwt_dn-1.5)/bdamp)
+                     Kl_fc = Kl_dn * bdamp * exp(-(zwt_dn-1.5)/bdamp)
                   ELSE
-                     Ks_fc = raniso * Ks_dn * ((1.5-zwt_dn) + bdamp)
+                     Kl_fc = Kl_dn * ((1.5-zwt_dn) + bdamp)
                   ENDIF
                ENDIF
 
-               ca = lenbdr * Ks_fc / theta_a_up / delp / area_up * deltime
-               cb = lenbdr * Ks_fc / theta_a_dn / delp / area_dn * deltime
+               ca = lenbdr * Kl_fc / theta_a_up / delp / area_up * deltime
+               cb = lenbdr * Kl_fc / theta_a_dn / delp / area_dn * deltime
 
-               xsubs_nb = (zsubs_up - zsubs_dn) * lenbdr * Ks_fc / (1+ca+cb) / delp
+               xsubs_nb = (zsubs_up - zsubs_dn) * lenbdr * Kl_fc / (1+ca+cb) / delp
 
                IF (.not. iam_lake) THEN
                   xsubs_nb = xsubs_nb / sum(hrus%agwt)
@@ -547,7 +554,7 @@ SUBROUTINE subsurface_flow (deltime)
 
          IF (allocated(theta_a_bsn)) deallocate(theta_a_bsn)
          IF (allocated(zwt_bsn    )) deallocate(zwt_bsn    )
-         IF (allocated(Ks_bsn     )) deallocate(Ks_bsn     )
+         IF (allocated(Kl_bsn     )) deallocate(Kl_bsn     )
 
       ENDIF
 
@@ -712,7 +719,7 @@ SUBROUTINE basin_neighbour_final ()
 
       IF (allocated(theta_a_nb)) deallocate(theta_a_nb)
       IF (allocated(zwt_nb    )) deallocate(zwt_nb    )
-      IF (allocated(Ks_nb     )) deallocate(Ks_nb     )
+      IF (allocated(Kl_nb     )) deallocate(Kl_nb     )
       IF (allocated(wdsrf_nb  )) deallocate(wdsrf_nb  )
       IF (allocated(agwt_nb   )) deallocate(agwt_nb   )
       IF (allocated(islake_nb )) deallocate(islake_nb )

main/LULCC/MOD_Lulcc_Driver.F90

@@ -27,6 +27,38 @@ SUBROUTINE LulccDriver (casename,dir_landdata,dir_restart,&
 !  07/2023, Wenzong Dong: porting to MPI version.
 !  08/2023, Wanyi Lin: add interface for Mass&Energy conserved scheme.
 !
+! Extra processes when adding a new variable and #define LULCC:
+!
+! 1. Save a copy of new variable (if called "var", save it to "var_") with 2 steps:
+!    1.1 main/LULCC/MOD_Lulcc_Vars_TimeVariables.F90's subroutine "allocate_LulccTimeVariables":
+!       allocate (var_(dimension))
+!    1.2 main/LULCC/MOD_Lulcc_Vars_TimeVariables.F90's subroutine "SAVE_LulccTimeVariables":
+!       var_ = var
+!
+! 2. Reassignment for the next year
+!    2.1 if used Same Type Assignment (SAT) scheme for variable recovery
+!       main/LULCC/MOD_Lulcc_Vars_TimeVariables.F90's subroutine "REST_LulccTimeVariables"
+!       var(np) = var_(np_)
+!
+!    2.2 if using Mass and Energy conservation (MEC) scheme for variable recovery
+!       2.2.1 main/LULCC/MOD_Lulcc_Vars_TimeVariables.F90's subroutine "REST_LulccTimeVariables":
+!          var(np) = var_(np_)
+!
+!       2.2.2 [No need for PFT/PC scheme] Mass and Energy conserve adjustment, add after line 519
+!             of MOD_Lulcc_MassEnergyConserve.F90.
+!
+!          o if variable should be mass conserved:
+!             var(:,np) = var(:,np) + var(:,frnp_(k))*lccpct_np(patchclass_(frnp_(k)))/sum_lccpct_np
+!
+!          o if variable should be energy conserved, take soil temperature "t_soisno" as an example: [May neeed extra calculation]
+!             t_soisno (1:nl_soil,np) = t_soisno (1:nl_soil,np) + &
+!                       t_soisno_(1:nl_soil,frnp_(k))*cvsoil_(1:nl_soil,k)*lccpct_np(patchclass_(frnp_(k)))/wgt(1:nl_soil)
+!             where cvsoil_ is the heat capacity, wgt is the sum of cvsoil_(1:nl_soil,k)*lccpct_np(patchclass_(frnp_(k))),
+!             which need to be calculated in advance.
+!
+! 3. Deallocate the copy of new variable in MOD_Lulcc_Vars_TimeVariables.F90's subroutine "deallocate_LulccTimeVariables":
+!    deallocate (var_)
+!
 !-----------------------------------------------------------------------
 
    USE MOD_Precision