diff --git a/SYNOP/STATS/fortran-programs/plots_for_one_station.f95 b/SYNOP/STATS/fortran-programs/plots_for_one_station.f95 index 0f5cfe0937a2eb8e5185e0acbfd423db400f5352..0b362641888b582feed8bd2bc93ee0383e65d62c 100644 --- a/SYNOP/STATS/fortran-programs/plots_for_one_station.f95 +++ b/SYNOP/STATS/fortran-programs/plots_for_one_station.f95 @@ -210,6 +210,7 @@ program plots_for_one_station integer,parameter :: n_colors=13 ! number of colours for GrADS output integer :: color(n_colors) ! colours for GrADS output + integer :: num_int_arr ! fix for Warning : Legacy Extension: REAL array index at (1) integer :: irec ! record number for GrADS output ! @@ -250,8 +251,8 @@ program plots_for_one_station ! f=miss open(unit=1,form='formatted',file=sim_file,status='old') - do while(.true.) -11 read(1,*,err=11,end=12)yyear,mmonth,dday,hhour,f1 + do while(.true.) +11 read(1,*,err=11,end=12) yyear,mmonth,dday,hhour,f1 if(l_round_6h)then call round_6h(yyear,mmonth,dday,hhour,year,month,day,hour) ! Rounding to nearest 6-h? else @@ -363,9 +364,13 @@ program plots_for_one_station write(1,'(A11)')'set cmark 3' write(1,'(A15)')'set digsiz 0.07' if(nyear.eq.1)then - write(1,'(A11,I2)')'set ccolor ',color(nint(1.+n_colors)/2.) + num_int_arr=nint(1.+n_colors)/2. + write(1,'(A11,I2)')'set ccolor ',color(num_int_arr) + !write(1,'(A11,I2)')'set ccolor ',color(nint(1.+n_colors)/2.) else - write(1,'(A11,I2)')'set ccolor ',color(nint(1+(i-1.)/(nyear-1.)*(n_colors-1.))) + num_int_arr=nint(1+(i-1.)/(nyear-1.)*(n_colors-1.)) + write(1,'(A11,I2)')'set ccolor ',color(num_int_arr) + !write(1,'(A11,I2)')'set ccolor ',color(nint(1+(i-1.)/(nyear-1.)*(n_colors-1.))) endif write(1,'(A7,I2,A1)')'d &0(z=',i,')' enddo diff --git a/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95 b/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95 index 5f129d6c6ad88f24eba9ee89859a2d7f6b43d0c2..46a8a608af26ddb7e8e1cd32fc327b9aac536c28 100644 --- a/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95 +++ b/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95 @@ -175,7 +175,7 @@ ! open(unit=1,form='formatted',status='old',file=infile) read(1,*) ! the first line in infile is a header - write(*,*)'Opened infile' + write(*,*)' Opened infile: ', infile do while(.true.) ! Loop continued until end of file reached if(l_code_in_char)then @@ -224,7 +224,7 @@ enddo 3 continue close(1) - write(*,*)'Read infile. n_line,n_station=',n_line,n_station + write(*,*)' Read infile: # of lines, # of stations: ',n_line,n_station ! if(n_line.ne.n_station*nhour_used)then ! write(*,*)'Something wrong !!!' diff --git a/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95~ b/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95~ deleted file mode 100644 index 743dc271290a7c694efb0f0624e0d6222cef2d5a..0000000000000000000000000000000000000000 --- a/SYNOP/STATS/fortran-programs/rank_histogram_summary_statistics.f95~ +++ /dev/null @@ -1,380 +0,0 @@ - program rank_histogram_summary_statistics -! -! Calculation of all-station rank histogram summary statistics: -! -! 1) distribution of p-values (in nbins_p bins) -! 2) Average frequencies in NQUANT+1 = 100 quantile bins -! (the quantile bin number is hard-coded!) -! -! The calculation is made for the UTC hours defined by -! the namelist parameters HOUR1, HOUR2, HOUR_STEP. It is assumed -! that the input file (INFILE) only includes data for these -! hours. -! -! Jouni Räisänen, University of Helsinki, August 2023 -! -!------------------------------------------------------------------------- -! -! INPUT FILE: -! -! infile: a text file including the rank histogram statistics for all stations -! -! The first line is a header. The other lines include: -! -! 1) station code (as a 3-character string, if l_code_in_char =.true. Else as an integer) -! 2) longitude -! 3) latitude -! 4) UTC hour (24 for 'all-UTC' statistics -! 5) total number of observations for this UTC hour -! 6) Mean squared deviation (MSD) from a flat quantile frequency histogram -! 7) p-value for MSD (= probability for getting at least as large MSD by chance) -! 8-107) frequencies in quantile space: 0-1%, 1-2% ... 98-99%, 99-100% -! -! Order of data assumed in the input file: -! -! line 1: header -! line 2: station 1, HOUR1 -! line 3: station 1, HOUR1 + HOUR_STEP -! line n = 1 + (HOUR2-HOUR1)/HOUR_STEP + 1 : station 2, HOUR2 -! line n+1: station2, HOUR1 -! etc. -!-------------------------------------------------------------------------- -! -! OUTPUT FILES: -! -! outfile_p_values: p-values at individual stations (This far: only the all-UTC p-values!) -! This file can be used for plotting the p-values on a map in python -! -! The first line is a header for eventual conversion to ODB. The other lines include -! -! 1) Station code as integer. If the original code was a 3-char string (GRUAN), -! the ordinal number of the station is written (for plotting in python) -! 2) longitude -! 3) latitude -! 4) p-value for MSD -!---------------------------------------------------------------------------- -! outfile_p_freq: frequency histogram of p-values (for different UTC hours separately) -! This file can be used for plotting a histogram of the p-value frequencies in python! -! -! The first line is a header. The others include -! -! 1) UTC hour -! 2-nbins_p+1) frequencies of p-values in nbins_p bins -!----------------------------------------------------------------------------- -! outfile_q = histogram of all-station-mean quantile bin frequencies (for different UTC hours separately) -! -! The first line is a header. The others include -! -! 1) UTC hour -! 2-101) quantile bin frequencies for 0-1%, 1-2% ... 99-100%. -!------------------------------------------------------------------------------ -! outfile_grads (only if grads_out=.true.): p-value frequencies and quantile bin frequencies in GrADS binary format -!------------------------------------------------------------------------------ -! text_output: frequencies of selected small p-values and some other data in free-text format, for each UTC hour separately. -!----------------------------------------------------------------------------------------------------------------------------- -! -! NAMELIST PARAMETERS: -! -! infile: see above -! outfile_p_values: see above -! outfile_p_freq: see above -! outfile_q: see above -! outfile_grads: see above -! text_output: see above -! nbins_p: number of p-value bins in outfile_p_freq -! grads_out: if .true. (.false.), outfile_grads is written (not written) -! -! hour1,hour2,hour_step: UTC hours in output = hour1, hour1+hour_step ... hour2 -! -! l_code_in_char: if .true., the station codes in infile are assumed to be 3-char string (else: integer) -! miss: missing value code. Default = -9.99e6. All values f with |f| >= |miss| are treated as missing. - - IMPLICIT NONE - INTEGER :: I,J ! loop variables - INTEGER,PARAMETER :: NHOUR=24 ! number of stations - INTEGER,PARAMETER :: NQUANT=99 ! number of quantiles from 1% to 99% - INTEGER,PARAMETER :: NPMAX=100 ! maximum bin number for frequency diagram of p-values - INTEGER,PARAMETER :: NSTATMAX=10000 ! maximum number of stations - REAL :: frequency_q1(nquant+1) ! quantile frequencies at a single station and UTC hour - REAL :: frequency_q(nquant+1,nhour+1) !all-station mean quantile frequencies for each UTC hour - REAL :: frequency_p(npmax,nhour+1) !frequencies of p-values for each UTC hours - REAL :: max_freq_p,max_freq_q ! maximum of p-value and quantile frequencies (for GrADS only) - REAL :: frequency_p01(nhour+1) ! frequency of p_values <= 0.1 % - REAL :: frequency_p1(nhour+1) ! same, <= 1 % - REAL :: frequency_p5(nhour+1) ! same, <= 5 % - character*2 :: number_of_bin ! two-digit code for quantiles for outfile_q and outfile_p_values headers - character*16 :: frequency_value ! frequency written in F16.6 format - character*160 :: infile,outfile_p_freq,outfile_p_values,outfile_q ! input and output files (see above) - character*160 :: outfile_grads,text_output ! output files (see above) - INTEGER,PARAMETER :: L_dataline=1700 ! maximum length of lines in text files - character*1700 :: dataline,headerline,emptyline ! string variables for input and output - - INTEGER :: nbins_p ! number of p-value bins in outfile_p_freq - INTEGER :: hour1,hour2,hour_step ! input / output UTC hours (see above) - INTEGER :: hour ! UTC hour - INTEGER :: n_line ! line count for lines read from infile - INTEGER :: n_station ! station count - INTEGER :: pbin ! index for p-value bin - - LOGICAL :: HOUR_USED(nhour+1) ! .true. for hours included in (hour1, hour1+hour_step,...,hour2) - INTEGER :: nhour_used ! number of UTC hours used - INTEGER :: station1_int ! the station code read from infile (RHARM, FMI) - CHARACTER*3 :: station1_char ! the station code read from infile (GRUAN) - LOGICAL :: l_code_in_char ! if .true., station codes in 3-char strings assumed - REAL :: longitude1,latitude1 !longitude and latitude - INTEGER :: station(nstatmax) !station codes of all stations - REAL :: longitude(nstatmax),latitude(nstatmax) ! station longitudes and latitudes - REAL :: p_value1, msd ! p-value and MSD read from infile - REAL :: p_value(nhour+1,nstatmax) ! p-values at individual stations for map - REAL :: total_number ! number of observations read from infile - - LOGICAL :: GRADS_OUT ! Output also as GrADS binaries, just for testing - INTEGER :: IREC ! record number for GrADS output -! - REAL :: MISS ! Missing value code - INTEGER :: valid_stations(nhour+1) ! stations with valid data - - NAMELIST/param/infile,outfile_p_values,outfile_p_freq,outfile_q,& - outfile_grads,text_output,& - nbins_p,grads_out,& - hour1,hour2,hour_step,& - l_code_in_char,miss - - miss=-9.99e6 ! default for missing value - READ(*,NML=PARAM) - valid_stations=0 -! -! Which UTC hours are used and what is their total number? -! - hour_used=.false. - nhour_used=0. - do i=hour1,hour2,hour_step - hour_used(i+1)=.true. - nhour_used=nhour_used+1 - enddo -! -! Empty data line -! - do i=1,L_dataline - emptyline(i:i)=' ' - enddo -!---------------------------------------------------------------------- -! -! Initialize counters etc. -! - frequency_p=0. - frequency_q=0. - frequency_p01=0. - frequency_p1=0. - frequency_p5=0. - n_station=0 - n_line=0 - p_value=miss -! -! Read the contents of the input file -! - open(unit=1,form='formatted',status='old',file=infile) - read(1,*) ! the first line in infile is a header - write(*,*)'Opened infile' - - do while(.true.) ! Loop continued until end of file reached - if(l_code_in_char)then - 1 read(1,*,err=1,end=3)& - station1_char,longitude1,latitude1,hour,total_number,msd,p_value1,(frequency_q1(i),i=1,nquant+1) - else - 2 read(1,*,err=2,end=3)& - station1_int,longitude1,latitude1,hour,total_number,msd,p_value1,(frequency_q1(i),i=1,nquant+1) - endif -! write(*,*)'n_station,Hour step:',n_station,hour_step - if(hour_used(hour+1))then ! Only use the statistics for the UTC hours define by hour1, hour2, hour_step - n_line=n_line+1 - if(mod(n_line,nhour_used).eq.1)then ! This assumes that all the required UTC hours are always included in infile - n_station=n_station+1 - if(l_code_in_char)then - station(n_station)=n_station - else - station(n_station)=station1_int - endif - longitude(n_station)=longitude1 - latitude(n_station)=latitude1 - endif - if(total_number.gt.0.and.abs(frequency_q1(1)).lt.abs(miss))then ! Only include stations with some valid data - p_value(hour+1,n_station)=p_value1 ! for map of p-values - valid_stations(hour+1)=valid_stations(hour+1)+1 - do i=1,nquant+1 - frequency_q(i,hour+1)=frequency_q(i,hour+1)+frequency_q1(i) ! update the quantile bin frequencies - enddo - pbin=min(1+(p_value1*nbins_p),real(nbins_p)) - frequency_p(pbin,hour+1)=frequency_p(pbin,hour+1)+1 ! update the p-value bin frequencies -! -! Frequencies of small p-values: -! - if(p_value1.le.0.001.and.abs(p_value1).lt.abs(miss))then - frequency_p01(hour+1)=frequency_p01(hour+1)+1. - endif - if(p_value1.le.0.01.and.abs(p_value1).lt.abs(miss))then - frequency_p1(hour+1)=frequency_p1(hour+1)+1. - endif - if(p_value1.le.0.05.and.abs(p_value1).lt.abs(miss))then - frequency_p5(hour+1)=frequency_p5(hour+1)+1. - endif - endif - - endif - enddo -3 continue - close(1) - write(*,*)'Read infile. n_line,n_station=',n_line,n_station - -! if(n_line.ne.n_station*nhour_used)then -! write(*,*)'Something wrong !!!' -! stop -! endif -!--------------------------------------------------------------------- -! Divide all the frequencies by the number of stations -! - do hour=hour1,hour2,hour_step - do i=1,nquant+1 - frequency_q(i,hour+1)=frequency_q(i,hour+1)/valid_stations(hour+1) - enddo - do i=1,nbins_p - frequency_p(i,hour+1)=frequency_p(i,hour+1)/valid_stations(hour+1) - enddo - frequency_p01(hour+1)=frequency_p01(hour+1)/valid_stations(hour+1) - frequency_p1(hour+1)=frequency_p1(hour+1)/valid_stations(hour+1) - frequency_p5(hour+1)=frequency_p5(hour+1)/valid_stations(hour+1) - enddo -! -! Find the normalized maximum frequencies for GrADS axis limits -! - max_freq_p=0 - max_freq_q=0 - do i=1,nbins_p - if(frequency_p(i,nhour+1).gt.max_freq_p/nbins_p)then - max_freq_p=frequency_p(i,nhour+1)*nbins_p - endif - enddo - max_freq_p=1.02*max_freq_p - - do i=1,nquant+1 - if(frequency_q(i,nhour+1).gt.max_freq_q/(nquant+1))then - max_freq_q=frequency_q(i,nhour+1)*(nquant+1) - endif - enddo - max_freq_q=1.02*max_freq_q -! -!--------------------------------------------------------------------- -! -! Write the text output file: -! 1) number of stations -! 2) frequencies of selected small p-values - - open(unit=1,form='formatted',file=text_output) - write(1,'(A20,I3)')'Number of stations: ',n_station - write(1,'(A24,I3)')'Number of p-value bins: ',nbins_p - if(grads_out)then - write(1,'(A23,F6.3)')'Axis_limit_p_for_GrADS ',max_freq_p - write(1,'(A23,F6.3)')'Axis_limit_q_for_GrADS ',max_freq_q - endif - do hour=hour1,hour2,hour_step - write(1,'(A22,I2.2,A6,F5.3)')'Frequency, p <= 0.001 ',hour,' UTC: ',& - frequency_p01(hour+1) - write(1,'(A22,I2.2,A6,F5.3)')'Frequency, p <= 0.01 ',hour,' UTC: ',& - frequency_p1(hour+1) - write(1,'(A22,I2.2,A6,F5.3)')'Frequency, p <= 0.05 ',hour,' UTC: ',& - frequency_p5(hour+1) - enddo - close(1) -! -! Write the average quantile bin frequencies and p-value bin frequencies -! to ODB compatible text files - -!-------------------------------------------------------------------------- -! -! Open the quantile bin frequency output file and write its first line. -! - open(unit=2,form='formatted',status='unknown',file=outfile_q) - headerline=emptyline - headerline='hour@hdr:integer ' - do j=0,nquant - write(number_of_bin,'(i2.2)')j - headerline=trim(headerline)//' f'//number_of_bin//'@body:real' - enddo - write(2,*)trim(headerline) -! -! Write the data lines for each UTC hour used -! - do hour=hour1,hour2,hour_step - j=hour+1 - write(dataline,'(I6)')& - hour - do i=1,nquant+1 - write(frequency_value,'(F16.6)')frequency_q(i,j) - dataline=trim(dataline)//frequency_value - enddo - write(2,*)trim(dataline) - enddo - close(2) -! -! Write the file giving the p-values for individual stations -! - open(unit=2,form='formatted',status='unknown',file=outfile_p_values) - headerline=' station@hdr:integer longitude@hdr:real latitude@hdr:real p_value@body:real' - write(2,*)trim(headerline) - do i=1,n_station - write(2,'(7X,I6,3F16.6)')& ! this far: just fwrite the station ordinal number, in case - ! the station code was a character string - station(i),longitude(i),latitude(i),p_value(nhour+1,i) - enddo - close(2) -! -! Open the p_value frequency output file and write its first line. -! - open(unit=2,form='formatted',status='unknown',file=outfile_p_freq) - headerline=emptyline - headerline='hour@hdr:integer ' - do j=1,nbins_p - write(number_of_bin,'(i2.2)')j-1 - headerline=trim(headerline)//' p'//number_of_bin//'@body:real' - enddo - write(2,*)trim(headerline) -! -! Write the data lines for each UTC hour used -! - do hour=hour1,hour2,hour_step - j=hour+1 - write(dataline,'(I6)')& - hour - do i=1,nbins_p - write(frequency_value,'(F16.6)')frequency_p(i,j) - dataline=trim(dataline)//frequency_value - enddo - write(2,*)trim(dataline) - enddo - close(2) -! -!------------------------------------------------------------------------------- -! -! Open the file for GrADS output and write its contents, for easier visualisation? -! - if(grads_out)then - open(unit=11,form='unformatted',file=outfile_grads,access='DIRECT',& - recl=4,status='unknown') -! - do hour=hour1,hour2,hour_step - j=(hour-hour1)/hour_step+1 - do i=1,nbins_p - irec=(j-1)*(nbins_p+nquant+1)+i - write(11,rec=irec)frequency_p(i,hour+1) - enddo - do i=1,nquant+1 - irec=(j-1)*(nbins_p+nquant+1)+nbins_p+i - write(11,rec=irec)frequency_q(i,hour+1) - write(*,*)hour,i,frequency_q(i,hour+1) - enddo - enddo - close(11) - - endif ! if (grads_out) - - END program rank_histogram_summary_statistics diff --git a/SYNOP/STATS/fortran-programs/rank_histograms_one_station.exe b/SYNOP/STATS/fortran-programs/rank_histograms_one_station.exe deleted file mode 100755 index 68af780d96994881fe44a9f93119733f5a2562a7..0000000000000000000000000000000000000000 Binary files a/SYNOP/STATS/fortran-programs/rank_histograms_one_station.exe and /dev/null differ diff --git a/SYNOP/STATS/fortran-programs/rank_histograms_one_station.f95 b/SYNOP/STATS/fortran-programs/rank_histograms_one_station.f95 index c73168c0dff7d4dda4d61f01838b755aab9438e9..2540ae9b0a0018c1acfb1c82f815b5a670be2f41 100644 --- a/SYNOP/STATS/fortran-programs/rank_histograms_one_station.f95 +++ b/SYNOP/STATS/fortran-programs/rank_histograms_one_station.f95 @@ -186,12 +186,10 @@ if(l_code_in_char)then ! The err=1 specifier ensures that header lines are skipped 1 read(1,*,err=1,end=3)& - station_char,longitude,latitude,month,day,hour,& - (q1(i),i=1,nquant) + station_char,longitude,latitude,month,day,hour,(q1(i),i=1,nquant) else 2 read(1,*,err=2,end=3)& - station_int,longitude,latitude,month,day,hour,& - (q1(i),i=1,nquant) + station_int,longitude,latitude,month,day,hour,(q1(i),i=1,nquant) endif do i=1,nquant quantiles(i,month,day,hour+1)=q1(i) @@ -217,7 +215,7 @@ open(unit=1,form='formatted',status='old',file=infile) do while(.true.) -11 read(1,*,err=11,end=12)yyear,mmonth,dday,hhour,value1 +11 read(1,*,err=11,end=12) yyear,mmonth,dday,hhour,value1 ! ! Rounding of hours? ! @@ -305,12 +303,10 @@ do while(.true.) if(l_code_in_char)then 21 read(1,*,err=21,end=23)station_char,longitude,latitude,realization,& - (msd_bootstrap(j+1),j=hour1,hour2,hour_step),& - msd_bootstrap(nhour+1) + (msd_bootstrap(j+1),j=hour1,hour2,hour_step),msd_bootstrap(nhour+1) else 22 read(1,*,err=22,end=23)station_int,longitude,latitude,realization,& - (msd_bootstrap(j+1),j=hour1,hour2,hour_step),& - msd_bootstrap(nhour+1) + (msd_bootstrap(j+1),j=hour1,hour2,hour_step),msd_bootstrap(nhour+1) endif ! ! Update the p-value counters diff --git a/SYNOP/STATS/produce_rank_histograms_all_stations.sh b/SYNOP/STATS/produce_rank_histograms_all_stations.sh index b8b064ed220d56bbe06e87465d24a4a50c866ae9..29ea15fff3e0339fb277d1d6a3771fa834d88272 100755 --- a/SYNOP/STATS/produce_rank_histograms_all_stations.sh +++ b/SYNOP/STATS/produce_rank_histograms_all_stations.sh @@ -40,6 +40,7 @@ echo "Bash version ${BASH_VERSION}..." # 1. Variable code # variable=$1 +echo " Meteorological variable code = $1" # # The following variables are included: # @@ -168,6 +169,8 @@ echo " List of synop stations (with geographical coordinates): $station_list" # Compile the Fortran program that calculates the rank histograms # echo " Compiling fortran-program to calculate rank histograms ..." +echo " fortran-programs/rank_histograms_one_station.f95" + gfortran fortran-programs/rank_histograms_one_station.f95 -o rank_histograms_one_station # ################################################################# @@ -186,6 +189,7 @@ while [ ${line_number} -le `expr $n_lines` ] do head -`expr ${line_number}` ${station_list} | tail -1 > input.txt read station longitude latitude < input.txt +echo " **********" echo " Synop station ID, longitude, latitude: ${station} ${longitude} ${latitude}" #echo ${station} @@ -193,17 +197,16 @@ echo " Synop station ID, longitude, latitude: ${station} ${longitude} ${latitu #echo ${latitude} # ################################################################ -# # Select the simulation data for the station (mimicked this far by observations!). # When simulation data are available, a means to retrieve the correct data points from it must be designed +echo " Selecting the simulation data (mimicked his far by observations) for synop station: ${station}" # odb_command="odb sql 'select year,month,day,hour,value where station=${station} and (year>=${year1} and year<=${year2}) and (hour=0 or hour=6 or hour=12 or hour=18) and variable=${variable}' -i ${sim_file} -o sim_data" eval ${odb_command} # ############################################################### -# # Select the quantiles for the station: -# +echo " Selecting the quantiles for synop station: ${station}" # 00, 06, 12 and 18 UTC are picked separately, since # select \* does not allow for parentheses (very strange) # @@ -224,9 +227,9 @@ cat quantile_selection_* > quantile_selection odb sql select \* where station=${station} -i ${msd_bootstrap_file} -o msd_bootstrap_selection # ############################################################### -# # Produce the rank histogram for one station. # Include data for 00, 06, 12 and 18 UTC. +echo " Producing the rank histogram for one synop station: ${station}" # ############################################################## # diff --git a/SYNOP/STATS/produce_standard_plots_all_stations.sh b/SYNOP/STATS/produce_standard_plots_all_stations.sh index 4398a3658afd03107c83092ff709590a25cae3ed..b45052afc7d1cf291b3bf5d5b0758fb915873a19 100755 --- a/SYNOP/STATS/produce_standard_plots_all_stations.sh +++ b/SYNOP/STATS/produce_standard_plots_all_stations.sh @@ -1,5 +1,7 @@ #!/bin/bash echo "Bash version ${BASH_VERSION}..." +echo "Python version" +which python ########################################################################## # # Produce the following plots: @@ -33,8 +35,10 @@ echo "Bash version ${BASH_VERSION}..." # their relevant parts are converted to text format for processing. # # Execution (e.g): ./produce_standard_plot_all_stations 39 2020 2020 1 12 -# +# ORIGINAL: # Jouni Räisänen, August 2023 +# MODIFIED: +# Alexander Mahura, Sep-Oct 2023 # ################################################################## # @@ -43,7 +47,8 @@ echo "Bash version ${BASH_VERSION}..." # 1. Variable code # variable=$1 -# +echo " Meteorological variable code = $1" + # The following variables are included: # # 91 'total amount of clouds' @@ -168,11 +173,13 @@ fi # and is available as a .txt file. # station_list=list_of_stations_${variable}.txt +echo " List of synop stations (with geographical coordinates): $station_list" # # ################################################################ # echo " Compiling the Fortran program needed for creating the figures ..." +echo " fortran-programs/plots_for_one_station.f95" # gfortran fortran-programs/plots_for_one_station.f95 -o plots_for_one_station # @@ -192,24 +199,26 @@ while [ ${line_number} -le `expr $n_lines` ] do head -`expr ${line_number}` ${station_list} | tail -1 > input.txt read station longitude latitude < input.txt +echo " **********" echo " Synop station ID, longitude, latitude: ${station} ${longitude} ${latitude}" # ################################################################ -# # Select the simulation data for the station (mimicked this far by observations!) # When simulation data are available, a means to retrieve the correct data points from it must be designed. +echo " Selecting the simulation data (mimicked his far by observations) for synop station: ${station}" # odb_command="odb sql 'select year,month,day,hour,value where station=${station} and (year>=${year1} and year<=${year2}) and (hour=0 or hour=6 or hour=12 or hour=18) and variable=${variable}' -i ${sim_file} -o sim_data" eval ${odb_command} + ################################################################ -# # Select the rank histogram data for the station +echo " Selecting rank histogram data for synop station: ${station}" # odb sql select \* where station=${station} -i ${rh_file} -o rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} # ############################################################### -# -# Select the quantiles for the station. +# Select the quantiles for the station +echo " Selecting the quantiles for synop station: ${station}" # # 00, 06, 12 and 18 UTC are picked separately, since # select \* does not allow for parentheses (very strange) @@ -245,8 +254,8 @@ cat quantile_selection_* > quantile_selection l_code_in_char=.false.,l_round_6h=.false. &END EOF -########################################################### -# + +################################################################ # Produce the figure in Python # cat < standard_plot.cfg @@ -268,26 +277,24 @@ rank_hists=rank_histogram_${variable}_${station}_${year1}${month1}-${year2}${mon [output] fig_name=${figure_dir}/standard_plot_${variable}_${station}_${year1}${month1}-${year2}${month2}_python.png EOF -######################## + +################################################################ # -#which python -echo " Calling python to plot quantiles rank histogram ..." -#echo Calling Python -#python3 python/plot_quantiles_rankhist.py standard_plot.cfg +echo " Calling python to plot quantiles rank histogram for synop station: ${station}" python3 python/plot_quantiles_rankhist.py standard_plot.cfg -#echo Called Python -####################### -# + +################################################################ # Remove unnecessary files -# +echo " Removing unnecessary temporary files ..." + rm input.txt rm vrange_* rm msd_and_p-value -rm msd_bootstrap_selection rm quantile_selection rm rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} rm sim_data -rm standard_plot_one_station +###rm msd_bootstrap_selection +###rm standard_plot_one_station rm time_series_commands rm time_series_${variable}_${station}_${year1}${month1}-${year2}${month2}.grads rm rank_histogram_${variable}_${station}_${year1}${month1}-${year2}${month2}.grads @@ -298,7 +305,11 @@ rm rank_histogram_${variable}_${station}_${year1}${month1}-${year2}${month2}.txt # ((line_number++)) done -################################################ + +################################################################ # The Fortran executable is not needed any more: -############################################## -rm plots_for_one_station +# +rm plots_for_one_station +rm quantile_selection_* +rm msd_and_p-value_* +rm coordinates diff --git a/SYNOP/STATS/summary_rank_histograms_all_stations.sh b/SYNOP/STATS/summary_rank_histograms_all_stations.sh index 7f038235c3f8327d5844f11ee8e1a0f8c727b24f..ce8e6e465d0e1e8c7f4d35b63addb6ac2048d7b6 100755 --- a/SYNOP/STATS/summary_rank_histograms_all_stations.sh +++ b/SYNOP/STATS/summary_rank_histograms_all_stations.sh @@ -1,5 +1,7 @@ #!/bin/bash echo "Bash version ${BASH_VERSION}..." +echo "Python version" +which python ################################################################## # # Calculation and plotting of summary statistics for quantile @@ -23,8 +25,10 @@ echo "Bash version ${BASH_VERSION}..." # (for the time range specified by the script arguments) # # Execution (e.g): ./summary_rank_histograms_all_stations 39 2020 2020 1 12 -# +# ORIGINAL: # Jouni Räisänen, August 2023 +# MODIFIED: +# Alexander Mahura, Sep-Oct 2023 # ################################################################## # @@ -33,6 +37,7 @@ echo "Bash version ${BASH_VERSION}..." # 1. Meteorological variable code # variable=$1 +echo " Meteorological variable code = $1" # # The following variables are included: # @@ -95,12 +100,11 @@ fi # ################################################################## # -# Compile the Fortran program that produces the rank histogram summary statistics +echo " Compiling the Fortran program that produces the rank histogram summary statistics ..." +echo " fortran-programs/rank_histogram_summary_statistics.f95" # gfortran fortran-programs/rank_histogram_summary_statistics.f95 -o rank_histogram_summary_statistics -echo "COMPILED - program that produces the rank histogram summary statistics" - # ################################################################## # @@ -108,7 +112,9 @@ echo "COMPILED - program that produces the rank histogram summary statistics" # # 1) Rank histogram directory and input name file name without extension # +echo " Directory for rank histogram ..." rh_dir=rank_histograms/${variable}_${year1}${month1}-${year2}${month2} +echo " $rh_dir" rh_file=${rh_dir}/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_all_stations rh_summary_file=${rh_dir}/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_summary # @@ -118,7 +124,10 @@ out_file=${rh_dir}/rh_summary_${variable}_${year1}${month1}-${year2}${month2} # # 3) Directory for figures # +echo " Directory for figures ..." figure_dir=figures +echo " $figure_dir" + # ################################################################## # @@ -129,7 +138,7 @@ odb sql select \* -i ${rh_file}.odb -o ${rh_file}.txt # ################################################################## # -# Calculate the rank histogram summary statistics +echo " Calculating the rank histogram summary statistics ..." # ./rank_histogram_summary_statistics <