Merge remote-tracking branch 'pengbinpeluo/agsys' into agsys

src/agsys/ctsm_interface/AgSysInterface.F90

@@ -156,7 +156,7 @@ subroutine AgSysDriver(this, num_pcropp, filter_pcropp, &
                      tt_in_phase       = this%agsys_inst%acc_thermal_time_in_phase_patch(p,:), &
                      days_after_phase  = this%agsys_inst%days_after_phase_patch(p,:), &
                      tt_after_phase    = this%agsys_inst%acc_thermal_time_after_phase_patch(p,:), &
-                     phase_target_tt   = this%agsys_inst%phase_target_thermal_time_path(p,:), &
+                     phase_target_tt   = this%agsys_inst%phase_target_thermal_time_patch(p,:), &
                      cumvd             = this%agsys_inst%acc_vernalization_days_patch(p))
              end if
 

src/agsys/ctsm_interface/AgSysParamReader.F90

@@ -77,11 +77,21 @@ subroutine ReadParams(crops)
           cultivar%rc_photoperiod_target_tt%y=[0._r8, 0._r8, 0._r8]
           cultivar%rc_photoperiod_target_tt%num_pts=3
         class is (agsys_crop_type_maize)
-          cultivar%tt_emerg_to_endjuv = 275._r8
-          cultivar%tt_flag_to_flower  = 1._r8
-          cultivar%tt_flower_to_maturity = 812._r8
-          cultivar%tt_flower_to_start_grain = 170._r8
-          cultivar%tt_maturity_to_ripe = 1._r8
+          cultivar%target_tt_emerg_to_endjuv = 275._r8
+          cultivar%target_tt_flag_to_flower  = 1._r8
+          cultivar%target_tt_flower_to_maturity = 812._r8
+          cultivar%target_tt_flower_to_start_grain = 170._r8
+          cultivar%target_tt_maturity_to_ripe = 1._r8
+
+          cultivar%leaf_init_rate = 23.2_r8
+          cultivar%leaf_no_seed = 6._r8
+          cultivar%leaf_no_at_emerg = 0.5_r8
+          cultivar%leaf_no_min = 5.0_r8
+          cultivar%leaf_no_max = 40._r8
+          cultivar%leaf_no_critical = 11._r8
+          cultivar%leaf_appearance_rate_early = 65._r8
+          cultivar%leaf_appearance_rate_late = 36._r8
+
           cultivar%potential_kernel_weight = 300._r8
           cultivar%leaf_no_dead_const = -0.025_r8
           cultivar%leaf_no_dead_slope = 0.00035_r8      

src/agsys/ctsm_interface/AgSysType.F90

@@ -83,7 +83,7 @@ module AgSysType
      real(r8), pointer, public :: acc_emerged_thermal_time_patch(:)  ! accumulated thermal time since emergence (deg-days)
      real(r8), pointer, public :: acc_thermal_time_in_phase_patch(:,:) ! accumulated thermal time in each phase (deg-days) [patch, phase]
      real(r8), pointer, public :: acc_thermal_time_after_phase_patch(:,:) ! accumulated thermal time after each phase (deg-days) [patch, phase]
-     real(r8), pointer, public :: phase_target_thermal_time_path(:,:)  !target thermal time to finish a phase (deg-days) [path, phase]
+     real(r8), pointer, public :: phase_target_thermal_time_patch(:,:)  !target thermal time to finish a phase (deg-days) [path, phase]
      real(r8), pointer, public :: acc_vernalization_days_patch(:) ! accumulated vernalization days (for crops with vernalization) (unit: days)
 
      real(r8), pointer, public :: h2osoi_liq_24hr_col(:,:)  ! 24-hour average h2osoi_liq (kg/m2), just over 1:nlevsoi

src/agsys/science/AgSysCropTypeGeneric.F90

@@ -38,6 +38,8 @@ module AgSysCropTypeGeneric
      procedure :: get_stress_phenol_inductive_phase
      procedure :: get_target_tt_photosensitive_phase
      procedure :: get_target_tt_inductive_phase
+     procedure :: update_target_tt_for_phases
+     procedure :: get_final_leaf_no
   end type agsys_crop_type_generic
 
   !!a container type to store all the cultivars of a given crop type
@@ -141,5 +143,19 @@ function get_target_tt_inductive_phase(this, cumvd) result(target_tt)
     write(iulog, *) "This call to ", subname, " is illegal"
     call endrun(msg="Illegal call!")
   end function get_target_tt_inductive_phase
-
+
+  subroutine update_target_tt_for_phases (this, env, das, current_stage_index, phase_target_tt)
+    class(agsys_crop_type_generic),         intent(in) :: this
+    type (agsys_environmental_inputs_type), intent(in) :: env
+    integer,  intent(in)                    :: das
+    integer,  intent(in)                    :: current_stage_index
+    real(r8), intent(inout)                 :: phase_target_tt(:)
+  end subroutine update_target_tt_for_phases
+
+  function get_final_leaf_no(this, phase_target_tt) result(final_leaf_no)
+    class(agsys_crop_type_generic), intent(in) :: this
+    real(r8), intent(in)                       :: phase_target_tt(:)
+    real(r8) :: final_leaf_no
+  end function get_final_leaf_no
+
 end module AgSysCropTypeGeneric

src/agsys/science/AgSysCropTypeMaize.F90

@@ -5,26 +5,34 @@ module AgSysCropTypeMaize
                                                phase_type_generic, &
                                                phase_type_germinating, phase_type_emerging, phase_type_photo_sensitive
   use AgSysEnvironmentalInputs,         only : agsys_environmental_inputs_type
-  use AgSysUtils,                       only : response_curve_type, interpolation
+  use AgSysUtils,                       only : response_curve_type, interpolation, bound
   use AgSysCropTypePhotoSensitive,      only : agsys_crop_type_photosensitive
   implicit none
 
   type, extends(agsys_crop_type_photosensitive), public :: agsys_crop_type_maize
     private
     !public data members
-    real(r8),                  public :: tt_emerg_to_endjuv
-    real(r8),                  public :: tt_emerg_to_fi
-    real(r8),                  public :: tt_flag_to_flower
-    real(r8),                  public :: tt_flower_to_start_grain
-    real(r8),                  public :: tt_flower_to_maturity
-    real(r8),                  public :: tt_maturity_to_ripe
+    real(r8),                  public :: target_tt_emerg_to_endjuv
+    real(r8),                  public :: target_tt_flag_to_flower
+    real(r8),                  public :: target_tt_flower_to_start_grain
+    real(r8),                  public :: target_tt_flower_to_maturity
+    real(r8),                  public :: target_tt_maturity_to_ripe
 
+    real(r8),                  public :: leaf_init_rate
+    real(r8),                  public :: leaf_no_seed
+    real(r8),                  public :: leaf_no_at_emerg
+    real(r8),                  public :: leaf_no_min
+    real(r8),                  public :: leaf_no_max
+    real(r8),                  public :: leaf_no_critical
+    real(r8),                  public :: leaf_appearance_rate_early
+    real(r8),                  public :: leaf_appearance_rate_late
     real(r8),                  public :: potential_kernel_weight
     real(r8),                  public :: leaf_no_dead_const
     real(r8),                  public :: leaf_no_dead_slope
   contains
     procedure :: init
-
+    procedure :: update_target_tt_for_phases
+    procedure :: get_final_leaf_no
   end type agsys_crop_type_maize
 contains
   subroutine init(this)
@@ -38,15 +46,23 @@ subroutine init(this)
     !!!initialize the parameters
     !!!TODO(pb, 2019-11-14) now we intialize the variables to 0, 
     !!!but later can initialize them into nan
-    this%tt_emerg_to_endjuv         = 0._r8
-    this%tt_emerg_to_fi             = 0._r8
-    this%tt_flag_to_flower          = 0._r8
-    this%tt_flower_to_start_grain   = 0._r8
-    this%tt_flower_to_maturity      = 0._r8
-    this%tt_maturity_to_ripe        = 0._r8
-    this%potential_kernel_weight = 0._r8
-    this%leaf_no_dead_const      = 0._r8
-    this%leaf_no_dead_slope      = 0._r8
+    this%target_tt_emerg_to_endjuv         = 0._r8
+    this%target_tt_flag_to_flower          = 0._r8
+    this%target_tt_flower_to_start_grain   = 0._r8
+    this%target_tt_flower_to_maturity      = 0._r8
+    this%target_tt_maturity_to_ripe        = 0._r8
+
+    this%leaf_init_rate             = 0._r8
+    this%leaf_no_seed               = 0._r8
+    this%leaf_no_at_emerg           = 0._r8
+    this%leaf_no_min                = 0._r8
+    this%leaf_no_max                = 0._r8
+    this%leaf_no_critical           = 0._r8
+    this%leaf_appearance_rate_early = 0._r8
+    this%leaf_appearance_rate_late  = 0._r8
+    this%potential_kernel_weight    = 0._r8
+    this%leaf_no_dead_const         = 0._r8
+    this%leaf_no_dead_slope         = 0._r8
 
     !!!initialize the phases
     this%phases%num_phases=11
@@ -63,4 +79,46 @@ subroutine init(this)
                             phase_type_generic, phase_type_generic, phase_type_generic, phase_type_generic, phase_type_generic]
 
   end subroutine init
+
+  subroutine update_target_tt_for_phases (this, env, das, current_stage_index, phase_target_tt)
+    class(agsys_crop_type_maize),           intent(in) :: this
+    type (agsys_environmental_inputs_type), intent(in) :: env
+    integer,  intent(in)                    :: das
+    integer,  intent(in)                    :: current_stage_index
+    real(r8), intent(inout)                 :: phase_target_tt(:)
+
+    real(r8) :: leaf_no_final
+    real(r8) :: target_tt_emerg_to_flag 
+    if ((das == 0) .and. (current_stage_index ==1)) then !!! on sowing, I put two triggers here to make sure it is on sowing date
+      phase_target_tt(1) = 0._r8  !!germination is determined by soil water in seeding layer, not thermal time
+      phase_target_tt(2) = this%p_shoot_lag+this%p_sowing_depth*this%p_shoot_rate
+      phase_target_tt(3) = this%target_tt_emerg_to_endjuv
+      phase_target_tt(4) = this%rc_photoperiod_target_tt%y(1)
+      phase_target_tt(6) = this%target_tt_flag_to_flower
+      phase_target_tt(7) = this%target_tt_flower_to_start_grain
+      phase_target_tt(9) = this%target_tt_flower_to_maturity*0.01  !!!TODO(pb, 2019-11-15) 0.01 can be a parameter
+      phase_target_tt(8) = this%target_tt_flower_to_maturity-phase_target_tt(7)-phase_target_tt(9)
+      phase_target_tt(10)= this%target_tt_maturity_to_ripe
+      phase_target_tt(11)= 0._r8  !!ripe to harvest, this is depend on farmer's behavior
+    else if ((current_stage_index >=4) .and. (current_stage_index<5)) then  !!photosensitive phase
+      phase_target_tt(4) = this%get_target_tt_photosensitive_phase(env)
+
+    else if (current_stage_index>=5) then   !!TODO(pb, 2019-11-15) we only need to do this at one step!
+      leaf_no_final=this%get_final_leaf_no(phase_target_tt)
+      target_tt_emerg_to_flag=(leaf_no_final-this%leaf_no_critical-this%leaf_no_at_emerg)*this%leaf_appearance_rate_early + &
+                           this%leaf_no_critical*this%leaf_appearance_rate_late
+      phase_target_tt(5)=target_tt_emerg_to_flag-phase_target_tt(3)-phase_target_tt(4)
+
+    end if 
+  end subroutine update_target_tt_for_phases
+
+  function get_final_leaf_no(this, phase_target_tt) result(final_leaf_no)
+    class(agsys_crop_type_maize), intent(in) :: this
+    real(r8), intent(in)                  :: phase_target_tt(:)
+    real(r8) :: final_leaf_no
+
+    final_leaf_no =(phase_target_tt(3)+phase_target_tt(4))/this%leaf_init_rate+this%leaf_no_seed
+    final_leaf_no = bound (final_leaf_no, this%leaf_no_min, this%leaf_no_max)
+  end function get_final_leaf_no
+
 end module AgSysCropTypeMaize

src/agsys/science/AgSysPhenology.F90

@@ -141,6 +141,7 @@ subroutine AgSysRunPhenology(crop, env, soil_cond, &
 
     !!!!
 
+    call crop%update_target_tt_for_phases(env, days_after_sowing, current_stage_index, phase_target_tt)
     tt=get_daily_tt(env%tair_max, env%tair_min, crop%rc_tair_tt)
     select case (crop%phases%phase_type(current_stage_index))
       case(phase_type_generic)
@@ -253,12 +254,17 @@ subroutine AgSysRunPhenology(crop, env, soil_cond, &
       current_stage = new_stage
     end if
 
-    if ((int(current_stage) >= crop%phases%num_phases) .or. (current_stage < 0.0_r8)) then
+    current_stage_index=int(current_stage)
+    if ((current_stage_index > crop%phases%num_phases) .or. (current_stage < 0.0_r8)) then
       write(iulog, *) "Illegal stage number.."
       call endrun(msg="Illegal stage number...")
     end if
-
-    days_after_sowing=days_after_sowing+1
+    if (crop_alive) then
+      days_after_sowing=days_after_sowing+1
+    end if
+    if (current_stage_index == crop%phases%num_phases) then
+      crop_alive = .false.
+    end if 
   end subroutine AgSysRunPhenology
 
   subroutine initialize_on_sowing(shoot_lag, shoot_rate, sowing_depth, target_tt)

src/agsys/science/AgSysRoot.F90

@@ -44,13 +44,18 @@ end subroutine get_soil_condition
   function get_sw_avail_ratio (env, root) result(sw_avail_ratio)
     type(agsys_environmental_inputs_type), intent(in) :: env
     type(agsys_root_type),                 intent(in) :: root
-    real(r8) :: sw_avail_ratio 
+    real(r8) :: sw_avail_ratio
+
+    !TODO(pb, 2019-11-15) do real calculation later using root and soil water profiles
+    sw_avail_ratio=1._r8 
   end function get_sw_avail_ratio
 
   function get_sw_avail_for_layer (env, depth) result(sw_avail)
     type(agsys_environmental_inputs_type), intent(in) :: env
     real(r8), intent(in) :: depth
     real(r8) :: sw_avail
+    !TODO(pb, 2019-11-15) do real calculation later using root and soil depth profiles
+    sw_avail = env%h2osoi_liq_24hr(2)  !! current use the second layer which is 0.02-0.06 m
   end function get_sw_avail_for_layer
 
   function get_sw_avail (env, root) result(sw_avail)