module time use utilities, only : dp, log_info, log_debug, log_progress implicit none ! total number of days of simulation integer, save :: nDaysSim ! number of years of simulations integer, save :: nbyr = 10 ! current years (YYYY) real year integer, save :: iyr = 2000 ! integer, save :: ieap ! current year starting with 1 integer, save :: iy = 0 ! current month integer, save :: mo = 0 ! current day integer, save :: ida = 0 ! number of days in the current year integer, save :: nd = 0 ! number of days integer, save :: ndsum = 0 ! ?? integer, save :: nt = 0 integer, save, dimension(13) :: nc = (/ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366/) ! integer, save :: iday ! integer, save :: inday ! counter of the current day of simulation (not of current year) integer, save :: daycounter = 0 namelist / TIME_PARAMETERS / & iyr, & nbyr contains subroutine time_initialise use utilities, only : is_leap_year use input, only : get_config_fid integer i read(get_config_fid(), TIME_PARAMETERS) ! *** count total number of days in simulation nDaysSim = 0 do i = iyr, iyr + nbyr - 1 if (is_leap_year(i) ) then nDaysSim = nDaysSim + 366 else nDaysSim = nDaysSim + 365 end if end do end subroutine time_initialise subroutine time_day_length(daylen, ida, mb, pit, ylc, yls) integer, intent(in) :: ida integer, intent(in) :: mb real(dp), intent(in) :: pit real(dp), dimension(:), intent(in) :: ylc real(dp), dimension(:), intent(in) :: yls real(dp), dimension(:), intent(out) :: daylen integer k real(dp) ch, h, sd, xi !**** CALC day length (from clgen) do k = 1, mb xi = ida sd = .4102 * sin((xi - 80.25) / pit) ch = - yls(k) * tan(sd) / ylc(k) if (ch .lt. 1.) then if (ch .le. - 1.) then h = 3.1416 else h = acos(ch) end if else h = 0. end if daylen(k) = 7.72 * h end do end subroutine time_day_length subroutine time_process_years integer year, month call log_info("time_process_years", 'SWIM starts simulation') do year = 1, nbyr call time_initialise_year(year) call log_info('time_process_years', 'Processing year', int=iyr) do month = 1, 12 call time_process_month(month) end do call time_finish_year() end do end subroutine time_process_years ! ---------------------------------------------------------------------- subroutine time_initialise_year(year) ! ---------------------------------------------------------------------- use input, only : read_real_column use crop, only : & arylda, & avylda, & mcrdb, & iyrrot, & nrotyrs use subbasin, only : & bRunoffdat, & discharge_input_file_id, & obs_discharge, & nqobs integer year, jj, k, skip ! iy - current year ! mo - current month ! nt - =1 if 365 days, =0 if 366 days ! nd - number of days in the current year iy = year mo = 1 nt = 1 if (mod(iyr, 4) .eq. 0) nt = 0 ida = 1 nd = 366 - nt ndsum = ndsum + nd !#### INITIALIZE CROPS: CALL INITCROP do jj = 1, mcrdb avylda(iy, jj) = 0. arylda(iy, jj) = 0. end do if (bRunoffdat) then ! read all columns by index (+1 for date column) ! skipping back up before reading except for 1st do k = 1, nqobs skip = -nd if (k == 1) skip = 0 call read_real_column(discharge_input_file_id, array=obs_discharge(:nd, k), index=k+1, skip=skip) end do end if ! count year of rotation (1-3) if (iyrrot < nrotyrs) then iyrrot = iyrrot + 1 else iyrrot = 1 end if if (iy == 1) iyrrot = 0 end subroutine time_initialise_year ! ---------------------------------------------------------------------- subroutine time_finish_year() ! ---------------------------------------------------------------------- use output, only : output_nashsutcliffe_efficiency, output_year use crop, only : istyr use subbasin, only : & runs, & bRunoffDat, & obs_discharge if (bRunoffDat) then call output_nashsutcliffe_efficiency(obs_discharge(:nd, 1), runs(:nd), 1, istyr, iy) end if call output_year call log_debug("time_finish_year", 'Completed year =', int=iyr) return end subroutine time_finish_year ! ---------------------------------------------------------------------- subroutine time_process_month(month) ! ---------------------------------------------------------------------- use output, only : output_month use vegetation, only : ialpha use nutrient, only : & degNgrw, & degNsub, & degNsur, & degPsur, & retNgrw, & retNsub, & retNsur, & retPsur use snow, only : bSnowModule use soil, only : psp, rtn use output, only : nsb, nvsub use crop, only : icc, mfe use evapotranspiration, only : pit, radiation_switch use catchment, only : bSubcatch integer, intent(in) :: month mo = month !ls** compute number of days of the month if (mo .eq. 2) then inday = nc(mo + 1) - nc(mo) - nt else inday = nc(mo + 1) - nc(mo) endif do iday = 1, inday call time_process_day(month, iday) end do call output_month(iyr, month) if (month .eq. 12) iyr = iyr + 1 call log_debug("time_process_month", "Completed month =", int=month) return end subroutine time_process_month subroutine time_process_day(mo1, iday) !**** CALLED IN MAIN !**** THIS SUBROUTINE COMPUTES ONE DAY !**** !**** CALLED: IN MAIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! ida - current day ! ieap - index for GRASS output in subbasin (yield) ! ieapu - index for GRASS output in subbasin (annual sums) ! xxswind - soil water index for basin ! xwysb - water yield for basin ! xnflow() - N flows for a chosen hydrotop !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ use catchment, only : bSubcatch, da, da9, subcatch_id, wy use crop, only : icc, mfe use erosion, only : & chc, & chxk, & conn, & cpp, & er, & xcklsp, & xnorg, & xnorgp, & xporg, & xpsedp use evapotranspiration, only : & canev, & ecal, & evapotranspiration_radiation, & humi, & omega, & pit, & ra, & radiation_switch, & snoev, & tmn, & ylc, & yls use groundwater, only : & additionalGwUptake, & revapst, & xet use hydrotope, only : smm #ifdef with_netcdf use input, only : input_nc_read_climate #endif use management, only : & bWAM_Module, & management_external_supply, & management_is_transfer_subbasin, & management_total_demand, & wam_d, & wam_y use nutrient, only : & degNgrw, & degNsub, & degNsur, & degPsur, & retNgrw, & retNsub, & retNsur, & retPsur, & xnflow, & yon, & yph use output, only : & area_tot_glacier, & area_tot_snow, & depth_ave_glacier, & depth_ave_snow, & ieapu, & nsb, & nvsub, & output_day use reservoir, only : & bRsvModule, & bRsvHydrograph, & reservoir_is_operational, & reservoir_process, & reservoir_subbasin, & rsvSubbasin, & xwysb, & xxswind use river, only : & accf, & river_route, & river_route_add, & river_transfer, & varoute use snow, only : & bSnowModule, & sml, & snowVal, & tmx, & ieapg use soil, only : pr, psp, rtn, xqd use subbasin, only : & aff, & dart, & flu, & icodes, & ihouts, & inum1s, & inum2s, & mb, & mhyd, & neap, & nqobs, & obs_discharge, & obssb, & precip, & qtl, & runs, & runsubbasin, & sbar, & sbp, & sda, & smq, & smsq, & sub, & subbasin_initialise_subbasin, & subbasin_read_climate, & subouthyd, & subp, & sumcn, & susb, & tx, & wysb, & xeo, & xpercn, & xqi, & xsep, & xssf, & xssfn, & xswind, & xyno3 use subbasin, only : & xysp, & yd, & bRunoffdat use vegetation, only : daylen integer, intent(in) :: mo1, iday integer ii, iik, k integer inum1, inum2 integer idum integer icode, icodep, ihout logical :: bRoute = .true. !#### RESERVOIR MODULE #### !ls** compute day number in the year if (mo1 .gt. 2) then ida = iday + nc(mo1) - nt else ida = iday + nc(mo1) endif daycounter = daycounter + 1 area_tot_snow = 0. depth_ave_snow = 0. !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then ieapg = 1 area_tot_glacier = 0. depth_ave_glacier = 0. end if !########################### xxswind = 0. xwysb = 0. ieap = 1 ieapu = 1 do ii = 1, 20 xnflow(ii) = 0. end do do k = 1, nsb sub(k) = 0. end do !#### CALL (NC_) READCLI - to read climate data daily #ifdef with_netcdf call input_nc_read_climate(flu, humi, mb, ra, subp, tmn, tmx, tx) #else call subbasin_read_climate(humi, mb, ra, subp, tmn, tmx, tx) #endif call time_day_length(daylen, ida, mb, pit, ylc, yls) !### VA !*********** For Radiation Data generated by Hargreaves Samani if (radiation_switch > 0) then call evapotranspiration_radiation(ida, mb, tmx) end if !### VA !ls ndmo(mo) = ndmo(mo) + 1 !ls dtot = dtot + 1. snoev = 0. pr = 0. !########################################################### START ROUTING !#### CALL subbasin, route, transfer, add do idum = 1, mhyd icode = icodes(idum) ihout = ihouts(idum) inum1 = inum1s(idum) inum2 = inum2s(idum) if (icode .gt. 0) then select case (icode) case (1) ! SUBBASIN command !################################# !#### WATER MANAGEMENT MODULE #### !################################# if (bWAM_Module) then if (ihout == 1) then wam_y(daycounter) = iyr wam_d(daycounter) = ida end if if (management_is_transfer_subbasin(inum1) ) then ! Summarise total inflow only from external sources and transmission losses ! Current day inflows from external sources are added to subbasins %inflow ! NOTE: Adding to routing variables sda(2, j) takes place in subroutine subbasin call management_external_supply(inum1, daycounter, ida, iyr) ! calculate total water demand of water user(s) ! NOTE: nothing is removed here, just computation of total demand call management_total_demand(inum1, daycounter, ida, iyr) ! Depending on whether the subbasin is a headwater or not ! inputs and outputs are added/removed in subroutines: ! 'subbasin' or 'add' end if end if !########################### !########################### !#### RESERVOIR MODULE #### !########################### if (bRsvModule) then ! if actual subbasin is reservoir skip subbasin call ! TODO: include res_active function if (rsvSubbasin(inum1) == 0 ) then if (subcatch_id(inum1) .ne. 0) then !write(*,*) 'Inside swim.f95::: olai(1, 1) = ', olai(1, 1) call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd) !write(*,*) 'Inside swim.f95::: olai(1, 1) After = ', olai(1,1) end if else if (reservoir_is_operational(iyr, ida, inum1) ) then call subbasin_initialise_subbasin(canev, sml, xcklsp, xet, xnorg, xnorgp, xporg, xpsedp, xqd, yon, yph) call reservoir_subbasin(inum1, da, aff, dart, flu, neap, precip, qtl, sbp, sda, smq, smsq, snowval, sub, subp, sumcn, susb, tx, varoute, wysb, xeo, xpercn, xqi, xsep, xssf, xssfn, xswind, xyno3, xysp, yd, yon, yph) else if (subcatch_id(inum1) .ne. 0) then call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd) end if end if ! calculate subbasin but only if subbasin is listed in subcatch.def ! for convenient subsetting of model eg. by subcatchment else if (subcatch_id(inum1) .ne. 0) then call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd) end if !########################### case (2) ! ROUTE command bRoute = .true. !########################### !#### RESERVOIR MODULE #### !########################### if (bRsvModule) then ! If inum1 is a reservoir subbasin then... if (bRsvHydrograph(ihout) ) then ! reservoir input is stored in varoute(2, inum2) and varoute(8, inum2) ! skip routing ! modify varoute(2, inum1), varoute(8, inum1) in reservoir functions ! DO NOT call route(icode, ihout, inum1, inum2) if (reservoir_is_operational(iyr, ida, inum1) ) then call reservoir_process(ihout, inum1, inum2, ecal, humi, ida, iyr, mo, nc, omega, ra, subp, tx, varoute) bRoute = .false. end if else bRoute = .true. !call route(icode, ihout, inum1, inum2) end if end if !########################### if (bRoute) then call river_route(ihout, inum1, inum2, chc, chxk, conn, cpp, da9, dart, er, flu, ida, iy, iyr, revapst, runs, sbar, sub, susb, xysp, yd, yon, yph) end if case (3) ! not implemented ! do nothing case (4) ! not implemented ihout = ihouts (idum - 1) icodep = icodes (idum - 1) call river_transfer() ! not implemented case (5) ! ADD command call river_route_add(bRunoffdat, ihout, inum1, inum2, additionalGwUptake, bWAM_Module, daycounter, ida, iyr, mb, nqobs, obssb, obs_discharge, runs, subouthyd, inum1s, bRsvModule, reservoir_is_operational(iyr, ida, inum2), rsvSubbasin) case default ! do nothing end select end if ! (icode .gt. 0) end do ! do idum = 1, mhyd !########################################################### END ROUTING !**** Calc wy - water yield wy = sub(8) + sub(9) - sub(10) + sub(15) + 1.e-20 !**** Correction - sediment yield sub(21) = sub(21) / da9 !**** Calculation of monthly sums do iik = 1, nvsub smm(iik) = smm(iik) + sub(iik) end do !**** Calc monthly water discharge if (bRunoffdat) accf(7) = accf(7) + obs_discharge(ida, 1) accf(8) = accf(8) + runs(ida) call output_day(iyr, mo1, iday) call log_progress('time_process_day', ida, nd) call log_debug("time_process_day", "Completed day =", int=iday) return end subroutine time_process_day ! subroutine time_process_day(year, month, day) ! use hydrotope, only: hydrotope_process, hydrotope_initialise_vars_timestep, n_hyrotopes ! use subbasin, only: subbasin_process, subbasin_initialise_timestep, n_subbasins ! integer, intent(in) :: year, month, day ! integer hydrotope, subbasin ! do hydrotope = 1, n_hyrotopes ! call hydrotope_process(hydrotope) ! end do ! do subbasin = 1, n_subbasins ! call subbasin_process(subbasin) ! end do ! call catchment_process(year, month, day) ! call output_day(year, day, month) ! end subroutine time_process_day end module time