      program read_and_interpolate_winds

!      include 'ecommons.h'
      include "../../etafcst_all/parmeta"
      include "../../etafcst_all/wave.inc"

!      parameter (rwest  = -8.0, reast  = 42.0)
!      parameter (rsout  = 29.0, rnort  = 48.0)
!      parameter (deltal =   .05)
!      parameter (deltaf =   .05)
!      parameter (nxp = (reast-rwest)/deltal+1)
!      parameter (nyp = (rnort-rsout)/deltaf+1)
!      parameter (init_gdsdir="../../data/prep/")
!      Send the output file to /vol1/meteo/wew.v0.1/wave/wam_in_out/Wind/nest1
      integer*4 nx,ny,n
      integer*4 imc,jmc,lmc
      integer gds(200)
      real awest, aeast, anorth, asouth, aresx, aresy
      character*64 init_gdsdir, indir, outdir, fname, infile, outfile
      character*64 eta_init
      character*11 indate
      character*5 cvar
c
c *** Input data variables.
c
      real*4,allocatable:: uw(:,:), vw(:,:), asm(:,:)
      real*4,allocatable:: uwam(:,:), vwam(:,:), smwam(:,:)
      real*4,allocatable:: ut(:,:), vt(:,:)
      real*4,allocatable:: sm(:,:), sice(:,:), sst(:,:)
      real*4,allocatable:: deta(:),aeta(:),eta1(:)
      real*4,allocatable:: ull(:,:), vll(:,:), smll(:,:)
c
!      common /wamgeometry/ rwest,reast,rnort,rsout,deltal,deltaf
      common /etageometry/ wbd,sbd,tlm0d,tph0d,dlmd,dphd
!      call eta_commons
c
!      do n=1,ninit
!         call get_sector_size(init_in(n))
!        write(6,*) 'call get_sector_size'
!        write(6,*) 'file= ', init_gdsdir
!      enddo
!      stop
      if (iargc().LT.3) THEN
         print*,'Usage:interp_coupling.x <INDIR> <OUTDIR> <yymmddhhfhr>'
         stop
      endif
c
      call getarg (1,indir)
      call getarg (2,outdir)
      call getarg (3,indate)
      infile=indir(1:index(indir,' ')-1)//'/'//indate(1:11)//'.ecm_wam'
      print*,infile
c
c *** Horizontally interpolate input data from native grid to ETA grid.
c
      init_gdsdir="../../../../data/prep/"    !change it in case of makefile
      n=index(init_gdsdir,' ')-1

      write(6,*) 'trying to open ',init_gdsdir(1:n)//"gdsinfo.ETA_ecm"

      open(unit=14,file=init_gdsdir(1:n)//"gdsinfo.ETA_ecm",
     &       form='unformatted', access='sequential')
      rewind 14
      read(14) GDS
      write(6,*) 'gds(1-14): ', (gds(I),I=1,14)
      close(14)

      nx=GDS(2)
      ny=GDS(3)
      asouth=GDS(4)/1000.
      awest=GDS(5)/1000.
      anorth=GDS(7)/1000.
      aeast=GDS(8)/1000.
      aresx=GDS(9)/1000.
      aresy=GDS(10)/1000.
      allocate (uw(nx,ny))
      allocate (vw(nx,ny))
      allocate (asm(nx,ny))
      print*,nx,ny,asouth,awest,anorth,aeast,aresx,aresy
!      nxp = (reast-rwest)/deltal+1
!      nyp = (rnort-rsout)/deltaf+1
      allocate (uwam(nxp,nyp))
      allocate (vwam(nxp,nyp))
      allocate (smwam(nxp,nyp))
      io = abs((awest-rwest))/aresx+1
      jo = abs((aeast-reast))/aresy+1
      iecm_to_iwam = aresx/deltal
!      print*,io,jo

!Mar      open (unit=15,file=infile,status='unknown',
!Mar     &      form='unformatted')
!Mar      read(15) ((uw(i,j), i=1,nx), j=1,ny)
!Mar      read(15) ((vw(i,j), i=1,nx), j=1,ny)
!Mar      read(15) ((asm(i,j), i=1,nx), j=1,ny)
!Mar      close(15)
c
!Mar      call interp2d (nx,ny,1,io,jo,iecm_to_iwam,nxp,nyp,uw,uwam)
!Mar      call interp2d (nx,ny,1,io,jo,iecm_to_iwam,nxp,nyp,vw,vwam)
!Mar      call interp2d (nx,ny,1,io,jo,iecm_to_iwam,nxp,nyp,asm,smwam)
c
!Mar      call minmax(nx,ny,uw,"Uecmwf:     ")
!Mar      call minmax(nxp,nyp,uwam,"Uecmwf->wam:")
!Mar      call minmax(nx,ny,vw,"Vecmwf:     ")
!Mar      call minmax(nxp,nyp,vwam,"Vecmwf->wam:")
c     WAM init file obtained from large scale model (ECMWF/GFS)
!      outfile=outdir(1:index(outdir,' ')-1)//'/'//indate(1:11)//'.wam'
!      print*,'Outfile: ',outfile
!      open (unit=16,file=outfile,status='unknown',
!     &      form='unformatted')
!      cvar(1:5)='u10  '
!      write(16) indate,cvar
!      write(16) ((uwam(i,j), i=1,nxp), j=1,nyp)
!      cvar(1:5)='v10  '
!      write(16) indate,cvar
!      write(16) ((vwam(i,j), i=1,nxp), j=1,nyp)
!      cvar(1:5)='smask'
!      write(16) indate,cvar
!      write(16) ((smwam(i,j), i=1,nxp), j=1,nyp)
!      close(16)
c
c *** Horizontally interpolate input data from native ETA grid to lat-lon grid.
c
      eta_init="../../../../eta/runs/"    !change it in case of makefile
      infile=eta_init(1:index(eta_init,' ')-1)//'/coupling.eta.file'
      print*, 'Access ',infile
      open (unit=17,file=infile,status='unknown',
     &      form='unformatted')
      read(17) imc,jmc,lmc
      print*, 'IM, JM, LM ',im,jm,lm
      allocate (ut(im,jm))
      allocate (vt(im,jm))
      allocate (sm(im,jm))
      allocate (sice(im,jm))
      allocate (sst(im,jm))
      allocate (deta(lm))
      allocate (eta1(lm-1))
      read(17) ut
      read(17) vt
      read(17) sm  !1-sea
      read(17) sice
      read(17) sst
      read(17) wbd,sbd,tlm0d,tph0d,dlmd,dphd,tlm0d
     .         ,deta,aeta,eta
      close(17)
!      print*,wbd,sbd,tlm0d,tph0d,dlmd,dphd,tlm0d
!      print*, deta,aeta,eta1
      print*,'NXP, NYP ',nxp,nyp

      allocate (ull(nxp,nyp))
      allocate (vll(nxp,nyp))
      allocate (smll(nxp,nyp))
      call imjm2ll(ut,vt,sm,ull,vll,smll)
!      print*,smll
      do i=1,nxp
      do j=1,nyp
         if (smll(i,j).ge.0.5) smll(i,j)=1.   !sea
         if (smll(i,j).lt.0.5) smll(i,j)=0.   !land
      enddo
      enddo
!      print*,smll
C
      call minmax(im,jm,ut,"Ueta:      ")
      call minmax(nxp,nyp,ull,"Ueta->wam:  ")
      call minmax(im,jm,vt,"Veta:      ")
      call minmax(nxp,nyp,vll,"Veta->wam:  ")
      call minmax(nxp,nyp,smll,"SMeta->wam: ")
c     WAM init file obtained from init and const eta files
      outfile=outdir(1:index(outdir,' ')-1)//'/'//indate(1:11)//
     &'.init.wam'
      print*,'Outfile: ',outfile
      open (unit=16,file=outfile,status='unknown',
     &      form='unformatted')
      cvar(1:5)='u10  '
      write(16) indate,cvar
      write(16) ((ull(i,j), i=1,nxp), j=1,nyp)
      cvar(1:5)='v10  '
      write(16) indate,cvar
      write(16) ((vll(i,j), i=1,nxp), j=1,nyp)
      cvar(1:5)='smask'
      write(16) indate,cvar
      write(16) ((smll(i,j), i=1,nxp), j=1,nyp)
      close(16)
      write(17) smll

 100  continue

      stop
      end

      subroutine interp2d(ncol,nrow,nlay,io,jo,nmesh,ncolf,nrowf,
     &                    cval,fval)
c
c     INTERP2D horizontally interpolates a coarse grid field to a
c     fine grid
c
c     Copyright 1996, 1997, 1998, 1999, 2000, 2001 
c     ENVIRON International Corporation
c          
c     Modifications:
c        none
c
c     Input arguments:
c        ncol              number of columns in parent grid
c        nrow              number of rows in parent grid 
c        nlay              number of layers in parent grid
c        io                starting i index for the fine grid
c        jo                starting j index for the fine grid
c        nmesh             mesh number
c        ncolf             number of columns in fine grid
c        nrowf             number of rows in fine grid
c        cval              cell centered value on coarse grid
c
c     Output arguments:
c        fval              cell centered value on fine grid
c
c     Subroutine called:
c        none
c
c     Called by:
c        STARTUP
c        INTRPDAT
c        INTRPCNC
c
      dimension cval(ncol,nrow,nlay),fval(ncolf,nrowf,nlay)
c
c-----Entry point
c
      do 50 k = 1,nlay
        jc = jo
        yc = float(jc - jo + 1) - 0.5
        do j = 1,nrowf
          ic = io
          xc = float(ic - io + 1) - 0.5
          yf = (j - 1.5)/float(nmesh)
          if (yc.le.yf) then
              jc = jc + 1
              yc = float(jc - jo + 1) - 0.5
          endif
          do i = 1,ncolf
            xf = (i - 1.5)/float(nmesh)
            if (xc.le.xf) then
              ic = ic + 1
              xc = float(ic - io + 1) - 0.5
            endif
            dcdx1 = (cval(ic,jc-1,k) - cval(ic-1,jc-1,k))
            dcdx2 = (cval(ic,jc,k) - cval(ic-1,jc,k))
            c1 = cval(ic,jc-1,k) - dcdx1*(xc - xf)
            c2 = cval(ic,jc,k) - dcdx2*(xc - xf)
            dcdy = (c2 - c1)
            fval(i,j,k) = c2 - dcdy*(yc - yf)
          enddo
        enddo
c
  50  continue
c
      return
      end
C
      SUBROUTINE imjm2ll(U,V,SM,ULL,VLL,SMLL)
C-----------------------------------------------------------------------
      INCLUDE "../../etafcst_all/parmeta"
      INCLUDE "../../etafcst_all/wave.inc"
      PARAMETER (IMT=2*IM-1,JMT=JM/2+1,LM1=LM-1,LP1=LM+1)
      PARAMETER (DTR = .01745329)
C-----------------------------------------------------------------------
      DIMENSION U(IM,JM),V(IM,JM),SM(IM,JM)
      DIMENSION ULL(NXP,NYP), VLL(NXP,NYP), SMLL(NXP,NYP)
      DIMENSION ISTH(NXP,NYP),JSTH(NXP,NYP),PPH(NXP,NYP),QQH(NXP,NYP)
C
!      common /wamgeometry/ rwest,reast,rnort,rsout,deltal,deltaf
      common /etageometry/ wbd,sbd,tlm0d,tph0d,dlmd,dphd
C
      print*,wbd,sbd,tlm0d,tph0d,dlmd,dphd
      DO J=NYP,1,-1
      DO I=1,NXP
         ALMD=(I-1)*DELTAL+RWEST
         APHD=(J-1)*DELTAF+RSOUT
         CALL TLL(ALMD,APHD,TLMD,TPHD,TPH0D,TLM0D)
         CALL PQIJSTH(TLMD,TPHD,DLMD,DPHD,WBD,SBD,IM,JM
     &                ,ISTH(I,J),JSTH(I,J),PPH(I,J),QQH(I,J))
      END DO
      END DO
C
      DO J=NYP,1,-1
      DO I=1,NXP
C                       !!!! geop points !!!!!!
         LISTH=ISTH(I,J)
         LJSTH=JSTH(I,J)

         LIWST=LISTH-MOD(LJSTH,2)
         LJWST=LJSTH+1

         LIEST=LISTH+MOD(LJSTH+1,2)
         LJEST=LJSTH+1

         LINRT=LISTH
         LJNRT=LJSTH+2

         IF(LISTH.GT.00.AND.LISTH.LT.IM.
     &      AND.LJSTH.GT.00.AND.LJSTH.LT.JM.
     &      AND.LIWST.GT.00.AND.LIWST.LT.IM.
     &      AND.LJWST.GT.00.AND.LJWST.LT.JM.
     &      AND.LIEST.GT.00.AND.LIEST.LT.IM.
     &      AND.LJEST.GT.00.AND.LJEST.LT.JM.
     &      AND.LINRT.GT.00.AND.LINRT.LT.IM.
     &      AND.LJNRT.GT.00.AND.LJNRT.LT.JM)THEN
            CALL BILIN(U(LISTH,LJSTH),U(LIWST,LJWST)
     &                ,U(LIEST,LJEST),U(LINRT,LJNRT)
     &                ,ULL(I,J),PPH(I,J),QQH(I,J))
            CALL BILIN(V(LISTH,LJSTH),V(LIWST,LJWST)
     &                ,V(LIEST,LJEST),V(LINRT,LJNRT)
     &                ,VLL(I,J),PPH(I,J),QQH(I,J))
            CALL BILIN(SM(LISTH,LJSTH),SM(LIWST,LJWST)
     &                ,SM(LIEST,LJEST),SM(LINRT,LJNRT)
     &                ,SMLL(I,J),PPH(I,J),QQH(I,J))
         ELSE
            print*,'H-points south outside :',LISTH,LJSTH
            print*,'H-points west outside :',LIWST,LJWST
            print*,'H-points north outside :',LINRT,LJNRT
            print*,'H-points east outside :',LIEST,LJEST
            ULL(I,J)=1.E+12
            VLL(I,J)=1.E+12
            SMLL(I,J)=1.E+12
         ENDIF
      END DO
      END DO
      RETURN
      END
C//////////////////////////////////////////////////////////////////
      SUBROUTINE PQIJSTH(TLMPT,TPHPT,DLMD,DPHD,WBD,SBD,IM,JM
     &                  ,ISTH,JSTH,P,Q)
      X1=(TLMPT-WBD)/DLMD
      Y1=(TPHPT-SBD)/DPHD
      X2=.50*( X1+Y1)
      Y2=.50*(-X1+Y1)+(IM-1)
      I2=INT(X2)
      J2=INT(Y2)
      P=X2-I2
      Q=Y2-J2
      JR=J2-(IM-1)
      ISTH=I2-JR+1
      ISTH=(ISTH+1)/2
      JSTH=I2+JR+1
      RETURN
      END
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C-----------------------------------------------------------------------
      SUBROUTINE TLL (ALMD,APHD,TLMD,TPHD,TPH0D,TLM0D)
C
C     Transformation from natural lat/lon coordinates (ALMD,APHD)
C                    to rotated coordinates (TLMD,TPHD) 
C
C-----------------------------------------------------------------------
      DTR=.01745329
C
      TPH0=TPH0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)
C
      RELM=(ALMD-TLM0D)*DTR
      SRLM=SIN(RELM)
      CRLM=COS(RELM)
      APH=APHD*DTR
      SPH=SIN(APH)
      CPH=COS(APH)
      CC=CPH*CRLM
      ANUM=CPH*SRLM
      DENOM=CTPH0*CC+STPH0*SPH
C
      TLMD=ATAN2(ANUM,DENOM)/DTR
      TPHD=ASIN(CTPH0*SPH-STPH0*CC)/DTR
C
      RETURN
      END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE BILIN(H00,H01,H10,H11,HN,P,Q)
C **********************************************************************
C *                                                                    *
C *  routine for bilinear interpolation                                *
C *                                                                    *
C **********************************************************************
      HN=H00+P*(H10-H00)+Q*(H01-H00)+P*Q*(H00-H10-H01+H11)
      RETURN
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE BILIN_S(H00,H01,H10,H11,
     &                   S00,S01,S10,S11,
     &                   HN,P,Q)
C **********************************************************************
C *                                                                    *
C *  routine for bilinear interpolation                                *
C *                                                                    *
C **********************************************************************
      HN=S00*H00+P*(S10*H10-S00*H00)+Q*(S01*H01-S00*H00)
     &          +P*Q*(S00*H00-S10*H10-S01*H01+S11*H11)
      RETURN
      END
C
      subroutine minmax(ii,jj,arin,vname)
      character*12 vname
      real arin(ii,jj)
      xmin=10000.
      xmax=-10000.
      do i=1,ii
      do j=1,jj
         if (arin(i,j).lt.xmin) xmin=arin(i,j)
         if (arin(i,j).gt.xmax) xmax=arin(i,j)
      enddo
      enddo
      print*,vname,xmin,xmax

      return
      end
