Evasinha/ehc/degree days

cpl/gcam_cpl_indices.F90

@@ -43,6 +43,9 @@ module gcam_cpl_indices
   integer, pointer, public ::index_x2z_Sl_hr(:)        ! total heterotrophic respiration
   integer, pointer, public ::index_x2z_Sl_npp(:)       ! net primary production
   integer, pointer, public ::index_x2z_Sl_pftwgt(:)    ! pft weights for each cell
+  integer, public ::index_x2z_Sl_forc_hdm = 0          ! human population density
+  integer, public ::index_x2z_Sl_hdd = 0               ! heating degree days
+  integer, public ::index_x2z_Sl_cdd = 0               ! cooling degree days
   integer, public ::nflds_x2z = 0
 
   !-----------------------------------------------------------------------
@@ -78,8 +81,10 @@ subroutine gcam_cpl_indices_init( )
          gcam_config, base_gcam_co2_file, base_gcam_lu_wh_file, &
          base_co2_surface_file, base_co2_shipment_file, base_co2_aircraft_file, &
          base_npp_file, base_hr_file, base_pft_file, &
+         base_hdd_file, base_cdd_file, &
          gcam2elm_co2_mapping_file, gcam2elm_luc_mapping_file, &
          gcam2elm_woodharvest_mapping_file, gcam2elm_cdensity_mapping_file, &
+         gcam2elm_degdays_mapping_file, &
          gcam_gridfile, elm2gcam_mapping_file, &
          gcam2glm_glumap, gcam2glm_baselu, gcam2glm_basebiomass, &
          country2grid_map, country2region_map, pop_iiasa_file, gdp_iiasa_file, &
@@ -89,8 +94,8 @@ subroutine gcam_cpl_indices_init( )
          pasture_addtreeonly, pasture_setherbfracrem, pasture_setavailtreefracrem, &         
          fdyndat_ehc, &
          read_scalars, scalar_source_dir, &
-         write_scalars, write_co2, &
-         elm_ehc_agyield_scaling, elm_ehc_carbon_scaling, ehc_eam_co2_emissions,&
+         write_scalars, read_hdd_cdd, write_hdd_cdd, write_co2, &
+         elm_ehc_agyield_scaling, elm_ehc_carbon_scaling, elm_ehc_hdd_cdd, ehc_eam_co2_emissions,&
          gcam_spinup, run_gcam
 
     nlfilename_iac = "gcam_in"
@@ -213,6 +218,11 @@ subroutine gcam_cpl_indices_set( )
 
     end do
 
+    ! Scalar per-gridcell lnd->iac field (no PFT loop needed)
+    index_x2z_Sl_forc_hdm = mct_avect_indexra(x2z, 'Sl_forc_hdm')
+    index_x2z_Sl_hdd      = mct_avect_indexra(x2z, 'Sl_hdd')
+    index_x2z_Sl_cdd      = mct_avect_indexra(x2z, 'Sl_cdd')
+
     ! iac -> atm
     ! Monthly sfc, low alt air, high alt air
     do m=1,12

cpl/iac_import_export.F90

@@ -33,6 +33,8 @@ subroutine iac_import(x2z, lnd2iac_vars)
     integer :: n,n1,p,g,i,j
     integer :: begg, endg
     character(len=32), parameter :: sub = '(iac_import)'
+    real(r8), parameter ::  dayspy   = 365._r8                ! days per year
+    real(r8), parameter ::  k2f      = 9.0_r8/5.0_r8          ! Kelvin-day to Fahrenheit-day conversion
 
     ! Gcam expects things in npp_m[lon][lat][pft] format, so we need
     ! to extract from the flattened column representation.
@@ -71,6 +73,15 @@ subroutine iac_import(x2z, lnd2iac_vars)
        end do ! global index g
     end do ! pft index p
 
+    ! Scalar per-gridcell fields
+    do g=iac_ctl%begg,iac_ctl%endg
+       i=iac_ctl%ilon(g)
+       j=iac_ctl%jlat(g)
+       lnd2iac_vars%forc_hdm(i,j) = 0._r8; if (index_x2z_Sl_forc_hdm > 0) lnd2iac_vars%forc_hdm(i,j) = x2z(index_x2z_Sl_forc_hdm, g)
+       lnd2iac_vars%hdd(i,j)      = 0._r8; if (index_x2z_Sl_hdd      > 0) lnd2iac_vars%hdd(i,j)      = x2z(index_x2z_Sl_hdd,g) * dayspy * k2f
+       lnd2iac_vars%cdd(i,j)      = 0._r8; if (index_x2z_Sl_cdd      > 0) lnd2iac_vars%cdd(i,j)      = x2z(index_x2z_Sl_cdd,g) * dayspy * k2f
+    end do ! global index g
+
   end subroutine iac_import
    !===============================================================================
 

iac/coupling/gcam_comp_mod.F90

@@ -236,6 +236,8 @@ subroutine gcam_init_mod(gcamo, gcamoemis,gcamoco2sfcjan, gcamoco2sfcfeb,   &
     gcam2elm_woodharvest_mapping_file(char_len+1:char_len+1) = c_null_char
     char_len = len_trim(gcam2elm_cdensity_mapping_file)
     gcam2elm_cdensity_mapping_file(char_len+1:char_len+1) = c_null_char
+    char_len = len_trim(gcam2elm_degdays_mapping_file)
+    gcam2elm_degdays_mapping_file(char_len+1:char_len+1) = c_null_char
     char_len = len_trim(base_co2_surface_file)
     base_co2_surface_file(char_len+1:char_len+1) = c_null_char
     char_len = len_trim(base_co2_shipment_file)
@@ -281,6 +283,7 @@ subroutine gcam_init_mod(gcamo, gcamoemis,gcamoco2sfcjan, gcamoco2sfcfeb,   &
          trim(gcam2elm_luc_mapping_file),&
          trim(gcam2elm_woodharvest_mapping_file),&
          trim(gcam2elm_cdensity_mapping_file),&
+         trim(gcam2elm_degdays_mapping_file),&
          trim(base_gcam_co2_file), trim(base_co2_surface_file),&
          trim(base_co2_shipment_file), trim(base_co2_aircraft_file),&
          iac_ctl%area, iac_ctl%nlon, iac_ctl%nlat, num_gcam_energy_regions, num_emiss_sectors, rr)
@@ -360,12 +363,14 @@ subroutine gcam_run_mod(gcamo, gcamoemis, gcamoco2sfcjan, gcamoco2sfcfeb,   &
 
 ! !LOCAL VARIABLES:
     integer :: ymd, tod, dt
-    integer :: i,j,wc,gs,cs,rs,ws,rr,ays
+    integer :: i,j,wc,gs,cs,rs,ws,rdd,wdd,rr,ays,dd
     character(len=256) :: scalar_source_dir_loc
     character(len=256) :: elm2gcam_mapping_file_loc 
     character(len=256) :: base_npp_file_loc
     character(len=256) :: base_hr_file_loc
     character(len=256) :: base_pft_file_loc
+    character(len=256) :: base_hdd_file_loc
+    character(len=256) :: base_cdd_file_loc
     character(len=*),parameter :: subname='(gcam_run_mod)'
 
 
@@ -412,6 +417,18 @@ subroutine gcam_run_mod(gcamo, gcamoemis, gcamoco2sfcjan, gcamoco2sfcfeb,   &
      ws = 0
   end if
 
+  ! for reading/writing hdd and cdd from/to diagnostic files
+  if ( read_hdd_cdd ) then
+     rdd = 1
+  else
+     rdd = 0
+  end if
+  if ( write_hdd_cdd ) then
+     wdd = 1
+  else
+     wdd = 0
+  end if
+
   ! for ag yield scaling
   if ( elm_ehc_agyield_scaling ) then
      ays = 1
@@ -426,21 +443,30 @@ subroutine gcam_run_mod(gcamo, gcamoemis, gcamoco2sfcjan, gcamoco2sfcfeb,   &
      cs = 0
   end if
 
+  ! for passing hdd and cdd from elm to ehc
+  if ( elm_ehc_hdd_cdd ) then
+     dd = 1
+  else
+     dd = 0
+  end if
+
   ! get some file names for scalars
   ! use local variables to avoid adding multiple null characters to the orig
   scalar_source_dir_loc=trim(scalar_source_dir)//c_null_char
   elm2gcam_mapping_file_loc=trim(elm2gcam_mapping_file)//c_null_char
   base_npp_file_loc=trim(base_npp_file)//c_null_char
   base_hr_file_loc=trim(base_hr_file)//c_null_char
   base_pft_file_loc=trim(base_pft_file)//c_null_char
+  base_hdd_file_loc=trim(base_hdd_file)//c_null_char
+  base_cdd_file_loc=trim(base_cdd_file)//c_null_char
 
   !  Call runcGCAM method of E3SM Interface 
   !  The yields and carbon density scalars are set within this function also
   call runcGCAM(ymd, gcamo, gcamoemis, trim(base_gcam_lu_wh_file), trim(base_gcam_co2_file), gs, &
-                iac_ctl%area, lnd2iac_vars%pftwgt, lnd2iac_vars%npp, lnd2iac_vars%hr, &
-                iac_ctl%nlon, iac_ctl%nlat, iac_ctl%npft, num_gcam_energy_regions, num_emiss_ctys, num_emiss_sectors, num_periods,&
-                elm2gcam_mapping_file_loc, iac_first_coupled_year, rs, scalar_source_dir_loc, ws, ays, cs,&
-                base_npp_file_loc, base_hr_file_loc, base_pft_file_loc, rr)
+                iac_ctl%area, lnd2iac_vars%pftwgt, lnd2iac_vars%npp, lnd2iac_vars%hr, lnd2iac_vars%hdd, lnd2iac_vars%cdd, lnd2iac_vars%forc_hdm, &
+                iac_ctl%landfrac, iac_ctl%nlon, iac_ctl%nlat, iac_ctl%npft, num_gcam_energy_regions, num_emiss_ctys, num_emiss_sectors, num_periods,&
+                elm2gcam_mapping_file_loc, iac_first_coupled_year, rs, scalar_source_dir_loc, ws, rdd, wdd, ays, cs, dd, &
+                base_npp_file_loc, base_hr_file_loc, base_pft_file_loc, base_hdd_file_loc, base_cdd_file_loc, rr)
 
   ! If co2 emissions need to be passed from GCAM to EAM, then call downscale CO2                                 
   if ( ehc_eam_co2_emissions ) then

iac/coupling/gcam_var_mod.F90

@@ -127,10 +127,13 @@ module gcam_var_mod
   character(len=256), public ::  base_npp_file   ! by pft, lat, lon, for land scalars
   character(len=256), public ::  base_hr_file    ! by pft, lat, lon, for land scalars
   character(len=256), public ::  base_pft_file   ! by pft, lat, lon, for land scalars
+  character(len=256), public ::  base_hdd_file   ! by lat, lon, for baseline hdd
+  character(len=256), public ::  base_cdd_file   ! by lat, lon, for baseline cdd
   character(len=256), public ::  gcam2elm_co2_mapping_file ! define gcam co2 out
   character(len=256), public ::  gcam2elm_luc_mapping_file ! def gcam lu out
   character(len=256), public ::  gcam2elm_woodharvest_mapping_file ! def gcam wh out
   character(len=256), public ::  gcam2elm_cdensity_mapping_file ! def gcam cdensity in
+  character(len=256), public ::  gcam2elm_degdays_mapping_file ! def gcam degday in
 
   ! Grid and mapping and initialization files
   character(len=256), public ::  gcam_gridfile ! definition of iac grid
@@ -165,6 +168,8 @@ module gcam_var_mod
   character(len=256), public :: scalar_source_dir
   logical, public :: write_scalars ! scalars will be written to a file.
   ! hr, area, pft weight) are passed from e3sm.
+  logical, public :: read_hdd_cdd ! if .FALSE., HDD/CDD are calculated; if .TRUE. HDD/CDD are read from file
+  logical, public :: write_hdd_cdd ! hdd cdd will be written to a file.
   logical, public :: write_co2 ! gridded co2 emissions will be
   ! written to a file (in addition to passed in code).
 
@@ -174,6 +179,8 @@ module gcam_var_mod
   ! land productivity from elm scale ag yield in gcam.
   logical, public :: elm_ehc_carbon_scaling ! if .TRUE., changes in
   ! land productivity from elm scale carbon density in gcam.
+  logical, public :: elm_ehc_hdd_cdd ! if .TRUE., HDD and CDD
+  ! passed from elm to gcam.
   logical, public :: ehc_eam_co2_emissions ! if .TRUE., energy system
   ! co2 is passed from gcam to eam.
 

iac/coupling/iac_data_mod.F90

@@ -61,6 +61,9 @@ module iac_data_mod
      real(r8), allocatable :: npp(:,:,:)
      real(r8), allocatable :: hr(:,:,:)
      real(r8), allocatable :: pftwgt(:,:,:)
+     real(r8), allocatable :: forc_hdm(:,:)       ! population density (lon,lat)
+     real(r8), allocatable :: hdd(:,:)
+     real(r8), allocatable :: cdd(:,:)
   end type lnd2iac_type
 
   type, public :: iac2lnd_type

iac/coupling/iac_init_mod.F90

@@ -59,8 +59,10 @@ subroutine iac_init(EClock)
          gcam_config, base_gcam_co2_file, base_gcam_lu_wh_file, &
          base_co2_surface_file, base_co2_shipment_file, base_co2_aircraft_file, &
          base_npp_file, base_hr_file, base_pft_file, &
+         base_hdd_file, base_cdd_file, &
          gcam2elm_co2_mapping_file, gcam2elm_luc_mapping_file, &
          gcam2elm_woodharvest_mapping_file, gcam2elm_cdensity_mapping_file, &
+         gcam2elm_degdays_mapping_file, &
          gcam_gridfile, elm2gcam_mapping_file, &
          gcam2glm_glumap, gcam2glm_baselu, gcam2glm_basebiomass, &
          country2grid_map, country2region_map, pop_iiasa_file, gdp_iiasa_file, &
@@ -70,8 +72,8 @@ subroutine iac_init(EClock)
          pasture_addtreeonly, pasture_setherbfracrem, pasture_setavailtreefracrem, &         
          fdyndat_ehc, &
          read_scalars, scalar_source_dir, &
-         write_scalars, write_co2, &
-         elm_ehc_agyield_scaling, elm_ehc_carbon_scaling, ehc_eam_co2_emissions, &
+         write_scalars, read_hdd_cdd, write_hdd_cdd, write_co2, &
+         elm_ehc_agyield_scaling, elm_ehc_carbon_scaling, elm_ehc_hdd_cdd, ehc_eam_co2_emissions, &
          gcam_spinup, run_gcam 
 
     nlfilename_iac = "gcam_in"
@@ -135,10 +137,13 @@ subroutine iac_init(EClock)
        write(iulog, '(A,A)') "base_npp_file = ", trim(base_npp_file )
        write(iulog, '(A,A)') "base_hr_file = ", trim(base_hr_file)
        write(iulog, '(A,A)') "base_pft_file = ", trim(base_pft_file )
+       write(iulog, '(A,A)') "base_hdd_file = ", trim(base_hdd_file )
+       write(iulog, '(A,A)') "base_cdd_file = ", trim(base_cdd_file )
        write(iulog, '(A,A)') "gcam2elm_co2_mapping_file = ", trim(gcam2elm_co2_mapping_file )
        write(iulog, '(A,A)') "gcam2elm_luc_mapping_file = ", trim(gcam2elm_luc_mapping_file)
        write(iulog, '(A,A)') "gcam2elm_woodharvest_mapping_file = ", trim(gcam2elm_woodharvest_mapping_file)
        write(iulog, '(A,A)') "gcam2elm_cdensity_mapping_file = ", trim(gcam2elm_cdensity_mapping_file)
+       write(iulog, '(A,A)') "gcam2elm_degdays_mapping_file = ", trim(gcam2elm_degdays_mapping_file)
 
        write(iulog,*) 'grid mapping and initialization files:'
        write(iulog, '(A,A)') "gcam_gridfile = ", trim(gcam_gridfile)
@@ -173,12 +178,27 @@ subroutine iac_init(EClock)
        write(iulog, '(A,L)') "read_scalars = ",read_scalars
        write(iulog, '(A,A)') "scalar_source_dir = ", trim(scalar_source_dir)
        write(iulog, '(A,L10)') "write_scalars = ",write_scalars
+       write(iulog, '(A,L10)') "read_hdd_cdd = ",read_hdd_cdd
+       write(iulog, '(A,L10)') "write_hdd_cdd = ",write_hdd_cdd
        write(iulog, '(A,L10)') "write_co2 = ",write_co2
        write(iulog, '(A,L10)') "elm_ehc_agyield_scaling = ", elm_ehc_agyield_scaling
        write(iulog, '(A,L10)') "elm_ehc_carbon_scaling = ", elm_ehc_carbon_scaling
+       write(iulog, '(A,L10)') "elm_ehc_hdd_cdd = ", elm_ehc_hdd_cdd
        write(iulog, '(A,L10)') "ehc_eam_co2_emissions = ", ehc_eam_co2_emissions
        write(iulog, '(A,L10)') "gcam_spinup = ",gcam_spinup
        write(iulog, '(A,L10)') "run_gcam = ",run_gcam
+       if (read_scalars .and. write_scalars) then
+          write(iulog,*) '('//trim(subname)//') WARNING: read_scalars and write_scalars are both true.'
+          write(iulog,*) '('//trim(subname)//') read_scalars takes precedence, so write_scalars is effectively false:'
+          write(iulog,*) '('//trim(subname)//') no scalar calculations or writes will be performed.'
+          write_scalars = .false.
+       end if
+       if (read_hdd_cdd .and. write_hdd_cdd) then
+          write(iulog,*) '('//trim(subname)//') WARNING: read_hdd_cdd and write_hdd_cdd are both true.'
+          write(iulog,*) '('//trim(subname)//') read_hdd_cdd takes precedence, so write_hdd_cdd is effectively false:'
+          write(iulog,*) '('//trim(subname)//') no HDD/CDD calculations or writes will be performed.'
+          write_hdd_cdd = .false.
+       end if
 
        !if (nsrest == nsrStartup .and. finidat_rtm /= ' ') then
        !   write(iulog,*) '   MOSART initial data   = ',trim(finidat_rtm)
@@ -238,6 +258,9 @@ subroutine iac_init(EClock)
     allocate(lnd2iac_vars%npp(iac_ctl%nlon,iac_ctl%nlat,iac_ctl%npft))
     allocate(lnd2iac_vars%hr(iac_ctl%nlon,iac_ctl%nlat,iac_ctl%npft))
     allocate(lnd2iac_vars%pftwgt(iac_ctl%nlon,iac_ctl%nlat,iac_ctl%npft))
+    allocate(lnd2iac_vars%forc_hdm(iac_ctl%nlon,iac_ctl%nlat))
+    allocate(lnd2iac_vars%hdd(iac_ctl%nlon,iac_ctl%nlat))
+    allocate(lnd2iac_vars%cdd(iac_ctl%nlon,iac_ctl%nlat))
 
     allocate(iac2lnd_vars%pct_pft(iac_ctl%nlon,iac_ctl%nlat,iac_ctl%npft))
     allocate(iac2lnd_vars%pct_pft_prev(iac_ctl%nlon,iac_ctl%nlat,iac_ctl%npft))