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

PROGRAM fbthin
   USE obs_fbm
   IMPLICIT NONE
   !
   ! Command line arguments for output file and input file
   !
#ifndef NOIARGTHINROTO
   INTEGER,EXTERNAL :: iargc
#endif
   INTEGER :: nargs
   CHARACTER(len=256) :: cdoutfile
   CHARACTER(len=256) :: cdinfile
   CHARACTER(len=256) :: cdtmp
   INTEGER :: nout,ninn,nadd,next,i,j,k
   LOGICAL :: lgrid
   !
   ! Feedback data
   !
   TYPE(obfbdata) :: fbdatain
   !
   ! Get number of command line arguments
   !
   nargs=IARGC()
   IF ((nargs /= 2)) THEN
      WRITE(*,'(A)')'Usage:'
      WRITE(*,'(A)')'fbthin inputfile outputfile'
      CALL abort()
   ENDIF
   CALL getarg(1,cdinfile)
   CALL getarg(2,cdoutfile)
   !
   ! Initialize feedback data
   !
   CALL init_obfbdata( fbdatain )
   !
   ! Read the file
   !
   CALL read_obfbdata( TRIM(cdinfile), fbdatain )
   !
   ! Do the thining
   !
   CALL fb_thin( fbdatain )
   !
   ! Write the file
   !
   CALL write_obfbdata( TRIM(cdoutfile), fbdatain )

CONTAINS
   
   SUBROUTINE fb_thin( fbdata )
      !
      ! Observation thinning stub
      !
      IMPLICIT NONE
      TYPE(obfbdata) :: fbdata
      REAL :: dlat = 1.0
      REAL :: dlon = 1.0
      LOGICAL :: lred = .TRUE.
      INTEGER, ALLOCATABLE :: thinobs(:,:), thinqc(:,:)
      REAL, ALLOCATABLE :: thindist(:,:), thingridlat(:), thingridlon(:,:)
      INTEGER :: nlon,nlat
      INTEGER :: ii,ij,iv,ilon,ilat,iqc,iobs,irej
      REAL :: zlon,zdlon,zrpi,zdist
      LOGICAL :: lexists
      NAMELIST/namthin/dlat,dlon,lred
      
      INQUIRE(file='namthin.in',exist=lexists)
      IF (lexists) THEN
         OPEN(10,file='namthin.in')
         READ(10,namthin)
         CLOSE(10)
      ENDIF
      WRITE(*,namthin)

      nlon=360/dlon
      nlat=180/dlat+1
      ALLOCATE(thinobs(nlon,nlat))
      ALLOCATE(thinqc(nlon,nlat))
      ALLOCATE(thindist(nlon,nlat))
      ALLOCATE(thingridlat(nlat))
      ALLOCATE(thingridlon(nlon,nlat))
      thinobs(:,:)=0
      thinqc(:,:)=0
      thindist(:,:)=0
      zrpi=ATAN(1.0)/45.
      DO ij=1,nlat
         thingridlat(ij)=(ij-1)*dlat-90.0
         IF (lred) THEN
            zdlon=MIN(dlon/MAX(COS(ABS(zrpi*thingridlat(ij))),0.0001),360.0)
         ELSE
            zdlon=dlon
         ENDIF
         DO ii=1,nlon
            thingridlon(ii,ij)=(ii-1)*zdlon
         ENDDO
      ENDDO
      irej=0
      iobs=0

      DO ii=1,fbdata%nobs

         ! Skip data with missing lon and lat.
         IF (fbdata%plam(ii)==fbrmdi) CYCLE
         IF (fbdata%pphi(ii)==fbrmdi) CYCLE

         ! Count number of observations with qcflag<=2.
         iqc=0
         DO iv=1,fbdata%nvar
            DO ij=1,fbdata%nlev
               IF(fbdata%ivlqc(ij,ii,iv)<=2) iqc=iqc+1
            ENDDO
         ENDDO

         ! Find ilat,ilon positions
         zlon=fbdata%plam(ii)
         IF (zlon>=360.0) zlon=zlon-360
         IF (zlon<0.0) zlon=zlon+360
         ! If reduced grid change zdlon accordingly
         IF (lred) THEN
            zdlon=MIN(dlon/MAX(COS(ABS(zrpi*fbdata%pphi(ii))),0.0001),360.0)
         ELSE
            zdlon=dlon
         ENDIF
         ilon=INT(zlon/zdlon)+1
         ilat=INT((fbdata%pphi(ii)+90.0)/dlat)+1
         iobs=iobs+1
         ! Check which observation to pick based on number of valid obs.
         IF (thinobs(ilon,ilat)>0) THEN
            irej=irej+1
            IF (iqc>=thinqc(ilon,ilat)) THEN
               zdist=distance(zrpi,zlon,fbdata%pphi(ii),&
                  & thingridlon(ilon,ilat),thingridlat(ilat))
               IF(iqc==thinqc(ilon,ilat)) THEN
                  IF (zdist<thindist(ilon,ilat)) THEN
                     fbdata%ivlqc(:,thinobs(ilon,ilat),:)=4
                     thinobs(ilon,ilat)=ii
                     thinqc(ilon,ilat)=iqc
                     thindist(ilon,ilat)=zdist
                  ELSE
                     fbdata%ivlqc(:,ii,:)=4
                  ENDIF
               ELSE
                  fbdata%ivlqc(:,thinobs(ilon,ilat),:)=4
                  thinobs(ilon,ilat)=ii
                  thinqc(ilon,ilat)=iqc
                  thindist(ilon,ilat)=zdist
               ENDIF
            ELSE
               fbdata%ivlqc(:,ii,:)=4
            ENDIF
         ELSE
            thinobs(ilon,ilat)=ii
            thinqc(ilon,ilat)=iqc
            thindist(ilon,ilat)=distance(zrpi,zlon,fbdata%pphi(ii),&
               & thingridlon(ilon,ilat),thingridlat(ilat))
         ENDIF
      ENDDO
      WRITE(*,*)'Observations considered = ',iobs
      WRITE(*,*)'Observations rejected   = ',irej

      DEALLOCATE(thinobs)
      DEALLOCATE(thinqc)
      DEALLOCATE(thindist)
      DEALLOCATE(thingridlat)
      DEALLOCATE(thingridlon)
      
   END SUBROUTINE fb_thin

   REAL FUNCTION distance( prpi, plon, plat, gridlon, gridlat )
      REAL :: prpi,plon,plat,gridlat,gridlon
      REAL :: zplat,zplon,zgridlat,zgridlon
      REAL :: za1,za2,zb1,zb2,zc1,zc2,zcos1,zcos2

      zplon=plon
      zgridlon=gridlon
      IF (zplon<-180) zplon=zplon+360.0
      IF (zplon>=180) zplon=zplon-360.0
      IF (zgridlon<-180) zgridlon=zgridlon+360.0
      IF (zgridlon>=180) zgridlon=zgridlon-360.0
      zplon=zplon*prpi
      zgridlon=zgridlon*prpi
      zplat=plat*prpi
      zgridlat=gridlat*prpi
      zcos1=COS(zplat)
      zcos2=COS(zgridlat)
      za1=SIN(zplat)
      za2=SIN(zgridlat)
      zb1=zcos1*COS(zplon)
      zb2=zcos2*COS(zgridlon)
      zc1=zcos1*SIN(zplon)
      zc2=zcos2*SIN(zgridlon)

      distance=ASIN(SQRT(1.0-(za1*za2+zb1*zb2+zc1*zc2)**2))
      
    END FUNCTION distance

END PROGRAM fbthin