add swr to options block of sfe (#16)

Author: skfogg

autotest/test_gwe_sfe_swr00.py

@@ -55,7 +55,7 @@
     [10.0],
 ]
 # shortwave radiation parameter values
-solr = 4788.08709  # unrealistically high to drive a 1 deg C rise in stream temperature
+solr = 47880870.9  # unrealistically high to drive a 1 deg C rise in stream temperature
 shd = 0.1
 swrefl = 0.03
 
@@ -435,7 +435,7 @@ def check_output(idx, test):
     df = pd.read_csv(fpth)
     calc_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy()
     # confirm stream width is 1.0 m
-    assert np.isclose(calc_strm_wid, 1.0, atol=1e-9), msg0
+    assert np.isclose(calc_strm_wid, 1.0, atol=1e-7), msg0
 
     # confirm that the energy added to the stream results in a 1 deg C rise in temp
     # temperature gradient
@@ -458,7 +458,7 @@ def check_output(idx, test):
     )
 
     assert np.isclose(df.loc[0, "RCH1_OUTFTEMP"], strm_temp + temp_rise, atol=1e-5), (
-        msg1
+        msg1 + "  " + str(df.loc[0, "RCH1_OUTFTEMP"])
     )
 
 

doc/mf6io/mf6ivar/dfn/gwe-sfe.dfn

@@ -286,6 +286,38 @@ tagged false
 longname file name of SHF information
 description defines a sensible heat flux (SHF) input file.  Records in the SHF file are used to calculate the flux of thermal energy in or out of a stream reach resulting from sensible heat exchange.
 
+block options
+name swr_filerecord
+type record swr6 filein swr6_filename
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
+block options
+name swr6
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname swr keyword
+description keyword to specify that record corresponds to the shortwave radiation (SWR) file.  The behavior of the SWR utility and a description of the input file is provided separately.
+
+block options
+name swr6_filename
+type string
+preserve_case true
+in_record true
+reader urword
+optional false
+tagged false
+longname file name of SWR information
+description defines a shortwave radiation (SWR) input file.  Records in the SWR file are used to calculate the flux of thermal energy in or out of a stream reach resulting from shortwave radiation heat flux.
+
 # --------------------- gwe sfe packagedata ---------------------
 
 block packagedata

src/Model/GroundWaterEnergy/gwe-sfe.f90

@@ -23,6 +23,7 @@
 ! -- SFE terms
 ! STRMBD-COND               idxbudsbcd    STRMBD-COND           ktf * wa / sbthk * (t_cell - t_feat)
 ! SENSIBLE HEAT FLUX        idxbudshf     SENS HEAT             cd * rho_a * C_p_a * wspd * (t_air - t_feat)
+! SHORTWAVE RADIATION       idxbudsWR     SHORTWAVE             (1 - shd) * (1 - swrefl) * solr
 
 ! -- terms from a flow file that should be skipped
 ! CONSTANT                  none          none                  none
@@ -49,6 +50,7 @@ module GweSfeModule
                           apt_process_obsID12
   use GweInputDataModule, only: GweInputDataType
   use SensHeatModule, only: ShfType, shf_cr
+  use ShortwaveModule, only: SwrType, swr_cr
   use MatrixBaseModule
   use InputOutputModule, only: openfile
   !
@@ -72,6 +74,7 @@ module GweSfeModule
     integer(I4B), pointer :: idxbudoutf => null() !< index of outflow terms in flowbudptr
 
     logical, pointer, public :: shf_active => null() !< logical indicating if a sensible heat flux object is active
+    logical, pointer, public :: swr_active => null() !< logical indicating if a shortwave radition heat flux object is active
 
     real(DP), dimension(:), pointer, contiguous :: temprain => null() !< rainfall temperature
     real(DP), dimension(:), pointer, contiguous :: tempevap => null() !< evaporation temperature
@@ -81,7 +84,10 @@ module GweSfeModule
     real(DP), dimension(:), pointer, contiguous :: rfeatthk => null() !< thickness of streambed material through which thermal conduction occurs
 
     type(ShfType), pointer :: shf => null() ! sensible heat flux (shf) object
+    type(SwrType), pointer :: swr => null() ! shortwave radiation heat flux (swr) object
+
     integer(I4B), pointer :: inshf => null() ! SHF (sensible heat flux utility) unit number (0 if unused)
+    integer(I4B), pointer :: inswr => null() ! SWR (shortwave radiation heat flux utility) unit number (0 if unused)
 
   contains
 
@@ -104,6 +110,7 @@ module GweSfeModule
     procedure :: sfe_outf_term
     procedure, private :: sfe_sbcd_term
     procedure, private :: sfe_shf_term
+    procedure, private :: sfe_swr_term
     procedure :: pak_df_obs => sfe_df_obs
     procedure :: pak_rp_obs => sfe_rp_obs
     procedure :: pak_bd_obs => sfe_bd_obs
@@ -231,6 +238,37 @@ subroutine sfe_options(this, option, found)
       !          would be from bnd_ar(), however, shf is only callable from
       !          sfe and not apt.
       call this%shf%ar()
+    case ('SWR6')
+      !
+      call this%parser%GetStringCaps(keyword)
+      if (trim(adjustl(keyword)) /= 'FILEIN') then
+        errmsg = 'SWR6 keyword must be followed by "FILEIN" '// &
+                 'then by filename.'
+        call store_error(errmsg)
+        call this%parser%StoreErrorUnit()
+      end if
+      if (this%swr_active) then
+        errmsg = 'Multiple SWR6 keywords detected in OPTIONS block. '// &
+                 'Only one SWR6 entry allowed for a package.'
+        call store_error(errmsg)
+      end if
+      this%swr_active = .true.
+      call this%parser%GetString(fname)
+      !
+      ! -- create sensible heat flux object
+      call openfile(this%inswr, this%iout, fname, 'SWR')
+      call swr_cr(this%swr, this%name_model, this%inswr, this%iout, this%ncv)
+      this%swr%inputFilename = fname
+      !
+      ! -- call _ar routine for swr sub-package
+      !    note: this is the best place to call the sub-package (swr) for
+      !          now because SFE does not override apt_ar(), and since sfe
+      !          does not run its own customized set of _ar() functionality
+      !          there is no opportunity to run swr%ar() from a more
+      !          natural point.  In other words, one place to call swr%ar
+      !          would be from bnd_ar(), however, swr is only callable from
+      !          sfe and not apt.
+      call this%swr%ar()
     case default
       !
       ! -- No options found
@@ -465,6 +503,16 @@ subroutine sfe_fc_expanded(this, rhs, ia, idxglo, matrix_sln)
         rhs(iloc) = rhs(iloc) + rhsval
       end do
     end if
+    ! -- Add shortwave radiation heat flux contribution
+    if (this%inswr /= 0) then
+      do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist
+        call this%sfe_swr_term(j, n1, n2, rrate, rhsval, hcofval)
+        iloc = this%idxlocnode(n1)
+        iposd = this%idxpakdiag(n1)
+        call matrix_sln%add_value_pos(iposd, hcofval)
+        rhs(iloc) = rhs(iloc) + rhsval
+      end do
+    end if
   end subroutine sfe_fc_expanded
 
   !> @ brief Add terms specific to sfr to the explicit sfe solve
@@ -547,6 +595,9 @@ function sfe_get_nbudterms(this) result(nbudterms)
     if (this%inshf /= 0) then
       nbudterms = nbudterms + 1
     end if
+    if (this%inswr /= 0) then
+      nbudterms = nbudterms + 1
+    end if
     !
   end function sfe_get_nbudterms
 
@@ -670,6 +721,27 @@ subroutine sfe_setup_budobj(this, idx)
         call this%budobj%budterm(idx)%update_term(n1, n2, q)
       end do
     end if
+    ! -- Shortwave radiation heat flux
+    if (this%inswr /= 0) then
+      text = '   SHORTWAVE-RADIATION'
+      idx = idx + 1
+      maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist
+      naux = 0
+      call this%budobj%budterm(idx)%initialize(text, &
+                                               this%name_model, &
+                                               this%packName, &
+                                               this%name_model, &
+                                               this%packName, &
+                                               maxlist, .false., .false., &
+                                               naux)
+      call this%budobj%budterm(idx)%reset(maxlist)
+      q = DZERO
+      do n = 1, maxlist
+        n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(n)
+        n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n)
+        call this%budobj%budterm(idx)%update_term(n1, n2, q)
+      end do
+    end if
   end subroutine sfe_setup_budobj
 
   !> @brief Copy flow terms into this%budobj
@@ -771,6 +843,17 @@ subroutine sfe_fill_budobj(this, idx, x, flowja, ccratin, ccratout)
         call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
       end do
     end if
+    ! -- Shortwave-radiation
+    if (this%inswr /= 0) then
+      idx = idx + 1
+      nlist = this%flowbudptr%budterm(this%idxbudgwf)%nlist
+      call this%budobj%budterm(idx)%reset(nlist)
+      do j = 1, nlist
+        call this%sfe_swr_term(j, n1, n2, q)
+        call this%budobj%budterm(idx)%update_term(n1, n2, q)
+        call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout)
+      end do
+    end if
   end subroutine sfe_fill_budobj
 
   !> @brief Allocate scalars specific to the streamflow energy transport (SFE)
@@ -792,7 +875,9 @@ subroutine allocate_scalars(this)
     call mem_allocate(this%idxbudiflw, 'IDXBUDIFLW', this%memoryPath)
     call mem_allocate(this%idxbudoutf, 'IDXBUDOUTF', this%memoryPath)
     call mem_allocate(this%shf_active, 'SHF_ACTIVE', this%memoryPath)
+    call mem_allocate(this%swr_active, 'SWR_ACTIVE', this%memoryPath)
     call mem_allocate(this%inshf, 'INSHF', this%memoryPath)
+    call mem_allocate(this%inswr, 'INSWR', this%memoryPath)
     !
     ! -- Initialize
     this%idxbudrain = 0
@@ -802,7 +887,9 @@ subroutine allocate_scalars(this)
     this%idxbudoutf = 0
     !
     this%shf_active = .false.
+    this%swr_active = .false.
     this%inshf = 0
+    this%inswr = 0
   end subroutine allocate_scalars
 
   !> @brief Allocate arrays specific to the streamflow energy transport (SFE)
@@ -837,6 +924,9 @@ subroutine sfe_allocate_arrays(this)
     if (this%inshf /= 0) then
       call this%shf%pbst_allocate_arrays()
     end if
+    if (this%inswr /= 0) then
+      call this%swr%pbst_allocate_arrays()
+    end if
   end subroutine sfe_allocate_arrays
 
   !> @brief Call Read and prepare routines for any active pbst subpackages
@@ -853,6 +943,10 @@ subroutine pbst_rp(this)
     if (this%inshf /= 0) then
       call this%shf%rp()
     end if
+    ! -- call shortwave radiation heat flux sub-package _rp() routine
+    if (this%inswr /= 0) then
+      call this%swr%rp()
+    end if
   end subroutine pbst_rp
 
   !> @brief Deallocate memory
@@ -868,6 +962,11 @@ subroutine sfe_da(this)
       call this%shf%da()
       deallocate (this%shf)
     end if
+    ! -- SWR (shortwave radiation heat flux)
+    if (this%inswr /= 0) then
+      call this%swr%da()
+      deallocate (this%swr)
+    end if
     !
     ! -- Deallocate scalars
     call mem_deallocate(this%idxbudrain)
@@ -877,7 +976,9 @@ subroutine sfe_da(this)
     call mem_deallocate(this%idxbudoutf)
     !
     call mem_deallocate(this%shf_active)
+    call mem_deallocate(this%swr_active)
     call mem_deallocate(this%inshf)
+    call mem_deallocate(this%inswr)
     !
     ! -- Deallocate time series
     call mem_deallocate(this%temprain)
@@ -1090,6 +1191,36 @@ subroutine sfe_shf_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
     if (present(hcofval)) hcofval = DZERO
   end subroutine sfe_shf_term
 
+  !> @brief Shortwave Radiation (SWR) term
+  !<
+  subroutine sfe_swr_term(this, ientry, n1, n2, rrate, rhsval, hcofval)
+    ! -- dummy
+    class(GweSfeType) :: this
+    integer(I4B), intent(in) :: ientry
+    integer(I4B), intent(inout) :: n1
+    integer(I4B), intent(inout) :: n2
+    real(DP), intent(inout), optional :: rrate
+    real(DP), intent(inout), optional :: rhsval
+    real(DP), intent(inout), optional :: hcofval
+    ! -- local
+    real(DP) :: shrtwvheat
+    real(DP) :: strmtemp
+    integer(I4B) :: auxpos
+    real(DP) :: sa !< surface area of stream reach, different than wetted area
+    !
+    n1 = this%flowbudptr%budterm(this%idxbudevap)%id1(ientry)
+    ! -- For now, there is only 1 aux variable under 'EVAPORATION'
+    auxpos = this%flowbudptr%budterm(this%idxbudevap)%naux
+    sa = this%flowbudptr%budterm(this%idxbudevap)%auxvar(auxpos, ientry)
+    !
+    strmtemp = this%xnewpak(n1)
+    call this%swr%swr_cq(n1, shrtwvheat)
+    !
+    if (present(rrate)) rrate = shrtwvheat * sa
+    if (present(rhsval)) rhsval = -rrate
+    if (present(hcofval)) hcofval = DZERO
+  end subroutine sfe_swr_term
+
   !> @brief Observations
   !!
   !! Store the observation type supported by the APT package and override
@@ -1171,6 +1302,11 @@ subroutine sfe_df_obs(this)
     !    for sens-heat-flux observation type.
     call this%obs%StoreObsType('shf', .true., indx)
     this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID
+    !
+    ! -- Store obs type and assign procedure pointer
+    !    for shortwave-radiation-flux observation type.
+    call this%obs%StoreObsType('swr', .true., indx)
+    this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID
   end subroutine sfe_df_obs
 
   !> @brief Process package specific obs
@@ -1202,6 +1338,8 @@ subroutine sfe_rp_obs(this, obsrv, found)
       call this%rp_obs_byfeature(obsrv)
     case ('SHF')
       call this%rp_obs_byfeature(obsrv)
+    case ('SWR')
+      call this%rp_obs_byfeature(obsrv)
     case default
       found = .false.
     end select
@@ -1249,6 +1387,10 @@ subroutine sfe_bd_obs(this, obstypeid, jj, v, found)
       if (this%iboundpak(jj) /= 0) then
         call this%sfe_shf_term(jj, n1, n2, v)
       end if
+    case ('SWR')
+      if (this%iboundpak(jj) /= 0) then
+        call this%sfe_swr_term(jj, n1, n2, v)
+      end if
     case default
       found = .false.
     end select