Fix crop phenology bugs with Gregorian calendar

src/biogeochem/CNPhenologyMod.F90

@@ -61,6 +61,7 @@ module CNPhenologyMod
   public :: SeasonalCriticalDaylength ! Critical day length needed for Seasonal decidious offset
   public :: get_swindow
   public :: was_sown_in_this_window
+  public :: PlantDate_to_PlantJday
 
   ! !PRIVITE MEMBER FIUNCTIONS:
   private :: CNPhenologyClimate             ! Get climatological everages to figure out triggers for Phenology
@@ -1792,7 +1793,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
     ! handle CN fluxes during the phenological onset                       & offset periods.
 
     ! !USES:
-    use clm_time_manager , only : get_prev_calday, get_curr_days_per_year, is_beg_curr_year
+    use clm_time_manager , only : get_prev_calday, get_prev_days_per_year, is_beg_curr_year
     use clm_time_manager , only : get_average_days_per_year
     use clm_time_manager , only : get_prev_date
     use clm_time_manager , only : is_doy_in_interval, is_end_curr_day
@@ -1921,11 +1922,15 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
          )
 
       ! get time info
-      dayspyr = get_curr_days_per_year()
+      dayspyr = get_prev_days_per_year()
       avg_dayspyr = get_average_days_per_year()
       jday    = get_prev_calday()
       call get_prev_date(kyr, kmo, kda, mcsec)
 
+      if (is_beg_curr_year()) then
+         call UpdateSowingWindows()
+      end if
+
       if (use_fertilizer) then
        ndays_on = 20._r8 ! number of days to fertilize
       else
@@ -2210,7 +2215,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
             end if
 
             ! days past planting may determine harvest
-            idpp = DaysPastPlanting(idop(p), jday)
+            idpp = DaysPastPlanting(idop(p))
 
             ! onset_counter initialized to zero when .not. croplive
             ! offset_counter relevant only at time step of harvest
@@ -2477,6 +2482,53 @@ subroutine CropPhase(bounds, num_pcropp, filter_pcropp, &
   end subroutine CropPhase
 
 
+  function PlantDate_to_PlantJday(plantdate) result(jday)
+    !
+    ! !DESCRIPTION:
+    ! Converts a plantdate from parameter file (e.g., mn[NS]Hplantdate) to a jday
+    ! that's actually used in the model (e.g., minplantjday).
+    !
+    ! !USES:
+    use clm_time_manager, only: get_calday, get_prev_date
+    ! !ARGUMENTS
+    integer, intent(in) :: plantdate
+    !
+    ! !LOCAL VARIABLES
+    integer :: kyr, kmo, kda, mcsec
+    !
+    ! Return value
+    integer :: jday
+
+    ! Get year to add to plantdate
+    call get_prev_date(kyr, kmo, kda, mcsec)
+
+    jday = int( get_calday(10000*kyr + plantdate, 0 ) )
+  end function PlantDate_to_PlantJday
+
+
+  subroutine UpdateSowingWindows()
+    !
+    ! !DESCRIPTION:
+    ! Updates sowing windows with dates for this year.
+    !
+    ! !USES:
+    use pftconMod, only: npcropmin, npcropmax
+    !
+    ! !LOCAL VARIABLES
+    integer :: n
+
+    do n = npcropmin, npcropmax
+      if (pftcon%is_pft_known_to_model(n)) then
+         minplantjday(n, inNH) = PlantDate_to_PlantJday(pftcon%mnNHplantdate(n))
+         maxplantjday(n, inNH) = PlantDate_to_PlantJday(pftcon%mxNHplantdate(n))
+
+         minplantjday(n, inSH) = PlantDate_to_PlantJday(pftcon%mnSHplantdate(n))
+         maxplantjday(n, inSH) = PlantDate_to_PlantJday(pftcon%mxSHplantdate(n))
+      end if
+   end do
+  end subroutine UpdateSowingWindows
+
+
   !-----------------------------------------------------------------------
   subroutine CropPhenologyInit(bounds)
     !
@@ -2507,15 +2559,7 @@ subroutine CropPhenologyInit(bounds)
     ! Convert planting dates into julian day
     minplantjday(:,:) = huge(1)
     maxplantjday(:,:) = huge(1)
-    do n = npcropmin, npcropmax
-       if (pftcon%is_pft_known_to_model(n)) then
-          minplantjday(n, inNH) = int( get_calday( pftcon%mnNHplantdate(n), 0 ) )
-          maxplantjday(n, inNH) = int( get_calday( pftcon%mxNHplantdate(n), 0 ) )
-
-          minplantjday(n, inSH) = int( get_calday( pftcon%mnSHplantdate(n), 0 ) )
-          maxplantjday(n, inSH) = int( get_calday( pftcon%mxSHplantdate(n), 0 ) )
-       end if
-    end do
+    call UpdateSowingWindows()
 
     ! Figure out what hemisphere each PATCH is in
     do p = bounds%begp, bounds%endp
@@ -2723,33 +2767,27 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, &
   end subroutine PlantCrop
 
   !-----------------------------------------------------------------------
-  function DaysPastPlanting(idop, jday_in)
+  function DaysPastPlanting(idop)
     ! !USES:
     use clm_time_manager, only : get_prev_calday, get_curr_days_per_year
+    use clm_varcon      , only : secspday
     !
     ! !ARGUMENTS:
     integer,           intent(in) :: idop ! patch day of planting
-    integer, optional, intent(in) :: jday_in ! julian day of the year
     !
     ! !LOCAL VARIABLES
     integer :: DaysPastPlanting
     integer :: jday
 
-    ! Must use separate jday_in and jday because we can't redefine an intent(in)
-    ! variable, even if it wasn't provided in the function call.
-    if (present(jday_in)) then
-       jday = jday_in
-    else
-       ! Use prev instead of curr to avoid jday=1 in last timestep of year
-       jday = get_prev_calday()
-    end if
+    ! Use prev instead of curr to avoid jday=1 in last timestep of year
+    jday = get_prev_calday()
 
     if (jday >= idop) then
        DaysPastPlanting = jday - idop
     else
-      ! As long as crops have at most a 365-day growing season, using get_curr_days_per_year()
-      ! should give the same result of this function as using get_prev_days_per_year().
-       DaysPastPlanting = jday - idop + get_curr_days_per_year()
+       ! get_curr_days_per_year() or get_prev_days_per_year() would only differ in the last timestep
+       ! of the year, but in that case this line is not reached.
+       DaysPastPlanting = jday - idop + get_curr_days_per_year(offset = -365*int(secspday))
     end if
 
   end function DaysPastPlanting

src/biogeochem/test/CNPhenology_test/CMakeLists.txt

@@ -1,5 +1,6 @@
 set (pfunit_sources
-  test_CNPhenology.pf)
+  test_CNPhenology.pf
+  test_CNPhenology_gregorian.pf)
 
 add_pfunit_ctest(CNPhenology
   TEST_SOURCES "${pfunit_sources}"

src/biogeochem/test/CNPhenology_test/test_CNPhenology.pf

@@ -501,4 +501,118 @@ contains
 
   end subroutine test_was_sown_in_this_window_sameday
 
+  @Test
+  subroutine test_DaysPastPlanting_mar25_feb9(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! Feb. 9 in a non-leap year
+    call unittest_timemgr_set_curr_date(1, 2, 9, this%dtime)
+
+    idop = 84  ! Mar. 25
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(321, idpp)
+
+  end subroutine test_DaysPastPlanting_mar25_feb9
+
+  @Test
+  subroutine test_DaysPastPlanting_feb9_mar25(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! Mar. 25 in a non-leap year
+    call unittest_timemgr_set_curr_date(1, 3, 25, this%dtime)
+
+    idop = 40  ! Feb. 9
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(44, idpp)
+
+  end subroutine test_DaysPastPlanting_feb9_mar25
+
+  @Test
+  subroutine test_DaysPastPlanting_feb9_EODdec31(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! Last timestep of Dec. 31 in a non-leap year
+    call unittest_timemgr_set_curr_date(3, 1, 1, 0)
+
+    idop = 40  ! Feb. 9
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(325, idpp)
+
+  end subroutine test_DaysPastPlanting_feb9_EODdec31
+
+  @Test
+  subroutine test_DaysPastPlanting_feb9_BODjan1(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! First timestep of Jan. 1 in a non-leap year following
+    ! another non-leap year
+    call unittest_timemgr_set_curr_date(3, 1, 1, this%dtime)
+
+    idop = 40  ! Feb. 9
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(326, idpp)
+
+  end subroutine test_DaysPastPlanting_feb9_BODjan1
+
+  @Test
+  subroutine test_DaysPastPlanting_feb9_mar25_leap(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! Mar. 25 in a leap year, although that doesn't matter in
+    ! this test module because there are no leap days
+    call unittest_timemgr_set_curr_date(4, 3, 25, this%dtime)
+
+    idop = 40  ! Feb. 9
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(44, idpp)
+
+  end subroutine test_DaysPastPlanting_feb9_mar25_leap
+
+  @Test
+  subroutine test_DaysPastPlanting_feb9_jan2_lastYrLeap(this)
+    class(TestCNPhenology), intent(inout) :: this
+    integer :: idop, idpp
+
+    ! Jan. 2 in the year AFTER a leap year, although that doesn't
+    ! matter in this test module because there are no leap days
+    call unittest_timemgr_set_curr_date(5, 1, 2, this%dtime)
+
+    idop = 40  ! Feb. 9
+    idpp = DaysPastPlanting(idop)
+
+    @assertEqual(327, idpp)
+
+  end subroutine test_DaysPastPlanting_feb9_jan2_lastYrLeap
+
+  @Test
+  subroutine test_PlantDate_to_PlantJday_jan1(this)
+    class(TestCNPhenology), intent(inout) :: this
+
+    @assertEqual(1, PlantDate_to_PlantJday(101))
+  end subroutine test_PlantDate_to_PlantJday_jan1
+
+  @Test
+  subroutine test_PlantDate_to_PlantJday_mar25(this)
+    class(TestCNPhenology), intent(inout) :: this
+
+    @assertEqual(84, PlantDate_to_PlantJday(325))
+  end subroutine test_PlantDate_to_PlantJday_mar25
+
+  @Test
+  subroutine test_PlantDate_to_PlantJday_dec31(this)
+    class(TestCNPhenology), intent(inout) :: this
+
+    @assertEqual(365, PlantDate_to_PlantJday(1231))
+  end subroutine test_PlantDate_to_PlantJday_dec31
+
 end module test_CNPhenology