C$$ MAIN PROGRAM DOCUMENTATION BLOCK C C MAIN PROGRAM: BUFR_TRANWINDSAT C PRGMMR: KEYSER ORG: NP22 DATE: 2014-01-20 C C ABSTRACT: READS IN WINDSAT SCATTEROMETER DATA FROM EDR FILES C GENERATED BY FNMOC AND NESDIS AND REFORMATS THEM INTO BUFR SUCH C THAT THEY CAN THEN BE INGESTED INTO THE NCEP BUFR DATA BASE VIA THE C PROGRAM BUFR_TRANJB. C C PROGRAM HISTORY LOG: C 2006-01-20 S. BENDER - ORIGINAL AUTHOR -- BASED ON PROGRAM WRITTEN C BY T. SMITH TO READ WINDSAT EDRs IN BINARY C FORMAT. C 2006-02-09 S. BENDER - NEW DATA FORMAT FROM NAVY, MODIFIED TO ADD C ALL 4 WINDS, + "SELECTED AMBIGUITY" VALUE C BASED ON WIND VECTOR ERROR. C 2006-07-06 D. KEYSER - ACCOUNT FOR NAVY ERROR IN STORING MODEL WIND C DIRECTION (ADD 360 IF LESS THAN ZERO - IF C NAVY FIXES THIS, CODE WILL STILL WORK C PROPERLY. C 2006-07-21 D. KEYSER - CORRECTED WORD LENGTH OF VARAIBLES PASSED C INTO W3LIB ROUTINE W3FC06 (MUST BE REALS OF C MACHINE WORD LENGTH, HAD BEEN HARDWIRED TO C REAL*4), WAS CAUSING WIND COMPONENTS TO BE C GARBAGE WHEN MISMATCH IN ARGUMENT WORD LENGTH C BETWEEN THIS CODE AND W3FC06 C 2006-08-02 D. KEYSER - REDEFINED RANGE OF VALIDITY FOR PARAMETERS C READ FROM EDR FILE BASED ON GROSS VALUE C CHECKS AND MISSING VALUE OF 255 FOR INTEGER*1 C VALUES (PREVIOUSLY ONLY INVALID/MISSING VALUE C FOR REALS AND INTEGERS OF ALL SIZES IN CHECKS C WAS -9999) C 2007-02-05 D. KEYSER - GENERALIZED TO ALSO HANDLE NESDIS WINDSAT C DATA (IN SAME FORMAT AS NAVY WINDSAT DATA); C WRITES NAVY AND WINDSAT DATA TO MESSAGE TYPE C BASED ON VARIABLE "TANKFILE" IN EXECUTING C SCRIPT, INSTEAD OF HARDWIRED TO NC012138 C (CURRENTLY NAVY GOES TO NCO12138 AND NESDIS C GOES TO NC012139) C 2007-05-25 D. KEYSER - MINOR CHANGES TO PREPARE FOR OPERATIONAL C IMPLEMENTATION C 2014-01-20 D. KEYSER - MINOR CHANGES C C USAGE: C INPUT FILES: C UNIT 05 - STANDARD INPUT. W3TRNARG PARSES ARGUMENTS FROM C - STANDARD INPUT. C UNIT 11 - WINDSAT EDR FORMAT ORBITAL FILE. C UNIT 20 - BUFR TABLE FILE CONTAINING BUFR TABLES A, B, AND C D (FOR UNIT 51). C C OUTPUT FILES: C UNIT 06 - PRINTOUT C UNIT 51 - POINTS TO THE OUTPUT BUFR FILE. TRANJB WILL PLACE C THE BUFR MESSAGES INTO THE PROPER TANKS. C C SUBPROGRAMS CALLED: C SYSTEM: - GETENV C LIBRARY: C W3LIB - W3TAGB W3TRNARG W3TAGE ERREXIT C BUFRLIB - DATELEN OPENBF W3MOVDAT W3FC06 OPENMB UFBINT C - WRITSB CLOSBF C C EXIT STATES: C COND = 0 - SUCCESSFUL RUN C = 1 - UNABLE TO PARSE INPUT ARGUMENTS IN W3TRNARG C = 253 - NO REPORTS WRITTEN OUT C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: IBM SP C C$$$ program BUFR_TRANWINDSAT implicit none c WindSat binary record variables c ------------------------------- real*8 xjd2000_8 ! seconds since 1/1/2000, 12UTC real*4 xlatitude_4 +, xlongitude_4 +, scanangle_4 ! angular scan position of the data +, eia_4 ! earth incidence angle +, caa_4 ! compass azimuth angle +, sst_4 ! sea sfc temperature +, vapor_4 ! retr. integr. columnar water vpr (mm) +, cloud_4 ! retr. integr. columnar cld liquid (mm) +, wspd_4(0:3) ! wind spds (up to 4 ambiguities avail) +, wdir_4(0:3) ! wind dirs (up to 4 ambiguities avail) ! ** NOTE! Raw windsat wdir convention ! is oceanographic; need to ! convert to metr. convention ! (90 deg means blowing FROM ! the east, not TOWARDS) +, chisq_4(0:3) ! chi-squared (metric used to rank ! ambiguities) +, xmod_wspd_4 ! model wind speed (probably at 10 m) +, xmod_wdir_4 ! model wind dir (probably at 10 m) +, rainrate_4 ! retr. rain rate integer*4 iscannum_4 ! num of BAPTA spins since start of ! data file +, iSDRrecnum_4 ! corresp. SDR rec# for a given EDR rec# +, iSDRQCflag_4 ! flag table detailing processing or QC ! errors +, iEDRQCflag1_4 ! flag table used to descr. data quality +, iEDRQCflag2_4 ! retrieved quality flag reserved for ! future use integer*2 idcnum_2 ! downcount number +, isfctype_2 ! WindSat surface type +, numambig_2 ! number of ambiguities (wind vectors ! for your choosing) +, iselambig_2 ! selected ambiguity (wind vector that ! won the prize) integer*1 iphiErr_1(0:3)! est. wind dir retr. error covariance +, issterr_1 ! est. retrieval error covariance for ! 1st rank ambiguity +, iwspderr_1 ! est. retrieval error covariance for ! 1st rank ambiguity +, ivaporerr_1 ! est. retrieval error covariance for ! 1st rank ambiguity +, iclouderr_1 ! est. retrieval error covariance for ! 1st rank ambiguity integer irec ! record counter +, iread ! good record read counter character*8 tlflag ! needed by W3TRNARG character*12 subdir,tankid ! needed by W3TRNARG character*80 appchr ! needed by W3TRNARG character*80 tankfile ! read from parent script, ! indicates BUFR message type integer inlun ! lun connected to input WindSat EDR file parameter(inlun = 11) integer ioutlun ! lun connected to output BUFR file parameter(ioutlun = 51) integer idxlun ! lun connected to external BUFR table parameter(idxlun = 20) character*8 msgtyp ! BUFR msg type (NC012138 for WindSat from ! Navy, NC012139 for WindSat from NESDIS) integer msgdate ! BUFR msg date (YYYYMMDDHH) +, idat(8) ! arrays used to convert time in xjd2000_8 +, jdat(8) ! format to calendar date/time real rinc(5) ! ditto real wdir(0:3) ! wdir_4 converted to machine byte size +, wspd(0:3) ! wspd_4 converted to machine byte size +, xmod_wdir ! xmod_wdir_4 converted to machine byte size +, xmod_wspd ! xmod_wspd_4 converted to machine byte size character*80 nemstg ! string of mnemonics representing ! values to be encoded into BUFR subset real*8 arr_sl_8(16,1) ! array used to interface with BUFRLIB ! routine ufbint for single level data real*8 arr_8(4,4) ! array used to interface with BUFRLIB ! routine ufbint for replicated data real*8 arr_iswv_8(1,2) ! array used to encode ISWV into BUFR ! file integer nlevstd ! number of levels stored by BUFRLIB ! software. should be =1, since ! WindSat data is single level data, ! not profile data. integer nrptswrt ! no. of reports written to BUFR file. real*8 bmiss_8 ! missing value flag for BUFRLIB parameter(bmiss_8 = 10E10) integer iveambig ! index (0:3) of ambiguity selected by ! use of wind vector error calculation real prevvecerr ! "previous" vector error +, currvecerr ! "current" vector error integer i ! counter integer jdate,kdate,LAPCHR,IERR,ltnkid,LSUBDR !used by W3TRNARG integer ierr_rec,nskip_qc,nskip_wind,nskip_mod_wind ! counters for skipped rpts integer npe0(0:3),nwmsg(0:3) real xmwdir ! meteorological wind direction +, xmod_u ! model wind u-component +, xmod_v ! model wind v-component +, u ! wind u-component +, v ! wind v-component integer nsstemsg, nreqvmsg ! counters for reports with missing ! SST, REQV c ------------------------------------------------ c Start program. c -------------- C----------------------------------------------------------------------- CALL W3TAGB('BUFR_TRANWINDSAT',2014,0020,0061,'NP22') CALL GETENV('TANKFILE',TANKFILE) PRINT *, ' ' PRINT *, ' ==> Welcome to BUFR_TRANWINDSAT -- Version 01/20/2014' PRINT *, ' ' IF(TANKFILE.EQ.'b012/xx138') THEN PRINT *,' -- processing WindSat data from Navy (NOGAPS model ', $ 'first guess)' ELSE IF(TANKFILE.EQ.'b012/xx139') THEN PRINT *,' -- processing WindSat data from NESDIS (GFS model ', $ 'first guess)' END IF PRINT *, ' ' CALL W3TRNARG(SUBDIR,LSUBDR,TANKID,LTNKID,APPCHR,LAPCHR,TLFLAG, . JDATE,KDATE,IERR) IF(IERR.NE.0) THEN WRITE(6,'('' UNABLE TO PARSE ARGS TO TRANSLATION ROUTINE - '', 1 '' RETURN CODE = '',I5)') IERR CALL W3TAGE('BUFR_TRANWINDSAT') CALL ERREXIT(IERR) ENDIF C----------------------------------------------------------------------- open(inlun,form='unformatted',recl=136,access='direct') call datelen(10) ! YYYYMMDDHH instead of YYMMDDHH in BUFRLIB call openbf(ioutlun,'OUT',idxlun) irec = 0 iread= 0 nrptswrt=0 nskip_qc=0 nskip_wind=0 nskip_mod_wind=0 nreqvmsg=0 nsstemsg=0 ierr_rec=0 npe0=0 nwmsg=0 IF(TANKFILE(1:1).EQ.'b') then msgtyp = 'NC'//TANKFILE(2:4)//TANKFILE(8:10) ! BUFR msg type ELSE msgtyp = 'NC012138' END IF c Read data out of binary file and write it out to BUFR c ----------------------------------------------------- do irec = irec + 1 read(inlun,rec=irec,err=99) + xjd2000_8 +, xlatitude_4 +, xlongitude_4 +, scanangle_4 +, eia_4 +, caa_4 +, iscannum_4 +, idcnum_2 +, isfctype_2 +, iSDRQCflag_4 +, iSDRrecnum_4 +, issterr_1 +, iwspderr_1 +, ivaporerr_1 +, iclouderr_1 +, sst_4 +, vapor_4 +, cloud_4 +, numambig_2 +, iselambig_2 +, wspd_4 +, wdir_4 +, chisq_4 +, xmod_wspd_4 +, xmod_wdir_4 +, iEDRQCflag1_4 +, iEDRQCflag2_4 +, rainrate_4 +, iphiErr_1 iread = iread + 1 wdir = wdir_4 wspd = wspd_4 xmod_wdir = xmod_wdir_4 xmod_wspd = xmod_wspd_4 c Per Jim J., boot obs where all 4 winds are missing or if the QC c values are missing (iEDRQCflag1_4, iphiErr_1, chisq_4). c --------------------------------------------------------------- if + ((wdir(0).lt.0..and.wspd(0).lt.0.) .and. + (wdir(1).lt.0..and.wspd(1).lt.0.) .and. + (wdir(2).lt.0..and.wspd(2).lt.0.) .and. + (wdir(3).lt.0..and.wspd(3).lt.0.)) then nskip_wind = nskip_wind + 1 cycle ! skip report if all 4 winds are missing endif if + ((iEDRQCflag1_4.lt.0) .and. +((iphiErr_1(0).lt.0.or.iphiErr_1(0).ge.255).and.chisq_4(0).lt.0.) + .and. +((iphiErr_1(1).lt.0.or.iphiErr_1(1).ge.255).and.chisq_4(1).lt.0.) + .and. +((iphiErr_1(2).lt.0.or.iphiErr_1(2).ge.255).and.chisq_4(2).lt.0.) + .and. +((iphiErr_1(3).lt.0.or.iphiErr_1(3).ge.255).and.chisq_4(3).lt.0.)) + then nskip_qc = nskip_qc + 1 cycle ! skip report if all QC indicators are missing endif c Convert xjd2000_8 value (seconds since 1/1/00 12UTC) to calendar date c and time. c ------------------------------------------------------------------ idat(1) = 2000 ! starting year idat(2) = 1 ! starting month idat(3) = 1 ! starting day idat(4) = 0 ! starting time zone idat(5) = 12 ! starting hour idat(6) = 0 ! starting minute idat(7) = 0 ! starting seconds idat(8) = 0 ! starting milliseconds rinc(1) = 0 ! increment day rinc(2) = 0 ! increment hour rinc(3) = 0 ! increment minute rinc(4) = xjd2000_8 ! increment seconds rinc(5) = 0 ! increment milliseconds call w3movdat(rinc,idat,jdat) msgdate = jdat(1)*1000000 + jdat(2)*10000 + + jdat(3)*100 + jdat(5) c Convert wind values from oceanograhpic conventions to meteorological c conventions (affects wdir(0:3), xmod_wdir). c Must also account for possible Navy error in storing xmod_wdir by c adding 360 if xmod_wdir less than ZERO - if Navy fixes this, code C will still work properly c -------------------------------------------------------------------- do i = 0,3 if(wdir(i).ge.0.) then xmwdir = mod((int(wdir(i))+180),360) + + wdir(i) - int(wdir(i)) wdir(i) = xmwdir endif enddo if(xmod_wdir.ge.-360.) then if(xmod_wdir.lt.0.) xmod_wdir = xmod_wdir + 360. if(xmod_wdir.ge.0..and.xmod_wdir.le.360.) then xmwdir = mod((int(xmod_wdir)+180),360) + + xmod_wdir - int(xmod_wdir) xmod_wdir = xmwdir else xmod_wdir = -9999. endif endif C Keyser: added skip if model wind is missing C ------------------------------------------- if(xmod_wdir.lt.0..and.xmod_wspd.lt.0.) then nskip_mod_wind = nskip_mod_wind + 1 cycle ! skip report if model wind is missing endif c Determine best ambiguity based on wind vector error. c ---------------------------------------------------- call w3fc06(xmod_wdir,xmod_wspd,xmod_u,xmod_v) ! get model u/v ! components prevvecerr = 99999. ! initial value iveambig = 7 ! initially, set to missing. do i=0,3 call w3fc06(wdir(i),wspd(i),u,v) currvecerr = ((u-xmod_u)**2 + (v-xmod_v)**2)**0.5 if(currvecerr.lt.prevvecerr) then cc print *, currvecerr, prevvecerr, iveambig, i, cc + u, xmod_u, v, xmod_v iveambig = i prevvecerr = currvecerr endif enddo ! i=0,3 c Write report to standard out. c ----------------------------- cc print 888, msgdate, jdat(6), jdat(7), isfctype_2, issterr_1, cc + iwspderr_1, ivaporerr_1, iclouderr_1, sst_4, cc + vapor_4, cloud_4, xmod_wdir, wdir(0:3), xmod_wspd, cc + wspd(0:3), chisq_4(0:3), iphiErr_1(0:3) 888 format(i10,i2.2,i2.2,1x,'SFC TYPE: ', i2,1x, + 'SSTER: ', i3,1x,'WSPDER: ',i3,1x, + 'VPRER: ', i3,1x,'CLDER: ',i3,1x, + 'SST: ',f8.2,1x, 'VPR: ', f8.2,1x, + 'CLD: ',f8.2,1x, 'MOD WDIR: ',f8.2,1x, + 'WDIRS: ',4(f8.2,1x),'MOD WSPD: ',f8.2,1x, + 'WSPDS: ',4(f8.2,1x),'CHSQS: ',4(f8.2,1x), + 'PHERRS: ', 4(i3,1x)) c Open BUFR message. c ------------------ call openmb(ioutlun,msgtyp,msgdate) c Store values in BUFRLIB memory in anticipation of writing BUFR subset. c ---------------------------------------------------------------------- arr_sl_8 = bmiss_8 arr_sl_8(1,1) = jdat(1) arr_sl_8(2,1) = jdat(2) arr_sl_8(3,1) = jdat(3) arr_sl_8(4,1) = jdat(5) arr_sl_8(5,1) = jdat(6) arr_sl_8(6,1) = jdat(7) arr_sl_8(7,1) = xlatitude_4 arr_sl_8(8,1) = xlongitude_4 nemstg = 'YEAR MNTH DAYS HOUR MINU SECO CLAT CLON' call ufbint(ioutlun,arr_sl_8,8,1,nlevstd,nemstg) arr_sl_8 = bmiss_8 arr_sl_8(1,1) = 283. ! SAID = 283 (proposed by J. Ator) if(isfctype_2.ge.0 .and. isfctype_2.le.7) then arr_sl_8(2,1) = isfctype_2 endif if(issterr_1.ge.0 .and. issterr_1.lt.255) then arr_sl_8(3,1) = float(issterr_1)*0.05 ! convert to K per ! windsat doc else nsstemsg = nsstemsg + 1 endif cc ! count those which are = -58 cc 2/24 if(issterr_1.eq.-58) then cc nsstemsg = nsstemsg + 1 cc endif if(iwspderr_1.ge.0 .and. iwspderr_1.lt.255) then arr_sl_8(4,1) = float(iwspderr_1)*0.05 ! convert to m/s per ! windsat doc endif if(ivaporerr_1.ge.0 .and. ivaporerr_1.lt.255) then arr_sl_8(5,1) = float(ivaporerr_1)*0.05 ! convert to mm per ! windsat doc endif if(iclouderr_1.ge.0 .and. iclouderr_1.lt.255) then arr_sl_8(6,1) = float(iclouderr_1)*0.002 ! convert to mm per ! windsat doc endif if(sst_4.gt.200.) then arr_sl_8(7,1) = sst_4 endif if(vapor_4.ge.0.) then arr_sl_8(8,1) = vapor_4 endif if(cloud_4.ge.0.) then arr_sl_8(9,1) = cloud_4 endif if(rainrate_4.ge.0.) then arr_sl_8(10,1) = rainrate_4/3600. ! Navy units = mm/hr, ! REQV = mm/sec else nreqvmsg = nreqvmsg + 1 endif if(xmod_wdir.ge.0.) then arr_sl_8(11,1) = xmod_wdir endif if(xmod_wspd.ge.0.) then arr_sl_8(12,1) = xmod_wspd endif if(iEDRQCflag1_4.ge.0) then arr_sl_8(13,1) = iEDRQCflag1_4 endif nemstg = 'SAID WSST SSTE SPDE VPRE CLDE SST1 '// + 'MRWVC MRLWC REQV MWD10 MWS10 WSEQC1' call ufbint(ioutlun,arr_sl_8,13,1,nlevstd,nemstg) c Add the "selected ambiguity" values, first the originaly that came c with the dataset and the second, the one selected by lowest wind c vector error. c ------------------------------------------------------------------ arr_iswv_8 = bmiss_8 if(iselambig_2.ge.0 .and. iselambig_2.le.3) then arr_iswv_8(1,1) = iselambig_2 endif arr_iswv_8(1,2) = iveambig call ufbint(ioutlun,arr_iswv_8,1,2,nlevstd,'ISWV') c Add 4 sets of winds. c -------------------- arr_8 = bmiss_8 do i = 0,3 if(wdir(i).ge.0.) then arr_8(1,i+1) = wdir(i) endif if(wspd(i).ge.0.) then arr_8(2,i+1) = wspd(i) endif if(chisq_4(i).ge.0.) then arr_8(3,i+1) = chisq_4(i) endif if(iphiErr_1(i).ge.0 .and. iphiErr_1(i).lt.255) then arr_8(4,i+1) = float(iphiErr_1(i))*0.2 ! windsat doc ! says X by ! 0.2 to convert ! to degrees if(iphiErr_1(i).eq.0) then npe0(i) = npe0(i) + 1 endif endif enddo c Write values to memory. nemstg = 'WD10 WS10 CHSQ PHER' call ufbint(ioutlun,arr_8,4,4,nlevstd,nemstg) c Write BUFR subset to output file. c --------------------------------- cc do i = 0,3 cc if(wdir(i).lt.0..and.wspd(i).lt.0.) then cc nwmsg(i) = nwmsg(i) + 1 cc endif cc enddo call writsb(ioutlun) nrptswrt = nrptswrt + 1 cycle 99 continue ! ERR in read comes here print *, "#### Error reading record ",irec ierr_rec = ierr_rec + 1 enddo 20 continue ! ALL DONE (input E-O-F comes here) print * print *,'*** PROCESSING ENDED NORMALLY ***' print * print *, 'TOTAL REPORTS READ FROM EDR FILE: ', iread print *, 'TOTAL NUMBER OF REPORTS WRITTEN TO BUFR FILE: ', + nrptswrt print * print *, 'TOTAL NUMBER OF REPORTS SKIPPED: ', + ierr_rec+nskip_wind+nskip_qc+nskip_mod_wind print *,'... number of reports that could not be read: ', + ierr_rec print *,'... number of reports with all four winds missing: ', + nskip_wind print *,'... number of reports with model wind missing: ', + nskip_mod_wind print *,'... number of reports with all QC indicators missing: ', + nskip_qc print * do i=0,3 print *, 'num rpts w/ iphiErr_1(',i,')=0:',npe0(i) enddo cc do i = 0,3 cc print *, 'num rpts w/ wind(',i,')=missing:',nwmsg(i) cc enddo print * print *, 'Number of reports with missing rain rate: ', nreqvmsg print * print *, 'Number of reports with missing SST: ', nsstemsg print * print *,'*** PROCESSING ENDED NORMALLY ***' print * c Close input file. c ----------------- close(inlun) c Close BUFR file. c ---------------- call closbf(ioutlun) IF(nrptswrt.EQ.0) THEN WRITE(6,2003) 2003 FORMAT(' NO REPORTS PROCESSED -- DISABLING ALL SUBSEQUENT ', 1 'PROCESSING.') CALL W3TAGE('BUFR_TRANWINDSAT') CALL ERREXIT(253) ENDIF CALL W3TAGE('BUFR_TRANWINDSAT') stop 901 format('time: ',i2.2,':',i2.2,':',i2.2,1x, + 'lat: ',f7.2,1x,'lon: ',f7.2,1x, + 'sfc type: ',a12,1x,'nambig: ',i8) 902 format('vapor ',f7.2,1x,'cloud ',f7.2,1x, + 'nambig ',i8,1x,'sel ',i8,1x, + 'wspd ',f6.2,1x,'wdir ',f6.2,1x, + 'qsc ',f6.3,1x, + 'mwsp ',f6.2,1x,'mwdir ',f6.2,1x, + 'EDRQCflag: ',i10) end