! Module m_bufr_rw ! ! Module to read and write observation data in bufr format. ! For ! Module to read and write observation data in bufr format. ! ! This is explicitly written to be compatible with datasets used at the GMAO ! for input into GSI (April 2008 versions) ! ! Initial Code provided by Meta Sienkiewicz (Jan. 2008) ! Code extensively modified by Ronald Errico (Jan. 2008) ! Height read/write for conventional data added by Ronald Errico (Sept. 2008) ! ! program history log ! 2008-02-xx errio - initial code ! 2008-09-18 yang - modify to include MSU use m_kinds implicit none ! private public read_check_idate public read_obs_test public read_write_obs_tq public read_write_obs_uv public read_write_obs_airs public read_write_obs_amsu_hirs ! contains ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! subroutine read_check_idate (idate,idatetime) ! ! Check that date on file and passed as argument agree ! The assumption here is that idate is of the form YYMMDDHH, but that ! idatetime is YYYYMMDDHH. ! integer, parameter :: ikind8=8 integer :: idate integer(ikind8) :: idatetime ! logical :: ldiff integer(ikind8) :: time1(4), time2(4) integer(ikind8) :: id integer :: i integer(ikind8), parameter :: fifty=50_ikind8 ! ! Unpack idate id=idate time1(1)=id/1000000 id=id-time1(1)*1000000 time1(2)=id/10000 id=id-time1(2)*10000 time1(3)=id/100 time1(4)=id-time1(3)*100 if (time1(1) < fifty) then time1(1)=2000+time1(1) else time1(1)=1900+time1(1) endif ! ! Unpack idatetime time2(1)=idatetime/1000000 id=idatetime-time2(1)*1000000 time2(2)=id/10000 id=id-time2(2)*10000 time2(3)=id/100 time2(4)=id-time2(3)*100 ! ! Compare values ldiff=.false. do i=1,4 if (time1(i) /= time2(i)) then ldiff=.true. endif enddo ! if (ldiff) then print *,' ' print *,'Error: unpacked idate and cdatetime differ:' print ('(a,i12)'),' idate on file =',idate print ('(a,i12)'),' cdatetime as argument =',idatetime stop endif ! end subroutine read_check_idate ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! Subroutine read_obs_test (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields, & nobs,obs_info,obs_levels,obs_values) ! ! Fill "fake" observation data for testing interpolation software. ! implicit none ! integer :: obs_ndim1 integer :: obs_ndim2 integer :: obs_ndim3 integer :: obs_ndim4 integer :: obs_nlevs integer :: obs_nfields integer :: nobs ! real(rkind3) :: obs_info(obs_ndim4) real(rkind3) :: obs_levels(obs_ndim1,obs_ndim3) real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! integer :: nt, nf, k, n, i real(rkind3) :: ft, ff, fk, fy, fx ! obs_info(5)=1. ! obs type obs_info(6)=1.e5 ! station pressure obs_info(7)=100. ! station elevation obs_nlevs=0 ! if obs not included in list below ! if (nobs==1) then obs_nlevs=1 obs_nfields=2 obs_info(2)=64. ! obs lon obs_info(3)=42. ! obs lat obs_info(4)=2. ! obs time obs_levels(1,1)=23000. ! p of obs obs_values(:,:)=nobs elseif (nobs==2) then obs_nlevs=3 obs_nfields=2 obs_info(2)=-120. ! obs lon obs_info(3)=-42. ! obs lat obs_info(4)=-2. ! obs time obs_levels(1,1)=1720000. obs_levels(2,1)=72000. obs_levels(3,1)=12000. obs_values(:,:)=nobs elseif (nobs==3) then obs_nlevs=1 obs_nfields=2 obs_info(2)=370. ! obs lon obs_info(3)=95. ! obs lat obs_info(4)=2. ! obs time obs_levels(1,1)=720. obs_values(:,:)=nobs elseif (nobs==4) then obs_nlevs=3 obs_nfields=2 obs_info(2)=0. ! obs lon obs_info(3)=-89. ! obs lat obs_info(4)=-2. ! obs time obs_levels(1,1)=82000. obs_levels(2,1)=52000. obs_levels(3,1)= 2200. obs_values(:,:)=nobs endif ! print *,' ' print *,' ' print ('(a,3i7,8f10.2)'),' OBS:',nobs,obs_nlevs,obs_nfields,obs_info if (obs_nlevs > 0) then do i=1,obs_nlevs print ('(i4,f10.2,15f12.5)'),i,obs_levels(i,1), & obs_values(i,1:obs_nfields) enddo endif ! end subroutine read_obs_test ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! Subroutine read_write_obs_tq (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierror,luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values, & l_z_present,lread) ! ! Read/write temperature, moisture, or pressure obs ! in NCEP .prepbufr (BUFR) format ! use m_parameters implicit none ! integer :: obs_ndim1 integer :: obs_ndim2 integer :: obs_ndim3 integer :: obs_ndim4 integer :: obs_nlevs integer :: obs_nfields integer :: nobs integer :: ierror integer :: luin integer :: luout ! logical :: ltest logical :: lread logical :: lsim ! indicates that these obs are to be simulated ! and thus read values can be over-written logical :: l_z_present ! indicates that heights are present with data ! real(rkind3) :: bmiss real(rkind3) :: obs_info(obs_ndim4) real(rkind3) :: obs_levels(obs_ndim1,obs_ndim3) real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! ! local variables logical, save :: l_z_same ! true if station elev=z obs for 1-level obs. logical, save :: l_z_zero ! true if sation elev =0 for 1-level obs. real(rkind3) :: bmiss_kind3 ! integer :: i ! loop counters integer :: ntq ! # of valid tq levels in a report integer :: nz ! # of non-missing z values real(8), save :: obsdat(13,255) ! obs profile from bufr file real(8), save :: drift(3,255) ! balloon drift real(8), save :: hdr(6) ! location and metadata real(8), save :: bout(4,255) ! output event real(8), save :: boutp(5,255) ! output event for pevstr integer, save :: indm(255) ! index of obs in vertical integer, save :: levtq, klev, levd, llev integer levdch character(len=52) hdstr, driftstr, tqdstr character(len=32) pevstr, tevstr, qevstr, zevstr data hdstr /'SID XOB YOB DHR TYP ELV ' / data driftstr /'XDR YDR HRDR ' / data tqdstr /'POB PQM ZOB ZQM TOB TQM QOB QQM CAT PPC ZPC TPC QPC '/ data pevstr /'POB PQM PPC PRC CAT ' / data tevstr /'TOB TQM TPC TRC ' / data qevstr /'QOB QQM QPC QRC ' / data zevstr /'ZOB ZQM ZPC ZRC ' / ! ! _OB indicates observation value, _QM indicates quality mark ! _RC indicates "reason code"; i.e., the reason for a particular quality mark ! _PC indicates "program code"; i.e., how obs values were computed ! Here, the _QM are copied from the original since (1) the quality-accepted obs ! count for the simulated data will then be similar to the original data ! and (2) any "missing data" created by the simulation program (e.g., ! below ground values in a raob report) will still be detectable as a ! missing value flag. ! The _PC values are copied from the original in case the DAS expects some ! particular values (in the GSI used at time of this development, _PC values ! are not used, but this can change...) ! real(rkind3), parameter :: c2kelvin=273.15 real(rkind3), parameter :: pscale=1.e2 real(rkind3), parameter :: qscale=1.e-6 ! ! bmiss_kind3=bmiss ! ! Check if this is a read or write if (lread) then nobs=nobs+1 ! ! Read information from the input file needed for interpolation call ufbint(luin,hdr,6,1,klev,hdstr) call ufbint(luin,obsdat,13,255,levtq,tqdstr) ! ! Check for levels with observations (non-missing pressure and ! one of T, q, and/or surface level). ! Also check if some observed heights are present ntq = 0 nz = 0 do i = 1,levtq if (obsdat(1,i)=180) .or. (nint(hdr(5))<=189) ) ) then l_z_zero=.true. else l_z_zero=.false. endif ! ! copy into main program arrays if (ntq > 0) then ! found tq info, continue processing do i=1,ntq ! if (obsdat(5,i) 300.) .and. & (obsdat(5,i)= 180) .and. (nint(hdr(5)) <= 189) ) ) ) then obs_levels(i,1)=-obs_levels(i,1) ! flag to indicate sfc. value endif ! if (l_z_present) then obs_levels(i,2)=obsdat(3,indm(i)) ! z obs value else obs_levels(i,2)=bmiss_kind3 endif enddo endif ! ! ntq=0 means no T,Q,Z data found obs_nlevs=ntq obs_nfields=2 ! 2 because both t and q fields obs_info(1)=hdr(1) ! id obs_info(2)=hdr(2) ! lon of obs obs_info(3)=hdr(3) ! lat of obs obs_info(4)=hdr(4) ! time of obs obs_info(5)=hdr(5) ! type of obs obs_info(6)=hdr(6) ! height of surface obs_info(7:obs_ndim4)=0. ! (not used) ! ! End of read of observation buffer ! else ! ! This is a call to write observation buffer ! Only write report if there are valid observations in it. ! if (obs_nlevs > 0) then ! ! Change surface elevation to NR elevation in both header and line of report do i=1,obs_nlevs if (nint(obsdat(9,indm(i))) == 0) then obsdat(3,indm(i))=obs_info(6) ! NR surface height hdr(6)=obs_info(6) ! NR surface height endif enddo ! ! For single-level reports: set to station elev if surface report, unless ! original value of station elev was 0, in which case the 0 value is preserved. ! Also, if station elev and obs z were originally the same, set the former to ! the latter computed from the nature run, unless both were originally 0. ! if (obs_nlevs==1) then if ((nint(hdr(5))>=181) .or. (nint(hdr(5))<=189) ) then if (l_z_zero) then hdr(6)=zero_k3 ! preserve 0 if original value was 0 obs_info(6)=zero_k3 else hdr(6)=obs_info(6) ! reset to sfc z in nature run endif endif if (l_z_same) then if (l_z_zero) then obs_levels(1,2)=zero_k3 ! hdr(6) already set to 0 above. else hdr(6)=obs_levels(1,2) ! e.g., for AIRCAR obs_info(6)=obs_levels(1,2) endif endif endif ! ! write out p event boutp = bmiss do i = 1,obs_nlevs if ( (obs_levels(i,1) == bmiss_kind3) ) then boutp(4,i) = 13. ! reason code else boutp(1,i) = obs_levels(i,1)/pscale boutp(2,i) = obsdat(2,indm(i)) ! copy quality mark boutp(4,i) = 0. ! reason code endif boutp(3,i) = obsdat(10,indm(i)) ! copy program code boutp(5,i) = obsdat(9,indm(i)) ! copy CAT enddo call ufbint(luout,boutp,5,obs_nlevs,llev,pevstr) if (llev .ne. obs_nlevs) then print *,'incompatible level number returned in write of pevstr' ierror=ierror+1 endif ! ! write out t event bout = bmiss do i = 1,obs_nlevs if (obs_values(i,1) 3) ) then ! write as T bout(1,i) = obs_values(i,1)-c2kelvin else ! write as Tv bout(1,i) = obs_values(i,1) * & (1. + ratio4R*obs_values(i,2)) - c2kelvin endif bout(2,i) = obsdat(6,indm(i)) ! copy quality mark bout(4,i) = 0. ! reason code else bout(4,i)=13. ! reason code endif bout(3,i) = obsdat(12,indm(i)) ! copy program code enddo call ufbint(luout,bout,4,obs_nlevs,llev,tevstr) if (llev .ne. obs_nlevs) then print *,'incompatible level number returned in write of tevstr' ierror=ierror+1 endif ! ! write out q event bout = bmiss do i = 1,obs_nlevs if ( (obs_values(i,2) 0) then do i=1,obs_nlevs print ('(i4,2f10.2,15f12.5)'),i,obs_levels(i,1:2), & obs_values(i,1:obs_nfields) enddo endif endif ! test on ltest ! end subroutine read_write_obs_tq ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! Subroutine read_write_obs_uv (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_nlevs,obs_nfields,nobs, & ierror,luin,luout,ltest,lsim,bmiss, & obs_info,obs_levels,obs_values, & l_z_present,lread) ! ! Read/write wind observations in NCEP .prepbufr (BUFR) format ! implicit none ! integer :: obs_ndim1 integer :: obs_ndim2 integer :: obs_ndim3 integer :: obs_ndim4 integer :: obs_nlevs integer :: obs_nfields integer :: nobs integer :: ierror integer :: luin integer :: luout ! logical :: ltest logical :: lread logical :: lsim ! indicates that these obs are to be simulated ! and thus read values can be over-written logical :: l_z_present ! indicates that heights are present with data ! real(8) :: bmiss real(rkind3) :: obs_info(obs_ndim4) real(rkind3) :: obs_levels(obs_ndim1,obs_ndim3) real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! ! local variables logical, save :: l_z_same ! true if station elev=z of 1-level observation logical, save :: l_z_zero ! true if sation elev =0 for 1-level obs. real(rkind3) :: bmiss_kind3 ! integer i ! loop counters integer nw ! # of wind levels in a report integer nz ! # of height values in a report real(8), save :: obsdat(11,255) ! obs profile from bufr file real(8), save :: drift(3,255) ! balloon drift real(8), save :: hdr(6) ! location and metadata real(8), save :: uvout(5,255) ! output wind array integer, save :: indw(255) ! index of wind obs in vertical real(8), save :: boutp(5,255) ! output for pevstr real(8), save :: boutz(4,255) ! output for zevstr integer, save :: levw, klev, levd, llev integer levdch real(rkind3), parameter :: pscale=1.e2 character(len=44) hdstr, driftstr, wndstr, uevstr, pevstr, zevstr data hdstr /'SID XOB YOB DHR TYP ELV ' / data driftstr /'XDR YDR HRDR ' / data wndstr /'POB PQM UOB VOB WQM CAT ZOB ZQM PPC WPC ZPC ' / data pevstr /'POB PQM PPC PRC CAT ' / data uevstr /'UOB VOB WQM WPC WRC ' / data zevstr /'ZOB ZQM ZPC ZRC ' / ! ! _OB indicates observation value, _QM indicates quality mark ! _RC indicates "reason code"; i.e., the reason for a particular quality mark ! _PC indicates "program code"; i.e., how obs values were computed ! Here, the _QM are copied from the original since (1) the quality-accepted obs ! count for the simulated data will then be similar to the original data ! and (2) any "missing data" created by the simulation program (e.g., ! below ground values in a raob report) will still be detectable as a ! missing value flag. ! The _PC values are copied from the original in case the DAS expects some ! particular values (in the GSI used at time of this development, _PC values ! are not used, but this can change...) ! ! bmiss_kind3=bmiss ! ! Check if this is a read or write if (lread) then nobs=nobs+1 ! Read information from the input file needed for interpolation call ufbint(luin,hdr,6,1,klev,hdstr) call ufbint(luin,obsdat,11,255,levw,wndstr) ! ! Check for levels with wind observations (and non-missing pressure) ! Also check if some observed heights are present nw = 0 nz = 0 do i = 1,levw if (obsdat(1,i)=280) .or. (nint(hdr(5))<=289) ) ) then l_z_zero=.true. else l_z_zero=.false. endif ! ! copy into main program arrays if (nw > 0) then ! found wind info, continue processing do i=1,nw obs_values(i,1)=obsdat(3,indw(i)) obs_values(i,2)=obsdat(4,indw(i)) ! ! If the intention is to simulate (i.e., over-write) observations and ! this obs value is either explicitly indicated as a surface value (e.g., ! in a multilevel RAOB report) or implied surface value in a surface-type ! report, flag the p value by making it negative so it will be changed later. obs_levels(i,1)=obsdat(1,indw(i))*pscale ! change hPa to Pa if (lsim .and. ( (nint(obsdat(6,indw(i))) == 0) .or. & ( ( nw == 1) .and. ( (nint(hdr(5)) == 281) .or. & ((nint(hdr(5)) >= 284) .and. (nint(hdr(5)) <= 289)) ) ) )) then obs_levels(i,1)=-obs_levels(i,1) ! flag to indicate surface value endif ! if (l_z_present) then obs_levels(i,2)=obsdat(7,indw(i)) else obs_levels(i,2)=bmiss_kind3 endif ! enddo endif ! ! nw=0 means no data of type wndstr found obs_nlevs=nw obs_nfields=2 ! 2 because both u and v fields obs_info(1)=hdr(1) ! id obs_info(2)=hdr(2) ! lon of obs obs_info(3)=hdr(3) ! lat of obs obs_info(4)=hdr(4) ! time of obs obs_info(5)=hdr(5) ! type of obs obs_info(6)=hdr(6) ! height of surface obs_info(7)=1000.*100 ! ps in pascals (not used for this data type) ** ! ! End of read of observation buffer ! else ! ! This is a call to write observation buffer ! ! Write the changes to wind obs into the output buffer ! Only write report if there are valid observations in it. ! if (obs_nlevs > 0) then ! ! First change station elevation as necessary to agree with specifications ! in the NCEP .prepbufr files if (obs_nlevs==1) then if ((nint(hdr(5))==281) .or. (nint(hdr(5))==287) ) then hdr(6)=obs_info(6) ! reset to sfc z in nature run elseif (l_z_zero) then hdr(6)=zero_k3 ! preserve 0 if original sfc value was 0 obs_info(6)=zero_k3 if (l_z_same) obs_levels(1,2)=zero_k3 endif if (l_z_same) then ! e.g., for AIRCAR, SATWND hdr(6)=obs_levels(1,2) obs_info(6)=obs_levels(1,2) endif if (nint(hdr(5))==285) then ! set to 10m for QKSWND hdr(6)=10._rkind3 obs_info(6)=10. endif endif ! ! Write out p event ! (note special setting for SSMI and QKSCAT Winds in NCEP .prebufr files) boutp = bmiss do i = 1,obs_nlevs if (obs_levels(i,1) == bmiss_kind3) then boutp(4,i) = 13. ! reason code elseif ((nint(hdr(5))==283) .or. (nint(hdr(5))==285)) then boutp(1,i) = 1013._rkind3 ! case for SPSSMI and QKSWND obs_levels(i,1)=101300. boutp(2,i) = obsdat(2,indw(i)) boutp(4,i) = 0. else boutp(1,i) = obs_levels(i,1)/pscale boutp(2,i) = obsdat(2,indw(i)) ! copy quality mark boutp(4,i) = 0. ! reason code endif boutp(3,i) = obsdat(9,indw(i)) ! copy program code boutp(5,i) = obsdat(6,indw(i)) ! copy CAT enddo call ufbint(luout,boutp,5,obs_nlevs,llev,pevstr) if (llev .ne. obs_nlevs) then print *,'incompatible level number returned in write of pevstr' ierror=ierror+1 endif ! ! Write out u and v event uvout = bmiss do i = 1,obs_nlevs if ( (obs_values(i,1) == bmiss_kind3) .or. & (obs_values(i,2) == bmiss_kind3) ) then uvout(5,i) = 13. ! reason code else uvout(1,i) = obs_values(i,1) uvout(2,i) = obs_values(i,2) uvout(3,i) = obsdat(5,indw(i)) ! copy quality mark uvout(5,i) = 0. ! reason code endif uvout(4,i) = obsdat(10,indw(i)) ! copy program code enddo ! call ufbint(luout,uvout,5,obs_nlevs,llev,uevstr) if (llev .ne. obs_nlevs) then print *,'incompatible level number returned in write of uvout' ierror=ierror+1 endif ! ! write out z event if (l_z_present) then boutz = bmiss do i = 1,obs_nlevs if (obs_levels(i,2) == bmiss_kind3) then boutz(4,i) = 13. ! reason code else boutz(1,i) = obs_levels(i,2) boutz(2,i) = obsdat(8,indw(i)) ! copy quality mark boutz(4,i) = 0. ! reason code endif boutz(3,i) = obsdat(11,indw(i)) ! copy program code enddo call ufbint(luout,boutz,4,obs_nlevs,llev,zevstr) if (llev .ne. obs_nlevs) then print *,'incompatible level number returned in write of zevstr' ierror=ierror+1 endif endif ! check on l_z_present ! ! reset drift information for radiosonde and pibal obs (120, 220, 221) if (hdr(5)==120..or.hdr(5)==220..or.hdr(5)==221.) then call ufbint(luin,drift,3,255,levd,driftstr) if (levd .ne. levw) then print *,'incompatible level number returned in read of drift' ierror=ierror+1 endif ! Here we (re)set the drift information on the data levels to match the ! location/time in the header. (If we were using the 'drift' we would ! need to index by 'indw'.) do i = 1,obs_nlevs drift(1,i)=hdr(2) drift(2,i)=hdr(3) drift(3,i)=hdr(4) enddo ! call ufbint(luout,drift,3,obs_nlevs,levdch,driftstr) if (obs_nlevs .ne. levdch) then print *,'incompatible level number returned in write of drift' ierror=ierror+1 endif endif ! test on obs_type hdr(5) ! Finally, write out the header information. call ufbint(luout,hdr,6,1,klev,hdstr) ! ! Finished updates; so write output buffer to output file for this report ! call writsb(luout) ! endif ! test on obs_nlevs ! endif ! test on whether to read or write observation buffer ! ! Optionally print some info if testing mode if (ltest) then print *,' ' if (lread) then print *,' Read observation buffer' else print *,' Write observation buffer' endif print ('(a,3i7,8f10.2)'),' OBS:',nobs,obs_nlevs,obs_nfields,obs_info if (obs_nlevs > 0) then do i=1,obs_nlevs print ('(i4,2f10.2,15f12.5)'),i,obs_levels(i,1:2), & obs_values(i,1:obs_nfields) enddo endif endif ! test on ltest ! end subroutine read_write_obs_uv ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! Subroutine read_write_obs_airs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,nobs, & ierror,luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,lread) ! ! Read/write bufr file of AIRS obs ! Each record is assumed to contain possible information for ! AIRS, AMSUA, and HSB data, although the latter is missing. ! use m_kinds ! precision specification module for real variables implicit none ! integer :: obs_ndim1 integer :: obs_ndim2 integer :: obs_ndim3 integer :: obs_ndim4 integer :: obs_n_channels integer :: obs_n_data integer :: nobs integer :: ierror integer :: luin integer :: luout ! logical :: ltest logical :: lread logical :: lsim ! not currently used in this subroutine ! real(8) :: bmiss real(rkind3) :: obs_info(obs_ndim4) real(rkind3) :: obs_channels(obs_ndim1,obs_ndim3) real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! ! local variables ! integer, parameter :: n_airs_channels=281 integer, parameter :: n_amsu_channels=15 integer, parameter :: n_hsb_channels=5 integer, parameter :: n_aquaspot=25 integer, parameter :: ndval=4 ! number of data values per channel integer, parameter :: nhead=12 ! number of header values integer, parameter :: max_ch=305 ! maximum number of channels integer :: i, i1, i2 real(8), save :: aquaspot(n_aquaspot) ! aqua sat. header information real(8), save :: hdr(nhead,3) ! instrument header information real(8), save :: obsdat(ndval,max_ch) ! obs profile from bufr file real(8), save :: bout(ndval,max_ch) ! output event real(rkind3) :: bmiss3 integer :: khead, kchan, khead_sum, kchan_sum character(len=70) :: hdstr ! header string character(len=32) :: chstr ! channel string data hdstr /'SIID CLATH CLONH YEAR MNTH DAYS HOUR MINU SECO SAZA FOVN BEARAZ' / data chstr /'CHNM LOGRCW ACQF TMBR '/ ! ! In the hearder string, clat and clon should always be the 2nd and 3rd values ! In the hearder string, YEAR etc. should always be values 4-9 ! ! Check if this is a read or write if (lread) then nobs=nobs+1 ! Read header information for 1 report ! First read "AQUASPOT" (header) information (This is not used in ! the simulation software, but is saved to maintain the proper header) call ufbseq(luin,aquaspot,n_aquaspot,1,khead,'AQUASPOT') ! ! read header information for 1 report call ufbrep(luin,hdr,nhead,3,khead,hdstr) if (khead /= 3) then print *,' ' print *,' Error: khead=',khead print *,'AIRS file does not have 3 instruments in each report' stop endif ! ! Copy header and "spot" information obs_n_data=2 obs_info(1:nhead)=hdr(1:nhead,1) obs_info(nhead+1:2*nhead)=hdr(1:nhead,2) obs_info(2*nhead+1:3*nhead)=hdr(1:nhead,3) obs_info(3*nhead+1:3*nhead+n_aquaspot)=aquaspot(1:n_aquaspot) ! ! Read channels for AIRS call ufbseq(luin,obsdat,ndval,n_airs_channels,kchan,'AIRSCHAN') kchan_sum=kchan do i=1,n_airs_channels obs_values(i,1)=obsdat(4,i) obs_values(i,2)=obsdat(3,i) obs_channels(i,1)=obsdat(1,i) obs_channels(i,2)=obsdat(2,i) enddo ! ! Read channels for AMSU-A call ufbseq(luin,obsdat,ndval,n_amsu_channels,kchan,'AMSUCHAN') kchan_sum=kchan_sum+kchan i1=n_airs_channels do i=1,n_amsu_channels i1=i1+1 obs_values(i1,1)=obsdat(4,i) obs_values(i1,2)=obsdat(3,i) obs_channels(i1,1)=obsdat(1,i) obs_channels(i1,2)=obsdat(2,i) enddo ! ! Read channels for HSB call ufbseq(luin,obsdat,ndval,n_hsb_channels,kchan,'HSBCHAN') kchan_sum=kchan_sum+kchan do i=1,n_hsb_channels i1=i1+1 obs_values(i1,1)=obsdat(4,i) obs_values(i1,2)=obsdat(3,i) obs_channels(i1,1)=obsdat(1,i) obs_channels(i1,2)=obsdat(2,i) enddo ! obs_n_channels=i1 ! ! Check that all required data has been read if (kchan_sum /= obs_n_channels) then print *,' ' print *,' Error in read of airs obs:' print *,' obs_n_channels=',obs_n_channels,' /= ',kchan_sum, & '=kchan_sum' print *,'n_airs_channels=',n_airs_channels print *,'n_amsu_channels=',n_amsu_channels print *,'n_hsb_channels=',n_hsb_channels stop endif ! ! End of read of observation buffer ! else ! ! This is a call to write observation buffer ! if (obs_n_channels > 0) then bmiss3=bmiss/1.1 ! ! Write out AIRS information ! First, copy header to correct locations in 'SPOT' sequences do i=1,3 i1=(i-1)*nhead hdr(1,i)=obs_info(1+i1) do i2=2,7 hdr(i2,i)=obs_info(2+i2+i1) enddo hdr( 8,i)=obs_info( 2+i1) hdr( 9,i)=obs_info( 3+i1) hdr(10,i)=obs_info(10+i1) hdr(11,i)=obs_info(12+i1) hdr(12,i)=obs_info(11+i1) enddo aquaspot(1:n_aquaspot)=obs_info(3*nhead+1:3*nhead+n_aquaspot) ! ! Write out "AQUASPOT" information call ufbseq(luout,aquaspot,n_aquaspot,1,khead,'AQUASPOT') ! ! Write out AIRS information do i=1,n_airs_channels if (obs_values (i,1) < bmiss3) then bout(4,i)=obs_values(i,1) else bout(4,i)=bmiss endif bout(3,i)=obs_values(i,2) bout(1,i)=obs_channels(i,1) bout(2,i)=obs_channels(i,2) enddo call drfini(luout,n_airs_channels, 1,'(AIRSCHAN)') call ufbseq(luout,hdr(1:nhead,1),nhead,1,khead,'AIRSSPOT') call ufbseq(luout,bout,ndval,n_airs_channels,kchan,'AIRSCHAN') kchan_sum=kchan ! ! Write out AMSU information i1=n_airs_channels do i=1,n_amsu_channels i1=i1+1 if (obs_values (i1,1) < bmiss3) then bout(4,i)=obs_values(i1,1) else bout(4,i)=bmiss endif bout(3,i)=obs_values(i1,2) bout(1,i)=obs_channels(i1,1) bout(2,i)=obs_channels(i1,2) enddo call ufbseq(luout,hdr(1:nhead,2),nhead,1,khead,'AMSUSPOT') call ufbseq(luout,bout,ndval,n_amsu_channels,kchan,'AMSUCHAN') kchan_sum=kchan+kchan_sum ! ! Write out HSB information do i=1,n_hsb_channels i1=i1+1 if (obs_values (i1,1) < bmiss3) then bout(4,i)=obs_values(i1,1) else bout(4,i)=bmiss endif bout(3,i)=obs_values(i1,2) bout(1,i)=obs_channels(i1,1) bout(2,i)=obs_channels(i1,2) enddo call ufbseq(luout,hdr(1:nhead,3),nhead,1,khead,'HSBSPOT') call ufbseq(luout,bout,ndval,n_hsb_channels,kchan,'HSBCHAN') kchan_sum=kchan+kchan_sum ! ! Check that all required data has been read if (kchan_sum /= obs_n_channels) then print *,' ' print *,' Error in write of airs observation report:' print *,' obs_n_channels=',obs_n_channels,' and kchan_sum=', & kchan_sum,' and i1=',i1 print *,'n_airs_channels=',n_airs_channels print *,'n_amsu_channels=',n_amsu_channels print *,'n_hsb_channels=',n_hsb_channels stop endif ! ! finished updates; so write output buffer to output file for this report ! call writsb(luout) ! endif ! test on obs_n_channels ! endif ! test on whether to read or write observation buffer ! ! Optionally print some info if testing mode if (ltest) then print *,' ' if (lread) then print *,' Read observation buffer:' else print *,' Write observation buffer:' endif print ('(a,2i8)'),' OBS number, channels:',nobs,obs_n_channels print ('(a,10f10.3)'),' OBS AQUASPOT( 1:10):',aquaspot( 1:10) print ('(a,10f10.3)'),' OBS AQUASPOT(11:20):',aquaspot(11:20) print ('(a,10f10.3)'),' OBS AQUASPOT(21:25):',aquaspot(21:25) print ('(a,12f10.3)'),' OBS HDR INFO AIRS:',obs_info(1:nhead) print ('(a,12f10.3)'),' OBS HDR INFO AMSU:', & obs_info(nhead+1:2*nhead) print ('(a,12f10.3)'),' OBS HDR INFO HSB:', & obs_info(2*nhead+1:3*nhead) if (obs_n_channels > 0) then do i=1,obs_n_channels print ('(i4,4f15.5)'),i,obs_channels(i,1:2),obs_values(i,1:2) enddo endif endif ! test on ltest ! end subroutine read_write_obs_airs ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! Subroutine read_write_obs_amsu_hirs (obs_ndim1,obs_ndim2,obs_ndim3, & obs_ndim4,obs_n_channels,obs_n_data,nobs, & ierror,luin,luout,ltest,lsim,bmiss, & obs_info,obs_channels,obs_values,lread) ! ! Read/write bufr file of HIRS, AMSU,and MSU obs ! implicit none ! integer :: obs_ndim1 integer :: obs_ndim2 integer :: obs_ndim3 integer :: obs_ndim4 integer :: obs_n_channels integer :: obs_n_data integer :: nobs integer :: ierror integer :: luin integer :: luout ! logical :: ltest logical :: lread logical :: lsim ! not currently used in this subroutine ! real(8) :: bmiss real(rkind3) :: obs_info(obs_ndim4) real(rkind3) :: obs_channels(obs_ndim1,obs_ndim3) real(rkind3) :: obs_values(obs_ndim1,obs_ndim2) ! ! local variables ! integer, parameter :: ndval=2 ! number of data values per channel integer, parameter :: nhead=16 ! number of header values integer, parameter :: max_ch=30 ! maximum number of channels integer :: i real(8), save :: obsdat(ndval,max_ch) ! obs profile from bufr file real(8), save :: hdr(nhead) ! location and metadata real(8), save :: bout(ndval,max_ch) ! output event integer :: khead, kchan character(len=80) :: hdstr ! header string character(len=10) :: chstr ! channel string data hdstr /'SIID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SAZA FOVN SOZA LSQL SAID HOLS HMSL ' / data chstr /'CHNM TMBR'/ ! order matter here ! In the hearder string, clat and clon should always be the 2nd and 3rd values ! In the hearder string, YEAR etc. should always be values 4-9 ! ! Check if this is a read or write if (lread) then nobs=nobs+1 ! ! read information from the input file needed for interpolation call ufbint(luin,hdr,nhead,1,khead,hdstr) call ufbrep(luin,obsdat,ndval,max_ch,obs_n_channels,chstr) ! ! copy into main program arrays do i=1,obs_n_channels obs_channels(i,1)=obsdat(1,i) obs_values(i,1)=obsdat(2,i) enddo ! obs_n_data=1 obs_info(1:nhead)=hdr(1:nhead) ! ! End of read of observation buffer ! else ! ! This is a call to write observation buffer ! if (obs_n_channels > 0) then ! ! Write out channel observations bout=bmiss do i=1,obs_n_channels bout(2,i)=obs_values(i,1) bout(1,i)=obs_channels(i,1) enddo ! call ufbrep(luout,bout,ndval,obs_n_channels,kchan,chstr) if (obs_n_channels /= kchan) then print *,'incompatible level number returned in write ', & 'of channel data' ierror=ierror+1 endif ! ! Write out the header information hdr(1:nhead)=obs_info(1:nhead) call ufbint(luout,hdr,nhead,1,khead,hdstr) ! ! finished updates; so write output buffer to output file for this report ! call writsb(luout) ! endif ! test on obs_n_channels ! endif ! test on whether to read or write observation buffer ! ! Optionally print some info if testing mode if (ltest) then print *,' ' if (lread) then print *,' Read observation buffer' else print *,' Write observation buffer' endif print ('(a,2i8)'),' OBS number, channels:',nobs,obs_n_channels print ('(a,12f10.3)'),' OBS HDR:',hdr if (obs_n_channels > 0) then do i=1,obs_n_channels print ('(i4,4f15.5)'),i,obs_channels(i,1:2),obs_values(i,1:2) enddo endif endif ! test on ltest ! end subroutine read_write_obs_amsu_hirs ! ! end module m_bufr_rw ! ! ! X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X ! logical function check_type(subset,type) character(len=*) subset character(len=*) type ! Code for making data selection of subtypes given the data type requested: ! subset is the data subset type name read from the BUFR files ! type is the data type specified in the main program argument list ! character(len=8), parameter :: temp_types(7) = (/ 'ADPUPA', & 'RASSDA', 'SFCSHP', 'ADPSFC', 'MSONET', 'AIRCAR', 'AIRCFT' /) ! character(len=8), parameter :: wind_types(13) = (/ 'SYNDAT', & 'ADPUPA','PROFLR','VADWND','AIRCFT','AIRCAR','SATWND', & 'SFCSHP','SPSSMI','ADPSFC','QKSWND','ERS1DA','MSONET' /) ! character(len=8), parameter :: hirs2_types(1) = 'NC021021' character(len=8), parameter :: hirs3_types(1) = 'NC021025' character(len=8), parameter :: airs_types(1) = 'NC021250' character(len=8), parameter :: amsua_types(1) = 'NC021023' character(len=8), parameter :: amsub_types(1) = 'NC021024' character(len=8), parameter :: msu_types(1) = 'NC021022' ! ! note for surface data may also want to include ADPUPA (upperair) ! since that also has sfc obs where CAT=0. The types listed ! below are ones that are -only- surface data character(len=8), parameter :: sfc_types(6) = (/ 'SFCSHP', & 'SPSSMI','ADPSFC','QKSWND','ERS1DA','MSONET' /) character(len=8), parameter :: sst_info(3) = (/ 'ADPUPA', & & 'ADPSFC', 'SFCSHP' /) select case ( type ) case ('MASS_') check_type = any( subset == temp_types ) case ('WIND_') check_type = any( subset == wind_types ) case ('HIRS2') check_type = any( subset == hirs2_types ) case ('HIRS3') check_type = any( subset == hirs3_types ) case ('AIRS_') check_type = any( subset == airs_types ) case ('AMSUA') check_type = any( subset == amsua_types ) case ('AMSUB') check_type = any( subset == amsub_types ) case ('_MSU_') check_type = any( subset == msu_types ) end select return end ! !