basic/bufr_append_sample.f900000644001460100027340000000201211631730152014522 0ustar mhurap4program bufr_append_sample ! ! sample of appending one observation into bufr file ! implicit none character(80):: hdstr='XOB YOB DHR' character(80):: obstr='TOB' real(8) :: hdr(3),obs(1,1) character(8) subset integer :: unit_out=10,unit_table=20 integer :: idate,iret ! set data values hdr(1)=85.0;hdr(2)=50.0;hdr(3)=0.2 obs(1,1)=300.0 idate=2008120101 ! YYYYMMDDHH subset='ADPSFC' ! surface land reports ! get bufr table from existing bufr file open(unit_table,file='prepobs_prep_app.bufrtable') open(unit_out,file='sample.bufr',status='old',form='unformatted') call openbf(unit_out,'IN',unit_out) call dxdump(unit_out,unit_table) call closbf(unit_out) ! append open(unit_out,file='sample.bufr',status='old',form='unformatted') call datelen(10) call openbf(unit_out,'APN',unit_table) call openmb(unit_out,subset,idate) call ufbint(unit_out,hdr,3,1,iret,hdstr) call ufbint(unit_out,obs,1,1,iret,obstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program basic/bufr_decode_sample.f900000644001460100027340000000164411631730401014505 0ustar mhurap4program bufr_decode_sample ! ! example of reading observations from bufr ! implicit none character(80):: hdstr='XOB YOB DHR' character(80):: obstr='TOB' real(8) :: hdr(3),obs(1,10) integer :: ireadmg,ireadsb character(8) subset integer :: unit_in=10 integer :: idate,iret,num_message,num_subset ! decode open(unit_in,file='sample.bufr',action='read',form='unformatted') call openbf(unit_in,'IN',unit_in) call datelen(10) num_message=0 msg_report: do while (ireadmg(unit_in,subset,idate) == 0) num_message=num_message+1 num_subset = 0 write(*,'(I10,I4,a10)') idate,num_message,subset sb_report: do while (ireadsb(unit_in) == 0) num_subset = num_subset+1 call ufbint(unit_in,hdr,3,1 ,iret,hdstr) call ufbint(unit_in,obs,1,10,iret,obstr) write(*,'(2I5,4f8.1)') num_subset,iret,hdr,obs(1,1) enddo sb_report enddo msg_report call closbf(unit_in) end program basic/bufr_encode_sample.f900000644001460100027340000000150411631521447014522 0ustar mhurap4program bufr_encode_sample ! ! example of writing one value into a bufr file ! implicit none character(80):: hdstr='XOB YOB DHR' character(80):: obstr='TOB' real(8) :: hdr(3),obs(1,1) character(8) subset integer :: unit_out=10,unit_table=20 integer :: idate,iret ! set data values hdr(1)=75.;hdr(2)=30.;hdr(3)=-0.1 obs(1,1)=287.15 idate=2008120100 ! YYYYMMDDHH subset='ADPUPA' ! upper-air reports ! encode open(unit_table,file='table_prepbufr.txt') open(unit_out,file='sample.bufr',action='write' & ,form='unformatted') call datelen(10) call openbf(unit_out,'OUT',unit_table) call openmb(unit_out,subset,idate) call ufbint(unit_out,hdr,3,1,iret,hdstr) call ufbint(unit_out,obs,1,1,iret,obstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program basic/decode_temperature.f900000644001460100027340000000054411660514103014542 0ustar mhurap4program bufr_decode_temperature implicit none real(8) :: obs character(8) subset integer :: idate,iret ! decode open(10,file=' t.bufr',action='read',form='unformatted') call openbf(10,'IN',10) call readmg(10,subset,idate,iret) call readsb(10,iret) call ufbint(10,obs,1,10,iret,'TOB') write(*,*) obs call closbf(10) end program basic/encode_temperature.f900000644001460100027340000000056011660514103014552 0ustar mhurap4program bufr_encode_temperature implicit none real(8) :: obs integer :: iret obs=10.15 ! encode open(20,file='bufrtable.txt') open(10,file='t.bufr',action='write',form='unformatted') call openbf(10,'OUT',20) call openmb(10, 'ADPUPA', 08120100) call ufbint(10,obs,1,1,iret,'TOB') call writsb(10) call closmg(10) call closbf(10) end program convention/prepbufr_append_retrieve.f900000644001460100027340000000376211555115555017104 0ustar mhurap4program prepbufr_append_retrieve ! ! append a retrieved data into prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=200 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) ! get bufr table from existing bufr file open(unit_table,file='prepobs_prep_app.bufrtable') open(unit_out,file='prepbufr',status='old',form='unformatted') call openbf(unit_out,'IN',unit_out) call dxdump(unit_out,unit_table) call closbf(unit_out) ! ! write observation into prepbufr file ! open(unit_out,file='prepbufr',status='old',form='unformatted') call datelen(10) call openbf(unit_out,'APN',unit_table) idate=2010050700 ! cycle time: YYYYMMDDHH subset='SATWND' ! upper-air (raob, drops) reports call openmb(unit_out,subset,idate) ! set headers hdr=10.0e10 c_sid='A114424Z'; hdr(1)=rstation_id hdr(2)=199.6; hdr(3)=13.9; hdr(4)=-1.0; hdr(6)=934.0; hdr(7)=255.0 ! set obs, qcf, oer for wind hdr(5)=251 ! report type: NESDIS VISIBLE CLOUD DRIFT ! (ALL LEVELS) (GOES) - u, v obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=906.0;obs(4,1)=934.0;obs(5,1)=-11.4;obs(6,1)=-3.3;obs(8,1)=6.0 qcf(1,1)=2.0 ;qcf(4,1)=2.0 ;qcf(5,1)=1.0 oer(5,1)=3.8 nlvl=1 ! encode wind obs call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr) call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr) call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program convention/prepbufr_append_surface.f900000644001460100027340000000463111555115555016703 0ustar mhurap4program prepbufr_append_surface ! ! append a surface observation into prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=1 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 character(8) :: c_sid real(8) :: rstation_id equivalence(rstation_id,c_sid) ! get bufr table from existing bufr file open(unit_table,file='prepobs_prep_app.bufrtable') open(unit_out,file='prepbufr',status='old',form='unformatted') call openbf(unit_out,'IN',unit_out) call dxdump(unit_out,unit_table) call closbf(unit_out) ! ! write observation into prepbufr file ! open(unit_out,file='prepbufr',status='old',form='unformatted') call datelen(10) call openbf(unit_out,'APN',unit_table) idate=2010050700 ! cycle time: YYYYMMDDHH subset='ADPSFC' ! surface land (SYNOPTIC, METAR) reports call openmb(unit_out,subset,idate) ! set headers hdr=10.0e10 c_sid='72408'; hdr(1)=rstation_id hdr(2)=284.8; hdr(3)=39.9; hdr(4)=-0.1; hdr(6)=9.0 ! set obs, qcf, oer for wind hdr(5)=281 ! report type obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=1008.6; obs(5,1)=4.0; obs(6,1)=-4.7; obs(8,1)=6.0 qcf(1,1)=2.0 ; qcf(5,1)=2.0 oer(5,1)=1.6 ! encode wind obs call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,mxlv,iret,obstr) call ufbint(unit_out,oer,mxmn,mxlv,iret,oestr) call ufbint(unit_out,qcf,mxmn,mxlv,iret,qcstr) call writsb(unit_out) ! set obs, qcf, oer for temperature and moisture hdr(5)=181 ! report type obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=1008.6;obs(2,1)=4925.0;obs(3,1)=25.1;obs(4,1)=9.0;obs(8,1)=0.0 qcf(1,1)=2.0 ;qcf(2,1)=2.0 ;qcf(3,1)=2.0 ;qcf(4,1)=2.0 oer(1,1)=0.5 ;oer(2,1)=0.6 ;oer(3,1)=1.5 ! encode temperature and moisture call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,mxlv,iret,obstr) call ufbint(unit_out,oer,mxmn,mxlv,iret,oestr) call ufbint(unit_out,qcf,mxmn,mxlv,iret,qcstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program convention/prepbufr_append_upperair.f900000644001460100027340000000576211555115555017110 0ustar mhurap4program prepbufr_append_upperair ! ! append a upper air observation into prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=200 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) ! get bufr table from existing bufr file open(unit_table,file='prepobs_prep_app.bufrtable') open(unit_out,file='prepbufr',status='old',form='unformatted') call openbf(unit_out,'IN',unit_out) call dxdump(unit_out,unit_table) call closbf(unit_out) ! ! write observation into prepbufr file ! open(unit_out,file='prepbufr',status='old',form='unformatted') call datelen(10) call openbf(unit_out,'APN',unit_table) idate=2010050700 ! cycle time: YYYYMMDDHH subset='ADPUPA' ! upper-air (raob, drops) reports call openmb(unit_out,subset,idate) ! set headers hdr=10.0e10 c_sid='71823'; hdr(1)=rstation_id hdr(2)=286.3; hdr(3)=53.8; hdr(4)=0.0; hdr(6)=307.0 ! set obs, qcf, oer for wind hdr(5)=220 ! report type: sounding obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=500.0; obs(5,1)=-2.0;obs(6,1)=0.7 ;obs(8,1)=1.0 qcf(1,1)=2.0 ; qcf(5,1)=2.0 oer(5,1)=2.5 obs(1,2)=432.5; obs(4,2)=6401.;obs(5,2)=-7.1;obs(6,2)=-1.2;obs(8,2)=4.0 qcf(1,2)=2.0 ; qcf(4,2)=2.0 ;qcf(5,2)= 2.0 oer(5,2)=2.6 nlvl=2 ! encode wind obs call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr) call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr) call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr) call writsb(unit_out) ! set obs, qcf, oer for temperature and moisture hdr(5)=120 ! report type: sounding obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=825.0;obs(2,1)=3672.0;obs(3,1)=0.8 ; obs(8,1)=2.0 qcf(1,1)=2.0 ;qcf(2,1)=2.0 ;qcf(3,1)=2.0 ; oer(2,1)=1.2 ;oer(3,1)=1.3 obs(1,2)=700.0;obs(2,2)=1157.0;obs(3,2)=-7.3;obs(4,2)=2841.;obs(8,2)=1.0 qcf(1,2)=2.0 ;qcf(2,2)=2.0 ;qcf(3,2)=2.0 ;qcf(4,2)=2.0 oer(2,2)=1.4 ;oer(3,2)=1.0 obs(1,3)=623.0;obs(2,3)=254.0 ;obs(3,3)=-12.9; ;obs(8,3)=2.0 qcf(1,3)=2.0 ;qcf(2,3)=2.0 ;qcf(3,3)=2.0 oer(2,3)=1.5 ;oer(3,3)=1.0 nlvl=3 ! encode temperature and moisture call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr) call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr) call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program convention/prepbufr_decode_all.f900000644001460100027340000000336311555115555016000 0ustar mhurap4program prepbufr_decode_all ! ! read all observations out from prepbufr. ! read bufr table from prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=250 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) INTEGER :: ireadmg,ireadsb character(8) :: subset integer :: unit_in=10,unit_table=24,idate,nmsg,ntb character(8) :: c_sid real(8) :: rstation_id equivalence(rstation_id,c_sid) integer :: i,k,iret ! ! open(unit_table,file='prepobs_prep.bufrtable') open(unit_in,file='prepbufr',form='unformatted',status='old') call openbf(unit_in,'IN',unit_in) call dxdump(unit_in,unit_table) call datelen(10) nmsg=0 msg_report: do while (ireadmg(unit_in,subset,idate) == 0) nmsg=nmsg+1 ntb = 0 write(*,*) write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate sb_report: do while (ireadsb(unit_in) == 0) ntb = ntb+1 call ufbint(unit_in,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_in,obs,mxmn,mxlv,iret,obstr) call ufbint(unit_in,oer,mxmn,mxlv,iret,oestr) call ufbint(unit_in,qcf,mxmn,mxlv,iret,qcstr) rstation_id=hdr(1) write(*,*) write(*,'(2I10,a14,8f14.1)') ntb,iret,c_sid,(hdr(i),i=2,8) DO k=1,iret write(*,'(i3,a10,9f14.1)') k,'obs=',(obs(i,k),i=1,9) write(*,'(i3,a10,9f14.1)') k,'oer=',(oer(i,k),i=1,7) write(*,'(i3,a10,9f14.1)') k,'qcf=',(qcf(i,k),i=1,7) ENDDO enddo sb_report enddo msg_report call closbf(unit_in) end program convention/prepbufr_encode_surface.f900000644001460100027340000000433611555115555016673 0ustar mhurap4program prepbufr_encode_surface ! ! write a surface observation into prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=1 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 character(8) :: c_sid real(8) :: rstation_id equivalence(rstation_id,c_sid) ! ! write observation into prepbufr file ! open(unit_table,file='prepobs_prep.bufrtable',action='read') open(unit_out,file='prepbufr',action='write',form='unformatted') call datelen(10) call openbf(unit_out,'OUT',unit_table) idate=2010050700 ! cycle time: YYYYMMDDHH subset='ADPSFC' ! surface land (SYNOPTIC, METAR) reports call openmb(unit_out,subset,idate) ! set headers hdr=10.0e10 c_sid='KTKI'; hdr(1)=rstation_id hdr(2)=263.4; hdr(3)=33.2; hdr(4)=-0.1; hdr(6)=179.0 ! set obs, qcf, oer for wind hdr(5)=287 ! report type obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=985.2; obs(5,1)=-2.8; obs(6,1)=-7.7; obs(8,1)=6.0 qcf(1,1)=2.0 ; qcf(5,1)=2.0 oer(5,1)=1.6 ! encode wind obs call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,mxlv,iret,obstr) call ufbint(unit_out,oer,mxmn,mxlv,iret,oestr) call ufbint(unit_out,qcf,mxmn,mxlv,iret,qcstr) call writsb(unit_out) ! set obs, qcf, oer for temperature and moisture hdr(5)=187 ! report type obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=985.2;obs(2,1)=12968.0;obs(3,1)=31.3;obs(4,1)=179.0;obs(8,1)=0.0 qcf(1,1)=2.0 ;qcf(2,1)=2.0 ;qcf(3,1)=2.0 ;qcf(4,1)=2.0 oer(1,1)=0.5 ;oer(2,1)=0.6 ;oer(3,1)=2.3 ! encode temperature and moisture call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,mxlv,iret,obstr) call ufbint(unit_out,oer,mxmn,mxlv,iret,oestr) call ufbint(unit_out,qcf,mxmn,mxlv,iret,qcstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program convention/prepbufr_encode_upperair.f900000644001460100027340000000615711555115555017075 0ustar mhurap4program prepbufr_encode_upperair ! ! write a upper air observation into prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=200 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) ! ! write observation into prepbufr file ! open(unit_table,file='prepobs_prep_app.bufrtable',action='read') open(unit_out,file='prepbufr',action='write',form='unformatted') call datelen(10) call openbf(unit_out,'OUT',unit_table) idate=2010050700 ! cycle time: YYYYMMDDHH subset='ADPUPA' ! upper-air (raob, drops) reports call openmb(unit_out,subset,idate) ! set headers hdr=10.0e10 c_sid='72293'; hdr(1)=rstation_id hdr(2)=242.9; hdr(3)=32.9; hdr(4)=0.0; hdr(6)=134.0 ! set obs, qcf, oer for wind hdr(5)=220 ! report type: sounding obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=998.0; obs(5,1)=4.6 ;obs(6,1)=2.2 ;obs(8,1)=3.0 qcf(1,1)=2.0 ; qcf(5,1)=2.0 oer(5,1)=2.3 obs(1,2)=850.0; obs(5,2)=2.0 ;obs(6,2)=-1.7;obs(8,2)=1.0 qcf(1,2)=2.0 ; qcf(5,2)=2.0 oer(5,2)=2.6 obs(1,3)=700.0; obs(5,3)=12.1;obs(6,3)=-4.4;obs(8,3)=1.0 qcf(1,3)=2.0 ; qcf(5,3)=2.0 oer(5,3)=2.5 nlvl=3 ! encode wind obs call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr) call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr) call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr) call writsb(unit_out) ! set obs, qcf, oer for temperature and moisture hdr(5)=120 ! report type: sounding obs=10.0e10;qcf=10.0e10;oer=10.0e10 obs(1,1)=998.0;obs(2,1)=8112.0;obs(3,1)=22.3;obs(4,1)=134.0;obs(8,1)=0.0 qcf(1,1)=2.0 ;qcf(2,1)=2.0 ;qcf(3,1)=2.0 ;qcf(4,1)=2.0 oer(1,1)=0.7 ;oer(2,1)=0.7 ;oer(3,1)=1.4 obs(1,2)=925.0;obs(2,2)=6312.0;obs(3,2)=14.1;obs(4,2)=779.0;obs(8,2)=1.0 qcf(1,2)=2.0 ;qcf(2,2)=2.0 ;qcf(3,2)=2.0 ;qcf(4,2)=2.0 oer(2,2)=0.9 ;oer(3,2)=1.5 obs(1,3)=850.0;obs(2,3)=2161.0;obs(3,3)=14.8;obs(4,3)=1493.;obs(8,3)=1.0 qcf(1,3)=2.0 ;qcf(2,3)=2.0 ;qcf(3,3)=2.0 ;qcf(4,3)=2.0 oer(2,3)=1.1 ;oer(3,3)=1.4 obs(1,4)=700.0;obs(2,4)=2131.0;obs(3,4)=9.2 ;obs(4,4)=3118.;obs(8,4)=1.0 qcf(1,4)=2.0 ;qcf(2,4)=2.0 ;qcf(3,4)=2.0 ;qcf(4,4)=2.0 oer(2,4)=1.4 ;oer(3,4)=1.0 nlvl=4 ! encode temperature and moisture call ufbint(unit_out,hdr,mxmn,1 ,iret,hdstr) call ufbint(unit_out,obs,mxmn,nlvl,iret,obstr) call ufbint(unit_out,oer,mxmn,nlvl,iret,oestr) call ufbint(unit_out,qcf,mxmn,nlvl,iret,qcstr) call writsb(unit_out) call closmg(unit_out) call closbf(unit_out) end program GSI_BUFR_interface/read_airs.f900000755001460100027340000000462311660514141015041 0ustar mhurap4program read_airs ! . . . . ! subprogram: read_airs read bufr format airs data use kinds, only: r_kind,r_double,i_kind ! Number of channels for sensors in BUFR integer(i_kind),parameter :: n_airschan = 281 !--- 281 subset ch out of 2378 ch for AIRS integer(i_kind),parameter :: n_amsuchan = 15 integer(i_kind),parameter :: n_hsbchan = 4 integer(i_kind),parameter :: n_totchan = n_amsuchan+n_airschan+n_hsbchan+1 ! BUFR file sequencial number character(len=512) :: table_file integer(i_kind) :: lnbufr = 10 integer(i_kind) :: lnbufrtab = 11 integer(i_kind) :: irec,isub,next ! Variables for BUFR IO real(r_double),dimension(2) :: aquaspot real(r_double),dimension(12,3) :: allspot real(r_double),dimension(n_totchan) :: allchan character(len=8) :: subset character(len=80) :: allspotlist integer(i_kind) :: iret, ireadmg,ireadsb logical :: airstab character(len=80) :: infile infile='airsbufr' allspotlist='SIID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ FOVN' ! Open BUFR file open(lnbufr,file=infile,form='unformatted') ! Open BUFR table table_file = 'airs_bufr.table' ! make table file name inquire(file=table_file,exist=airstab) if (airstab) then open(lnbufrtab,file=trim(table_file)) call openbf(lnbufr,'IN',lnbufrtab) else call openbf(lnbufr,'IN',lnbufr) endif call datelen(10) ! Big loop to read data file next=0 do while(ireadmg(lnbufr,subset,idate)>=0) next=next+1 read_loop: do while (ireadsb(lnbufr)==0) ! Read AIRSSPOT , AMSUSPOT and HSBSPOT call ufbrep(lnbufr,allspot,12,3,iret,allspotlist) if(iret /= 3) cycle read_loop ! Check that number of airs channel equals n_airschan ! only done until they match for one record and ndata is updated ! Read AIRSCHAN or AMSUCHAN or HSBCHAN call ufbrep(lnbufr,allchan,1,n_totchan,iret,'TMBR') if( iret /= n_totchan)then write(6,*)'READ_AIRS: ### ERROR IN READING ', ' BUFR DATA:', & iret, ' CH DATA IS READ INSTEAD OF ',n_totchan cycle read_loop endif ! Read AQUASPOT call ufbint(lnbufr,aquaspot,2,1,iret,'SOZA SOLAZI') enddo read_loop enddo call closbf(lnbufr) ! Close bufr file end program read_airs GSI_BUFR_interface/read_bufrtovs.f900000755001460100027340000000570011660514141015752 0ustar mhurap4program read_bufrtovs ! subprogram: read_bufrtovs read bufr tovs 1b data !$$$ use kinds, only: r_kind,r_double,i_kind implicit none ! Declare passed variables character(len=20) :: infile ! Declare local parameters integer(i_kind),parameter:: n1bhdr=13 integer(i_kind),parameter:: n2bhdr=14 ! Declare local variables logical hirs,msu,amsua,amsub,mhs,hirs4,hirs3,hirs2,ssu character(14):: obstype character(14):: infile2 character(8) subset character(80) hdr1b,hdr2b integer(i_kind) ireadsb,ireadmg,next integer(i_kind) i,j,k,llll,llb,lll integer(i_kind) iret,idate,nchanl,n integer(i_kind) lnbufr real(r_double),allocatable,dimension(:):: data1b8 real(r_double),dimension(n1bhdr):: bfr1bhdr real(r_double),dimension(n2bhdr):: bfr2bhdr !************************************************************************** obstype='mhs' infile='bufrtovs' hirs2 = obstype == 'hirs2' hirs3 = obstype == 'hirs3' hirs4 = obstype == 'hirs4' hirs = hirs2 .or. hirs3 .or. hirs4 msu= obstype == 'msu' amsua= obstype == 'amsua' amsub= obstype == 'amsub' mhs = obstype == 'mhs' ssu = obstype == 'ssu' if ( hirs ) then nchanl=19 else if ( msu ) then nchanl=4 else if ( amsua ) then nchanl=15 else if ( amsub ) then nchanl=5 else if ( mhs ) then nchanl=5 else if ( ssu ) then nchanl=3 endif llb=1 lll=1 lnbufr=10 allocate(data1b8(nchanl)) ! Big loop over standard data feed and possible ears data do llll=llb,lll ! Set bufr subset names based on type of data to read ! Open unit to satellite bufr file infile2=infile ! if(llll == 2)then ! infile2=trim(infile)//'ears' ! if(amsua .and. kidsat >= 200 .and. kidsat <= 207)go to 500 ! end if ! Reopen unit to satellite bufr file call closbf(lnbufr) open(lnbufr,file=infile2,form='unformatted',status = 'old',err = 500) call openbf(lnbufr,'IN',lnbufr) ! Loop to read bufr file next=0 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0) next=next+1 write(*,*) subset,idate read_loop: do while (ireadsb(lnbufr)==0) ! Read header record. (llll=1 is normal feed, 2=EARS data) hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS' call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b) write(*,*) (bfr1bhdr(i),i=1,n1bhdr) hdr2b ='SAZA SOZA BEARAZ SOLAZI' call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b) write(*,*) (bfr2bhdr(i),i=1,4) ! Read data record. Increment data counter if(llll == 1)then call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') else call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') end if write(*,*) (data1b8(i),i=1,iret) ! End of bufr read loops enddo read_loop enddo read_subset call closbf(lnbufr) 500 continue end do end program read_bufrtovs GSI_BUFR_interface/read_gps.f900000755001460100027340000000661411660514141014676 0ustar mhurap4program read_gps ! subprogram: read_gps read in and reformat gps data use kinds, only: r_kind,i_kind,r_double implicit none ! Declare passed variables character(len=80) :: infile ! Declare local parameters integer(i_kind),parameter:: maxlevs=500 ! Declare local variables character(10) nemo character(80) hdr1a character,dimension(8):: subset integer(i_kind) lnbufr,i,k,m,maxobs,ireadmg,ireadsb,said,ptid integer(i_kind) idate integer(i_kind) iret,levs,levsr,nreps_ROSEQ1 integer(i_kind),parameter:: mxib=31 ! integer(i_kind) ibit(mxib),nib ! logical lone logical ref_obs real(r_kind) qfro integer(i_kind),parameter:: n1ahdr=10 real(r_double),dimension(n1ahdr):: bfr1ahdr real(r_double),dimension(50,maxlevs):: data1b real(r_double),dimension(50,maxlevs):: data2a real(r_double),dimension(maxlevs):: nreps_this_ROSEQ2 data lnbufr/10/ data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU' / data nemo /'QFRO'/ !*********************************************************************************** infile='gpsbufr' ref_obs=.true. ! Open file for input, then read bufr data open(lnbufr,file=infile,form='unformatted') call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) if (iret/=0) goto 1010 ! Big loop over the bufr file do while(ireadmg(lnbufr,subset,idate)==0) read_loop: do while(ireadsb(lnbufr)==0) ! Read/decode data in subset (profile) ! Extract header information call ufbint(lnbufr,bfr1ahdr,n1ahdr,1,iret,hdr1a) call ufbint(lnbufr,qfro,1,1,iret,nemo) ! if (said == 4) then ! Gras ! call upftbv(lnbufr,nemo,qfro,mxib,ibit,nib) ! ! if(lone) then ! write(6,*)'READ_GPS: bad profile said=',said,'ptid=',ptid,& ! ' SKIP this report' ! cycle read_loop ! endif ! endif ! Read bending angle information ! Get the number of occurences of sequence ROSEQ2 in this subset ! (will also be the number of replications of sequence ROSEQ1), nreps_ROSEQ1 ! Also determine the number of replications of sequence ROSEQ2 nested inside ! each replication of ROSEQ1, ! nreps_this_ROSEQ2(1:nreps_ROSEQ1) - currently = 3 frequencies (L1, L2, zero) call ufbint(lnbufr,nreps_this_ROSEQ2,1,maxlevs,nreps_ROSEQ1,'{ROSEQ2}') ! Store entire contents of ROSEQ1 sequence (including contents of nested ROSEQ2 sequence) ! in array data1b call ufbseq(lnbufr,data1b,50,maxlevs,levs,'ROSEQ1') if(levs.ne.nreps_ROSEQ1) then write(6,*) 'READ_GPS: **WARNING** said,ptid=',said,ptid,& ' mismatch between sequence of ROSEQ1 and ROSEQ2 occurence',levs,nreps_ROSEQ1, & ' SKIP this report' cycle read_loop endif ! Check we have the same number of levels for ref and bending angle ! when ref_obs on to get lat/lon information call ufbseq(lnbufr,data2a,50,maxlevs,levsr,'ROSEQ3') ! refractivity if ((ref_obs).and.(levs/=levsr)) then write(6,*) 'READ_GPS: **WARNING** said,ptid=',said,ptid,& ' with gps_bnd levs=',levs,& ' and gps_ref levsr=',levsr,& ' SKIP this report' cycle read_loop endif enddo read_loop ! subsets enddo ! messages ! Close unit to input file 1010 continue call closbf(lnbufr) end program read_gps GSI_BUFR_interface/read_prepbufr.f900000755001460100027340000001420511660514141015725 0ustar mhurap4program read_prepbufr use kinds, only: r_single,r_kind,r_double,i_kind implicit none character(len=80) :: infile,obstype ! Declare local variables character(40) drift,hdstr,qcstr,oestr,sststr,satqcstr,levstr,hdstr2 character(40) metarcldstr,geoscldstr,metarvisstr,metarwthstr character(80) obstr character(10) date character(8) subset character(8) prvstr,sprvstr integer(i_kind) ireadmg,ireadsb,icntpnt,icntpnt2,icount,iiout integer(i_kind) lunin,i,maxobs,j,idomsfc,itemp,it29 integer(i_kind) metarcldlevs,metarwthlevs integer(i_kind) k,nmsg,kx, nreal,idate,iret,ncsave,levs integer(i_kind) ntb real(r_double) vtcd real(r_double),dimension(8):: hdr real(r_double),dimension(8,255):: drfdat,qcmark,obserr real(r_double),dimension(9,255):: obsdat real(r_double),dimension(8,1):: sstdat real(r_double),dimension(2,10):: metarcld real(r_double),dimension(1,10):: metarwth real(r_double),dimension(1,1) :: metarvis real(r_double),dimension(4,1) :: geoscld real(r_double),dimension(1):: satqc real(r_double),dimension(1,1):: r_prvstg,r_sprvstg real(r_double),dimension(1,255):: levdat real(r_double),dimension(255,20):: tpc real(r_double),dimension(2,255,20):: tobaux ! data statements data hdstr /'SID XOB YOB DHR TYP ELV SAID T29'/ data hdstr2 /'TYP SAID T29 SID'/ data obstr /'POB QOB TOB ZOB UOB VOB PWO CAT PRSS' / data drift /'XDR YDR HRDR '/ data sststr /'MSST DBSS SST1 SSTQM SSTOE '/ data qcstr /'PQM QQM TQM ZQM WQM NUL PWQ '/ data oestr /'POE QOE TOE NUL WOE NUL PWE '/ data satqcstr /'QIFN'/ data prvstr /'PRVSTG'/ data sprvstr /'SPRVSTG'/ data levstr /'POB'/ data metarcldstr /'CLAM HOCB'/ ! cloud amount and cloud base height data metarwthstr /'PRWE'/ ! present weather data metarvisstr /'HOVI'/ ! visibility data geoscldstr /'CDTP TOCC GCDTT CDTP_QM'/ ! NESDIS cloud products: cloud top pressure, temperature,amount logical tob,qob,uvob,spdob,sstob,pwob,psob logical metarcldobs,geosctpobs logical driftl,convobs data lunin / 13 / ! Initialize variables infile='prepbufr' nreal=0 satqc=0.0_r_kind obstype='t' tob = obstype == 't' uvob = obstype == 'uv' spdob = obstype == 'spd' psob = obstype == 'ps' qob = obstype == 'q' pwob = obstype == 'pw' sstob = obstype == 'sst' metarcldobs = obstype == 'mta_cld' geosctpobs = obstype == 'gos_ctp' convobs = tob .or. uvob .or. spdob .or. qob !------------------------------------------------------------------------ ! Open, then read date from bufr data call closbf(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) maxobs=0 nmsg=0 ntb = 0 msg_report: do while (ireadmg(lunin,subset,idate) == 0) nmsg=nmsg+1 loop_report: do while (ireadsb(lunin) == 0) ntb = ntb+1 ! Extract type information call ufbint(lunin,hdr,4,1,iret,hdstr2) kx=hdr(1) ! For the satellite wind to get quality information and check if it will be used if( kx ==243 .or. kx == 253 .or. kx ==254 ) then call ufbint(lunin,satqc,1,1,iret,satqcstr) if(satqc(1) < 85.0_r_double) cycle loop_report ! QI w/o fcst (su's setup endif ! Save information for next read if(ncsave /= 0) then call ufbint(lunin,levdat,1,255,levs,levstr) maxobs=maxobs+max(1,levs) end if end do loop_report enddo msg_report if (nmsg==0) goto 900 write(6,*)'READ_PREPBUFR: messages/reports = ',nmsg,'/',ntb !------------------------------------------------------------------------ ! Obtain program code (VTCD) associated with "VIRTMP" step if(tob)call ufbqcd(lunin,'VIRTMP',vtcd) ! loop over convinfo file entries; operate on matches !DTC comment out the loop loop_convinfo because we want to read all typies !DTC loop_convinfo: do nx=1, ntread call closbf(lunin) open(lunin,file=infile,form='unformatted') call openbf(lunin,'IN',lunin) call datelen(10) ! Big loop over prepbufr file ntb = 0 nmsg = 0 icntpnt=0 icntpnt2=0 loop_msg: do while (ireadmg(lunin,subset,idate)== 0) loop_readsb: do while(ireadsb(lunin) == 0) ! use msg lookup table to decide which messages to skip ! use report id lookup table to only process matching reports ntb = ntb+1 ! Extract type, date, and location information call ufbint(lunin,hdr,8,1,iret,hdstr) ! Balloon drift information available for these data !DTC driftl=kx==120.or.kx==220.or.kx==221 ! Extract data information on levels call ufbint(lunin,obsdat,9,255,levs,obstr) call ufbint(lunin,qcmark,8,255,levs,qcstr) call ufbint(lunin,obserr,8,255,levs,oestr) if(sstob)then sstdat=1.e11_r_kind call ufbint(lunin,sstdat,8,1,levs,sststr) else if(metarcldobs) then metarcld=1.e11_r_kind metarwth=1.e11_r_kind metarvis=1.e11_r_kind call ufbint(lunin,metarcld,2,10,metarcldlevs,metarcldstr) call ufbint(lunin,metarwth,1,10,metarwthlevs,metarwthstr) call ufbint(lunin,metarvis,1,1,iret,metarvisstr) if(levs /= 1 ) then write(6,*) 'READ_PREPBUFR: error in Metar observations, levs sould be 1 !!!' stop 110 endif else if(geosctpobs) then geoscld=1.e11_r_kind call ufbint(lunin,geoscld,4,1,levs,geoscldstr) endif ! If temperature ob, extract information regarding virtual ! versus sensible temperature ! if(tob) then ! call ufbevn(lunin,tpc,1,255,20,levs,'TPC') ! if (.not. twodvar_regional .or. .not.tsensible) then ! else !peel back events to store sensible temp in case temp is virtual ! call ufbevn(lunin,tobaux,2,255,20,levs,'TOB TQM') ! end if ! end if end do loop_readsb ! ! End of bufr read loop enddo loop_msg ! Close unit to bufr file call closbf(lunin) ! Normal exit !DTC enddo loop_convinfo! loops over convinfo entry matches 900 continue close(lunin) end program read_prepbufr radiance/bufr_decode_radiance.f900000644001460100027340000000501411555115555015466 0ustar mhurap4program bufr_decode_radiance ! ! read all radaince observations out from bufr. ! read bufr table from prepbufr file ! implicit none integer, parameter :: mxmn=35, mxlv=250 character(80):: hdstr= & 'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS' character(80):: hdr2b='SAZA SOZA BEARAZ SOLAZI' character(80):: obstr='TMBR' ! character(80):: obstr='TMBRST' real(8) :: hdr(mxmn),hdr2(mxmn),obs(mxlv) INTEGER :: ireadmg,ireadsb character(8) :: subset integer :: unit_in=10,unit_table=24,idate,nmsg,ntb integer,parameter:: max_sat_type=20 integer :: nsat_type(max_sat_type),nsat_num(max_sat_type) integer :: i,k,iret,ksatid,nchanl,num_sat_type,ii ! nchanl=15 nsat_num=0 nsat_type=0 ! open(unit_table,file='radiance.bufrtable') open(unit_in,file='1bamua',form='unformatted',status='old') call openbf(unit_in,'IN',unit_in) call dxdump(unit_in,unit_table) call datelen(10) nmsg=0 ntb = 0 num_sat_type = 0 msg_report: do while (ireadmg(unit_in,subset,idate) == 0) nmsg=nmsg+1 write(*,*) write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate sb_report: do while (ireadsb(unit_in) == 0) ntb = ntb+1 call ufbint(unit_in,hdr ,mxmn,1 ,iret,hdstr) call ufbint(unit_in,hdr2,mxmn,1 ,iret,hdr2b) call ufbrep(unit_in,obs ,1 ,nchanl,iret,obstr) ksatid=nint(hdr(1)) ! write(*,*) ! write(*,'(2I10,I14,13f8.1,)') ntb,iret,ksatid,(hdr(i),i=3,10),hdr(13), & ! (hdr2(i),i=1,4) ! write(*,'(a10,15f7.1)') 'obs=',(obs(i),i=1,iret) ! satellite type inventory if(num_sat_type == 0 ) then num_sat_type=1 nsat_type(num_sat_type)=ksatid nsat_num(num_sat_type)= 1 else ii=0 DO i=1,num_sat_type if(nsat_type(i) == ksatid) ii=i ENDDO if( ii > 0 .and. ii <=num_sat_type) then nsat_num(ii)=nsat_num(ii) + 1 else num_sat_type=num_sat_type+1 if(num_sat_type > max_sat_type) then write(*,*) 'Error: too many satellite types' write(*,*) 'Need to increase max_sat_type' stop 1234 endif nsat_type(num_sat_type)=ksatid nsat_num(num_sat_type)=1 endif endif enddo sb_report enddo msg_report call closbf(unit_in) write(*,*) 'message=',nmsg,' subset=',ntb DO i=1,num_sat_type write(*,'(i4,2i10)') i,nsat_type(i),nsat_num(i) ENDDO end program