SUBROUTINE getdata (whme, usme, nofile)

use params
#ifdef ice
use ice1
#endif
use gridpar
use wind
use testo
use stat
use map
implicit none

REAL   , INTENT(OUT):: whme(ngx,ngy), usme(ngx,ngy)
LOGICAL, INTENT(OUT):: nofile

INTEGER,PARAMETER :: ISAT = 33
#ifdef ice
INTEGER :: MASK_ICE(NGX,NGY)
#endif

REAL    :: CDG(NGX,NGY)
CHARACTER (LEN=12) :: BEGINDT, ENDDT, CDATE
LOGICAL :: EOFD, FLAGGED
REAL    :: WHSE(NGX,NGY), USSE(NGX,NGY)
INTEGER :: NUMBWH(NGX,NGY),NUMBUS(NGX,NGY)
INTEGER :: I, J, IG, IJ, IX, IY, NN, NNMIN, ICOUNTD, ICOUNTB, ICOUNTF, ICOUNTT
REAL    :: FRIC_VEL, WIND_SPEED, USSE2, WHSE2, RLAT, RLON, TLAT, TLON, FWM, SWH, &
           TAU, WS, XI

!     OPEN DATA FILE
OPEN (ISAT, FILE = 'SATDATA', STATUS = 'OLD', ERR = 400)

GOTO 500
400 CONTINUE
NOFILE = .TRUE.
RETURN
500 CONTINUE

NOFILE = .FALSE.
!*    1. FIX TIME WINDOW EQUAL TO ONE WIND TIME STEP.
!        --------------------------------------------
BEGINDT = IDTPRO
CALL INCDATE (BEGINDT, - IDELWO / 2)
ENDDT = BEGINDT
CALL INCDATE (ENDDT, IDELWO)
WRITE (IU06, * ) ' GRDATA : SELECTING ALTIMETER MEASUREMENTS '
WRITE (IU06,  * ) '          PERIOD FROM ', BEGINDT, ' TO ', ENDDT
!     FIX THE MINIMUM NUMBER OF OBSERVATIONS PER GRID BOX REQUIRED
!     IN ORDER TO ACCEPT THE DATA AT THE GRIB POINTS
IF (XDELLA.GE.0.5) THEN
   NNMIN = 4
ELSE
   NNMIN = 1
ENDIF
WRITE (IU06,  * ) '   GRIDDED DATA WILL BE KEPT WHEN THE NUMBER OF'
WRITE (IU06,  * ) '   OBSERVATIONS PER GRID BOX IS LARGER THAN ',NNMIN - 1
! --- INITIALIZE ARRAYS
DO J = 1, NGY
   DO I = 1, NGX
      WHME (I, J) = 0.
      USME (I, J) = 0.
      WHSE (I, J) = 0.
      USSE (I, J) = 0.
      NUMBWH (I, J) = 0
      NUMBUS (I, J) = 0
   ENDDO
ENDDO
! --- FIND GRIDDED DRAG COEFFICIENT
DO J = 1, NGY
   DO I = 1, NGX
      CDG (I, J) = - 99.
   ENDDO
ENDDO
DO IG = 1, IGL
   DO IJ = IJS (IG), IJL (IG)
      IX = IXLG (IJ, IG)
      IY = KXLT (IJ, IG)
      FRIC_VEL = USOLD (IJ, IG)
      WIND_SPEED = U10OLD (IJ, IG)
      TAU = FRIC_VEL**2 + 0.0001
      CDG (IX, IY) = TAU / (WIND_SPEED**2 + 0.01)
   ENDDO
ENDDO
#ifdef ice
! --- RECONSTRUCT ICE MASK
DO IG = 1, IGL
   DO IJ = IJS (IG), IJL (IG)
      IX = IXLG (IJ, IG)
      IY = KXLT (IJ, IG)
      MASK_ICE (IX, IY) = ICEMASK (IJ, IG)
      IF (MASK_ICE (IX, JY) .EQ.0) CDG (IX, JY) = - 99.
   ENDDO
ENDDO
#endif
! ---- READ MEASUREMENTS
 3000 CONTINUE
flagged = .TRUE.
CALL READSAT (ISAT, CDATE, RLAT, RLON, SWH, WS, flagged, EOFD)
IF (EOFD) GOTO 4000
ICOUNTD = ICOUNTD+1
!     CHECK DATE.
IF (CDATE (1:12) .LT.BEGINDT) then
   icountb = icountb + 1
   GOTO 3000
ENDIF
IF (CDATE (1:12) .LE.ENDDT) THEN
!     CHECK QUALITY.
   IF (flagged) THEN
      icountf = icountf + 1
      GOTO 3000
   ENDIF
!     DATE IS INSIDE THE TIME WINDOW and record flagged as reliable.
!     THE INDECES ARE COMPUTED
   J = NINT ( (RLAT - AMOSOP) / XDELLA) + 1
   XI = MOD (RLON - AMOWEP + 720., 360.)
   I = NINT (XI / XDELLO) + 1

   IF (IPER.EQ.1.AND.I.EQ.NGX + 1) I = 1
   IF (I.GT.NGX.OR.I.LT.1) GOTO 3000

   IF (J.GT.NGY.OR.J.LT.1) GOTO 3000
#ifdef ice
   IF (MASK_ICE (I, J) .EQ.0) GOTO 3000
#endif

   IF (CDG (I, J) .LT.0.) GOTO 3000

   ICOUNTT = ICOUNTT + 1
   IF (SWH.GE.0.) THEN
      WHME (I, J) = WHME (I, J) + SWH
      WHSE (I, J) = WHSE (I, J) + SWH**2
      NUMBWH (I, J) = NUMBWH (I, J) + 1
      IF (WS.GE.0.) THEN
         USME (I, J) = USME (I, J) + WS
         USSE (I, J) = USSE (I, J) + WS**2
         NUMBUS (I, J) = NUMBUS (I, J) + 1
      ENDIF
   ENDIF
   GOTO 3000
ENDIF

4000 CONTINUE

CLOSE (ISAT)

WRITE (IU06, * ) 'MEASUREMENTS READ:', ICOUNTT
! --- DETERMINE MEAN VALUES AND DEVIATIONS
DO J = 1, NGY
   DO I = 1, NGX
      NN = NUMBWH (I, J)
      IF (NN.GE.NNMIN) THEN
!           MEAN VALUES
         WHME (I, J) = WHME (I, J) / REAL (NN)
!           COMPUTE THE STANDARD DEVIATION.
         IF (NN.GE.2) THEN
            WHSE2 = (WHSE(I, J) - WHME(I, J) **2 * REAL(NN)) / (REAL(NN) - 1.)
            WHSE (I, J) = SQRT (MAX (WHSE2, 0.) )
         ENDIF
      ENDIF
   ENDDO
ENDDO
DO J = 1, NGY
   DO I = 1, NGX
      NN = NUMBUS (I, J)
      IF (NN.GE.NNMIN) THEN
!           MEAN VALUES
         USME (I, J) = USME (I, J) / REAL (NN)
!           COMPUTE THE STANDARD DEVIATION.
         IF (NN.GE.2) THEN
            USSE2 = (USSE(I, J) - USME(I, J) **2 * REAL(NN)) / (REAL(NN) - 1.)
            USSE (I, J) = SQRT(MAX (USSE2, 0.) )
         ENDIF
      ENDIF
   ENDDO
ENDDO
! --- FLAG MEASUREMENTS
DO J = 1, NGY
   DO I = 1, NGX
! --- FLAG MEASURED WAVE HEIGHT
      NN = NUMBWH (I, J)
      IF (NN.LT.NNMIN) THEN
!           IF THERE ARE NO MEASUREMENTS
         WHME (I, J) = - 3.
      ENDIF
!           IF SCATTER IS TOO LARGE
      IF (NN.GE.2) THEN
         FWM = 0.5 * WHME (I, J)
         IF (WHSE (I, J) .GT.FWM) WHME (I, J) = - 3.
      ENDIF
! --- FLAG MEASURED FRICTION VELOCITY
      NN = NUMBUS (I, J)
      IF (NN.LT.NNMIN) THEN
!           IF THERE ARE NO MEASUREMENTS
         USME (I, J) = - 3.
      ENDIF
!           IF SCATTER IS TOO LARGE
      IF (NN.GE.2) THEN
         FWM = 0.5 * USME (I, J)
         IF (USSE (I, J) .GT.FWM) USME (I, J) = - 3.
      ENDIF
!         IF THERE IS NO MEASURED WAVE HEIGHT
      IF (WHME (I, J) .LT.0.) USME (I, J) = - 3.
   ENDDO
ENDDO
!*     OBTAIN FRICTION VELOCITIES FROM OBSERVED WINDS AND THE DRAG.
!      ------------------------------------------------------------
DO J = 1, NGY
   DO I = 1, NGX
      IF (USME(I,J).GT.0.) USME(I, J) = SQRT(ABS(CDG(I,J))) * USME (I, J)
      IF (CDG (I, J) .LT.0.0) THEN
         USME (I, J) = - 3.
         WHME (I, J) = - 3.
      ENDIF
   ENDDO
ENDDO
DO J = 1, NGY
   DO I = 1, NGX
      IF (USME (I, J) .GT.0.) THEN
         TLON = (I - 1) * XDELLO + AMOWEP
         TLAT = (J - 1) * XDELLA + AMOSOP
         WRITE(IU06,*) 'GRIDDED VALUE FOR POINT:',TLON,TLAT,WHME(I,J),USME(I, J)
      ENDIF
   ENDDO
ENDDO

RETURN
END SUBROUTINE getdata
