Dynamical lakes and new lake cover map

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