! FILE crop.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine crpmd(j, je, k, n) hydrotop ! subroutine operat(j, je, k, n) crpmd ! subroutine growth(j, je, k, n) crpmd module crop use utilities, only : dp, path_max_length, log_info, log_debug, log_warn, log_error implicit none ! number of crops in crop.csv integer, save :: mcrdb = 0 ! current rotation year (land management, 1 - 3) integer, save :: iyrrot = 0 ! number of rotation years integer :: nrotyrs = 3 ! number of management operations during one year integer, save, dimension(:), allocatable :: mgt_nop ! day of operation read from landman.csv integer, save, dimension(:, :), allocatable :: mgt_idop ! operation code: 1 - planting, 2 - harvest & kill read from landman.csv integer, save, dimension(:, :), allocatable :: mgt_iopc ! crop number read from landman.csv integer, save, dimension(:, :), allocatable :: mgt_ncrp ! day of fertilization read from landman.csv integer, save, dimension(:, :), allocatable :: mgt_idfe ! N fertilizer kg/ha read from landman.csv real(dp), save, dimension(:, :), allocatable:: mgt_fen ! org N fertilizer kg/ha read from landman.csv real(dp), save, dimension(:, :), allocatable:: mgt_feno ! P fertilizer kg/ha read from landman.csv real(dp), save, dimension(:, :), allocatable:: mgt_fep ! management IDs in landman.csv integer, save, dimension(:), allocatable :: mgt_id ! management / rotation year in landman.csv integer, save, dimension(:), allocatable :: mgt_yr ! land use id read from landman.csv integer, save, dimension(:), allocatable :: mgt_lu_id ! Number of land management operations of current hydrotop integer, save :: cur_nop = 0 ! mop = max number of agric. operations within one year integer :: mop = 7 ! ida, to calc ndgro - number of growth days integer, save :: idayx = -99 ! number of growth days integer, save :: ndgro ! day to write crop yield for GIS output integer, save :: ndpri ! actual evapotranspiration, mm real(dp), save :: actual real(dp), save :: sdn = 0. ! sum N stress days real(dp), save :: sdp = 0. ! sum P stress days ! code for the day length effect in crop: 0 - without, 1 - with integer, save :: idlef ! switch code to print from crop() in old code, left to keep some aggregation alive integer, save :: icrop = 1 ! integer, save :: istyr ! crop number (database) integer, save, dimension(:, :), allocatable :: nucr ! vegetation cover, kg / ha real(dp), save, dimension(:, :), allocatable :: cva ! vegetation index, = 1 if vegetation is growing integer, save, dimension(:, :), allocatable :: igro ! harvest index heat unit real(dp), save, dimension(:, :), allocatable :: huharv ! harvest index real(dp), save, dimension(:, :), allocatable :: hia ! harvest index, adjusted real(dp), save, dimension(:, :), allocatable :: hiad ! fraction of root weight real(dp), save, dimension(:, :), allocatable :: rwt ! fresh organic N from residue in a layer, kg / ha real(dp), save, dimension(:, :, :), allocatable :: fon ! fresh org N, kg / ha real(dp), save, dimension(:, :, :), allocatable :: fop ! crop yield, kg / ha real(dp), save, dimension(:, :), allocatable :: yld ! crop yield for subbasin and soil, kg / ha real(dp), save, dimension(:, :), allocatable :: ylda ! actual transp. by plants, mm real(dp), save, dimension(:, :), allocatable :: swh ! potent. transp. by plants, mm real(dp), save, dimension(:, :), allocatable :: swp ! N uptake accumulated, kg / ha real(dp), save, dimension(:, :), allocatable :: snup ! P uptake in hydrotop, kg / ha real(dp), save, dimension(:, :), allocatable :: spup ! av yld per sub, soil, crop & aryld(j, k, icr): frac. area real(dp), save, dimension(:, :, :), allocatable :: avyld ! fraction of area by crop per sub, soil real(dp), save, dimension(:, :, :), allocatable :: aryld ! av. yld per soil, crop, kg / ha real(dp), save, dimension(:, :), allocatable :: avylds ! fraction of area by crop per soil real(dp), save, dimension(:, :), allocatable :: arylds ! day of operation integer, save, dimension(:), allocatable :: idop ! operation code: 1 - planting, 2 - harvest & kill integer, save, dimension(:), allocatable :: iopc ! crop number integer, save, dimension(:), allocatable :: ncrp ! day of fertilization integer, save, dimension(:), allocatable :: idfe ! amount of min N fertilizers applied, kg N / ha real(dp), save, dimension(:), allocatable :: fen ! amount of P fertilizers applied, kg P / ha real(dp), save, dimension(:), allocatable :: fep ! amount of org N fertilizers applied, kg N / ha real(dp), save, dimension(:), allocatable :: feno ! harvest index for crop (from database) real(dp), save, dimension(:), allocatable :: hi ! maximum root depth, mm real(dp), save, dimension(:), allocatable :: rdmx ! potential heat units required for maturity of crop integer, save, dimension(:), allocatable :: hun ! av yld per crop, kg / ha real(dp), save, dimension(:), allocatable :: avyldc ! fraction of area by crop real(dp), save, dimension(:), allocatable :: aryldc ! av yld per year, crop, kg / ha real(dp), save, dimension(:, :), allocatable :: avylda ! fraction of area by crop per year real(dp), save, dimension(:, :), allocatable :: arylda integer :: iww = 45 integer :: iwb = 36 integer :: iwr = 42 integer :: isba = 22 integer :: ipo = 20 integer :: icc = 51 ! mfe = max number of fertilizations integer :: mfe = 7 ! crop number integer, save, dimension(:), allocatable :: icnum ! crop name (abbreviation) character(len = 4), save, dimension(:), allocatable :: cnam ! complex number: fraction of grow. season, max corresp. LAI real(dp), save :: dlp1 = 0 ! complex number: fraction of grow. season, max corresp. LAI real(dp), save :: dlp2 = 0 ! fraction of nitrogen in yield, kg N / kg yield real(dp), save, dimension(:), allocatable :: cnyld ! fraction of phosphorus in yield, kg P / kg yield real(dp), save, dimension(:), allocatable :: cpyld ! 2nd point on radiation use efficiency curve real(dp), save, dimension(:), allocatable :: pt2 integer, save, dimension(4) :: k1 = (/ 9, 98, 915, 92/) integer, save, dimension(4) :: k6 = (/ 51, 78, 648, 0/) integer, save, dimension(4) :: k8 = (/ 20, 90, 215, 31/) integer, save, dimension(4) :: k9 = (/ 320, 73, 631, 49/) ! biomass - energy ratio for crop, kg m2 MJ - 1 ha - 1 d - 1 real(dp), save, dimension(:), allocatable :: be ! index for dormant period integer, save, dimension(:, :), allocatable :: nclc ! index for dormant period integer, save, dimension(:, :), allocatable :: avyldrot ! index for dormant period integer, save, dimension(:, :), allocatable :: aryldrot ! index for dormant period integer, save, dimension(:, :), allocatable :: irotup ! index for dormant period integer, save, dimension(:, :), allocatable :: iccup ! crop rotation type of the hydrotope integer, save, dimension(:), allocatable :: ihydRot ! fertility class of the hydrotope integer, save, dimension(:), allocatable :: ihydFert ! hydrotope's crop taken from ihydRot integer, save, dimension(:), allocatable :: ihydRotCrp character(len=path_max_length), save :: landmgtdat = 'landmgt.csv' ! file containing land management operations ! init. CN for cropland, cond 2 real(dp), save :: cnum2 = 1. ! total number of management operation years (lines in landman.csv) integer, save :: mgt_tot = 0 integer :: crop_input_file_id, crop_management_input_file_id ! Output variable ids integer :: & crop_yield_output_id = 0 character(len=path_max_length) :: crop_input_file = "crop.csv" character(len=path_max_length) :: crop_management_input_file = "crop_management.csv" namelist / CROP_PARAMETERS / & cnum2, & icc, & idlef, & ipo, & isba, & iwb, & iwr, & iww, & mfe, & mop, & nrotyrs contains subroutine crop_initialise(iyr, mb, meap, ms, nbyr) use input, only: input_open_file, input_count_rows, get_config_fid use output, only: output_register_hydrotope_var integer, intent(in) :: iyr, mb, meap, ms, nbyr istyr = iyr read(get_config_fid(), CROP_PARAMETERS) crop_yield_output_id = output_register_hydrotope_var("crop_yield", .false.) crop_input_file_id = input_open_file(crop_input_file) crop_management_input_file_id = input_open_file(crop_management_input_file) ! **** count number of managment operation years in landmgt.csv mgt_tot = input_count_rows(crop_management_input_file_id, .true.) ! **** count number of crops in crop.dat mcrdb = input_count_rows(crop_input_file_id, .true.) call crop_allocate(mb, mcrdb, meap, ms, nbyr) call crop_read_management_input end subroutine crop_initialise subroutine crop_read_input(bn1, bn2, bn3, bnu1, bnu2, bp1, bp2, bp3, bpu1, bpu2, cvm) !**** PURPOSE: THIS SUBROUTINE READS CROP PARAMETERS from crop.dat !**** CALLED IN: MAIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! ald1(ic) = shape parameter for the LAI developement equation ! ald2(ic) = shape parameter for the LAI developement equation ! almn(ic) = minimum Leaf Area Index (for forest and natural vegetation) ! be(ic) = biomass-energy ratio ! blai(ic) = maximum Leaf Area Index for crop ! bn1(ic) = nitrogen uptake parameter #1: normal fraction of N ! in crop biomass at emergence, kg N/kg biomass ! bn2(ic) = nitrogen uptake parameter #2: normal fraction of N ! in crop biomass at 0.5 maturity, kg N/kg biomass ! bn3(ic) = nitrogen uptake parameter #3: normal fraction of N ! in crop biomass at maturity, kg N/kg biomass ! bnu1(ic) = shape parameter to calc optimal N fraction in crop biomass ! bnu2(ic) = shape parameter to calc optimal N fraction in crop biomass ! bp1(ic) = phosphorus uptake parameter #1: normal fraction of P ! in crop biomass at emergence, kg P/kg biomass ! bp2(ic) = phosphorus uptake parameter #2: normal fraction of P ! in crop biomass at 0.5 maturity, kg P/kg biomass ! bp3(ic) = phosphorus uptake parameter #3: normal fraction of P ! P in crop biomass at maturity, kg P/kg biomass ! bpu1(ic) = shape parameter to calc optimal P fraction in crop biomass ! bpu2(ic) = shape parameter to calc optimal P fraction in crop biomass ! cnyld(ic) = fraction of nitrogen in yield, kg N/kg yield ! cpyld(ic) = fraction of phosphorus in yield, kg P/kg yield ! cvm(ic) = minimum value of C factor for water erosion ! dlai(ic) = fraction of growing season when leaf area declines ! dlp1 = complex number: fraction of grow. season, max corresp. LAI ! dlp2 = complex number: fraction of grow. season, max corresp. LAI ! hi(ic) = harvest index ! hun(ic) = potential heat units required for maturity of crop ! icnum(ic) = crop number ! ilcc(ic) = land cover category: ! (1) annual crop (row crop / small grain) ! (2) annual winter crop (row crop / small grain) ! (3) perennial (grass, brush, urban, water) ! (4) woods ! (5) annual legumes (row crop) ! (6) annual winter legumes (row crop) ! (7) perennial legumes (grass) ! pt2(ic) = 2nd point on radiation use efficiency curve: ! complex number: The value to the left of (.) is a CO2 ! atm. concentration higher than the ambient (in units of ! microliters CO2/liter air). The value to the right of (.) ! is the corresponding biomass-energy ratio divided by 100 ! rdmx(ic) = maximum root depth (m), then converted to mm ! sla(ic) = specific leaf area, m2/kg ! tb(ic) = base temperature for plan growth, degrees C ! to(ic) = optimal temperature for plant growth, degrees C ! >>>>> ! >>>>> STATIC PARAMETERS ! b1 = intermediate parameter ! b2 = intermediate parameter ! b3 = intermediate parameter ! ic = count parameter ! titl = text ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use vegetation, only : & vegetation_s_curve_parameters, & ald1, & ald2, & blai, & almn, & ilcc, & sla, & to, & tb, & dlai use input, only : & read_integer_column, & read_string_column, & read_real_column, & input_count_rows real(dp), dimension(:), intent(inout) :: bn1 real(dp), dimension(:), intent(inout) :: bn2 real(dp), dimension(:), intent(inout) :: bn3 real(dp), dimension(:), intent(inout) :: bnu1 real(dp), dimension(:), intent(inout) :: bnu2 real(dp), dimension(:), intent(inout) :: bp1 real(dp), dimension(:), intent(inout) :: bp2 real(dp), dimension(:), intent(inout) :: bp3 real(dp), dimension(:), intent(inout) :: bpu1 real(dp), dimension(:), intent(inout) :: bpu2 real(dp), dimension(:), intent(inout) :: cvm real(dp), dimension(:), allocatable :: dlp1_array real(dp), dimension(:), allocatable :: dlp2_array integer ic real(dp) b1, b2, b3 real(dp) x3, x4 x3 = .5 x4 = 1. allocate(dlp1_array(mcrdb)) allocate(dlp2_array(mcrdb)) call read_integer_column(crop_input_file_id, "icnum", icnum, 0) call read_string_column(crop_input_file_id, "cpnm", cnam, "") call read_real_column(crop_input_file_id, "be", be, 0.0_dp) call read_real_column(crop_input_file_id, "hi", hi, 0.0_dp) call read_real_column(crop_input_file_id, "blai", blai, 0.0_dp) call read_real_column(crop_input_file_id, "dlp1", dlp1_array, 0.0_dp) call read_real_column(crop_input_file_id, "dlp2", dlp2_array, 0.0_dp) call read_real_column(crop_input_file_id, "bn1", bn1, 0.0_dp) call read_real_column(crop_input_file_id, "bn2", bn2, 0.0_dp) call read_real_column(crop_input_file_id, "bn3", bn3, 0.0_dp) call read_real_column(crop_input_file_id, "bp1", bp1, 0.0_dp) call read_real_column(crop_input_file_id, "bp2", bp2, 0.0_dp) call read_real_column(crop_input_file_id, "bp3", bp3, 0.0_dp) call read_real_column(crop_input_file_id, "cnyld", cnyld, 0.0_dp) call read_real_column(crop_input_file_id, "cpyld", cpyld, 0.0_dp) call read_real_column(crop_input_file_id, "rdmx", rdmx, 0.0_dp) call read_real_column(crop_input_file_id, "cvm", cvm, 0.0_dp) call read_real_column(crop_input_file_id, "pt2", pt2, 0.0_dp) call read_integer_column(crop_input_file_id, "hun", hun, 0) ! Vegetation variables call read_real_column(crop_input_file_id, "almn", almn, 0.0_dp) call read_integer_column(crop_input_file_id, "ird", ilcc, 0) call read_real_column(crop_input_file_id, "sla", sla, 0.0_dp) call read_real_column(crop_input_file_id, "to", to, 0.0_dp) call read_real_column(crop_input_file_id, "tb", tb, 0.0_dp) call read_real_column(crop_input_file_id, "dlai", dlai, 0.0_dp) !**** READ 5 - crop parameter database (unit 5-crop.dat) do ic = 1, mcrdb dlp1 = dlp1_array(ic) dlp2 = dlp2_array(ic) cvm(ic) =log(cvm(ic)) rdmx(ic) = rdmx(ic) * 1000. !**** RECALCULATE: dlp1, dlp2 need to be broken into two variables !**** Isolate number to left of decimal: b1, b2 - fraction of the grow. season b1 = .01 * Int(dlp1) b2 = .01 * Int(dlp2) !**** Isolate number to right of decimal, dlp1, 2 is now max LAI corr. to b1, 2 dlp1 = dlp1 - Int(dlp1) dlp2 = dlp2 - Int(dlp2) !#### determine shape parameters for the LAI developement equation call vegetation_s_curve_parameters(dlp1, dlp2, b1, b2, ald1(ic), ald2(ic)) !**** Calc nitrogen uptake parameters b1 = bn1(ic) - bn3(ic) b2 = 1. - (bn2(ic) - bn3(ic)) / b1 b3 = 1. - .00001 / b1 !#### CALL ASCRV call vegetation_s_curve_parameters(b2, b3, x3, x4, bnu1(ic), bnu2(ic)) !**** CORRECT bp3 & Calc phosphorus uptake parameters if (bp2(ic) - bp3(ic) < .0001) bp3(ic) = .75 * bp3(ic) b1 = bp1(ic) - bp3(ic) b2 = 1. - (bp2(ic) - bp3(ic)) / b1 b3 = 1. - .00001 / b1 !#### CALL ASCRVg if (b2 > 0. .AND. b2 < 1.) then call vegetation_s_curve_parameters(b2, b3, x3, x4, bpu1(ic), bpu2(ic)) else call log_error("crop_read_input", & "Value of b2 in call vegetation_s_curve_parameters() is causing an error "// & "(ic, b2, b1, bp2, bp3):", int=ic, reals=(/b2, b1, bp2, bp3/)) end if end do deallocate(dlp1_array) deallocate(dlp2_array) call log_debug("crop_read_input", 'Crop parameters READ!') end subroutine crop_read_input subroutine crop_read_management_input !------------------------------------------------------------------------------- ! Author : stefan.liersch@pik-potsdam.de ! Date : 2009-11-25 ! modified: 2009-12-08 ! ! PURPOSE : THIS SUBROUTINE READS CROP MANAGEMENT DATA FROM FILE: landman.csv ! ! CALLED : from program main !------------------------------------------------------------------------------- use input, only : read_integer_column, read_real_column character(len=80) :: colno integer :: j = 0 call read_integer_column(crop_management_input_file_id, "manag_id", mgt_id, 0) call read_integer_column(crop_management_input_file_id, "nop", mgt_nop, 0) call read_integer_column(crop_management_input_file_id, "year", mgt_yr, 0) call read_integer_column(crop_management_input_file_id, "lu_id", mgt_lu_id, 0) do j = 1, mgt_nop(1) write(colno, '(I5)') j colno = trim(adjustl(colno)) ! call read_integer_column(1818, "idop"//colno, temp(1, j), 0) call read_integer_column(crop_management_input_file_id, "idop"//colno, mgt_idop(:, j), 0) call read_integer_column(crop_management_input_file_id, "iopc"//colno, mgt_iopc(:, j), 0) call read_integer_column(crop_management_input_file_id, "ncrp"//colno, mgt_ncrp(:, j), 0) call read_integer_column(crop_management_input_file_id, "idfe"//colno, mgt_idfe(:, j), 0) call read_real_column(crop_management_input_file_id, "fen"//colno, mgt_fen(:, j), 0.0_dp) call read_real_column(crop_management_input_file_id, "feno"//colno, mgt_feno(:, j), 0.0_dp) call read_real_column(crop_management_input_file_id, "fep"//colno, mgt_fep(:, j), 0.0_dp) end do close(crop_management_input_file_id) end subroutine crop_read_management_input subroutine crop_allocate(mb, mcrdb, meap, ms, nbyr) integer, intent(in) :: mb, mcrdb, meap, ms, nbyr allocate(aryld(mb, ms, mcrdb)) allocate(arylda(nbyr, mcrdb)) allocate(aryldc(mcrdb)) allocate(aryldrot(23, mcrdb)) allocate(arylds(ms, mcrdb)) allocate(avyld(mb, ms, mcrdb)) allocate(avylda(nbyr, mcrdb)) allocate(avyldc(mcrdb)) allocate(avyldrot(23, mcrdb)) allocate(avylds(ms, mcrdb)) allocate(be(mcrdb)) allocate(cnam(mcrdb)) ! crop name (4 letters) allocate(cnyld(mcrdb)) allocate(cpyld(mcrdb)) allocate(cva(mb, meap)) allocate(fen(mfe)) allocate(feno(mfe)) allocate(fep(mfe)) allocate(fon(mb, meap, 2)) allocate(fop(mb, meap, 2)) allocate(hi(mcrdb)) allocate(hia(mb, meap)) allocate(hiad(mb, meap)) allocate(huharv(mb, meap)) allocate(hun(mcrdb)) allocate(iccup(mb, meap)) allocate(icnum(mcrdb)) allocate(idfe(mfe)) allocate(idop(mop)) allocate(igro(mb, meap)) allocate(ihydFert(mb * meap)) allocate(ihydRot(mb * meap)) allocate(ihydRotCrp(mb * meap)) allocate(iopc(mop)) allocate(irotup(mb, meap)) allocate(mgt_fen(mgt_tot, mop)) allocate(mgt_feno(mgt_tot, mop)) allocate(mgt_fep(mgt_tot, mop)) allocate(mgt_idfe(mgt_tot, mop)) allocate(mgt_idop(mgt_tot, mop)) allocate(mgt_iopc(mgt_tot, mop)) allocate(mgt_ncrp(mgt_tot, mop)) allocate(mgt_nop(mgt_tot)) allocate(nclc(mb, meap)) allocate(ncrp(mop)) allocate(nucr(mb, meap)) allocate(pt2(mcrdb)) allocate(rdmx(mcrdb)) allocate(rwt(mb, meap)) allocate(snup(mb, meap)) allocate(spup(mb, meap)) allocate(swh(mb, meap)) allocate(swp(mb, meap)) allocate(yld(mb, meap)) allocate(ylda(mb, ms)) allocate(mgt_id(mgt_tot)) allocate(mgt_lu_id(mgt_tot)) allocate(mgt_yr(mgt_tot)) aryld = 0. arylda = 0. arylda = 0. aryldc = 0. aryldrot = 0 arylds = 0. avyld = 0. avylda = 0. avylda = 0. avyldc = 0. avyldrot = 0 avylds = 0. be = 0. cnyld = 0. cnyld = 0. cpyld = 0. cpyld = 0. cva = 0. fen = 0. feno = 0. fep = 0. fon = 0. fop = 0. hi = 0. hia = 0. hiad = 0. huharv = 0. hun = 0 iccup = 0 icnum = 0 idfe = 0 idop = 0 igro = 0 ihydRot = 0 ihydRot = 0 ihydRot = 0 iopc = 0 irotup = 0 mgt_fen = 0. mgt_fen = 0. mgt_feno = 0. mgt_feno = 0. mgt_fep = 0. mgt_fep = 0. mgt_idfe = 0 mgt_idfe = 0 mgt_idop = 0 mgt_idop = 0 mgt_iopc = 0 mgt_iopc = 0 mgt_ncrp = 0 mgt_ncrp = 0 mgt_nop = 0 nclc = 0 ncrp = 0 nucr = 0 pt2 = 0. rdmx = 0. rwt = 0. snup = 0. spup = 0. swh = 0. swp = 0. yld = 0. ylda = 0. mgt_id = 0 mgt_lu_id = 0 mgt_yr = 0 end subroutine crop_allocate subroutine crop_deallocate deallocate(aryld) deallocate(avyld) deallocate(arylda) deallocate(aryldc) deallocate(aryldrot) deallocate(arylds) deallocate(avylda) deallocate(avyldc) deallocate(avyldrot) deallocate(avylds) deallocate(be) deallocate(cnam) deallocate(cnyld) deallocate(cpyld) deallocate(cva) deallocate(fen) deallocate(feno) deallocate(fep) deallocate(fon) deallocate(fop) deallocate(hi) deallocate(hia) deallocate(hiad) deallocate(huharv) deallocate(hun) deallocate(iccup) deallocate(icnum) deallocate(idfe) deallocate(idop) deallocate(igro) deallocate(iopc) deallocate(irotup) deallocate(nclc) deallocate(ncrp) deallocate(nucr) deallocate(pt2) deallocate(rdmx) deallocate(rwt) deallocate(snup) deallocate(spup) deallocate(swh) deallocate(swp) deallocate(yld) deallocate(ylda) end subroutine crop_deallocate !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine crop_initialise_hydrotope(nsub, nhru, iy, mstruc) !------------------------------------------------------------------------------- ! Author : stefan.liersch@pik-potsdam.de ! Date : 2009-11-25 ! modified: 2009-12-08 ! ! PURPOSE : THIS SUBROUTINE INITIALIZES CROP MANAGEMENT DATA AT HRU LEVEL ! ! CALLED : from subroutine hydrotop !------------------------------------------------------------------------------- !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! fen(if) = amount of min N fertilizers applied, kg N/ha ! feno(if) = amount of org N fertilizers applied, kg N/ha ! fep(if) = amount of P fertilizers applied, kg P/ha ! icc = index for cover crop corr. number in crop database ! idfe(if) = day of fertilization ! idop(io) = day of operation ! iopc(io) = operation code: 1 - planting, 2 - harvest & kill ! 3 - harvest, 4 - kill integer, intent(in) :: iy integer, dimension(:, :, :), intent(in) :: mstruc integer, intent(in) :: nsub, nhru integer :: i = 0, j = 0 ! number of years of rotation (3) + one initialization year (year 0) integer :: rotation integer :: management rotation = nrotyrs + 1 ! find the correct index for cropland management ! this depends on the mgt_id and current rotation year management = mstruc(nsub, nhru, 4) i = management * rotation - rotation + 1 if (iy > 1) i = i + iyrrot ! iyrrot is calculated in mainpro cur_nop = mgt_nop(i) ! hydrotop - specific number of mgt operations do j = 1, mgt_nop(i) idop(j) = mgt_idop(i, j) iopc(j) = mgt_iopc(i, j) ncrp(j) = mgt_ncrp(i, j) idfe(j) = mgt_idfe(i, j) fen(j) = mgt_fen(i, j) feno(j) = mgt_feno(i, j) fep(j) = mgt_fep(i, j) end do end subroutine crop_initialise_hydrotope subroutine crop_process(j, je, k, n, wet, disLastDay, additionalGwUptake, avt, bWAM_Module, dart, daycounter, es, fc, flu, frar, humi, ida, iy, iyr, mstruc, nbyr, nn, nveg, pit, ra, sbar, sep, ste, tmn, tx, uap, ylc, yls, z, bSnowModule, tmit) !**** PURPOSE: THIS SUBROUTINE CALCULATES DAILY POTENTIAL & ACTUAL GROWTH ! OF TOTAL PLANT BIOMASS AND ROOTS AND CALCULATES LEAF AREA INDEX. ! IT ADJUSTS DAILY BIOMASS TO WATER, TEMP. & NUTR. STRESS. !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! actual = actual evapotranspiration, mm ! alai(j, je) = leaf area index ! cva(j, je) = vegetation cover, kg/ha ! dm(j, je) = total biomass, kg/ha ! ep = plant transpiration, mm ! es = soil evaporation, mm ! g(j, je) = fraction of heat units to maturity accumulated ! icrop = switch code to print from crop() ! icrsb = number of subbasin to print from crop(), if icrop = 1 ! icrso = number of soil to print from crop(), if icrop = 1 ! ida = current day ! igro(j, je) = vegetation index, =1 if vegetation is growing ! rd(j, je) = root depth, mm ! rsd(j, je, 2)= crop residue in two upper soil layers, kg/ha ! ts = temperature stress factor ! uap = P uptake in hydrotope, kg/ha ! ws(j, je) = water stress factor ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use vegetation, only : vegetation_water_stress, olai, dm, ws, rsd, rd, g, alai, ts, ep real(dp), dimension(:), intent(inout) :: additionalGwUptake real(dp), dimension(:), intent(inout) :: disLastDay real(dp), dimension(:), intent(in) :: avt logical, intent(in) :: bWAM_Module real(dp), dimension(:), intent(in) :: dart integer, intent(in) :: daycounter real(dp), intent(in) :: es real(dp), dimension(:, :), intent(in) :: fc real(dp), dimension(:), intent(in) :: flu real(dp), dimension(:, :), intent(in) :: frar real(dp), dimension(:), intent(in) :: humi integer, intent(in) :: ida integer, intent(in) :: iy integer, intent(in) :: iyr integer, dimension(:, :, :), intent(in) :: mstruc integer, intent(in) :: nbyr integer, intent(in) :: nn integer, dimension(:, :), intent(in) :: nveg real(dp), intent(in) :: pit real(dp), dimension(:), intent(in) :: ra real(dp), dimension(:), intent(in) :: sbar real(dp), intent(inout) :: sep real(dp), dimension(:, :, :), intent(inout) :: ste real(dp), dimension(:), intent(in) :: tmn real(dp), dimension(:), intent(in) :: tx real(dp), intent(out) :: uap real(dp), dimension(:), intent(in) :: ylc real(dp), dimension(:), intent(in) :: yls real(dp), dimension(:, :), intent(in) :: z logical, intent(in) :: bSnowModule real(dp), intent(in) :: tmit integer j, je, n, k, wet uap = 0. ts = 0. !**** CALC vegetation cover cva(j, je) = .8 * dm(j, je) + rsd(j, je, 1) !#### CALL OPERAT call crop_operation(j, je, k, alai, dm, frar, g, ida, iy, olai, rd, rsd, ws) !#### CALL WSTRESS to COMPUTE WATER STRESS if (igro(j, je) .ge. 1) call vegetation_water_stress(j, je, k, n, wet, disLastDay, additionalGwUptake, bWAM_Module, dart, daycounter, fc, frar, humi, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, rdmx, sbar, sep, ste, z) actual = ep + es !#### CALL GROWTH call crop_growth(avt, bSnowModule, flu, frar, ida, j, je, n, nn, nveg, pit, ra, tmit, tmn, tx, ylc, yls) return end subroutine crop_process subroutine crop_operation(j, je, k, alai, dm, frar, g, ida, iy, olai, rd, rsd, ws) !**** PURPOSE: TO DEFINE PLANT OPERATIONS !**** CALLED IN: CRPMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alai(j, je) = leaf area index ! aryld(j, k, icr) = fraction of area by crop per sub, soil ! arylda(iy, icr) = fraction of area by crop per year ! aryldc(icr) = fraction of area by crop ! arylds(k, icr) = fraction of area by crop per soil ! avyld(j, k, icr) = av. yld per sub, soil, crop, kg/ha ! avylda(iy, icr) = av yld per year, crop, kg/ha ! avyldc(icr) = av yld per crop, kg/ha ! avylds(k, icr) = av. yld per soil, crop, kg/ha ! cva(j, je) = vegetation cover, kg/ha ! dm(j, je) = total biomass, kg/ha ! disLastDay = discharge of last day, m3/s ! fon(j, je, l) = fresh org N, kg/ha ! fop(j, je, l) = fresh org P, kg/ha ! frar(j, je) = fractional area of hydrotope in subbasin ! g(j, je) = fraction of heat units to maturity accumulated ! hi(icr) = harvest index for crop (database), for maize & potat. ! hia(j, je) = harvest index ! hiad(j, je) = harvest index, adjusted ! huharv(j, je) = harvest index heat unit ! icc = index for cover crop corr. number in crop database ! ida = current day ! idayx = par = ida, to calc ndgro - number of growth days ! idop(5, iop) = day of operation ! igro(j, je) = vegetation index, =1 if yes ! iopc(5, iop) = opeartion code: 1 - planting, ... ! ipo = index for potatoes corr. number in crop database ! isba = index for s. barley corr. number in crop databas ! istyr = starting year ! iwb = index for w. barley corr. number in crop databas ! iwr = index for w. rye corr. number in crop databas ! iww = index for w. wheat corr. number in crop databas ! iy = current year as counter (1, ..., nbyr) ! ncrp(iop) = crop number ! ndgro = number of growth days ! ndpri = day to write crop yield for GIS output ! nucr(j, je) = crop number (database) ! olai(j, je) = alai(j, je) - leaf area index ! rsd(j, je, 2) = residue, kg/ha ! rwt(j, je) = fraction of root weight ! sbar(j) = subbasin area, m2 ! snup(j, je) = N uptake, kg/ha ! spup(j, je) = P uptake, kg/ha ! swh(j, je) = actual transp. by plants, mm ! swp(j, je) = potent. transp. by plants, mm ! ws(j, je) = water stress ! yld(j, je) = crop yield, kg/ha ! ylda(j, k) = crop yield for subbasin and soil, kg/ha ! >>>>> ! >>>>> STATIC PARAMETERS ! icr = local par ! ii = local par ! ioper = local par ! xx = local par ! yield = current yield ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ use output, only: output_store_hydrotope_value !**** Include common parameters real(dp), dimension(:, :), intent(inout) :: alai real(dp), dimension(:, :), intent(inout) :: dm real(dp), dimension(:, :), intent(in) :: frar real(dp), dimension(:, :), intent(inout) :: g integer, intent(in) :: ida integer, intent(in) :: iy real(dp), dimension(:, :), intent(out) :: olai real(dp), dimension(:, :), intent(inout) :: rd real(dp), dimension(:, :, :), intent(inout) :: rsd real(dp), dimension(:, :), intent(inout) :: ws integer j, je, k integer icr, ii, ioper real(dp) xx, yield !*********************************************************** START IF (IGRO=0) !**** CHECK if day of planting, then goto 10 - planting if (igro(j, je) .eq. 0) then do ii = 1, mop if (ida .eq. idop(ii)) then nucr(j, je) = ncrp(ii) ioper = iopc(ii) if (ioper .eq. 1) then igro(j, je) = 1 g(j, je) = 0. dm(j, je) = 0.01 snup(j, je) = 0. spup(j, je) = 0. swh(j, je) = 0. swp(j, je) = 0. huharv(j, je) = 0. hia(j, je) = 0. olai(j, je) = 0. rwt(j, je) = 0. idayx = - 99 ndgro = 0 EXIT end if endif end do endif !*********************************************************** END IF (IGRO=0) !*********************************************************** START IF(IGRO=1) !**** CHECK if day of harvest and kill if (igro(j, je) .eq. 1) then do ii = 1, mop if (ida .eq. idop(ii)) then nucr(j, je) = ncrp(ii) ioper = iopc(ii) if (ioper .eq. 2) then !**** CALC HARVEST AND KILL igro(j, je) = 0 if (hiad(j, je) .gt. hi(nucr(j, je))) hiad(j, je) = hi(nucr(j, je)) if (hiad(j, je) .gt. 1.) hiad(j, je) = 1. !**** CALC residue & fresh org N and P (no residue, fon, fop for cover crop) if (nucr(j, je) .ne. icc) then rsd(j, je, 1) = 0.25 * (1. - rwt(j, je)) * (1. - hiad(j, je)) * dm(j, je) & + rsd(j, je, 1) rsd(j, je, 2) = 0.75 * (1. - rwt(j, je)) * (1. - hiad(j, je)) * dm(j, je) & + rsd(j, je, 2) if (rsd(j, je, 1) .le. 0.) rsd(j, je, 1) = 1.e-6 if (rsd(j, je, 2) .le. 0.) rsd(j, je, 2) = 1.e-6 fon(j, je, 1) = rsd(j, je, 1) * .008 fon(j, je, 2) = rsd(j, je, 2) * .008 fop(j, je, 1) = rsd(j, je, 1) * .0011 fop(j, je, 2) = rsd(j, je, 2) * .0011 endif !**** CALC yield ! ATTN! Harvest index used for grains: hia(), for maize, potatoes: hi() ! No yield for cover crop (icc) if (nucr(j, je) .eq. iww .or. nucr(j, je) .eq. iwb .or. & nucr(j, je) .eq. iwr .or. nucr(j, je) .eq. isba) then yield = 0.85 * dm(j, je) * hia(j, je) else if (nucr(j, je) .eq. ipo) then yield = 1.00 * dm(j, je) * hi(nucr(j, je)) else yield = 0.85 * dm(j, je) * hi(nucr(j, je)) endif if (nucr(j, je) .eq. icc) yield = 0. ylda(j, k) = yield yld(j, je) = yld(j, je) + yield call output_store_hydrotope_value(crop_yield_output_id, j, je, yield) dm(j, je) = 0. rd(j, je) = 0. ws(j, je) = 1. alai(j, je) = 0. hia(j, je) = 0. cva(j, je) = rsd(j, je, 1) + rsd(j, je, 2) g(j, je) = 0. idayx = 0 icr = nucr(j, je) !**** CALC Day to write crop yield for GIS output (except cover crop) if (icr .ne. icc) ndpri = ida + 3 !**** CALC average yield ! avyld(j, k, icr) av yld per sub, soil, crop & aryld(j, k, icr): frac. area ! avylds(k, icr) av yld per soil, crop & arylds(k, icr): frac. area ! avyldc(icr) av yld per crop & aryldc(icr): frac. area ! avylda(iy, icr) av yld per year, crop & arylda(iy, icr): frac. area if (icrop == 1) then avyld(j, k, icr) = avyld(j, k, icr) + ylda(j, k) * frar(j, je) / 100. aryld(j, k, icr) = aryld(j, k, icr) + frar(j, je) endif avylds(k, icr) = avylds(k, icr) + ylda(j, k) * frar(j, je) / 100. arylds(k, icr) = arylds(k, icr) + frar(j, je) avyldc(icr) = avyldc(icr) + ylda(j, k) * frar(j, je) / 100. aryldc(icr) = aryldc(icr) + frar(j, je) avylda(iy, icr) = avylda(iy, icr) + ylda(j, k) * frar(j, je) / 100. arylda(iy, icr) = arylda(iy, icr) + frar(j, je) return !**** END CALC HARVEST AND KILL end if if (ioper .eq. 3) then !**** CALC HARVEST ONLY - CUTTING, NO KILL if (hiad(j, je) .gt. hi(nucr(j, je))) & hiad(j, je) = hi(nucr(j, je)) yield = (1. - rwt(j, je)) * dm(j, je) * hiad(j, je) yld(j, je) = yld(j, je) + yield xx = dm(j, je) dm(j, je) = dm(j, je) - yield alai(j, je) = alai(j, je) * dm(j, je) / xx g(j, je) = g(j, je) * dm(j, je) / xx return !**** END CALC HARVEST ONLY end if if (ioper .eq. 4) then !**** CALC KILL igro(j, je) = 0 dm(j, je) = 0. ws(j, je) = 1. alai(j, je) = 0. cva(j, je) = 0. g(j, je) = 0. rd = 0. !**** END CALC KILL end if endif end do endif !*********************************************************** END IF (IGRO=1) return end subroutine crop_operation subroutine crop_growth(avt, bSnowModule, flu, frar, ida, j, je, n, nn, nveg, pit, ra, tmit, tmn, tx, ylc, yls) !**** PURPOSE: TO SIMULATE PLANT GROWTH !**** CALLED IN: CRPMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! actual = actual evapotranspiration, mm ! alai(j, je) = leaf area index ! ald1(icr) = shape parameter for the LAI developement equation ! ald2(icr) = shape parameter for the LAI developement equation ! be(icr) = biomass-energy ratio for crop, kg m2 MJ-1 ha-1 d-1 ! blai(icr) = max LAI for crop ! dlai(icr) = fraction of season, when LAI declines ! dm(j, je) = total biomass, kg/ha ! flu(j) = fraction of subbasin area in the basin ! frar(j, je) = fractional areas of hydrotope in subbasin ! g(j, je) = fraction of heat units to maturity accumulated, 1/1 ! hi(icr) = harvest index for crop (from database) ! hia(j, je) = harvest index ! hiad(j, je) = harvest index, adjusted ! huharv(j, je) = harvest index heat unit ! hun(icr) = potential heat units required for maturity of crop ! ialpha = switch parameter for CO2 effect on net photosynth. ! ida = current day ! idayx = ida, to calc ndgro - number of growth days ! idlef = code for the day length effect in crop 0/1 ! igro(j, je) = vegetation index, =1 if yes ! ilcc(icr) = land cover category --> readcrp ! ndgro = number of growth days ! nucr(j, je) = crop number (database) ! olai(j, je) = alai(j, je) - leaf area index !block pit = 58.13 ! potentl = potential evapotranspiration, mm ! ra(j) = solar radiation in subbasin j, J/cm^2 ! rd(j, je) = root depth, mm ! rdmx(icr) = maximum root depth, mm ! rwt(j, je) = fraction of root weight ! sdn = sum N stress days ! sdp = sum P stress days ! sdt = sum temp stress days ! sdw = sum water stress days ! strsn = N stress factor ! strsp = P stress factor ! swh(j, je) = actual transp. by plants, mm ! swp(j, je) = potent. transp. by plants, mm ! tb(icr) = base temp. for plant growth, degree C ! ts = temp. stress ! tsav(j, je) = temp. stress, accumulated ! tx(j) = average daily temp., degree C ! ws(j, je) = water stress ! wsav(j, je) = water stress, accumulated ! ylc(j) = cos(lat()/clt), lat() - lat, clt=57.296 ! yls(j) = sin(lat()/clt), lat() - lat, clt=57.296 ! >>>>> ! >>>>> STATIC PARAMETERS ! beadj = adjusted on CO2 biomass-energy ratio ! ch = local par ! dayl = day length ! ddm = delta dm ! delg = delta g ! deltalai = delta LAI ! dlef = local par ! duma = local par ! f = fraction of plant's maximum leaf area index corresponding to ! to a given fraction of potential heat units for plant ! ff = delta f for the day ! h = local par ! heat = local par ! par = local par ! reg = local par ! sd = local par ! tgx = local par ! x1 = local par ! x2 = local par ! x3 = local par ! x4 = local par ! xi = local par ! xinc = local par ! xy = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use vegetation, only : & vegetation_adjust_energy_ratio, & vegetation_s_curve, & vegetation_temperature_stress, & vegetation_water_stress, & ald1, & wsav, & sdt, & ialpha, & blai, & olai, & tsav, & ald2, & dm, & tb, & g, & alai, & ilcc, & ws, & to, & dlai, & potentl, & sdw, & rd, & ts use nutrient, only : & nutrient_phosphorus_uptake, & nutrient_nitrate_uptake, & strsp, & strsn real(dp), dimension(:), intent(in) :: avt real(dp), dimension(:), intent(in) :: flu real(dp), dimension(:, :), intent(in) :: frar integer, intent(in) :: ida integer, intent(in) :: nn integer, dimension(:, :), intent(in) :: nveg real(dp), intent(in) :: pit real(dp), dimension(:), intent(in) :: ra real(dp), dimension(:), intent(in) :: tmn real(dp), dimension(:), intent(in) :: tx real(dp), dimension(:), intent(in) :: ylc real(dp), dimension(:), intent(in) :: yls logical, intent(in) :: bSnowModule real(dp), intent(in) :: tmit integer j, je, n real(dp) beadj, ch, dayl, ddm, delg, deltalai, dlef, duma, f, ff real(dp) h, heat, par, reg, sd, tgx, x1, x2, x3, x4, xi, xinc, xy x1 = 0.05 x2 = 0.9 x3 = 10.05 x4 = 50.90 !ToDo: dlef = 0. ! ??? ! dlef is not initialized if g(j, je) >=0.5 .AND. idlef.NOT.0 dlef = 0. !*********************************************************** START IF (IGRO>=1) !**** CALC DAILY INCREASE IN HEAT UNITS delg & ACCUMULATED HEAT UNITS g() if (igro(j, je) .ge. 1) then if (idayx .ne. ida .and. g(j, je) .gt. 0.) then ndgro = ndgro + 1 idayx = ida end if !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then delg = (tmit - tb(nucr(j, je))) / hun(nucr(j, je)) else delg = (tx(j) - tb(nucr(j, je))) / hun(nucr(j, je)) end if !########################### !**** CALC day length (from clgen) xi = ida sd = .4102 * sin((xi - 80.25) / pit) ch = - yls(j) * tan(sd) / ylc(j) if (ch .lt. 1.) then if (ch .le. - 1.) then h = 3.1416 else h = acos(ch) end if else h = 0. end if dayl = 7.72 * h !**** Correction of delg on day length (Jan Gr�fe) if (g(j, je) .lt. 0.5) then dlef = 0.25 + 0.75 * (dayl - 8.) / 8. if (dlef .gt. 1.) dlef = 1. if (dlef .lt. 0.25) dlef = 0.25 endif if (idlef .eq. 0) then delg = delg * 1. else delg = delg * dlef endif if (delg .lt. 0.) delg = 0. g(j, je) = g(j, je) + delg if (g(j, je) .gt. 1.) g(j, je) = 1. !*********************************************************** START IF (G<=1) !**** GROWTH SEASON if (g(j, je) .le. 1.) then !#### CALL TSTRESS - CALC TEMPERATURE STRESS !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then tgx = tmit - tb(nucr(j, je)) else tgx = tx(j) - tb(nucr(j, je)) end if !########################### if (tgx .le. 0.) then ts = 0. else call vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx) end if !**** COMPUTE ROOT DEPTH if (ilcc(nucr(j, je)) .le. 2 .or. ilcc(nucr(j, je)) .eq. 4 & .or. ilcc(nucr(j, je)) .eq. 5 & .OR.nucr(j, je) == icc ) then ! cover crop included by S. Liersch rd(j, je) = 2.5 * g(j, je) * rdmx(nucr(j, je)) if (rd(j, je) .gt. rdmx(nucr(j, je))) rd(j, je) = rdmx(nucr(j, je)) else rd(j, je) = rdmx(nucr(j, je)) endif !**** COMPUTE potential increase in biomass - ddm ! STANDARD VERSION: WITHOUT CO2 FERTILIZATION par = .005 * ra(j) * (1. - exp(- .65 * (alai(j, je) + .05))) ddm = be(nucr(j, je)) * par if (ddm .lt. 0.) ddm = 0. !**** COMPUTE potential increase in biomass - ddm ! VERSION for CLIMATE IMPACT ASSESSMENT: WITH CO2 FERTILIZATION !#### CALL ADJUSTBE to adjust be if ialpha = 1 if (ialpha .eq. 1) then call vegetation_adjust_energy_ratio(j, je, beadj, be, nucr, tx) ddm = beadj * par if (ddm .lt. 0.) ddm = 0. endif !#### CALL NUPTAKE & PUPTAKE - CALC N AND P UPTAKE call nutrient_nitrate_uptake(j, je, nucr(j, je), dm, flu, frar, g, ida, nn, snup) call nutrient_phosphorus_uptake(j, je, nucr(j, je), dm, flu, frar, g, ida, nn, spup) !**** CALC crop growth regulating factor: reg strsn = 1. strsp = 1. reg = amin1(real(ws(j, je), 4), real(ts, 4), real(strsn, 4), real(strsp, 4)) if (reg .lt. 0.) reg = 0. if (reg .gt. 1.) reg = 1. tsav(j, je) = tsav(j, je) + ts wsav(j, je) = wsav(j, je) + ws(j, je) !**** CALC biomass dm() & root weight fraction rwt() dm(j, je) = dm(j, je) + ddm * reg rwt(j, je) = (.4 - .2 * g(j, je)) !**** CALC harvest index under favourable conditions - hia() ! & real(dp) harvest index - hiad() f = g(j, je) / (g(j, je) + & exp(ald1(nucr(j, je)) - ald2(nucr(j, je)) * g(j, je))) ff = f - huharv(j, je) huharv(j, je) = f if (igro(j, je) .gt. 0) then if (g(j, je) .gt. 0.5) then swh(j, je) = swh(j, je) + actual swp(j, je) = swp(j, je) + potentl endif heat = 100. * (1. - (dlai(nucr(j, je)) - g(j, je))) !#### CALL SCURVE call vegetation_s_curve(x1, x2, x3, x4) hia(j, je) = hi(nucr(j, je)) * 100 * g(j, je) / & (100 * g(j, je) + exp(11.1 - 10. * g(j, je))) xinc = 100. * swh(j, je) / (swp(j, je) + 1.e-6) if (xinc .lt. 0.) xinc = 0. duma = xinc / (xinc + exp(x1 - x2 * xinc)) hiad(j, je) = hia(j, je) * duma else endif !**** CALC SUM STRESS DAYS sdw = sdw + (1. - ws(j, je)) * flu(j) * frar(j, je) sdt = sdt + (1. - ts) * flu(j) * frar(j, je) sdn = sdn + (1. - strsn) * flu(j) * frar(j, je) sdp = sdp + (1. - strsp) * flu(j) * frar(j, je) !**** COMPUTE LEAF AREA INDEX - alai() if (g(j, je) .le. dlai(nucr(j, je))) then if (alai(j, je) .gt. blai(nucr(j, je))) & alai(j, je) = blai(nucr(j, je)) xy = sqrt(reg) deltalai = ff * blai(nucr(j, je)) * & (1. - exp(5. * (alai(j, je) - blai(nucr(j, je))))) * xy alai(j, je) = alai(j, je) + deltalai if (alai(j, je) .gt. blai(nucr(j, je))) & alai(j, je) = blai(nucr(j, je)) olai(j, je) = alai(j, je) else alai(j, je) = 16. * olai(j, je) * (1. - g(j, je)) ** 2 if (alai(j, je) .gt. olai(j, je)) alai(j, je) = olai(j, je) end if end if !*********************************************************** END IF (G<=1) end if !*********************************************************** END IF (IGRO>=1) return end subroutine crop_growth !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine crop_yield_output(j, jek, k, ieap, ms, ndgro, tsav, wsav, ylda, is_cropland) !**** PURPOSE: Write crop yield for GRASS output !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! ieap = counter ! ndgro = number days of growth - calc. in crop ! tsav(j, je) = temp stress, sum for the growth period ! wsav(j, je) = water stress, sum for the growth period ! ylda(j, k) = crop yield, kg/ha ! >>>>> ! >>>>> STATIC PARAMETERS ! mts = mean temp stress for the growth period, in % ! mws = mean water stress for the growth period, in % ! myld = ylda(j, k)/100 - yield in kg/dt, integer ! xyld = crop yield, dt/ha ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters integer, intent(inout) :: ieap integer, intent(in) :: ms integer, intent(in) :: ndgro real(dp), dimension(:, :), intent(inout) :: tsav real(dp), dimension(:, :), intent(inout) :: wsav real(dp), dimension(:, :), intent(in) :: ylda logical, intent(in) :: is_cropland integer :: j, jek, k, mts, mws, myld real(dp) :: xyld if (is_cropland .and. k .lt. ms .and. ndgro .gt. 0) then xyld = ylda(j, k) / 100. + 0.5 myld = int(xyld) mws = int(wsav(j, jek) / ndgro * 100) mts = int(tsav(j, jek) / ndgro * 100) else myld = 0 mws = 0 mts = 0 end if ieap = ieap + 1 wsav(j, jek) = 0. tsav(j, jek) = 0. end subroutine crop_yield_output end module crop