diff --git a/SATELLITE/STATST/fortran-programs/plots_for_one_location.f95 b/SATELLITE/STATST/fortran-programs/plots_for_one_location.f95 new file mode 100644 index 0000000000000000000000000000000000000000..4d5e293ce02f70e189e7d749d298e12795d684bd --- /dev/null +++ b/SATELLITE/STATST/fortran-programs/plots_for_one_location.f95 @@ -0,0 +1,654 @@ +program plots_for_one_station + ! + ! For producing files that can be used to generate "one-station standard plots" + ! in python and in GrADS. However, these plots need to be generated separately, + ! i.e., python or GrADS is not called directly from this program. + ! + ! The following plots can be produced, depending on the details of the + ! python or GrADS script: + ! + ! 1) The time series of simulated data at one station location, + ! plotted against the background of the observed quantiles. + ! Separate plots for 00, 06, 12 and 18 UTC. The plots are + ! organized so that all the simulation data are projected + ! on a single annual cycle (normal year, leap days excluded) + ! regardless of the length of the simulation. + ! + ! 2) A quantile rank histogram, showing the normalized frequencies + ! for 100 percentile bins from 0-1% to 99-100%. The Mean + ! Square Deviation and its p-value are written as text in this plot. + ! + ! NB hard-coding: the program currently gives output only for 00, 06, 12 and 18 UTC. + ! + ! ---------------------------------------------------------------------- + ! Jouni Räisänen, University of Helsinki, August 2023 + !----------------------------------------------------------------------- + ! + ! INPUT FILES (all as text files retrieved from .odb format): + ! ------------------------------------------------------------ + ! + ! sim_file = simulated data from the station location of interest + ! + ! Each line of this file includes: + ! + ! 1) year + ! 2) month + ! 3) day + ! 4) hour + ! 5) data value + ! + ! quantile_file = observed quantiles for this station + ! + ! Each line of this file includes: + ! + ! 1) station code + ! 2) longitude + ! 3) latitude + ! 4) month + ! 5) day + ! 6) hour + ! 7-105) quantiles from 1% to 99% + ! + ! rank_histogram_file = frequencies and MSD and P-values for plotting the rank histogram + ! + ! Each line of this file includes: + ! + ! 1) station code + ! 2) longitude + ! 3) latitude + ! 4) UTC hour + ! 5) total number of observation for the UTC hour + ! 6) Mean squared deviation (MSD) of the quantile bin rank histogram wrt a flat histogram + ! 7) p-value of the MSD value, relative to the bootstrap sampling distribution + ! 8-107) frequencies for quantile bins 0-1%, 1-2%, ... ,98-99%, 99-100% + ! + ! ---------------------------------------------------------------------------------------- + ! + ! OUTPUT FILES in ODB compatible text format (for python) + ! + ! txt_data_time_series: individual data values from sim_file + ! + ! The first line is a header. The other lines include: + ! + ! 1) station code + ! 2) longitude + ! 3) latitude + ! 4) year + ! 5) day of year (counted from Jan 1) + ! 6) UTC hour + ! 7) data value + ! + ! txt_data_quantiles: selected quantiles from the oberved distribution as a function of UTC hour and time of year + ! + ! The first line is a header. The other lines include: + ! + ! 1) station code + ! 2) longitude + ! 3) latitude + ! 4) day of year (counted from Jan 1) + ! 5) UTC hour + ! 6-12) quantiles 1, 10, 25, 50, 75, 90 and 99% + ! + ! txt_data_rank_histogram = quantile bin frequencies (separately for each UTC hour) + ! + ! The first line is a header. The other lines include: + ! + ! 1) station code + ! 2) longitude + ! 3) latitude + ! 4) UTC hour + ! 5) Mean squared deviation (MSD) of the quantile bin rank histogram wrt a flat histogram + ! 6) p-value of the MSD value, relative to the bootstrap sampling distribution + ! 7-106) frequencies for quantile bins 0-1%, 1-2%, ... ,98-99%, 99-100% + ! + ! ---------------------------------------------------------------------------------------- + ! + ! OUTPUT FILES in GrADS binary format (for GrADS) + ! + ! grads_data_time_series = individual data values from sim_file + ! + values of selected quantiles as a fuction + ! of UTC hour and time of the year. + ! + ! grads_data_rank_histogram = quantile bin frequencies (separately for each UTC hour) + ! + ! + ! In addition, the program writes a few auxilliary texts files with hard-coded names + ! for use in GrADS scripts: + ! + ! "vrange_commands" : ranges needed on the y axis of plots + ! "time_series_commands" : for plotting different years in time series with different colors + ! "coordinates" : for including coordinate values in GrADS plot + ! "msd_and_p-value_00" : for including 00 UTC MSD and p-value in GrADS plot + ! "msd_and_p-value_06" : for including 06 UTC MSD and p-value in GrADS plot + ! "msd_and_p-value_12" : for including 12 UTC MSD and p-value in GrADS plot + ! "msd_and_p-value_18" : for including 18 UTC MSD and p-value in GrADS plot + ! "msd_and_p-value" : for including all-UTC MSD and p-value in GrADS plot + ! + !------------------------------------------------------------------------------------------ + ! + ! NAMELIST PARAMETERS: + ! + ! a) input files (see above) + ! + ! sim_file + ! quantile_file + ! rank_histogram_file + ! + ! b) output files in ODB compatible text format (see above) + ! + ! txt_data_time_series + ! txt_data_rank_histogram + ! txt_data_quantiles + ! + ! c) output files in GrADS binary format + ! + ! grads_data_time_series + ! grads_data_rank_histogram + ! + ! d) other + ! + ! year1,month1 = beginning of the period in time series output + ! year2,month2 = end of the period in time series output + ! include_rank_histograms: if .true. (.false.), rank histogram data is included (excluded) in output + ! l_code_in_char : if .true., (.false.) station codes are assumed to be 3-character strings (integers) + ! l_round_6h : if .true., UTC hours are rounded to the closesest even 6-h (00, 06, 12 or 18 UTC). + ! if .false., data for non-matching UTC hours are excluded from output + ! + ! miss = missing value code. Default = -9.99e6. All data f with |f| >= |miss| are treated as missing + ! + !------------------------------------------------------------------------------------------- + ! + implicit none + character*160 :: sim_file ! see above + character*160 :: quantile_file ! see above + character*160 :: rank_histogram_file ! see above + character*160 :: grads_data_time_series ! see above + character*160 :: grads_data_rank_histogram ! see above + character*160 :: txt_data_time_series ! see above + character*160 :: txt_data_quantiles ! see above + character*160 :: txt_data_rank_histogram ! see above + integer :: year1,year2 ! first and last year of simulation + integer :: month1,month2 ! first and last month of simulation + integer :: nyear ! number of years = year2-year1+1 + ! + ! Data structures for input data + ! + integer,parameter :: nhour=1 ! only 12 UTC used used for AMSU-A!!! + integer,parameter :: nmax_years=200 ! maximum length of a simulation in years + integer,parameter :: ndays_year=365 ! number of days in a non-leap year + + real :: f(nmax_years,ndays_year,nhour) ! the simulated values + real :: f1 ! a single simulated value + ! + integer,parameter :: nquant=99 ! number of quantiles + real :: quant(nquant,ndays_year,nhour) ! quantile values from quantile_file + real :: quant1(nquant) ! quantile values for a single day and UTC hour + ! + real :: freq(nquant+1,nhour+1) ! quantile bin frequencies from the rank histogram file + real :: freq1(nquant+1) ! frequencies for a single UTC hour + real :: msd(nhour+1) ! MSD value from the rank histogram file + real :: msd1 ! a single MSD value + real :: p_value(nhour+1) ! p-values from the rank histogram files + real :: p_value1 ! a single p-value + + integer :: yyear,mmonth,dday,hhour ! time as read from sim_file + integer :: year,month,day,hour ! time after eventual rounding to nearest 6-h + character*3 :: station_char ! station code as 3-char string + integer :: station_int ! station code as integer + logical :: l_code_in_char ! .true. for station codes as string + logical :: l_round_6h ! .true. for rounding time to nearest 6-h + real :: ntot ! total number of observations as read from rank_histogram_file + real :: lon,lat ! longitude and latitude + integer :: i,j,k ! loop variables + integer :: day_of_year ! day count from Jan 1 (in a non-leap year) + + real :: miss ! missing value code + + real :: fmin(nhour),fmax(nhour) ! minimum and maximum values for time series plots + real :: fmin_down(nhour),fmax_up(nhour) ! same, but extended outward + real :: freqmax(nhour+1) ! largest frequency for the rank histogram + + 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 + + integer :: nlines, io +! +! For text output in odb-compatible format: +! + INTEGER, PARAMETER :: L_dataline=1700 ! max length of input or output text line + character*1700 :: dataline,headerline,emptyline ! for reading / writing text files + character*2 :: number_of_bin ! quantile bin numbers for header in txt_data_quantiles + character*16 :: frequency_value ! frequency value written as character string +! + logical :: include_rank_histograms ! true for including rank histograms + + namelist/param/sim_file,quantile_file,rank_histogram_file,& + grads_data_time_series,grads_data_rank_histogram,& + year1,year2,month1,month2,& + txt_data_time_series,txt_data_rank_histogram,txt_data_quantiles,& + include_rank_histograms,l_code_in_char,l_round_6h,miss + + data color/9,14,4,11,5,13,3,10,7,12,8,2,6/ +! +! default values +! + miss=-9.99e6 + l_code_in_char=.false. + l_round_6h=.false. + include_rank_histograms=.true. +! +! Read the namelist and count the number of years for output +! + read(*,nml=param) + nyear=year2-year1+1 + + !************************************************************************************************ + ! Read the data from the input files. Only data for the period year1/month1-year2/month2 + ! will be included in the time series output, even if other data exist. + ! + ! 1) time series + ! + + f=miss + print*,' Reading sim_file: ',sim_file + open(unit=1,form='formatted',file=sim_file,status='old') + do while(.true.) + !do i=1,nlines +11 read(1,*,err=11,end=12) yyear,mmonth,dday,hhour,f1 +!11 read(1,*,err=11,end=12) station_char,yyear,mmonth,dday,f1,hhour +! read(1,*) station_char,yyear,mmonth,dday,f1,hhour + !print*,yyear,mmonth,dday,f1,hhour + if(l_round_6h)then + call round_6h(yyear,mmonth,dday,hhour,year,month,day,hour) ! Rounding to nearest 6-h? + else + year=yyear + month=mmonth + day=dday + hour=12 !hhour + endif + !print*,i,yyear,mmonth,dday,f1,hhour +! print*,'hep1' + !print*,1+hour/6 + if((year.gt.year1.or.(year.eq.year1.and.month.ge.month1)).and. & ! Is the time within the + (year.lt.year2.or.(year.eq.year2.and.month.le.month2))) then ! wanted range? + !print*,'hep' + if(.not.(month.eq.2.and.day.eq.29))then ! leap days are skipped from plotting + f(year-year1+1,day_of_year(month,day),1)=f1 ! hard-coding to 6-hourly resolution + endif + endif + !print*,'hep2',f1, f(year-year1+1,day_of_year(month,day),1) + !print*,i,yyear,mmonth,dday,f1,hhour + enddo +12 continue + close(1) + ! + ! 2) quantiles + ! + quant=miss + print*,' Reading quantile file: ',quantile_file + open(unit=1,form='formatted',file=quantile_file,status='old') + do while(.true.) + if(l_code_in_char)then +21 read(1,*,err=21,end=23)station_char,lon,lat,month,day,hour,(quant1(i),i=1,nquant) + !print*,station_char,lon,lat,month,day,hour + else +22 read(1,*,err=22,end=23)station_int,lon,lat,month,day,hour,(quant1(i),i=1,nquant) + endif + if((hour.eq.0.or.hour.eq.6.or.hour.eq.12.or.hour.eq.18)& + .and.(.not.(month.eq.2.and.day.eq.29))) then ! leap days are skipped from plotting + do i=1,nquant + quant(i,day_of_year(month,day),1)=quant1(i) + enddo + endif + enddo +23 continue + close(1) + ! + ! 3) quantile bin frequencies and related MSD and p-value statistics + ! +! open(unit=2,form='formatted',status='old',file=rank_histogram_file) +! nlines = 0 +! DO +! READ(2,*,iostat=io) +! IF (io/=0) EXIT +! nlines = nlines + 1 +! ENDDO +! close(2) + + freq=miss + if(include_rank_histograms)then + print*,' Reading rank histograms: ',rank_histogram_file + open(unit=1,form='formatted',file=rank_histogram_file,status='old') + do while(.true.) + !do j=2,nlines + if(l_code_in_char)then + !print*,'hep0' +31 read(1,*,err=31,end=33)station_char,lon,lat,hour,ntot,msd1,p_value1,(freq1(i),i=1,nquant+1) + else +32 read(1,*,err=32,end=33)station_int,lon,lat,hour,ntot,msd1,p_value1,(freq1(i),i=1,nquant+1) + read(1,*)station_int,lon,lat,hour,ntot,msd1,p_value1,(freq1(i),i=1,nquant+1) + endif + !if((hour.eq.0.or.hour.eq.6.or.hour.eq.12.or.hour.eq.18.or.hour.eq.24))then + !print*,'hep1' + if((hour.eq.12))then + !print*,'hep1',hour,msd1,p_value1 + msd(1)=msd1 + p_value(1)=p_value1 + do i=1,nquant+1 + freq(i,1)=freq1(i) + enddo + !print*,'hep2',msd(1),p_value(1) + endif + !print*,'hep2' + enddo +33 continue + close(1) + endif + +!********************************************************************************************** +! Find the minimum and maximum values of f (or the 1st and 99th percentile) +! for defining the vranges for the GrADS time series plots + + do j=1,nhour + fmin(j)=quant(1,1,j) + fmax(j)=quant(nquant,1,j) + do i=1,ndays_year + if(quant(1,i,j).lt.fmin(j))fmin(j)=quant(1,i,j) + if(quant(nquant,i,j).gt.fmax(j))fmax(j)=quant(nquant,i,j) + do k=1,nyear + if(f(k,i,j).lt.fmin(j).and.abs(f(k,i,j)).lt.abs(miss))fmin(j)=f(k,i,j) + if(f(k,i,j).gt.fmax(j).and.f(k,i,j).lt.abs(miss))fmax(j)=f(k,i,j) + enddo + enddo + fmin_down(j)=fmin(j)-0.05*(fmax(j)-fmin(j)) + fmax_up(j)=fmax(j)+0.05*(fmax(j)-fmin(j)) + enddo +! +! Find the largest frequency from the rank histograms +! + freqmax=0. + do j=1,nhour+1 + do i=1,nquant+1 + if(freq(i,j).gt.freqmax(j))freqmax(j)=freq(i,j) + enddo + enddo +! +! Write the commands to the auxiliary script files for Grads +! + open(unit=1,form='formatted',file='vrange_commands') + write(1,'(A11,2F10.3)')'set vrange ',fmin_down(1),fmax_up(1) + write(1,'(A11,2F10.3)')'set vrange ',fmin_down(2),fmax_up(2) + write(1,'(A11,2F10.3)')'set vrange ',fmin_down(3),fmax_up(3) + write(1,'(A11,2F10.3)')'set vrange ',fmin_down(4),fmax_up(4) + write(1,'(A13,F7.3)')'set vrange 0 ',1.1*(nquant+1)*freqmax(1) ! conversion to relative frequencies + write(1,'(A13,F7.3)')'set vrange 0 ',1.1*(nquant+1)*freqmax(2) ! conversion to relative frequencies + write(1,'(A13,F7.3)')'set vrange 0 ',1.1*(nquant+1)*freqmax(3) ! conversion to relative frequencies + write(1,'(A13,F7.3)')'set vrange 0 ',1.1*(nquant+1)*freqmax(4) ! conversion to relative frequencies + write(1,'(A13,F7.3)')'set vrange 0 ',1.1*(nquant+1)*freqmax(5) ! conversion to relative frequencies + close(1) + + open(unit=1,form='formatted',file='time_series_commands') + do i=1,nyear + write(1,'(A12)')'set cthick 4' + write(1,'(A12)')'set cstyle 0' + write(1,'(A11)')'set cmark 3' + write(1,'(A15)')'set digsiz 0.07' + if(nyear.eq.1)then + 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 + 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 + close(1) + + open(unit=1,form='formatted',file='coordinates') + write(1,'(F8.3)')lon + write(1,'(F8.3)')lat + close(1) + + if(include_rank_histograms)then +! open(unit=1,form='formatted',file='msd_and_p-value_00') +! write(1,'(A15)')'set strsiz 0.17' +! write(1,'(A18)')'set string 1 l 4 0' +! write(1,'(A25,F5.3)')'draw string 1.3 &0 MSD=',msd(1) +! write(1,'(A29,F5.3)')'draw string 3.0 &0 p-value=',p_value(1) +! close(1) +! open(unit=1,form='formatted',file='msd_and_p-value_06') +! write(1,'(A15)')'set strsiz 0.17' +! write(1,'(A18)')'set string 1 l 4 0' +! write(1,'(A25,F5.3)')'draw string 1.3 &0 MSD=',msd(2) +! write(1,'(A29,F5.3)')'draw string 3.0 &0 p-value=',p_value(2) +! close(1) +! open(unit=1,form='formatted',file='msd_and_p-value_12') +! write(1,'(A15)')'set strsiz 0.17' +! write(1,'(A18)')'set string 1 l 4 0' +! write(1,'(A25,F5.3)')'draw string 1.3 &0 MSD=',msd(3) +! write(1,'(A29,F5.3)')'draw string 3.0 &0 p-value=',p_value(3) +! close(1) +! open(unit=1,form='formatted',file='msd_and_p-value_18') +! write(1,'(A15)')'set strsiz 0.17' +! write(1,'(A18)')'set string 1 l 4 0' +! write(1,'(A25,F5.3)')'draw string 1.3 &0 MSD=',msd(4) +! write(1,'(A29,F5.3)')'draw string 3.0 &0 p-value=',p_value(4) +! close(1) + open(unit=1,form='formatted',file='msd_and_p-value') + write(1,'(A15)')'set strsiz 0.17' + write(1,'(A18)')'set string 1 l 4 0' + write(1,'(A25,F5.3)')'draw string 1.3 &0 MSD=',msd(1) !msd(5) + write(1,'(A29,F5.3)')'draw string 3.0 &0 p-value=',p_value(1) !p_value(5) + close(1) + endif +! +!*************************************************************** +! Write the GrADS binary files needed for the plots + + open(unit=1,form='unformatted',file=grads_data_time_series,access='DIRECT',& + recl=4,status='unknown') + do j=1,nhour + do i=1,ndays_year + irec=((i-1)+(j-1)*ndays_year)*(nyear+7) + do k=1,nyear + write(1,rec=irec+k)f(k,i,j) ! time series + enddo + write(1,rec=irec+nyear+1)quant(1,i,j) ! selected quantile + write(1,rec=irec+nyear+2)quant(10,i,j) + write(1,rec=irec+nyear+3)quant(25,i,j) + write(1,rec=irec+nyear+4)quant(50,i,j) + write(1,rec=irec+nyear+5)quant(75,i,j) + write(1,rec=irec+nyear+6)quant(90,i,j) + write(1,rec=irec+nyear+7)quant(99,i,j) + enddo + enddo + close(1) + + if(include_rank_histograms)then + open(unit=1,form='unformatted',file=grads_data_rank_histogram,access='DIRECT',& + recl=4,status='unknown') + do j=1,nhour+1 + irec=(nquant+1)*(j-1) + do i=1,nquant+1 + write(1,rec=irec+i)freq(i,j) !quantile bin frequencies + enddo + enddo + close(1) + endif + +!*************************************************************** +! Write the time series, quantiles and rank histogram data to text files + +! +! Empty data line +! + do i=1,L_dataline + emptyline(i:i)=' ' + enddo +! +! Time series. 365 days per year, 4 times per day, from year1 to year2 +! + open(unit=1,form='formatted',file=txt_data_time_series) + headerline=emptyline + if(l_code_in_char)then + headerline='station@hdr:string longitude@hdr:real latitude@hdr:real year@hdr:integer'//& + ' day_of_year@hdr:integer hour@hdr:integer value@body:real' + else + headerline='station@hdr:integer longitude@hdr:real latitude@hdr:real year@hdr:integer'//& + ' day_of_year@hdr:integer hour@hdr:integer value@body:real' + endif + write(1,*)trim(headerline) + + do year=year1,year2 + do j=1,365 + !do hour=0,18,6 + do hour=12,12 + !print*,'###########################################' + !print*,year-year1+1,j,1+hour/6 !f(1,1,1) + !print*,f(year-year1+1,j,1+hour/6),miss + !print*,f(:,j,1) + if(abs(f(year-year1+1,j,1)).lt.abs(miss))then + if(l_code_in_char)then + write(1,'(A3,2F16.6,1X,I4,1X,I3,1X,I2,F16.6)')station_char,lon,lat,year,j,hour,f(year-year1+1,j,1) + else + write(1,'(I6,2F16.6,1X,I4,1X,I3,1X,I2,F16.6)')station_int,lon,lat,year,j,hour,f(year-year1+1,j,1) + endif + else + if(l_code_in_char)then + write(1,'(A3,2F16.6,1X,I4,1X,I3,1X,I2,A16)')station_char,lon,lat,year,j,hour,' NaN' + else + write(1,'(I6,2F16.6,1X,I4,1X,I3,1X,I2,A16)')station_int,lon,lat,year,j,hour,' NaN' + endif + endif + enddo + enddo + enddo + close(1) +! +! Selected quantiles of the observed distribution. 365 days per year, 4 time per day. +! + open(unit=1,form='formatted',file=txt_data_quantiles) + headerline=emptyline + if(l_code_in_char)then + headerline='station@hdr:string longitude@hdr:real latitude@hdr:real day_of_year@hdr:integer hour@hdr:integer'//& + ' q01@body:real q10@body:real q25@body:real q50@body:real q75@body:real q90@body:real q99@body:real' + else + headerline='station@hdr:integer longitude@hdr:real latitude@hdr:real day_of_year@hdr:integer hour@hdr:integer'//& + ' q01@body:real q10@body:real q25@body:real q50@body:real q75@body:real q90@body:real q99@body:real' + endif + write(1,*)trim(headerline) + do year=year1,year2 + do j=1,365 + !do hour=0,18,6 + do hour=12,12 + if(l_code_in_char)then + write(1,'(A3,2F16.6,1X,I3,1X,I2,7F16.6)')station_char,lon,lat,j,hour,& + quant(1,j,1),quant(10,j,1),quant(25,j,1),quant(50,j,1),& + quant(75,j,1),quant(90,j,1),quant(99,j,1) + else + write(1,'(I6,2F16.6,1X,I3,1X,I2,7F16.6)')station_int,lon,lat,j,hour,& + quant(1,j,1),quant(10,j,1),quant(25,j,1),quant(50,j,1),& + quant(75,j,1),quant(90,j,1),quant(99,j,1) + endif + enddo + enddo + enddo + close(1) +! +! Rank histogram of frequencies. 100 bins from 0-1% to 99-100. Before that, the MSD and p-values + ! + if(include_rank_histograms)then + open(unit=1,form='formatted',file=txt_data_rank_histogram) + headerline=emptyline + if(l_code_in_char)then + headerline='station@hdr:string longitude@hdr:real latitude@hdr:real hour@hdr:integer' + else + headerline='station@hdr:integer longitude@hdr:real latitude@hdr:real hour@hdr:integer' + endif + headerline=trim(headerline)//' msd@body:real p_value@body:real' + do j=0,nquant + write(number_of_bin,'(i2.2)')j + headerline=trim(headerline)//' f'//number_of_bin//'@body:real' + enddo + write(1,*)trim(headerline) + + !do hour=0,24,6 + do hour=12,12 + if(l_code_in_char)then + write(dataline,'(A3,2F16.6,1X,I2,2F16.6)')& + station_char,lon,lat,hour,msd(1),p_value(1) + else + write(dataline,'(I6,2F16.6,1X,I2,2F16.6)')& + station_int,lon,lat,hour,msd(1),p_value(1) + endif + + do i=1,nquant+1 + write(frequency_value,'(F16.6)')freq(i,1) + dataline=trim(dataline)//frequency_value + enddo + write(1,*)trim(dataline) + enddo + + close(1) + endif + + end program plots_for_one_station + + integer function day_of_year(month,day) +! +! Number of day 'day' in the 'month':th month from the beginning of a non-leap year +! + implicit none + integer :: year,month,day + integer,parameter :: nmon=12 + integer :: days_before(nmon) ! this is for normal years + integer :: dnumber + data days_before / 0,31,59,90,120,151,181,212,243,273,304,334 / + + day_of_year=days_before(month)+day + return + end function day_of_year + +! +! rounding the time to the nearest full 6 hours +! +subroutine round_6h(year,month,day,hour,year1,month1,day1,hour1) +implicit none +integer :: year,month,day,hour ! time before rounding +integer :: year1,month1,day1,hour1 ! time after rounding +integer :: ndays_month ! number of days in a month +ndays_month=31 +if(month.eq.4.or.month.eq.6.or.month.eq.9.or.month.eq.11)ndays_month=30 +if(month.eq.2)then + if((mod(year,4).eq.0.and.mod(year,100).gt.0).or.(mod(year,400).eq.0))then + ndays_month=29 + else + ndays_month=28 + endif +endif +year1=year +month1=month +day1=day +hour1=6*((hour+3)/6) +if(hour.ge.21)then ! Hour was rounded forward to 00 UTC -> increment day + hour1=0 + day1=day+1 + if(day1.gt.ndays_month)then ! increment month + day1=1 + month1=month1+1 + if(month1.gt.12)then ! increment year + month1=1 + year1=year+1 + endif + endif +endif +return +end subroutine round_6h + + + diff --git a/SATELLITE/STATST/fortran-programs/rank_histogram_sumstats_locations.f95 b/SATELLITE/STATST/fortran-programs/rank_histogram_sumstats_locations.f95 new file mode 100644 index 0000000000000000000000000000000000000000..e4fe3456acaa2ab60e5929bf7bb625e2d567849a --- /dev/null +++ b/SATELLITE/STATST/fortran-programs/rank_histogram_sumstats_locations.f95 @@ -0,0 +1,381 @@ + 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: ', 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 + !print*,station1_char,longitude1,latitude1,hour,total_number,msd,p_value1 +! 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: # of lines, # of stations: ',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/SATELLITE/STATST/fortran-programs/rank_histograms_one_location.f95 b/SATELLITE/STATST/fortran-programs/rank_histograms_one_location.f95 new file mode 100644 index 0000000000000000000000000000000000000000..acee2b1aae02a89a28cbe2b2b51cb9f216004d45 --- /dev/null +++ b/SATELLITE/STATST/fortran-programs/rank_histograms_one_location.f95 @@ -0,0 +1,618 @@ + program rank_histograms_one_station +! +! Produce the rank histogram (i.e., the frequencies at which the simulated or +! observed data at one station fall between different quantiles of a pre-computed +! quantile distribution (below 1%, within 1-2%, ... above 99%) +! +! Jouni Räisänen, July 2015 +! +!---------------------------------------------------------------------- +! Input files (text format) +!---------------------------------------------------------------------- +! +! INFILE: simulated or observed time series at the station location. +! The first line is a header. The other lines include: +! +! 1) year +! 2) month +! 3) day +! 4) hour +! 5) data value +! +! QUANTILE_FILE: quantile values, separately for each day of the year! +! The first line is a header. The other lines include: +! +! 1) station code +! 2) longitude +! 3) latitude +! 4) month +! 5) day +! 6) hour +! 7-105) The 99 quantiles from 1 % to 99 % +! +! if (l_code_in_char) is .true., the station code is assumed to be a 3-char string. +! Otherwise, it is assumed to be an integer +! +! MSD_BOOTSTRAP_FILE: MSD (mean squared deviation) distribution of quantile frequencies +! from bootstrap tests. +! +! The first 1st line is a header. The other lines include +! +! 1) station code +! 2) longitude +! 3) latitude +! 4) number of realization +! 5 to N-1) MSD values for HOUR1, HOUR1 + HOUR_STEP ... HOUR2 +! N) MSD calculated from frequencies averaged over all the UTC hours within +! HOUR1, HOUR1 + HOUR_STEP ... HOUR2 +! +! if (l_code_in_char) is .true., the station code is assumed to be a 3-char string. +! Otherwise, it is assumed to be an integer +! +!---------------------------------------------------------------------- +! Output files: +!---------------------------------------------------------------------- +! +! OUTFILE: frequencies of quantiles. +! The first 1st line is a header. The other lines include +! +! 1) station code +! 2) longitude +! 3) latitude +! 4) hour +! 5) number of obervations for this (UTC) hour +! 6) MSD for the quantile frequencies +! 7) p-value of MSD (=fraction of bootstrap realizations with larger or equal MSD) +! 8-107) quantile frequencies (on scale 0-1) from 0-1 % ... 99-100% +! +! if (l_code_in_char) is .true., the station code is assumed to be a 3-char string. +! Otherwise, it is assumed to be an integer. +! +! UTC hours included: HOUR1, HOUR1 + HOUR_STEP ... HOUR2 + all-UTC-statistics coded with HOUR=24 +! +! OUTFILE_GRADS: frequencies of quantiles in GrADS binary output (if grads_out=.true.) +! +! The resulting GraDS binary file includes 7 variables (1 level for 1-6, 100 levels for 7): +! +! 1) longitude +! 2) latitude +! 3) hour +! 4) number of obervations for this (UTC) hour +! 5) MSD for the quantile frequencies +! 6) p-value of MSD (=fraction of bootstrap realizations with larger or equal MSD) +! 7) quantile frequencies (on scale 0-1) from 0-1 % ... 99-100% +! +!--------------------------------------------------------------- +! Namelist parameters +!--------------------------------------------------------------- +! +! NAMELIST/param/infile,outfile,quantile_file,msd_bootstrap_file,& +! hour1,hour2,hour_step,year1,year2,month1,month2,grads_out,outfile_grads,& +! l_code_in_char,l_round_hour,miss +! +! infile : file name for simulated or observed time series (see above) +! quantile_file : file name for pre-computed quantiles (see above) +! msd_bootstrap_file : file name for pre-computed bootstrap statistics (see above) +! outfile : name of output file (see above) +! outfile_grads : output file in GrADS binary format (see above) +! grads_out : if .true., outfile_grads is written +! l_code_in_char: if .true., 3-character station codes assumed (else integer) +! +! hour1 = first UTC hour analysed +! hour2 = last UTC hour analysed +! hour_step = step of analyzed UTC hours +! l_round_hour : if .true., hours are rounded to the nearest analysis hour. +! if .false., data for non-matching hours is ignored +! +! year1,month1 = beginning of the analysis period +! year2,month2 = end of the analysis period +! +! miss = missing value code. All input data f with abs(f) >= abs(miss) is treated as missing +! +!--------------------------------------------------------------- +! + IMPLICIT NONE + INTEGER :: I,J,K,Q + INTEGER,PARAMETER :: NMON=12, NDAY=31, NHOUR=24, NQUANT=99 + REAL :: quantiles(nquant,nmon,nday,nhour) ! array for quantile values + REAL :: q1(nquant) ! quantiles for one day + hour + REAL :: frequency(nquant+1,nhour+1) ! frequency of quantiles. NHOUR+1 for the all-UTC frequencies + REAL :: expected_freq ! expected frequency of quantiles (1/100) + REAL :: MSD(nhour+1) ! MSD statistics per UTC hour + REAL :: p_value(nhour+1) ! p values per UTC hour + REAL :: msd_bootstrap(nhour+1) ! msd values for one bootstrap realization + INTEGER :: n_bootstrap ! total number of bootstrap realizations + INTEGER :: realization ! number of bootstrap realization read from msd_boostrap_file + REAL :: value1 ! individual data values read from the input file + INTEGER :: ntimes_hour(nhour) ! total number of times with data for an UTC hour + INTEGER :: ntimes_hour_tot ! total number of times with data for all UTC hours + CHARACTER*160 :: infile,outfile,quantile_file,msd_bootstrap_file + CHARACTER*160 :: outfile_grads + INTEGER,PARAMETER :: L_dataline=1700 ! maximum length of lines in text files + character*1700 :: dataline,headerline,emptyline ! character strings for writing of output files + character*2 :: number_of_bin ! two-digit code for quantiles for outfile header + character*16 :: frequency_value ! frequency written in F16.6 format +! + INTEGER :: yyear,mmonth,dday,hhour ! time as read from the input file + INTEGER :: year,month,day,hour ! time after eventual rounding of hours + + INTEGER :: year1,year2,month1,month2 ! period used in calculations (see above) + INTEGER :: hour1,hour2,hour_step ! hours used in calculations (see above) + LOGICAL :: HOUR_USED(nhour) ! true if the UTC hour is used + INTEGER :: nhour_used ! number of UTC hours used + INTEGER :: station_int ! station code integer (for synop stations amd RHARM soundings) + CHARACTER*3 :: station_char ! station with characters (for GRUAN soundings) + LOGICAL :: l_code_in_char ! .true. for station code in characters + REAL :: longitude,latitude ! station longitude and latitude + + LOGICAL :: GRADS_OUT ! Output also as GrADS binaries, just for testing + INTEGER :: IREC ! record number for GrADS output + + LOGICAL :: L_ROUND_HOUR ! if .true., hour rounded to the nearest output hour +! + REAL :: MISS ! missing value code + integer :: nlines, io +! +!----------------------------------------------------------------------- +! + NAMELIST/param/infile,outfile,quantile_file,msd_bootstrap_file,& + hour1,hour2,hour_step,year1,year2,month1,month2,grads_out,outfile_grads,& + l_code_in_char,l_round_hour,miss + + MISS=-9.99e6 + READ(*,NML=PARAM) +! +! 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 +!---------------------------------------------------------------------- +! +! Read the contents of the quantile_file to array 'quantiles' +! + quantiles=miss + +! open(unit=2,form='formatted',status='old',file=quantile_file) + print*,'Reading quantilefile: ',quantile_file + +! nlines = 0 +! DO +! READ(2,*,iostat=io) +! IF (io/=0) EXIT +! nlines = nlines + 1 +! END DO +! close(2) +! print*,nlines + + open(unit=1,form='formatted',status='old',file=quantile_file) + + do while(.true.) + !do j=2,nlines !while(.true.) + !print*,j + if(l_code_in_char)then + !print*,j + ! 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(k),k=1,nquant) + !read(1,*) station_char,longitude,latitude,month,day,hour,(q1(k),k=1,nquant) + !read(1,*) station_char,longitude + !print*,q1(10) + else +2 read(1,*,err=2,end=3) station_int,longitude,latitude,month,day,hour,(q1(i),i=1,nquant) + !read(1,*,err=2,end=3)& + ! 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) + enddo +! copy quantiles from 28 February to 29 February + if(month.eq.2.and.day.eq.28)then + do i=1,nquant + quantiles(i,month,day+1,hour+1)=q1(i) + enddo + endif + enddo +3 continue + close(1) +! +!---------------------------------------------------------------------- +! +! Read the station observations and count the absolute quantile bin frequencies +! + frequency=0 + ntimes_hour=0 + ntimes_hour_tot=0 + +! open(unit=2,form='formatted',status='old',file=infile) + print*,'Reading infile: ',infile + +! nlines = 0 +! DO +! READ(2,*,iostat=io) +! IF (io/=0) EXIT +! nlines = nlines + 1 +! END DO +! close(2) + + open(unit=1,form='formatted',status='old',file=infile) + + do while(.true.) +! do i=2,nlines !while(.true.) +!11 read(1,*,err=11,end=12)station_char,yyear,mmonth,dday,value1,hhour +11 read(1,*,err=11,end=12)yyear,mmonth,dday,hhour,value1 +! read(1,*)station_char,yyear,mmonth,dday,value1,hhour + !print*,yyear,mmonth,dday,hhour,value1 +! +! Rounding of hours? +! + if(l_round_hour)then + call round_hour(yyear,mmonth,dday,hhour,year,month,day,hour,& + hour1,hour2,hour_step) + else + year=yyear + month=mmonth + day=dday + hour=12 !hhour + endif +! +! If the hour is used and the date is within th analyzed period, +! update the quantile bin frequencies + + !print*,hour_used + hour=12 + !print*,hour+1 + if(hour_used(hour+1))then + !print*,'#########################' + if((year.gt.year1.or.(year.eq.year1.and.month.ge.month1)).and. & + (year.lt.year2.or.(year.eq.year2.and.month.le.month2)).and. & + abs(value1).lt.abs(miss))then + ntimes_hour(hour+1)=ntimes_hour(hour+1)+1 + ntimes_hour_tot=ntimes_hour_tot+1 + call update_frequencies & + (frequency(1,hour+1),value1,quantiles(1,month,day,hour+1),nquant,miss) + endif + endif + + enddo + station_char='a00' +12 continue + close(1) +! +!---------------------------------------------------------------------- +! +! Convert absolute frequencies to relative frequencies +! + !print*,frequency(:,:) + do hour=hour1,hour2,hour_step + do i=1,nquant+1 + if(ntimes_hour(hour+1).gt.0)then + !print*,'hep' + frequency(i,hour+1)=frequency(i,hour+1)/ntimes_hour(hour+1) + else + frequency(i,hour+1)=miss + endif + enddo + enddo +! +!-------------------------------------------------------------------------------- +! +! Calculation of MSD. Because it can be assumed that model data are always available, +! all UTC hours get the same weight in the calculation of the all-UTC MSD. +! + expected_freq=1./(nquant+1.) + !print*,frequency ! -->> missing + do j=hour1,hour2,hour_step + if(frequency(1,j+1).eq.miss)then + msd(j+1)=miss + else + do i=1,nquant+1 + msd(j+1)=msd(j+1)+((frequency(i,j+1)/expected_freq-1.)**2.)/(nquant+1.) + enddo + endif + enddo + do i=1,nquant+1 + do j=hour1,hour2,hour_step + frequency(i,nhour+1)=& + frequency(i,nhour+1)+frequency(i,j+1)/nhour_used + enddo + msd(nhour+1)=msd(nhour+1)+((frequency(i,nhour+1)/expected_freq-1.)**2.)/(nquant+1.) + enddo + do j=hour1,hour2,hour_step + if(frequency(1,j+1).eq.miss)then + do i=1,nquant+1 + frequency(i,nhour+1)=miss + enddo + msd(nhour+1)=miss + endif + enddo + !print*,msd + +!---------------------------------------------------------------------- +! +! Read the bootstrap MSD:s and calculate the fractions of them that exceed the actual value +! + n_bootstrap=0 + p_value=0. + +! open(unit=2,form='formatted',status='old',file=msd_bootstrap_file) + print*,'Reading msd_bootstrap_file: ',msd_bootstrap_file + +! nlines = 0 +! DO +! READ(2,*,iostat=io) +! IF (io/=0) EXIT +! nlines = nlines + 1 +! END DO +! close(2) + + open(unit=1,form='formatted',status='old',file=msd_bootstrap_file) + do while(.true.) + !do i=2,nlines ! 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) + !read(1,*) station_char,longitude,latitude,realization,& + ! (msd_bootstrap(j+1),j=hour1,hour2,hour_step),& + ! msd_bootstrap(nhour+1) + !print*,station_char,longitude,latitude,realization !OK + !print*,(msd_bootstrap(j+1),j=hour1,hour2,hour_step),msd_bootstrap(nhour+1) !OK + 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) + !read(1,*) station_int,longitude,latitude,realization,& + + endif +! +! Update the p-value counters +! + n_bootstrap=n_bootstrap+1 + do hour=hour1,hour2,hour_step + if(msd_bootstrap(hour+1).ge.msd(hour+1))then + p_value(hour+1)=p_value(hour+1)+1. + endif + enddo + if(msd_bootstrap(nhour+1).ge.msd(nhour+1))then + p_value(nhour+1)=p_value(nhour+1)+1. + endif + enddo +23 continue +! +! Convert the p_values from absolute counts to relative frequencies +! + do hour=hour1,hour2,hour_step + p_value(hour+1)=p_value(hour+1)/n_bootstrap + enddo + p_value(nhour+1)=p_value(nhour+1)/n_bootstrap +! +!-------------------------------------------------------------------------- +! +! Open the output file and write its header line in ODB compatible text format +! + open(unit=2,form='formatted',status='unknown',file=outfile) + headerline=emptyline + if(l_code_in_char)then + headerline='station@hdr:string longitude@hdr:real latitude@hdr:real hour@hdr:integer' + else + headerline='station@hdr:integer longitude@hdr:real latitude@hdr:real hour@hdr:integer' + endif + headerline=trim(headerline)//' total_number@body:real' + headerline=trim(headerline)//' msd@body:real' + headerline=trim(headerline)//' p_value@body:real' +! +! quantile bin frequency variable names: f00, f01 ... f99 +! + 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 + !print*,msd + !print*,ntimes_hour + if(l_code_in_char)then + write(dataline,'(A3,2F16.6,1X,I2,I7,2F16.6)')& + station_char,longitude,latitude,hour,ntimes_hour(j),msd(j),p_value(j) + else + write(dataline,'(I6,2F16.6,1X,I2,I7,2F16.6)')& + station_int,longitude,latitude,hour,ntimes_hour(j),msd(j),p_value(j) + endif + + do i=1,nquant+1 + write(frequency_value,'(F16.6)')frequency(i,j) + dataline=trim(dataline)//frequency_value + enddo + write(2,*)trim(dataline) + enddo +! +! Write the data line for the all-UTC statistics +! + if(l_code_in_char)then + write(dataline,'(A3,2F16.6,1X,I2,I7,2F16.6)')& + station_char,longitude,latitude,24,ntimes_hour_tot,msd(nhour+1),p_value(nhour+1) + else + write(dataline,'(I6,2F16.6,1X,I2,I7,2F16.6)')& + station_int,longitude,latitude,24,ntimes_hour_tot,msd(nhour+1),p_value(nhour+1) + endif + do i=1,nquant+1 + write(frequency_value,'(F16.6)')frequency(i,nhour+1) + dataline=trim(dataline)//frequency_value + enddo + write(2,*)trim(dataline) + + close(2) +! +!--------------------------------------------------------------------------------- +! +! Open the file for GrADS output and write its contents, for visualisation in GrADS? +! + if(grads_out)then + open(unit=11,form='unformatted',file=outfile_grads,access='DIRECT',& + recl=4,status='unknown') +! +! Again, first each UTC hour separately +! + do hour=hour1,hour2,hour_step + j=(hour-hour1)/hour_step + write(11,rec=1+j*(nquant+6))longitude + write(11,rec=2+j*(nquant+6))latitude + write(11,rec=3+j*(nquant+6))real(ntimes_hour(hour+1)) + write(11,rec=4+j*(nquant+6))msd(hour+1) + write(11,rec=5+j*(nquant+6))p_value(hour+1) + do i=1,nquant+1 + write(11,rec=5+j*(nquant+6)+i)frequency(i,hour+1) + enddo + enddo + ! + ! The all-UTC statistics + ! + write(11,rec=1+(j+1)*(nquant+6))longitude + write(11,rec=2+(j+1)*(nquant+6))latitude + write(11,rec=3+(j+1)*(nquant+6))real(ntimes_hour_tot) + write(11,rec=4+(j+1)*(nquant+6))msd(nhour+1) + write(11,rec=5+(j+1)*(nquant+6))p_value(nhour+1) + do i=1,nquant+1 + write(11,rec=5+(j+1)*(nquant+6)+i)frequency(i,nhour+1) + enddo + + close(11) + + endif ! if (grads_out) + + END program rank_histograms_one_station +! +!---------------------------------------------------------------- +! + subroutine update_frequencies & + (frequency,f,quantiles,nquant,miss) +! +! Find the location of the data value (f) within the +! quantile value table (quantiles, ascending order) +! and update the table of absolute quantile bin +! frequencies (frequency) +! +! frequency (in/out) = quantile bin frequencies +! f (in) = individual data value +! quantiles (in) = quantile values +! nquant (in) = number of quantiles +! miss (in) = missing value code + + implicit none + integer :: nquant + real :: quantiles(nquant) + real :: frequency(nquant+1) + real :: f + integer :: i1,i2,i,ind + real :: miss +! +! Quantiles should never be missing, but if they are, +! this should concern all quantiles simultaneously. +! Therefore, check just the 1st quantile + + if(abs(quantiles(1)).ge.abs(miss))then + write(*,*)'Quantiles missing!?' + stop + endif +! +! Find the position of f (ind) in the quantile table +! + if(f.lt.quantiles(1))then + ind=1 + else + if(f.ge.quantiles(nquant))then + ind=nquant+1 + else + i1=1 + i2=nquant + do while (i2.gt.i1+1) + i=(i1+i2)/2 + if(f.ge.quantiles(i))then + i1=(i1+i2)/2 + else + i2=(i1+i2)/2 + endif + enddo + ind=i1+1 + endif + endif +! +! Update the frequency table +! + frequency(ind)=frequency(ind)+1. + + return + end subroutine update_frequencies + +!---------------------------------------------------------------------- + +subroutine round_hour(yyear,mmonth,dday,hhour,year,month,day,hour,& + hour1,hour2,hour_step) +! +! Rounding the hour to the nearest hour within (hour1, hour1 + hour_step, ... hour2) +! +implicit none +integer :: yyear,mmonth,dday,hhour ! year, month, day and hour before rounding (IN) +integer :: year,month,day,hour ! year, month, day and hour after rounding (OUT) +integer :: hour1,hour2,hour_step ! target hours: hour1, hour1 + hour_step ... hour2 +integer :: ndays_month(12) ! number of days per month +integer :: hour_index +ndays_month=31 +ndays_month(4)=30 +ndays_month(6)=30 +ndays_month(9)=30 +ndays_month(11)=30 +ndays_month(2)=28 +if((mod(yyear,4).eq.0.and.mod(yyear,100).gt.0).or.(mod(yyear,400).eq.0))then + ndays_month(2)=29 +endif +year=yyear +month=mmonth +day=dday +! +! round the hour to the nearest output hour +! +hour_index=nint((hhour-hour1+0.)/hour_step) +hour=hour1+hour_step*hour_index + if(hour.ge.24)then ! hhour was rounded forward to next day + hour=hour1 + day=dday+1 + if(day.gt.ndays_month(month))then + day=1 + month=month+1 + if(month.gt.12)then + month=1 + year=yyear+1 + endif + endif + endif + if(hour.lt.0)then ! Hhour was rounded backward to previous day + hour=hour2 + day=dday-1 + if(day.eq.0)then + month=month-1 + if(month.eq.0)then + month=12 + year=yyear-1 + endif + day=ndays_month(month) + endif + endif + +return +end subroutine round_hour + + + diff --git a/SATELLITE/STATST/list_of_all_amsua_locations.txt b/SATELLITE/STATST/list_of_all_amsua_locations.txt new file mode 100644 index 0000000000000000000000000000000000000000..2dbd42037c1acfc54928098f15a96b6f1ea1c991 --- /dev/null +++ b/SATELLITE/STATST/list_of_all_amsua_locations.txt @@ -0,0 +1 @@ +a00 45.0 70.0 diff --git a/SATELLITE/STATST/produce_rank_histograms_all_locations.sh b/SATELLITE/STATST/produce_rank_histograms_all_locations.sh new file mode 100755 index 0000000000000000000000000000000000000000..05927d124ed18850d7d19a8c809873d06383d120 --- /dev/null +++ b/SATELLITE/STATST/produce_rank_histograms_all_locations.sh @@ -0,0 +1,306 @@ +#!/bin/bash +echo "Bash version ${BASH_VERSION}..." +#source ../set_env.sh +################################################################## +# +# Calculation of quantile space rank histograms and their +# summary statistics (MSD + p-value) for all stations in the +# RHARM data sample +# +# The following need to be given as arguments: +# +# - variable +# - first year of raw data +# - last year of raw data +# - first month of raw data (default 1) +# - last month of raw data (default 12) +# +# As the input for this script, the following files are needed: +# +# 1) List of stations +# 2) Raw simulation data at station coordinates +# (but here, the station observations are used as a surrogate of this file) +# 3) File including quantiles as a function of time of year +# 4) File including the rank histogram bootstrap MSD values +# +# Important notes: +# +# 0) It is assumed that all the files are in ODB format. However, +# their relevant parts are converted to text format for processing. +# 1) The list of station coordinates can be retrieved from the raw +# ODB file, but this part is commented here. +# 2) The rank histogram bootstrap MSD values must be available +# for the same length of simulation as given as arguments to this script. +# +# Execution (e.g): ./produce_rank_histograms_all_stations 2010 2010 1 12 +# +# Jouni Räisänen, July 2023 +# +################################################################## +# +# Arguments: +# +# 0. Variable +# +variable=$1 # btmp for AMSU-A Channel 5 +echo " Variable number: 194" +echo " Variale name: btmp" +echo " Variable AMSU-A channel: $variable" +# +# 1.-2: First and last year +# +year1=$2 +year2=$3 +let nyears=year2-year1+1 +# +# 3.-4: First and last month +# +month1="${4:-1}" +month2="${5:-12}" +# +# Find the length of the simulation period. +# It this is 12 months or more, it is rounded to the nearest integer +# number of years. Otherwise, the number of months is recorded. +# +let nmonths=12*nyears+month2-month1-11 +echo 'Number of months' $nmonths +# +n=$( echo "$nmonths / 12" | bc -l ) +nyears_round=$( echo "($n + 0.5) / 1" | bc ) +echo 'Number of years rounded' ${nyears_round} +# +if [ $nmonths -ge 12 ] +then + period=${nyears_round}yr +else + period=${nmonths}mon +fi +echo $period +# +# Add zero to ahead of $month1 and $month2 for naming of files if needed? +# +if [ $month1 -le 9 ] +then + month1='0'${month1} +fi +if [ $month2 -le 9 ] +then + month2='0'${month2} +fi +# +################################################################## +###module load grads +# +# Add odb_api to $PATH +# +# On Puhti +#PATH="${PATH}":/projappl/project_2001011/odb_api/bin +# +# On Lumi +#PATH="${PATH}":/projappl/project_465000454/ama/software/odb_api/bin +# +# On Lumi (asked Devaraju Narayanappa, CSC) +# to make in /project/project_465000454/devaraju/modules/LUMI/23.03/C +# module use /project/project_465000454/devaraju/modules/LUMI/23.03/C +# +################################################################## +# +# ----- File names. Hard-coded, at least this far ----- +# +# NB: it would be more efficient to have ODB files for individual +# stations instead for the ones that include everything for all stations! +# +# 1) Raw data from observations (USED or NOT???) +#echo " Directory with raw data from observations (pre-computed) ..." +#obs_dir=RHARM_vertical_differences_odb +#obs_dir=/scratch/project_465000454/ama/STATDATARADSOUND/RHARM_vertical_differences_odb +#obs_dir= # make some climatological means here, NOT USED??? +#echo " $obs_dir" + +# 2) Raw data from model simulation (this far, mimicked by observations) +echo " Directory with data from model simulations ..." +#sim_dir=RHARM_vertical_differences_odb +#sim_dir=/scratch/project_465000454/ama/STATDATARADSOUND/RHARM_vertical_differences_odb +#sim_dir=$final_output_dir +#sim_dir=/scratch/project_465000454/lautuppi/AMSUA_demo_ODB/test_odb_files #LT +sim_dir=/scratch/project_465000454/ama/AMSUA_demo_ODB/test_odb_files #AM +echo " $sim_dir" + +# 3) Pre-computed quantiles +echo " Directory with pre-computed quantiles ..." +#quant_dir="quantiles" +#quant_dir=/scratch/project_465000454/ama/STATDATARADSOUND/quantiles +#quant_dir=$scriptdir/STATST/pre-computed_data/quantiles #LT +quant_dir=/scratch/project_465000454/ama/STATDATASATELLITE/quantiles #AM +echo " $quant_dir" + +# 4) Pre-computed bootstrap MSD statistics +echo " Directory with pre-computed bootstrap MSD statistics ..." +#bootstrap_dir=bootstrap_statistics/${variable}_${period} +#bootstrap_dir=/scratch/project_465000454/ama/STATDATARADSOUND/bootstrap_statistics/${variable}_${period} +#quant_dir=$scriptdir/STATST/pre-computed_data/bootstrap_statistics #LT +bootstrap_dir=/scratch/project_465000454/ama/STATDATASATELLITE/bootstrap_statistics #AM +echo " $bootstrap_dir" + +# 5) Directory for results +echo " Directory for output results ..." +outdir=rank_histograms/${variable}_${year1}${month1}-${year2}${month2} + +#[ ! -d ${outdir} ] && mkdir ${outdir} +echo " $outdir" +if test -d "$outdir"; then + echo $outdir exists +else + echo $outdir does not exist + echo " Creating directory: $outdir" + mkdir -p $outdir +fi + +# 6) Directory for temporary results +echo " Directory for temporary output results ..." +outtmpdir=tmp_output + +echo " $outtmpdir" +if test -d "$outtmpdir"; then + echo $outtmpdir exists +else + echo $outtmpdir does not exist + echo " Creating directory: $outtmpdir" + mkdir -p $outtmpdir +fi + +################################################################## +# It is assumed that the list of stations has been pre-produced +# +echo " (Note) List includes all locations-areas with obs.data ..." +station_list=list_of_all_amsua_locations.txt +echo " location-area ID, longitude, latitude" +echo " List of selected locations-areas (with geographical coordinates): $station_list" + +################################################################ +echo " Compiling fortran-program to calculate rank histograms ..." +echo " fortran-programs/rank_histograms_one_location.f95" + +gfortran fortran-programs/rank_histograms_one_location.f95 -o rank_histograms_one_location.exe + +################################################################# +# +n_lines=$(cat ${station_list} | wc | awk NF=1) +line_number=1 +# +# Loop over all stations +# +while [ ${line_number} -le `expr $n_lines` ] +do +# head -`expr ${line_number}` ${station_list} | tail -1 > input.txt +# read station longitude latitude < input.txt + head -`expr ${line_number}` ${station_list} | tail -1 > input.txt + read station longitude latitude < input.txt + #read station longitude latitude n00 n12 < input.txt +echo " **********" +echo " Location-area ID, longitude, latitude: ${station} ${longitude} ${latitude}" + +################################################################# +# File names +# +#obs_file=${obs_dir}/station_${station}_var_2_diff_H_2_p_85000.odb # Not used??? + +echo " Paths to files with simulation, quantiles, and MSD bootstrap data ..." +sim_file=${sim_dir}/amsua_c5_2000.odb #odb +echo " sim_file: $sim_file" +quant_file=${quant_dir}/70N_45E_c5bt_quantiles.txt #odb +echo " quant_file: $quant_file" +msd_bootstrap_file=${bootstrap_dir}/70N_45E_c5bt_MSDstats_$year1.txt #odb +echo " msd_bootstrap_file: $msd_bootstrap_file" + +# Produce daily means +echo " Calling python script (daymean.py) to produce daily mean ..." +#LT odb sql 'select year,month,day,hour,modvalue' -i ${sim_file} -o tmp_output/sim_data_hourly +#LT python3 python/daymean.py tmp_output/sim_data_hourly tmp_output/sim_data +# AM updated below +odb sql 'select year,month,day,hour,modvalue' -i ${sim_file} -o $outtmpdir/sim_data_hourly +python3 python/daymean.py $outtmpdir/sim_data_hourly $outtmpdir/sim_data + +################################################################ +# Select the simulation data for the location-area +# (mimicked this far by observations!). +echo " Selecting the simulation data for AMSU-A location-area ID: ${station}" + +#odb_command="odb sql 'select year,month,day,hour,value where station=${station} and (year>=${year1} and year<=${year2}) and (hour=21 or hour=22 or hour=23 or hour=0 or hour=1 or hour=2 or hour=9 or hour=10 or hour=11 or hour=12 or hour=13 or hour=14)' -i ${sim_file} -o sim_data" +#odb_command="odb sql 'select year,month,day,hour,value where station=${station} and (year>=${year1} and year<=${year2})' -i ${sim_file} -o sim_data" +#eval ${odb_command} +# +############################################################### +# Select the quantiles for the nearest location-area +echo " Selecting the quantiles for AMSU-A location-area ID: ${station}" +# +#test -f quantile_selection && rm quantile_selection +#odb sql select \* -i ${quant_file} -o quantile_selection +# +################################################################ +# +# Get the distribution of pre-computed Mean Square Deviations from +# the msd_bootstrap_file. Note that these values depend on the length +# of the time period (included in the name of msd_bootstrap_file) +# and (to a smaller extent) the selected station. +# +#test -f msd_bootstrap_selection && rm msd_bootstrap_selection +#odb sql select \* -i ${msd_bootstrap_file} -o msd_bootstrap_selection +#test -f $sim_dir/70N_45E_c5bt_MSDstats_$year1.txt && rm $sim_dir/70N_45E_c5bt_MSDstats_$year1.txt + +############################################################### +# Produce the rank histogram for one location-area +echo " Producing the rank histogram for one AMSU-A location-area ID: ${station}" + +############################################################### +echo " Checking existence of file: rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} ..." +outfile=${outdir}/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} + +echo " $outfile" +if test -f "$outfile"; then + echo $outfile exists +else + echo $outfile does not exist + echo " Creating file: $outfile" + touch $outfile +fi + +# +./rank_histograms_one_location.exe < ${outdir}/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_all_locations.odb + +###################################################################### +echo " Removing unnecessary files ..." + +#rm input.txt +#rm msd_bootstrap_selection +#rm quantile_selection +#rm sim_data + +rm rank_histograms_one_location.exe + diff --git a/SATELLITE/STATST/produce_standard_plots_all_locations.sh b/SATELLITE/STATST/produce_standard_plots_all_locations.sh new file mode 100755 index 0000000000000000000000000000000000000000..a2c0db8b60f66caf03d4802fa4f631fe153e76ab --- /dev/null +++ b/SATELLITE/STATST/produce_standard_plots_all_locations.sh @@ -0,0 +1,319 @@ +#!/bin/bash +echo "Bash version ${BASH_VERSION}..." +echo "Python version" +which python +########################################################################## +# +# Produce the following plots: +# +# (a) individual data values (from a model simulation or a real station) +# against the background of quantiles from the observed distribution +# for the same station +# (b) rank histogram for the same station +# +# for all stations in the RHARM data sample, for T(2 m) - T(850 hPa) +# +# The following need to be given as arguments: +# +# - first year of raw data +# - last year of raw data +# - first month of raw data (default 1) +# - last month of raw data (default 12) +# +# As the input for this script, the following files are needed: +# +# 1) List of stations +# 2) Raw simulation data at station coordinates +# (but here, the station observations are used as a surrogate of this file) +# 3) File including quantiles as a function of time of year +# 4) File including the rank histogram bootstrap MSD values +# +# Important notes: +# +# 1) It is assumed that all the files are in ODB format. However, +# their relevant parts are converted to text format for processing. +# 2) The rank histogram bootstrap MSD values must be available +# for the same length of simulation as given as arguments to this script. +# +# Execution (e.g): ./produce_standard_plot_all_stations 2010 2010 +# +# ORIGINAL: +# Jouni Räisänen, July 2023 +# MODIFIED: +# Alexander Mahura, Sep-Oct-Nov 2023 +# Lauri Tuppi, Jan-Feb-Mar 2024 +# Alexander Mahura, Mar-Apr 2024 +# +################################################################## +# +# Arguments: +# +# The Variable is always the same +# +variable=$1 # btmp for AMSU-A Channel 5 +echo " Variable number: 194" +echo " Variale name: btmp" +echo " Variable AMSU-A channel: $variable" + +# +# 1.-2: First and last year +# +year1=$2 +year2=$3 +let nyears=year2-year1+1 +# +# 3.-4: First and last month +# +month1="${4:-1}" +month2="${5:-12}" +# +# Find the length of the simulation period. +# It this is 12 months or more, it is rounded to the nearest integer +# number of years. Otherwise, the number of months is recorded. +# +let nmonths=12*nyears+month2-month1-11 +echo 'Number of months' $nmonths +# +n=$( echo "$nmonths / 12" | bc -l ) +nyears_round=$( echo "($n + 0.5) / 1" | bc ) +echo 'Number of years rounded' ${nyears_round} +# +if [ $nmonths -ge 12 ] +then + period=${nyears_round}yr +else + period=${nmonths}mon +fi +echo $period +# +# Add zero to ahead of $month1 and $month2 for naming of files if needed? +# +if [ $month1 -le 9 ] +then + month1='0'${month1} +fi +if [ $month2 -le 9 ] +then + month2='0'${month2} +fi +# +################################################################## +# +# Add odb_api to $PATH +# On Puhti +#PATH="${PATH}":/projappl/project_2001011/odb_api/bin +# On Lumi +#export PATH="/projappl/project_465000454/ama/software/odb_api/bin:$PATH" +# +# On Lumi (asked Devaraju Narayanappa, CSC) +# to make in /project/project_465000454/devaraju/modules/LUMI/23.03/C +# module use /project/project_465000454/devaraju/modules/LUMI/23.03/C +# module load odb_api/0.18.1-cpeCray-23.03.lua +# module load python-climatedt/3.11.3-cpeCray-23.03.lua +# +# +################################################################## +# +# ----- Directory names. Hard-coded, at least this far ----- +# +# 1) Raw data from observations +#echo " Directory with raw data from observations (pre-computed) ..." +#obs_dir=RHARM_vertical_differences_odb +#obs_dir=/scratch/project_465000454/ama/STATDATARADSOUND/RHARM_vertical_differences_odb +#echo " $obs_dir" + +# 2) Raw data from model simulation (this far, mimicked by observations) +echo " Directory with data from model simulations (this far, mimicked by observatios) ..." +#sim_dir=RHARM_vertical_differences_odb +#sim_dir=/scratch/project_465000454/lautuppi/AMSUA_demo_ODB/test_odb_files #LT +sim_dir=/scratch/project_465000454/ama/AMSUA_demo_ODB/test_odb_files #AM +echo " sim_dr: $sim_dir" + +# 3) Pre-computed quantiles +echo " Directory for pre-computed quantiles ..." +#quant_dir=quantiles +#quant_dir=$scriptdir/STATST/pre-computed_data/quantiles +quant_dir=/scratch/project_465000454/ama/STATDATASATELLITE/quantiles #AM +echo " quant_dir: $quant_dir" + +# 3) Pre-computed rank histogram data +#rh_dir=rank_histograms/${variable}_${year1}${month1}-${year2}${month2} +echo " Directory for computed rank histogram data ..." +rh_dir=rank_histograms/${variable}_${year1}${month1}-${year2}${month2} +echo " rh_dir: $rh_dir" + +# 4) Directory for figures +echo " Directory for figures ..." +figure_dir=figures/standard_plots_${variable}_${year1}${month1}-${year2}${month2} +echo " $figure_dir" +if test -d "$figure_dir"; then + echo $figure_dir exists +else + echo $figure_dir does not exist + echo " Creating directory: $figure_dir" + mkdir -p $figure_dir +fi + +################################################################## +# It is assumed that the list of stations has been pre-produced +# +#echo " (Note) List includes all stations with obs/data points ..." +#station_list=list_of_all_radsound_stations.txt +#echo " (Note) List includes stations with at least 1000 obs/data points ..." +station_list=list_of_all_amsua_locations.txt +echo " locationID, longitude, latitude" +echo " List of AMSU-A locations (with geographical coordinates): $station_list" + +################################################################ +# +echo " Compiling the Fortran program needed for creating the figures ..." +echo " fortran-programs/plots_for_one_location.f95" + + gfortran fortran-programs/plots_for_one_location.f95 -o plots_for_one_location.exe + +################################################################# +# +# Count the number of lines in the station list +# +n_lines=$(cat ${station_list} | wc | awk NF=1) +# +# Skip the first line which contains no station_ID +# +line_number=1 #2 +# +# Loop over all stations +# +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 " AMSU-A location-areaID, longitude, latitude: ${station} ${longitude} ${latitude}" +# +#sim_file=${sim_dir}/station_${station}_var_2_diff_H_2_p_85000.odb +sim_file=tmp_output/sim_data +#quant_file=${quant_dir}/quantiles_${station}_T2-T850.odb +quant_file=${quant_dir}/${station}_c5bt_quantiles.txt #odb +rh_file=${rh_dir}/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2}.odb +# +################################################################ +# Select the simulation data for the station (mimicked this far by observations!) +echo " Selecting the simulation data for AMSU-A location: ${station}" +# +#odb_command="odb sql 'select year,month,day,hour,value where (year>=${year1} and year<=${year2})' -i ${sim_file} -o sim_data" +#eval ${odb_command} + +################################################################ +# Select the rank histogram data for the station +echo " Selecting rank histogram data for AMSU-A location: ${station}" +# +#echo $rh_file +#echo tmp_output/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} +#odb sql select \* where station=${station} -i ${rh_file} -o tmp_output/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} +odb ls $rh_file | grep ${station} > tmp_output/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} +# +############################################################### +# Select the quantiles for the station +# 00 and 12 UTC are picked separately, since +# select \* does not allow for parentheses (very strange) +echo " Selecting the quantiles for AMSU-A location: ${station}" +# +#rm quantile_selection +#odb sql select \* -i ${quant_file} -o quantile_selection +#rm quantile_selection_* +#odb sql select \* where hour=0 and station=${station} -i ${quant_file} -o quantile_selection_00 +#odb sql select \* where hour=12 and station=${station} -i ${quant_file} -o quantile_selection_12 +#cat quantile_selection_* > quantile_selection +# +################################################################ +# +# Run the fortran progam that prepares the plotting +# +############################################################### +# +echo $quant_file +echo tmp_output/rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} +#echo +./plots_for_one_location.exe < tmp_output/standard_plot.cfg +[data] +variable=${variable} +station_id=${station} +longitude=${longitude} +latitude=${latitude} +year_beg=${year1} +year_end=${year2} +month_beg=${month1} +month_end=${month2} + +[input] +sim_data=tmp_output/time_series_${variable}_${station}_${year1}${month1}-${year2}${month2}.txt +quantile_sel=tmp_output/quantiles_${variable}_${station}.txt +rank_hists=tmp_output/rank_histogram_${variable}_${station}_${year1}${month1}-${year2}${month2}.txt + +[output] +fig_name=${figure_dir}/standard_plot_${variable}_${station}_${year1}${month1}-${year2}${month2}.png +EOF + +################################################################ +# +echo " Calling python to plot quantiles rank histogram for AMSU-A location-area: ${station}" +python3 python/plot_quantiles_rankhist_amsua.py tmp_output/standard_plot.cfg + +################################################################ +# Remove unnecessary files +echo " Removing unnecessary temporary files ..." + +if [ 0 -eq 1 ]; then + rm input.txt + rm vrange_* + rm msd_and_p-value + rm quantile_selection + rm rank_histograms_${variable}_${station}_${year1}${month1}-${year2}${month2} + rm sim_data + ##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 + rm standard_plot.cfg + rm time_series_${variable}_${station}_${year1}${month1}-${year2}${month2}.txt + rm quantiles_${variable}_${station}.txt + rm rank_histogram_${variable}_${station}_${year1}${month1}-${year2}${month2}.txt +fi +# +((line_number++)) +done + +################################################################ +# The Fortran executable is not needed any more: +# +exit +rm rank_histogram.ctl +rm time_series.ctl +rm quantile_selection_* +rm msd_and_p-val* +rm plots_for_one_station +rm coordinates diff --git a/SATELLITE/STATST/python/daymean.py b/SATELLITE/STATST/python/daymean.py new file mode 100644 index 0000000000000000000000000000000000000000..1151f143f8e09c339fb96d0930fea6464b22097e --- /dev/null +++ b/SATELLITE/STATST/python/daymean.py @@ -0,0 +1,26 @@ +import numpy as np +import sys + +#hourly_data_file = '/scratch/project_465000454/lautuppi/AMSUA_demo_ODB/amsua_c5_2000_test.txt' +hourly_data_file = sys.argv[1] + +#daily_data_file = 'sim_data' +daily_data_file = sys.argv[2] + +hourly_data = np.loadtxt(hourly_data_file, skiprows=1) + +shape = np.shape(hourly_data) +#print(shape) +daily_data = np.zeros((int(shape[0]/24),int(shape[1]))) + +for i in range(shape[0]): # set metadata + if i%24 == 0: + #print(hourly_data[i,:]) + daily_data[int(i/24),0] = int(hourly_data[i,0]) + daily_data[int(i/24),1] = int(hourly_data[i,1]) + daily_data[int(i/24),2] = int(hourly_data[i,2]) + daily_data[int(i/24),4] = np.mean(hourly_data[i:i+23,4]) + +#print(daily_data[:,4]) +#daily_data +np.savetxt(daily_data_file,daily_data,fmt=['%d', '%d', '%d', '%d', '%.8f']) diff --git a/SATELLITE/STATST/python/plot_p_values_map_amsua.py b/SATELLITE/STATST/python/plot_p_values_map_amsua.py new file mode 100644 index 0000000000000000000000000000000000000000..ae5fac6dc31cf43a10357c00f5e7cd2fd361275d --- /dev/null +++ b/SATELLITE/STATST/python/plot_p_values_map_amsua.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python + +''' +Written (original non-python) by Jouni Räisänen, 2023-10-07 +Written (python) & Updates by Madeleine Ekblom, 2023-04-06 +Updates by Lauri Tuppi, 2024-02/03-... +Updates by Alexander Mahura, 2024-03/04-... + +LT: Plots p values on a map restricted to an area limited by lons (19, 32) and lats (59, 71) +AM: Plots p values from quantile bin rank histogram analysis to global maps (scale 0-1, only low values indicate significance). +P-values are plotted for ? which exactly time UTC? + +Example: +LT: +$ python3 plot_p_values_map.py p_values_as_text_file +AM: +python3 plot_p_values_map_rank_hist_amsua.py p_values_as_text_file timestring figname + +p_values_as_text_file is a text file containing "satellite area ID" (stationID in SYNOP and RADSOUND), longitude, latitude, +and p value with one header row, timestring indicates the simulation period used. It is only used in the map headers. + +''' + +import numpy as np +import matplotlib.pyplot as plt +import sys +import cartopy.crs as ccrs +from matplotlib.lines import Line2D + +def read_data(p_values_file): + ''' + Description: Reads in file containing p values at different stations + Input: file containing p values with "satellite area ID" (sid), longitude (lon), latitude (lat), and p value (p). Header = 1 row. + Output: structure numpy array with sid, lon, lat, p + ''' + + #p_values_data = np.loadtxt(p_values_file, skiprows=1, dtype={'names': ('sid', 'lon', 'lat', 'p'), 'formats':('S3', 'f4', 'f4', 'f4')}) #AM + p_values_data = np.loadtxt(p_values_file, skiprows=1, dtype={'names': ('sid', 'lon', 'lat', 'p'), 'formats':('i4', 'f4', 'f4', 'f4')}) #LT + + return p_values_data + +def plot_p_values(p_values): + ''' + Description: plot p values on a map restricted to an area limited by lons (19, 32) and lats (59, 71) + Input: numpy array + Output: saves a figure in the running directory + ''' + + # define lat and lon for global map + lon_min=-180 + lon_max=180 + lat_min=-90 + lat_max=90 + + # define colour and limits for the dots + colors = ['darkviolet', 'blue', 'skyblue', 'lime', 'yellow', 'orange', 'red', 'grey'] + limits = [0,0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,1] + + # set up figure + #fig = plt.figure(figsize=(12,4.5)) #LT + fig = plt.figure(figsize=(12,11)) #AM + + # plotting details: # AM + fs1=20 # fontsize for title + fs2=12 # fontsize for labels, subplot titles + fs3=10 # fontsize for figure explanation + + # add subplot with map, for 12 UTC data /as plotted daily mean values/ + ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree()) + ax.coastlines() + ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) + ax.gridlines(draw_labels=True) + + # basic plot, not changing colors of p: + ax.scatter(p_values['lon'], p_values['lat'], c=p_values['p']) + + #loop over limits: + legend_elements = [] + for n in range(8): + # find indices for p-values + p_mod = 0.5+1.0001*(p_values['p']-0.5) # AM + #p_ind = ((p_values['p'] > limits[n]) & (p_values['p'] <= limits[n+1])) + p_ind = ((p_mod > limits[n]) & (p_mod <= limits[n+1])) # AM + # scatter plot with dots coloured based on p-values + ax.scatter(p_values['lon'][p_ind], p_values['lat'][p_ind], c=colors[n]) + # add to legend_elements: colour and limit + legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1]))) + + # add legend to lower left corner of map + ax.legend(handles=legend_elements, loc='lower left') # AM + #ax.legend(handles=legend_elements, loc='center left', bbox_to_anchor=(1.1, 0.5)) # LT + + # add subplot title + ax.set_title('Daily Mean', fontsize=fs2, fontweight='bold') + + # explanation for figure: AM + text_expl = r'$\bf{p-value}$'+': probability that the MSD statistics for the quantile rank \nhistogram would be exceeded by pure chance when comparing \na perfect model simulation with observations.' + # add text box with explanation to the bottom of the figure: AM + ax.text(0.2, 1.25, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes) + + # explanation of brightness temperature ... AM + #text_btmp= r'$\bf{Brightness\ temperature}$'+': measurement of radiance of microwave\nradiation traveling upward from top of the atmosphere to satellite\n(expressed in units of temperature of equivalent black body)' + text_btmp= r'$\bf{Brightness\ temperature}$'+': Measure of upwelling microwave radiation\nat top of the atmosphere (expressed in units of temperature of equivalent black body)' + # add text box with explanation to the bottom of the figure + ax.text(0.2, -0.25, text_btmp, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes) + + # add figure title: AM + satellite_and_channel='AMSU-A Ch.5 : ' + fig.suptitle(satellite_and_channel + 'SATELLITE areas:\np-values for brightness temperature (btmp)\nquantiles in ' + timestring, fontsize=fs1) #, fontweight='bold') + + # save figure + plt.savefig(figname, dpi=300) #AM + #plt.savefig('p_values.png', dpi=300) #LT + + +def main(p_values_file): + ''' + Main: (1) reads in p values, and (2) plots p-values on a global map + ''' + p_values_data = read_data(p_values_file) + + plot_p_values(p_values_data) + +if __name__=='__main__': + p_values_file = sys.argv[1] + timestring = sys.argv[2] + figname = sys.argv[3] #AM + + main(p_values_file) diff --git a/SATELLITE/STATST/python/plot_quantiles_rankhist_amsua.py b/SATELLITE/STATST/python/plot_quantiles_rankhist_amsua.py new file mode 100644 index 0000000000000000000000000000000000000000..116389b8a3cee1d3f5f0afd2b399a5d8b2ac9793 --- /dev/null +++ b/SATELLITE/STATST/python/plot_quantiles_rankhist_amsua.py @@ -0,0 +1,358 @@ +#!/usr/bin/env python + +''' +Originally written (for GRADs) by Jouni Räisänen (script for plotting quantiles and rank histogram) +Written in python by Madeleine Ekblom, 2023-05-10 +Updates by Lauri Tuppi, 2024-02/03-... +Updates by Alexander Mahura, 2024-03/04-... + +Example: +$ python3 plot_quantiles_rankhist.py example.cfg + +example.cfg is a text file containing: +[data] +variable=39 +station_id=102019 +longitude=23.576000 +latitude=68.602997 +year_beg=2010 +year_end=2012 +month_beg=01 +month_end=12 + +[input] +sim_data=time_series_102019_2010-2012.txt +quantile_sel=quantiles_102019.txt +rank_hists=rank_histogram_102019_2010-2012.txt + +[output] +fig_name=standard_plot_102019_2010-2012_python.png + +''' + +import numpy as np +import matplotlib.pyplot as plt +import sys +import configparser +import datetime +from matplotlib.dates import drange, MonthLocator, DateFormatter +import matplotlib.lines as mlines +import matplotlib.patches as mpatches +import matplotlib.ticker as ticker + +config = configparser.ConfigParser() + +def plot_data2(quantiles_data, rank_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end,station_id, savefig=False): + ''' + quantiles data give the lines and the time series data give the dots + ''' + # set up date array consisting of all days in beg year: + days = np.arange(1,366) +# day1 = datetime.date(year_beg, 1, 1) +# day2 = datetime.date(year_beg+1, 1, 1) + day1 = datetime.date(2016, 1, 1) + day2 = datetime.date(2016, 12, 31) + delta = datetime.timedelta(days=1) + dates = drange(day1, day2, delta) + + # arrays for quantile names, axes titles, hours to plot + #q_names = ['q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'] + sub_titles = ['Daily Mean'] + hours = [12] + + # calculate number of years: + nyears = year_end - year_beg + 1 + + # set up figure + fig = plt.figure(figsize=(10,10), layout='constrained') + #spec = fig.add_gridspec(2,5) + spec = fig.add_gridspec(3,5) + fig.suptitle(fig_title, fontsize=20) + + # counter + c = 0 + + # plot quantiles/time series for daily means: +# for i in range(1): +# for j in range(1): + # set up axis: title, xaxis + ax = fig.add_subplot(spec[0,0:5]) + ax.set_title(sub_titles[c]) + ax.set_xlim(dates[0], dates[-1]) + ax.xaxis.set_major_locator(MonthLocator()) + ax.xaxis.set_major_formatter(DateFormatter('%b')) + + ####################################################### + # legend elemnents for quantile data, one legend for model and one legend for observations + exp_id='a0h3' + color='blue' + alpha = [0.5, 0.3, 0.1] + linecolor = 'red' + ndigits = 3 + fs2 = 14 + fs3 = 12 + fs4 = 10 + scatter_size = 5 + legend_elements1 = [mlines.Line2D([],[], marker='.', linestyle='', color='k', label='Model data (' + exp_id + ')')] + legend_elements2= [mpatches.Patch(color=color, alpha=alpha[0], label='25-50%'), + mpatches.Patch(color=color, alpha=alpha[1], label='10-25%'), + mpatches.Patch(color=color, alpha=alpha[2], label='1-10%'), + mlines.Line2D([],[], color=linecolor, label='50%'), + mpatches.Patch(color=color, alpha=alpha[0], label='50-75%'), + mpatches.Patch(color=color, alpha=alpha[1], label='75-90%'), + mpatches.Patch(color=color, alpha=alpha[2], label='90-99%')] + + # text box for empty data: + text_no_data = "No data available" + + fig.suptitle(fig_title, fontsize=18) #20) + # counter + c = 0 + + # set up axis: title, xaxis, x- and y-labels + ax.set_title(sub_titles[c], fontsize=fs2) + ax.set_xlim(dates[0], dates[-1]) + ax.xaxis.set_major_locator(MonthLocator()) + ax.xaxis.set_major_formatter(DateFormatter('%b')) + ax.set_xlabel('Time', fontsize=fs2) + #ax.set_ylabel('$\Delta$T = 2t - t850 ($^\circ$C)', fontsize=fs2) + ax.set_ylabel('Brightness temperature, btmp ($^\circ$K)', fontsize=fs2) + + # quantile data: + # find quantile data for daily means + # qh_ind = (quantiles_data['hour'][:] == hours[c]) + + # get quantile data, qXX_data = quantile XX data at hour qh_ind +# q01_data = quantiles_data[q_names[0]][qh_ind] +# q10_data = quantiles_data[q_names[1]][qh_ind] +# q25_data = quantiles_data[q_names[2]][qh_ind] +# q50_data = quantiles_data[q_names[3]][qh_ind] +# q75_data = quantiles_data[q_names[4]][qh_ind] +# q90_data = quantiles_data[q_names[5]][qh_ind] +# q99_data = quantiles_data[q_names[6]][qh_ind] + q01_data = quantiles_data[:365,3]#-273.15 + q10_data = quantiles_data[:365,4]#-273.15 + q25_data = quantiles_data[:365,5]#-273.15 + q50_data = quantiles_data[:365,6]#-273.15 + q75_data = quantiles_data[:365,7]#-273.15 + q90_data = quantiles_data[:365,8]#-273.15 + q99_data = quantiles_data[:365,9]#-273.15 + + # solid line at q50: + ax.plot(dates, q50_data, c=linecolor, ls='-') + + # shading will be the most opaque close to 50 and more transparent close to 1/99 + + # fill area btw q25-q50, q50-q75: alpha = alpha[0] + ax.fill_between(dates, q25_data[:], q50_data[:], color=color, alpha=alpha[0]) + ax.fill_between(dates, q50_data[:], q75_data[:], color=color, alpha=alpha[0]) + + # fill area btw q10-q25, q75-q90 + ax.fill_between(dates, q10_data[:], q25_data[:], color=color, alpha=alpha[1]) + ax.fill_between(dates, q75_data[:], q90_data[:], color=color, alpha=alpha[1]) + # fill area btw q01-q10, q90-q99 + ax.fill_between(dates, q01_data[:], q10_data[:], color=color, alpha=alpha[2]) + ax.fill_between(dates, q90_data[:], q99_data[:], color=color, alpha=alpha[2]) + + # plot time series data: + # find time series where hour = 0, 12 + th_ind = (time_series_data['hour'][:] == hours[c]) + t_data = time_series_data['value'][th_ind] # all years + + # check if data available for given station and time: check that the length of array t_data > 0 + if (len(t_data) > 0): + for n in range(nyears): + # plot data with black, small dots: + ax.scatter(dates, t_data[n*365:(n+1)*365], s=scatter_size, c='k', marker='.') + else: + ax.text(0.5, 0.5, text_box, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # increase counter by 1 + c = c + 1 + ########################################## + + + # plot rank histogram data: + ax2 = fig.add_subplot(spec[1,:]) + #print("#######################") + #print(rank_data) + #ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[1,4]) + ' p-value:' + str(rank_data[1,5])) + #ax2.bar(np.arange(100), rank_data[1,6:]*100) + #ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[3]) + ' p-value:' + str(rank_data[4])) + ax2.set_title(('Rank histogram\n' + r'$\bf{MSD}$: ' + str(rank_data[3]) +' '+ r'$\bf{ p-value}$: '+ str(rank_data[4])), fontsize=fs2) + ax2.bar(np.arange(100), rank_data[5:]*100) + ax2.set_xlim(0,100) + ax2.axhline(y=1) + ax2.set_ylabel('Relative frequency', fontsize=fs2) + ax2.set_xlabel('Percentile', fontsize=fs2) + + ax3 = fig.add_subplot(spec[2,1:2]) + ax3.axis('off') + # add two legends: model data, and observations (with two columns) + #rank_hist_mean = r'$\bf{MSD}$: ' + str(np.round(rank_data[4], ndigits)) +' '+ r'$\bf{ p-value}$: ' + str(np.round(rank_data[5], ndigits)) + # comment, because !!! doublication of info on top of graph + #rank_hist_mean = r'$\bf{MSD}$: ' + str(np.round(rank_data[3], ndigits)) +' '+ r'$\bf{ p-value}$: ' + str(np.round(rank_data[4], ndigits)) + #ax3.text(0.25,0.65, rank_hist_mean, bbox={'pad':10, 'facecolor':'white'}, ha='center', fontsize=fs3) # 0.1, 0.9 + #### + leg1 = ax3.legend(loc='upper center', handles=legend_elements1, fontsize=fs4) + leg2 = ax3.legend(loc='lower center', handles=legend_elements2, ncols=2, fontsize=fs4, title='Quantiles for observations', title_fontsize=fs4) + ax3.add_artist(leg1) + ax3.add_artist(leg2) + + ax4 = fig.add_subplot(spec[2,2:4]) + ax4.axis('off') + ###rank_hist_mean = r'$\bf{MSD}$: ' + str(np.round(rank_data[4], ndigits)) +' '+ r'$\bf{ p-value}$: ' + str(np.round(rank_data[5], ndigits)) + # add text boxes with rank histogram data + ###ax4.text(0.1,0.9, rank_hist_mean, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # figure explanation: + #expl_data = 'Mean Square Deviation (' + r'$\bf{MSD}$' + ') relative to flat histogram\n' + r'$\bf{p-value}$' + ' of MSD relative to sampling distribution\n' + expl_data = 'Mean Square Deviation (' + r'$\bf{MSD}$' + ') relative to flat histogram\n' + r'$\bf{p-value}$' + ' of MSD relative to sampling distribution\n' + expl_quantile='Quantiles are estimated based on observations\n in the area ' + station_id + ' (period: 2010-2020)\n' + expl_model='Model data from experiment (' + exp_id + ')\n interpolated to centre location of the area ' + station_id + + fig_text = expl_data + expl_quantile + expl_model + ax4.text(0.2,0.5, fig_text, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs4) # 0.1, 0.2 + + # explanation of brightness temperature ... AM + expl_btmp= r'$\bf{Brightness\ temperature}$'+': measurement of radiance\nof microwave radiation traveling upward from top\nof the atmosphere to satellite (expressed in units\nof temperature of nequivalent black body)' + expl_btmp= r'$\bf{Brightness\ temperature}$'+': Measure of upwelling\nmicrowave radiation at top of the atmosphere\n(expressed in units of temperature of\nequivalent black body)' + # add text box with explanation to the bottom of the figure + ax4.text(0.2,0.1, expl_btmp, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs4) #, transform=ax.transAxes) + + if savefig: + plt.savefig(fig_name, dpi=300) + + +def read_quantiles_data(quantile_sel_file): + print(quantile_sel_file) + ''' + Reads quantile data from text file and returns a structured numpy array + Input: text file + Output:a structured numpy array: sid, lon, lat, day_of_year, hour, q01, q10, q25, q50, q75, q90, q99 + ''' + # Header line contains: + #station@hdr:integer longitude@hdr:real latitude@hdr:real day_of_year@hdr:integer hour@hdr:integer q01@body:real q10@body:real q25@body:real q50@body:real q75@body:real q90@body:real q99@body:real + + #quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, usecols=[0,1,2,3,4] dtype={'names':('sid', 'lon', 'lat', 'day_of_year', 'hour', 'q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'), 'formats': ('S3', 'f4', 'f4', 'i4', 'i4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4')}) + #quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, usecols=[1,2,3,4,15,30,55,80,95,104], dtype={'names':('lon', 'lat', 'hour', 'q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99'), 'formats': ('f4', 'f4', 'i4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4')}) + #quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, usecols=[1,2,3,6,15,30,55,80,95,104]) + quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, usecols=[0,1,3,4,5,6,7,8,9,10]) + #quantile_data = np.loadtxt(quantile_sel_file, skiprows=1, usecols=[1,2,4,5,6,7,8,9,10,11]) + #print(np.shape(quantile_data)) + return quantile_data + +def read_rank_data(rank_hists_file): + ''' + Reads rank histogram data binned into a 100 bins: 0-1%, 1-2%, ..., 99-100% + Input: text file + Output: numpy array (unstructured) first columns contain station id, longitude, latitude, hour, msd, p_value, and then follows the bins f00 = 0-1%, ..., f99=99-100% + ''' + + #Header line contains: + # station@hdr:integer longitude@hdr:real latitude@hdr:real hour@hdr:integer msd@body:real p_value@body:real f00@body:real f01@body:real f02@body:real f03@body:real f04@body:real f05@body:real f06@body:real f07@body:real f08@body:real f09@body:real f10@body:real f11@body:real f12@body:real f13@body:real f14@body:real f15@body:real f16@body:real f17@body:real f18@body:real f19@body:real f20@body:real f21@body:real f22@body:real f23@body:real f24@body:real f25@body:real f26@body:real f27@body:real f28@body:real f29@body:real f30@body:real f31@body:real f32@body:real f33@body:real f34@body:real f35@body:real f36@body:real f37@body:real f38@body:real f39@body:real f40@body:real f41@body:real f42@body:real f43@body:real f44@body:real f45@body:real f46@body:real f47@body:real f48@body:real f49@body:real f50@body:real f51@body:real f52@body:real f53@body:real f54@body:real f55@body:real f56@body:real f57@body:real f58@body:real f59@body:real f60@body:real f61@body:real f62@body:real f63@body:real f64@body:real f65@body:real f66@body:real f67@body:real f68@body:real f69@body:real f70@body:real f71@body:real f72@body:real f73@body:real f74@body:real f75@body:real f76@body:real f77@body:real f78@body:real f79@body:real f80@body:real f81@body:real f82@body:real f83@body:real f84@body:real f85@body:real f86@body:real f87@body:real f88@body:real f89@body:real f90@body:real f91@body:real f92@body:real f93@body:real f94@body:real f95@body:real f96@body:real f97@body:real f98@body:real f99@body:real + + #rank_hists = np.loadtxt(rank_hists_file, skiprows=1, usecols=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105]) + rank_hists = np.loadtxt(rank_hists_file, skiprows=1) + + return rank_hists + +def read_time_series_data(sim_data_file): + ''' + Reads time series data + Input: text file + Output: structured numpy array: sid, lon, lat, year, day_of_year, hour, value + ''' + + # Header line contains + # station@hdr:integer longitude@hdr:real latitude@hdr:real year@hdr:integer day_of_year@hdr:integer hour@hdr:integer value@body:real + + #time_series_data = np.loadtxt(sim_data_file, skiprows=1, usecols=[1,2,3,4,5,6], dtype={'names': ('lon', 'lat', 'year', 'day_of_year', 'hour', 'value'), 'formats': ('f4', 'f4', 'i4', 'i4', 'i4', 'f4')}) + time_series_data = np.loadtxt(sim_data_file, skiprows=1, dtype={'names': ('lon', 'lat', 'year', 'day_of_year', 'hour', 'value'), 'formats': ('f4', 'f4', 'i4', 'i4', 'i4', 'f4')}) + + return time_series_data + +def main(config_file): + ''' + Main: + reads config files + reads input data: time_series_data, quantiles_data, rank_hist_data + calls plotting function and saves figure with figure name given in config file + ''' + ## Read from config file ## + config.read(config_file) + + # data to be plotted + variable = config['data']['variable'] + station_id = config['data']['station_id'] + longitude = config['data']['longitude'] + latitude = config['data']['latitude'] + year_beg = int(config['data']['year_beg']) + year_end = int(config['data']['year_end']) + month_beg = config['data']['month_beg'] + month_end = config['data']['month_end'] + + # input files + sim_file = config['input']['sim_data'] + quantile_file = config['input']['quantile_sel'] + rank_hist_file = config['input']['rank_hists'] + + # output files + fig_name = config['output']['fig_name'] + + ## Read input data ## + time_series_data = read_time_series_data(sim_file) + quantiles_data = read_quantiles_data(quantile_file) + rank_hist_data = read_rank_data(rank_hist_file) + + ## Plot data ## +# if (month_beg < 10): +# str_month_beg = '0' + str(month_beg) +# else: +# str_month_beg = str(month_beg) + +# if (month_end < 10): +# str_month_end = '0' + str(month_end) +# else: +# str_month_end = str(month_end) + + # round lon, lat to three digits: + ndigits = 3 + lon = np.round(float(longitude), ndigits) + lat = np.round(float(latitude), ndigits) + + # define lon/lat text strings for title of the form lat: XX N/S, lat: XX E/W + if (lat > 0): + # if lat positive -> N + lat_text = ' (lat: ' + str(lat) + '$^\circ$N, ' + else: + # if negative -> S + lat_text = ' (lat ' + str(-lat) + '$^\circ$S, ' + + if (0 <= lon) and (lon <= 180): + # lon btw 0 an 180 -> E + lon_text = ' lon: ' + str(lon) + '$^\circ$E)' + else: + if lon < 0: + # lon < 0 --> - lon W + lon_text = ' lon: ' + str(-lon) + '$^\circ$W)' + else: + # 180 < lon < 360 --> - (lon - 360) W + lon = -(lon - 360) + lon_text = ' lon: ' + str(lon) + '$^\circ$W)' + + satellite_and_channel='AMSU-A Ch.5 : ' + fig_title = satellite_and_channel + 'SATELLITE area: ' + station_id + lat_text + lon_text + '\nModel data from ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end #AM + + #fig_title = 'AMSU-A : SATELLITE area: ' + station_id + ', ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + '\nLat=' + latitude + ' Lon=' + longitude #LT + plot_data2(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, station_id, True) + + +if __name__=='__main__': + if (len(sys.argv) < 2): + sys.exit("Error: config file must be added as 2nd argument") + elif (len(sys.argv) > 2): + sys.exit("Error: only add config file as argument") + config_file = sys.argv[1] + + main(config_file) diff --git a/SATELLITE/STATST/python/plot_rankhist_sum_all_locations_amsua.py b/SATELLITE/STATST/python/plot_rankhist_sum_all_locations_amsua.py new file mode 100644 index 0000000000000000000000000000000000000000..fb3a137e4ba07c70e3f397216e4fe24976ab84a0 --- /dev/null +++ b/SATELLITE/STATST/python/plot_rankhist_sum_all_locations_amsua.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python + +''' +Originally written (for GRADs) by Jouni Räisänen (script for plotting rank histogram summart statistics) +Written in python by Madeleine Ekblom, 2023-09-12 +Updates by Lauri Tuppi, 2024-02/03-... +Updates by Alexander Mahura, 2024-03/04-... + +Example: +$ python3 plot_rank_hist_sum_all_stations.py rh_summary_file nstat p01 p1 p5 max_freq_p max_freq_q + +''' + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) +import sys + +def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q): + ''' + Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q + Description: Plots p and q frequencies and saves the figure with figname + ''' + # add more variables and units when needed + variable_shortnames={'5': 'btmp'} # clarify: ? is it number 5 in sat.data + variable_units={'5': '$^\circ$K'} # clarify: ? is it in deg C or deg K + + # variables to add, shortnames and units need to be checked + # total amount of clouds, 91 + # sea level pressure, 108 + # 1-hour precipitation, 80 + # relative humidity, 58 + # 10-minute precipitation intensity, 999 + # snow depth, 71 + # 2m temp, 39 --> 2t, degrees Celcius + # dew point temperature, 40 + # visibility, 6 + # wind direction, 111 + # max wind gust in last 10 mins,261 + # surface wind speed, 221 + + variable='5' #'39' 2t for synop data # ! change to be an input value, as string + + # title + satellite_and_channel='AMSU-A Ch.5 : ' + title=satellite_and_channel + 'SATELLITE areas:\n Rank histogram summary statistics\n for brightness temperature (' + variable_shortnames[variable] + ')' + + # text box information + text_box = 'Number of areas: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01: ' + str(p1) + '\nFrequency, p<0.05: ' + str(p5) + + # these can also be input arguments if the number of bins are not constant + number_of_p_bins=20 + number_of_q_bins=100 + + # set up figure + fig = plt.figure(figsize=(8,12)) + gs = fig.add_gridspec(3, hspace=0.5) + ax1, ax2, ax3 = gs.subplots(sharex=False, sharey=False) + + # plotting details + fs1=20 # fontsize for title + fs2=14 # fontsize for labels, subplot titles, textbox + fs3=10 # fontsize for text explanation + + # Plot p values + ax1.set_title('Normalized p-value frequencies', fontsize=fs2) + ax1.set_xlabel('p-value', fontsize=fs2) + ax1.set_ylabel('Relative frequency', fontsize=fs2) + ax1.set_ylim([0,max_freq_p]) + ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) # add +0.5 to center the bars + ax1.xaxis.set_major_locator(MultipleLocator(2)) + ax1.xaxis.set_minor_locator(MultipleLocator(1)) + p_xlabels = ax1.get_xticks() + ax1.set_xticks(p_xlabels, p_xlabels/number_of_p_bins) + ax1.set_xlim([0,number_of_p_bins]) + + # Plot q values + ax2.set_title('Normalized quantile frequencies', fontsize=fs2) + ax2.set_xlabel('Quantile (%)', fontsize=fs2) + ax2.set_ylabel('Relative frequency', fontsize=fs2) + ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) # add +0.5 to center the bars + ax2.xaxis.set_major_locator(MultipleLocator(10)) + ax2.xaxis.set_minor_locator(MultipleLocator(5)) + ax2.set_ylim([0,max_freq_q]) + ax2.set_xlim([0,number_of_q_bins]) + + # remove axis from ax3 + ax3.axis('off') + + # add text box with additional infromation: + ax3.text(0.3, 0.8, text_box, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs2) + + # explanation for figure + text_pvalue=r'$\bf{Normalized\ p-values}$'+': Frequency distribution of p-values at\nall individual stations, normalized to give an expected value of 1.\n' + text_quantile=r'$\bf{Normalized\ quantiles}$'+': Frequency distribution of model data\nrelative to quantiles of observations, averaged over all stations\nand divided by the expected value of 0.01.' + text_expl = text_pvalue + text_quantile + # add text box with explanation to the bottom of the figure + ax3.text(0.15,0.3, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # explanation of brightness temperature ... AM + text_btmp= r'$\bf{Brightness\ temperature}$'+': Measure of upwelling microwave radiation\nat top of the atmosphere (expressed in units of temperature\nof equivalent black body)' + #text_btmp= r'$\bf{Brightness\ temperature}$'+': Measurement of radiance of microwave\nradiation traveling upward from top of the atmosphere to satellite\n(expressed in units of temperature of equivalent black body)' + # add text box with explanation to the bottom of the figure + ax3.text(0.15,-0.1, text_btmp, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3) + + # add title to figure + fig.suptitle(title, fontsize=fs1) + + # save figure + figname='rank_hist_sumstats.png' + plt.savefig(figname, dpi=300) + +def read_data(rh_summary_file): + ''' + Input: path/to/rh_summary_file + Output: p_freq, q_freq + Description: reads in p frequencies and q frequencies + ''' + + # first line contains a header and first column contains hour + p_freq_all = np.loadtxt(rh_summary_file + '_p-freq.txt', skiprows=1) + q_freq_all = np.loadtxt(rh_summary_file + '_q-freq.txt', skiprows=1) + + # p_freq and q_freq contain data 00, 06, 12, 18, 24 (the last is used in these plots) + p_freq = p_freq_all[-1,1:] + q_freq = q_freq_all[-1,1:] + + return p_freq, q_freq + +if __name__=='__main__': + rh_summary_file = sys.argv[1] + number_of_stations = sys.argv[2] + p01 = sys.argv[3] + p1 = sys.argv[4] + p5 = sys.argv[5] + max_freq_p = float(sys.argv[6]) + max_freq_q = float(sys.argv[7]) + + p_freq, q_freq = read_data(rh_summary_file) + plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q) + diff --git a/SATELLITE/STATST/summary_rank_histograms_all_locations.sh b/SATELLITE/STATST/summary_rank_histograms_all_locations.sh new file mode 100755 index 0000000000000000000000000000000000000000..d3d2d191deb7fceef3879e6f3d650a56b59ddb7e --- /dev/null +++ b/SATELLITE/STATST/summary_rank_histograms_all_locations.sh @@ -0,0 +1,250 @@ +#!/bin/bash +echo "Bash version ${BASH_VERSION}..." +echo "Python version" +which python +################################################################## +# +# Calculation and plotting of summary statistics for quantile +# space rank histograms: +# +# - a map of p-values +# - bar plot of p-value distribution + rank histogram +# averaged over all stations +# +# The following need to be given as arguments: +# +# - variable +# - first year of raw data +# - last year of raw data +# - first month of raw data (default 1) +# - last month of raw data (default 12) +# +# As the input for this script, the following files are needed: +# +# 1) Rank histogram statistics files for all available stations +# (for the time range specified by the script arguments) +# +# Execution (e.g): ./summary_rank_histograms_all_stations 2010 2010 1 12 + +# ORIGINAL: +# Jouni Räisänen, Aug 2023 +# MODIFIED: +# Alexander Mahura, Sep-Oct-Nov 2023 +# Lauri Tuppi, Jan-Feb-Mar 2024 +# Alexander Mahura, Mar-Apr 2024 +# +################################################################## +# +# Arguments: +# +# 0. Variable code (not needed: always T2-T850) +# +variable=$1 # btmp for AMSU-A Channel 5 +echo " Variable number: 194" +echo " Variale name: btmp" +echo " Variable AMSU-A channel: $variable" + +# 1.-2: First and last year +# +year1=$2 +year2=$3 +let nyears=year2-year1+1 +# +# 3.-4: First and last month +# +month1="${4:-1}" +month2="${5:-12}" +# +# Add zero to ahead of $month1 and $month2 for naming of files if needed? +# +if [ $month1 -le 9 ] +then + month1='0'${month1} +fi +if [ $month2 -le 9 ] +then + month2='0'${month2} +fi +# +################################################################## +# +# Add odb_api to $PATH +# +# On Puhti +#PATH="${PATH}":/projappl/project_2001011/odb_api/bin +# On Lumi +#export PATH="/projappl/project_465000454/ama/software/odb_api/bin:$PATH" + +# On Lumi +#PATH="${PATH}":/projappl/project_465000454/ama/software/odb_api/bin +# +# On Lumi (asked Devaraju Narayanappa, CSC) +# to make in /project/project_465000454/devaraju/modules/LUMI/23.03/C +# module use /project/project_465000454/devaraju/modules/LUMI/23.03/C +# module load odb_api/0.18.1-cpeCray-23.03.lua +# module load python-climatedt/3.11.3-cpeCray-23.03.lua +# +################################################################## +# +echo " Compiling the Fortran program that produces the rank histogram summary statistics ..." +echo " fortran-programs/rank_histogram_sumstats_locations.f95" +# +gfortran fortran-programs/rank_histogram_sumstats_locations.f95 -o rank_histogram_sumstats_locations.exe + +################################################################## +# ----- File names. Hard-coded, at least this far ----- +# +# NB: it would be more efficient to have ODB files for individual +# stations instead for the ones that include everything for all stations! +# +# 1) Rank histogram directory and input name file name without extension +echo " Directory for rank histograms ..." +#rh_dir=/scratch/project_2001011/RHARM_to_Alexander_120923/rank_histograms/${variable}_${year1}${month1}-${year2}${month2} +rh_dir=rank_histograms/${variable}_${year1}${month1}-${year2}${month2} +echo " $rh_dir" +echo " Input file with rank histograms for all stations ..." +rh_file=${rh_dir}/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_all_locations +echo " $rh_file" +echo " Input file with sammy rank histograms ..." +#rh_summary_file=${rh_dir}/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_summary +rh_summary_file=tmp_output/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_summary +echo " $rh_summary_file" + +# 2) Name of output file(s) without extension +echo " Name of output file without extension ..." +out_file=${rh_dir}/rh_summary_${variable}_${year1}${month1}-${year2}${month2} +echo " $out_file" + +# 3) Directory for figures +echo " Directory for figures ..." +figure_dir=figures +echo " $figure_dir" + +################################################################## +echo " Converting the all-station ODB format rank histogram file to txt ..." +echo " (it is format for reading in Fortran)" +# +#odb sql select \* -i ${rh_file}.odb -o ${rh_file}.txt +odb sql select \* -i ${rh_file}.odb -o tmp_output/rank_histograms_${variable}_${year1}${month1}-${year2}${month2}_all_locations.txt +# +################################################################## +echo " Calculating the rank histogram summary statistics ..." +# +./rank_histogram_sumstats_locations.exe < line_1 +read v1 v2 v3 number_of_stations < line_1 +echo ' number_of_stations:' $number_of_stations + +head -2 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 number_of_p_bins < line_1 +echo ' number_of_p_bins:' $number_of_p_bins + +head -3 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 max_freq_p < line_1 +#echo 'max_freq_p:' $max_freq_p +head -4 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 max_freq_q < line_1 +#echo 'max_freq_q:' $max_freq_q +echo " max_freq_p: $max_freq_p ; max_freq_q: $max_freq_q" + +head -5 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p01_00 < line_1 +#echo 'p01_00:' $p01_00 +head -6 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p1_00 < line_1 +#echo 'p1_00:' $p1_00 +head -7 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p5_00 < line_1 +#echo 'p5_00:' $p5_00 +echo " p01_00: $p01_00 ; p1_00: $p1_00 ; p5_00: $p5_00" + +head -8 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p01_12 < line_1 +#echo 'p01_12:' $p01_12 +head -9 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p1_12 < line_1 +#echo 'p1_12:' $p1_12 +head -10 ${rh_summary_file}_text_output | tail -1 > line_1 +read v1 v2 v3 v4 v5 v6 p5_12 < line_1 +#echo 'p5_12:' $p5_12 +echo " p01_12: $p01_12 ; p1_12: $p1_12 ; p5_12: $p5_12" + +numbers_of_p_bins_up=$( echo "(${number_of_p_bins} + 0.5)" | bc ) +echo $( echo "(${number_of_p_bins} + 0.5)" | bc ) > line_1 +read number_of_p_bins_up < line_1 +rm line_1 + +###################################################################### +# +echo "====================================================================" +echo " (1) Plotting - P-Values Summary on 2D map for Daily Mean btmp ..." +# +###################################################################### + +python_script=python/plot_p_values_map_amsua.py + +figname=${figure_dir}/p-values_map_amsua_${variable}_${year1}${month1}-${year2}${month2}.png +echo " Output file: $figname" + +python3 ${python_script} ${rh_summary_file}_p-values.txt ${year1}${month1}-${year2}${month2} ${figname} + +##################################################################### +# +echo " (2) Plotting - Rank Histogram Summary Statistics for Daily Mean btmp ..." +# +##################################################################### + +utc_term="Daily Mean" +echo " $utc_term - btmp (194) - AMSU-A channel : ${variable}" + +python_script=python/plot_rankhist_sum_all_locations_amsua.py + +echo " number_of_stations : ${number_of_stations}" +echo " p01_12, p1_12, p5_12, max_freq_p, max_freq_q : ${p01_12} ${p1_12} ${p5_12} ${max_freq_p} ${max_freq_q}" + +python3 ${python_script} ${rh_summary_file} ${number_of_stations} ${p01_12} ${p1_12} ${p5_12} ${max_freq_p} ${max_freq_q} ${utc_term} ${figname} + +mv rank_hist_sumstats.png ${figure_dir}/rankhist_sumstats_amsua_${variable}_${year1}${month1}-${year2}${month2}.png +echo " Output file for Daily Mean: ${figure_dir}/rankhist_sumstats_amsua_${variable}_${year1}${month1}-${year2}${month2}.png" + +###################################################################### +# Delete files that are not needed any more +rm -rf tmp_output +rm -rf *.exe + +if [ 0 -eq 1 ]; then + rm rank_histogram_sumstats_locations.exe + rm tmp_output/coordinates + rm tmp_output/input.txt + rm tmp_output/msd_and_p-value + rm plots_for_one_location.exe + rm tmp_output/sim_data + rm tmp_output/sim_data_hourly + rm tmp_output/time_series_commands + rm tmp_output/vrange_commands + rm rank_histograms_bootstrap.exe + rm rank_histograms_one_location.exe + rm tmp_output/quantiles* + rm tmp_output/rank_histogram_* + rm tmp_output/rank_histograms_* + rm tmp_output/standard_plot.cfg + rm tmp_output/time_series* +fi diff --git a/SATELLITE/STATST/test_stats.bash b/SATELLITE/STATST/test_stats.bash new file mode 100755 index 0000000000000000000000000000000000000000..7cc7a597245428b84091c95d94eabaf98c3e705f --- /dev/null +++ b/SATELLITE/STATST/test_stats.bash @@ -0,0 +1,47 @@ +#!/bin/bash + +# This script is supposed to emulate running ../amsua_stats.sh $channel $b_yy $b_yy $b_mm_start $b_mm_end for one time. + +source ../set_env.sh a + +channel=5 +b_yy_start=2000 +b_yy_end=2000 +b_mm_start=01 +b_mm_end=12 + +echo "channel, b_yy_start, b_yy_end, b_mm_start, b_mm_end : $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end" + +echo "*****************************************************************" +echo " STATS - PRODUCING QUANTILE RANK HISTOGRAM STATISTICS AND PLOTS" +echo " ON EXAMPLE OF AMSU-A BRIGHTNESS TEMPERATURE: channel 5" +echo " AT SELECTED LOCATIONS" +echo "*****************************************************************" + +#cd STATST +pwd + +echo "*****************************************************************" +echo " STEP STATST-1 - PRODUCE RANK HISTOGRAMS FOR THE SELECTED LOCATIONS" +echo "*****************************************************************" + +./produce_rank_histograms_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "*****************************************************************" +echo " STEP STATST-2 - PRODUCE STANDARD PLOTS FOR EACH LOCATION" +echo "*****************************************************************" + +./produce_standard_plots_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "*****************************************************************" +echo " STEP STATST-3 - PRODUCE SUMMARY RANK HISTOGRAMS FOR ALL LOCATIONS" +echo "*****************************************************************" + +./summary_rank_histograms_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "*****************************************************************" +echo " ALL STATST-# STEP COMPLETED" +echo "*****************************************************************" + +exit + diff --git a/SATELLITE/amsua_mod.sh b/SATELLITE/amsua_mod.sh new file mode 100755 index 0000000000000000000000000000000000000000..5d865aeb55d5b9398a276d00d6b72585833a2bff --- /dev/null +++ b/SATELLITE/amsua_mod.sh @@ -0,0 +1,131 @@ +#!/bin/bash + +# See main_amsua.sh for construction of command line arguments. +cdate=${1?} #YYYYMMDDHH +midpoint_time=${2?} #HHMMSS + +echo "Entering amsua_mod.sh of AMSU-A workflow." +#source set_env.sh + +# EUMETSAT NWP-SAF Radiance Simulator satellite observation operator control block. +# This part writes running instructions for Radiance Simulator. + +cdate_s2=$(echo $cdate | cut -c 1-8) + +if [ $process_mode == 'rf' ]; then + obs_datafile=$amsua_input/amsua_obs_$cdate_s2$midpoint_time.dat + model_datafile=$gsv_depo_processed/amsu-a_input_${cdate_s2}_limited.grb +elif [ $process_mode == 'cli' ]; then + obs_datafile=$scriptdir/other_files/amsua_locations_sample.dat + model_datafile=$gsv_depo/$cdate/amsu-a_input_${cdate_s2}_limited.grb +fi + +# This is for RadSim 3.2 +cat > $amsua_input/radsim_cfg_$cdate.nl << EOF +&radsim_nl +obs_datafile = '$obs_datafile' +model_datafile = '$model_datafile' +model_filetype = 1 +rttov_coeffs_dir = '$rttov_coefs/rttov13pred54L' +rttov_hydrotable_dir = '$rttov_coefs/hydrotable' +rttov_sccld_dir = '$rttov_coefs/cldaer_ir' +rttov_coeffs_type = '.dat' +rttov_coeffs_options = '' !'_o3co2' +emiss_atlas_dir = '$rttov_coefs/emis_data' +brdf_atlas_dir = '$rttov_coefs/brdf_data' +platform = 'metop' +inst = 'amsua' +satid = 2 +channels = 5 +nthreads = 1 +temporal_data = .false. +enable_footprints = .false. +run_scatt = .false. +output_mode = 1 +seaice_threshold = 0.2 +ozone_data = .false. +clw_data = .false. +mw_clw_scheme = 2 +addsolar = .false. +ir_addclouds = .false. +ircloud_ice_scheme = 1 +fastem_version = 6 +do_lambertian = .false. +do_nlte_correction = .false. +output_dir = '$amsua_output' +output_file = 'mod_amsua_$cdate.nc' +write_radiances = .false. +write_profiles = .false. +write_latlon = .false. +write_emiss = .false. +write_brdf = .false. +write_trans = .false. +write_tjac = .false. +write_qjac = .false. +write_o3jac = .false. +write_tskinjac = .false. +write_wind10mjac = .false. +write_emissjac = .false. +/ +EOF +#source $radsim_run $amsua_input/radsim_cfg_$cdate.nl +export LD_LIBRARY_PATH=/opt/cray/pe/netcdf/4.9.0.7/gnu/9.1/lib:/lib:/scratch/project_465000454/devaraju/SW/gcc/eccodes-2.32.0/lib:/lib:/lib:$LD_LIBRARY_PATH +ulimit -s unlimited 2>/dev/null +radsim.exe $amsua_input/radsim_cfg_$cdate.nl + +############################################################################### +# Converting amsu-a mod data file into odb (with import command) +# Opt1: saving/adding records to ODB (for mod data at synop stations locations) +# Opt2: saving adding records to ODB (for odb already having obs synop data) +# +# modify-update (& unify) according to different types of obs data: +# SYNOP (only at ground), +# RADSOUND (at 3D - lat/lon/alt) +# SATELLITE (at 2D - area within lat/lon) + +mkdir $amsua_output/tmp_$cdate +echo " Converting/adding/saving records with amsu-a mod data into ODB ..." + +amsua_output_txt1=$amsua_output/tmp_$cdate/tmp_file1_$cdate.txt +amsua_output_txt2=$amsua_output/tmp_$cdate/tmp_file2_$cdate.txt +python3 convert_netcdf2txt.py $amsua_output/mod_amsua_$cdate.nc $amsua_output_txt1 +rm $amsua_output/mod_amsua_$cdate.nc + +b_yy=$( echo $cdate | cut -c 1-4 ) +b_mm=$( echo $cdate | cut -c 5-6 ) +b_dd=$( echo $cdate | cut -c 7-8 ) +b_hh=$( echo $cdate | cut -c 9-10 ) +b_min="00" +b_sec="00" + +# .............................................................................. +# AMSU-A data ................................................................... +# Defining header for adding mod amsua data to odb file ... +#header_mod_amsua_data="year@descr:integer month@descr:integer day@descr:integer hour@descr:integer minute@descr:integer second@descr:integer longitude@descr:real latitude@descr:real value@descr:real variable@descr:integer station@descr:string" +if [ $process_mode == 'rf' ]; then + header_mod_amsua_data="year@descr:integer month@descr:integer day@descr:integer hour@descr:integer minute@descr:integer second@descr:integer longitude@descr:real latitude@descr:real solar_zenith@descr:real zenith@descr:real lsm@descr:real obsvalue@descr:real modvalue@descr:real scanline@descr:real scanpos@descr:real satellite_identifier@descr:integer" +elif [[ $process_mode == 'cli' || $process_mode == 'dummy' ]]; then + header_mod_amsua_data="year@descr:integer month@descr:integer day@descr:integer hour@descr:integer minute@descr:integer second@descr:integer longitude@descr:real latitude@descr:real zenith@descr:real modvalue@descr:real" +fi + +echo $header_mod_amsua_data >> $amsua_output_txt2 + +nlines=$(wc -l < $amsua_output_txt1) +for (( line=1; line<=$nlines; line++ )); do + amsua_record=$( sed -n ${line}p $amsua_output_txt1 ) + echo $b_yy ' ' $b_mm ' ' $b_dd ' ' $b_hh ' ' $b_min ' ' $b_sec ' ' $amsua_record >> $amsua_output_txt2 +done + +tmp_mod_amsua_odb=$amsua_output/tmp_$cdate/tmp_mod_amsua_odb.odb +# Adding data to odb file ... + +echo " Path and Name of odb-file with saved records:" +echo " $mod_amsua_odb" +odb import -d ' ' $amsua_output_txt2 $tmp_mod_amsua_odb +cat $tmp_mod_amsua_odb >> $mod_amsua_odb + +rm -rf $amsua_output/tmp_$cdate +rm -f $amsua_input/radsim_cfg_$cdate.nl +rm -f $amsua_input/amsua_obs_$cdate_s2$midpoint_time.dat + +echo "Exiting amsua_mod.sh of AMSU-A workflow." diff --git a/SATELLITE/amsua_obs.sh b/SATELLITE/amsua_obs.sh new file mode 100755 index 0000000000000000000000000000000000000000..5498cddf1ac09b6b8e654c4a2152d9c4456bdb58 --- /dev/null +++ b/SATELLITE/amsua_obs.sh @@ -0,0 +1,57 @@ +#!/bin/bash + +# See main_amsua.sh for construction of command line arguments. +cdate=${1?} #YYYYMMDDHH +midpoint_time=${2?} #HHMMSS + +echo "Entering amsua_obs.sh of AMSU-A workflow." +source set_env.sh b + +# Selection of observation metadata (latitude, longitude, ...) +# In the case of reality-following simulation it is possible to make one-to-one comparison with actual observations, so select AMSU-A metadata matching the area selected above. +# In the case of free-running climate simulation, just use a file containing manually set points. + +#if [ 1 -eq 1 ]; then +#start_minutes=$(hhmmss_to_minutes "$start_time") +#end_minutes=$((start_minutes + time_range)) +#end_time=$(minutes_to_hhmmss "$end_minutes") +#mid_minutes=$((start_minutes + time_range / 2)) +#midpoint_time=$(minutes_to_hhmmss "$mid_minutes") +mid_minutes=$(hhmmss_to_minutes "$midpoint_time") +start_minutes=$((mid_minutes - time_range / 2)) +end_minutes=$((mid_minutes + time_range / 2)) + +start_time=$(minutes_to_hhmmss "$start_minutes") +#echo $start_minutes $start_time '#############################' +end_time=$(minutes_to_hhmmss "$end_minutes") + +time_step=-1 +cdate_long1=$(compute_end_date "$cdate" "$time_step") +cdate_s1=$(echo $cdate_long1 | cut -c 1-8) +cdate_s2=$(echo $cdate | cut -c 1-8) + +#cdate1=$cdate # decided to give the dates as command line arguments +#cdate_s1=$(echo $cdate1 | cut -c 1-8) +#cdate2=`exec $calendar $cdate1 + 1` +#cdate_s2=$(echo $cdate2 | cut -c 1-8) + +input=$amsua_obs_database/$cdate_s1.dat #20161130.dat +#outf=$output/filtered_metadata/amsua_obs_$cdate_s2$midpoint_time.dat +outf=$amsua_input/amsua_obs_$cdate_s2$midpoint_time.dat +#amsua_obs_20161202000000.dat +if [[ $process_mode == 'rf' || $process_mode == 'dummy' ]]; then + echo $input + echo $start_time + echo $time_range + echo $outf + python3 filter_data.py $input $start_time $time_range $outf + if [ 1 -eq 1 ]; then # prepend info for RadSim. + nobs=$(wc -l < $outf) + echo -e "2\n$nobs\n5\n1 2 15 5 4 7 8\n$(cat $outf)" > $outf + fi + start_time=$end_time +elif [ $process_mode == 'cli' ]; then + test -f $outf || cp $scriptdir/other_files/amsua_locations_sample.dat $outf +fi + +echo "Exiting amsua_obs.sh of AMSU-A workflow." diff --git a/SATELLITE/amsua_stats.sh b/SATELLITE/amsua_stats.sh new file mode 100755 index 0000000000000000000000000000000000000000..4da20ccf4f766b651f8ee549df6be45690b0c776 --- /dev/null +++ b/SATELLITE/amsua_stats.sh @@ -0,0 +1,74 @@ +#!/bin/bash +echo "Bash version ${BASH_VERSION}..." +######################################################################## +# Author: Alexander Mahura : 2023-08-22 +# updated by Lauti Tuppi : 2024-02/03-... +# updated by Alexander Mahura : 2024-03/04-... +######################################################################## + +# PRE-SETUP (TESTING FOR LIMITED SATELLITE ODB DATASET & LIMITED GEOGRAPHICAL AREA) +# FOR BRIGHTNESS TEMPERATURE (btmp) +# OTHER METEOROLOGICAL PARAMETER (FROM SATELLITE AVAILABLE CAN BE LISTED) CAN BE ADDED +# +# Producing quantile rank histogram statistics and plots (using limited open odb datasets) +# +# INPUT: +# start & end year, start & end month, variable name and number, channel numberfor AMSU-A + +b_yy_start=$1 +b_yy_end=$2 +b_mm_start=$3 +b_mm_end=$4 +varMETname=$5 +varMETnum=$6 +channel=$7 + +#channel=5 +#b_yy_start=2000 +#b_yy_end=2000 +#b_mm_start=01 +#b_mm_end=12 + +echo " b_yy_start, b_yy_end, b_mm_start, b_mm_end : $b_yy_start $b_yy_end $b_mm_start $b_mm_end" +echo " varMETname, varMETnum, channel : $varMETname $varMETnum $channel" + +echo "*****************************************************************" +echo " STATST - PRODUCING QUANTILE RANK HISTOGRAM STATISTICS AND PLOTS" +echo " ON EXAMPLE OF BRIGHTNESS TEMPERATURE: AMSU-A Channel 5" +echo " AT SELECTED LOCATIONS/ AREAS" +echo "*****************************************************************" + +cd STATST +pwd +#t exit + +echo "***************************************************************************" +echo " STEP STATST-1 - PRODUCE RANK HISTOGRAMS FOR THE SELECTED LOCATIONS/ AREAS " +echo "***************************************************************************" + +#./produce_rank_histograms_all_stations.sh 39 2020 2020 1 12 +#./produce_rank_histograms_all_locations.sh $varMETnum $b_yy_start $b_yy_end $b_mm_start $b_mm_end +./produce_rank_histograms_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "****************************************************************" +echo " STEP STATST-2 - PRODUCE STANDARD PLOTS FOR EACH LOCATION/ AREA " +echo "****************************************************************" + +#./produce_standard_plots_all_stations.sh 39 2020 2020 1 12 +#./produce_standard_plots_all_locations.sh $varMETnum $b_yy_start $b_yy_end $b_mm_start $b_mm_end +./produce_standard_plots_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "**************************************************************************" +echo " STEP STATST-3 - PRODUCE SUMMARY RANK HISTOGRAMS FOR ALL LOCATIONS/ AREAS " +echo "**************************************************************************" + +#./summary_rank_histograms_all_stations.sh 39 2020 2020 1 12 +#./summary_rank_histograms_all_locations.sh $varMETnum $b_yy_start $b_yy_end $b_mm_start $b_mm_end +./summary_rank_histograms_all_locations.sh $channel $b_yy_start $b_yy_end $b_mm_start $b_mm_end + +echo "*****************************************************************" +echo " ALL STATST-# STEP COMPLETED" +echo "*****************************************************************" + +exit + diff --git a/SATELLITE/convert_netcdf2txt.py b/SATELLITE/convert_netcdf2txt.py new file mode 100644 index 0000000000000000000000000000000000000000..7086f25ea8b7c459a62e93dc4dcfec7628ed6e1d --- /dev/null +++ b/SATELLITE/convert_netcdf2txt.py @@ -0,0 +1,31 @@ +import netCDF4 +from netCDF4 import Dataset +import numpy as np +import sys + +fname=sys.argv[1] #'/scratch/project_2001011/Lauri/AMSUA_demo_wf/amsua_output/mod_amsua_2016120102.nc' +outf=sys.argv[2] + +infile=Dataset(fname,'r') + +bt=infile.variables['bt'][:] +lat=infile.variables['lat'][:] +lon=infile.variables['lon'][:] +satzen=infile.variables['satzen'][:] +solzen=infile.variables['solzen'][:] +lsm=infile.variables['lsm'][:] + +shape=np.shape(bt) +ncols=3+(shape[0]*shape[1]) +data_array=np.zeros((shape[1],ncols)) + +for i in range(shape[1]): + data_array[i,0]=lon[i] + data_array[i,1]=lat[i] + data_array[i,2]=satzen[i] + for j in range(shape[0]): + data_array[i,j+3]=bt[j,i] + +#print(data_array) + +np.savetxt(outf,data_array) diff --git a/SATELLITE/convert_pl2ml.py b/SATELLITE/convert_pl2ml.py new file mode 100755 index 0000000000000000000000000000000000000000..19dccc3a1f2f7f8fafa06adf29e497086831d336 --- /dev/null +++ b/SATELLITE/convert_pl2ml.py @@ -0,0 +1,199 @@ +#!/usr/bin/env python + +''' +Written by Madeleine Ekblom, June 2023 + +Description: +Convert data on pressure levels to model levels and write model-level data to netcdf + +Set python environment (contains e.g., metview and metview-python): +export PATH="/projappl/project_2001011/madeleine/python_envs/pymet/bin:$PATH" + +Example: +python3 convert_oifs_pl2mpl.py input_pl input_surf vct_dat + +''' + +import sys +import numpy as np +#import cfgrib +import xarray as xr +import pandas as pd +import scipy.interpolate as intp +#import metview as mv +from netCDF4 import Dataset + +def interpolate_pl2ml(input_pl, input_surf, vct_dat, output_file, lat_S, lat_N, lon_W, lon_E): + ''' + Reads grib files, interpolates data to model levels, writes interpolated data to file + Input: input_ml: grib file on regular grid and pressure levels (model levels and surface levels), vct_dat: spherical harmonics, output_file (/path/to/netcdf/file) lat_S, lat_N, lon_W, lon_E + Output: writes to file interpolated data as netcdf + ''' + + # Read in data on pressure levels: + #pl_data = mv.read(input_pl) + pl_data = Dataset(input_pl, mode='r') + + # Read in data on surface levels: + #surf_data = mv.read(input_surf) + surf_data = Dataset(input_surf, mode='r') + + # Convert Fieldsets to xarray Dataset: + #pl_xa = pl_data.to_dataset() + #surf_xa = surf_data.to_dataset() + + # Get pressure levels, latitude, and longitude + #p = pl_xa.coords['isobaricInhPa'] + #lats = pl_xa.coords['latitude'].values + #lons = pl_xa.coords['longitude'].values + #time = pl_xa.coords['time'].values + #p = pl_data.variables['plev'][:] + p = pl_data.variables['level'][:] + lats = pl_data.variables['lat'][:] + lons = pl_data.variables['lon'][:] + time = pl_data.variables['time'][:] + #print('#################',time_t) + + # Get specific humidity, cloud liquid water content, temperature + q0 = pl_data.variables['q'] + q = q0[0,:,:,:] + t0 = pl_data.variables['t'] + t = t0[0,:,:,:] + #print(t[:,50,50]) # OK + clwc0 = pl_data.variables['clwc'] + clwc = clwc0[0,:,:,:] +# q = pl_data.variables['Q'][:].squeeze() +# t = pl_data.variables['T'][:].squeeze() +# clwc = pl_data.variables['CLWC'][:].squeeze() + + # Get logarithmic surface pressure and surface pressure + #lnsp = surf_xa.data_vars['lnsp'] + sp = surf_data.variables['sp'] # (Pa) +# sp = surf_data.variables['SP'][:] + + + # Read in vct for model level conversion + vct_dat = np.loadtxt(vct_dat, skiprows=1, dtype={'names':('k', 'vct_a','vct_b'), 'formats':('i4', 'float64', 'float64')}) + vct_a = vct_dat['vct_a'] + vct_b = vct_dat['vct_b'] + k = vct_dat['k'] + nlevs = np.max(k) + + # Create list of model levels: + #time = np.ones(1) + mlevs = np.arange(1,nlevs+1) + nlats = len(lats) + nlons = len(lons) + + ### Vertical interpolation ### + # Calculate pressure on half-levels: p_half(levs, lats,lons) = a(levs) + b(levs)*sp(lats, lons) + p_half = np.zeros((nlevs+1, nlats, nlons)) + check = sp + #print(check[:,50]) + for n in range(nlevs+1): + p_half[n,:,:] = (vct_a[n] + vct_b[n]*sp)/100 # (sp in Pa -> divide by 100 to get hPa) + + # Calculate pressure on model levels: p_model = ( p_half_above + p_half below )/2 + p_ml_temp = (p_half[:-1,:,:] + p_half[1:,:,:])/2 + p_ml = np.zeros((1, nlevs, nlats, nlons)) + p_ml[0,:,:,:]=p_ml_temp + + # Use regular grid interpolator to form interpolation function: + print(np.shape(q)) + q_interp = intp.RegularGridInterpolator((p,lats,lons), q, bounds_error=False, fill_value=None) + clwc_interp = intp.RegularGridInterpolator((p,lats,lons), clwc, bounds_error=False, fill_value=None) + t_interp = intp.RegularGridInterpolator((p,lats,lons), t, bounds_error=False, fill_value=None) + + # Form new grid points: p_ml x lats x lons + # lats = nlats x 1 + # lons = nlons x 1 + # p_ml = nlevs x nlats x nlons + + # Find closest indices corresponding to lat_south, lat_north, lon_west, lon_east + lat_S_ind = np.argmin(np.abs(lats - lat_S)) + lat_N_ind = np.argmin(np.abs(lats - lat_N)) + lon_W_ind = np.argmin(np.abs(lons - lon_W)) + lon_E_ind = np.argmin(np.abs(lons - lon_E)) + + # Create zero numpy arrays for saving interpolated data + q_intp = np.zeros((1, nlevs, nlats, nlons)) + clwc_intp = np.zeros((1, nlevs, nlats, nlons)) + t_intp = np.zeros((1, nlevs, nlats, nlons)) + + # typeOfLevel=hybrid + c1 = 1 # counter + + # only loop over restricted area from west to east and north to south: (lon_W, lat_N) --> (lon_E, lat_S) + if (lat_N_ind>lat_S_ind): + aaa=lat_N_ind + bbb=lat_S_ind + lat_N_ind=bbb + lat_S_ind=aaa + for lon_ind in range(lon_W_ind, lon_E_ind): # --> nlons, if looping over all grid points + #for lon_ind in range(nlons): + c2=1 + for lat_ind in range(lat_N_ind, lat_S_ind): # --> nlats, if looping over all grid points + #for lat_ind in range(nlats): + col = np.array((p_ml[0,:,lat_ind,lon_ind], lats[lat_ind]*np.ones(nlevs), lons[lon_ind]*np.ones(nlevs))).T + for lev_ind in range(nlevs): + res_q = q_interp(col[lev_ind,:]) + res_clwc = clwc_interp(col[lev_ind,:]) + res_t = t_interp(col[lev_ind,:]) + q_intp[0,lev_ind,lat_ind,lon_ind] = res_q + clwc_intp[0,lev_ind,lat_ind,lon_ind] = res_clwc + t_intp[0,lev_ind,lat_ind,lon_ind] = res_t + print('lon: ',c1,' lat: ',c2) + c2=c2+1 + c1=c1+1 + # Create xarray Dataset for interpolated data: + q_xr = xr.DataArray(q_intp, dims=("time", "level", "latitude", "longitude"), coords={"time": time, "level": mlevs, "latitude": lats, "longitude": lons}) + clwc_xr = xr.DataArray(clwc_intp, dims=("time", "level", "latitude", "longitude"), coords={"time": time, "level": mlevs, "latitude":lats, "longitude":lons}) + t_xr = xr.DataArray(t_intp, dims=("time", "level", "latitude", "longitude"), coords={"time": time, "level": mlevs, "latitude": lats, "longitude": lons}) + p_xr = xr.DataArray(p_ml, dims=("time", "level", "latitude", "longitude"), coords={"time": time, "level": mlevs, "latitude": lats, "longitude": lons}) + + data_intp = xr.Dataset(dict(q=q_xr, clwc=clwc_xr, t=t_xr, pres=p_xr)) + + data_intp.attrs['gridType'] = 'regular_ll' + data_intp.attrs['typeOfLevel'] = 'hybrid' + + # Metadata need to be added to dataset/provided to the function before/when converting to fieldset (how?) + + # key/value + typeOfGrid='regular_ll' + typeOfLevel='hybrid' + paramID=[133, 246, 130, 54] + units=[] + shortName=['q', 'clwc', 't', 'pres'] + name=[] + + # Save interpolated to netcdf: + data_intp.to_netcdf(output_file) + +def main(input_pl, input_surf, vct_dat, output_file, lat_S, lat_N, lon_W, lon_E): + ''' + Input: path/to/grib file: pressure level inputs on regular grid, path/to/grib file: surface level inputs on regular grid, vct_dat containing information related to model levels (A, B) + Calls interpolation from pressure to model level for a given grid box (lat_N, lat_S, lon_W, lonE) + ''' + + interpolate_pl2ml(input_pl, input_surf, vct_dat, output_file, lat_S, lat_N, lon_W, lon_E) + + +if __name__=='__main__': + # check number of input arguments: + if len(sys.argv) < 9: + sys.exit("Error: path/to/grib/file (regular grid, pressure levels) path/to/grib/file (regular grid, surface levels), vct_data and path/to/output/netcdf must be passed as command line arguments" +) + elif len(sys.argv) > 9: + sys.exit("Error: too many command line arguments passed. Only four required: path/to/grib/file (regular grid, pressure levels), path/to/grib/file (regular grid, surface levels), vct_data, and /path/to/output/netfcdf ") + + input_pl = sys.argv[1] + input_surf = sys.argv[2] + vct_data = sys.argv[3] + output_file = sys.argv[4] + + lat_S = float(sys.argv[5]) + lat_N = float(sys.argv[6]) + lon_W = float(sys.argv[7]) + lon_E = float(sys.argv[8]) + + main(input_pl, input_surf, vct_data, output_file, lat_S, lat_N, lon_W, lon_E) diff --git a/SATELLITE/filter_data.py b/SATELLITE/filter_data.py new file mode 100755 index 0000000000000000000000000000000000000000..913978d3a818c84c5376df64f6c167562df6eba8 --- /dev/null +++ b/SATELLITE/filter_data.py @@ -0,0 +1,83 @@ +import sys +from datetime import datetime, timedelta +# A more advanced version of AMSU-A observation filtering. +# This version is also capable of handling cases around the midnight. + +def parse_time(indate, time_str): + # This makes sure that the time is always HHMMSS. + # Then the time string is converted into datetime object. + if len(time_str)==1: + time_str = '00000' + time_str + elif len(time_str)==2: + time_str = '0000' + time_str + elif len(time_str)==3: + time_str = '000' + time_str + elif len(time_str)==4: + time_str = '00' + time_str + elif len(time_str)==5: + time_str = '0' + time_str + elif len(time_str)>6: + print("Error! The time_str cannot have more than 6 characters! HHMMSS, please.") + sys.exit(1) + date_time_str = indate + ' ' + time_str + return datetime.strptime(date_time_str, "%Y%m%d %H%M%S") + +def filter_data(input_file, start_time, end_time, output_file, mode): + # This is the actual data selection routine. + # All obs between start time and end time are written into a new file. + lat_S = 69.5 + lat_N = 70.5 + lon_W = 43.5 + lon_E = 46.5 + write_obs = "T" + indate=input_file[-12:-4] + with open(input_file, "r") as f_in, open(output_file, mode) as f_out: + counter = 0 + for line in f_in: + columns = line.split() + time = parse_time(indate,columns[1]) +# if start_time <= time < end_time: +# counter = counter + 1 +# line_out = columns[2]+' '+columns[3]+' '+columns[4]+' '+columns[5]+' '+columns[6]+' '+columns[8]+' '+columns[9]+'\n' +# f_out.write(line_out) + if start_time <= time < end_time and lon_W <= float(columns[2]) < lon_E and lat_S <= float(columns[3]) < lat_N: + counter = counter + 1 + if write_obs == "T": + line_out = columns[1]+' '+ columns[2]+' '+columns[3]+' '+columns[7]+'\n' + elif write_obs == "F": + line_out = columns[2]+' '+columns[3]+' '+columns[4]+' '+columns[5]+' '+columns[6]+' '+columns[8]+' '+columns[9]+'\n' + f_out.write(line_out) + + print("Number of selected observations: ", counter) + print("Filtered data has been written to", output_file) + +if __name__ == "__main__": + if len(sys.argv) != 5: + print("Usage: python script.py ") + sys.exit(1) + + # Read the command line arguments. + input_file = sys.argv[1] + start_time_str = sys.argv[2] + time_range = int(sys.argv[3]) + output_file = sys.argv[4] + + # Get the start date from the name of the file. + # Also extract the input directory from the full path. + start_date_str=input_file[-12:-4] + indir=input_file[:-12] + start_time = parse_time(start_date_str, start_time_str) + end_time = start_time + timedelta(minutes=time_range) + mode="w" # "w" for writing + + filter_data(input_file, start_time, end_time, output_file, mode) + + end_date_str = end_time.strftime("%Y%m%d") + # This is a special treatment for obs windows extending over midnight. + if start_date_str != end_date_str: + print('Discontinuity of time at midnight.') + input_file = indir + end_date_str + '.dat' + print('-> Collecting the remaining observations from the file corresponding to the next day (', input_file, ') and appending them to the output.') + mode = "a" # "a" for appending + filter_data(input_file, start_time, end_time, output_file, mode) + diff --git a/SATELLITE/gsv_mod_data.sh b/SATELLITE/gsv_mod_data.sh new file mode 100755 index 0000000000000000000000000000000000000000..c84bd577d8763845257eee88d930a31e42b9f69a --- /dev/null +++ b/SATELLITE/gsv_mod_data.sh @@ -0,0 +1,170 @@ +#!/bin/bash + +echo "Entering gsv_mod_data.sh of AMSU-A workflow." +cdate_s2=${1?} +time_s=${2?} + +source set_env.sh b + +hour=$(echo $time_s | cut -c 1-2) + +test -d $gsv_depo/$cdate_s2$hour || mkdir $gsv_depo/$cdate_s2$hour +yaml_file_pl=$gsv_depo/$cdate_s2$hour/request_pl.yaml +cat > $yaml_file_pl << EOF +class: d1 # For all our runs (for the moemnt) +dataset: climate-dt # For all our runs +activity: CMIP6 # For the moment, maybe it will change in the future +realization: 1 # For current data set to 1 +generation: 1 # For current data set to 1 +experiment: hist # For historical runs +stream: clte # For all our climateDT variables +model: IFS-NEMO # Depends on the model +type: fc # For all our datasets +expver: a0h3 # Depends on the expid of the experiment +param: [t,q,clwc] # [t] - for RADSOUND& [t, q, clwc] - for SATELLITE +date: "$cdate_s2/to/$cdate_s2" # "19900101/to/19971231" # "19910501/to/19910531" +time: "$time_s/to/$time_s/by/0100" # "0000/to/2300/by/0100" - for SATELLITE +levtype: pl # sfc/pl/o2d/o3d/sol/hl -> Check on confluence table +resolution: standard # Standard or high +levelist: [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10, 5 ,1] #- for SATELLITE +grid: $grid +method: nn +EOF +python3 gsv_retriever.py $yaml_file_pl $gsv_depo/$cdate_s2$hour + +yaml_file_sfc=$gsv_depo/$cdate_s2$hour/request_sfc.yaml +cat > $yaml_file_sfc << EOF +class: d1 # For all our runs (for the moemnt) +dataset: climate-dt # For all our runs +activity: CMIP6 # For the moment, maybe it will change in the future +realization: 1 # For current data set to 1 +generation: 1 # For current data set to 1 +experiment: hist # For historical runs +stream: clte # For all our climateDT variables +model: IFS-NEMO # Depends on the model +type: fc # For all our datasets +expver: a0h3 # Depends on the expid of the experiment +param: [2t,2d,skt,sp,10u,10v] # +date: "$cdate_s2/to/$cdate_s2" # "19900101/to/19971231" # "19900101/to/19900131" +time: "$time_s/to/$time_s/by/0100"# "0000/to/2300/by/1200" - for RADSOUND "0000/to/2300/by/0100" - for SYNOP and SATELLITE +levtype: sfc # sfc/pl/o2d/o3d/sol/hl -> Check on confluence table +resolution: standard # Standard or high +grid: $grid +method: nn +EOF +python3 gsv_retriever.py $yaml_file_sfc $gsv_depo/$cdate_s2$hour + +if [ $time_s == '0000' ]; then + yaml_file_o2d=$gsv_depo/$cdate_s2$hour/request_o2d.yaml + cat > $yaml_file_o2d << EOF + class: d1 # For all our runs (for the moemnt) + dataset: climate-dt # For all our runs + activity: CMIP6 # For the moment, maybe it will change in the future + realization: 1 # For current data set to 1 + generation: 1 # For current data set to 1 + experiment: hist # For historical runs + stream: clte # For all our climateDT variables + model: IFS-NEMO # Depends on the model + type: fc # For all our datasets + expver: a0h3 # Depends on the expid of the experiment + param: [263001] # ["vg_siconc"] # required in SATELLITE Part (Note: 263XXX are ocean variables from NEMO) + date: "$cdate_s2/to/$cdate_s2" # "19900101/to/19971231" + time: "0000" # "0000" - since outputting ocean variables as daily averages + levtype: o2d # sfc/pl/o2d/o3d/sol/hl -> Check on confluence table (Note, NEMO: levtype "o2d" (= Ocean 2D) + resolution: standard # Standard or high + grid: $grid + method: nn +EOF + python3 gsv_retriever.py $gsv_depo/$cdate_s2$hour/request_o2d.yaml $gsv_depo/$cdate_s2$hour +else + cp $gsv_depo/${cdate_s2}00/avg_siconc_$res.nc $gsv_depo/$cdate_s2$hour/ +fi + +cdo merge $gsv_depo/$cdate_s2$hour/t_$res.nc $gsv_depo/$cdate_s2$hour/q_$res.nc $gsv_depo/$cdate_s2$hour/clwc_$res.nc $gsv_depo/$cdate_s2$hour/pl.nc + +lsm_static=$scriptdir/other_files/lsm_$res.grb1 #nc +zsurf_static=$scriptdir/other_files/z_surf_$res.grb1 #nc +u10=$gsv_depo/$cdate_s2$hour/10u_$res.nc +u10c=$gsv_depo/$cdate_s2$hour/10u_${res}_cdo_copy.nc +v10=$gsv_depo/$cdate_s2$hour/10v_$res.nc +v10c=$gsv_depo/$cdate_s2$hour/10v_${res}_cdo_copy.nc +t2=$gsv_depo/$cdate_s2$hour/2t_$res.nc +t2c=$gsv_depo/$cdate_s2$hour/2t_${res}_cdo_copy.nc +d2=$gsv_depo/$cdate_s2$hour/2d_$res.nc +d2c=$gsv_depo/$cdate_s2$hour/2d_${res}_cdo_copy.nc +skt=$gsv_depo/$cdate_s2$hour/skt_$res.nc +sktc=$gsv_depo/$cdate_s2$hour/skt_${res}_cdo_copy.nc +sp=$gsv_depo/$cdate_s2$hour/sp_$res.nc +spc=$gsv_depo/$cdate_s2$hour/sp_${res}_cdo_copy.nc +ci=$gsv_depo/$cdate_s2$hour/avg_siconc_$res.nc +cic=$gsv_depo/$cdate_s2$hour/avg_siconc_${res}_cdo_copy.nc + +cdo -f nc copy $u10 $u10c +cdo -f nc copy $v10 $v10c +cdo -f nc copy $t2 $t2c +cdo -f nc copy $d2 $d2c +cdo -f nc copy $skt $sktc +cdo -f nc copy $sp $spc +cdo -f nc copy $ci $cic + +############################################################################################ +# Convert pressure levels to "look-a-like" model levels, because Radiance Simulator does not accept pressure levels. + +# Output file (as netcdf): +intp_file=$gsv_depo/$cdate_s2$hour/ml.nc + +# vct_dat: get A and B coefficients of target vertical grid +vct_dat=$scriptdir/other_files/vct_91.dat + +# run python script: convert_pl2ml.py +python3 convert_pl2ml.py $gsv_depo/$cdate_s2$hour/pl.nc $sp $vct_dat ${intp_file} $lat_S $lat_N $lon_W $lon_E + +$scriptdir/nc2grib.exe < A demo workflow without real AMSU-A processing. + +echo "====================================================================" +echo " DestinE Digital Twin --- OBSALL Apps for AMSU-A OBS DATA:" +echo " TESTING/ RUNNING FOR: " +echo " -- LIMITED GEOGRAPHICAL AREA: BARENTS SEA" +echo " -- FOR AMSU-A CHANNEL 5" +echo "====================================================================" +echo " " + +echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" +echo " STEP #0 : PRE-CHEKING EXISTANCE OF NECESSARY DIRECTORIES ..." +echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" +pwd + +source set_env.sh a +source set_env.sh b + +echo " Directory for pre-processing modelled data extracted with gsv ..." +echo " amsua_process_root: $amsua_process_root" + +export gsv_depo=$amsua_process_root/gsv_depo +echo " Checking that a storage area for GSV files exists ..." +if test -d "$gsv_depo"; then + echo $gsv_depo exists already. +else + echo $gsv_depo does not exist yet. + echo " Creating directory $gsv_depo" + mkdir -p $gsv_depo +fi + +export timestamp_files=$amsua_process_root/timestamp_files +echo " Checking that a storage area for GSV files exists ..." +if test -d "$timestamp_files"; then + echo $timestamp_files exists already. +else + echo $timestamp_files does not exist yet. + echo " Creating directory $timestamp_files" + mkdir -p $timestamp_files +fi + +#export gsv_depo_processed=$gsv_depo/processed +#echo "Checking that a place for GSV pre-processing exists." +#if test -d "$gsv_depo_processed"; then +# echo $gsv_depo_processed exists already. +#else +# echo $gsv_depo_processed does not exist yet. +# echo " Creating directory $gsv_depo_processed" +# mkdir $gsv_depo_processed +#fi + +export amsua_input=$amsua_process_root/amsua_metadata +echo " Checking that a place for observation metadata input files for Radiance Simulator exists ..." +if test -d "$amsua_input"; then + echo $amsua_input exists already. +else + echo $amsua_input does not exist yet. + echo " Creating directory $amsua_input" + mkdir $amsua_input +fi + +export amsua_output=$amsua_process_root/amsua_output +echo " Checking that a place for intermediate output of Radiance Simulator exists ..." +if test -d "$amsua_output"; then + echo $amsua_output exists already. +else + echo $amsua_output does not exist yet. + echo " Creating directory $amsua_output" + mkdir $amsua_output +fi + +echo " Checking that a place for simulated AMSU-A obs in ODB form exists ..." +export final_output_dir +if test -d "$final_output_dir"; then + echo $final_output_dir exists already. +else + echo $final_output_dir does not exist yet. + echo " Creating directory $final_output_dir" + mkdir $final_output_dir +fi + +echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" +echo " ..... STEP 0 - COMPLETED ....." +echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" +#exit 1 + +##test -f nc2grib.exe || ftn -I$MY_ENV/include -L$MY_ENV/lib64 -leccodes_f90 nc2grib.f90 -o nc2grib.exe +#exit 1 +echo " Compiling nc2grib.exe for converting from NetCDF to GRIB format ... " +ftn -I$MY_ENV/include -L$MY_ENV/lib64 -leccodes_f90 nc2grib.f90 -o nc2grib.exe + +tmp_file_tstamps=$timestamp_files/file_with_tstamps_$sdate.txt +rm -f $tmp_file_tstamps +#ls GSVMODDATA/ > $tmp_file_tstamps +##ls $gsv_depo_processed > $tmp_file_tstamps +#echo '2000010100' > $tmp_file_tstamps +#tail -n 1 "$tmp_file_tstamps" | tee >(wc -c | xargs -I {} truncate "$tmp_file_tstamps" -s -{}) +#exit 1 + +echo " Checking sdate, edate, cur_date ..." +echo " sdate: $sdate" +echo " edate: $edate" + +cur_date=$sdate +while [ $cur_date -le $edate ] + do + echo " Adding date: " $cur_date "to the list of steps." + echo $cur_date >> $tmp_file_tstamps + cur_date=$(compute_end_date "$cur_date" "$step_h") + done + +echo " cur_date: $cur_date" + +nrec_file=$( sed -n '$=' $tmp_file_tstamps) +echo "nrec_file: $nrec_file" +for (( nnrec=1; nnrec<=$nrec_file; nnrec++ )); do + b_yy=$( sed -n ${nnrec}p $tmp_file_tstamps | cut -c 1-4 ) + b_mm=$( sed -n ${nnrec}p $tmp_file_tstamps | cut -c 5-6 ) + b_dd=$( sed -n ${nnrec}p $tmp_file_tstamps | cut -c 7-8 ) + b_hh=$( sed -n ${nnrec}p $tmp_file_tstamps | cut -c 9-10 ) + b_min="00" + b_sec="00" + #echo "b_yy, b_mm, b_dd, b_hh : $b_yy $b_mm $b_dd $b_hh" + #tail -n +2 "$tmp_file_tstamps" > "$tmp_file_tstamps.tmp" && mv "$tmp_file_tstamps.tmp" "$tmp_file_tstamps" + # see line above (it is added at the bottom of the script) + + echo " " + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " STEP #1 : PRE-PROCESSING MOD DATA EXTRACTED WITH GSV_INTERFACE ..." + echo " SCRIPT --- gsv_mod_data.sh" $b_yy$b_mm$b_dd $b_hh$b_min + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + ./gsv_mod_data.sh $b_yy$b_mm$b_dd $b_hh$b_min + + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " ..... STEP 1 - COMPLETED ....." + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + echo "====================================================================" + echo " START CALCULATIONS FOR TIME-SLICE: $b_yy $b_mm $b_dd $b_hh " + echo "====================================================================" + + echo " " + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " STEP #2 : PREPARATION OF AMSU-A OBS METADATA ..." + echo " SCRIPT --- amsua_obs.sh" + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + ./amsua_obs.sh $b_yy$b_mm$b_dd$b_hh $b_hh$b_min$b_sec + + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " ..... STEP 2 - COMPLETED ....." + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " STEP #3 : RUNNING RADIATION SIMULATOR OBSERVATION OPERATOR TO " + echo " SIMULATE OBSERVATIONS AT SELECTED POSITIONS" + echo " AND ADDING THE SIMULATED OBSERVATIONS DATA TO ODB ..." + echo " SCRIPT --- amsua_mod.sh" + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + ./amsua_mod.sh $b_yy$b_mm$b_dd$b_hh $b_hh$b_min$b_sec + + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + echo " ..... STEP 3 - COMPLETED ....." + echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" + + echo "====================================================================" + echo " END CALCULATIONS FOR TIME-SLICE: $b_yy $b_mm $b_dd $b_hh " + echo "====================================================================" + + # Removing record with completed timestamp from file ... + # in order to start next time-slice calculations +# tail -n +2 "$tmp_file_tstamps" > "$tmp_file_tstamps.tmp" && mv "$tmp_file_tstamps.tmp" "$tmp_file_tstamps" + + echo " " + echo "====================================================================" + echo " Checking if STATS should be run ..." + echo " STATS - PRODUCING QUANTILE RANK HISTOGRAM STATISTICS AND PLOTS" + echo " SCRIPT --- amsua_stats.sh" + echo " STATS is run on a monthly basis - i.e., done on 2nd day of a month" + echo " ... BUT NOW : $b_yy $b_mm $b_dd $b_hh" + + #if [ $b_dd == "31" ] && [ $b_hh == "00" ]; then + if [ $b_dd == "01" ] && [ $b_hh == "00" ]; then # for testing 2000-02-01 + #b_mm_start="1" + #b_mm_end="12" + + # btmp (https://codes.ecmwf.int/grib/param-db/194) - brightness temperature + varMETname="btmp" + varMETnum="194" + b_yy="2000" #AM #b_yy="2016" #LT + channel="5" # AMSU-A channel number + b_mm_start="01" + b_mm_end="12" + + echo " ... CALCULATING NOW : $b_yy $b_mm $b_dd $b_hh" + echo "====================================================================" + echo " b_yy, b_mm, b_dd, b_hh, b_mm_start, b_mm_end : $b_yy $b_mm $b_dd $b_hh $b_mm_start $b_mm_end" + echo " b_yy, b_yy, b_mm_start, b_mm_end, varMETname, varMETnum, channel : $b_yy $b_yy $b_mm_start $b_mm_end $varMETname $varMETnum $channel" + + #./radsound_stats.sh $varMETnum $b_yy $b_yy $b_mm_start $b_mm_end + #./amsua_stats.sh $b_yy $b_yy $b_mm_start $b_mm_end + ./amsua_stats.sh $b_yy $b_yy $b_mm_start $b_mm_end $varMETname $varMETnum $channel + fi + + #echo " Checking size of file with time stamps: $tmp_file_tstamps ..." + +# file_actualsize=$(wc -c <"$tmp_file_tstamps") +# if [ $file_actualsize == "0" ]; then +# echo size is equal to zero bytes +# echo "====================================================================" +# echo " AMSU-A : END LOOP OVER ALL AVAILABLE TIME-SLICES " +# echo "====================================================================" +# rm $tmp_file_tstamps +# exit 1 +# else + # echo size is over zero bytes + # fi + +done + +#rm nc2grib.exe +exit 1 diff --git a/SATELLITE/nc2grib.f90 b/SATELLITE/nc2grib.f90 new file mode 100644 index 0000000000000000000000000000000000000000..f827b93b5162f4346b611148535c89e16756dd96 --- /dev/null +++ b/SATELLITE/nc2grib.f90 @@ -0,0 +1,408 @@ +program nc2grib + use eccodes + use netcdf + + implicit none + + ! Variables to store dimensions and data + integer :: ncid, lon_dimid, lat_dimid, time_dimid, nx, ny, nl, nt, varid + integer :: i, j, k, level_dimid, level_len + real, allocatable :: t(:,:,:), clwc(:,:,:), q(:,:,:), p(:,:,:) + real, allocatable :: u10(:,:), v10(:,:), d2(:,:), t2(:,:), ci(:,:), skt(:,:), sp(:,:), lsm(:,:), z(:,:) + real, allocatable :: lons(:), lats(:), levs(:), afull(:), bfull(:), coords(:,:) + character(len=160) :: file_name, vname, indir, res, infile3d, indir_static, outfile, ofname + character(len=160) :: coordsfile, grib_samples, outdir + character(len=160) :: file_u10,file_v10,file_d2,file_t2,file_ci,file_skt,file_sp,file_lsm,file_z + integer :: status, xtype, ierr + + namelist/param/infile3d,file_u10,file_v10,file_d2,file_t2,file_ci,file_skt,file_sp, & + file_lsm,file_z,coordsfile,grib_samples,res,outdir,outfile + read(*,nml=param) + + ofname= trim(outdir)//trim(outfile) + + call griddims(infile3d,nt,nl,ny,nx) ! This is the order for Python data + !call griddims(infile3d,nt,nx,ny,nl) ! This is the order for CDO 3D data? + print*,"nc2grib.f90: Got dimensions:",nt,nl,ny,nx + + ! Allocate memory for data + allocate(lons(nx)) + allocate(lats(ny)) + allocate(levs(nl)) + allocate(t(nx, ny, nl)) + allocate(clwc(nx, ny, nl)) + allocate(q(nx, ny, nl)) + allocate(p(nx, ny, nl)) + allocate(u10(nx,ny)) + allocate(v10(nx,ny)) + allocate(d2(nx,ny)) + allocate(t2(nx,ny)) + allocate(ci(nx,ny)) + allocate(skt(nx,ny)) + allocate(sp(nx,ny)) + allocate(lsm(nx,ny)) + allocate(z(nx,ny)) + allocate(coords(3,nl)) + allocate(afull(nl),bfull(nl)) + print*,"nc2grib.f90: Allocation done!" + + call read_vert_coords(coordsfile,nl,coords) + afull=coords(2,:) + bfull=coords(3,:) + + call read_nc_3D(infile3d,lons,lats,levs,t,q,clwc,p,nt,nl,ny,nx) + print*,'nc2grib.f90: Got 3D data!' + call read_nc_2D(file_u10,file_v10,file_d2,file_t2,file_ci,file_skt,file_sp,res,u10,v10,d2,t2,ci,skt,sp,ny,nx) + print*,'nc2grib.f90: Got 2D data!' + call read_nc_2D_static(file_lsm,file_z,res,lsm,z,ny,nx) + print*,'nc2grib.f90: Got 2D static data!' + print*,"nc2grib.f90: All reading done!" + call write_fields(ofname,grib_samples,nx,ny,nl,afull,bfull,t,p,q,clwc, & + u10,v10,d2,t2,ci,skt,sp,lsm,z,ierr) + print*,'nc2grib.f90: Writing done!' + ! Deallocate memory + deallocate(t, clwc, q, p, u10, v10, d2, t2, ci, skt, sp, lsm, z) + +contains + + SUBROUTINE griddims(infile,nt,nl,ny,nx) + IMPLICIT NONE + INTEGER, INTENT(OUT) :: nx, ny, nl, nt + INTEGER :: ncid + CHARACTER(LEN=160), INTENT(IN) :: infile + CHARACTER(LEN=160) :: xname, yname, tname, lname + !Fortran 90 netCDF introduction + !Open netCDF file + !:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_open(infile, nf90_nowrite, ncid)) + !Inquire about the dimensions + !:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_inquire_dimension(ncid,1,tname,nt)) + CALL check(nf90_inquire_dimension(ncid,2,lname,nl)) + CALL check(nf90_inquire_dimension(ncid,3,yname,ny)) + CALL check(nf90_inquire_dimension(ncid,4,xname,nx)) + !Close netCDF file + !:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_close(ncid)) + END SUBROUTINE griddims + + SUBROUTINE read_nc_3D(infile,lons,lats,levs,t,q,clwc,p,nt,nl,ny,nx) + IMPLICIT NONE + REAL, DIMENSION(nx), INTENT(OUT) :: lats + REAL, DIMENSION(ny), INTENT(OUT) :: lons + REAL, DIMENSION(nl), INTENT(OUT) :: levs + REAL, DIMENSION(nt) :: times + REAL, DIMENSION(nx,ny,nl), INTENT(OUT) :: t,q,clwc,p + INTEGER, INTENT(IN) :: nx, ny, nl, nt + INTEGER, DIMENSION(4) :: dimids + INTEGER :: ncid, xtype, ndims, varid, i + CHARACTER(LEN=160), INTENT(IN) :: infile + CHARACTER(LEN=160) :: xname, yname, vname + !Open netCDF file + !:-------:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_open(infile, nf90_nowrite, ncid)) + !Get the values of the coordinates and put them in xpos & ypos + !:-------:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_inquire_variable(ncid,1,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,times)) ! 0. OK? + CALL check(nf90_inquire_variable(ncid,2,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,levs)) ! OK + CALL check(nf90_inquire_variable(ncid,3,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,lats)) ! not OK, but no big deal. + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,lons)) ! OK + + !Get the data + !:-------:-------:-------:-------:-------:-------:-------:-------: + print*,'nc2grib.f90: Start reading 3D data.' + CALL check(nf90_inquire_variable(ncid,5,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,q)) + print*,' Q(lowest level), min=',minval(q(:,:,91)),'max=',maxval(q(:,:,91)) + CALL check(nf90_inquire_variable(ncid,6,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,clwc)) + print*,' CLWC(lowest level), min=',minval(clwc(:,:,91)),'max=',maxval(clwc(:,:,91)) + CALL check(nf90_inquire_variable(ncid,7,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,t)) + print*,' T(lowest level), min=',minval(t(:,:,91)),'max=',maxval(t(:,:,91)) + CALL check(nf90_inquire_variable(ncid,8,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,p)) + print*,' P(lowest level), min=',minval(p(:,:,91)),'max=',maxval(p(:,:,91)) + + !Close netCDF file + !:-------:-------:-------:-------:-------:-------:-------:-------: + CALL check(nf90_close(ncid)) + END SUBROUTINE read_nc_3D + + subroutine read_nc_2D(file_u10,file_v10,file_d2,file_t2,file_ci,file_skt,file_sp, & + res,u10,v10,d2,t2,ci,skt,sp,ny,nx) + implicit none + character(len=160), intent(in) :: res, file_u10, file_v10, file_d2, file_t2, file_ci, file_skt, file_sp + INTEGER, INTENT(IN) :: nx, ny + REAL, DIMENSION(nx,ny), INTENT(OUT) :: u10,v10,d2,t2,ci,skt,sp + INTEGER, DIMENSION(3) :: dimids + INTEGER :: ncid, xtype, ndims, varid + CHARACTER(LEN=160) :: vname + + print*,'nc2grib.f90: Start reading 2D data.' + CALL check(nf90_open(file_u10, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,u10)) + CALL check(nf90_close(ncid)) + print*,' 10U, min=',minval(u10(:,:)),'max=',maxval(u10(:,:)) + + CALL check(nf90_open(file_v10, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,v10)) + CALL check(nf90_close(ncid)) + print*,' 10V, min=',minval(v10(:,:)),'max=',maxval(v10(:,:)) + + CALL check(nf90_open(file_d2, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,d2)) + CALL check(nf90_close(ncid)) + print*,' 2D, min=',minval(d2(:,:)),'max=',maxval(d2(:,:)) + + CALL check(nf90_open(file_t2, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,t2)) + CALL check(nf90_close(ncid)) + print*,' 2T, min=',minval(t2(:,:)),'max=',maxval(t2(:,:)) + + CALL check(nf90_open(file_ci, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,ci)) + CALL check(nf90_close(ncid)) + print*,' CI, min=',minval(ci(:,:)),'max=',maxval(ci(:,:)) + + CALL check(nf90_open(file_skt, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,skt)) + CALL check(nf90_close(ncid)) + print*,' SKT, min=',minval(skt(:,:)),'max=',maxval(skt(:,:)) + + CALL check(nf90_open(file_sp, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,sp)) + CALL check(nf90_close(ncid)) + print*,' SP, min=',minval(sp(:,:)),'max=',maxval(sp(:,:)) + + end subroutine read_nc_2D + + subroutine read_nc_2D_static(file_lsm,file_z,res,lsm,z,ny,nx) + implicit none + CHARACTER(LEN=160), INTENT(IN) :: file_lsm, file_z, res + INTEGER, INTENT(IN) :: nx, ny + REAL, DIMENSION(nx,ny), INTENT(OUT) :: lsm, z + INTEGER, DIMENSION(3) :: dimids + INTEGER :: ncid, xtype, ndims, varid + CHARACTER(LEN=160) :: vname + + print*,'nc2grib.f90: Start reading 2D static data.' + CALL check(nf90_open(file_lsm, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,lsm)) + CALL check(nf90_close(ncid)) + print*,' LSM, min=',minval(lsm(:,:)),'max=',maxval(lsm(:,:)) + + CALL check(nf90_open(file_z, nf90_nowrite, ncid)) + CALL check(nf90_inquire_variable(ncid,4,vname,xtype,ndims,dimids)) + CALL check(nf90_inq_varid(ncid,vname,varid)) + CALL check(nf90_get_var(ncid,varid,z)) + CALL check(nf90_close(ncid)) + print*,' Z_SURF, min=',minval(z(:,:)),'max=',maxval(z(:,:)) + + end subroutine read_nc_2D_static + + SUBROUTINE check(istatus) + IMPLICIT NONE + INTEGER, INTENT (IN) :: istatus + IF (istatus /= nf90_noerr) THEN + write(*,*) TRIM(ADJUSTL(nf90_strerror(istatus))) + END IF + END SUBROUTINE check + + subroutine write_fields(ofname,samplefile,nx,ny,nlev,afull,bfull,t,p,q,clwc, & + u10,v10,d2,t2,ci,skt,sp,lsm,z,ierr) + character(len=160) :: ofname,samplefile + integer :: nx, ny, nlev, ierr, k + real :: afull(nlev+1), bfull(nlev+1) + real :: t(nx,ny,nlev),p(nx,ny,nlev),q(nx,ny,nlev),clwc(nx,ny,nlev) + real :: u10(nx,ny),v10(nx,ny),d2(nx,ny),t2(nx,ny),ci(nx,ny),skt(nx,ny),sp(nx,ny),lsm(nx,ny),z(nx,ny) + integer :: infilecode, outfilecode, igrib_in, igrib_out + real, allocatable :: ci_v(:) + real :: a1 + + allocate(ci_v(nx*ny)) + + call codes_open_file(infilecode,samplefile,'r') + call codes_open_file(outfilecode,ofname,'w') + call codes_grib_new_from_file(infilecode,igrib_in) + call codes_clone(igrib_in,igrib_out) + + ! 3D temperature + do k=1,nlev + call codes_set(igrib_out,'values',pack(t(:,:,k),mask=.true.)) + call codes_set(igrib_out,'level',k) + call codes_set(igrib_out,'indicatorOfParameter',130) + call codes_set(igrib_out,'paramId',130) + call codes_set(igrib_out,'shortName','t') + call codes_write(igrib_out,outfilecode) + enddo + + ! 3D pressure + do k=1,nlev + call codes_set(igrib_out,'values',pack(p(:,:,k),mask=.true.)) + call codes_set(igrib_out,'level',k) + call codes_set(igrib_out,'indicatorOfParameter',54) + call codes_set(igrib_out,'paramId',54) + call codes_set(igrib_out,'shortName','pres') + call codes_write(igrib_out,outfilecode) + enddo + + ! 3D specific humidity + do k=1,nlev + call codes_set(igrib_out,'values',pack(q(:,:,k),mask=.true.)) + call codes_set(igrib_out,'level',k) + call codes_set(igrib_out,'indicatorOfParameter',133) + call codes_set(igrib_out,'paramId',133) + call codes_set(igrib_out,'shortName','q') + call codes_write(igrib_out,outfilecode) + enddo + + ! 3d cloud liquid water content + do k=1,nlev + call codes_set(igrib_out,'values',pack(clwc(:,:,k),mask=.true.)) + call codes_set(igrib_out,'level',k) + call codes_set(igrib_out,'indicatorOfParameter',246) + call codes_set(igrib_out,'paramId',246) + call codes_set(igrib_out,'shortName','clwc') + call codes_write(igrib_out,outfilecode) + enddo + + ! 2D surface pressure + call codes_set(igrib_out,'values',pack(sp,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',134) + call codes_set(igrib_out,'paramId',134) + call codes_set(igrib_out,'shortName','sp') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D skin temperature + call codes_set(igrib_out,'values',pack(skt,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',235) + call codes_set(igrib_out,'paramId',235) + call codes_set(igrib_out,'shortName','skt') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D sea ice cover + ci_v=pack(ci,mask=.true.) + do i = 1,size(ci_v) + if (isnan(ci_v(i))) then + ci_v(i)=9999.0 + endif + enddo + call codes_set(igrib_out,'values',ci_v) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',31) + call codes_set(igrib_out,'paramId',31) + call codes_set(igrib_out,'shortName','ci') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D 10m U wind + call codes_set(igrib_out,'values',pack(u10,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',165) + call codes_set(igrib_out,'paramId',165) + call codes_set(igrib_out,'shortName','10u') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D 10m V wind + call codes_set(igrib_out,'values',pack(v10,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',166) + call codes_set(igrib_out,'paramId',166) + call codes_set(igrib_out,'shortName','10v') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D 2m temperature + call codes_set(igrib_out,'values',pack(t2,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',167) + call codes_set(igrib_out,'paramId',167) + call codes_set(igrib_out,'shortName','2t') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D 2m dewpoint temperature + call codes_set(igrib_out,'values',pack(d2,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',168) + call codes_set(igrib_out,'paramId',168) + call codes_set(igrib_out,'shortName','2d') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D land-sea mask + call codes_set(igrib_out,'values',pack(lsm,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',172) + call codes_set(igrib_out,'paramId',172) + call codes_set(igrib_out,'shortName','lsm') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + ! 2D surface geopotential (orography) + call codes_set(igrib_out,'values',pack(z,mask=.true.)) + call codes_set(igrib_out,'level',0) + call codes_set(igrib_out,'indicatorOfParameter',129) + call codes_set(igrib_out,'paramId',129) + call codes_set(igrib_out,'shortName','z') + call codes_set(igrib_out,'typeOfLevel','surface') + call codes_write(igrib_out,outfilecode) + + call codes_release(igrib_out) + call codes_release(igrib_in) + call codes_close_file(infilecode) + call codes_close_file(outfilecode) + + end subroutine write_fields + + subroutine read_vert_coords(coordsfile,nlev,coords) + character(len=160) :: coordsfile, line + integer :: nlev, i + real,dimension(3,nlev+1) :: coords + open(unit=7, file=coordsfile) + i = 1 + do + read(7,'(A)',end=11) line + line = adjustl(line) + if (line(1:1) == '#') cycle + read(line,*) coords(:,i) + i = i +1 + enddo + 11 close(7) + end subroutine read_vert_coords + +end program nc2grib diff --git a/SATELLITE/other_files/README_other_files.txt b/SATELLITE/other_files/README_other_files.txt new file mode 100644 index 0000000000000000000000000000000000000000..4613fe19930bf7e56586fe7e8e0e8c66eb871272 --- /dev/null +++ b/SATELLITE/other_files/README_other_files.txt @@ -0,0 +1,28 @@ +### Some info about the other files required by AMSU-A workflow. ### + +*** amsua_locations_sample.dat +- If the process mode is "cli", a static user defined list of locations is used. The file contains lon-lat points where AMSU-A irradiances are computed. + +*** lsm_0.25.grib2 +- Land-sea mask of IFS on 0.25 x 0.25 degree grid. Pole points are included. + +*** lsm_r180x90.nc +- Land-sea mask of IFS interpolated on r180x90 (2 x 2 deg) grid and converted to netcdf with CDO. Pole points are excluded. This can be used as input for the AMSU-A input file construction with nc2grib.f90. + +*** lsm_r1440x720.nc +- Land-sea mask of IFS interpolated on r1440x720 (0.25 x 0.25 deg) grid and converted to netcdf with CDO. Pole points are excluded. This can be used as input for the AMSU-A input file construction with nc2grib.f90. + +*** sample_r1440x720.grib1 +- A sample file from which the grib entries are cloned when constructing the AMSU-A input file with nc2grib.f90. + +*** vct_91.dat +- IFS A- and B-coefficients for constructing 3D hybrid sigma grid. + +*** z_surf_0.25.grib2 +- Surface orography of IFS on 0.25 x 0.25 degree grid. Pole points are included. + +*** z_surf_r180x90.nc +- Surface orography of IFS interpolated on r180x90 (2 x 2 deg) grid and converted to netcdf with CDO. Pole points are excluded. This can be used as input for the AMSU-A input file construction with nc2grib.f90. + +*** z_surf_r1440x720.nc +- Surface orography of IFS interpolated on r1440x720 (0.25 x 0.25 deg) grid and converted to netcdf with CDO. Pole points are excluded. This can be used as input for the AMSU-A input file construction with nc2grib.f90. diff --git a/SATELLITE/other_files/amsua_locations_sample.dat b/SATELLITE/other_files/amsua_locations_sample.dat new file mode 100644 index 0000000000000000000000000000000000000000..0a30f48966dec062efb8580d96ebaced6af5820f --- /dev/null +++ b/SATELLITE/other_files/amsua_locations_sample.dat @@ -0,0 +1,5 @@ +2 +1 +5 +1 2 15 5 4 +45.000000 70.000000 120.000000 0.0000000 0.000000 diff --git a/SATELLITE/other_files/lsm_0.25.grib2 b/SATELLITE/other_files/lsm_0.25.grib2 new file mode 100644 index 0000000000000000000000000000000000000000..608d55adb9cb38b203273bf22711de239a325134 Binary files /dev/null and b/SATELLITE/other_files/lsm_0.25.grib2 differ diff --git a/SATELLITE/other_files/lsm_r1440x720.nc b/SATELLITE/other_files/lsm_r1440x720.nc new file mode 100644 index 0000000000000000000000000000000000000000..43f44e2750926bf9acafdb24b56652fbd5da53ad Binary files /dev/null and b/SATELLITE/other_files/lsm_r1440x720.nc differ diff --git a/SATELLITE/other_files/lsm_r180x90.nc b/SATELLITE/other_files/lsm_r180x90.nc new file mode 100644 index 0000000000000000000000000000000000000000..8b7acd2f7245278210d56b36c3007a576adc8011 Binary files /dev/null and b/SATELLITE/other_files/lsm_r180x90.nc differ diff --git a/SATELLITE/other_files/sample_r1440x720.grib1 b/SATELLITE/other_files/sample_r1440x720.grib1 new file mode 100644 index 0000000000000000000000000000000000000000..f8fbdbad4324250751ed3154aec4e6f84c307073 Binary files /dev/null and b/SATELLITE/other_files/sample_r1440x720.grib1 differ diff --git a/SATELLITE/other_files/vct_91.dat b/SATELLITE/other_files/vct_91.dat new file mode 100644 index 0000000000000000000000000000000000000000..85d6f4e23129e85afacac814a60d5c2d318e1604 --- /dev/null +++ b/SATELLITE/other_files/vct_91.dat @@ -0,0 +1,93 @@ +# k vct_a(k) [Pa] vct_b(k) [] + 0 0.00000000000000000 0.00000000000000000 + 1 2.00004005432128906 0.00000000000000000 + 2 3.98083209991455078 0.00000000000000000 + 3 7.38718605041503906 0.00000000000000000 + 4 12.90831947326660156 0.00000000000000000 + 5 21.41361236572265625 0.00000000000000000 + 6 33.95285797119140625 0.00000000000000000 + 7 51.74660110473632812 0.00000000000000000 + 8 76.16765594482421875 0.00000000000000000 + 9 108.71556091308593750 0.00000000000000000 + 10 150.98602294921875000 0.00000000000000000 + 11 204.63745117187500000 0.00000000000000000 + 12 271.35650634765625000 0.00000000000000000 + 13 352.82449340820312500 0.00000000000000000 + 14 450.68579101562500000 0.00000000000000000 + 15 566.51922607421875000 0.00000000000000000 + 16 701.81335449218750000 0.00000000000000000 + 17 857.94580078125000000 0.00000000000000000 + 18 1036.16650390625000000 0.00000000000000000 + 19 1237.58544921875000000 0.00000000000000000 + 20 1463.16394042968750000 0.00000000000000000 + 21 1713.70959472656250000 0.00000000000000000 + 22 1989.87438964843750000 0.00000000000000000 + 23 2292.15551757812500000 0.00000000000000000 + 24 2620.89843750000000000 0.00000000000000000 + 25 2976.30224609375000000 0.00000000000000000 + 26 3358.42578125000000000 0.00000000000000000 + 27 3767.19604492187500000 0.00000000000000000 + 28 4202.41650390625000000 0.00000000000000000 + 29 4663.77636718750000000 0.00000000000000000 + 30 5150.85986328125000000 0.00000000000000000 + 31 5663.15625000000000000 0.00000000000000000 + 32 6199.83935546875000000 0.00000000000000000 + 33 6759.72705078125000000 0.00000000000000000 + 34 7341.46972656250000000 0.00000027240000122 + 35 7942.92626953125000000 0.00001391160003550 + 36 8564.62402343750000000 0.00005466720176628 + 37 9208.30566406250000000 0.00013136409688741 + 38 9873.56054687500000000 0.00027888480690308 + 39 10558.88183593750000000 0.00054838409414515 + 40 11262.48437500000000000 0.00100013439077884 + 41 11982.66210937500000000 0.00170107535086572 + 42 12713.89746093750000000 0.00276471930555999 + 43 13453.22558593750000000 0.00426704855635762 + 44 14192.00976562500000000 0.00632216688245535 + 45 14922.68554687500000000 0.00903499033302069 + 46 15638.05371093750000000 0.01250826288014650 + 47 16329.56054687500000000 0.01685957796871662 + 48 16990.62304687500000000 0.02218864485621452 + 49 17613.28125000000000000 0.02861034870147705 + 50 18191.02929687500000000 0.03622690960764885 + 51 18716.96875000000000000 0.04514613375067711 + 52 19184.54492187500000000 0.05547422915697098 + 53 19587.51367187500000000 0.06731618940830231 + 54 19919.79687500000000000 0.08077730238437653 + 55 20175.39453125000000000 0.09596405923366547 + 56 20348.91601562500000000 0.11297900974750519 + 57 20434.15820312500000000 0.13193482160568237 + 58 20426.21875000000000000 0.15293352305889130 + 59 20319.01171875000000000 0.17609108984470367 + 60 20107.03125000000000000 0.20152013003826141 + 61 19785.35742187500000000 0.22931487858295441 + 62 19348.77539062500000000 0.25955447554588318 + 63 18798.82226562500000000 0.29199340939521790 + 64 18141.29687500000000000 0.32632938027381897 + 65 17385.59570312500000000 0.36220258474349976 + 66 16544.58593750000000000 0.39920476078987122 + 67 15633.56640625000000000 0.43690633773803711 + 68 14665.64550781250000000 0.47501641511917114 + 69 13653.21972656250000000 0.51327973604202271 + 70 12608.38378906250000000 0.55145847797393799 + 71 11543.16699218750000000 0.58931714296340942 + 72 10471.31054687500000000 0.62655889987945557 + 73 9405.22265625000000000 0.66293358802795410 + 74 8356.25292968750000000 0.69822359085083008 + 75 7335.16455078125000000 0.73222380876541138 + 76 6353.92089843750000000 0.76467949151992798 + 77 5422.80273437500000000 0.79538476467132568 + 78 4550.21582031250000000 0.82418543100357056 + 79 3743.46435546875000000 0.85095041990280151 + 80 3010.14697265625000000 0.87551838159561157 + 81 2356.20263671875000000 0.89776724576950073 + 82 1784.85461425781250000 0.91765093803405762 + 83 1297.65612792968750000 0.93515706062316895 + 84 895.19354248046875000 0.95027381181716919 + 85 576.31414794921875000 0.96300709247589111 + 86 336.77236938476562500 0.97346603870391846 + 87 162.04342651367187500 0.98223811388015747 + 88 54.20833587646484375 0.98915296792984009 + 89 6.57562780380249023 0.99420416355133057 + 90 0.00316000008024275 0.99763011932373047 + 91 0.00000000000000000 1.00000000000000000 diff --git a/SATELLITE/other_files/z_surf_0.25.grib2 b/SATELLITE/other_files/z_surf_0.25.grib2 new file mode 100644 index 0000000000000000000000000000000000000000..a14645fbe45e6859d89a8066df3eff0fe2d023b3 Binary files /dev/null and b/SATELLITE/other_files/z_surf_0.25.grib2 differ diff --git a/SATELLITE/other_files/z_surf_r1440x720.nc b/SATELLITE/other_files/z_surf_r1440x720.nc new file mode 100644 index 0000000000000000000000000000000000000000..9adb6964eecbba86a11db244562ccf28038e9ea4 Binary files /dev/null and b/SATELLITE/other_files/z_surf_r1440x720.nc differ diff --git a/SATELLITE/other_files/z_surf_r180x90.nc b/SATELLITE/other_files/z_surf_r180x90.nc new file mode 100644 index 0000000000000000000000000000000000000000..75d5319ee5fa74f7e374ee277b9762824d288aaa Binary files /dev/null and b/SATELLITE/other_files/z_surf_r180x90.nc differ diff --git a/SATELLITE/set_env.sh b/SATELLITE/set_env.sh new file mode 100755 index 0000000000000000000000000000000000000000..6fe01415ebf967a1584cf89840c35c5e3fd815af --- /dev/null +++ b/SATELLITE/set_env.sh @@ -0,0 +1,314 @@ +#!/bin/bash + +part=${1?} + +################################################################### + +echo "#######################################################" +echo "### SATELLITE PART ###" +echo "### START Part a: ###" +echo "### Define variables and load modules ###" +echo "#######################################################" + +if [ $part == 'a' ]; then + + ### The place where you have put the AMSU-A workflow scripts. ### + #scriptdir=/projappl/project_465000454/lautuppi/for_Alexander/AMSUA_workflow + scriptdir=/projappl/project_465000454/ama/AMSUA + + ### Set some variables that define where you want the workflow to run and process data. ### + # You can change these to be whatever you like. + + # Your scratch directory for data. + ##scratch=/scratch/project_465000454/lautuppi #LT + echo " Scratch directory with required data for SATELLITE Part ..." + scratch=/scratch/project_465000454/ama + echo " scratch: $scratch" + + # writing oputput into dirs + # /scratch/project_465000454/lautuppi/AMSUA_demo_wf --- for temporary output + # /scratch/project_465000454/lautuppi/AMSUA_demo_ODB --- for final output + + # GSV data is prepared for AMSU-A computation in here. + # Also other temporary intermediate files appear here. + echo " GSV data prepared for AMSU-A Part computation ... writing temporary output to dir ..." + amsua_process_root=$scratch/AMSUA_demo_wf + echo " amsua_process_root: $amsua_process_root" + + # A place for output ODB files and output ODB file name. + echo " Path to saving output ODB files for AMSU-A Part..." + final_output_dir=$scratch/AMSUA_demo_ODB + echo " final_output_dir: $final_output_dir" + echo " Name of output ODB file for saving model data for satellite Part ..." + mod_amsua_odb=$final_output_dir/amsua_mod_data_2000.odb + echo " mod_amsua_odb: $mod_amsua_odb" + + ### Nudged storyline simulations require actual observations that are stored here. ### + echo " Nudged storyline simulations require actual observations ..." + amsua_obs_database=/scratch/project_465000454/lautuppi/AMSU-A_observations/channel_5/daily_txt + echo " amsua_obs_database: $amsua_obs_database" + + # IFS-NEMO simulations - experiment name. + echo " Experiment name (from IFS-NEMO simulations) ..." + exp=a0h3 + echo " expid: $exp" + + # Start and end dates for a chunk to be processed; length of a GSV timestep in hours; option to remove intermediate files on the fly. + sdate=2000020100 # Please always start at 00UTC + edate=2000020101 #2000010500 + + step_h=1 + remove_files=1 # 1 means remove files. + + # Select an area for processing. + echo " Selected geographical area for processing ..." + lat_S=69.5 + lat_N=70.5 + lon_W=43.5 + lon_E=46.5 + echo " lat_S= $lat_S lat_N= $lat_N lon_W= $lon_W lon_E= $lon_W" + + # Resolution for GSV interface output. + # Note that res has to be defined manually, and grid and res must match. + # Only regular lon-lat is acceptable. + echo " Selected resolution for gsv extracted modeled data ..." + grid='0.25/0.25' + res=r1440x720 + echo " $grid deg/deg" + + # Set processing mode: "cli" for free-running climate simulation and "rf" for the storylines simulations. + process_mode=cli + + # Set path to coefficients required by Radiance Simulator + echo " Path to coefficients/constants required by Radiance Simulator software ..." + rttov_coefs=/projappl/project_465000454/lautuppi/software/radiance_simulator/v3.2/RTTOV/rtcoef_rttov13 + echo " rttov_coefs: $rttov_coefs" + + # FDB5 configuration file corresponding the selected experiment. + #FDB5_CONFIG_FILE=/projappl/project_465000454/experiments/a07q/fdb/config.yaml + echo " Path to configuration yaml-file corresponding to selected expid ..." + FDB5_CONFIG_FILE=/projappl/project_465000454/experiments/$exp/config.yaml + echo " FDB5_CONFIG_FILE: $FDB5_CONFIG_FILE" + + # GSV Configuration files + echo " Path to GSV configuration files ..." + GSV_WEIGHTS_PATH=/scratch/project_465000454/igonzalez/gsv_weights + GSV_TEST_FILES=/scratch/project_465000454/igonzalez/gsv_test_files + GRID_DEFINITION_PATH=/scratch/project_465000454/igonzalez/grid_definitions + echo " GSV_WEIGHTS_PATH: $GSV_WEIGHTS_PATH" + echo " GRID_DEFINITION_PATH: $GRID_DEFINITION_PATH" + + # Put the GSV interface tool on Python search path. + echo " Path to GSV interface tool ..." + GSV_SRC_DIR=/projappl/project_465000454/lautuppi/gsv_interface/ + echo " GSV_SRC_DIR: $GSV_SRC_DIR" + PYTHONPATH=$GSV_SRC_DIR:$PYTHONPATH + + # Other paths required for nc2grib.f90 file conversion tool. + echo " Paths required for conversion file from NetCDF to GRIB format (nc2grib.f90) ..." + MY_ENV=/projappl/project_465000454/jlento/PrgEnv-gnu/8.4.0 + echo " $MY_ENV" + PATH=$MY_ENV/bin:$PATH + CPATH=$MY_ENV/include:$CPATH + LIBRARY_PATH=$MY_ENV/lib64:$LIBRARY_PATH + LD_LIBRARY_PATH=$MY_ENV/lib64:$LD_LIBRARY_PATH + + #### Not sure if this stuff will be needed ############ + ############# Perhaps useful in the storylines mode. ########## + # These times should correspond to the GSV timestep. cdate_start is the time t-30min in YYYYMMDDHH, cdate_mid is the time t, cdate_end is time t+30min. start_time is t-30 min in HHMMSS and time range is the timestep length. + # cdate_mid=2000010100 #2016120100 + #cdate_end=2016120100 + #start_time=233000 +# time_range=60 +# +echo " Exporting needed environment variables (i.e., global variables accessible by all processes) ..." + export scratch scriptdir amsua_process_root amsua_obs_database final_output_dir mod_amsua_odb mother_GSV_file cdate_mid time_range process_mode radsim_run rttov_coefs remove_files grid res FDB5_CONFIG_FILE GSV_WEIGHTS_PATH GSV_TEST_FILES GRID_DEFINITION_PATH PYTHONPATH MY_ENV PATH CPATH LIBRARY_PATH LD_LIBRARY_PATH + export sdate edate step_h lat_S lat_N lon_W lon_E + + # Load modules etc. + echo " Loading required modules ..." + module use /project/project_465000454/devaraju/modules/LUMI/23.03/C + module load LUMI/23.03 + module load partition/C + module load ecCodes/2.32.0-cpeCray-23.03 + module load pyfdb/0.0.2-cpeCray-23.03 + module load python-climatedt/3.11.3-cpeCray-23.03.lua + module load odb_api/0.18.1-cpeCray-23.03.lua + module load cray-hdf5/1.12.2.3 + module load cray-netcdf/4.9.0.3 + module load rttov/13.2 + module load radsim/3.2 + module load fdb/5.11.94-cpeCray-23.03 + module load eckit/1.25.0-cpeCray-23.03 + module load metkit/1.11.0-cpeCray-23.03 + module load python-climatedt/3.11.7-cpeCray-23.03 + module load PrgEnv-gnu + #### module load python-climatedt/3.11.3-cpeCray-23.03.lua ### older version + #### clarify pyfdb & fdb modules loaded ??? +fi + +echo "#######################################################" +echo "### SATELLITE PART ###" +echo "### END Part a: ###" +echo "### Define a set of costomized bash functions ###" +echo "#######################################################" + +#################################################################### + +echo "#######################################################" +echo "### SATELLITE PART ###" +echo "### START Part b: ###" +echo "### Define a set of costomized bash functions ###" +echo "#######################################################" + +if [ $part == 'b' ]; then + + # Function to convert HHMMSS to minutes since midnight + echo " Function to convert HHMMSS to minutes since midnight ..." + + hhmmss_to_minutes() { + local hhmmss=$1 + local hh=${hhmmss:0:2} + local mm=${hhmmss:2:2} + local ss=${hhmmss:4:2} + echo $((10#$hh * 60 + 10#$mm + 10#$ss / 60)) + } + + # Function to convert minutes since midnight to HHMMSS + echo " Function to convert minutes since midnight to HHMMSS ..." + + minutes_to_hhmmss() { + local minutes=$1 + if [ $minutes -ge 0 ]; then + local hh=$((minutes / 60)) + if [ $hh -ge 24 ]; then + local hh=$((hh - 24)) + fi + local mm=$((minutes % 60)) + elif [ $minutes -lt 0 ]; then + local hh=$((23 + minutes / 60)) + if [ $hh -lt 0 ]; then + local hh=$((hh + 24)) + fi + local mm=$((-1 * minutes % 60)) + fi + printf "%02d%02d00\n" $hh $mm + } + + # Function to check if a year is a leap year + echo " Function to check if a year is a leap year ..." + + is_leap_year() { + local year=$1 + if ((year % 4 == 0 && year % 100 != 0)) || ((year % 400 == 0)); then + return 0 + else + return 1 + fi + } + + # Function to calculate the number of days in a given month of a specific year + echo " Function to calculate the number of days in a given month of a specific year ..." + + days_in_month() { + local month=$1 + local year=$2 + month_c1=$(echo $end_month | cut -c 1 ) + if [ $month_c1 -eq 0 ]; then + month=$(echo $month | cut -c 2 ) + fi + case $month in + 1|3|5|7|8|10|12) echo 31 ;; # January, March, May, July, August, October, December + 4|6|9|11) echo 30 ;; # April, June, September, November + 2) if is_leap_year "$year"; then + echo 29 # February in a leap year + else + echo 28 # February in a non-leap year + fi + ;; + *) echo "Invalid month: $month"; exit 1 ;; + esac + } + + # Function to compute the end date based on initial date (YYYYMMDDHH) and time step (hours) + echo " Function to compute the end date based on initial date (YYYYMMDDHH) and time step (hours) ..." + + compute_end_date() { + local initial_date=$1 + local time_step=$2 + + local year=${initial_date:0:4} + local month=${initial_date:4:2} + local day=${initial_date:6:2} + local hour=${initial_date:8:2} + + # Calculate the end date + echo " Calculate the end date ..." + + local end_year=$year + local end_month=$month + local end_day=$day + local end_hour=$((10#$hour + 10#$time_step)) + + # Handle negative time steps + echo "Handle negative time steps ..." + + while ((end_hour < 0)); do + end_hour=$((end_hour + 24)) + end_day=$((10#$end_day - 1)) + if ((end_day == 0)); then + end_month=$((10#$end_month - 1)) + if ((end_month == 0)); then + end_month=12 + end_year=$((end_year - 1)) + fi + end_day=$(days_in_month $end_month $end_year) + fi + done + + # Handle positive time steps + echo " Handle positive time steps ..." + + while ((end_hour >= 24)); do + end_hour=$((end_hour - 24)) + end_day=$((10#$end_day + 1)) + local days_in_current_month=$(days_in_month $end_month $end_year) + if ((end_day > days_in_current_month)); then + end_day=1 + end_month=$((10#$end_month + 1)) + if ((end_month > 12)); then + end_month=1 + end_year=$((end_year + 1)) + fi + fi + done + + # Format end date + echo " Format end date ..." + + if [ $end_month -lt 10 ]; then + end_month="0"$(echo $end_month) + fi + if [ ${#end_month} -gt 2 ]; then + end_month=$(echo $end_month | cut -c 2-3 ) + fi + if [ $end_day -lt 10 ]; then + end_day="0"$(echo $end_day) + fi + if [ ${#end_day} -gt 2 ]; then + end_day=$(echo $end_day | cut -c 2-3 ) + fi + if [ $end_hour -lt 10 ]; then + end_hour="0"$(echo $end_hour ) + fi + echo $end_year$end_month$end_day$end_hour + } +fi + +echo "#######################################################" +echo "### SATELLITE PART ###" +echo "### END Part b: ###" +echo "### Define a set of costomized bash functions ###" +echo "#######################################################" + +################################################################