!***************************************************
!
!      subroutine tide_harmcompo
!
!***************************************************
!
!     objet    : corps de la subroutine de calcul des niveaux
!                de la surface libre a la limite ouverte (marine)
!                pour un modele de rang 0 par
!                composition harmonique. Cette composition harmonique
!                basee sur les ondes du spectre de maree. Il y a 6 possibilites
!                de selection ; 3 sont specifiques a des configurations de
!                modeles :
!                 - key_mdn : modele de jcs Manche-Mer du Nord (11 km)
!                 - key_gdg  : modele de pl golfe de Gascogne (11 km)
!                 - par defaut: modele de pg Islande-Portugal  (5 km)
!                Deux donnent la possibilite de couvrir tout le monde :
!                 - key_tide_fes2004
!                 - key_tide_fes2012
!                Une donne la possibilite de couvrir tout le plateau
!                Manche Atlantique Francais :
!                 - 1
!                Le spectre peut etre reduit a la seule onde M2 (tide_nbharmcp=1)
!                et le cycle Morte-Eau Vive-Eau reproduit a l aide des
!                coefficients de maree (eventuellement constant pour
!                simuler des marees schematiques).
!
!     appelant         :
!     appeles          :
!     canal d ecriture :
!     canal de lecture :
!                      :
!     cle cpp          : 1, key_tide_fes2004
!                      :
!     variables in     :
!                      :
!     variables out    :
!                      :
!     variables inout  :
!                      :
!     dernieres modifs : 14/06/2004
!&E
!&E ** History :
!&E       !          Original code
!&E       !  2009-03 (V. Garnier)  removing of xeuv argument
!&E       !  2009-10 (A. Fiandrino)  add harmonic composition for lagoon
!&E
!                      :
!***********************************************************
!


!***************************************************
!
       subroutine tide_harmcompo
!
!***************************************************
!
!     objet    : calcul des niveaux de la surface libre a la
!                limite pour un modele de rang 0 par composition
!                harmonique (type Simon). Composition harmonique
!                basee sur 8 ondes (M2, N2, S2, K2, O1, K1, P1, Q1).
!
!     appelant         : main
!     appeles          : tide_astroarg, tide_calcoef
!     canal d ecriture :
!     canal de lecture : 7
!                      :
!     cle cpp          :
!                      :
!     variables in     : t, h0, methom, parcoef, omegat
!                      :
!     variables out    :
!                      :
!     variables inout  :
!&E
!&E ** History :
!&E       !          Original code
!&E       !  2009-03 (V. Garnier)  removing of xeuv argument
!&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E       !  2013-06 (S. Petton, F. Dumas) Restart : initialisation of the water depth
!&E                                        one day before to make sure all the time steps
!&E                                        used for the temporal initialisation are correct
!&E
!
!*********************************************************************
!
USE parameters
USE ionc4, only : ionc4_openr,ionc4_corres,ionc4_flimax,ionc4_fljmax,&
                  ionc4_read_lat,ionc4_read_lon,ionc4_read_subxy,ionc4_read_subxy
use comvars2d
use commar
use comaree
implicit none
!
      integer,parameter                   :: nbfes=114
      character(len=100)                  :: nomlu
      integer                             :: i,j,n,nuu,l,tide_nbharmcp2,kount,nifes,njfes,k
      integer                             :: imarmin,imarmax,jmarmin,jmarmax,ixmar,jxmar,cdf_id
      real(kind=rlg)                      :: xt,g0fes,phi0fes,dgfes,dphifes
      real(kind=rlg),dimension(:),pointer :: lonmar,latmar
      real(kind=rlg)                      :: long,latit
      real(kind=riosh),dimension(:,:),pointer       :: wra,wrg
      real(kind=rsh),dimension(:,:,:),pointer       :: wrp,wip
      real(kind=riosh),parameter                    :: valmanq_tide = 999.
      CHARACTER(LEN=19)  :: tool_sectodat 
!
      common /fes99/ g0fes,phi0fes,dgfes,dphifes,wrp,wip
      common /fes99/ imarmin,imarmax,jmarmin,jmarmax
!
!************************************************************
!     Executable statements
!************************************************************
!
!   initialisation
!---------------------
!
      if(nbouc==0) then
       tide_nbharmcp=nbfes
       l_presence(:)=.false.
       nu(:)=int(valmanq)
       a0sud(:,:)=valmanq
       ph0sud(:,:)=valmanq
       a0nord(:,:)=valmanq
       ph0nord(:,:)=valmanq
       a0ouest(:,:)=valmanq
       ph0ouest(:,:)=valmanq
       a0est(:,:)=valmanq
       ph0est(:,:)=valmanq
!
! initialisation du tableau des constituants harmoniques
!-------------------------------------------------
!
       call tide_data
       do i=2,nbfes+1
         nu(i-1)=i
       enddo
        l_presence(4:115)=.true.
!
!
! lecture des conditions limites en phase et en amplitude
!--------------------------------------------------------
!
        IF (l_tide_surge) file_tide_harmcp = file_tide_harmcp_surge
        call ionc4_openr(file_tide_harmcp)
        call ionc4_corres(file_tide_harmcp,cdf_id)
        ixmar=ionc4_flimax(cdf_id)
        jxmar=ionc4_fljmax(cdf_id)
        allocate(latmar(1:jxmar),lonmar(1:ixmar))
        call ionc4_read_lat(file_tide_harmcp,latmar,1,jxmar)
        call ionc4_read_lon(file_tide_harmcp,lonmar,1,ixmar)
!
        do i=1,ixmar-1
        if(gwac >= lonmar(i) .and. gwac < lonmar(i+1)) imarmin=i

        if(geac >= lonmar(i) .and. geac < lonmar(i+1)) imarmax=i+1	

       enddo
!
       do j=1,jxmar-1
        if(fisac >= latmar(j) .and. fisac < latmar(j+1)) jmarmin=j

        if(finac >= latmar(j) .and. finac < latmar(j+1)) jmarmax=j+1

       enddo
!
       print*,'TIDES',imarmin,imarmax,jmarmin,jmarmax
       allocate(wrp(imarmin:imarmax,jmarmin:jmarmax,nbfes))
       allocate(wip(imarmin:imarmax,jmarmin:jmarmax,nbfes))
       allocate(wra(imarmin:imarmax,jmarmin:jmarmax))
       allocate(wrg(imarmin:imarmax,jmarmin:jmarmax))
!
       do n=1,nbfes
       nomlu='amp_'//trim(nommar(nu(n)))
       call ionc4_read_subxy(file_tide_harmcp,nomlu,wra,&
	                      imarmin,imarmax,jmarmin,jmarmax,1,1)
       nomlu='phas_'//trim(nommar(nu(n)))
       call ionc4_read_subxy(file_tide_harmcp,nomlu,wrg,&
                      imarmin,imarmax,jmarmin,jmarmax,1,1)
         do j=jmarmin,jmarmax
          do i=imarmin,imarmax
           if (wra(i,j) /= valmanq_tide .and. wrg(i,j) /= valmanq_tide) then
            wrp(i,j,n)= wra(i,j)*cos(-wrg(i,j)*pi/180)
            wip(i,j,n)=-wra(i,j)*sin(-wrg(i,j)*pi/180)
           else
            wrp(i,j,n)=valmanq
            wip(i,j,n)=valmanq
           endif
          enddo
         enddo
       end do
!
        g0fes=lonmar(imarmin)
	phi0fes=latmar(jmarmin)
	dgfes=lonmar(2)-lonmar(1)
	dphifes=latmar(2)-latmar(1)
!
!
! interpolation sur les points des frontieres ouvertes
!--------------------------------------------------------
!
!do j=jmin,jmax
        do j=7,160

! ce qui est vrai pour le curvi devrait l etre aussi pour le cartesien, on peut
! donc enlever ce qui est apres le else
         long=lon2d_g(imin+31,j)
         latit=lat2d_g(imin+31,j)
	 if (long>valmanq_curv) &

         call interpolfes(long,latit,nbfes,a0ouest(1,j),ph0ouest(1,j))
!#ifdef 1
!         long=lon2d_g(imax,j)
!         latit=lat2d_g(imax,j)
!	 if (long>valmanq_curv) &
!#else
!         long=gwac+(imax-1)*dgac
!         latit=fisac+(j-1)*dfiac
!#endif
!         call interpolfes(long,latit,nbfes,a0est(1,j),ph0est(1,j))
        end do
!!
!        do i=imin,imax
!#ifdef 1
!         long=lon2d_g(i,jmin+1)
!         latit=lat2d_g(i,jmin+1)
!	 if (long>valmanq_curv) &
!#else
!         long=gwac+(i-1)*dgac
!         latit=fisac
!#endif
!         call interpolfes(long,latit,nbfes,a0sud(1,i),ph0sud(1,i))
!#ifdef 1
!         long=lon2d_g(i,jmax)
!         latit=lat2d_g(i,jmax)
!	 if (long>valmanq_curv) &
!#else
!         long=gwac+(i-1)*dgac
!         latit=fisac+(jmax-1)*dfiac
!#endif
!         call interpolfes(long,latit,nbfes,a0nord(1,i),ph0nord(1,i))
!        end do
      deallocate(wip,wrp,latmar,lonmar)
!
!
! on reordonne les ondes pour decrire le spectre de facon compatible
! avec la liste SIMON
!-------------------------------------------------
!
      kount=1
      do while(kount/=0)
      kount=0
      do n=1,tide_nbharmcp-1
       if(nu(n) > nu(n+1)) then
        kount=kount+1
        xt=nu(n)
        nu(n)=nu(n+1)
        nu(n+1)=xt
       end if
      end do
      end do
!
      if (l_tide_M2harm) then  ! On ne conserve que l onde M2
         tide_nbharmcp=1 
         nu(1) = im2
      end if
!
      end if

      if (l_init_restart .and. nbouc == 0) then
! initialisation des hauteurs reduites pour J-1
         cdate = tool_sectodat(t-86400.0_rsh)
         CALL tool_decompdate(cdate,ijour,imois,ian,iheure,iminu,isec)
         call tide_predic
         cdate = tool_sectodat(t)
         CALL tool_decompdate(cdate,ijour,imois,ian,iheure,iminu,isec)
      endif
      call tide_predic
!
      end

!****************************************************************************
!
       subroutine interpolfes(long,latit,nbfes,a0,ph0)
!
!***************************************************************************
!
!     objet    : calcul des niveaux de la surface libre a la
!                limite pour un modele de rang 0 par composition
!                harmonique (type Simon). Composition harmonique
!                basee sur 8 ondes (M2, N2, S2, K2, O1, K1, P1, Q1).
!
!     appelant         : main
!     appeles          : tide_astroarg, tide_calcoef
!     canal d ecriture :
!     canal de lecture : 7
!                      :
!     cle cpp          :
!                      :
!     variables in     : t, h0, methom, parcoef, omegat
!                      :
!     variables out    :
!                      :
!     variables inout  :
!&E ** History :
!&E       !  19XX (F. Dumas) Original code
!&E       !  2011-10 (F. Dumas) Integration of a bilinear interpolation for high resolution grids
!&E                  1
!&E                  The two types of interpolation are not compatible yet
!
!*********************************************************************
!
USE parameters
use commar
implicit none
!
      integer                                :: nbfes,i0,j0,i,j,n
      integer                                :: imarmin,imarmax,jmarmin,jmarmax
      real(kind=rlg)                         :: ptot,g0fes,phi0fes,dgfes,dphifes,long,latit,pds

      real(kind=rlg)                         :: var1,var2,xm,ym,xdp,ydp

      real(kind=rsh),dimension(nbfes)                  :: sre,sim
      real(kind=rsh),dimension(nmarmax)                :: a0,ph0
      real(kind=rsh),dimension(:,:,:),pointer          :: wrp,wip
     
      common /fes99/ g0fes,phi0fes,dgfes,dphifes,wrp,wip
      common /fes99/ imarmin,imarmax,jmarmin,jmarmax
!
!     Executable statements
!************************************************************
!
! Assignement of sre and sim

       sre(:)=0.0_rsh
       sim(:)=0.0_rsh
!
       i0=int((long-g0fes)/dgfes)+imarmin
       j0=int((latit-phi0fes)/dphifes)+jmarmin


       xm=(long-g0fes)/dgfes+imarmin
       ym=(latit-phi0fes)/dphifes+jmarmin
       xdp=xm-i0
       ydp=ym-j0

       do n=1,nbfes
        var1=wrp(i0+1,j0  ,n)*xdp+(1.0-xdp)*wrp(i0,j0  ,n)
        var2=wrp(i0+1,j0+1,n)*xdp+(1.0-xdp)*wrp(i0,j0+1,n)
        sre(n)=var2*ydp+(1.0-ydp)*var1
        var1=wip(i0+1,j0  ,n)*xdp+(1.0-xdp)*wip(i0,j0  ,n)
        var2=wip(i0+1,j0+1,n)*xdp+(1.0-xdp)*wip(i0,j0+1,n)
        sim(n)=var2*ydp+(1.0-ydp)*var1
       enddo

!
       do n=1,nbfes

         a0(nu(n))=sqrt(sre(n)*sre(n)+sim(n)*sim(n))
         if(sre(n)==0.0_rsh .and. sim(n)==0.0_rsh) then
          ph0(nu(n))=0.0_rsh
         else
          ph0(nu(n))=atan2(sim(n),sre(n))*180.0_rsh/pi
         end if
         if(ph0(nu(n)) < 0.0_rsh) ph0(nu(n))=ph0(nu(n))+360.0_rsh

       enddo   
!
      end

