Merge pull request #363 from zhangsp8/master

Bugs fixed before release.

Makefile

@@ -61,6 +61,8 @@ OBJS_SHARED =    \
 				  MOD_SrfdataRestart.o         \
 				  MOD_ElmVector.o              \
 				  MOD_HRUVector.o              \
+				  MOD_MeshFilter.o             \
+				  MOD_RegionClip.o             \
 				  MOD_Urban_Const_LCZ.o
 
 ${OBJS_SHARED} : %.o : %.F90 ${HEADER}
@@ -80,8 +82,6 @@ OBJS_MKSRFDATA = \
 				  Aggregation_TopographyFactors.o   \
 				  Aggregation_Urban.o               \
 				  Aggregation_SoilTexture.o         \
-				  MOD_MeshFilter.o                  \
-				  MOD_RegionClip.o                  \
 				  MKSRFDATA.o
 
 $(OBJS_MKSRFDATA) : %.o : %.F90 ${HEADER} ${OBJS_SHARED}

main/HYDRO/MOD_Hydro_SoilWater.F90

@@ -159,13 +159,12 @@ SUBROUTINE get_water_equilibrium_state ( &
    END SUBROUTINE get_water_equilibrium_state
 
    ! --- soil water movement ---
-   SUBROUTINE soil_water_vertical_movement (         &
-         nlev,  dt,   sp_zc, sp_zi,   is_permeable,  &
-         porsl, vl_r, psi_s, hksat,   nprm,   prms,  &
-         porsl_wa,                                   &
-         qgtop, etr,  rootr, rootflux, rsubst,  qinfl,         &
-         ss_dp, zwt,  wa,    ss_vliq, smp,    hk  ,  &
-         tolerance )
+   SUBROUTINE soil_water_vertical_movement (                      &
+         nlev,       dt,    sp_zc,  sp_zi,    is_permeable,  porsl,    &
+         vl_r,       psi_s, hksat,  nprm,     prms,          porsl_wa, &
+         qgtop,      etr,   rootr,  rootflux, rsubst,        qinfl,    &
+         ss_dp,      zwt,   wa,     ss_vliq,  smp,           hk,       &
+         tolerance,  wblc)
 
    !=======================================================================
    ! this is the main subroutine to execute the calculation of
@@ -214,9 +213,11 @@ SUBROUTINE soil_water_vertical_movement (         &
 
    real(r8), intent(in) :: tolerance
 
+   real(r8), intent(out) :: wblc
+
    ! Local variables
    integer  :: lb, ub, ilev, izwt
-   real(r8) :: sumroot, deficit, wexchange
+   real(r8) :: sumroot, deficit, etrdef, wexchange
    real(r8) :: dp_m1, psi, vliq, zwtp, air
    logical  :: is_sat
 
@@ -230,9 +231,7 @@ SUBROUTINE soil_water_vertical_movement (         &
    integer  :: lbc_typ_sub
    real(r8) :: lbc_val_sub
 
-#ifdef CoLMDEBUG
-   real(r8) :: w_sum_before, w_sum_after, wblc
-#endif
+   real(r8) :: w_sum_before, w_sum_after
 
    real(r8) :: tol_q, tol_z, tol_v, tol_p
 
@@ -249,7 +248,6 @@ SUBROUTINE soil_water_vertical_movement (         &
       ! water table location
       izwt = findloc_ud(zwt >= sp_zi, back=.true.)
 
-#ifdef CoLMDEBUG
       ! total water mass
       w_sum_before = ss_dp
       DO ilev = 1, nlev
@@ -265,7 +263,6 @@ SUBROUTINE soil_water_vertical_movement (         &
          ENDIF
       ENDDO
       w_sum_before = w_sum_before + wa
-#endif
 
       ! transpiration
       IF(.not. DEF_USE_PLANTHYDRAULICS)THEN
@@ -275,15 +272,17 @@ SUBROUTINE soil_water_vertical_movement (         &
             WHERE (is_permeable)
                etroot = etr * max(rootr, 0.) / sumroot
             END WHERE
-            deficit = 0.
+            etrdef = 0.
          ELSE
-            deficit = etr*dt
+            etrdef = etr*dt
          ENDIF
       ELSE
-         deficit = 0.
+         etrdef = 0.
          etroot(:) = rootflux
       ENDIF
 
+      deficit = etrdef
+
       DO ilev = 1, izwt-1
          IF (is_permeable(ilev)) THEN
 
@@ -413,7 +412,6 @@ SUBROUTINE soil_water_vertical_movement (         &
 
       qinfl = qgtop - (ss_dp - dp_m1)/dt
 
-#ifdef CoLMDEBUG
       ! total water mass
       w_sum_after = ss_dp
       DO ilev = 1, nlev
@@ -430,13 +428,14 @@ SUBROUTINE soil_water_vertical_movement (         &
       ENDDO
       w_sum_after = w_sum_after + wa
 
-      wblc = w_sum_after - (w_sum_before + (qgtop - sum(etroot) - rsubst) * dt)
+      wblc = w_sum_after - (w_sum_before + (qgtop - sum(etroot) - rsubst) * dt - etrdef)
 
       IF (abs(wblc) > tolerance) THEN
-         write(*,*) 'soil_water_vertical_movement balance error: ', wblc
-         write(*,*) w_sum_after, w_sum_before, qgtop, etr, rsubst, is_permeable(1), ss_dp
+         write(*,*) 'soil_water_vertical_movement balance error: ', wblc, ' in mm.'
+         write(*,*) 'qtop: ', qgtop, 'etr: ', etr, 'rsubst: ', rsubst, 'surf dep: ', ss_dp
+         write(*,*) 'permeable (1-10): ', is_permeable
+         write(*,*) 'vliq (1-10): ', ss_vliq
       ENDIF
-#endif
 
       DO ilev = 1, nlev
          IF (ilev < izwt) THEN

main/MOD_Eroot.F90

@@ -62,7 +62,7 @@ SUBROUTINE eroot (nl_soil,trsmx0,porsl, &
    real(r8), intent(in) :: sc_vgm   (1:nl_soil)
    real(r8), intent(in) :: fc_vgm   (1:nl_soil)
 #endif
-   real(r8), intent(in) :: psi0(1:nl_soil)         ! saturated soil suction (mm) (NEGATIVE)
+   real(r8), intent(in) :: psi0(1:nl_soil)         ! saturated soil suction (cm) (NEGATIVE)
    real(r8), intent(in) :: rootfr(1:nl_soil)       ! fraction of roots in a layer,
    real(r8), intent(in) :: dz_soisno(1:nl_soil)    ! layer thickness (m)
    real(r8), intent(in) :: t_soisno(1:nl_soil)     ! soil/snow skin temperature (K)
@@ -77,7 +77,7 @@ SUBROUTINE eroot (nl_soil,trsmx0,porsl, &
    real(r8) roota             ! accumulates root resistance factors
    real(r8) rresis(1:nl_soil) ! soil water contribution to root resistance
    real(r8) s_node            ! vol_liq/porosity
-   real(r8) smpmax            ! wilting point potential in mm
+   real(r8) smpmax            ! wilting point potential in cm
    real(r8) smp_node          ! matrix potential
 
    integer i                  ! loop counter

main/MOD_Forcing.F90

@@ -27,6 +27,7 @@ MODULE MOD_Forcing
    USE MOD_MonthlyinSituCO2MaunaLoa
    USE MOD_Vars_Global, only : pi
    USE MOD_OrbCoszen
+   USE MOD_UserDefFun
 
    IMPLICIT NONE
 
@@ -873,7 +874,7 @@ SUBROUTINE read_forcing (idate, dir_forcing)
          IF (p_is_worker) THEN
             DO j = 1, numpatch
                   a = forc_swrad(j)
-                  IF (isnan(a)) a = 0
+                  IF (isnan_ud(a)) a = 0
                   calday = calendarday(idate)
                   sunang = orb_coszen (calday, patchlonr(j), patchlatr(j))
                   IF (sunang.eq.0) THEN

main/MOD_ForcingDownscaling.F90

@@ -20,6 +20,7 @@ MODULE MOD_ForcingDownscaling
    USE MOD_Namelist
    USE MOD_Const_Physical
    USE MOD_Vars_Global
+   USE MOD_UserDefFun
 
    IMPLICIT NONE
 
@@ -657,7 +658,7 @@ SUBROUTINE downscale_shortwave( &
          tcf_type(i) = (1+cos(slp_type_c(i)))/2-svf
          IF (tcf_type(i)<0) tcf_type(i) = 0
 
-         IF (isnan(alb)) THEN
+         IF (isnan_ud(alb)) THEN
             refl_swrad_type(i) = -1.0e36
          ELSE
             IF ((balb<0).or.(balb>1)) balb = 0

main/MOD_SoilSnowHydrology.F90

@@ -658,6 +658,7 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
    logical  :: is_permeable(1:nl_soil)
    real(r8) :: dzsum, dz
    real(r8) :: icefracsum, fracice_rsub, imped
+   real(r8) :: wblc
 
 #ifdef CROP
    integer  :: ps, pe
@@ -967,7 +968,7 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
          eff_porosity(1:nl_soil), theta_r(1:nl_soil), psi0(1:nl_soil), hksati(1:nl_soil), &
          nprms, prms(:,1:nl_soil), porsl(nl_soil),     &
          qgtop, etr, rootr(1:nl_soil), rootflux(1:nl_soil), rsubst, qinfl, &
-         wdsrf, zwtmm, wa, vol_liq(1:nl_soil), smp(1:nl_soil), hk(1:nl_soil), 1.e-3)
+         wdsrf, zwtmm, wa, vol_liq(1:nl_soil), smp(1:nl_soil), hk(1:nl_soil), 1.e-3, wblc)
 
       ! update the mass of liquid water
       DO j = nl_soil, 1, -1
@@ -1001,6 +1002,22 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
       wice_soisno(1) = max(0., wice_soisno(1) + (qfros_soil-qsubl_soil) * deltim)
 ENDIF
 
+      ! water imbalance mainly due to insufficient liquid water for evapotranspiration
+      IF (wblc > 0.) THEN
+         DO j = 1, nl_soil
+            IF (wice_soisno(j) > wblc) THEN
+               wice_soisno(j) = wice_soisno(j) - wblc
+               wblc = 0.
+               EXIT
+            ELSE
+               wblc = wblc - wice_soisno(j)
+               wice_soisno(j) = 0.
+            ENDIF
+         ENDDO
+
+         IF (wblc > 0.) wa = wa - wblc
+      ENDIF
+
 #ifndef CatchLateralFlow
       IF (.not. is_dry_lake) THEN
          IF (wdsrf > pondmx) THEN
@@ -1074,7 +1091,11 @@ SUBROUTINE WATER_VSF (ipatch, patchtype,is_dry_lake, lb, nl_soil, deltim   ,&
 
 
          IF (lb >= 1) THEN
-            wetwat = wdsrf + wa + wetwat + (gwat - etr + qsdew + qfros - qsubl) * deltim
+            IF (.not.DEF_SPLIT_SOILSNOW) THEN
+               wetwat = wdsrf + wa + wetwat + (gwat - etr + qsdew + qfros - qsubl) * deltim
+            ELSE
+               wetwat = wdsrf + wa + wetwat + (gwat - etr + qsdew_soil + qfros_soil - qsubl_soil) * deltim
+            ENDIF
          ELSE
             wetwat = wdsrf + wa + wetwat + (gwat - etr) * deltim
          ENDIF

mkinidata/MOD_Initialize.F90

@@ -236,6 +236,8 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
    real(r8)           :: filval
    type(block_data_real8_2d)  :: fsatmax_grid, fsatdcf_grid
    type(spatial_mapping_type) :: map_simtop_para
+
+   logical :: use_soiltext
    ! 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; &    
@@ -342,12 +344,19 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, &
       ENDIF
 #endif
 
-      CALL soiltext_readin (dir_landdata, lc_year)
-      IF (p_is_worker) THEN
-         IF (numpatch > 0) THEN
-            DO ipatch = 1, numpatch
-               BVIC(ipatch)=BVIC_USDA(soiltext(ipatch))
-            ENDDO
+#ifdef CatchLateralFlow
+      use_soiltext = .true.
+#else
+      use_soiltext = (DEF_Runoff_SCHEME == 3) ! only for Simple VIC
+#endif
+      IF (use_soiltext) THEN
+         CALL soiltext_readin (dir_landdata, lc_year)
+         IF (p_is_worker) THEN
+            IF (numpatch > 0) THEN
+               DO ipatch = 1, numpatch
+                  BVIC(ipatch)=BVIC_USDA(soiltext(ipatch))
+               ENDDO
+            ENDIF
          ENDIF
       ENDIF
 

mksrfdata/Aggregation_SoilTexture.F90

@@ -71,7 +71,7 @@ SUBROUTINE Aggregation_SoilTexture ( &
 #endif
 
 #ifdef SinglePoint
-      IF (USE_SITE_soilparameters) THEN
+      IF (USE_SITE_soilparameters .or. (DEF_Runoff_SCHEME /= 3)) THEN
          RETURN
       ENDIF
 #endif

mksrfdata/MOD_SingleSrfdata.F90

@@ -276,12 +276,13 @@ SUBROUTINE read_surface_data_single (fsrfdata, mksrfdata)
          CALL ncio_read_serial (fsrfdata, 'soil_alpha_vgm        ', SITE_soil_alpha_vgm        )
          CALL ncio_read_serial (fsrfdata, 'soil_L_vgm            ', SITE_soil_L_vgm            )
          CALL ncio_read_serial (fsrfdata, 'soil_n_vgm            ', SITE_soil_n_vgm            )
-#else
-         !SITE_soil_theta_r(:) = 0.
 #endif
          CALL ncio_read_serial (fsrfdata, 'soil_BA_alpha         ', SITE_soil_BA_alpha         )
          CALL ncio_read_serial (fsrfdata, 'soil_BA_beta          ', SITE_soil_BA_beta          )
-         CALL ncio_read_serial (fsrfdata, 'soil_texture          ', SITE_soil_texture          )
+
+         IF (DEF_Runoff_SCHEME == 3) THEN ! for Simple VIC
+            CALL ncio_read_serial (fsrfdata, 'soil_texture       ', SITE_soil_texture          )
+         ENDIF
       ENDIF
 
       IF (DEF_USE_BEDROCK) THEN
@@ -463,7 +464,10 @@ SUBROUTINE read_urban_surface_data_single (fsrfdata, mksrfdata, mkrun)
 #endif
          CALL ncio_read_serial (fsrfdata, 'soil_BA_alpha         ', SITE_soil_BA_alpha         )
          CALL ncio_read_serial (fsrfdata, 'soil_BA_beta          ', SITE_soil_BA_beta          )
-         CALL ncio_read_serial (fsrfdata, 'soil_texture          ', SITE_soil_texture          )
+
+         IF (DEF_Runoff_SCHEME == 3) THEN ! for Simple VIC
+            CALL ncio_read_serial (fsrfdata, 'soil_texture       ', SITE_soil_texture          )
+         ENDIF
       ENDIF
 
       IF (DEF_USE_BEDROCK) THEN
@@ -643,9 +647,11 @@ SUBROUTINE write_surface_data_single (numpatch, numpft)
       CALL ncio_write_serial (fsrfdata, 'soil_BA_beta ', SITE_soil_BA_beta , 'soil')
       CALL ncio_put_attr     (fsrfdata, 'soil_BA_alpha', 'source', source)
       CALL ncio_put_attr     (fsrfdata, 'soil_BA_beta ', 'source', source)
-
-      CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture)
-      CALL ncio_put_attr     (fsrfdata, 'soil_texture ', 'source', source)
+
+      IF (DEF_Runoff_SCHEME == 3) THEN ! for Simple VIC
+         CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture)
+         CALL ncio_put_attr     (fsrfdata, 'soil_texture ', 'source', source)
+      ENDIF
 
       IF(DEF_USE_BEDROCK)THEN
          CALL ncio_write_serial (fsrfdata, 'depth_to_bedrock', SITE_dbedrock)
@@ -842,8 +848,10 @@ SUBROUTINE write_urban_surface_data_single (numurban)
       CALL ncio_put_attr     (fsrfdata, 'soil_BA_alpha', 'source', source)
       CALL ncio_put_attr     (fsrfdata, 'soil_BA_beta ', 'source', source)
 
-      CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture)
-      CALL ncio_put_attr     (fsrfdata, 'soil_texture ', 'source', source)
+      IF (DEF_Runoff_SCHEME == 3) THEN ! for Simple VIC
+         CALL ncio_write_serial (fsrfdata, 'soil_texture ', SITE_soil_texture)
+         CALL ncio_put_attr     (fsrfdata, 'soil_texture ', 'source', source)
+      ENDIF
 
       IF(DEF_USE_BEDROCK)THEN
          CALL ncio_write_serial (fsrfdata, 'depth_to_bedrock', SITE_dbedrock)

share/MOD_Namelist.F90

@@ -1074,6 +1074,11 @@ SUBROUTINE read_namelist (nlfile)
          DEF_HIST_mode = 'one'
 #endif
 
+         IF (DEF_simulation_time%timestep > 3600.) THEN
+            write(*,*) '                  *****                  '
+            write(*,*) 'Warning: timestep should be less than or equal to 3600 seconds.'
+            CALL CoLM_Stop ()
+         ENDIF
 
 ! ===============================================================
 ! ----- Macros&Namelist conflicts and dependency management -----

share/MOD_SPMD_Task.F90

@@ -218,7 +218,7 @@ SUBROUTINE divide_processes_into_groups (ngrp)
       ELSE
          p_is_io     = .false.
          p_is_worker = .false.
-         p_my_group  = -1
+         p_my_group  = p_np_glb
       ENDIF
 
       ! 3. Construct IO communicator and address book.