PROGRAM fbstat
   USE obs_fbm
   USE fbaccdata
   USE coords
   USE omonainfo
   IMPLICIT NONE
   TYPE(obfbdata) :: fbdata
   CHARACTER(len=256) :: filename
   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji
#ifndef NOIARGCPROTO
   INTEGER,EXTERNAL :: iargc
#endif
   REAL,DIMENSION(:),ALLOCATABLE :: zlev
   INTEGER :: nmlev, nfiles
   LOGICAL :: lexists,lomona,ltext
   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdat3d
   REAL, ALLOCATABLE, DIMENSION(:,:) :: zdat2d
   INTEGER,DIMENSION(1) :: itime
   INTEGER :: inidate,icurdate,loopno,ityp,iloopno
   INTEGER :: nvar,nadd
   CHARACTER(len=4) :: expver
   CHARACTER(len=7) :: cltyp
   CHARACTER(len=128) :: cdfmthead,cdfmtbody
   LOGICAL :: lnear,linner,linnerp,linnerini,lpassive,lhistogram,lfound
   INTEGER :: nqc
   CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: cname,caddname
   CHARACTER(len=20) :: carea
   INTEGER :: jstabox, jendbox
   NAMELIST/namfbstat/ltext,lomona,nmlev,inidate,icurdate,loopno,&
      &               expver,lnear,linner,lpassive,lhistogram, &
      &               zhistmax,zhistmin,zhiststep,zcheck,carea

   ltext=.TRUE.
   lomona=.FALSE.
   nmlev=31
   inidate=19010101
   icurdate=19010116
   loopno=-1
   expver='test'
   lnear=.TRUE.
   linner=.FALSE.
   lpassive=.FALSE.
   lhistogram=.FALSE.
   zhistmin(:)=-10.0
   zhistmax(:)=10.0
   zhiststep(:)=0.1
   zcheck(:)=10.0
   nqc=2
   carea='all'
   INQUIRE(file='namfbstat.in',exist=lexists)
   IF (lexists) THEN
      OPEN(10,file='namfbstat.in')
      READ(10,namfbstat)
      CLOSE(10)
      WRITE(*,namfbstat)
   ENDIF
   IF (iargc()==0) THEN
      WRITE(*,*)'Usage:'
      WRITE(*,*)'fbstat [-nmlev] <filenames>'
      CALL abort
   ENDIF
   jfirst=1
   DO ji=1,2
      CALL getarg(jfirst,filename)
      IF (filename=='-42') THEN
         nmlev=42
         jfirst=jfirst+1
      ELSEIF(filename=='-31') THEN
         nmlev=31
         jfirst=jfirst+1
      ELSEIF(filename=='-1') THEN
         nmlev=1
         lnear=.TRUE.
         jfirst=jfirst+1
      ELSEIF(filename=='-q') THEN
         jfirst=jfirst+1
         CALL getarg(jfirst,filename)
         READ(filename,'(I4)')nqc
         IF ((nqc<0).OR.(nqc>4)) THEN
            WRITE(*,*)'Quality control option (-q) should be 1 to 4'
            CALL abort
         ENDIF
         jfirst=jfirst+1
      ENDIF
   END DO
   nfiles=iargc()
   
   IF (carea/='all') THEN
      IF (lomona) THEN
         WRITE(*,*)'For omona files carea has to be all'
         CALL abort
      ENDIF
      lfound=.FALSE.
      DO jbox=1,nbox
         IF (TRIM(carea)==TRIM(cl_boxes(jbox))) THEN
            jstabox=jbox
            jendbox=jbox
            lfound=.TRUE.
         ENDIF
      ENDDO
      IF (.NOT.lfound) THEN
         WRITE(*,*)'Area not found'
         CALL abort
      ENDIF
   ELSE
      jstabox=1
      jendbox=nbox
   ENDIF
   PRINT *,'jstabox,jendbox=',jstabox,jendbox

   IF (.NOT.lnear) nmlev=nmlev-1

   ALLOCATE(&
      & zlev(nmlev) &
      & )
   IF(lnear) THEN
      CALL getlevs(nmlev,zlev)
   ELSE
      CALL getlevsmean(nmlev,zlev)
   ENDIF

   DO jfile=jfirst, nfiles
      CALL getarg(jfile,filename)
      WRITE(*,*)'Handling file : ',TRIM(filename)
      CALL read_obfbdata(TRIM(filename),fbdata)
      IF (jfile==jfirst) THEN
         nvar=fbdata%nvar
         nadd=fbdata%nadd
         IF (lhistogram) THEN
            IF (nvar>maxvars) THEN
               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
               WRITE(*,*)'if you want histograms'
               CALL abort
            ENDIF
            DO jvar = 1, nvar
               hist(jvar)%npoints=(zhistmax(jvar)-zhistmin(jvar))&
                  &              /zhiststep(jvar)+1
               WRITE(*,*)'Number of points in histogram      = ',&
                  &  hist(jvar)%npoints
               WRITE(*,*)'Size of histogram array (elements) = ',&
                  &  hist(jvar)%npoints*nmlev*nbox*nadd
               ALLOCATE(&
                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nbox,nadd) &
                  & )
               hist(jvar)%nhist(:,:,:,:)=0
            ENDDO
         ENDIF
         ALLOCATE(&
            & inum(nmlev,nbox,nadd,nvar),  &
            & zmean(nmlev,nbox,nadd,nvar), &
            & zrms(nmlev,nbox,nadd,nvar),  &
            & zsdev(nmlev,nbox,nadd,nvar), &
            & cname(nvar),                 &
            & caddname(nadd)               &
            & )
         DO jvar = 1, nvar
            cname(jvar) = fbdata%cname(jvar)
         END DO
         DO jadd = 1, nadd
            caddname(jadd) = fbdata%caddname(jadd)
         END DO
         inum(:,:,:,:)=0
         zmean(:,:,:,:)=0.0
         zrms(:,:,:,:)=0.0
         zsdev(:,:,:,:)=0.0
      ELSE
         IF (fbdata%nvar/=nvar) THEN
            WRITE(*,*)'Different number of nvar ',fbdata%nvar,' in ',&
               & TRIM(filename)
            CALL abort
         ENDIF
         IF (fbdata%nadd/=nadd) THEN
            WRITE(*,*)'Different number of nadd ',fbdata%nadd,' in ',&
               & TRIM(filename)
            CALL abort
         ENDIF
      ENDIF
      CALL fb_stat(fbdata,jstabox,jendbox,nmlev,zlev,lnear,nqc,lhistogram)
      CALL dealloc_obfbdata(fbdata)
   ENDDO

   DO jvar=1, nvar
      DO jadd=1, nadd
         DO jbox=jstabox, jendbox
            DO jlev = 1, nmlev
               IF ( inum(jlev,jbox,jadd,jvar) > 0 ) THEN
                  zmean(jlev,jbox,jadd,jvar) = &
                     & zmean(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar)
                  zrms(jlev,jbox,jadd,jvar) = &
                     & SQRT(zrms(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar))
                  zsdev(jlev,jbox,jadd,jvar) = &
                     & SQRT(MAX(zrms(jlev,jbox,jadd,jvar)**2-zmean(jlev,jbox,jadd,jvar)**2,0.0))
               ELSE
                  zmean(jlev,jbox,jadd,jvar) = fbrmdi
                  zrms(jlev,jbox,jadd,jvar) = fbrmdi
                  zsdev(jlev,jbox,jadd,jvar) = fbrmdi
               ENDIF
            ENDDO
         ENDDO
      ENDDO
   ENDDO

   IF (ltext) THEN
      DO jvar=1, nvar
         DO jadd=1, nadd
            DO jbox=jstabox, jendbox
               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
                  &                  TRIM(caddname(jadd)),'_',&
                  &                  TRIM(cl_boxes(jbox)),'.dat'
               OPEN(10,file=TRIM(filename))
               DO jlev = 1, nmlev
                  WRITE(10,'(F16.7,2I12,3F17.10)') zlev(jlev), &
                     & jlev, inum(jlev,jbox,jadd,jvar), &
                     & zmean(jlev,jbox,jadd,jvar), &
                     & zrms(jlev,jbox,jadd,jvar), &
                     & zsdev(jlev,jbox,jadd,jvar)
               ENDDO
               CLOSE(10)
            ENDDO
         ENDDO
      ENDDO
   ENDIF

   IF (lhistogram) THEN
      DO jvar=1, nvar
         DO jadd=1, nadd
            DO jbox=jstabox, jendbox
               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
                  &                  TRIM(caddname(jadd)),'_',&
                  &                  TRIM(cl_boxes(jbox)),'_histogram.dat'
               OPEN(10,file=TRIM(filename))
               WRITE(10,'(A10,100F10.2)')'#',(zlev(jlev),jlev=1,nmlev)
               DO ji=1,hist(jvar)%npoints
                  WRITE(10,'(F10.2,100I10)') &
                     & zhistmin(jvar)+(ji-1)*zhiststep(jvar), &
                     & (hist(jvar)%nhist(ji,jlev,jbox,jadd),jlev=1,nmlev)
               ENDDO
               CLOSE(10)
            ENDDO
         ENDDO
      ENDDO
   ENDIF

   IF (lomona) THEN

      IF (nmlev>1) THEN
         ALLOCATE(zdat3d(nmlev,nbox,1))
      ELSE
         ALLOCATE(zdat2d(nbox,1))
      ENDIF

      cl_expnam=expver
      WRITE(cl_date,'(I8.8)')inidate
      i_dp = nmlev
      itime=icurdate
      linnerp=.TRUE.
      iloopno = loopno
      linnerini = linner
      DO jvar = 1, nvar
         linner = linnerini
         loopno = iloopno
         SELECT CASE (TRIM(cname(jvar)))
         CASE('POTM')
            cl_var = 'votemper'
         CASE('PSAL')
            cl_var='vosaline'
         CASE('SLA')
            cl_var='sossheig'
         CASE('SST')
            cl_var='sosstsst'
         END SELECT
         DO jadd = 1, nadd
            linner = (caddname(jadd)(1:3)=='Hxa').OR.linner
            IF (lpassive) THEN
               ityp=145
            ELSE
               IF (linner) THEN
                  linnerp=.TRUE.
                  ityp=144
                  IF (jadd>1) loopno=loopno+1
               ELSE
                  ityp=142
                  IF (.NOT.linnerp) THEN
                     IF (jadd>1) loopno=loopno+1
                  ENDIF
               ENDIF
            ENDIF
            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
            CALL obs_variable_att(cltyp)
            IF (nmlev>1) THEN
               zdat3d(:,:,1) = zmean(:,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
            ELSE
               zdat2d(:,1) = zmean(1,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
            ENDIF
            IF (lpassive) THEN
               ityp=245
            ELSE
               IF (linner) THEN
                  ityp=244
               ELSE
                  ityp=242
               ENDIF
            ENDIF
            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
            CALL obs_variable_att(cltyp)
            IF (nmlev>1) THEN
               zdat3d(:,:,1) = zrms(:,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
            ELSE
               zdat2d(:,1) = zrms(1,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
            ENDIF
            IF (lpassive) THEN
               ityp=345
            ELSE
               IF (linner) THEN
                  ityp=344
               ELSE
                  ityp=342
               ENDIF
            ENDIF
            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
            CALL obs_variable_att(cltyp)
            IF (nmlev>1) THEN
               zdat3d(:,:,1) = zsdev(:,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
            ELSE
               zdat2d(:,1) = zsdev(1,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
            ENDIF
            IF (lpassive) THEN
               ityp=445
            ELSE
               IF (linner) THEN
                  ityp=444
               ELSE
                  ityp=442
               ENDIF
            ENDIF
            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
            CALL obs_variable_att(cltyp)
            IF (nmlev>1) THEN
               zdat3d(:,:,1) = inum(:,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
            ELSE
               zdat2d(:,1) = inum(1,:,jadd,jvar)
               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
                  &                    cl_boxes,REAL(fbrmdi))
            ENDIF
            linner=.FALSE.
         ENDDO
      ENDDO

      IF (nmlev>1) THEN
         DEALLOCATE(zdat3d)
      ELSE
         DEALLOCATE(zdat2d)
      ENDIF
      
   ENDIF

CONTAINS

   SUBROUTINE fb_stat(fbdata,nstabox, nendbox, nmlev,zlev,lnear,kqc,lhist)
      USE fbaccdata
      USE coords
      TYPE(obfbdata) :: fbdata
      INTEGER :: nstabox,nendbox
      INTEGER :: nmlev
      REAL :: zlev(nmlev)
      LOGICAL :: lnear
      INTEGER :: kqc
      LOGICAL :: lhist
      INTEGER :: jlev, jobs, jvar, klev,jlev2,ih
      REAL :: zarea(4),zlat,zlon,zdiff,zdiff2
      
      DO jbox = nstabox, nendbox
         CALL coord_area(cl_boxes(jbox),zarea)
         DO jobs = 1, fbdata%nobs
            zlat = fbdata%pphi(jobs)
            zlon = fbdata%plam(jobs)
            IF (zlon<0) zlon=zlon+360
            IF (zlon>360) zlon=zlon-360
            IF ( ( zlat .GE. zarea(3) ) .AND. &
               & ( zlat .LE. zarea(4) ) .AND. &
               & ( ( ( zlon .GE. zarea(1) ) .AND. &
               &     ( zlon .LE. zarea(2) ) ) .OR. &
               &   ( ( zarea(2) .LE. zarea(1) ) .AND. &
               &     ( zlon .GE. zarea(1) ) .AND. &
               &     ( zlon .LE. 360        ) ) .OR. &
               &   ( ( zarea(2)   .LE. zarea(1) ) .AND. &
               &     ( zlon .GE. 0          ) .AND. &
               &     ( zlon .LE. zarea(2) ) ) ) ) THEN
            
               DO jlev = 1, fbdata%nlev
                  DO jvar = 1, fbdata%nvar
                     IF (nmlev==1) THEN
                        klev=1
                     ELSE
                        IF (lnear) THEN
                           zdiff=ABS(fbdata%pdep(jlev,jobs)-zlev(1))
                           klev=1
                           DO jlev2=2,nmlev
                              zdiff2=ABS(fbdata%pdep(jlev,jobs)-zlev(jlev2))
                              IF (zdiff2<zdiff) THEN
                                 klev=jlev2
                                 zdiff=zdiff2
                              ENDIF
                           ENDDO
                        ELSE
                           klev = fbdata%iobsk(jlev,jobs,jvar)-1
                        ENDIF
                        IF ( klev > nmlev ) THEN
                           DO jadd = 1, fbdata%nvar
                              IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar))<9000 ) THEN
                                 WRITE(*,*)'Error in fb_stat'
                                 WRITE(*,*)'Increase nmlev to at least ',klev
                                 klev=nmlev
                                 CALL abort
                              ENDIF
                           ENDDO
                        ENDIF
                     ENDIF
                     IF ( fbdata%ivlqc(jlev,jobs,jvar) > kqc ) CYCLE 
                     IF (( klev > 0 ).AND. &
                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
                        DO jadd = 1, fbdata%nadd
                           IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar)) < 9000 ) THEN
                              zdiff = ( fbdata%padd(jlev,jobs,jadd,jvar) - &
                                 &      fbdata%pob(jlev,jobs,jvar) )
                              inum(klev,jbox,jadd,jvar) = inum(klev,jbox,jadd,jvar) + 1
                              zmean(klev,jbox,jadd,jvar) = zmean(klev,jbox,jadd,jvar) + &
                                 & zdiff
                              zrms(klev,jbox,jadd,jvar) = zrms(klev,jbox,jadd,jvar) + &
                                 & zdiff * zdiff
                              IF (ABS(zdiff)>zcheck(jvar)) THEN
                                 WRITE(*,*)'Departure outside check range ',&
                                    & TRIM(fbdata%cname(jvar)),' entry ',&
                                    & fbdata%caddname(jadd)
                                 WRITE(*,*)'Depar = ',zdiff
                                 WRITE(*,*)'Check = ',zcheck(jvar)
                                 WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
                                 WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
                                 WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
                                 WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
                                 WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
                                 WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
                                 WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
                                 WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
                              ENDIF
                              IF (lhist) THEN
                                 ih=NINT((zdiff-zhistmin(jvar))/zhiststep(jvar))+1
                                 IF ((ih>=1).AND.(ih<=hist(jvar)%npoints)) THEN
                                    hist(jvar)%nhist(ih,klev,jbox,jadd) = &
                                       hist(jvar)%nhist(ih,klev,jbox,jadd) +1
                                 ELSE
                                    WRITE(*,*)'Histogram value outside range for ',&
                                       & TRIM(fbdata%cname(jvar)),' entry ',&
                                       & fbdata%caddname(jadd)
                                    WRITE(*,*)'Value = ',zdiff
                                    WRITE(*,*)'Range = ',zhistmin(jvar),zhistmax(jvar)
                                    WRITE(*,*)'Step  = ',zhiststep(jvar)
                                    WRITE(*,*)'Index = ',ih
                                    WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
                                    WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
                                    WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
                                    WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
                                    WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
                                    WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
                                    WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
                                    WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
                                 ENDIF
                              ENDIF
                           ENDIF
                        ENDDO
                     ENDIF
                  ENDDO
               ENDDO
            ENDIF
         ENDDO
      ENDDO
      
   END SUBROUTINE fb_stat

   SUBROUTINE getlevsmean(nmlev,zlev)
      IMPLICIT NONE
      INTEGER :: nmlev
      REAL,DIMENSION(nmlev) :: zlev
      REAL,DIMENSION(nmlev+1) :: ztmp
      INTEGER :: i

      zlev(:)=9999.9
      CALL getlevs(nmlev+1,ztmp)
      DO i=1,nmlev
         zlev(i)=0.5*(ztmp(i)+ztmp(i+1))
      ENDDO

   END SUBROUTINE getlevsmean

   SUBROUTINE getlevs(nmlev,zlev)
      IMPLICIT NONE
      INTEGER :: nmlev
      REAL,DIMENSION(nmlev) :: zlev
      
      zlev(:)=9999.9

      IF (nmlev==42) THEN
         zlev(1)=5.02159
         zlev(2)=15.07854
         zlev(3)=25.16046
         zlev(4)=35.27829
         zlev(5)=45.44776
         zlev(6)=55.69149
         zlev(7)=66.04198
         zlev(8)=76.54591
         zlev(9)=87.27029
         zlev(10)=98.31118
         zlev(11)=109.8062
         zlev(12)=121.9519
         zlev(13)=135.0285
         zlev(14)=149.4337
         zlev(15)=165.7285
         zlev(16)=184.6975
         zlev(17)=207.4254
         zlev(18)=235.3862
         zlev(19)=270.5341
         zlev(20)=315.3741
         zlev(21)=372.9655
         zlev(22)=446.8009
         zlev(23)=540.5022
         zlev(24)=657.3229
         zlev(25)=799.5496
         zlev(26)=967.9958
         zlev(27)=1161.806
         zlev(28)=1378.661
         zlev(29)=1615.291
         zlev(30)=1868.071
         zlev(31)=2133.517
         zlev(32)=2408.583
         zlev(33)=2690.780
         zlev(34)=2978.166
         zlev(35)=3269.278
         zlev(36)=3563.041
         zlev(37)=3858.676
         zlev(38)=4155.628
         zlev(39)=4453.502
         zlev(40)=4752.021
         zlev(41)=5050.990
         zlev(42)=5350.272
      ELSEIF (nmlev==31) THEN
         zlev(1)=4.999938
         zlev(2)=15.00029
         zlev(3)=25.00176
         zlev(4)=35.00541
         zlev(5)=45.01332
         zlev(6)=55.0295
         zlev(7)=65.06181
         zlev(8)=75.12551
         zlev(9)=85.25037
         zlev(10)=95.49429
         zlev(11)=105.9699
         zlev(12)=116.8962
         zlev(13)=128.6979
         zlev(14)=142.1953
         zlev(15)=158.9606
         zlev(16)=181.9628
         zlev(17)=216.6479
         zlev(18)=272.4767
         zlev(19)=364.303
         zlev(20)=511.5348
         zlev(21)=732.2009
         zlev(22)=1033.217
         zlev(23)=1405.698
         zlev(24)=1830.885
         zlev(25)=2289.768
         zlev(26)=2768.242
         zlev(27)=3257.479
         zlev(28)=3752.442
         zlev(29)=4250.401
         zlev(30)=4749.913
         zlev(31)=5250.227
      ELSEIF (nmlev==1) THEN
         zlev(1)=0.0
      ELSE
         WRITE(*,*) 'Unknown number of levels'
         CALL abort
      ENDIF

   END SUBROUTINE getlevs

END PROGRAM fbstat