! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ##
!              ### Irrigation event with irrigation amount calculation based on soil type ###
! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ## ! ##

! Working in the interface mo_nwp_sfc_interface.f90

USE mo_time_config,         ONLY: time_config ! JL
USE mtime,                  ONLY: datetime, newDatetime, getDayOfYearFromDateTime, deallocateDatetime  ! JL

    REAL(wp) :: s_fc (nproma) ! JL
    REAL(wp) :: s_pwp (nproma) ! JL

    ! JL IRRIGATION VARIABLES
    REAL(wp) :: crop_mo_t (nproma) !Fraction
    REAL(wp) :: crop_irr_t (nproma) ! JL Fraction post-flooding or irrigated croplands
    INTEGER :: crop_irr_tile (nproma) ! JL Tile post-flooding or irrigated croplands
    REAL(wp) :: lu_class_fraction (nproma) ! Landuse class fraction

    INTEGER :: end_hour !UTC

    INTEGER :: irr_on_day, &
               irr_on_hour, &
               irr_on, &
               irr_start_hour, & ! UTC
               irr_num_hours, &
               irr_freq ! Days
    REAL(wp) :: irr_amount (nproma) ! for FC calculations
    INTEGER ::  m_styp_irr ! JL soil type index
    REAL(wp) :: d_frac !  depleted fraction range between 0.5 and 0.8
    REAL(wp) :: root_d_irr ! root depth for irrigation [mm]
    REAL(wp) :: constants_irrigation (nproma)
    INTEGER :: day_difference, days_since_start

    REAL(wp) :: IRRIGATION_CHANNEL (nproma)         ! JL Irrigation channel in mm

    TYPE(datetime)                       :: mtime_now
    INTEGER  :: mtime_now_day
    INTEGER  :: mtime_now_year
    CHARACTER(4) :: year_string
    TYPE(datetime)                       :: mtime_irri_start
    TYPE(datetime)                       :: mtime_irri_end
    INTEGER  :: mtime_irri_start_day, mtime_irri_end_day


    ! IRR EVENT JL
    irr_start_hour = 5
    irr_num_hours = 5 !24h irrigation for the first free simulations
    irr_freq = 12
    d_frac = 0.65
    root_d_irr = 810 ! mm until 5 soil layer in TERRA
    ! Dates current, start, end
    mtime_now = time_config%tc_current_date
    mtime_now_day = getDayOfYearFromDateTime(mtime_now)
    mtime_now_year = mtime_now%date%year ! Year can change..
    ! Convert the integer year to a string
    WRITE(year_string, '(I4)') mtime_now_year
    mtime_irri_start = newDatetime(TRIM(year_string) // "-05-01T00:00:00")
    mtime_irri_end = newDatetime(TRIM(year_string) // "-09-01T00:00:00")
    mtime_irri_start_day = getDayOfYearFromDateTime(mtime_irri_start)
    mtime_irri_end_day = getDayOfYearFromDateTime(mtime_irri_end)
    ! Irrigation frequency
    day_difference = mtime_now_day - mtime_irri_start_day
    days_since_start = MOD(day_difference, irr_freq)
    ! irr_on_day
    IF (mtime_irri_start_day.lt.mtime_irri_end_day) THEN
            IF (mtime_now_day.ge.mtime_irri_start_day .AND. mtime_now_day.lt.mtime_irri_end_day) THEN
                    IF (days_since_start == 0) THEN
                            irr_on_day = 1
                    ELSE
                            irr_on_day = 0
                    ENDIF
            ELSE
                    irr_on_day = 0
            ENDIF
    ELSE
            irr_on_day = 0
    ENDIF

    ! irr_on_hour
    end_hour = irr_start_hour + irr_num_hours
    if( end_hour.gt.23) end_hour = end_hour - 24
    IF (irr_on_day == 1) THEN
            IF (mtime_now%time%hour >= irr_start_hour .AND. mtime_now%time%hour < end_hour) THEN
                    irr_on_hour = 1
                          !PRINT *, "mtime_now HOURS is ", mtime_now%time%hour

            ELSE
                    irr_on_hour = 0
            ENDIF
    ELSE
            irr_on_hour = 0
    ENDIF
    ! irrigation event irr_on
    IF (irr_on_hour.EQ.1 .AND. irr_on_day.EQ.1) THEN
            irr_on = 1
    ELSE
            irr_on = 0
    ENDIF



          ! JL ADDING IRRIGATION WATER TO THE SYSTEM
          ! Soil type is the same for 
          m_styp_irr = soiltyp_t(ic) ! Current soil type
          s_fc(ic) = cfcap (m_styp_irr)              ! field capacity
          s_pwp(ic) = cpwp (m_styp_irr)              ! plant wilting point
          ! Irrigation amount based on texture characteristics, depletion frac and root depth
          irr_amount(ic) = (s_fc(ic) - s_pwp(ic)) * d_frac * root_d_irr ! [mm]

          constants_irrigation(ic) = irr_amount(ic)*0.000277778_wp/irr_num_hours ! [mm/s]

          IF (irr_on.eq.1 .AND. lc_class_t(ic) == ext_data%atm%i_lc_crop_irrig) THEN
             ! Here irrigation fraction (as tile) from GlobCover is considered
             ! tcall_sfc_jg is my dt in TERRA, units in seconds
            IRRIGATION_CHANNEL(ic) = tcall_sfc_jg*crop_irr_tile(ic)*constants_irrigation(ic) ! In mm
          ELSE
            IRRIGATION_CHANNEL(ic) = 0._wp
          ENDIF

        ! JL Adding IRRIGATION_CHANNEL to prr_gsp, in kg/m2s if there is no precipitation
          IF (lc_class_t(ic) == ext_data%atm%i_lc_crop_irrig .AND. prr_gsp_t(ic).LE.0.00001_wp) THEN
            prr_gsp_t(ic) = prr_gsp_t(ic) + IRRIGATION_CHANNEL(ic)
          ELSE
            prr_gsp_t(ic) = prr_gsp_t(ic)
          ENDIF