program read_airs implicit none integer :: i,ii,j,k,kk,l,l2,nr,iuf,rec, recl integer :: num,num_pres_1, num_pres_2,num_index_t, num_index_mr integer :: analysis_time, time_window integer, parameter :: leng=337, nrec=12150, hr_bufr=24 integer :: hr_airs integer,dimension(5):: idate5 integer,dimension(5):: iadate integer :: nmind,gstime integer :: num_message real :: tdiff,twind,sstime integer, parameter :: mxmn=35, mxlv=62 integer :: mxlv2 character(80):: hdstr='SID XOB YOB DHR TYP ELV SAID T29' character(80):: obstr='POB QOB TOB ZOB UOB VOB PWO CAT PRSS' character(80):: qcstr='PQM QQM TQM ZQM WQM NUL PWQ ' character(80):: oestr='POE QOE TOE NUL WOE NUL PWE ' real(8) :: hdr(mxmn),obs(mxmn,mxlv),qcf(mxmn,mxlv),oer(mxmn,mxlv) character(8) :: subset integer :: unit_out=10,unit_table=20,idate,iret,nlvl character(8) :: c_sid real(8) :: rstation_id equivalence(rstation_id,c_sid) ! real*4 prof(leng) ! real*4 prof_all(nrec,leng,hr_airs) ! real*4 selected_prof_all(100000,leng) real*4 pres(101) ! character*200 airs_file_path_and_name(hr_airs) character*200 prepbufr_file_path_and_name(hr_bufr) REAL, ALLOCATABLE, DIMENSION(:) :: prof REAL, ALLOCATABLE, DIMENSION(:,:,:) :: prof_all REAL, ALLOCATABLE, DIMENSION(:,:) :: selected_prof_all CHARACTER*200, ALLOCATABLE,DIMENSION(:) :: airs_file_path_and_name character(4) :: year_c character(2) :: month_c, day_c, hour_c character(10) :: idate_c real :: year, month, day,hour twind=1.5 ! read in dimension of hr_airs, i.e, how many airs data files open(9000,file='number_of_airs_files.txt', form='formatted') read(9000,*) hr_airs close(9000) ALLOCATE(prof(leng)) ALLOCATE(prof_all(nrec,leng,hr_airs)) ALLOCATE(airs_file_path_and_name(hr_airs)) write(*,*) hr_airs ! open airs sfov sounding file name list open(9000,file='InputFileNames_airs.path',form='formatted') do L=1,hr_airs read(9000,*) airs_file_path_and_name(L) enddo close(9000) iuf = 90 do L=1,hr_airs open(iuf,file=airs_file_path_and_name(L),recl=leng,access='direct', & status='old',convert='big_endian') write(*,*) airs_file_path_and_name(L) do 1000 nr=1,nrec read(iuf,rec=nr) prof prof_all(nr,:,L) = prof(:) 1000 continue close(iuf) enddo ! read in pressure open(unit=20,file='pres.dat') do i=1,101 read(20,*) pres(i) enddo close(20) !! read in one day prepbufr files. ! open prepbufr file name list open(9000,file='InputFileNames_prepbufr.path',form='formatted') do L=1,hr_bufr read(9000,*) prepbufr_file_path_and_name(L) enddo close(9000) ! this is the first loop , loop over bufr files in one day DO L=1,hr_bufr write(*,*) L ! perform data selection, hour by hour, per prepbufr file ! read in year, month, day, hour (sring) year_c =prepbufr_file_path_and_name(L)(91:94) month_c=prepbufr_file_path_and_name(L)(95:96) day_c =prepbufr_file_path_and_name(L)(97:98) hour_c =prepbufr_file_path_and_name(L)(107:108) !convert string to real value read(year_c,*) year read(month_c,*) month read(day_c,*) day read(hour_c,*) hour write(*,*) year write(*,*) month write(*,*) day write(*,*) hour !! read in analysis time iadate(1)=year iadate(2)=month iadate(3)=day iadate(4)=hour iadate(5)=0.0 ! get gstime relative analysis time call w3fs21(iadate,gstime) ! write(*,*) gstime ! write(*,*) iadate ! first get the total number of observations num=0 do L2=1,hr_airs ! loop over airs data files do i=1,nrec !!just select the data with clear sky scenes ! if (prof_all(I,50,L2) .NE. -9999.0 ) then if ( prof_all(I,332,L2) .NE. -1. .and. prof_all(I,332,L2) .NE. 1 & .and. prof_all(I,333,L2) .NE. -1. .and. prof_all(I,333,L2) .NE. 1 & .and. prof_all(I,334,L2) .NE. -1. .and. prof_all(I,334,L2) .NE. 1 & .and. prof_all(I,335,L2) .NE. -1. .and. prof_all(I,335,L2) .NE. 1 & .and. prof_all(I,336,L2) .NE. -1. .and. prof_all(I,336,L2) .NE. 1 & .and. prof_all(I,337,L2) .NE. -1. .and. prof_all(I,337,L2) .NE. 4 ) then idate5(1)= prof_all(I,322,L2) !year idate5(2)= prof_all(I,323,L2) !month idate5(3)= prof_all(I,324,L2) !day idate5(4)= prof_all(I,325,L2) !hour idate5(5)= prof_all(I,326,L2) !minute call w3fs21(idate5,nmind) sstime = nmind + prof_all(I,327,L2)/60. tdiff=(sstime-gstime)/60. !! write(*,*) tdiff !! select data in the time window if (abs(tdiff) < twind ) then num=num+1 endif endif enddo enddo ! if there is obervation in the time window, then append them to the prepbufr file if ( num .GT. 0 ) then allocate(selected_prof_all(num,leng)) write(*,*) 'total selected observations =', num num=0 do L2=1,hr_airs ! loop over airs data files do i=1,nrec !!just select the data with clear sky scenes and the best QC marks (could be changed if using more data) if ( prof_all(I,332,L2) .NE. -1. .and. prof_all(I,332,L2) .NE. 1 & .and. prof_all(I,333,L2) .NE. -1. .and. prof_all(I,333,L2) .NE. 1 & .and. prof_all(I,334,L2) .NE. -1. .and. prof_all(I,334,L2) .NE. 1 & .and. prof_all(I,335,L2) .NE. -1. .and. prof_all(I,335,L2) .NE. 1 & .and. prof_all(I,336,L2) .NE. -1. .and. prof_all(I,336,L2) .NE. 1 & .and. prof_all(I,337,L2) .NE. -1. .and. prof_all(I,337,L2) .NE. 4 ) then idate5(1)= prof_all(I,322,L2) !year idate5(2)= prof_all(I,323,L2) !month idate5(3)= prof_all(I,324,L2) !day idate5(4)= prof_all(I,325,L2) !hour idate5(5)= prof_all(I,326,L2) !minute call w3fs21(idate5,nmind) sstime = nmind + prof_all(I,327,L2)/60. tdiff=(sstime-gstime)/60. !! write(*,*) tdiff !! select data in the time window if (abs(tdiff) < twind ) then num=num+1 selected_prof_all(num,:) = prof_all(I,:,L2) endif endif enddo enddo write(*,*) 'total seleted observations =', num open(unit=100,file='name_of_append_prepbufr_files.txt',form='formatted',access='append') write(100,'(A)') prepbufr_file_path_and_name(L) close(100) do i=1,num do j=1,leng if ( selected_prof_all(i,j) .EQ. -9999. ) then selected_prof_all(i,j) = 10.0e10 endif enddo enddo ! get bufr table from existing bufr file open(unit_table,file='prepbufr.table') ! open(unit_out,file='prepbufr',status='old',form='unformatted') open(unit_out,file=prepbufr_file_path_and_name(L),status='old',form='unformatted') write(*,*) prepbufr_file_path_and_name(L) call openbf(unit_out,'IN',unit_out) call dxdump(unit_out,unit_table) call closbf(unit_out) close(unit_out) ! !just print some sample data ! do i=1,101 ! write(*,*) pres(i),selected_prof_all(1,2+i)-273.15 ! enddo ! do i=1,101 ! write(*,*) pres(i),selected_prof_all(1,103+i)*1000. ! enddo ! ! apend observation into prepbufr file ! ! open(unit_out,file='prepbufr',status='old',form='unformatted') open(unit_out,file=prepbufr_file_path_and_name(L),status='old',form='unformatted') ! write(*,*) prepbufr_file_path_and_name(L) call datelen(10) call openbf(unit_out,'APN',unit_table) ! idate=2010050906 ! cycle time: YYYYMMDDHH idate_c= year_c//month_c//day_c//hour_c read(idate_c,*) idate write(*,*) idate subset='SATEMP' ! satellite derived t/q reports num_message = INT(num/1000.) + 1 write(*,*) 'number of messages=',num_message ! do ii=1, num_message call openmb(unit_out,subset,idate) do i=1,num ! do i= (ii-1)*1000+1,ii*1000 ! if ( i .gt. num) then ! write(*,*) 'exit, i=', i ! exit ! endif write(*,*) i ! call openmb(unit_out,subset,idate) ! get offset time idate5(1)= selected_prof_all(I,322) !year idate5(2)= selected_prof_all(I,323) !month idate5(3)= selected_prof_all(I,324) !day idate5(4)= selected_prof_all(I,325) !hour idate5(5)= selected_prof_all(I,326) !minute call w3fs21(idate5,nmind) sstime = nmind + selected_prof_all(I,327)/60. tdiff=(sstime-gstime)/60. ! set headers hdr=10.0e10 c_sid='airssfov'; hdr(1)=rstation_id hdr(2)=selected_prof_all(i,2); hdr(3)=selected_prof_all(i,1); hdr(4)=tdiff; hdr(6)=934.0; hdr(7)=255.0 ! set obs, qcf, oer for AIRS SFOV sounding hdr(5)=140 ! report typ obs=10.0e10;qcf=10.0e10;oer=10.0e10 kk=0 do k=1,90 if ( pres(101-k+1) > 50 .and. pres(101-k+1) < 1000 .and. selected_prof_all(i,103-k+1) .NE. 10.0e10 ) then ! just choose the data in 100 mb - 800 mb period or other period ) kk=kk+1 obs(1,kk) = pres(101-k+1) ! if ( selected_prof_all(i,103-k+1) .EQ. 10.0e10 ) then ! write (*,*) 'haidao missing data' ! write (*,*) k ! write(*,*) selected_prof_all(i,103-k+1) ! obs(1,kk) = 10.0e10 ! else ! obs(1,kk) = pres(101-k+1) ! flip over, so from bottom to top ! endif ! for water vapor, convert mixing ratio to specific humidity obs(2,kk)= selected_prof_all(i,204-k+1)/(1+selected_prof_all(i,204-k+1))*1000.;obs(3,kk)=selected_prof_all(i,103-k+1)-273.15 qcf(1,kk)=2.0 ;qcf(2,kk)=2.0 ;qcf(3,kk)=2.0 oer(1,kk)=0.5 ;oer(2,kk)=0.6 ;oer(3,kk)=1.5 endif enddo write(*,*) kk mxlv2=kk write(*,*) 'mxlv2', mxlv2 ! append temperature and moisture call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,mxlv2,iret,obstr) call ufbint(unit_out,oer,mxmn,mxlv2,iret,oestr) call ufbint(unit_out,qcf,mxmn,mxlv2,iret,qcstr) call writsb(unit_out) ! call closmg(unit_out) enddo !doloop for append each obs to prepbufr file call closmg(unit_out) call closbf(unit_out) !close bf close(unit_out) !close bufr file close(unit_table) ! close bufr table deallocate(selected_prof_all) endif ! if num .GT. 0 enddo ! doloop for hr_bufr, the biggest doloop end program