      character*12 cdate,date_start,date_end
      include 'wam.h'                 
      DIMENSION WHT(NGX,NGY), WDIR(NGX,NGY), WMP(NGX,NGY)
      DIMENSION WINDSP(NGX,NGY), WINDDIR(NGX,NGY), WPP(NGX,NGY)
      DIMENSION USTAR(NGX,NGY), Z0(NGX,NGY), ROUGH(NGX,NGY)
      DIMENSION CDG(NGX,NGY),CD(NGX,NGY),UST(NGX,NGY),VST(NGX,NGY)
      DIMENSION DIR(NGX,NGY),PPER(NGX,NGY),WAVE_STRESS(NGX,NGY)
      DIMENSION PERM(NGX,NGY),BETA(NGX,NGY),BETAG(NGX,NGY)
      DIMENSION UDIR(NGX,NGY),VDIR(NGX,NGY),AUX(NGX,NGY),FSM(NGX,NGY)
      DIMENSION TAUWG(NGX,NGY),TAU_REL(NGX,NGY),SIG_WHT(NGX,NGY)
      DIMENSION TAUWG_X(NGX,NGY),TAUWG_Y(NGX,NGY)
      DIMENSION TAUWGHF_X(NGX,NGY),TAUWGHF_Y(NGX,NGY)
      DIMENSION TAUW_X(NGX,NGY),TAUW_Y(NGX,NGY)
      DIMENSION TAUWHF_X(NGX,NGY),TAUWHF_Y(NGX,NGY)
      DIMENSION USTOKES(NGX,NGY),VSTOKES(NGX,NGY)
      DIMENSION MASK(NGX,NGY),WDISS(NGX,NGY),DISS(NGX,NGY)
      character dayc*2,monc*2,yrc*2,kechar*3

        pi = 3.1415927 

c I (1--->NGX): W(-7)------->E(42)
c J (1--->NGY): N(47.25)------->S(30.25)

c DATE START: 01/01/2000 21:00 UTC
      open(10,file='filedate')
         read(10,'(3(1x,a2))') monc,dayc,yrc
      close(10)
      date_start=yrc//monc//dayc//'000000'
      date_end = date_start
      call incdate(date_start,4*3600)

      call incdate(date_end,7*24*3600)

c  ------- SPECIAL VALUE IS 9999.
      spval = 9999.
c WHT --> WAVE HEIGHT FIELD 
c WDIR --> WAVE DIRECTION FIELD (METEOROLOGICAL CONVENTION)
c WMP --> WAVE MEAN PERIOD FIELD
c WPP --> WAVE PEAK PERIOD FIELD
c WINDSP --> WIND SPEED FIELD
c WINDDIR --> WIND DIRECTION FIELD (METEOROLOGICAL CONVENTION)

      OPEN(140,file='MASKA_MED')
       DO I = 1,NGX
        READ(140,'(1000i1)') (MASK(I,J),J=1,NGY)
       END DO
      CLOSE(140)


        do i = 1,ngx
        do j = 1,ngy
         fsm(i,j)=float(mask(i,j))
        end do
        end do


      cdate = date_start

      do ii = 1,1

      CALL  READ_MAP_FILES(CDATE,WHT,WDIR,WMP,WPP,
     + WINDSP,WINDDIR,CDG,TAUWG,UST,VST,ROUGH,WDISS,
     + TAUWG_X,TAUWG_Y,BETAG)

      end do

50    continue


 900  FORMAT(8e12.5)

       close(65)
       close(60)



      stop
      end
c ---------------------------------------------------------------

      SUBROUTINE INCDATE (DATE, ISHIFT)

C ----------------------------------------------------------------------
C
C**** *INCDATE* - TO UPDATE DATE TIME GROUP
C
C     L. BERTOTTI, P.JANSSEN.
C
C     H. GUNTHER   ECMWF  NOVEMBER 1989    NEGATIVE INCREMENTS.
C
C*    PURPOSE.
C     --------
C
C       UPDATING DATE TIME GROUP.
C
C**   INTERFACE.
C     ----------
C
C       *CALL* *INCDATE (IIDATE,ISHIFT)*
C         *IIDATE*  INTEGER - DATE TIME GROUP (YYMMDDHHMM)
C         *ISHIFT* INTEGER - TIME INCREMENT IN SECONDS, WHERE
C                            ABS (ISHIFT) HAS TO BE LESS THEN 1 YEAR
C
C     METHOD.
C     -------
C
C       NONE.
C
C     EXTERNALS.
C     ----------
C
C       NONE.
C
C     REFERENCES.
C     -----------
C
C       NONE.
C
C ----------------------------------------------------------------------
C
      CHARACTER*12 DATE
      DIMENSION MON(12)
C
      DATA MON /31,28,31,30,31,30,31,31,30,31,30,31/
C
C ----------------------------------------------------------------------
C
C*    1.0 SPLITE DATE TIME GROUP INTO MINUTE, HOUR, DAY, MONTH, YEAR.
C         -----------------------------------------------------------
C
 1000 CONTINUE
      READ(DATE,5) IYEAR,MONTH,IDAY,IHOUR,MINUT,ISEC
      IF (MOD(IYEAR,4).EQ.0) THEN
         MON(2) = 29
      ELSE
         MON(2) = 28
      ENDIF

C
C ----------------------------------------------------------------------
C
C*    2.0 ADD SECONDS AND UPDATE DATE AND TIME.
C         -------------------------------------
      ISEC=ISEC+MOD(ISHIFT,60)
      IF (ISEC.GE.60) THEN
        MINUT = MINUT+ISEC/60
        ISEC = ISEC-(ISEC/60)*60
      ELSE IF (ISEC.LT.0) THEN
        MINUT = MINUT +(ISEC-59)/60
        ISEC = ISEC-((ISEC-59)/60)*60
      END IF

C
 2000 CONTINUE
      MINUT=MINUT+ISHIFT/60
C
C     2.1 POSITIVE SHIFT GREATER THAN 1 MINUTE.
C
      IF (MINUT.GE.60) THEN
         IHOUR = IHOUR + MINUT/60
         MINUT = MINUT - (MINUT/60)*60
         IF (IHOUR.GE.24) THEN
            IDAY = IDAY + IHOUR/24
            IHOUR = IHOUR - (IHOUR/24)*24
            IF (IDAY.GT.MON(MONTH)) THEN
 1300          CONTINUE
               IDAY=IDAY-MON(MONTH)
               MONTH=MONTH+1
               IF(MONTH.EQ.13) THEN
                  MONTH = 1
                  IYEAR=MOD(IYEAR+1,100)
                  IF (MOD(IYEAR,4).EQ.0) THEN
                     MON(2)=29
                  ELSE
                     MON(2)=28
                  ENDIF
               END IF
               IF(IDAY.GT.MON(MONTH)) GO TO 1300
            END IF
         END IF
      ELSE IF (MINUT.LT.0) THEN
C
C     2.2 NEGATIVE SHIFT.
C
         IHOUR = IHOUR + (MINUT-59)/60
         MINUT = MINUT - ((MINUT-59)/60)*60
         IF (IHOUR.LT.0) THEN
            IDAY = IDAY + (IHOUR-23)/24
            IHOUR = IHOUR - ((IHOUR-23)/24)*24
            IF (IDAY.LT.1) THEN
 1400          CONTINUE
               MONTH=MONTH-1
               IF(MONTH.EQ.0) THEN
                  MONTH = 12
                  IYEAR=MOD(IYEAR+99,100)
                  IF (MOD(IYEAR,4).EQ.0) THEN
                     MON(2)=29
                  ELSE
                     MON(2)=28
                  ENDIF
               END IF
               IDAY=IDAY+MON(MONTH)
               IF(IDAY.LT.1) GO TO 1400
            END IF
         END IF
      END IF
C
C ----------------------------------------------------------------------
C
C*    3.0 COMPOSE NEW DATE TIME GROUP.
C         ----------------------------
C
 3000 CONTINUE
      WRITE(DATE,5) IYEAR,MONTH,IDAY,IHOUR,MINUT,ISEC
C
    5 FORMAT(6I2.2)
C
      RETURN
      END

      SUBROUTINE READ_MAP_FILES(IIDATE,WAVEHT,WAVEDIR,WAVEMFR,WAVEPFR,
     +  WIND10,WINDDIR,CDG,TAUWG,UST,VST,ROUGH,WDISS,
     + TAUWG_X,TAUWG_Y,BETAG)
C ----------------------------------------------------------------------
C
C ----------------------------------------------------------------------
C
      include 'wam.h'

C
      CHARACTER*12 IIDATE,JJDATE,DATE_START,DATE_END
      CHARACTER MAPFILE*72,DIRFILE*57
      LOGICAL IEOF, FFLAG(8),FFLAGS(4)
      DIMENSION WAVEHT(NGX,NGY), WAVEDIR(NGX,NGY), WAVEMFR(NGX,NGY)
      DIMENSION WINDSPD(NGX,NGY), WINDDIR(NGX,NGY), WAVEPFR(NGX,NGY)
      DIMENSION CDG(NGX,NGY), TAUWG(NGX,NGY),WIND10(NGX,NGY)
      DIMENSION TAUWG_X(NGX,NGY),TAUWG_Y(NGX,NGY),BETAG(NGX,NGY)
      DIMENSION TAUWGHF_X(NGX,NGY),TAUWGHF_Y(NGX,NGY)
      DIMENSION WINDSEA_DIR(NGX,NGY),UST(NGX,NGY),VST(NGX,NGY)
      DIMENSION MASK(NGX,NGY),ROUGH(NGX,NGY),WDISS(NGX,NGY)

C
C ----------------------------------------------------------------------
C
C*    1. DATA HEADER FROM WAVE MODEL OUTPUT (SUB. OUTGRID)
C        -------------------------------------------------
C

c     WRITE(IIDATE,5) IYR,IMON,IDAY,IHR,0,0
C
    5 FORMAT(6I2.2)

       NX = NGX
       NY = NGY

       IUNIT = 80

C ----------------------------------------------------------------------
C
C*    2. DATA FROM WAVE MODEL OUTPUT (SUB OUTGRID).
C        ------------------------------------------
c---------------------------------------------------
      

      DIRFILE=
     + '/vol1/meteo/gkorres/WAM_CPL_ATMOS/wam_in_out/swamp_nest1/'
      MAPFILE=DIRFILE//'MAP'//IIDATE(1:12)
      OPEN(IUNIT,FILE='MAP120105010000',FORM='UNFORMATTED')

      READ(IUNIT, END=3000) JJDATE, DNGX, DNGY,
     1                      AMOWEP, AMOSOP, AMOEAP, AMONOP


      READ(IUNIT, END=3000) (FFLAG(J),J=1,8)
        READ(IUNIT, END=3000) ((WAVEHT(I,K) ,I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WAVEDIR(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WAVEMFR(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WINDSPD(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WINDDIR(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WAVEPFR(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((CDG(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((TAUWG(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((UST(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((VST(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((ROUGH(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((WDISS(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((TAUWG_X(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((TAUWG_Y(I,K),I=1,NX),K=1,NY)
        READ(IUNIT, END=3000) ((BETAG(I,K),I=1,NX),K=1,NY)
        CLOSE(IUNIT)


        amax = -99
c       pp = 0.0
c       cd_m = 0.0
        do i = 1,nx
        do j = 1,ny
c       if (cdg(i,j).gt.amax) then
        if (betag(i,j).gt.amax) then
c       if (rough(i,j).gt.amax) then
c           amax = cdg(i,j)
            amax = betag(i,j)
c           amax = rough(i,j)
c           ii = i
c           jj = j
        end if
        if (waveht(i,j).gt.0.001) then
        pp = pp + 1.
c       cd_m = cd_m + cdg(i,j)
        cd_m = cd_m + betag(i,j)
c       cd_m = cd_m + rough(i,j)
        end if
        end do
        end do
c       print *,'Roughness max and mean values: ', amax , cd_m/pp
C
C ----------------------------------------------------------------------
C
C*    3. END OF FILE
C        -----------
C

 3000 CONTINUE
      CLOSE(IUNIT)

C
	 DO 3200 K=1,NY
	 DO 3200 I=1,NX
	    IF (WAVEHT(I,K).LT.0.00001) THEN
	       WAVEDIR(I,K) = 9999.
	       WAVEHT(I,K) = 9999.
	       WAVEMFR(I,K) = 9999.
	       WAVEPFR(I,K) = 9999.
	    ELSE
c       IF (WAVEDIR(I,K).LT.0.) WAVEDIR(I,K) = WAVEDIR(I,K)+360.
	       WAVEMFR(I,K) = 1./WAVEMFR(I,K)
	       WAVEPFR(I,K) = 1./WAVEPFR(I,K)
c              WAVEDIR(I,K) = mod((wavedir(i,k)+180.),360.)
	    ENDIF
 3200    CONTINUE

              print *,WAVEDIR
              stop
C
C*    3.3.1 CONVERT FRICTION VELOCITIES TO U10.
C           -----------------------------------
C
	 DO 3300 K=1,NY
	 DO 3300 I=1,NX
	    IF (WAVEHT(I,K).LT.0.00001) THEN
	       WIND10(I,K) = 9999.
	    ELSE
        WIND10(I,K) = WINDSPD(I,K)
	    ENDIF
 3300    CONTINUE
C
	 DO 3500 K=1,NY
	 DO 3500 I=1,NX
	    IF (WAVEHT(I,K).LT.0.00001) THEN
	       WINDDIR(I,K) = 9999.
	    ELSE
	       IF (WINDDIR(I,K).LT.0.) WINDDIR(I,K) = WINDDIR(I,K)+360.
	    ENDIF
 3500    CONTINUE


 7000    CONTINUE

 8000    CONTINUE


         RETURN
         END
C
      SUBROUTINE SMOOTH(A,FSM,IN,JN,NITS,IM,JM,B)
C-------------------------------------------------------------------
C     THIS ROUTINE SMOOTHS DATA WITH A FIVE POINT LAPLACIAN FILTER.
C-------------------------------------------------------------------
      DIMENSION A(IN,JN),B(IN,JN),FSM(IN,JN)
      DO 100 N=1,NITS
      do j=1,jm
      do i=1,im
      B(i,j)=A(i,j)
      end do
      end do
      DO 200 J=2,JM-1
      DO 200 I=2,IM-1
      SMFAC=FSM(I+1,J)+FSM(I,J-1)+FSM(I-1,J)+FSM(I,J+1)+1.E-10
      IF (SMFAC.LT.1.) THEN
      B(I,J) = A(I,J)
      GOTO 200
      END IF
      B(I,J)=A(I,J)+(.5/SMFAC)
     1               *(A(I+1,J)*FSM(I+1,J)+A(I,J-1)*FSM(I,J-1)
     2                +A(I-1,J)*FSM(I-1,J)+A(I,J+1)*FSM(I,J+1)
     3                -SMFAC*A(I,J))
  200 CONTINUE
      do j=1,jm
      do i=1,im
      a(i,j)=b(i,j)*fsm(i,j)
      end do
      end do
  100 CONTINUE
      do i=1,im
      A(i,1)=A(i,2)
      end do
      do j=1,jm
      a(im,j)=a(im-1,j)
      end do
      RETURN
      END

