Added lake landunit initialisation for gridcells which will grow lake

src/biogeophys/LakeTemperatureMod.F90

@@ -1000,9 +1000,10 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
 
        end do
     end do
-	write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts
+	! write(iulog,*)'Energy content of lake after calculating lake temperature (J/m²)', ncvts
 
-    lake_heat(c) = ncvts(c) 
+	! IV: currently commented out: caused crash. To do: look at this part of the code!!!
+    ! lake_heat(c) = ncvts(c) 
 
     do j = -nlevsno + 1, nlevgrnd
        do fc = 1, num_lakec

src/biogeophys/TotalWaterAndHeatMod.F90

@@ -748,7 +748,7 @@ subroutine ComputeHeatLake(bounds, num_lakec, filter_lakec, &
                latent_heat_liquid = latent_heat_liquid_lake(c))
     end do
 
-     write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c)
+    ! write(iulog,*) 'lake heat (J/m^2)', heat_lake(c)+latent_heat_liquid(c)
 
 ! Add lake heat here if wanted to incorporate
     do fc = 1, num_lakec

src/main/clm_initializeMod.F90

@@ -13,7 +13,7 @@ module clm_initializeMod
   use clm_varctl      , only : is_cold_start, is_interpolated_start
   use clm_varctl      , only : iulog
   use clm_varctl      , only : use_lch4, use_cn, use_cndv, use_c13, use_c14, use_fates
-  use clm_instur      , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft, wt_glc_mec, topo_glc_mec
+  use clm_instur      , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft, wt_glc_mec, topo_glc_mec, haslake
   use perf_mod        , only : t_startf, t_stopf
   use readParamsMod   , only : readParameters
   use ncdio_pio       , only : file_desc_t
@@ -161,7 +161,7 @@ subroutine initialize1( )
     allocate (fert_cft     (begg:endg, cft_lb:cft_ub       ))
     allocate (wt_glc_mec  (begg:endg, maxpatch_glcmec))
     allocate (topo_glc_mec(begg:endg, maxpatch_glcmec))
-
+    allocate (haslake      (begg:endg                      ))
     ! Read list of Patches and their corresponding parameter values
     ! Independent of model resolution, Needs to stay before surfrd_get_data
 
@@ -244,7 +244,7 @@ subroutine initialize1( )
     ! Some things are kept until the end of initialize2; urban_valid is kept through the
     ! end of the run for error checking.
 
-    deallocate (wt_lunit, wt_cft, wt_glc_mec)
+    deallocate (wt_lunit, wt_cft, wt_glc_mec, haslake)
 
     call t_stopf('clm_init1')
 
@@ -258,6 +258,7 @@ subroutine initialize2( )
     ! CLM initialization - second phase
     !
     ! !USES:
+
     use shr_orb_mod           , only : shr_orb_decl
     use shr_scam_mod          , only : shr_scam_getCloseLatLon
     use seq_drydep_mod        , only : n_drydep, drydep_method, DD_XLND

src/main/clm_varsur.F90

@@ -40,6 +40,9 @@ module clm_instur
 
   ! subgrid glacier_mec sfc elevation
   real(r8), pointer :: topo_glc_mec(:,:) 
+
+  ! whether we have lake to initialise in each grid cell
+  logical , pointer :: haslake(:)
   !-----------------------------------------------------------------------
 
 end module clm_instur

src/main/subgridMod.F90

@@ -39,7 +39,8 @@ module subgridMod
   public :: subgrid_get_info_glacier_mec
   public :: subgrid_get_info_crop
   public :: crop_patch_exists ! returns true if the given crop patch should be created in memory
-
+  public :: lake_landunit_exists ! returns true if the lake landunit should be created in memory
+
   ! !PRIVATE MEMBER FUNCTIONS:
   private :: subgrid_get_info_urban
 
@@ -393,10 +394,10 @@ subroutine subgrid_get_info_lake(gi, npatches, ncols, nlunits)
     character(len=*), parameter :: subname = 'subgrid_get_info_lake'
     !-----------------------------------------------------------------------
 
-    ! We currently do NOT allow the lake landunit to expand via dynamic landunits, so we
-    ! only need to allocate space for it where its weight is currently non-zero.
-
-    if (wt_lunit(gi, istdlak) > 0.0_r8) then
+    ! We do allow the lake landunit to expand via dynamic landunits, so we
+    !  need to allocate space for where it is known that the lake unit will grow.
+    
+    if (lake_landunit_exists(gi) ) then
        npatches = 1
        ncols = 1
        nlunits = 1
@@ -485,7 +486,6 @@ subroutine subgrid_get_info_crop(gi, npatches, ncols, nlunits)
     !-----------------------------------------------------------------------
 
     npatches = 0
-
     do cft = cft_lb, cft_ub
        if (crop_patch_exists(gi, cft)) then
           npatches = npatches + 1
@@ -559,6 +559,45 @@ function crop_patch_exists(gi, cft) result(exists)
 
   end function crop_patch_exists
 
+!-----------------------------------------------------------------------
+  function lake_landunit_exists(gi) result(exists)
+    !
+    ! !DESCRIPTION:
+    ! Returns true if a land unit for lakes should be created in memory
+    ! which is defined for gridcells which will grow lake, given by haslake
+    ! 
+    ! !USES:
+    use dynSubgridControlMod , only : get_do_transient_lakes
+    use clm_instur           , only : haslake
+    !
+    ! !ARGUMENTS:
+    logical :: exists  ! function result
+    integer, intent(in) :: gi  ! grid cell index
+    !
+    ! !LOCAL VARIABLES:
+
+    character(len=*), parameter :: subname = 'lake_landunit_exists'
+    !-----------------------------------------------------------------------
+
+    if (get_do_transient_lakes()) then
+       ! To support dynamic landunits, we initialise a lake land unit in each grid cell in which there are lakes. 
+       ! This is defined by the haslake variable
+
+       if (haslake(gi)) then
+            exists = .true.
+       else
+            exists = .false.
+       end if
+
+    else 
+        ! For a run without transient lakes, only allocate memory for lakes actually present in run)
+        if (wt_lunit(gi, istdlak) > 0.0_r8) then
+            exists = .true.
+        else
+            exists = .false.
+        end if
+    end if
 
+  end function lake_landunit_exists
 
 end module subgridMod

src/main/subgridWeightsMod.F90

@@ -301,7 +301,7 @@ logical function is_active_l(l, glc_behavior)
     ! Determine whether the given landunit is active
     !
     ! !USES:
-    use landunit_varcon, only : istsoil, istice_mec, isturb_MIN, isturb_MAX
+    use landunit_varcon, only : istsoil, istice_mec, isturb_MIN, isturb_MAX, istdlak
     !
     ! !ARGUMENTS:
     implicit none
@@ -361,7 +361,15 @@ logical function is_active_l(l, glc_behavior)
        if (lun%itype(l) == istsoil) then
           is_active_l = .true.
        end if
-
+
+       ! Set all lake land units to active
+       ! By doing this, lakes are also run virtually in grid cells which will grow
+       ! lakes during the transient run. 
+
+       if (lun%itype(l) == istdlak) then
+            is_active_l = .true.
+        end if
+
     end if
 
   end function is_active_l

src/main/surfrdMod.F90

@@ -284,10 +284,12 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
     !    o real % abundance PFTs (as a percent of vegetated area)
     !
     ! !USES:
-    use clm_varctl  , only : create_crop_landunit
-    use fileutils   , only : getfil
-    use domainMod   , only : domain_type, domain_init, domain_clean
-    use clm_instur  , only : wt_lunit, topo_glc_mec
+    use clm_varctl          , only : create_crop_landunit
+    use fileutils           , only : getfil
+    use domainMod           , only : domain_type, domain_init, domain_clean
+    use clm_instur          , only : wt_lunit, topo_glc_mec
+    use dynSubgridControlMod, only : get_flanduse_timeseries
+    use clm_varctl          , only : fname_len
     !
     ! !ARGUMENTS:
     integer,          intent(in) :: begg, endg      
@@ -304,8 +306,11 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
     logical           :: readvar              ! true => variable is on dataset
     real(r8)          :: rmaxlon,rmaxlat      ! local min/max vars
     type(file_desc_t) :: ncid                 ! netcdf id
+    type(file_desc_t) :: ncid_dynuse          ! netcdf id for landuse timeseries file
     logical           :: istype_domain        ! true => input file is of type domain
     logical           :: isgrid2d             ! true => intut grid is 2d 
+    character(len=fname_len)  :: fdynuse              ! landuse.timeseries filename
+
     character(len=32) :: subname = 'surfrd_get_data'    ! subroutine name
     !-----------------------------------------------------------------------
 
@@ -409,6 +414,33 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
        write(iulog,*) 'Successfully read surface boundary data'
        write(iulog,*)
     end if
+
+
+    ! IV: add here call to subroutine to read in lake mask (necessary for initialisation of dynamical lakes)
+    ! First open landuse.timeseries file for this. 
+
+     if (masterproc) then
+       write(iulog,*) 'Attempting to read landuse.timeseries data .....'
+       if (lfsurdat == ' ') then
+          write(iulog,*)'fdynuse must be specified'
+          call endrun(msg=errMsg(sourcefile, __LINE__))
+       endif
+    endif   
+
+
+    ! open landuse_timeseries file 
+    fdynuse = get_flanduse_timeseries()
+
+    call getfil(fdynuse, locfn, 0 )
+
+    call ncd_pio_openfile (ncid_dynuse, trim(locfn), 0)    
+
+    ! read the lakemask
+    call surfrd_lakemask(begg, endg, ncid_dynuse)
+
+    ! close landuse_timeseries file again 
+    call ncd_pio_closefile(ncid_dynuse)
+
 
   end subroutine surfrd_get_data
 
@@ -620,7 +652,7 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, cftsize, natpft_size )
     wt_nat_patch(begg:,natpft_lb:natpft_size-1+natpft_lb) = array2D(begg:,:)
     deallocate( array2D )
 
-  end subroutine surfrd_cftformat
+  end subroutine surfrd_cftformat      
 
 !-----------------------------------------------------------------------
   subroutine surfrd_pftformat( begg, endg, ncid )
@@ -795,4 +827,48 @@ subroutine surfrd_veg_dgvm(begg, endg)
 
   end subroutine surfrd_veg_dgvm
 
+  !-----------------------------------------------------------------------
+  subroutine surfrd_lakemask(begg, endg, ncid)
+    !
+    ! !DESCRIPTION:
+    ! Reads the lake mask, indicating where lakes are and will grow
+    ! of the landuse.timeseries file.
+    ! Necessary for the initialisation of the lake land units
+    !
+    ! !USES:
+     use clm_instur      , only : haslake
+    ! !ARGUMENTS:
+    integer,           intent(in)    :: begg, endg  
+    type(file_desc_t), intent(inout) :: ncid 
+    ! 
+    !
+    ! !LOCAL VARIABLES:
+    logical          :: readvar
+    integer ,pointer :: haslake_id(:)
+    !
+    !
+    character(len=*), parameter :: subname = 'surfrd_lakemask'
+    !
+    !-----------------------------------------------------------------------
+
+
+    allocate(haslake_id(begg:endg))
+
+    call ncd_io(ncid=ncid, varname='HASLAKE'  , flag='read', data=haslake_id, &
+           dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) call endrun( msg=' ERROR: HASLAKE is not on landuse.timeseries file'//errMsg(sourcefile, __LINE__))
+
+    where (haslake_id == 1.)
+        haslake = .true.
+    elsewhere
+        haslake = .false.
+    end where
+
+
+    deallocate(haslake_id)
+
+
+  end subroutine surfrd_lakemask
+
+
 end module surfrdMod