c --- Univ. of Athens - Ocean Group - G.Korres 07/03/2001

      INCLUDE 'wind.h'
c -------------------------------------------------------------
      PARAMETER (NX=NXP , NY = NYP)
      PARAMETER (NFIELDS = 2, NRAD = 2)
c ---- frequency of meteo updates (in hours)
      PARAMETER (MET_FREQ = 1 )

      character cdate*11, cvar*5 , dump*80
      character mdate*16, mvar*5 
      dimension mask_wam(nx,ny),mask_meteo(nxp,nyp)
      dimension aux(nx,ny),store_1(nx,ny,nfields),
     + store_2(nx,ny,nrad),rain_prev(nx,ny)
      dimension amask(nx , ny),uc(nx,ny),vc(nx,ny)
      integer hours_to_integrate,starting_hour,write_hours
      character*10 starting_date
      logical first
      real*8 time,time0

c --- FIELDS FILE NAMES
      character indf*3,fname1*7,chara*1,directory_fields*80
      character fname*80 , filename*20
c --- RADIATION FILE NAMES
      character inds*3,directory_radiation*80
      character sfile*12
      DATA FNAME1/'TIELDS.'/
      data first/.true./


      NAMELIST/METEO/HOURS_TO_INTEGRATE,STARTING_DATE,
     1 STARTING_HOUR, DIRECTORY_FIELDS


       print *,nx,ny
c ---- read in namelist
      READ(5,METEO)
      WRITE(*,METEO)

c --- define here the location of meteo mask
          open(150,file='ATMOS_MASK20')
           do j = ny,1,-1
            read(150,'(1001i4)') (mask_meteo(i,j),i=1,nx)
           end do
          close(150)
           do i = 1,nx
            do j = 1,ny
             amask(i,j) = float(mask_meteo(i,j))
            end do
           end do
           close(150)


      ldf = 0
      do ii = 1,80
      chara = directory_fields(ii:ii)
      if (chara.eq.'') goto 10
      ldf = ldf + 1
      end do
10    continue





      CALL TRANSLATE(STARTING_DATE,IYR_N,IMO_N,IDAY_N)
      ihour_n = starting_hour
      call findtime0(time0,year,ihour_n,iday_n,imo_n,iyr_n)

      iyr_n = iyr_n-2000
      time = time0
      print *,'IYR_N:',iyr_n



c -----------------------------------------------------------------
c construct inp.prm file here as if everything is to fail!!!
c -----------------------------------------------------------------
        MODEL_RUN = 0  ! no run


c ---- do some checks on starting_hour and hours_to_integrate

c -----------------------------------------------------------------

        last_f = hours_to_integrate/met_freq + 1



        iyear = iyr_n
        imon = imo_n
        iday = iday_n
        ihour = ihour_n 

c ---- do the fields files (u_10,v_10)
        index_f = 0
        index_f = starting_hour - 18      
        index_f = starting_hour 
        do 500 nsteps = 1, last_f

      print *,'*******************************'
c     print *,iyr,imo_n,iday_n,ihour_n

        CALL CLOCK(IYEAR,IMON,IDAY,IHOUR,IMIN,FILENAME)


c       fname = directory_fields(1:ldf)//filename        
        write(indf,'(i3)') index_f
        if (index_f .lt. 10) indf(1:2) = '00'
        if (index_f .lt. 100) indf(1:1) = '0'
        index_f = index_f + met_freq

        fname = directory_fields(1:ldf)//fname1//indf

        write(*,*) 'OPENING FILE:',fname
        open(40,file=fname(1:ldf+20),form='unformatted',
     +                              status='old',err=1000)

        do 450 k = 1,nfields
        call read2d_40(cdate,nx,ny,aux,cvar,iproc)


c --- if reading for any reason is not correct ...exit
         if (iproc .eq. 0) goto 1000

c ---- initial check if fields.00met_freq have the correct date.
 

666     continue

         do i = 1,nx
         do j = 1,ny
         store_1(i,j,k) = aux(i,j)
         end do
         end do


450     continue

c ---------------------------------------------------------

      call met_clock(year,time,time0,iyy,imo,iday,ihour,imin,
     + isec)

        write(sfile,'(i2,i2,i2,i2)') iyy,imo,iday,ihour
        if (iyy.lt.10) sfile(1:1) = '0'
        if (imo.lt.10) sfile(3:3) = '0'
        if (iday.lt.10) sfile(5:5) = '0'
        if (ihour.lt.10) sfile(7:7) = '0'

        sfile(9:12) = '.dat'

        open(50,file=sfile,form='unformatted')
        call write2d_50(nx,ny,nfields,store_1,amask,uc,vc,
     1  rwest,rsout,deltal,deltaf)
        close(50)

        time = time + (met_freq * 3600.)/86400.


        close(40)

        ihour = ihour + met_freq
500     continue


        close(50)


c--------------------------------------------------------------------


        goto 1500
        model_run = 1

c --------------------------------------------------------------------


1000     continue

c ---- come here if errors occur
         open(44,file='norun1')
         close(44)

1500     continue

         stop
         end


        subroutine read2d_40(cdate,nx,ny,arr,cvar,iproc)
        DIMENSION ARR(NX,NY)
        CHARACTER CDATE*11, CVAR*5
        character*3 kechar

        iproc = 0
        read(40,err=800) cdate, cvar
        print *,cdate,cvar
        read(40,err=800) ((arr(i,j), i = 1, nx),j = 1, ny)
        iproc = 1
800     continue

        RETURN
        END
 


        subroutine write2d_50(nx,ny,nf,arr,amask,uc,vc,
     1   wlon,slat,dx,dy)
        parameter ( xland = 9999.)
        DIMENSION ARR(NX,NY,NF),UC(NX,NY),VC(NX,NY),AMASK(NX,NY)


        do i = 1,nx
        do j = 1,ny
        uc(i,j) = arr(i,j,1)
        vc(i,j) = arr(i,j,2)
c----------------------------------------
c----------------------------------------
        if (amask(i,j).eq.0.) then
        uc(i,j) = xland
        vc(i,j) = xland
        end if
        end do
        end do


        write(50) uc,vc

        RETURN
        END


        
      
      SUBROUTINE TRANSLATE(TITLE,IYR_N,IMO_N,IDAY_N)
      CHARACTER TITLE*10,SYM(0:9)*1
      CHARACTER*1 AUX
      DIMENSION IHOLD(3)
      INTEGER QUANT(4)
      DATA QUANT/1,10,100,1000/
      DATA SYM/'0','1','2','3','4','5','6','7','8','9'/

      K = 1
      L = 1
        DO N = 1,3
        IHOLD(N) = 0
        END DO
      DO 200 I = 10,1,-1
      AUX = TITLE(I:I)
      IF (AUX.EQ.'/') THEN
        L = 1
        K = K+1
        GOTO 200
      END IF
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 100
       END IF
       END DO
100    CONTINUE
       IHOLD(K) = IHOLD(K) + MULT*QUANT(L)
       L = L + 1
200    CONTINUE

       IYR_N = IHOLD(1)
       IMO_N = IHOLD(2)
       IDAY_N = IHOLD(3)

c      print *,iyr_n,imo_n,iday_n
       RETURN
       END


      SUBROUTINE TRANSLATE_MET(TITLE,IYR,IMO,IDAY,IHOUR)
      CHARACTER*16 TITLE,SYM(0:9)*1
      CHARACTER*1 AUX
      DIMENSION IHOLD(3)
      INTEGER QUANT(2)
      DATA QUANT/10,1/
      DATA SYM/'0','1','2','3','4','5','6','7','8','9'/


c --- first read which day
      L = 1
      IDAY = 0
      DO 200 I = 1,2
      AUX = TITLE(I:I)
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 100
       ELSEIF (AUX.EQ.'') THEN
       MULT = 0
       END IF
       END DO
100    CONTINUE
       IDAY = IDAY + MULT*QUANT(L)
       L = L + 1
200    CONTINUE
c --- then which month
      L = 1
      IMO = 0
      DO 300 I = 4,5
      AUX = TITLE(I:I)
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 250
       ELSEIF (AUX.EQ.'') THEN
       MULT = 0
       END IF
       END DO
250    CONTINUE
       IMO = IMO  + MULT*QUANT(L)
       L = L + 1
300    CONTINUE

c --- then which year
      L = 1
      IYR = 0
      DO 400 I = 7,8
      AUX = TITLE(I:I)
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 350
       ELSEIF (AUX.EQ.'') THEN
       MULT = 0
       END IF
       END DO
350    CONTINUE
       IYR = IYR  + MULT*QUANT(L)
       L = L + 1
400    CONTINUE

c --- then which hour
      L = 1
      IHOUR = 0
      DO 500 I = 10,11
      AUX = TITLE(I:I)
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 450
       ELSEIF (AUX.EQ.'') THEN
       MULT = 0
       END IF
       END DO
450    CONTINUE
       IHOUR = IHOUR  + MULT*QUANT(L)
       L = L + 1
500    CONTINUE
      L = 1
      DO 600 I = 14,15
      AUX = TITLE(I:I)
       DO M = 0,9
       IF (AUX.EQ.SYM(M)) THEN
       MULT = M
       GOTO 550
       ELSEIF (AUX.EQ.'') THEN
       MULT = 0
       END IF
       END DO
550    CONTINUE
       IHOUR = IHOUR  + MULT*QUANT(L)
       L = L + 1
600    CONTINUE



        RETURN
        END


      subroutine findtime0(time0,year,iihour,iiday,iimo,iiyr)

      real*8 time0
c -- Univ. of Athens - Ocean
      dimension month_o(12),month_lp(12),month(0:12)
      data month_o/31,59,90,120,151,181,212,243,273,304,334,365/
      data month_lp/31,60,91,121,152,182,213,244,274,305,335,366/


      year=float(iiyr)

      year_days = 365.
      if (mod(year,4.).eq.0.) year_days = 366.

      month(0) = 0
      do k = 1,12
      month(k) = month_o(k)
      end do

      if (year_days.eq.366.) then
       do k = 1,12
       month(k) = month_lp(k)
       end do
      end if
C
C
      time0=float(month(iimo-1)+iiday-1) + float(iihour)/24.

      return
      end

      subroutine met_clock(year,time,time0,iyr,imo,iday,
     +     ihour,imin,isec)

c --- Univ. of Athens - Ocean

      real*8 time,time0,hh,xmm,xsec
      dimension month_o(12),month_lp(12),month(0:12),m_days(12)
      dimension m_days_o(12),m_days_lp(12)
      data month_o/31,59,90,120,151,181,212,243,273,304,334,365/
      data month_lp/31,60,91,121,152,182,213,244,274,305,335,366/
      data m_days_o/31,28,31,30,31,30,31,31,30,31,30,31/
      data m_days_lp/31,29,31,30,31,30,31,31,30,31,30,31/
      DAYI=1.E0/86400.E0
  

      iadd_yr = 0
      year_days = 365.
      if (mod(year,4.).eq.0.) year_days = 366.


      month(0) = 0
      do k = 1,12
      month(k) = month_o(k)
      m_days(k) = m_days_o(k)
      end do
      if (year_days.eq.366.) then
       do k = 1,12
       month(k) = month_lp(k)
       m_days(k) = m_days_lp(k)
       end do
      end if


       sub_days = 0
c     if (time.ge.year_days+1) then
      if (time.ge.year_days) then
       sub_days = year_days
       iadd_yr = 1
      end if
      dd = mod(time,year_days)

      hh = dd - int(dd)
      hh = hh * 24.
      xmm = (time-sub_days)/dayi - int(dd)*86400.  - int(hh)*60.*60.
      xmm = xmm/60.
      if (xmm.ge.60.) then
      hh = hh + 1.
      xmm = 0.
      end if

      xsec = (xmm - int(xmm))*60.

      time_elaps = time - time0
      hh_p = time_elaps * 24
      ihour_p = int(hh_p)

      iday = int(dd) + 1
      ihour = int(hh)
      imin = int(xmm)
      isec = int(xsec)

      do k = 1,12
      if (iday.le.month(k)) then
      imo = k
      goto 6000
      end if
      end do

6000  continue
      iyr = year - 1900. + iadd_yr
      if (year.ge.2000.) iyr = year - 2000. + iadd_yr
      iday = iday - month(k-1)
      if (iyr.eq.100) iyr = 0

      write(*,600) time,iday,imo,iyr,ihour,imin,isec

600   format(f7.2,3x,i2,'/',i2,'/',i2,3x,i4,i4,i4)



      return
      end

      function gust(airsea,u10)
      grav = 9.81
      ZWND = 10.
      STABSH =    1.4
      STABOF =   -0.01
      CNEG   =   -0.1
      CPOS   =    0.1
      FNEG   = -150.
      FPOS   =  150.
             STAB0  = ZWND * GRAV / 273.
c
               STAB   = STAB0 * AIRSEA / MAX(5.,U10)**2
               STAB   = MAX ( -1. , MIN ( 1. , STAB ) )
c
               THARG1 = MAX ( 0. , FNEG*(STAB-STABOF))
               THARG2 = MAX ( 0. , FPOS*(STAB-STABOF))
               COR1   = CNEG * TANH(THARG1)
               COR2   = CPOS * TANH(THARG2)
c
               ASF = SQRT ( (1.+COR1+COR2)/STABSH )
               gust = 1./ASF

        return
        end

      SUBROUTINE CLOCK(IYEAR,IMON,IDAY,IHOUR,IMIN,FILENAME)
      CHARACTER FILENAME*20
      DIMENSION MONTHDAYS(12)

       IF (IHOUR.GE.24) THEN
       IHOUR = IHOUR - 24
       IDAY = IDAY + 1
       END IF
        CALL LEAP_YEAR(IYEAR,MONTHDAYS)
       IF (IDAY.GT.MONTHDAYS(IMON)) THEN
        IDAY = 1
        IMON = IMON + 1
       END IF
       IF (IMON.EQ.13) THEN
       IMON = 1
       IYEAR = IYEAR + 1
       END IF

c --- LOCAL TIME CONVERTED TO UTC
 

c ---- CONSTRUCT FILE NAME
       WRITE(FILENAME,'(6x,i4,1x,i2,1x,i2,1x,i2)') IYEAR,IDAY,IMON,IHOUR
       IF (IMON.LT.10) FILENAME(15:15)='0'
       IF (IDAY.LT.10) FILENAME(12:12)='0'
       IF (IHOUR.LT.10) FILENAME(18:18)='0'
       FILENAME(1:6) = 'wind1-'
       FILENAME(11:11)='-'
       FILENAME(14:14)='-'
       FILENAME(17:17)='-'
       FILENAME(20:20)='h'


      RETURN
      END
      SUBROUTINE LEAP_YEAR(IYEAR,MONTHDAYS)
      dimension m_days_o(12),m_days_lp(12),monthdays(12)

      data m_days_o/31,28,31,30,31,30,31,31,30,31,30,31/
      data m_days_lp/31,29,31,30,31,30,31,31,30,31,30,31/

      do k = 1,12
      monthdays(k) = m_days_o(k)
      end do

      if (mod(iyear,4).eq.0) then
      do k = 1,12
      monthdays(k) = m_days_lp(k)
      end do
      endif

      return
      end

