Merge pull request #349 from CoLM-SYSU/master

Update from master

.github/workflows/build_CoLM_gnu.yml

@@ -22,9 +22,10 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, ubuntu-20.04]
+        os: [ubuntu-latest, ubuntu-22.04]
         mpi: ['mpich', 'openmpi', 'intelmpi']
         toolchain:
+          - {compiler: gcc, version: 14}
           - {compiler: gcc, version: 13}
           - {compiler: gcc, version: 12}
           - {compiler: gcc, version: 11}
@@ -34,10 +35,17 @@ jobs:
        #   - os: ubuntu-latest
        #     toolchain: {compiler: gcc, version: 12}
         exclude:
-          - os: ubuntu-20.04
-            toolchain: {compiler: gcc, version: 12}
+          - os: ubuntu-22.04
+            toolchain: {compiler: gcc, version: 14}
+          - mpi: 'mpich'
+            toolchain: {compiler: gcc, version: 9}
+            # mpich compiling FCFLAG has --fallow-type-mismatch, and will cause gfortran 9 failing
           - os: ubuntu-latest
             toolchain: {compiler: gcc, version: 9}
+          - os: ubuntu-latest
+            toolchain: {compiler: gcc, version: 10}
+          - os: ubuntu-latest
+            toolchain: {compiler: gcc, version: 11}
 
     steps:
      - uses: actions/checkout@v4
@@ -49,9 +57,9 @@ jobs:
        with:
           compiler: ${{ matrix.toolchain.compiler }}
           version: ${{ matrix.toolchain.version }} 
-     - name: Install netcdf-fortran library
+     - name: Install netcdf-fortran and blas library
        shell: bash -l {0}
-       run: sudo apt update && sudo apt install -y netcdf-bin libnetcdf-dev libnetcdff-dev
+       run: sudo apt update && sudo apt install -y netcdf-bin libnetcdf-dev libnetcdff-dev "liblapack*" "libblas*"
      - name: Test mpi
        #run: mpif90 -v
        run: which mpif90

CaMa/src/MOD_CaMa_Vars.F90

@@ -539,7 +539,7 @@ SUBROUTINE flux_map_and_write_2d_cama (is_hist, &
       IF (itime_in_file == 1) THEN
          CALL ncio_put_attr (file_hist, varname, 'long_name', longname)
          CALL ncio_put_attr (file_hist, varname, 'units', units)
-         CALL ncio_put_attr (file_hist, varname, 'missing_value',spval)
+         CALL ncio_put_attr (file_hist, varname, 'missing_value',real(real(spval,kind=JPRM),kind=8))
       ENDIF
 
    END SUBROUTINE flux_map_and_write_2d_cama

CaMa/src/MOD_CaMa_colmCaMa.F90

@@ -28,7 +28,7 @@ MODULE MOD_CaMa_colmCaMa
    USE CMF_CTRL_OUTPUT_MOD,       only: CMF_OUTPUT_INIT,CMF_OUTPUT_END,NVARSOUT,VAROUT
    USE YOS_CMF_INPUT,             only: NXIN, NYIN, DT,DTIN,IFRQ_INP,LLEAPYR,NX,NY,RMIS,DMIS
    USE MOD_Precision,             only: r8,r4
-   USE YOS_CMF_INPUT ,            only: LROSPLIT,LWEVAP,LWINFILT
+   USE YOS_CMF_INPUT ,            only: LROSPLIT,LWEVAP,LWINFILT,CSETFILE
    USE YOS_CMF_MAP,               only: D1LON, D1LAT
    USE YOS_CMF_INPUT,             only: WEST,EAST,NORTH,SOUTH
 
@@ -73,6 +73,11 @@ SUBROUTINE colm_CaMa_init
 #ifdef USEMPI
       CALL mpi_barrier (p_comm_glb, p_err)
 #endif
+
+      IF (p_is_master) THEN
+
+         CSETFILE = DEF_CaMa_Namelist
+
          !Namelist handling
          CALL CMF_DRV_INPUT
          !get the time information from colm namelist
@@ -181,14 +186,33 @@ SUBROUTINE colm_CaMa_init
             END SELECT
          ENDDO
 
+      ENDIF 
+
       !Broadcast the variables to all the processors
       CALL mpi_bcast (NX      ,   1, MPI_INTEGER,   p_address_master, p_comm_glb, p_err) ! number of grid points in x-direction of CaMa-Flood
       CALL mpi_bcast (NY      ,   1, MPI_INTEGER,   p_address_master, p_comm_glb, p_err) ! number of grid points in y-direction of CaMa-Flood
       CALL mpi_bcast (IFRQ_INP ,   1, MPI_INTEGER,  p_address_master, p_comm_glb, p_err) ! input frequency of CaMa-Flood (hour)
       CALL mpi_bcast (LWEVAP ,   1, MPI_LOGICAL,  p_address_master, p_comm_glb, p_err)   ! switch for inundation evaporation
       CALL mpi_bcast (LWINFILT ,   1, MPI_LOGICAL,  p_address_master, p_comm_glb, p_err) ! switch for inundation re-infiltration
-      CALL mpi_bcast (real(D1LAT,kind=8)    ,   1, MPI_REAL8,   p_address_master, p_comm_glb, p_err) ! 
-      CALL mpi_bcast (real(D1LON,kind=8)    ,   1, MPI_REAL8,   p_address_master, p_comm_glb, p_err)  !    
+
+      IF (.not. allocated(D1LAT))  allocate (D1LAT(NY))
+      IF (.not. allocated(D1LON))  allocate (D1LON(NX))
+
+#ifdef SinglePrec_CMF
+      CALL mpi_bcast (D1LAT, NY, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (D1LON, NX, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (SOUTH,  1, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (NORTH,  1, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (WEST ,  1, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (EAST ,  1, MPI_REAL4, p_address_master, p_comm_glb, p_err) !
+#else
+      CALL mpi_bcast (D1LAT, NY, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (D1LON, NX, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (SOUTH,  1, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (NORTH,  1, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (WEST ,  1, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+      CALL mpi_bcast (EAST ,  1, MPI_REAL8, p_address_master, p_comm_glb, p_err) !
+#endif
 
       !allocate the data structure for cama
       CALL gcama%define_by_center (D1LAT,D1LON,real(SOUTH,kind=8), real(NORTH,kind=8), real(WEST,kind=8), real(EAST,kind=8)) !define the grid for cama

CaMa/src/yos_cmf_input.F90

@@ -27,7 +27,7 @@ MODULE YOS_CMF_INPUT
    character(LEN=256)              :: CSETFILE                !! input namelist file name
 
    DATA LLOGOUT       /.TRUE./
-   DATA CLOGOUT       /'../run/log_CaMa.txt'/
+   DATA CLOGOUT       /'./log_CaMa.txt'/
    DATA CSETFILE      /'../run/cama_flood.nml'/
 
    !================================================

main/HYDRO/MOD_Catch_LateralFlow.F90

@@ -272,7 +272,7 @@ SUBROUTINE lateral_flow (deltime)
          CALL mpi_allreduce (MPI_IN_PLACE, toldis, 1, MPI_REAL8, MPI_SUM, p_comm_worker, p_err)
 #endif
          IF (p_iam_worker == 0) THEN
-            write(*,'(A,F10.5,A,ES10.3,A,ES10.3,A)') 'Total surface water error: ', dtolw, &
+            write(*,'(A,F10.2,A,ES10.3,A,ES10.3,A)') 'Total surface water error: ', dtolw, &
                '(m^3) in area ', landarea, '(m^2), discharge ', toldis, '(m^3)' 
          ENDIF
 
@@ -282,7 +282,7 @@ SUBROUTINE lateral_flow (deltime)
          CALL mpi_allreduce (MPI_IN_PLACE, dtolw,  1, MPI_REAL8, MPI_SUM, p_comm_worker, p_err)
 #endif
          IF (p_iam_worker == 0) THEN
-            write(*,'(A,F10.5,A,ES10.3,A)') 'Total ground  water error: ', dtolw, &
+            write(*,'(A,F10.2,A,ES10.3,A)') 'Total ground  water error: ', dtolw, &
                '(m^3) in area ', landarea, '(m^2)'
          ENDIF
       ENDIF

main/HYDRO/MOD_Catch_RiverLakeFlow.F90

@@ -454,21 +454,22 @@ SUBROUTINE river_lake_flow (dt)
                momen_riv_ta(:) = momen_riv_ta(:) + momen_riv(:) * dt_this
                discharge   (:) = discharge   (:) + hflux_fc (:) * dt_this
             ENDIF
+
+            DO i = 1, nbasin
+               IF (lake_id(i) > 0) THEN ! for lakes
+                  hs = basin_hru%substt(i)
+                  he = basin_hru%subend(i)
+                  DO j = hs, he
+                     wdsrf_hru(j) = max(wdsrf_bsn(i) - (lakes(i)%depth(1) - lakes(i)%depth0(j-hs+1)), 0.)
+                     wdsrf_hru_ta(j) = wdsrf_hru_ta(j) + wdsrf_hru(j) * dt_this
+                  ENDDO
+               ENDIF
+            ENDDO
 
             dt_res = dt_res - dt_this
 
          ENDDO
 
-         DO i = 1, nbasin
-            IF (lake_id(i) > 0) THEN ! for lakes
-               hs = basin_hru%substt(i)
-               he = basin_hru%subend(i)
-               DO j = hs, he
-                  wdsrf_hru(j) = max(wdsrf_bsn(i) - (lakes(i)%depth(1) - lakes(i)%depth0(j-hs+1)), 0.)
-               ENDDO
-            ENDIF
-         ENDDO
-
          wdsrf_bsn_prev(:) = wdsrf_bsn(:)
 
          IF (allocated(wdsrf_bsn_ds )) deallocate(wdsrf_bsn_ds )

main/HYDRO/MOD_Catch_RiverLakeNetwork.F90

@@ -682,6 +682,9 @@ SUBROUTINE river_lake_network_init ()
 
                   IF (inb <= 0) THEN
                      outletwth(ibasin) = 0
+                     IF (riverdown(ibasin) > 0) THEN
+                        outletwth(ibasin) = 90.
+                     ENDIF
                   ELSE
                      outletwth(ibasin) = elementneighbour(ibasin)%lenbdr(inb)
                   ENDIF

main/MOD_Const_LC.F90

@@ -596,7 +596,7 @@ MODULE MOD_Const_LC
    real(r8), parameter, dimension(N_land_classification) :: ck0_igbp &
       = (/3.95, 3.95, 3.95, 3.95, 3.95, 3.95, &
           3.95, 3.95, 3.95, 3.95, 3.95, 3.95, &
-          3.95, 3.95,   0.,   0.,   0.  /)
+          3.95, 3.95, 3.95, 3.95, 3.95  /)
 !end plant hydraulic parameters
 
    ! lambda for WUE stomata model

main/MOD_Const_PFT.F90

@@ -1604,7 +1604,7 @@ MODULE MOD_Const_PFT
 
    ! shape-fitting parameter for vulnerability curve (-)
    real(r8), parameter :: ck_p(0:N_PFT+N_CFT-1) &
-      = (/  0.,  3.95, 3.95,  3.95, 3.95,  3.95, 3.95, 3.95&
+      = (/3.95,  3.95, 3.95,  3.95, 3.95,  3.95, 3.95, 3.95&
          ,3.95,  3.95, 3.95,  3.95, 3.95,  3.95, 3.95, 3.95&
 #ifdef CROP
          ,3.95,  3.95, 3.95,  3.95, 3.95,  3.95, 3.95, 3.95&

main/MOD_SoilSnowHydrology.F90

@@ -683,6 +683,7 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
    type(cell_data_struct) :: cell
    real(r8) :: wliq_soisno_tmp(1:nl_soil)
 
+   real(r8), parameter :: e_ice=6.0      !soil ice impedance factor
 
 !=======================================================================
 ! [1] update the liquid water within snow layer and the water onto soil
@@ -1020,6 +1021,16 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
       ENDIF
 #endif
 
+      DO j = 1, nl_soil
+         IF(t_soisno(j) <= tfrz) THEN
+            ! consider impedance factor
+            vol_ice(j) = max(min(porsl(j), wice_soisno(j)/(dz_soisno(j)*denice)), 0.)
+            icefrac(j) = vol_ice(j)/porsl(j)
+            imped = 10.**(-e_ice*icefrac(j))
+            hk(j) = imped * hk(j)
+         ENDIF
+      ENDDO
+
 #ifndef CatchLateralFlow
       err_solver = (sum(wliq_soisno(1:))+sum(wice_soisno(1:))+wa+wdsrf) - w_sum &
          - (gwat-etr-rsur-rsubst)*deltim
@@ -1044,7 +1055,8 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
 
 #if(defined CoLMDEBUG)
       IF(abs(err_solver) > 1.e-3)THEN
-         write(6,'(A,E20.5,I0)') 'Warning (WATER_VSF): water balance violation', err_solver,landpatch%eindex(ipatch)
+         write(6,'(A,E20.5,A,I0)') 'Warning (WATER_VSF): water balance violation', err_solver, & 
+            ' in element ', landpatch%eindex(ipatch)
       ENDIF
       IF (any(wliq_soisno < -1.e-3)) THEN
          write(6,'(A,10E20.5)') 'Warning (WATER_VSF): negative soil water', wliq_soisno(1:nl_soil)

mkinidata/MOD_Initialize.F90

@@ -725,8 +725,15 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
       IF (DEF_USE_SoilInit) THEN
 
          fsoildat = DEF_file_SoilInit
+
          IF (p_is_master) THEN
             inquire (file=trim(fsoildat), exist=use_soilini)
+            IF (use_soilini) THEN
+               write(*,'(/,2A)') 'Use soil water content, soil temperature and water table depth ' &
+                  // 'to initialize soil state from file ', trim(fsoildat)
+            ELSE
+               write(*,*) 'No initial data for soil state from ', trim(fsoildat)
+            ENDIF
          ENDIF
 #ifdef USEMPI
          CALL mpi_bcast (use_soilini, 1, MPI_LOGICAL, p_address_master, p_comm_glb, p_err)
@@ -1083,19 +1090,25 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
 
 
       ! for SOIL Water INIT by using water table depth
-      fwtd = trim(DEF_dir_runtime) // '/wtd.nc'
-      IF (p_is_master) THEN
-         inquire (file=trim(fwtd), exist=use_wtd)
-         IF (use_soilini) use_wtd = .false.
-         IF (use_wtd) THEN
-            write(*,'(/, 2A)') 'Use water table depth and derived equilibrium state ' &
-               // ' to initialize soil water content: ', trim(fwtd)
-         ENDIF
-      ENDIF
+      use_wtd = (.not. use_soilini) .and. DEF_USE_WaterTableInit
+
+      IF (use_wtd) THEN
 
+         fwtd = DEF_file_WaterTable 
+
+         IF (p_is_master) THEN
+            inquire (file=trim(fwtd), exist=use_wtd)
+            IF (use_wtd) THEN
+               write(*,'(/, 2A)') 'Use water table depth and derived equilibrium state ' &
+                  // ' to initialize soil water content from file ', trim(fwtd)
+            ELSE
+               write(*,*) 'No initial data for water table depth from ', trim(fwtd)
+            ENDIF
+         ENDIF
 #ifdef USEMPI
-      CALL mpi_bcast (use_wtd, 1, MPI_LOGICAL, p_address_master, p_comm_glb, p_err)
+         CALL mpi_bcast (use_wtd, 1, MPI_LOGICAL, p_address_master, p_comm_glb, p_err)
 #endif
+      ENDIF
 
       IF (use_wtd) THEN
 

run/forcing/CRUJRA.nml

@@ -18,7 +18,7 @@
    DEF_forcing%endyr              = 2022    ! end year of forcing data
    DEF_forcing%endmo              = 12      ! end month of forcing data
    DEF_forcing%dtime              = 21600 21600 21600 21600 21600 21600 21600 21600
-   DEF_forcing%offset             =   0    0     0    10800    0     0     0   10800
+   DEF_forcing%offset             = 10800 10800 10800 10800 10800 10800 0 10800
    DEF_forcing%nlands             = 1       ! land grid number in 1d
 
    DEF_forcing%leapyear           = .false. ! leapyear calendar

share/MOD_Namelist.F90

@@ -203,6 +203,9 @@ MODULE MOD_Namelist
    logical :: DEF_USE_CN_INIT   = .false.
    character(len=256) :: DEF_file_cn_init  = 'null'
 
+   logical :: DEF_USE_WaterTableInit = .false.
+   character(len=256) :: DEF_file_WaterTable = 'null'
+
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! ----- Part 9: LULCC related ------
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -301,6 +304,9 @@ MODULE MOD_Namelist
    !2: To allow monthly ndep data to be read in
    integer :: DEF_NDEP_FREQUENCY = 1
 
+   ! ----- CaMa-Flood -----
+   character(len=256) :: DEF_CaMa_Namelist = 'null'
+
    ! ----- lateral flow related -----
    logical :: DEF_USE_EstimatedRiverDepth     = .true.
    character(len=256) :: DEF_ElementNeighbour_file = 'null'
@@ -993,9 +999,14 @@ SUBROUTINE read_namelist (nlfile)
       DEF_USE_CN_INIT,                        &
       DEF_file_cn_init,                       &
 
+      DEF_USE_WaterTableInit,                 &
+      DEF_file_WaterTable,                    &
+
       DEF_file_snowoptics,                    &
       DEF_file_snowaging ,                    &
 
+      DEF_CaMa_Namelist,                      &
+
       DEF_ElementNeighbour_file,              &
 
       DEF_DA_obsdir,                          &
@@ -1305,6 +1316,15 @@ SUBROUTINE read_namelist (nlfile)
 #endif
 #endif
 
+! ----- Soil water and temperature Initialization ----- Namelist conflicts
+
+         IF (DEF_USE_SoilInit .and. DEF_USE_WaterTableInit) THEN
+            write(*,*) '                  *****                  '
+            write(*,*) 'If both DEF_USE_SoilInit and DEF_USE_WaterTableInit are .TRUE.,   '
+            write(*,*) 'initial value of water table depth is read from DEF_file_SoilInit,'
+            write(*,*) 'instead of DEF_file_WaterTable (which is useless in this CASE).   '
+         ENDIF 
+
 ! ----- dynamic lake run ----- Macros&Namelist conflicts and dependency management
 
 #ifndef CATCHMENT
@@ -1477,10 +1497,15 @@ SUBROUTINE read_namelist (nlfile)
       CALL mpi_bcast (DEF_USE_CN_INIT                        ,1   ,mpi_logical   ,p_address_master ,p_comm_glb ,p_err)
       CALL mpi_bcast (DEF_file_cn_init                       ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
 
+      CALL mpi_bcast (DEF_USE_WaterTableInit                 ,1   ,mpi_logical   ,p_address_master ,p_comm_glb ,p_err)
+      CALL mpi_bcast (DEF_file_WaterTable                    ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
+
       CALL mpi_bcast (DEF_USE_SNICAR                         ,1   ,mpi_logical   ,p_address_master ,p_comm_glb ,p_err)
       CALL mpi_bcast (DEF_file_snowoptics                    ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
       CALL mpi_bcast (DEF_file_snowaging                     ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
 
+      CALL mpi_bcast (DEF_CaMa_Namelist                      ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
+
       CALL mpi_bcast (DEF_ElementNeighbour_file              ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)
 
       CALL mpi_bcast (DEF_DA_obsdir                          ,256 ,mpi_character ,p_address_master ,p_comm_glb ,p_err)

share/MOD_SPMD_Task.F90

@@ -331,11 +331,12 @@ SUBROUTINE CoLM_stop (mesg)
 
    IMPLICIT NONE
    character(len=*), optional :: mesg
+   integer :: errorcode
 
       IF (present(mesg)) write(*,*) trim(mesg)
 
 #ifdef USEMPI
-      CALL mpi_abort (p_comm_glb, p_err)
+      CALL mpi_abort (p_comm_glb, errorcode, p_err)
 #else
       STOP
 #endif