Dynamic lakes - without tools changes

bld/CLMBuildNamelist.pm

@@ -2372,6 +2372,7 @@ sub setup_logic_dynamic_subgrid {
 
    setup_logic_do_transient_pfts($opts, $nl_flags, $definition, $defaults, $nl);
    setup_logic_do_transient_crops($opts, $nl_flags, $definition, $defaults, $nl);
+   setup_logic_do_transient_lakes($opts, $nl_flags, $definition, $defaults, $nl);
    setup_logic_do_harvest($opts, $nl_flags, $definition, $defaults, $nl);
 
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'reset_dynbal_baselines');
@@ -2537,6 +2538,63 @@ sub setup_logic_do_transient_crops {
    }
 }
 
+sub setup_logic_do_transient_lakes {
+   #
+   # Set do_transient_lakes default value, and perform error checking on do_transient_lakes
+   #
+   # Assumes the following are already set in the namelist (although it's okay
+   # for them to be unset if that will be their final state):
+   # - flanduse_timeseries
+   #
+   my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
+
+   my $var = 'do_transient_lakes';
+
+   # cannot_be_true will be set to a non-empty string in any case where
+   # do_transient_lakes should not be true; if it turns out that
+   # do_transient_lakes IS true in any of these cases, a fatal error will be
+   # generated
+   my $cannot_be_true = "";
+
+   my $n_dom_pfts = $nl->get_value( 'n_dom_pfts' );
+   my $n_dom_landunits = $nl->get_value( 'n_dom_landunits' );
+   my $toosmall_soil = $nl->get_value( 'toosmall_soil' );
+   my $toosmall_crop = $nl->get_value( 'toosmall_crop' );
+   my $toosmall_glacier = $nl->get_value( 'toosmall_glacier' );
+   my $toosmall_lake = $nl->get_value( 'toosmall_lake' );
+   my $toosmall_wetland = $nl->get_value( 'toosmall_wetland' );
+   my $toosmall_urban = $nl->get_value( 'toosmall_urban' );
+
+   if (string_is_undef_or_empty($nl->get_value('flanduse_timeseries'))) {
+      $cannot_be_true = "$var can only be set to true when running a transient case (flanduse_timeseries non-blank)";
+   }
+
+   if (!$cannot_be_true) {
+      # Note that, if the variable cannot be true, we don't call add_default
+      # - so that we don't clutter up the namelist with variables that don't
+      # matter for this case
+      add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var);
+   }
+
+   # Make sure the value is false when it needs to be false - i.e., that the
+   # user hasn't tried to set a true value at an inappropriate time.
+
+   if (&value_is_true($nl->get_value($var)) && $cannot_be_true) {
+      $log->fatal_error($cannot_be_true);
+   }
+
+   # if do_transient_lakes is .true. and any of these (n_dom_* or toosmall_*)
+   # are > 0 or collapse_urban = .true., then give fatal error
+   if (&value_is_true($nl->get_value($var))) {
+      if (&value_is_true($nl->get_value('collapse_urban'))) {
+         $log->fatal_error("$var cannot be combined with collapse_urban");
+      }
+      if ($n_dom_pfts > 0 || $n_dom_landunits > 0 || $toosmall_soil > 0 || $toosmall_crop > 0 || $toosmall_glacier > 0 || $toosmall_lake > 0 || $toosmall_wetland > 0 || $toosmall_urban > 0) {
+         $log->fatal_error("$var cannot be combined with any of the of the following > 0: n_dom_pfts > 0, n_dom_landunit > 0, toosmall_soil > 0._r8, toosmall_crop > 0._r8, toosmall_glacier > 0._r8, toosmall_lake > 0._r8, toosmall_wetland > 0._r8, toosmall_urban > 0._r8");
+      }
+   }
+}
+
 sub setup_logic_do_harvest {
    #
    # Set do_harvest default value, and perform error checking on do_harvest

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -2965,6 +2965,7 @@ lnd/clm2/surfdata_map/landuse.timeseries_ne30np4_hist_16pfts_Irrig_CMIP6_simyr18
 <!-- (some other defaults for dynamic subgrid   -->
 <!-- are set in the BuildNamelist code)         -->
 <!-- =========================================  -->
+<do_transient_lakes>.false.</do_transient_lakes>
 <reset_dynbal_baselines>.false.</reset_dynbal_baselines>
 
 <!-- =========================================  -->

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2426,6 +2426,12 @@ If TRUE, apply transient crops from flanduse_timeseries file.
 (Only valid for transient runs, where there is a flanduse_timeseries file.)
 </entry>
 
+<entry id="do_transient_lakes" type="logical" category="physics"
+       group="dynamic_subgrid" valid_values="" >
+If TRUE, apply transient lakes from flanduse_timeseries file.
+(Only valid for transient runs, where there is a flanduse_timeseries file.)
+</entry>
+
 <entry id="do_harvest" type="logical" category="physics"
        group="dynamic_subgrid" valid_values="" >
 If TRUE, apply harvest from flanduse_timeseries file.

src/biogeophys/BalanceCheckMod.F90

@@ -20,6 +20,7 @@ module BalanceCheckMod
   use SoilHydrologyType  , only : soilhydrology_type
   use SurfaceAlbedoType  , only : surfalb_type
   use WaterStateType     , only : waterstate_type
+  use LakestateType      , only : lakestate_type
   use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type
   use WaterDiagnosticType, only : waterdiagnostic_type
   use Wateratm2lndType   , only : wateratm2lnd_type
@@ -122,7 +123,7 @@ end function GetBalanceCheckSkipSteps
   !-----------------------------------------------------------------------
   subroutine BeginWaterBalance(bounds, &
        num_nolakec, filter_nolakec, num_lakec, filter_lakec, &
-       water_inst, soilhydrology_inst, &
+       water_inst, soilhydrology_inst, lakestate_inst, &
        use_aquifer_layer)
     !
     ! !DESCRIPTION:
@@ -136,6 +137,7 @@ subroutine BeginWaterBalance(bounds, &
     integer                   , intent(in)    :: num_lakec            ! number of column lake points in column filter
     integer                   , intent(in)    :: filter_lakec(:)      ! column filter for lake points
     type(water_type)          , intent(inout) :: water_inst
+    type(lakestate_type)      , intent(in)    :: lakestate_inst
     type(soilhydrology_type)  , intent(in)    :: soilhydrology_inst
     logical                   , intent(in)    :: use_aquifer_layer    ! whether an aquifer layer is used in this run
     !
@@ -150,6 +152,7 @@ subroutine BeginWaterBalance(bounds, &
             num_nolakec, filter_nolakec, &
             num_lakec, filter_lakec, &
             soilhydrology_inst, &
+            lakestate_inst, &
             water_inst%bulk_and_tracers(i)%waterstate_inst, &
             water_inst%bulk_and_tracers(i)%waterdiagnostic_inst, &
             water_inst%bulk_and_tracers(i)%waterbalance_inst, &
@@ -161,7 +164,8 @@ end subroutine BeginWaterBalance
   !-----------------------------------------------------------------------
   subroutine BeginWaterBalanceSingle(bounds, &
        num_nolakec, filter_nolakec, num_lakec, filter_lakec, &
-       soilhydrology_inst, waterstate_inst, waterdiagnostic_inst, waterbalance_inst, &
+       soilhydrology_inst, lakestate_inst, waterstate_inst, & 
+       waterdiagnostic_inst, waterbalance_inst, &
        use_aquifer_layer)
     !
     ! !DESCRIPTION:
@@ -175,6 +179,7 @@ subroutine BeginWaterBalanceSingle(bounds, &
     integer                   , intent(in)    :: num_lakec            ! number of column lake points in column filter
     integer                   , intent(in)    :: filter_lakec(:)      ! column filter for lake points
     type(soilhydrology_type)  , intent(in)    :: soilhydrology_inst
+    type(lakestate_type)      , intent(in)    :: lakestate_inst
     class(waterstate_type)    , intent(inout) :: waterstate_inst
     class(waterdiagnostic_type), intent(in)   :: waterdiagnostic_inst
     class(waterbalance_type)  , intent(inout) :: waterbalance_inst
@@ -210,7 +215,7 @@ subroutine BeginWaterBalanceSingle(bounds, &
          water_mass = begwb(bounds%begc:bounds%endc))
 
     call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
-         waterstate_inst, &
+         waterstate_inst, lakestate_inst, &
          subtract_dynbal_baselines = .false., &
          water_mass = begwb(bounds%begc:bounds%endc))
 

src/biogeophys/LakeHydrologyMod.F90

@@ -164,7 +164,7 @@ subroutine LakeHydrology(bounds, &
          t_soisno             =>  temperature_inst%t_soisno_col         , & ! Output: [real(r8) (:,:) ]  snow temperature (Kelvin)             
          dTdz_top             =>  temperature_inst%dTdz_top_col         , & ! Output: [real(r8) (:)   ]  temperature gradient in top layer K m-1] !TOD 
          snot_top             =>  temperature_inst%snot_top_col         , & ! Output: [real(r8) (:)   ]  snow temperature in top layer [K]  !TODO
-         t_sno_mul_mss        => temperature_inst%t_sno_mul_mss_col     , & ! Output: [real(r8) (:)   ]  col snow temperature multiplied by layer mass, layer sum (K * kg/m2) 
+         t_sno_mul_mss        =>  temperature_inst%t_sno_mul_mss_col     , & ! Output: [real(r8) (:)   ]  col snow temperature multiplied by layer mass, layer sum (K * kg/m2) 
 
          begwb                =>  b_waterbalance_inst%begwb_col             , & ! Input:  [real(r8) (:)   ]  water mass begining of the time step    
          endwb                =>  b_waterbalance_inst%endwb_col             , & ! Output: [real(r8) (:)   ]  water mass end of the time step         
@@ -649,7 +649,7 @@ subroutine LakeHydrology(bounds, &
     ! Determine ending water balance and volumetric soil water
 
     call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, &
-         b_waterstate_inst, &
+         b_waterstate_inst, lakestate_inst, &
          subtract_dynbal_baselines = .false., &
          water_mass = endwb(bounds%begc:bounds%endc))
 

src/biogeophys/LakeTemperatureMod.F90

@@ -134,6 +134,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
     type(energyflux_type)  , intent(inout) :: energyflux_inst
     type(temperature_type) , intent(inout) :: temperature_inst
     type(lakestate_type)   , intent(inout) :: lakestate_inst
+
     !
     ! !LOCAL VARIABLES:
     real(r8), parameter :: p0 = 1._r8                                      ! neutral value of turbulent prandtl number
@@ -247,7 +248,8 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
          t_grnd          =>   temperature_inst%t_grnd_col          , & ! Input:  [real(r8) (:)   ]  ground temperature (Kelvin)             
          t_soisno        =>   temperature_inst%t_soisno_col        , & ! Output: [real(r8) (:,:) ]  soil (or snow) temperature (Kelvin)   
          t_lake          =>   temperature_inst%t_lake_col          , & ! Output: [real(r8) (:,:) ]  col lake temperature (Kelvin)             
-
+         lake_heat       =>   temperature_inst%lake_heat           , & ! Output: [real(r8) (:) ]    col lake heat (J/m²) 
+
          beta            =>   lakestate_inst%betaprime_col         , & ! Output: [real(r8) (:)   ]  col effective beta: sabg_lyr(p,jtop) for snow layers, beta otherwise
          lake_icefrac    =>   lakestate_inst%lake_icefrac_col      , & ! Output: [real(r8) (:,:) ]  col mass fraction of lake layer that is frozen
          lake_icefracsurf => lakestate_inst%lake_icefracsurf_col   , & ! Output: [real(r8) (:,:) ]  col mass fraction of surface lake layer that is frozen
@@ -999,8 +1001,13 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
           ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
                    + cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) 
           fin(c) = fin(c) + phi(c,j)
+
        end do
     end do
+	! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts
+
+	! IV: currently commented out: caused crash. To do: look at this part of the code!!!
+    ! lake_heat(c) = ncvts(c) 
 
     call waterstatebulk_inst%CalculateTotalH2osno(bounds, num_lakec, filter_lakec, &
          caller = 'LakeTemperature-2', &
@@ -1022,6 +1029,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
     end do
 
 
+
     ! Check energy conservation.
     do fp = 1, num_lakep
        p = filter_lakep(fp)

src/biogeophys/TemperatureType.F90

@@ -117,6 +117,9 @@ module TemperatureType
      real(r8), pointer    :: fact_col              (:,:) ! used in computing tridiagonal matrix
      real(r8), pointer    :: c_h2osfc_col          (:)   ! heat capacity of surface water
 
+     ! lake heat
+     real(r8), pointer    :: lake_heat              (:)  ! total heat of lake water (J/m²)
+
    contains
 
      procedure, public  :: Init
@@ -279,6 +282,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%fact_col                 (begc:endc, -nlevsno+1:nlevgrnd)) ; this%fact_col                 (:,:) = nan
     allocate(this%c_h2osfc_col             (begc:endc))                      ; this%c_h2osfc_col             (:)   = nan
 
+    allocate(this%lake_heat                (begc:endc))                      ; this%lake_heat                (:)   = nan
   end subroutine InitAllocate
 
   !------------------------------------------------------------------------
@@ -620,6 +624,11 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp )
             ptr_patch=this%t_veg10_night_patch, default='inactive')
     endif
 
+    ! add lake heat history field here
+    this%lake_heat(begc:endc) = spval
+    call hist_addfld1d (fname='LAKE_HEAT', units='J/m^2', &
+            avgflag='A', long_name='Heat content of gridcell lake water', &
+            ptr_col=this%lake_heat, default='active')
 
   end subroutine InitHistory