SUBROUTINE update_monthly_update_vars()

    USE pd_dates_vars
    USE pd_update_vars
    USE pd_report_vars

    USE pd_var_tables

    IMPLICIT NONE

    include "w1clig.inc"        ! precip
    include "p1werm.inc"        ! mntime (maximum # of time steps/day)
    include "w1wind.inc"        ! awu(mntime), awudmx
    include "w1pavg.inc"        ! awdair
    include "m1sim.inc"         ! ntstep (actual # of time steps/day)

    include "h1et.inc"          ! ah0drat (dryness ratio)

    include "h1db1.inc"         ! ahzsnd(s) snow depth in mm

    include "erosion/m2geo.inc" ! imax, jmax, ix, jy of simulation grid
    include "erosion/e2erod.inc"! egt, egtcs, egtss, egt10

    INTEGER :: s                ! local variable (subregion)
    INTEGER :: i,j              ! local loop variables
    INTEGER :: y                ! local loop variables
    INTEGER :: ngdpt            ! number of simulation grid datapoints

    INTEGER :: cnt              ! number of simulation grid datapoints
    REAL    :: sum_salt

    REAL :: we

    ! Threshold value for determining erosive wind energy (m/s)
    REAL, PARAMETER :: wind_energy_thresh = 8.0
    ! Threshold value for determining protective snow depth (mm)
    REAL, PARAMETER :: snow_depth_thresh = 20.0
    ! Threshold value for determining erosion loss and deposition regions
    REAL, PARAMETER :: eros_thresh = 0.025

    ! Flag to specify whether we have experienced an erosion event or not.
    ! It is set in the "Salt_Loss2" section and used in the "Trans_Cap"
    ! and "Sheltered" code sections.
    LOGICAL :: Have_Erosion

    Have_Erosion = .FALSE.   ! Initialize for each invocation of routine


    !variables summed for period
    monthly_update(Precipi)%val = monthly_update(Precipi)%val + awzdpt
    monthly_update(Precipi)%cnt = monthly_update(Precipi)%cnt + 1

    ngdpt = (imax-1) * (jmax-1)
    DO i = 1, imax-1 
       DO j = 1, jmax-1 
          monthly_update(Eros_loss)%val = monthly_update(Eros_loss)%val +     &
                                       egt(i,j)/ngdpt
          monthly_update(Salt_loss)%val = monthly_update(Salt_loss)%val +     &
                                       (egt(i,j) - egtss(i,j))/ngdpt
          monthly_update(Susp_loss)%val = monthly_update(Susp_loss)%val +     &
                                       egtss(i,j)/ngdpt
          monthly_update(PM10_loss)%val = monthly_update(PM10_loss)%val +     &
                                       egt10(i,j)/ngdpt
       END DO
    END DO
    monthly_update(Eros_loss)%cnt = monthly_update(Eros_loss)%cnt + 1
    monthly_update(Salt_loss)%cnt = monthly_update(Salt_loss)%cnt + 1
    monthly_update(Susp_loss)%cnt = monthly_update(Susp_loss)%cnt + 1
    monthly_update(PM10_loss)%cnt = monthly_update(PM10_loss)%cnt + 1

    ! Determine the saltation loss and area from region that generated it
    sum_salt = 0.0; cnt = 0
    DO i = 1, imax-1 
       DO j = 1, jmax-1 
          IF ((egt(i,j) - egtss(i,j)) < -eros_thresh) THEN
             sum_salt = sum_salt + egt(i,j) - egtss(i,j)
             cnt = cnt + 1
          END IF
       END DO
    END DO
    IF (cnt /= 0) THEN
       Have_Erosion = .TRUE.  !We have erosion occuring, set flag for use later
       monthly_update(Salt_loss2)%val =                                 &
            monthly_update(Salt_loss2)%val + sum_salt/cnt
       monthly_update(Salt_loss2)%cnt =                                 &
            monthly_update(Salt_loss2)%cnt + 1

       CALL run_ave (monthly_update(Salt_loss2_area), (ix*jy)*cnt/10000.0, 1)
       CALL run_ave (monthly_update(Salt_loss2_frac), REAL(cnt)/ngdpt, 1)
    END IF

    ! Determine the saltation deposition and area from region that generated it
    sum_salt = 0.0; cnt = 0
    DO i = 1, imax-1 
       DO j = 1, jmax-1 
          IF ((egt(i,j) - egtss(i,j)) > eros_thresh) THEN
             sum_salt = sum_salt + egt(i,j) - egtss(i,j)
             cnt = cnt + 1
          END IF
       END DO
    END DO
    IF (cnt /= 0) THEN
       monthly_update(Salt_dep2)%val =                                  &
            monthly_update(Salt_dep2)%val + sum_salt/cnt
       monthly_update(Salt_dep2)%cnt =                                  &
            monthly_update(Salt_dep2)%cnt + 1

       CALL run_ave (monthly_update(Salt_dep2_area), (ix*jy)*cnt/10000.0, 1)
       CALL run_ave (monthly_update(Salt_dep2_frac), REAL(cnt)/ngdpt, 1)
    END IF

IF (Have_Erosion == .TRUE.) THEN      !We have erosion, so compute TC, etc.

    ! Determine the region under saltation transport capacity
    ! (heavy erosion flux rates but zero net erosion and deposition)
    ! Note: Currently this assumes no barrier protection - LEW
    !       Also, the "sheltered area" is set equal to the "transport capacity area"
    cnt = 0
    DO i = 1, imax-1 
       DO j = 1, jmax-1 
          IF ( ABS((egt(i,j) - egtss(i,j))) < eros_thresh) THEN
             cnt = cnt + 1
          END IF
       END DO
    END DO
!   IF (cnt /= 0) THEN
       ! Note: We currently don't have a way of computing the Flux Rate over a grid area
       CALL run_ave (monthly_update(Trans_Cap_area), (ix*jy)*cnt/10000.0, 1)
       CALL run_ave (monthly_update(Trans_Cap_frac), REAL(cnt)/ngdpt, 1)
!Assume Sheltered region is same as Transport Region for now....
       CALL run_ave (monthly_update(Sheltered_area), (ix*jy)*cnt/10000.0, 1)
       CALL run_ave (monthly_update(Sheltered_frac), REAL(cnt)/ngdpt, 1)
!   END IF

END IF  !Have_Erosion flag

! Sum boundary losses  (ave value per boundary grid point)
    DO i = 0, imax 
      ! Note that egt contains creep+saltation not total soil loss on boundary
      monthly_update(Salt_1)%val = monthly_update(Salt_1)%val +         &
                                   egt(i,0)/(imax-1)
      monthly_update(Salt_3)%val = monthly_update(Salt_3)%val +         &
                                   egt(i,jmax)/(imax-1)
      monthly_update(Susp_1)%val = monthly_update(Susp_1)%val +         &
                                   egtss(i,0)/(imax-1)
      monthly_update(Susp_3)%val = monthly_update(Susp_3)%val +         &
                                   egtss(i,jmax)/(imax-1)
      monthly_update(PM10_1)%val = monthly_update(PM10_1)%val +         &
                                   egt10(i,0)/(imax-1)
      monthly_update(PM10_3)%val = monthly_update(PM10_3)%val +         &
                                   egt10(i,jmax)/(imax-1)
    END DO
    monthly_update(Salt_1)%cnt = monthly_update(Salt_1)%cnt + 1
    monthly_update(Salt_3)%cnt = monthly_update(Salt_3)%cnt + 1
    monthly_update(Susp_1)%cnt = monthly_update(Susp_1)%cnt + 1
    monthly_update(Susp_3)%cnt = monthly_update(Susp_3)%cnt + 1
    monthly_update(PM10_1)%cnt = monthly_update(PM10_1)%cnt + 1
    monthly_update(PM10_3)%cnt = monthly_update(PM10_3)%cnt + 1

    DO j = 0, jmax 
      ! Note that egt contains creep+saltation not total soil loss on boundary
      monthly_update(Salt_2)%val = monthly_update(Salt_2)%val +         &
                                   egt(0,j)/(jmax-1)
      monthly_update(Salt_4)%val = monthly_update(Salt_4)%val +         &
                                   egt(imax,j)/(jmax-1)
      monthly_update(Susp_2)%val = monthly_update(Susp_2)%val +         &
                                   egtss(0,j)/(jmax-1)
      monthly_update(Susp_4)%val = monthly_update(Susp_4)%val +         &
                                   egtss(imax,j)/(jmax-1)
      monthly_update(PM10_2)%val = monthly_update(PM10_2)%val +         &
                                   egt10(0,j)/(jmax-1)
      monthly_update(PM10_4)%val = monthly_update(PM10_4)%val +         &
                                   egt10(imax,j)/(jmax-1)
    END DO
    monthly_update(Salt_2)%cnt = monthly_update(Salt_2)%cnt + 1
    monthly_update(Salt_4)%cnt = monthly_update(Salt_4)%cnt + 1
    monthly_update(Susp_2)%cnt = monthly_update(Susp_2)%cnt + 1
    monthly_update(Susp_4)%cnt = monthly_update(Susp_4)%cnt + 1
    monthly_update(PM10_2)%cnt = monthly_update(PM10_2)%cnt + 1
    monthly_update(PM10_4)%cnt = monthly_update(PM10_4)%cnt + 1

    !variables running averaged for period
    we = 0.0
    IF (awudmx > wind_energy_thresh) THEN
       DO i = 1, ntstep
          IF (awu(i) > wind_energy_thresh) THEN
            we = we + 0.5*awdair*(awu(i)**2) * (awu(i) - wind_energy_thresh) *        &
               (86400.0/ntstep) * (0.001)    ! (s/day) and (J/kJ)
          END IF
       END DO
    END IF
    CALL run_ave(monthly_update(Wind_energy), we, 1)

    CALL run_ave(monthly_update(Dryness_ratio), ah0drat, 1)

    s = 1  !currently have only one subregion
    ! Note that the 20mm depth should be a global parameter
    ! It is currently stuck in erosion.for as a local parameter there
    IF (ahzsnd(s) > snow_depth_thresh) THEN
       CALL run_ave(monthly_update(Snow_cover), 1.0, 1)
    ELSE
       CALL run_ave(monthly_update(Snow_cover), 0.0, 1)
    END IF

END SUBROUTINE update_monthly_update_vars


SUBROUTINE update_monthly_report_vars(cur_month, cur_year, nrot_years)

    USE pd_dates_vars
    USE pd_update_vars
    USE pd_report_vars

    USE pd_var_tables

    IMPLICIT NONE

    INTEGER, INTENT (IN) :: cur_month
    INTEGER, INTENT (IN) :: cur_year
    INTEGER, INTENT (IN) :: nrot_years

    INTEGER :: i, y    ! local loop variables


    !variables averaged for reporting period
    DO y=0,nrot_years
       IF ( (y == 0) .or. (mod(cur_year-1,nrot_years)+1 == y) )  THEN
          DO i=1,Max_cli_vars
             CALL run_ave (monthly_report(i,cur_month,y),          &
                 monthly_update(i)%val, 1 )
          END DO
          DO i=Max_cli_vars+1,PM10_4
             CALL run_ave (monthly_report(i,cur_month,y),          &
                 monthly_update(i)%val, 1 )
          END DO
          CALL run_ave (monthly_report(Salt_loss2,cur_month,y),          &
                 monthly_update(Salt_loss2)%val, 1 )
          CALL run_ave (monthly_report(Salt_dep2,cur_month,y),          &
                 monthly_update(Salt_dep2)%val, 1 )
          CALL run_ave (monthly_report(Trans_Cap,cur_month,y),          &
                 monthly_update(Trans_Cap)%val, 1 )

          ! If we have saltating loss, add the area and fraction info
          IF (monthly_update(Salt_loss2)%cnt > 0) THEN 
             CALL run_ave (monthly_report(Salt_loss2_area,cur_month,y),          &
                 monthly_update(Salt_loss2_area)%val, 1 )
             CALL run_ave (monthly_report(Salt_loss2_frac,cur_month,y),          &
                 monthly_update(Salt_loss2_frac)%val, 1 )
             CALL run_ave (monthly_report(Salt_dep2_area,cur_month,y),          &
                 monthly_update(Salt_dep2_area)%val, 1 )
             CALL run_ave (monthly_report(Salt_dep2_frac,cur_month,y),          &
                 monthly_update(Salt_dep2_frac)%val, 1 )
             CALL run_ave (monthly_report(Trans_Cap_area,cur_month,y),          &
                 monthly_update(Trans_Cap_area)%val, 1 )
             CALL run_ave (monthly_report(Trans_Cap_frac,cur_month,y),          &
                 monthly_update(Trans_Cap_frac)%val, 1 )
             CALL run_ave (monthly_report(Sheltered_area,cur_month,y),          &
                 monthly_update(Sheltered_area)%val, 1 )
             CALL run_ave (monthly_report(Sheltered_frac,cur_month,y),          &
                 monthly_update(Sheltered_frac)%val, 1 )
          END IF
       END IF
    END DO

    ! reset update_vars
    DO i=1,Max_monthly_vars
       monthly_update(i)%cnt = 0
       monthly_update(i)%val = 0.0
    END DO
    monthly_dates(cur_month,0)%ey = monthly_dates(cur_month,0)%ey + 1

END SUBROUTINE update_monthly_report_vars