add heatcond1d_cn_dense for densening grids inside

src/suews/src/suews_ctrl_driver.f95

@@ -1525,7 +1525,7 @@ SUBROUTINE SUEWS_cal_Main( &
             flag_converge = .FALSE.
          ELSE
             flag_converge = .TRUE.
-            ! PRINT *, 'Iteration done in', i_iter, ' iterations'
+            PRINT *, 'Iteration done in', i_iter, ' iterations'
             ! PRINT *, ' qh_residual: ', qh_residual, ' qh_resist: ', qh_resist
             ! PRINT *, ' dif_qh: ', ABS(qh_residual - qh_resist)
             ! PRINT *, ' abs. dif_tsfc: ', dif_tsfc_iter
@@ -1535,11 +1535,11 @@ SUBROUTINE SUEWS_cal_Main( &
          i_iter = i_iter + 1
          ! force quit do-while loop if not convergent after 100 iterations
          IF (Diagnose == 1 .AND. i_iter == max_iter) THEN
-            ! PRINT *, 'Iteration did not converge in', i_iter, ' iterations'
+            PRINT *, 'Iteration did not converge in', i_iter, ' iterations'
             ! PRINT *, ' qh_residual: ', qh_residual, ' qh_resist: ', qh_resist
             ! PRINT *, ' dif_qh: ', ABS(qh_residual - qh_resist)
-            ! PRINT *, ' Ts_iter: ', Ts_iter, ' TSfc_C: ', TSfc_C
-            ! PRINT *, ' abs. dif_tsfc: ', dif_tsfc_iter
+            PRINT *, ' Ts_iter: ', Ts_iter, ' TSfc_C: ', TSfc_C
+            PRINT *, ' abs. dif_tsfc: ', dif_tsfc_iter
             ! exit
          END IF
 

src/suews/src/suews_phys_ehc.f95

@@ -116,6 +116,7 @@ SUBROUTINE heatcond1d_CN(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)
       REAL(KIND(1D0)), INTENT(out) :: Qs, Tsfc
       LOGICAL, INTENT(in) :: bctype(2) ! if true, use surrogate flux as boundary condition
       REAL(KIND(1D0)) :: alpha, T_lw, T_up
+      REAL(KIND(1D0)) :: Qs_acc
 
       LOGICAL, INTENT(in) :: debug
       INTEGER :: i, n
@@ -155,10 +156,12 @@ SUBROUTINE heatcond1d_CN(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)
       dt_remain = dt
       dt_step_cfl = 0.002*MINVAL(dx**2/(k/rhocp))
       !PRINT *, 'dt_step_cfl: ', dt_step_cfl
+
+      Qs_acc = 0.0    ! accumulated heat storage
       DO WHILE (dt_remain > 1E-10)
          dt_step = MIN(dt_step_cfl, dt_remain)
-         T_tmp(1) = T_up
-         T_tmp(n) = T_lw
+         !T_tmp(1) = T_up
+         !T_tmp(n) = T_lw
          ! set the tridiagonal matrix for 1D heat conduction solver based on Crank-Nicholson method
          DO i = 1, n
             IF (i == 1) THEN
@@ -194,6 +197,7 @@ SUBROUTINE heatcond1d_CN(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)
          ! solve the tridiagonal matrix using Thomas algorithm
          CALL thomas_triMat(vec_lw, vec_diag, vec_up, vec_rhs, n, T_tmp)
          dt_remain = dt_remain - dt_step
+         Qs_acc = Qs_acc + (T_up - T_tmp(1)) * k(1) / (dx(1)*0.5) * dt_step
       END DO
 
       T_out = T_tmp
@@ -219,8 +223,142 @@ SUBROUTINE heatcond1d_CN(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)
       ! ------considering there might be fluxes going out from the lower boundary
       Qs = (T_up - T_out(1))*k(1)/(dx(1)*0.5)
       ! Qs = (T_out(1) - T_out(2)) * k(1) / dx(1)
+      ! Qs = Qs_acc / dt
    END SUBROUTINE heatcond1d_CN
 
+
+   SUBROUTINE heatcond1d_CN_dense(T, Qs, Tsfc, dx, dt, k, rhocp, bc, bctype, debug)
+      REAL(KIND(1D0)), INTENT(inout) :: T(:)
+      REAL(KIND(1D0)), INTENT(in) :: dx(:), dt, k(:), rhocp(:), bc(2)
+      REAL(KIND(1D0)), INTENT(out) :: Qs, Tsfc
+      LOGICAL, INTENT(in) :: bctype(2) ! if true, use surrogate flux as boundary condition
+      REAL(KIND(1D0)) :: alpha, T_lw, T_up
+      REAL(KIND(1D0)) :: Qs_acc
+
+      LOGICAL, INTENT(in) :: debug
+      INTEGER :: i, ids_new, n, n_tmp
+      INTEGER :: ratio
+      REAL(KIND(1D0)), ALLOCATABLE :: T_tmp(:), k_itf(:)
+      REAL(KIND(1D0)), ALLOCATABLE :: T_in(:), T_out(:)
+      REAL(KIND(1D0)), ALLOCATABLE :: k_tmp(:), rhocp_tmp(:), dx_tmp(:)
+      REAL(KIND(1D0)), ALLOCATABLE :: vec_lw(:), vec_up(:), vec_diag(:), vec_rhs(:)
+
+      REAL(KIND(1D0)) :: dt_remain
+      REAL(KIND(1D0)) :: dt_step
+      REAL(KIND(1D0)) :: dt_step_cfl
+
+      ratio = 3   ! ratio between the number of layers in the used dense and input coarse grids
+      n = SIZE(T) 
+      n_tmp = n * ratio
+
+      ALLOCATE (T_tmp(1:n_tmp)) ! temporary temperature array
+      ALLOCATE (T_in(1:n_tmp)) ! initial temperature array
+      ALLOCATE (k_itf(1:n_tmp - 1)) ! thermal conductivity at interfaces
+      ALLOCATE (T_out(1:n)) ! output temperature array
+      ALLOCATE (k_tmp(1:n_tmp))
+      ALLOCATE (rhocp_tmp(1:n_tmp))
+      ALLOCATE (dx_tmp(1:n_tmp))
+
+      ALLOCATE (vec_lw(1:n_tmp - 1))
+      ALLOCATE (vec_up(1:n_tmp - 1))
+      ALLOCATE (vec_diag(1:n_tmp))
+      ALLOCATE (vec_rhs(1:n_tmp))
+
+      alpha = 0.5
+      T_up = bc(1)
+      T_lw = bc(2)
+
+      ! save initial temperatures
+      DO i = 1, n
+         T_in((i - 1)*ratio + 1:i*ratio) = T(i)
+         T_tmp((i - 1)*ratio + 1:i*ratio) = T(i)
+         k_tmp((i - 1)*ratio + 1:i*ratio) = k(i)
+         rhocp_tmp((i - 1)*ratio + 1:i*ratio) = rhocp(i)
+         dx_tmp((i - 1)*ratio + 1:i*ratio) = dx(i) / ratio
+      END DO
+
+      ! calculate the depth-averaged thermal conductivity
+      DO i = 1, n_tmp - 1
+         k_itf(i) = (k_tmp(i)*k_tmp(i + 1)*(dx_tmp(i) + dx_tmp(i + 1)))/(k_tmp(i)*dx_tmp(i + 1) + k_tmp(i + 1)*dx_tmp(i))
+      END DO
+
+      dt_remain = dt
+      dt_step_cfl = 0.01*MINVAL(dx_tmp**2/(k_tmp/rhocp_tmp))
+      !PRINT *, 'dt_step_cfl: ', dt_step_cfl
+
+      Qs_acc = 0.0    ! accumulated heat storage
+      DO WHILE (dt_remain > 1E-10)
+         dt_step = MIN(dt_step_cfl, dt_remain)
+         !T_tmp(1) = T_up
+         !T_tmp(n) = T_lw
+         ! set the tridiagonal matrix for 1D heat conduction solver based on Crank-Nicholson method
+         DO i = 1, n_tmp
+            IF (i == 1) THEN
+               vec_up(i) = (1.0 - alpha) * k_itf(i) / (0.5*(dx_tmp(i + 1) + dx_tmp(i)))
+               vec_diag(i) = -(1.0 - alpha)* k_tmp(1) / (0.5*(dx_tmp(i))) &
+                             - (1.0 - alpha) * k_itf(i)/(0.5*(dx_tmp(i + 1) + dx_tmp(i))) &
+                             - rhocp_tmp(i) * dx_tmp(i) / dt_step
+               vec_rhs(i) = -rhocp_tmp(i) * dx_tmp(i) / dt_step * T_tmp(i) &
+                            - alpha * k_tmp(1)/(0.5 * dx_tmp(i)) * (T_up - T_tmp(i)) &
+                            + alpha * k_itf(i)/(0.5 * (dx_tmp(i + 1) + dx_tmp(i)))*(T_tmp(i) - T_tmp(i + 1)) &
+                            - (1.0 - alpha) * k_tmp(1) / (0.5*(dx_tmp(i)))*T_up
+            ELSE IF (i == n_tmp) THEN
+               vec_lw(i - 1) = (1.0 - alpha) * k_itf(i - 1)/(0.5*(dx_tmp(i - 1) + dx_tmp(i)))
+               vec_diag(i) = -(1.0 - alpha) * k_itf(i - 1)/(0.5*(dx_tmp(i - 1) + dx_tmp(i))) &
+                             - (1.0 - alpha) * k_tmp(n)/(0.5*(dx_tmp(i))) &
+                             - rhocp_tmp(i) * dx_tmp(i) / dt_step
+               vec_rhs(i) = -rhocp_tmp(i) * dx_tmp(i) / dt_step * T_tmp(i) &
+                            - alpha * k_itf(i - 1) / (0.5*(dx_tmp(i - 1) + dx_tmp(i)))*(T_tmp(i - 1) - T_tmp(i)) &
+                            + alpha * k_tmp(n) / (0.5*(dx_tmp(i)))*(T_tmp(i) - T_lw) &
+                            - (1.0 - alpha) * k_tmp(n) / (0.5*(dx_tmp(i))) * T_lw
+            ELSE
+               vec_lw(i - 1) = (1.0 - alpha) * k_itf(i - 1)/(0.5*(dx_tmp(i - 1) + dx_tmp(i)))
+               vec_up(i) = (1.0 - alpha) * k_itf(i)/(0.5*(dx_tmp(i + 1) + dx_tmp(i)))
+               vec_diag(i) = -(1.0 - alpha)*k_itf(i - 1)/(0.5*(dx_tmp(i - 1) + dx_tmp(i))) &
+                             - (1.0 - alpha)*k_itf(i)/(0.5*(dx_tmp(i + 1) + dx_tmp(i))) &
+                             - rhocp_tmp(i)*dx_tmp(i)/dt_step
+               vec_rhs(i) = -rhocp_tmp(i)*dx_tmp(i)/dt_step*T_tmp(i) &
+                            - alpha*k_itf(i - 1)/(0.5*(dx_tmp(i - 1) + dx_tmp(i)))*(T_tmp(i - 1) - T_tmp(i)) &
+                            + alpha*k_itf(i)/(0.5*(dx_tmp(i + 1) + dx_tmp(i)))*(T_tmp(i) - T_tmp(i + 1))
+            END IF
+         END DO
+
+         ! solve the tridiagonal matrix using Thomas algorithm
+         CALL thomas_triMat(vec_lw, vec_diag, vec_up, vec_rhs, n, T_tmp)
+         dt_remain = dt_remain - dt_step
+         Qs_acc = Qs_acc + (T_up - T_tmp(1)) * k_tmp(1) / (dx_tmp(1) * 0.5) * dt_step
+      END DO
+
+      DO i = 1, n
+         ids_new = (i - 1)*ratio + (ratio - 1) / 2 + 1
+         T_out(i) = T_tmp(ids_new)
+      END DO
+      Tsfc = T_out(1)
+      T = T_out
+
+      ! new way for calcualating heat storage
+      ! Qs = SUM( &
+      !      (([bc(1), T_out(1:n - 1)] + T_out)/2. & ! updated temperature
+      !       -([bc(1), T_in(1:n - 1)] + T_in)/2) & ! initial temperature
+      !      *rhocp*dx/dt)
+      ! Qs = SUM( &
+      !      (T_out - T_in) & ! initial temperature
+      !      *rhocp*dx/dt)
+      ! ---Here we use the outermost surface temperatures to calculate
+      ! ------the heat flux from the surface as the change of Qs for SEB
+      ! ------considering there might be fluxes going out from the lower boundary
+      Qs = (T_up - T_tmp(1)) * k_tmp(1)/(dx_tmp(1)*0.5)
+      ! Qs = (T_out(1) - T_out(2)) * k(1) / dx(1)
+      ! Qs = Qs_acc / dt
+      IF (debug) THEN
+         !PRINT *, "T_up: ", T_up, "T_lw: ", T_lw
+         !PRINT *, "T_out: ", T_out
+         !PRINT *, "T_in: ", T_in
+         !PRINT *, "Qs_last: ", Qs
+         !PRINT *, "Qs_acc_avg: ", Qs_acc / dt
+      END IF
+   END SUBROUTINE heatcond1d_CN_dense
+
 END MODULE heatflux
 
 MODULE EHC_module
@@ -251,7 +389,7 @@ SUBROUTINE EHC( &
       USE allocateArray, ONLY: &
          nsurf, ndepth, &
          PavSurf, BldgSurf, ConifSurf, DecidSurf, GrassSurf, BSoilSurf, WaterSurf
-      USE heatflux, ONLY: heatcond1d_vstep, heatcond1d_CN
+      USE heatflux, ONLY: heatcond1d_vstep, heatcond1d_CN, heatcond1d_CN_dense
 
       IMPLICIT NONE
       INTEGER, INTENT(in) :: tstep
@@ -448,13 +586,13 @@ SUBROUTINE EHC( &
                ! surface heat flux
                IF (i_group == 1) THEN
                   bc(1) = qg_roof(i_facet)
-                  !debug = .FALSE.
+                  debug = .True.
                ELSE IF (i_group == 2) THEN
                   bc(1) = qg_wall(i_facet)
-                  !debug = .FALSE.
+                  debug = .FALSE.
                ELSE IF (i_group == 3) THEN
                   bc(1) = QG_surf(i_facet)
-                  !debug = .True.
+                  debug = .FALSE.
                END IF
                ! bctype(1) = .TRUE.
                bc(1) = tsfc_cal(i_facet)
@@ -483,6 +621,7 @@ SUBROUTINE EHC( &
                ! CALL heatcond1d_ext( &
                !CALL heatcond1d_vstep( &
                CALL heatcond1d_CN( &
+               ! CALL heatcond1d_CN_dense( &
                   temp_cal(i_facet, :), &
                   QS_cal(i_facet), &
                   tsfc_cal(i_facet), &
@@ -496,16 +635,18 @@ SUBROUTINE EHC( &
                ! update temperature at all inner interfaces
                ! tin_cal(i_facet, :) = temp_all_cal
                ! IF (i_group == 3 .AND. i_facet == 3) THEN
-               ! PRINT *, i_facet, ' temp_cal after: ', temp_cal(i_facet, :)
+               IF (i_facet == 1) THEN
+                     PRINT *, i_facet, ' temp_cal after: ', temp_cal(i_facet, :)
+                     PRINT *, 'QS_cal after: ', QS_cal(i_facet)
                ! PRINT *, i_facet, 'tsfc_cal after: ', tsfc_cal(i_facet)
-               ! PRINT *, 'QS_cal after: ', QS_cal(i_facet)
+
                ! PRINT *, 'k_cal: ', k_cal(i_facet, 1:ndepth)
                ! PRINT *, 'cp_cal: ', cp_cal(i_facet, 1:ndepth)
                ! PRINT *, 'dz_cal: ', dz_cal(i_facet, 1:ndepth)
                ! PRINT *, 'bc: ', bc
                ! PRINT *, ''
 
-               ! END IF
+               END IF
             END IF
 
             IF (dz_cal(i_facet, 1) /= -999.0 .AND. use_heatcond1d_water) THEN
@@ -526,6 +667,7 @@ SUBROUTINE EHC( &
                ! CALL heatcond1d_ext( &
                !CALL heatcond1d_vstep( &
                CALL heatcond1d_CN( &
+               ! CALL heatcond1d_CN_dense( &
                   temp_cal(i_facet, :), &
                   QS_cal(i_facet), &
                   tsfc_cal(i_facet), &