program simple_graphics use module_netcdf_io use kwm_date_utilities implicit none character(len=256) :: validation_data = "bondville.val" character(len=256) :: ncfile1 character(len=256) :: ncfile2 character(len=256) :: executable_name integer, parameter :: dt = 30 integer, parameter :: iunit = 10 integer :: ierr character(len=12) :: startdate ! YYYYmmddhhmm character(len=12) :: enddate ! YYYYmmddhhmm integer :: ncid1 integer :: ncid2 integer :: idt integer :: iarg integer :: numarg character(len=256) :: arg integer, external :: iargc ! Default values, overridden as necessary by command-line arguments. startdate = "0000000000"//"00" enddate = "0000000000"//"00" ncfile1 = "-" ncfile2 = "-" open(iunit, file=trim(validation_data), status='old', form='formatted', action='read') numarg = iargc() iarg = 1 do while (iarg <= numarg) call getarg(iarg, arg) print*, 'arg = "'//trim(arg)//'"' iarg = iarg + 1 if ( (trim(arg) == "-help") .or. (trim(arg) == "-h") .or. (trim(arg) == "--help") ) then call getarg(0, executable_name) write(*,'(/,"Usage: ",A, " [-startdate ] [-enddate ] []")') trim(executable_name) write(*,'(/," Where is an optional starting date/time (YYYYMMDDHH).")') write(*,'( " The default is to start plotting from the first available time.")') write(*,'( " is an optional ending date/time (YYYYMMDDHH).")') write(*,'( " The default is to end plotting at the last available time.")') write(*,'( " is the NetCDF output file from the 1d driver.")') write(*,'( " is an optional the NetCDF output file from a second run of the 1D")') write(*,'( " driver, for comparison.")') write(*,*) stop else if (trim(arg) == "-startdate") then call getarg(iarg, startdate) iarg = iarg + 1 do idt = 1,12 if (startdate(idt:idt) == " ") startdate(idt:idt) = "0" enddo print*, 'startdate = ', startdate else if (trim(arg) == "-enddate") then call getarg(iarg, enddate) iarg = iarg + 1 do idt = 1,12 if (enddate(idt:idt) == " ") enddate(idt:idt) = "0" enddo print*, 'enddate = ', enddate else if (ncfile1 == "-") then ncfile1 = arg print*, 'ncfile1 = ', trim(ncfile1) else if (ncfile2 == "-") then ncfile2 = arg print*, 'ncfile2 = ', trim(ncfile2) endif enddo ncid1 = -1 if (ncfile1 /= "-") then call output_open_for_read(trim(ncfile1), ncid1) endif ncid2 = -1 if (ncfile2 /= "-") then call output_open_for_read(trim(ncfile2), ncid2) endif print*, 'ncid1, ncid2 = ', ncid1, ncid2 call gopks(6,0) call gopwk(1,20,1) call gacwk(1) call pcsetc("FC", "~") call pcseti("FN", 21) call gaseti("LTY", 1) call gscr(1, 0, 1.0, 1.0, 1.0) call gscr(1, 1, 0.0, 0.0, 0.0) call gscr(1, 2, 0.0, 0.0, 1.0) call gscr(1, 3, 0.8, 0.0, 0.0) call gscr(1, 4, 1.0, 0.0, 1.0) call gscr(1, 255, 0.8, 0.8, 0.8) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "H", "SHEAT", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "LE", "ETA", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "CH", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "RES", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "SHEAT", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "PRCP", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "SSOIL", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "EDIR", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ETT", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ET:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ET:Level2", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ET:Level3", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ET:Level4", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ESNOW", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "DRIP", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "DEW", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "ALBEDO", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "Z0", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "SHDFAC", startdate, enddate, dt, 2, 3, 4) call make_graphs(-1, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "-", "EMISSI", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "SKN_IRT", "T1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_02", "STC:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_04", "STC:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_08", "STC:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_16", "STC:Level1", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_32", "STC:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_16", "STC:Level2", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_32", "STC:Level2", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T_64", "STC:Level3", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "SM_05", "SMC:Level1", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "SM_20", "SMC:Level2", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "SM_60", "SMC:Level3", startdate, enddate, dt, 2, 3, 4) call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "LW_IN", "LWDN", startdate, enddate, dt, 2, 3, 4) !KWM call make_graphs(iunit, ncid1, ncid2, trim(ncfile1), trim(ncfile2), "T", "SFCTMP", startdate, enddate, dt, 2, 3, 4) call gdawk(1) call gclwk(1) call gclks() end program simple_graphics !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine make_graphs(iunit, ncid1, ncid2, ncfile1, ncfile2, vald_name, noah_name, startdate, enddate, dt, icol1, icol2, icol3) use kwm_date_utilities use module_netcdf_io implicit none integer, intent(in) :: iunit, ncid1, ncid2 character(len=*), intent(in) :: vald_name, noah_name character(len=*), intent(in) :: ncfile1, ncfile2 character(len=12), intent(inout) :: startdate, enddate integer, intent(in) :: dt integer, intent(in) :: icol1, icol2, icol3 real, pointer, dimension(:) :: avgd_vald, avgd_noah1, avgd_noah2, xavgd real, pointer, dimension(:) :: vald, xval real, pointer, dimension(:) :: noah1, x1 real, pointer, dimension(:) :: noah2, x2 integer :: lvl, indx integer :: noah_startindex, noah_endindex integer :: itime character(len=256) :: text integer :: idtm real :: xmin_full real :: xmax_full character(len=256) :: units1 character(len=256) :: description1 character(len=256) :: units2 character(len=256) :: description2 type frame_definition real :: left real :: right real :: bottom real :: top end type frame_definition type(frame_definition) :: box = frame_definition (0.15, 0.85, 0.3, 0.7) allocate(avgd_vald(0:1439)) allocate(avgd_noah1(0:1439)) allocate(avgd_noah2(0:1439)) allocate(xavgd(0:1439)) do itime = 0, 1439 xavgd(itime) = real(itime) enddo print*, 'Make graphs: "'//vald_name//'" "'//noah_name//'"' if (ncid1 > 0) then call read_nc_single_field(ncid1, noah_name, startdate, enddate, noah1, x1, units1, description1) endif if (ncid2 > 0) then call read_nc_single_field(ncid2, noah_name, startdate, enddate, noah2, x2, units2, description2) endif if (iunit > 0) then call read_vald_single_field(iunit, vald_name, vald, startdate, enddate, xval, dt) endif ! ! Number of minutes: ! print*, 'startdate = ', startdate print*, 'enddate = ', enddate call geth_idts(enddate, startdate, idtm) ! Our X-coordinate is in minutes. xmin_full = 0.0 xmax_full = real(idtm) call set(0., 1., 0., 1., 0., 1., 0., 1., 1) write(text,'("Time series from ",A4,2("-",A2)," to ",A4,2("-",A2))') & startdate(1:4), startdate(5:6), startdate(7:8), enddate(1:4), enddate(5:6), enddate(7:8) call pchiqu(0.5, 0.85, trim(text), 0.015, 0., 0.) if (iunit > 0) then call pcseti("CC", icol1) call pchiqu(0.33, 0.80, "Obs "//vald_name, 0.012, 0.0, 0.0) endif if (ncid1 > 0) then call pcseti("CC", icol2) call strrep(units1, "{", "~S~") call strrep(units1, "}", "~N~") call pchiqu(0.66, 0.80, "Model: "//ncfile1//": "//noah_name//" ("//trim(units1)//")", 0.012, 0.0, 0.0) call pchiqu(0.66, 0.78, trim(description1), 0.009, 0.0, 0.0) endif if (ncid2 > 0) then call pcseti("CC", icol3) call strrep(units2, "{", "~S~") call strrep(units2, "}", "~N~") call pchiqu(0.66, 0.74, "Model: "//ncfile2//": "//noah_name//" ("//trim(units2)//")", 0.012, 0.0, 0.0) call pchiqu(0.66, 0.72, trim(description2), 0.009, 0.0, 0.0) endif call pcseti("CC", 1) call pchiqu(box%left, box%bottom, "~V-3Q~"//startdate(1:10), 0.010, 0., 0.) call pchiqu(box%right, box%bottom, "~V-3Q~"//enddate(1:10), 0.010, 0., 0.) ! Draw model and obs call gasetr("XLO", 1.E5) ! Send the default X-axis labels somewhere else entirely. call drawlines(xmin_full, xmax_full, xval, vald, x1, noah1, x2, noah2, icol1, icol2, icol3) call frame() ! Compute average diurnal cycles. ! call average_diurnal_cycle(startdate, size(vald), dt, vald, noah1, noah2, avgd_vald, avgd_noah1, avgd_noah2) if (associated(vald)) then call average_diurnal_cycle(startdate, vald, xval, avgd_vald) else deallocate(avgd_vald) nullify(avgd_vald) endif if (associated(noah1)) then call average_diurnal_cycle(startdate, noah1, x1, avgd_noah1) else deallocate(avgd_noah1) nullify(avgd_noah1) endif if (associated(noah2)) then call average_diurnal_cycle(startdate, noah2, x2, avgd_noah2) else deallocate(avgd_noah2) nullify(avgd_noah2) endif ! Draw average diurnal cycle call set(0., 1., 0., 1., 0., 1., 0., 1., 1) call pchiqu(0.5, 0.88, "Average diurnal cycle", 0.015, 0., 0.) write(text,'("averaged over the period ",A4,2("-",A2)," to ",A4,2("-",A2))') & startdate(1:4), startdate(5:6), startdate(7:8), enddate(1:4), enddate(5:6), enddate(7:8) call pchiqu(0.5, 0.85, trim(text), 0.015, 0., 0.) if (iunit>0) then call pcseti("CC", icol1) call pchiqu(0.33, 0.80, "Obs "//vald_name, 0.012, 0.0, 0.0) endif if (ncid1>0) then call pcseti("CC", icol2) call strrep(units1, "{", "~S~") call strrep(units1, "}", "~N~") call pchiqu(0.66, 0.80, "Model: "//ncfile1//": "//noah_name//" ("//trim(units1)//")", 0.012, 0.0, 0.0) call pchiqu(0.66, 0.78, trim(description1), 0.009, 0.0, 0.0) endif if (ncid2>0) then call pcseti("CC", icol3) call strrep(units2, "{", "~S~") call strrep(units2, "}", "~N~") call pchiqu(0.66, 0.74, "Model: "//ncfile2//": "//noah_name//" ("//trim(units2)//")", 0.012, 0.0, 0.0) call pchiqu(0.66, 0.72, trim(description2), 0.009, 0.0, 0.0) endif call pcseti("CC", 1) call gasetc("XLF", "(I3)") call gasetr("XLO", 20.0) ! Back to default position for x-axis labels call drawlines(0.0, 1439.0, xavgd, avgd_vald, xavgd, avgd_noah1, xavgd, avgd_noah2, icol1, icol2, icol3) call frame() if (associated(noah1)) then deallocate(noah1) nullify(noah1) endif if (associated(avgd_noah1)) then deallocate(avgd_noah1) nullify(avgd_noah1) endif if (associated(noah2)) then deallocate(noah2) nullify(noah2) endif if (associated(avgd_noah2)) then deallocate(avgd_noah2) nullify(avgd_noah2) endif if (associated(vald)) then deallocate(vald) nullify(vald) endif if (associated(avgd_vald)) then deallocate(avgd_vald) nullify(avgd_vald) endif contains !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine read_vald_single_field(iunit, name, vald, startdate, enddate, x, dt) use kwm_date_utilities implicit none integer, intent(in) :: iunit character(len=*), intent(in) :: name real, dimension(:), pointer :: vald real, dimension(:), pointer :: x character(len=12), intent(in) :: startdate character(len=12), intent(in) :: enddate integer, intent(in) :: dt character(len=256), parameter :: rdfmt = '(A4,4(1x,A2), 31(F12.5))' character(len=12) :: readdate real :: w_speed real :: w_dir real :: t real :: rh real :: p real :: rg real :: par_in real :: par_out real :: rnet real :: ghf real :: rain real :: wet real :: skn_irt real :: t_02 real :: t_04 real :: t_08 real :: t_16 real :: t_32 real :: t_64 real :: u_bar real :: eddyuw real :: uprim2 real :: vprim2 real :: wprim2 real :: H real :: le real :: co2 real :: lw_in real :: sm_05 real :: sm_20 real :: sm_60 integer :: ierr integer :: itime integer :: idt integer :: idtm integer :: ntimes integer :: i call geth_idts(enddate, startdate, idt) ntimes = idt/dt + 1 allocate(vald(ntimes)) vald = -1.E36 allocate(x(ntimes)) x = -1.E36 itime = 0 READLOOP : do read(iunit, FMT=trim(rdfmt), iostat=ierr) readdate(1:4), readdate(5:6), readdate(7:8), readdate(9:10), readdate(11:12), & w_speed, w_dir, & T, RH, P, Rg, Par_in, Par_out, & rnet, GHF, RAIN, wet, SKN_IRT, & T_02, T_04, T_08, T_16, T_32, T_64, & u_bar, eddyuw, uprim2, vprim2, wprim2, & H, LE, CO2, LW_in, & sm_05, sm_20, sm_60 if (ierr < 0) then write(*,'("Hit the end of file.")') write(*,'("Unexpected error exit.")') stop endif if (ierr /= 0) then write(*,'("Error reading from data file.")') stop endif if (readdate < startdate) then cycle READLOOP else if (readdate > enddate) then exit READLOOP else itime = itime + 1 call geth_idts(readdate(1:12), startdate(1:12), idtm) x(itime) = real(idtm) if (name == "T") then if (T > -6998) vald(itime) = t + 273.15 else if (name == "H") then if (h > -6998) vald(itime) = h else if (name == "LE") then if (le > -6998) vald(itime) = le else if (name == "SKN_IRT") then if (skn_irt > -6998) vald(itime) = skn_irt + 273.15 else if (name == "T_02") then if (t_02 > -6998) vald(itime) = t_02 + 273.15 else if (name == "T_04") then if (t_04 > -6998) vald(itime) = t_04 + 273.15 else if (name == "T_08") then if (t_08 > -6998) vald(itime) = t_08 + 273.15 else if (name == "T_16") then if (t_16 > -6998) vald(itime) = t_16 + 273.15 else if (name == "T_32") then if (t_32 > -6998) vald(itime) = t_32 + 273.15 else if (name == "T_64") then if (t_64 > -6998) vald(itime) = t_64 + 273.15 else if (name == "LW_IN") then if (lw_in > -6998) vald(itime) = lw_in else if (name == "SM_05") then if (sm_05 > -6998) vald(itime) = sm_05 else if (name == "SM_20") then if (sm_20 > -6998) vald(itime) = sm_20 else if (name == "SM_60") then if (sm_60 > -6998) vald(itime) = sm_60 else print*, 'Variable name = "'//name//'"' stop "Variable not found in validation file" endif if (readdate >= enddate) exit READLOOP endif enddo READLOOP rewind(iunit) end subroutine read_vald_single_field !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine read_nc_single_field(ncid, name, startdate, enddate, noah, x, units, description) use module_netcdf_io use kwm_date_utilities implicit none integer, intent(in) :: ncid character(len=*), intent(in) :: name character(len=12), intent(inout) :: startdate character(len=12), intent(inout) :: enddate real, dimension(:), pointer :: noah real, dimension(:), pointer :: x character(len=256), intent(out) :: units character(len=256), intent(out) :: description integer :: noah_startindex integer :: noah_endindex character(len=12), pointer, dimension(:) :: noah_times integer :: itime integer :: lvl integer :: indx real, pointer, dimension(:,:) :: noah2d integer :: i integer :: idtm call get_from_netcdf_output(ncid, "Times", noah_times) if (startdate == "000000000000") then noah_startindex = 1 startdate = noah_times(noah_startindex) write(*,'("Setting start date to ", A)') startdate else noah_startindex = -1 do itime = 1, size(noah_times) if (noah_times(itime) == startdate) then noah_startindex = itime endif enddo endif if (enddate == "000000000000") then noah_endindex = size(noah_times) enddate = noah_times(noah_endindex) write(*,'("Setting end date to ", A)') enddate else noah_endindex = -1 do itime = noah_startindex, size(noah_times) if (noah_times(itime) == enddate) then noah_endindex = itime exit endif enddo endif indx = index(name,":Level") if (indx > 0) then read( name(indx+6:),*) lvl call get_from_netcdf_output(ncid, name(1:indx-1), noah2d, units, description, noah_startindex, noah_endindex) allocate(noah(size(noah2d,2))) noah = noah2d(lvl,:) deallocate(noah2d) nullify(noah2d) else call get_from_netcdf_output(ncid, name, noah, units, description, noah_startindex, noah_endindex) endif allocate(x(size(noah))) do i = 1, size(x) call geth_idts(noah_times(i+noah_startindex-1), startdate, idtm) x(i) = real(idtm) enddo deallocate(noah_times) nullify(noah_times) end subroutine read_nc_single_field !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine average_diurnal_cycle(startdate, val, xval, avgd_val) use kwm_date_utilities implicit none character(len=12) :: startdate real, dimension(:), pointer :: val real, dimension(:), pointer :: xval real, dimension(0:1439) :: avgd_val real :: xbf, xaf, vbf, vaf, fraction, vn integer :: indx, i, j, offm integer, dimension(0:1439) :: avgd_count avgd_count = 0 avgd_val = 0.0 call geth_idts(startdate, startdate(1:8)//"0000", offm) do i = 1, size(val)-1 xbf = xval(i) xaf = xval(i+1) vbf = val(i) vaf = val(i+1) if ( (xbf > -1.E25) .and. (xaf > -1.E25) .and. (vbf > -1.E25) .and. (vaf > -1.E25) ) then do j = nint(xbf), nint(xaf) fraction = (j-xbf)/(xaf-xbf) vn = ( vbf * (1.-fraction) ) + ( vaf * fraction ) indx = mod(j+offm,1440) avgd_val(indx) = avgd_val(indx) + vn avgd_count(indx) = avgd_count(indx) + 1 enddo endif enddo where (avgd_count > 0) avgd_val = avgd_val / real(avgd_count) elsewhere avgd_val = -1.E36 endwhere end subroutine average_diurnal_cycle #if 0 subroutine average_diurnal_cycle(startdate, ntimes, dt, val, noah1, noah2, avgd_val, avgd_noah1, avgd_noah2) use kwm_date_utilities implicit none character(len=12), intent(in) :: startdate integer, intent(in) :: ntimes integer, intent(in) :: dt real, dimension(:), pointer :: val real, dimension(:), pointer :: noah1 real, dimension(:), pointer :: noah2 real, dimension(:), pointer :: avgd_val real, dimension(:), pointer :: avgd_noah1 real, dimension(:), pointer :: avgd_noah2 real, allocatable, dimension(:) :: avgd_count integer :: itime character(len=12) :: nowdate integer :: idt if (associated(val) .and. associated(noah1) .and. associated(noah2)) then allocate(avgd_val(1440/dt)) allocate(avgd_noah1(1440/dt)) allocate(avgd_noah2(1440/dt)) allocate(avgd_count(1440/dt)) avgd_val = 0. avgd_noah1 = 0. avgd_noah2 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (val(itime) > -1.E25) then avgd_val(idt) = avgd_val(idt) + val(itime) avgd_noah1(idt) = avgd_noah1(idt) + noah1(itime) avgd_noah2(idt) = avgd_noah2(idt) + noah2(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_val = avgd_val / avgd_count avgd_noah1 = avgd_noah1 / avgd_count avgd_noah2 = avgd_noah2 / avgd_count elsewhere avgd_val = -1.E36 avgd_noah1 = -1.E36 avgd_noah2 = -1.E36 endwhere elseif (associated(val) .and. associated(noah1)) then allocate(avgd_val(1440/dt)) allocate(avgd_noah1(1440/dt)) allocate(avgd_count(1440/dt)) avgd_val = 0. avgd_noah1 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (val(itime) > -1.E25) then avgd_val(idt) = avgd_val(idt) + val(itime) avgd_noah1(idt) = avgd_noah1(idt) + noah1(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_val = avgd_val / avgd_count avgd_noah1 = avgd_noah1 / avgd_count elsewhere avgd_val = -1.E36 avgd_noah1 = -1.E36 endwhere elseif (associated(val) .and. associated(noah2)) then allocate(avgd_val(1440/dt)) allocate(avgd_noah2(1440/dt)) allocate(avgd_count(1440/dt)) avgd_val = 0. avgd_noah2 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (val(itime) > -1.E25) then avgd_val(idt) = avgd_val(idt) + val(itime) avgd_noah2(idt) = avgd_noah2(idt) + noah2(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_val = avgd_val / avgd_count avgd_noah2 = avgd_noah2 / avgd_count elsewhere avgd_val = -1.E36 avgd_noah2 = -1.E36 endwhere elseif (associated(noah1) .and. associated(noah2)) then allocate(avgd_noah1(1440/dt)) allocate(avgd_noah2(1440/dt)) allocate(avgd_count(1440/dt)) avgd_noah1 = 0. avgd_noah2 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (noah1(itime) > -1.E25) then avgd_noah1(idt) = avgd_noah1(idt) + noah1(itime) avgd_noah2(idt) = avgd_noah2(idt) + noah2(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_noah1 = avgd_noah1 / avgd_count avgd_noah2 = avgd_noah2 / avgd_count elsewhere avgd_noah1 = -1.E36 avgd_noah2 = -1.E36 endwhere elseif (associated(val)) then allocate(avgd_val(1440/dt)) allocate(avgd_count(1440/dt)) avgd_val = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (val(itime) > -1.E25) then avgd_val(idt) = avgd_val(idt) + val(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_val = avgd_val / avgd_count elsewhere avgd_val = -1.E36 endwhere elseif (associated(noah1)) then allocate(avgd_noah1(1440/dt)) allocate(avgd_count(1440/dt)) avgd_noah1 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (noah1(itime) > -1.E25) then avgd_noah1(idt) = avgd_noah1(idt) + noah1(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_noah1 = avgd_noah1 / avgd_count elsewhere avgd_noah1 = -1.E36 endwhere elseif (associated(noah2)) then allocate(avgd_noah2(1440/dt)) allocate(avgd_count(1440/dt)) avgd_noah2 = 0. avgd_count = 0. do itime = 1, ntimes call geth_newdate(nowdate, startdate, (itime-1)*dt) call geth_idts(nowdate, nowdate(1:8)//"0000", idt) idt = 1+(idt/dt) if (noah2(itime) > -1.E25) then avgd_noah2(idt) = avgd_noah2(idt) + noah2(itime) avgd_count(idt) = avgd_count(idt) + 1 endif enddo where (avgd_count > 0) avgd_noah2 = avgd_noah2 / avgd_count elsewhere avgd_noah2 = -1.E36 endwhere endif end subroutine average_diurnal_cycle #endif !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine drawlines(xmin, xmax, x1, val1, x2, val2, x3, val3, icol1, icol2, icol3) implicit none real, intent(in) :: xmin, xmax integer, intent(in) :: icol1, icol2, icol3 real, dimension(:), pointer:: val1, val2, val3 real, dimension(:), pointer:: x1, x2, x3 real :: ymin, ymax integer :: expn, nrange ymin = 1.E36 ymax = -1.E36 ! Compute reasonable min/max values for the Y axis if (associated(val1)) then ymin = min(minval(val1, mask=(val1>-1.E25)), ymin) ymax = max(maxval(val1, mask=(val1>-1.E25)), ymax) endif if (associated(val2)) then ymin = min(minval(val2, mask=(val2>-1.E25)), ymin) ymax = max(maxval(val2, mask=(val2>-1.E25)), ymax) endif if (associated(val3)) then ymin = min(minval(val3, mask=(val3>-1.E25)), ymin) ymax = max(maxval(val3, mask=(val3>-1.E25)), ymax) endif if (ymin < ymax) then call find_good_range(ymin, ymax, expn, nrange) else nrange = 2 expn = 0 ymin = ymin - 1.0 ymax = ymax + 1.0 endif call set (box%left, box%right, box%bottom, box%top, xmin, xmax, ymin, ymax, 1) call gasetc("YLF", '(F9.4)') call gasetc("XLF", '(F20.4)') call periml(0, 4, nrange, 0) if ( (ymin < 0) .and. (ymax > 0)) then call gsplci(255) call line(0.0, 0.0, 1.0, 0.0) call gsplci(1) endif if (associated(val1)) then call drawline(x1, val1, size(x1), icol1) endif if (associated(val2)) then call drawline(x2, val2, size(x2), icol2) endif if (associated(val3)) then call drawline(x3, val3, size(x3), icol3) endif end subroutine drawlines !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ end subroutine make_graphs !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine drawline(x, val, ntimes, icol) integer, intent(in) :: ntimes real, intent(in), dimension(ntimes) :: x real, intent(in), dimension(ntimes) :: val integer, intent(in) :: icol logical :: pd integer :: itimes pd = .FALSE. call gsplci(icol) do itime = 1, ntimes if (val(itime) > -1.E25) then if (.not. pd) then call frstpt(x(itime), val(itime)) pd = .TRUE. else call vector(x(itime), val(itime)) endif else call plotif(0., 0., 2) pd = .FALSE. endif enddo call plotif(0., 0., 2) pd = .FALSE. call gsplci(1) end subroutine drawline !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine find_good_range(ymin, ymax, expn, irange) ! ! Given an initial and (minimum and maximum of some field), ! return "good round numbers" for and . ! is the number of steps between our good round and good round ! is the power-of-ten scaling applied to get our round numbers. ! ! There are steps of size 10.0** between output values of and ! implicit none real, intent(inout) :: ymin real, intent(inout) :: ymax integer, intent(out) :: expn integer, intent(out) :: irange integer :: iymin integer :: iymax expn = floor(log10(ymax-ymin)) - 1 iymin = floor ( (10.0**(real(-expn))) * ymin ) iymax = ceiling ( (10.0**(real(-expn))) * ymax ) irange = (iymax-iymin) ! If we've got too many steps between iymin and iymax, ! then shift to the next power-of-ten for our scaling. if (irange > 20) then iymin = floor (real(iymin)*0.1) iymax = ceiling(real(iymax)*0.1) expn = expn + 1 irange = (iymax-iymin) endif ymin = real(iymin) * 10.0**(real(expn)) ymax = real(iymax) * 10.0**(real(expn)) end subroutine find_good_range !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------ subroutine strrep(string, rep1, rep2) ! In string 1, replace expression rep1 with expression rep2. ! if rep2 is shorter than rep1, fill the end of the string ! with blanks implicit none character(len=*) :: string, rep1, rep2 integer :: idx, inlen, len1, len2 do inlen = len(string) len1 = len(rep1) len2 = len(rep2) idx = index(string, rep1)-1 if (idx == -1) then return else string = string(1:idx)// rep2 // & string((idx+len1+1):inlen) // & " " endif enddo end subroutine strrep !------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------