C This file contains Fortran source code for a sample main program and two subroutines: C Subroutine DECODE_SCAT_ORBIT_V04 which reads in L2B data files C and places the data into common /SCAT_L2B/ C Subroutine QC_CHECK which is our suggested algorithm for C discarding observations contaminated by rain or having other problems. C !!! This program assumes you are maintaining the same directory structure !!! C !!! You still need to change the path to the downloaded files !!! C !!! see format statement 9002 below !!! c adapted to read both seawinds and quikscat data July 23, 2004. D.Smith c added comments Jan 1, 2007 D.Smith c changed to v04 May 2, 2011 D.Smith c fixed small error in call to subroutine Sep 23, 2011 D.Smith PROGRAM GET_SCAT_ORBIT_V04 PARAMETER(NCEL=76, NSCAN=1624) LOGICAL*4 LEXIST CHARACTER*100 FILENAME CHARACTER*21 ATIME INTEGER*1 NUDGE_NCEP,NUDGE_ECMWF COMMON /SCAT_L2B/ ATIME(NSCAN),PHI_TRACK(NSCAN),XLAT(NCEL,NSCAN), & XLON(NCEL,NSCAN),ICLASS(NCEL,NSCAN),NUM_AMBIG(NCEL,NSCAN), & ISELECT(NCEL,NSCAN),IRAIN_SCAT(NCEL,NSCAN),WINAL(4,NCEL,NSCAN), & PHIAL(4,NCEL,NSCAN),SOSAL(4,NCEL,NSCAN),WINDS(NCEL,NSCAN), & PHIWS(NCEL,NSCAN),WINGCM(NCEL,NSCAN),DIRGCM(NCEL,NSCAN), & RAD_RAIN(NCEL,NSCAN),MIN_DIFF(NCEL,NSCAN),NUDGE_NCEP,NUDGE_ECMWF PARAMETER(ILU=5) DO 100 IREV=4741,4741 !seawinds use 1867,1867 - for quikscat use 4741,4741 IDIR1=INT(IREV/1000)*1000 IDIR2=IDIR1+999 WRITE(FILENAME,9002) IDIR1,IDIR2,IREV 9002 FORMAT('MYDRIVE:\MYDIRECTORY\',i5.5,'to',i5.5,'\winvec_',I5.5,'_v04') !change this path as needed INQUIRE(FILE=FILENAME,EXIST=LEXIST) IF(.NOT.LEXIST) GOTO 100 WRITE(*,*) FILENAME, IREV OPEN(ILU,FILE=FILENAME,ACCESS="DIRECT",RECL=3738450) CALL DECODE_SCAT_ORBIT_V04(ILU) CLOSE(ILU) DO 60 ISCAN=1172,1176 !for seawinds verification use 1378,1382; for quikscat use 1172,1176 DO 50 ICEL=18,23 !for seawinds verification use 55,60; for quikscat use 18,23 C CALL QC_CHECK(ICEL,ISCAN, IBAD) C IF(IBAD.EQ.1) CYCLE if(iclass(icel,iscan).le.1) cycle write(*,'(2i5,7f8.2,i4,f8.2)') iscan,icel,xlat(icel,iscan), . xlon(icel,iscan),winal(iselect(icel,iscan),icel,iscan), . phial(iselect(icel,iscan),icel,iscan), . sosal(iselect(icel,iscan),icel,iscan), . wingcm(icel,iscan),dirgcm(icel,iscan), . irain_scat(icel,iscan),rad_rain(icel,iscan) 50 CONTINUE 60 CONTINUE 100 CONTINUE END C =========================================================================================== C =========================================================================================== C =========================================================================================== SUBROUTINE DECODE_SCAT_ORBIT_V04(ILU) C ILU is the unit number for the L2B data file C Data are read from ILU and put into common /SCAT_L2B/, which contains the following: C ATIME is the 21 character time string, format= YYYY-DDDTHH:MM:SS.sss C PHI_TRACK is the direction of the subtrack velocity vector relative to north C XLAT is the geodetic latitude C XLON is the east longitude (0-360) C ICLASS indicates the expected quality of the vector retrieval C ICLASS=0 denotes no retrieval was done (either no observations or only one flavor of observation) C ICLASS=1 denotes 2 flavors of observations in wind vector cell C ICLASS=2 denotes 3 flavors of observations in wind vector cell C ICLASS=3 denotes 4 flavors of observations in wind vector cell C We suggest using just cases for which ICLASS.GE.2 However, this will exclude all outer swath data C NUM_AMBIG is the number of ambiguites (0 to 4) C ISELECT is the selected ambiguity (0 to 4) C IRAIN_SCAT is the rain flag derived from the scatterometer measurements C IRAIN_SCAT=1 indicates rain C WINAL is the wind speed for the various ambiguities (10 meter surface winds) C PHIAL is the wind direction for the various ambiguties (oceanographic convention) C SOSAL is the normalized rms after-the-fit residual of the observation minus model sigma-nought C Large SOSAL values indicate the observations did not fit the geophysical model function. C We suggest discarding observations for which SOSAL.GT.1.9. C WINDS is the smoothed version of the selected wind, WINAL(ISELECT) C PHIWS is the smoothed version of the selected wind, PHIAL(ISELECT) C WINGCM is the general circulation model wind speed used for nudging (either NCEP or ECMWF) C DIRGCM is the general circulation model wind direction used for nudging (either NCEP or ECMWF) c c Two flags exist at the end of the swath data file. Each is INTEGER(1) and tell which GCM was used c for nudging (KGCM1=NCEP [0=used, 1=not used], KGCM2=ECMWF [0=used, 1= not used]) c if no GCM was available for a specific date, the orbit was not processed (this has not happened). c In general, ECMWF data are used for nudging in any reprocessed data and NCEP data are used in c all real-time processing. This is due to the fact that ECMWF data are available with a 3-4 day delay. c since the climate model data are trilinearly interpolated, it is possible that one of each field is used C RAD_RAIN is the TMI or SSMI columar rain rate (rain rate times rain column height, km*mm/hour) C C RAD_RAIN=-999 no TMI or SSMI rain avaliable C RAD_RAIN=0 no rain C RAD_RAIN=0.1 possible rain C RAD_RAIN=0.2 through 25.4 definite rain and the given value is the columnar rain rate C "no rain" means no rain was detected within +/- 50 km and +/- time limit given in MIN_DIFF. C "possible rain" means some rain was detected within +/- 50 km and +/- time limit given in MIN_DIFF. C "definite rain" means rain was detected within +/- 25 km and +/- time limit given in MIN_DIFF. C We suggest discarding observations for which RAD_RAIN.GT.0.15. C C MIN_DIFF is the difference in minutes between the scatterometer and radiometer measurements. A value C of 255 implies that no TMI or SSMI collocation could be found within 254 minutes of the SCAT observation. C C NUDGE_NCEP is a variable that shows whether NCEP data are used(0)/not used(1) in nudging scat data. NCEP data C are used for nudging during near real-time processing. C NUDGE_ECMWF shows whether ECMWF data are used(0)/not used(1) in nudging scat data. This occurs only C during reprocessing as ECWMF data are not timely enough (3-4 day delay) for real time processing. PARAMETER(NCEL=76, NSCAN=1624) CHARACTER(21) ATIME INTEGER(1) NUDGE_NCEP,NUDGE_ECMWF COMMON /SCAT_L2B/ ATIME(NSCAN),PHI_TRACK(NSCAN),XLAT(NCEL,NSCAN), & XLON(NCEL,NSCAN),ICLASS(NCEL,NSCAN),NUM_AMBIG(NCEL,NSCAN), & ISELECT(NCEL,NSCAN),IRAIN_SCAT(NCEL,NSCAN),WINAL(4,NCEL,NSCAN), & PHIAL(4,NCEL,NSCAN),SOSAL(4,NCEL,NSCAN),WINDS(NCEL,NSCAN), & PHIWS(NCEL,NSCAN),WINGCM(NCEL,NSCAN),DIRGCM(NCEL,NSCAN), & RAD_RAIN(NCEL,NSCAN),MIN_DIFF(NCEL,NSCAN),NUDGE_NCEP,NUDGE_ECMWF CHARACTER(1) ABUF0(NSCAN),ABUF1(2,NCEL,NSCAN),ABUF2(2,NCEL,NSCAN) CHARACTER(1) ABUF3(NCEL,NSCAN),ABUF4(NCEL,NSCAN) CHARACTER(1) ABUF5S(2,NCEL,NSCAN),ABUF6S(NCEL,NSCAN) CHARACTER(1) ABUF5(2,4,NCEL,NSCAN),ABUF6(4,NCEL,NSCAN) CHARACTER(1) ABUF7(4,NCEL,NSCAN), ABUF8(2,NCEL,NSCAN),ABUF9(NCEL,NSCAN) CHARACTER(1) ABUF10(NCEL,NSCAN),ABUF11(NCEL,NSCAN) READ(ILU,REC=1) ATIME,ABUF0,ABUF1,ABUF2,ABUF3,ABUF4,ABUF5S,ABUF6S,ABUF5, & ABUF6,ABUF7,ABUF8,ABUF9,ABUF10,ABUF11,NUDGE_NCEP,NUDGE_ECMWF DO ISCAN=1,NSCAN IVAL=ICHAR(ABUF0(ISCAN)) IF(IVAL.EQ.255) THEN PHI_TRACK(ISCAN)=-999 !NO PHI_TRACK AVAILABLE ELSE PHI_TRACK(ISCAN)=180.5 +IVAL ENDIF DO ICEL=1,NCEL ILAT1=ICHAR(ABUF1(1,ICEL,ISCAN)) ILAT2=ICHAR(ABUF1(2,ICEL,ISCAN)) ILON1=ICHAR(ABUF2(1,ICEL,ISCAN)) ILON2=ICHAR(ABUF2(2,ICEL,ISCAN)) XLAT(ICEL,ISCAN)=0.01*(256*ILAT1+ILAT2-9000) XLON(ICEL,ISCAN)=0.01*(256*ILON1+ILON2) IVAL=ICHAR(ABUF3(ICEL,ISCAN)) ICLASS(ICEL,ISCAN)=INT(IVAL/64) IVAL=IVAL-64*ICLASS(ICEL,ISCAN) ISELECT(ICEL,ISCAN)=INT(IVAL/8) NUM_AMBIG(ICEL,ISCAN)=IVAL-8*ISELECT(ICEL,ISCAN) IVAL=ICHAR(ABUF4(ICEL,ISCAN)) IRAIN_SCAT(ICEL,ISCAN)=INT(IVAL/64) !LOWER ORDER BITS FOR INTERNAL USE IF(ICLASS(ICEL,ISCAN).EQ.0) CYCLE IWIN1=ICHAR(ABUF5s(1,ICEL,ISCAN)) IWIN2=ICHAR(ABUF5s(2,ICEL,ISCAN)) WINDS(ICEL,ISCAN)=0.01*(256*IWIN1+IWIN2) IPHI=ICHAR(ABUF6S(ICEL,ISCAN)) PHIWS(ICEL,ISCAN)=1.5*IPHI DO IAL=1,NUM_AMBIG(ICEL,ISCAN) IWIN1=ICHAR(ABUF5(1,IAL,ICEL,ISCAN)) IWIN2=ICHAR(ABUF5(2,IAL,ICEL,ISCAN)) WINAL(IAL,ICEL,ISCAN)=0.01*(256*IWIN1+IWIN2) ! 10 m Surface Wind IPHI=ICHAR(ABUF6(IAL,ICEL,ISCAN)) PHIAL(IAL,ICEL,ISCAN)=1.5*IPHI ! oceanographic convention (direction wind blows to) ISOS=ICHAR(ABUF7(IAL,ICEL,ISCAN)) SOSAL(IAL,ICEL,ISCAN) =0.02*ISOS ENDDO IWIN1=ICHAR(ABUF8(1,ICEL,ISCAN)) IWIN2=ICHAR(ABUF8(2,ICEL,ISCAN)) WINGCM(ICEL,ISCAN)=0.01*(256*IWIN1+IWIN2) IPHI=ICHAR(ABUF9(ICEL,ISCAN)) DIRGCM(ICEL,ISCAN)=1.5*IPHI IRAD_RAIN=ICHAR(ABUF10(ICEL,ISCAN)) IF(IRAD_RAIN.EQ.255) THEN RAD_RAIN(ICEL,ISCAN)=-999. !NO TMI or SSMI AVAILABLE ELSE RAD_RAIN(ICEL,ISCAN)=0.1*IRAD_RAIN !0.1 KM MM IS A SPECIAL CODE FOR POSSIBLE RAIN ENDIF IMIN_DIFF=ICHAR(ABUF11(ICEL,ISCAN)) IF(IMIN_DIFF.EQ.255) THEN MIN_DIFF(ICEL,ISCAN)=999999. !NO TMI or SSMI AVAILABLE ELSE MIN_DIFF(ICEL,ISCAN)=IMIN_DIFF !MINUTES ENDIF ENDDO ENDDO RETURN END C =========================================================================================== C =========================================================================================== C =========================================================================================== SUBROUTINE QC_CHECK(ICEL,ISCAN, IBAD) C If the returned value of IBAD = 1, discard the wind vector cell. PARAMETER(NCEL=76, NSCAN=1624) CHARACTER(21) ATIME INTEGER(1) NUDGE_NCEP,NUDGE_ECMWF COMMON /SCAT_L2B/ ATIME(NSCAN),PHI_TRACK(NSCAN),XLAT(NCEL,NSCAN), & XLON(NCEL,NSCAN),ICLASS(NCEL,NSCAN),NUM_AMBIG(NCEL,NSCAN), & ISELECT(NCEL,NSCAN),IRAIN_SCAT(NCEL,NSCAN),WINAL(4,NCEL,NSCAN), & PHIAL(4,NCEL,NSCAN),SOSAL(4,NCEL,NSCAN),WINDS(NCEL,NSCAN), & PHIWS(NCEL,NSCAN),WINGCM(NCEL,NSCAN),DIRGCM(NCEL,NSCAN), & RAD_RAIN(NCEL,NSCAN),MIN_DIFF(NCEL,NSCAN),NUDGE_NCEP,NUDGE_ECMWF IBAD=0 IF(ICLASS(ICEL,ISCAN).LE.1) IBAD=1 IF(SOSAL(1,ICEL,ISCAN).GT.1.9) IBAD=1 IF(IRAIN_SCAT(ICEL,ISCAN).EQ.1) IBAD=1 IF(MIN_DIFF(ICEL,ISCAN).LE.30.AND.RAD_RAIN(ICEL,ISCAN).GT.0.15) IBAD=1 RETURN END