Skip to content

Commit

Permalink
BUGFIX: GDD methods now working and verified
Browse files Browse the repository at this point in the history
  • Loading branch information
smwesten-usgs committed Jun 10, 2020
1 parent 11ae49e commit 55709e8
Show file tree
Hide file tree
Showing 10 changed files with 15,563 additions and 100 deletions.
30 changes: 22 additions & 8 deletions src/growing_degree_day.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ module growing_degree_day

public :: GDD_BASE, GDD_MAX, GDD_RESET_DATE
public :: growing_degree_day_calculate, growing_degree_day_initialize
public :: modified_growing_degree_day_calculate

real (c_float), allocatable :: GDD_BASE(:)
real (c_float), allocatable :: GDD_MAX(:)
Expand Down Expand Up @@ -140,25 +141,38 @@ impure elemental subroutine growing_degree_day_calculate( gdd, tmean, order )
real (c_float), intent(in) :: tmean
integer (c_int), intent(in) :: order

! [ LOCALS ]
real (c_float) :: delta
real (c_float) :: bounded_mean_air_temp

associate( doy_to_reset_gdd => GDD_RESET_DATE( order ), &
gdd_max => GDD_MAX( order ), &
gdd_base => GDD_BASE( order ) )

if ( SIM_DT%iDOY == doy_to_reset_gdd ) gdd = 0.0_c_float

bounded_mean_air_temp = enforce_bounds(value=tmean, minval=gdd_base, maxval=gdd_max)

delta = bounded_mean_air_temp - gdd_base
gdd = gdd + delta
gdd = gdd + max(tmean, gdd_base) - gdd_base

end associate

end subroutine growing_degree_day_calculate

!--------------------------------------------------------------------------------------------------

impure elemental subroutine modified_growing_degree_day_calculate( gdd, tmin, tmax, order )

! [ ARGUMENTS ]
real (c_float), intent(inout) :: gdd
real (c_float), intent(in) :: tmin
real (c_float), intent(in) :: tmax
integer (c_int), intent(in) :: order

associate( doy_to_reset_gdd => GDD_RESET_DATE( order ), &
gdd_max => GDD_MAX( order ), &
gdd_base => GDD_BASE( order ) )

if ( SIM_DT%iDOY == doy_to_reset_gdd ) gdd = 0.0_c_float

gdd = gdd + max((max(tmin, gdd_base) + min(tmax,gdd_max))/2.0_c_float, gdd_base) - gdd_base

end associate

end subroutine modified_growing_degree_day_calculate

end module growing_degree_day
37 changes: 18 additions & 19 deletions src/growing_degree_day_baskerville_emin.F90
Original file line number Diff line number Diff line change
Expand Up @@ -136,14 +136,14 @@ end subroutine growing_degree_day_be_initialize

!--------------------------------------------------------------------------------------------------

elemental subroutine growing_degree_day_be_calculate( gdd, tmean, tmin, tmax, order )
subroutine growing_degree_day_be_calculate( gdd, tmean, tmin, tmax, order )

! [ ARGUMENTS ]
real (c_float), intent(inout) :: gdd
real (c_float), intent(in) :: tmean
real (c_float), intent(in) :: tmin
real (c_float), intent(in) :: tmax
integer (c_int), intent(in) :: order
real (c_float), intent(inout) :: gdd(:)
real (c_float), intent(in) :: tmean(:)
real (c_float), intent(in) :: tmin(:)
real (c_float), intent(in) :: tmax(:)
integer (c_int), intent(in) :: order(:)

! [ LOCALS ]
real (c_float) :: delta
Expand All @@ -152,42 +152,41 @@ elemental subroutine growing_degree_day_be_calculate( gdd, tmean, tmin, tmax, or
real (c_float) :: W
real (c_float) :: At
real (c_float) :: A
integer (c_int) :: indx

associate( doy_to_reset_gdd => GDD_RESET_DATE( order ), &
gdd_max => GDD_MAX( order ), &
gdd_base => GDD_BASE( order ) )
do indx=lbound(order,1),ubound(order,1)

if ( SIM_DT%iDOY == doy_to_reset_gdd ) gdd = 0.0_c_float
if ( SIM_DT%iDOY == GDD_RESET_DATE( order(indx) ) ) gdd(indx) = 0.0_c_float

tmax_l = min( tmax, gdd_max )
tmax_l = min( tmax(indx), GDD_MAX( order(indx) ) )

if ( tmax_l <= gdd_base ) then
if ( tmax_l <= GDD_BASE( order(indx) ) ) then

dd = 0.0_c_float

elseif ( tmin >= gdd_base ) then
elseif ( tmin(indx) >= GDD_BASE( order(indx) ) ) then

dd = min(tmean, gdd_max - gdd_base )
dd = min(( tmean(indx) - GDD_BASE( order(indx) )), (GDD_MAX( order(indx) ) - GDD_BASE( order(indx) ) ) )

else

W = ( tmax_l - tmin) / 2.0_c_float
W = ( tmax_l - tmin(indx)) / 2.0_c_float

At = ( gdd_base - tmean ) / W
At = ( GDD_BASE( order(indx) ) - tmean(indx) ) / W

if ( At > 1 ) At = 1.0_c_float
if ( At < -1 ) At = -1._c_float

A = asin( At )

dd = ( ( W * cos( A ) ) - ( ( gdd_base - tmean ) &
dd = ( ( W * cos( A ) ) - ( ( GDD_BASE( order(indx) ) - tmean(indx) ) &
* ( real( PI / 2._c_double, c_float ) - A ) ) ) / PI

endif

gdd = gdd + dd
gdd(indx) = gdd(indx) + dd

end associate
enddo

end subroutine growing_degree_day_be_calculate

Expand Down
59 changes: 45 additions & 14 deletions src/model_domain.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1251,22 +1251,38 @@ subroutine set_method_pointers_sub(this, sCmdText, argv_list )
.or. ( Method_Name .strapprox. "BE" ) &
.or. ( Method_Name .strapprox. "SINUSOIDAL" ) ) then

! this%init_GDD => model_initialize_GDD_be
! this%calc_GDD => model_calculate_GDD_be
! call LOGS%WRITE( "==> Growing degree-day (GDD) will be calculated as described " &
! //"in Baskerville and Emin (1969)", iLogLevel = LOG_ALL, lEcho = FALSE )
this%init_GDD => model_initialize_GDD_be
this%calc_GDD => model_calculate_GDD_be
call LOGS%WRITE( "==> Growing degree-day (GDD) will be calculated as described " &
//"in Baskerville and Emin (1969)", iLogLevel = LOG_ALL, lEcho = FALSE )
! this%init_GDD => model_initialize_GDD
! this%calc_GDD => model_calculate_GDD
! call LOGS%WRITE( "==> THIS OPTION IS TEMPORARILY DISABLED! Growing degree-day (GDD) will be calculated using " &
! //"simple averaging of TMAX and TMIN.", iLogLevel = LOG_ALL, lEcho = FALSE )

elseif( ( Method_Name .strapprox. "SIMPLE" ) &
.or. ( Method_Name .strapprox. "SIMPLE_GROWING_DEGREE_DAY" ) &
.or. ( Method_Name .strapprox. "SIMPLE_GROWING_DEGREE-DAY") ) then

this%init_GDD => model_initialize_GDD
this%calc_GDD => model_calculate_GDD
call LOGS%WRITE( "==> THIS OPTION IS TEMPORARILY DISABLED! Growing degree-day (GDD) will be calculated using " &
call LOGS%WRITE( "==> Growing degree-day (GDD) will be calculated using " &
//"simple averaging of TMAX and TMIN.", iLogLevel = LOG_ALL, lEcho = FALSE )

else
elseif( ( Method_Name .strapprox. "MODIFIED" ) &
.or. ( Method_Name .strapprox. "MODIFIED_GROWING_DEGREE-DAY" ) &
.or. ( Method_Name .strapprox. "MODIFIED_GROWING_DEGREE_DAY" ) ) then

this%init_GDD => model_initialize_GDD
this%calc_GDD => model_calculate_GDD
call LOGS%WRITE( "==> Growing degree-day (GDD) will be calculated using " &
//"simple averaging of TMAX and TMIN.", iLogLevel = LOG_ALL, lEcho = FALSE )
this%calc_GDD => model_calculate_modified_GDD
call LOGS%WRITE( "==> Modified growing degree-day (GDD) will be calculated using " &
//"a simple averaging of TMAX and TMIN.", iLogLevel = LOG_ALL, lEcho = FALSE )

else

call warn("Your control file specifies an unknown or unsupported GROWING DEGREE-DAY method.", &
lFatal = TRUE, iLogLevel = LOG_ALL, lEcho = TRUE )

endif

elseif ( sCmdText .containssimilar. "IRRIGATION" ) then
Expand Down Expand Up @@ -2654,7 +2670,22 @@ subroutine model_calculate_GDD( this )

end subroutine model_calculate_GDD

!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------

subroutine model_calculate_modified_GDD( this )

use growing_degree_day, only : modified_growing_degree_day_calculate

class (MODEL_DOMAIN_T), intent(inout) :: this

call modified_growing_degree_day_calculate( gdd=this%gdd, &
tmin=this%tmin, &
tmax=this%tmax, &
order=this%sort_order )

end subroutine model_calculate_modified_GDD

!--------------------------------------------------------------------------------------------------

subroutine model_initialize_GDD_be( this )

Expand Down Expand Up @@ -2684,10 +2715,10 @@ subroutine model_calculate_GDD_be( this )
class (MODEL_DOMAIN_T), intent(inout) :: this

call growing_degree_day_be_calculate( gdd=this%gdd, &
tmean=this%tmean, &
tmin=this%tmin, &
tmax=this%tmax, &
order=this%sort_order )
tmean=this%tmean, &
tmin=this%tmin, &
tmax=this%tmax, &
order=this%sort_order )

end subroutine model_calculate_GDD_be

Expand Down
127 changes: 127 additions & 0 deletions test/integration_tests/cs/central_sands_swb2.ctl
Original file line number Diff line number Diff line change
@@ -0,0 +1,127 @@
%% Central Sands, Wisconsin
%% Example of net infiltration calculation with crop water demand
%% and irrigation included in water budget

! Grid definition: projected coordinates are Wisconsin Transverse Mercator (83/91), meters
! nx ny xll yll resolution
GRID 400 346 545300 432200 90.0
BASE_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m

%% Define methods
-----------------

INTERCEPTION_METHOD BUCKET
EVAPOTRANSPIRATION_METHOD HARGREAVES
RUNOFF_METHOD CURVE_NUMBER
SOIL_MOISTURE_METHOD FAO-56_TWO_STAGE
PRECIPITATION_METHOD GRIDDED
GROWING_DEGREE_DAY_METHOD MODIFIED_GROWING_DEGREE_DAY
FOG_METHOD NONE
FLOW_ROUTING_METHOD NONE
IRRIGATION_METHOD FAO-56
ROOTING_DEPTH_METHOD DYNAMIC
CROP_COEFFICIENT_METHOD FAO-56
DIRECT_RECHARGE_METHOD NONE
SOIL_STORAGE_MAX_METHOD CALCULATED
AVAILABLE_WATER_CONTENT_METHOD GRIDDED

GROWING_SEASON 133 268 TRUE


%% define location, projection, and conversions for weather data
----------------------------------------------------------------

PRECIPITATION NETCDF prcp_Daymet_v3_%y.nc
PRECIPITATION_GRID_PROJECTION_DEFINITION +proj=lcc +lat_1=25.0 +lat_2=60.0 +lat_0=42.5 +lon_0=-100.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
PRECIPITATION_NETCDF_Z_VAR prcp
PRECIPITATION_SCALE_FACTOR 0.03937008
PRECIPITATION_MISSING_VALUES_CODE -9999.0
PRECIPITATION_MISSING_VALUES_OPERATOR <=
PRECIPITATION_MISSING_VALUES_ACTION zero

TMAX NETCDF tmax_Daymet_v3_%y.nc
TMAX_GRID_PROJECTION_DEFINITION +proj=lcc +lat_1=25.0 +lat_2=60.0 +lat_0=42.5 +lon_0=-100.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
TMAX_SCALE_FACTOR 1.8
TMAX_ADD_OFFSET 32.0
TMAX_MISSING_VALUES_CODE -9999.0
TMAX_MISSING_VALUES_OPERATOR <=
TMAX_MISSING_VALUES_ACTION mean

TMIN NETCDF tmin_Daymet_v3_%y.nc
TMIN_GRID_PROJECTION_DEFINITION +proj=lcc +lat_1=25.0 +lat_2=60.0 +lat_0=42.5 +lon_0=-100.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
TMIN_SCALE_FACTOR 1.8
TMIN_ADD_OFFSET 32.0
TMIN_MISSING_VALUES_CODE -9999.0
TMIN_MISSING_VALUES_OPERATOR <=
TMIN_MISSING_VALUES_ACTION mean

INITIAL_CONTINUOUS_FROZEN_GROUND_INDEX CONSTANT 100.0
UPPER_LIMIT_CFGI 83.
LOWER_LIMIT_CFGI 55.

%% specify location and projection for input GIS grids
------------------------------------------------------

FLOW_DIRECTION ARC_GRID d8_flow_direction.asc
FLOW_DIRECTION_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m

HYDROLOGIC_SOILS_GROUP ARC_GRID hydrologic_soils_group.asc
HYDROLOGIC_SOILS_GROUP_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m

LAND_USE ARC_GRID landuse.asc
LANDUSE_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m

AVAILABLE_WATER_CONTENT ARC_GRID available_water_capacity.asc
AVAILABLE_WATER_CONTENT_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m

IRRIGATION_MASK ARC_GRID irrigation_mask_from_cdl.asc
IRRIGATION_MASK_PROJECTION_DEFINITION +proj=tmerc +lat_0=0.0 +lon_0=-90.0 +k=0.9996 +x_0=520000 +y_0=-4480000 +datum=NAD83 +units=m


%% specify location and names for all lookup tables
---------------------------------------------------

LAND_USE_LOOKUP_TABLE Landuse_lookup_CDL.txt
IRRIGATION_LOOKUP_TABLE Irrigation_lookup_CDL.txt

%% initial conditions for soil moisture and snow storage amounts
%% may be specified as grids, but using a constant amount and
%% allowing the model to "spin up" for a year is also acceptable.

INITIAL_PERCENT_SOIL_MOISTURE CONSTANT 100.0
INITIAL_SNOW_COVER_STORAGE CONSTANT 2.0

# POTATOES
DUMP_VARIABLES COORDINATES 561167. 445224.

# RYE
DUMP_VARIABLES COORDINATES 546094., 438492.

# EVERGREEN FOREST
DUMP_VARIABLES COORDINATES 556791. 457569.

# WINTER WHEAT
DUMP_VARIABLES COORDINATES 568129. 458340.

# DEC FOREST
DUMP_VARIABLES COORDINATES 553927. 459454.

# SWEET CORN
DUMP_VARIABLES COORDINATES 555602. 434644.

# SOYBEANS
DUMP_VARIABLES COORDINATES 558663. 432949.

# SUGAR BEETS
DUMP_VARIABLES COORDINATES 570741. 445112.


OUTPUT ENABLE bare_soil_evaporation crop_et snowmelt soil_storage delta_soil_storage
OUTPUT ENABLE growing_season growing_degree_day

%% start and end date may be any valid dates in SWB version 2.0
%% remember to allow for adequate model spin up; running the
%% model for just a month or two will give questionable results

START_DATE 01/01/2012
END_DATE 12/31/2013
Loading

0 comments on commit 55709e8

Please sign in to comment.