SUBROUTINE update_yrly_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 yrly_update(Precipi)%val = yrly_update(Precipi)%val + awzdpt yrly_update(Precipi)%cnt = yrly_update(Precipi)%cnt + 1 ngdpt = (imax-1) * (jmax-1) DO i = 1, imax-1 DO j = 1, jmax-1 yrly_update(Eros_loss)%val = yrly_update(Eros_loss)%val + & egt(i,j)/ngdpt yrly_update(Salt_loss)%val = yrly_update(Salt_loss)%val + & (egt(i,j) - egtss(i,j))/ngdpt yrly_update(Susp_loss)%val = yrly_update(Susp_loss)%val + & egtss(i,j)/ngdpt yrly_update(PM10_loss)%val = yrly_update(PM10_loss)%val + & egt10(i,j)/ngdpt END DO END DO yrly_update(Eros_loss)%cnt = yrly_update(Eros_loss)%cnt + 1 yrly_update(Salt_loss)%cnt = yrly_update(Salt_loss)%cnt + 1 yrly_update(Susp_loss)%cnt = yrly_update(Susp_loss)%cnt + 1 yrly_update(PM10_loss)%cnt = yrly_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 !write(UNIT=6,FMT="(A,i5,A,i5)",ADVANCE="NO") 'ngdpt: ',ngdpt, ' Loss: ', cnt Have_Erosion = .TRUE. !We have erosion occuring, set flag for use later yrly_update(Salt_loss2)%val = & yrly_update(Salt_loss2)%val + sum_salt/cnt yrly_update(Salt_loss2)%cnt = & yrly_update(Salt_loss2)%cnt + 1 CALL run_ave (yrly_update(Salt_loss2_area), (ix*jy)*cnt/10000.0, 1) CALL run_ave (yrly_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 (Have_Erosion == .TRUE.) THEN !We have erosion so compute TC, etc. !write(UNIT=6,FMT="(A,i5)",ADVANCE="NO") " Dep:", cnt !ENDIF IF (cnt /= 0) THEN yrly_update(Salt_dep2)%val = & yrly_update(Salt_dep2)%val + sum_salt/cnt yrly_update(Salt_dep2)%cnt = yrly_update(Salt_dep2)%cnt + 1 CALL run_ave (yrly_update(Salt_dep2_area), (ix*jy)*cnt/10000.0, 1) CALL run_ave (yrly_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 !write(UNIT=6,FMT="(A,i5)",ADVANCE="YES") " TC: ", cnt !IF (cnt == 0) THEN !write(6,*) 'Have Erosion, but cnt is zero', Have_Erosion, cnt !ENDIF ! IF (cnt /= 0) THEN ! Note: We currently don't have a way of computing the Flux Rate over a grid area CALL run_ave (yrly_update(Trans_Cap_area), (ix*jy)*cnt/10000.0, 1) CALL run_ave (yrly_update(Trans_Cap_frac), REAL(cnt)/ngdpt, 1) !Assume Sheltered region is same as Transport Region for now as mentioned in Note above. CALL run_ave (yrly_update(Sheltered_area), (ix*jy)*cnt/10000.0, 1) CALL run_ave (yrly_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 yrly_update(Salt_1)%val = yrly_update(Salt_1)%val + egt(i,0)/(imax-1) yrly_update(Salt_3)%val = yrly_update(Salt_3)%val + egt(i,jmax)/(imax-1) yrly_update(Susp_1)%val = yrly_update(Susp_1)%val + egtss(i,0)/(imax-1) yrly_update(Susp_3)%val = yrly_update(Susp_3)%val + egtss(i,jmax)/(imax-1) yrly_update(PM10_1)%val = yrly_update(PM10_1)%val + egt10(i,0)/(imax-1) yrly_update(PM10_3)%val = yrly_update(PM10_3)%val + egt10(i,jmax)/(imax-1) END DO yrly_update(Salt_1)%cnt = yrly_update(Salt_1)%cnt + 1 yrly_update(Salt_3)%cnt = yrly_update(Salt_3)%cnt + 1 yrly_update(Susp_1)%cnt = yrly_update(Susp_1)%cnt + 1 yrly_update(Susp_3)%cnt = yrly_update(Susp_3)%cnt + 1 yrly_update(PM10_1)%cnt = yrly_update(PM10_1)%cnt + 1 yrly_update(PM10_3)%cnt = yrly_update(PM10_3)%cnt + 1 DO j = 0, jmax ! Note that egt contains creep+saltation not total soil loss on boundary yrly_update(Salt_2)%val = yrly_update(Salt_2)%val + egt(0,j)/(jmax-1) yrly_update(Salt_4)%val = yrly_update(Salt_4)%val + egt(imax,j)/(jmax-1) yrly_update(Susp_2)%val = yrly_update(Susp_2)%val + egtss(0,j)/(jmax-1) yrly_update(Susp_4)%val = yrly_update(Susp_4)%val + egtss(imax,j)/(jmax-1) yrly_update(PM10_2)%val = yrly_update(PM10_2)%val + egt10(0,j)/(jmax-1) yrly_update(PM10_4)%val = yrly_update(PM10_4)%val + egt10(imax,j)/(jmax-1) END DO yrly_update(Salt_2)%cnt = yrly_update(Salt_2)%cnt + 1 yrly_update(Salt_4)%cnt = yrly_update(Salt_4)%cnt + 1 yrly_update(Susp_2)%cnt = yrly_update(Susp_2)%cnt + 1 yrly_update(Susp_4)%cnt = yrly_update(Susp_4)%cnt + 1 yrly_update(PM10_2)%cnt = yrly_update(PM10_2)%cnt + 1 yrly_update(PM10_4)%cnt = yrly_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(yrly_update(Wind_energy), we, 1) CALL run_ave(yrly_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(yrly_update(Snow_cover), 1.0, 1) ELSE CALL run_ave(yrly_update(Snow_cover), 0.0, 1) END IF END SUBROUTINE update_yrly_update_vars SUBROUTINE update_yrly_report_vars(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) :: nrot_years INTEGER, INTENT (IN) :: cur_year 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 CALL run_ave (yrly_report(Precipi,y), & yrly_update(Precipi)%val, 1 ) CALL run_ave (yrly_report(Wind_energy,y), & yrly_update(Wind_energy)%val, 1 ) CALL run_ave (yrly_report(Dryness_ratio,y), & yrly_update(Dryness_ratio)%val, 1 ) CALL run_ave (yrly_report(Snow_cover,y), & yrly_update(Snow_cover)%val, 1 ) CALL run_ave (yrly_report(Eros_loss,y), & yrly_update(Eros_loss)%val, 1 ) CALL run_ave (yrly_report(Salt_loss,y), & yrly_update(Salt_loss)%val, 1 ) CALL run_ave (yrly_report(Susp_loss,y), & yrly_update(Susp_loss)%val, 1 ) CALL run_ave (yrly_report(PM10_loss,y), & yrly_update(PM10_loss)%val, 1 ) CALL run_ave (yrly_report(Salt_loss2,y), & yrly_update(Salt_loss2)%val, 1 ) CALL run_ave (yrly_report(Salt_dep2,y), & yrly_update(Salt_dep2)%val, 1 ) CALL run_ave (yrly_report(Trans_Cap,y), & yrly_update(Trans_Cap)%val, 1 ) ! If we have saltating loss, add the area and fraction info IF (yrly_update(Salt_loss2)%cnt > 0) THEN CALL run_ave (yrly_report(Salt_loss2_area,y), & yrly_update(Salt_loss2_area)%val, 1 ) !write(6,*) '1yrly_update(Salt_loss2_frac): ',y, yrly_update(Salt_loss2_frac)%val, & ! yrly_update(Salt_loss2_frac)%cnt !write(6,*) '1yrly_report(Salt_loss2_frac): ',y, yrly_report(Salt_loss2_frac,y)%val, & ! yrly_report(Salt_loss2_frac,y)%cnt CALL run_ave (yrly_report(Salt_loss2_frac,y), & yrly_update(Salt_loss2_frac)%val, 1 ) !write(6,*) '2yrly_report(Salt_loss2_frac): ',y, yrly_report(Salt_loss2_frac,y)%val, & ! yrly_report(Salt_loss2_frac,y)%cnt CALL run_ave (yrly_report(Salt_dep2_area,y), & yrly_update(Salt_dep2_area)%val, 1 ) CALL run_ave (yrly_report(Salt_dep2_frac,y), & yrly_update(Salt_dep2_frac)%val, 1 ) CALL run_ave (yrly_report(Trans_Cap_area,y), & yrly_update(Trans_Cap_area)%val, 1 ) !write(6,*) '1yrly_update(Trans_Cap_frac): ',y, yrly_update(Trans_Cap_frac)%val, & ! yrly_update(Trans_Cap_frac)%cnt !write(6,*) '1yrly_report(Trans_Cap_frac): ',y, yrly_report(Trans_Cap_frac,y)%val, & ! yrly_report(Trans_Cap_frac,y)%cnt CALL run_ave (yrly_report(Trans_Cap_frac,y), & yrly_update(Trans_Cap_frac)%val, 1 ) !write(6,*) '2yrly_report(Trans_Cap_frac): ',y, yrly_report(Trans_Cap_frac,y)%val, & ! yrly_report(Trans_Cap_frac,y)%cnt CALL run_ave (yrly_report(Sheltered_area,y), & yrly_update(Sheltered_area)%val, 1 ) CALL run_ave (yrly_report(Sheltered_frac,y), & yrly_update(Sheltered_frac)%val, 1 ) END IF ! per boundary grid pt CALL run_ave (yrly_report(Salt_1,y), & yrly_update(Salt_1)%val, 1 ) CALL run_ave (yrly_report(Salt_2,y), & yrly_update(Salt_2)%val, 1 ) CALL run_ave (yrly_report(Salt_3,y), & yrly_update(Salt_3)%val, 1 ) CALL run_ave (yrly_report(Salt_4,y), & yrly_update(Salt_4)%val, 1 ) CALL run_ave (yrly_report(Susp_1,y), & yrly_update(Susp_1)%val, 1 ) CALL run_ave (yrly_report(Susp_2,y), & yrly_update(Susp_2)%val, 1 ) CALL run_ave (yrly_report(Susp_3,y), & yrly_update(Susp_3)%val, 1 ) CALL run_ave (yrly_report(Susp_4,y), & yrly_update(Susp_4)%val, 1 ) CALL run_ave (yrly_report(PM10_1,y), & yrly_update(PM10_1)%val, 1 ) CALL run_ave (yrly_report(PM10_2,y), & yrly_update(PM10_2)%val, 1 ) CALL run_ave (yrly_report(PM10_3,y), & yrly_update(PM10_3)%val, 1 ) CALL run_ave (yrly_report(PM10_4,y), & yrly_update(PM10_4)%val, 1 ) END IF END DO ! reset update_vars DO i=1,Max_yrly_vars yrly_update(i)%cnt = 0 yrly_update(i)%val = 0.0 END DO yrly_dates(0)%ey = yrly_dates(0)%ey + 1 END SUBROUTINE update_yrly_report_vars