source: branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbsel.F90 @ 2893

PROGRAM fbsel
   USE obs_fbm
   USE date_utils
   IMPLICIT NONE
   TYPE(obfbdata) :: fbdatain,fbdataout
   CHARACTER(len=256) :: filenamein,filenameout,filenametmp,cnameout
#ifndef NOIARGCPROTO
   INTEGER,EXTERNAL :: iargc
#endif
   INTEGER,PARAMETER :: maxtyp=1023
   INTEGER,PARAMETER :: maxdates=20
   INTEGER :: nqc,ntyp,ndates,ninidate(maxdates),nenddate(maxdates)
   LOGICAL :: lsplitqc,lsplittyp
   INTEGER :: iqc,ityp,idate
   REAL :: maxlat,minlat,maxlon,minlon
   NAMELIST/namsel/nqc,ntyp,ndates,ninidate,nenddate,lsplitqc,lsplittyp, &
      &            maxlat,minlat,maxlon,minlon

   IF (iargc()/=2) THEN
      WRITE(*,*)'Usage:'
      WRITE(*,*)'fbsel <input filename> <output filename>'
      CALL abort
   ENDIF

   CALL getarg(1,filenamein)
   CALL getarg(2,filenameout)
   
   nqc=-1
   ntyp=-1
   ndates=1
   ninidate=19500101
   nenddate=21000101

   lsplitqc=.FALSE.
   lsplittyp=.FALSE.
   maxlat=1e+38
   minlat=-1e+38
   maxlon=1e+38
   minlon=-1e+38
   OPEN(10,file='namsel.in')
   READ(10,namsel)
   CLOSE(10)
   WRITE(*,namsel)

   CALL init_obfbdata(fbdatain)
   CALL init_obfbdata(fbdataout)

   WRITE(*,*)'Reading file : ',TRIM(filenamein)
   CALL read_obfbdata(TRIM(filenamein),fbdatain)
   WRITE(*,*)'Number of observations before selection = ',fbdatain%nobs
   DO idate=1,ndates
      IF (ndates==1) THEN
         cnameout=filenameout
      ELSE
         WRITE(cnameout,'(I2.2,2A)')idate,'_',TRIM(filenameout)
      ENDIF
      IF (lsplitqc) THEN
         DO iqc=1,3
            CALL fb_sel(fbdatain,fbdataout,iqc,ntyp, &
               &        ninidate(idate),nenddate(idate), &
               &        maxlat,minlat,maxlon,minlon)
            WRITE(filenametmp,'(A,I2.2,A,A)')'qc_',iqc,'_',TRIM(cnameout)
            IF (fbdataout%nobs>0) THEN
               WRITE(*,*)'QC selected = ',iqc
               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
            ENDIF
            CALL dealloc_obfbdata(fbdataout)
         ENDDO
      ELSEIF (lsplittyp) THEN
         DO ityp=0,maxtyp
            CALL fb_sel(fbdatain,fbdataout,nqc,ityp, &
               &        ninidate(idate),nenddate(idate), &
               &        maxlat,minlat,maxlon,minlon)
            WRITE(filenametmp,'(A,I4.4,A,A)')'typ_',ityp,'_',TRIM(cnameout)
            IF (fbdataout%nobs>0) THEN
               WRITE(*,*)'Type = ',ityp
               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
            ENDIF
            CALL dealloc_obfbdata(fbdataout)
         ENDDO
      ELSE
         CALL fb_sel(fbdatain,fbdataout,nqc,ntyp, & 
            &        ninidate(idate),nenddate(idate), &
            &        maxlat,minlat,maxlon,minlon)
         WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
         WRITE(*,*)'Writing file : ',TRIM(cnameout)
         CALL write_obfbdata(TRIM(cnameout),fbdataout)
         CALL dealloc_obfbdata(fbdataout)
      ENDIF
   ENDDO

CONTAINS

   SUBROUTINE fb_sel(fbdatain,fbdataout,nqc,ntyp,ninidate,nenddate,&
      &              maxlat,minlat,maxlon,minlon)
      TYPE(obfbdata) :: fbdatain,fbdataout
      INTEGER :: nqc,ntyp,ninidate,nenddate
      REAL :: maxlat,minlat,maxlon,minlon
      INTEGER :: jobs
      INTEGER :: iqc,ityp
      LOGICAL :: llvalid(fbdatain%nobs)
      INTEGER :: iyea,imon,iday
      REAL(KIND=dp) :: zjini,zjend

      iyea=ninidate/10000
      imon=ninidate/100-iyea*100
      iday=ninidate-iyea*10000-imon*100
      CALL greg2jul(0,0,0,iday,imon,iyea,zjini)
      iyea=nenddate/10000
      imon=nenddate/100-iyea*100
      iday=nenddate-iyea*10000-imon*100
      CALL greg2jul(0,0,0,iday,imon,iyea,zjend)
      DO jobs = 1, fbdatain%nobs
         llvalid(jobs)=.TRUE.
         IF (nqc/=-1) THEN
            CALL check_prof(fbdatain,jobs,iqc)
            llvalid(jobs)=(iqc==nqc).AND.llvalid(jobs)
         ENDIF
         IF (ntyp/=-1) THEN
            READ(fbdatain%cdtyp(jobs),'(I4)')ityp
            llvalid(jobs)=(ityp==ntyp).AND.llvalid(jobs)
         ENDIF
         IF (ninidate/=-1) THEN
            llvalid(jobs)=(fbdatain%ptim(jobs)>zjini).AND.llvalid(jobs)
         ENDIF
         IF (nenddate/=-1) THEN
            llvalid(jobs)=(fbdatain%ptim(jobs)<=zjend).AND.llvalid(jobs)
         ENDIF
         llvalid(jobs)=(fbdatain%pphi(jobs)<=maxlat).AND. &
            &          (fbdatain%pphi(jobs)>=minlat).AND. &
            &          (fbdatain%plam(jobs)<=maxlon).AND. &
            &          (fbdatain%plam(jobs)>=minlon).AND.llvalid(jobs)
         ! Add more checks here...
      ENDDO

      CALL subsamp_obfbdata(fbdatain,fbdataout,llvalid)

   END SUBROUTINE fb_sel

   SUBROUTINE check_prof(fbdata,iobs,iqc)
      
      TYPE(obfbdata) :: fbdata
      INTEGER :: iobs,iqc
      INTEGER :: i,ivar
      
      LOGICAL :: lpart,lfull

      lpart=.false.
      lfull=.true.
      DO ivar=1,fbdata%nvar
         DO i=1,fbdata%nlev
            IF ((fbdata%ivlqc(i,iobs,ivar)>2).AND.&
               &(fbdata%ivlqc(i,iobs,ivar)<9)) lpart = .TRUE.
            IF (fbdata%ivlqc(i,iobs,ivar)<=2) lfull = .FALSE.
         ENDDO
      ENDDO

      IF(lfull) THEN
         iqc=3
      ELSEIF (lpart) then
         iqc=2
      ELSE
         iqc=1
      ENDIF

   END SUBROUTINE check_prof

END PROGRAM fbsel