Abc lhf work (#23)

Author: skfogg

autotest/test_gwe_sfe_abc_lhf00.py

@@ -58,7 +58,7 @@
 ]
 # sensible and latent heat flux parameter values
 wspd = 20.0
-tatm = 270.87556216  # unrealistically high to drive a -1C change in stream temperature
+tatm = 270.79111 # unrealistically high to drive a -1C change in stream temperature
 # shortwave radiation parameter values
 solr = 47880870.9  # unrealistically high to drive a 1 deg C rise in stream temperature
 shd = 1.0  # 100% shade "turns off" solar flux
@@ -90,7 +90,7 @@
 perlen = [1]
 
 nouter, ninner = 1000, 300
-hclose, rclose, relax = 1e-3, 1e-4, 0.97
+hclose, rclose, relax = 1e-10, 1e-10, 0.97
 
 #
 # MODFLOW 6 flopy GWF object
@@ -455,11 +455,12 @@ def check_output(idx, test):
 
     # confirm that the energy added to the stream results in a -1C change in temp
     # temperature gradient
+    
     tgrad = tatm - strm_temp
     shf_ener_per_sqm = c_d * rhoa * Cpa * wspd * tgrad
     swr_ener_per_sqm = solr * (1 - shd) * (1 - swrefl)
     # latent calcs
-    L = (2499.64 - 2.51 * strm_temp) * 1000
+    L = (2499.64 - (2.51 * strm_temp)) * 1000
     e_w = 6.1275 * math.exp(17.2693882 * (strm_temp / (strm_temp + 273.16 - 35.86)))
     e_s = 6.1275 * math.exp(17.2693882 * (tatm / (tatm + 273.16 - 35.86)))
     e_a = rh / 100 * e_s
@@ -473,22 +474,25 @@ def check_output(idx, test):
     )
     # calculate expected temperature change based on energy transfer
     temp_change = ener_transfer / (surf_Q_in[idx][0] * Cpw * rhow)
+    
+    print(str(strm_temp + temp_change))
+
+    #new_strm_temp = strm_temp + temp_change
 
     fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv")
     assert os.path.isfile(fpth2)
     df2 = pd.read_csv(fpth2)
 
     # confirm 1 deg C decrease in temp
-    msg1 = (
-        "The calculated temperature change \
-        is: "
-        + str(temp_change)
+    msg1 = "Python temperature change is = " + str(temp_change)
+    msg2 = "MODFLOW temperature = " + str(df2.loc[0, "RCH1_OUTFTEMP"])
+    msg3 = "MODFLOW temperature change is " + str(
+        strm_temp - df2.loc[0, "RCH1_OUTFTEMP"]
     )
-    msg2 = "The temperature is:" + str(df2.loc[0, "RCH1_OUTFTEMP"])
 
     assert np.isclose(
         df2.loc[0, "RCH1_OUTFTEMP"], strm_temp + temp_change, atol=1e-6
-    ), msg2
+    ), msg2 + ". " + msg3 + ". " + msg1
 
 
 # - No need to change any code below

make/makefile

@@ -425,6 +425,8 @@ $(OBJDIR)/tsp-obs.o \
 $(OBJDIR)/tsp-mvt.o \
 $(OBJDIR)/tsp-adv.o \
 $(OBJDIR)/gwf-uzf.o \
+$(OBJDIR)/tsp-apt_weird_chngs.o \
+$(OBJDIR)/gwt-mst_updates.o \
 $(OBJDIR)/GwtDspOptions.o \
 $(OBJDIR)/gwf-tvs.o \
 $(OBJDIR)/GwfStorageUtils.o \

src/Model/GroundWaterEnergy/gwe-lhf.f90

@@ -139,16 +139,20 @@ subroutine lhf_cq(this, ifno, tstrm, rhow, lhflx)
     real(DP) :: sat_vap_ta
     real(DP) :: amb_vap_atm
     real(DP) :: evap_rate
+    real(DP) :: test_mult
     !
     ! -- calculate latent heat flux using mass transfer equation
     !! EQ BASED ON TEMPERATURES IN CELSIUS
     ! -- temperature dependent latent heat of vaporization:
-    latent_heat_vap = (2499.64 - (2.51*tstrm))*1000 ! tstrm in C
+    !test_mult = 2.51*tstrm
 
-    sat_vap_tw = 6.1275*exp(17.2693882*((tstrm)/(tstrm + DCTOK - 35.86)))
-    sat_vap_ta = 6.1275*exp(17.2693882*((this%tatm(ifno))/(this%tatm(ifno) + DCTOK - 35.86)))
+    latent_heat_vap = (2499.64_DP - (2.51_DP*tstrm))*1000.0_DP ! tstrm in C
+    !latent_heat_vap = 2449440.0
 
-    amb_vap_atm = (this%rh(ifno)/100)*sat_vap_ta
+    sat_vap_tw = 6.1275_DP*exp(17.2693882_DP*((tstrm)/(tstrm + DCTOK - 35.86_DP)))
+    sat_vap_ta = 6.1275_DP*exp(17.2693882_DP*((this%tatm(ifno))/(this%tatm(ifno) + DCTOK - 35.86_DP)))
+    
+    amb_vap_atm = (this%rh(ifno)/100.0_DP)*sat_vap_ta
 
     evap_rate = (this%wfint + this%wfslope*this%wspd(ifno))*(sat_vap_tw - amb_vap_atm)