!     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, &
    crop_input_file, &
    crop_management_input_file


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, 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(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, 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
    !      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, &
      nutrient_on, &
      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
      if (nutrient_on) call nutrient_nitrate_uptake( &
          j, je, nucr(j, je), dm, flu, frar, g, ida, nn, snup)
      if (nutrient_on) 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
