!------------------------------------------------------------
!
        subroutine tide_predic
   !&E
   !&E ** History :
   !&E       !          Original code
   !&E       !  2009-03 (V. Garnier)  removing of xeuv argument
   !&E                                initialization of obc_tide_ssh*
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
!
!------------------------------------------------------------

#include "toolcpp.h"
!
USE parameters
use comvars2d
use commar
use cometdyn
use comaree
   USE obccombine, ONLY : obc_tide_sshn,obc_tide_sshs,obc_tide_sshe,obc_tide_sshw
implicit none
!
      integer                   :: i,j,k,n,nol,nd0,ir0,ir1,ir2,nd,jourmem=-1
      integer,save :: nbouc_restart=0
      integer,dimension(8),save :: num
      integer,dimension(144)    :: rr0,rr1,rr2
      integer,dimension(144),PARAMETER :: ide=(/ &
        55555,  56555,  57555,  63655,  65455,  73555,  75555, 125755, &
       127555, 135655, 137455, 145555, 146555, 147555, 153655, 155655, &
       157455, 162556, 163555, 164555, 165555, 166554, 167555, 173655, &
       175455, 183555, 185555, 195455, 209655, 217755, 219555, 225655, &
       225855, 227655, 229455, 235555, 235755, 237555, 238554, 245655, &
       247455, 253555, 253755, 254556, 255555, 256554, 257555, 263655, &
       265455, 265655, 272556, 273555, 274554, 275555, 283455, 285455, &
       291555, 293555, 335655, 345555, 355555, 363555, 364555, 365555, &
       375555, 381555, 382555, 383555, 385555, 417755, 419555, 427655, &
       429455, 435555, 435755, 437555, 445655, 447455, 453555, 455555, &
       457555, 463655, 465455, 465655, 471555, 472556, 473555, 475555, &
       481655, 483455, 485455, 491555, 493555, 625655, 627655, 629455, &
       635555, 635755, 637555, 645655, 647455, 653555, 655555, 657555, &
       663655, 665455, 665655, 672556, 673555, 675555, 681655, 683455, &
       685455, 691555, 693555, 827655, 835755, 845655, 847455, 853555, &
       855555, 857555, 863655, 865455, 865655, 873555, 875555, 881655, &
       883455, 883655, 885455, 891555, 893555,1027655,1035755,1045655, &
      1055555,1063655,1073555,1075555,1083455,1083655,1091555,  85455/)
      real(kind=rsh)                               :: t0,tempis,tide_coef
!
      data rr1/ &
      0,0,0,-72,-66,72,0,189,189,189,188,189,0,0,227,0,0,0,0,0,  20,3*0 &
      ,0,0,0,0,-112,-75,-112,279,0,-75,-75,279,-37,-37,0,-37,-37,260,0,0 &
      ,-37,0,-37,-37,-37,-37,0,2,0,0,-37,-37,0,0,152, 61,-56,189,0,0,0,0 &
      ,0, 20, 20,-112,-112,-112,-112,186,-75,-112,-75,-75,223,-75,-75,-37, &
      -112,-37,261,-37,-37,-37,-37,-75,-75,0,0, &
      148,-149,-149,148,-112,-149,-112,-112,186,-112,-112,-75,-149,-75, &
      -75,-75,-75,-37,-112,-112,-37,-37, &
      -112,-149,-149,-149,148,-149,-149,-112,-149,-112,-112,-112,-75, &
      -149,-75,-149,-75,-75, &
      -224,-186,-186,-186,-149,-149,-149,-186,-112,-112,41/
!
      data rr2/ &
       3*0,-65,-65,-64,414, &
       6*0,218,0,-200,-219,3*0,-135,3*0,-198,-195,-640,-640, &
       17*0,-477,298,2*0,260,3*0,298,-37,260,-37,260, &
       0,20,3 *0,135,439,2*0,-135,-433, &
       11*0,298,0,-37,298,3*0,298,-37,-37,261,0,298, &
       10*0,297,0,-37,297,2*0,297,0,-37,260,0,297, &
       6*0,297,2*0,297,0,297,0,-37,297,260,0,297, &
       6*0,297,-37,297,0,0/
!
      data rr0/ &
       13*1000,4*-1000,3*1000,8*-1000,19*1000,2*-1000,14*1000, &
       2*-1000,2*1000,2*-1000,54*1000,-1000,19*1000,1000/
!
      rr0(53)=-1000
      rr1(54)=32
!
!
!  initialisation
!-----------------
!
      tempis =iheure*3600.0+iminu*60.0+isec
!
      if(nbouc == 0 .and. nbouc_restart==0) then
        nbouc_restart=1
!
! lecture des conditions aux limites
!-----------------------------------
! a faire une bonne fois pour toute
!----------------------------------
!
      do j=jmin,jmax
       num(:)=0
       rouest(:,:,j)=0.0_rsh
       gouest(:,:,j)=0.0_rsh
       rest(:,:,j)=0.0_rsh
       gest(:,:,j)=0.0_rsh
       nomd(:,:)=0
       fr1(:,:)=0.0_rsh
       fr2(:,:)=0.0_rsh
       nd0=0
       k=1
       n=1
       do i=1,tide_nbharmcp
        nol=ide(nu(i))
        ir0=rr0(nu(i))
        ir1=rr1(nu(i))
        ir2=rr2(nu(i))
        nd=nol/100000
        if(nd /= nd0) then
         num(k)=n-1
         k=nd+1
         if(nd >= 6) k=(nd+6)/2
         nd0=nd
         n=1
        end if
!
        rouest(n,k,j)=a0ouest(nu(i),j)*100.0_rsh
        gouest(n,k,j)=ph0ouest(nu(i),j)
        rest(n,k,j)=a0est(nu(i),j)*100.0_rsh
        gest(n,k,j)=ph0est(nu(i),j)
        nomd(n,k)=nol
        if(ir0 < 0) gouest(n,k,j)=gouest(n,k,j)+180.0_rsh
        if(ir0 < 0) gest(n,k,j)=gest(n,k,j)+180.0_rsh
        fr1(n,k)=ir1/float(ir0)
        fr2(n,k)=ir2/float(ir0)
        n=n+1
       end do
       num(k)=n-1
       hnouest(129,j)=-valmanq
       hnest(129,j)=-valmanq
      end do
!
      do i=imin,imax
       num(:)=0
       rnord(:,:,i)=0.0_rsh
       gnord(:,:,i)=0.0_rsh
       rsud(:,:,i)=0.0_rsh
       gsud(:,:,i)=0.0_rsh
       nd0=0
       k=1
       n=1
       do j=1,tide_nbharmcp
        nol=ide(nu(j))
        ir0=rr0(nu(j))
        ir1=rr1(nu(j))
        ir2=rr2(nu(j))
        nd=nol/100000
        if(nd /= nd0) then
         num(k)=n-1
         k=nd+1
         if(nd >= 6)k=(nd+6)/2
         nd0=nd
         n=1
        end if
        rsud(n,k,i)=a0sud(nu(j),i)*100.0_rsh
        gsud(n,k,i)=ph0sud(nu(j),i)
        rnord(n,k,i)=a0nord(nu(j),i)*100.0_rsh
        gnord(n,k,i)=ph0nord(nu(j),i)
        if(ir0 < 0) gsud(n,k,i)=gsud(n,k,i)+180.0_rsh
        if(ir0 < 0) gnord(n,k,i)=gnord(n,k,i)+180.0_rsh
        n=n+1
       end do
       num(k)=n-1
       hnsud(129,i)=-valmanq
       hnnord(129,i)=-valmanq
      end do
!
! choix du time zone
!----------------------
!
!     t0=-tu
!
!
! que represente k a ce stade et pour la suite
! lui donner un nom plus explicite et puis le sauver
!---------------------------------------------
!
!
      end if
!
! calcul des 128 hauteurs par jour
!---------------------------------
!
   t0=0.0_rsh

   if(ijour /= jourmem) then
     do i=liminp1m2,limaxp2
       IF_MPI (ljmin==jmin) THEN
         if(h0(i,jmin+1)/=-valmanq) call ma1(t0,ijour,imois,ian,hnsud(1,i),rsud(1,1,i),gsud(1,1,i),nomd,fr1,fr2,num)
       ENDIF_MPI
       IF_MPI (ljmax==jmax) THEN
         if(h0(i,jmax)/=-valmanq)   call ma1(t0,ijour,imois,ian,hnnord(1,i),rnord(1,1,i),gnord(1,1,i),nomd,fr1,fr2,num)
       ENDIF_MPI
     end do

     do j=ljminp1m2,ljmaxp2        !7,160            !ljminp1,ljmaxp1 !iman
if(j>=7 .and. j<=160) then
       IF_MPI (limin==imin) THEN
         if(h0(imin+31,j)/=-valmanq) call ma1(t0,ijour,imois,ian,hnouest(1,j),rouest(1,1,j),gouest(1,1,j),nomd,fr1,fr2,num)
       ENDIF_MPI
       IF_MPI (limax==imax) THEN
         if(h0(imax,j)/=-valmanq)   call ma1(t0,ijour,imois,ian,hnest(1,j),rest(1,1,j),gest(1,1,j),nomd,fr1,fr2,num)
       ENDIF_MPI
end if
     end do
   end if

! interpolation parmi les 128 hauteurs par jour
!----------------------------------------------


   IF (l_tide_M2harm) CALL tide_calcoef(tide_coef)

   IF_MPI (limin==imin) THEN
     obc_tide_sshs(imin)=0.0_rsh
     obc_tide_sshn(imin)=0.0_rsh
   ENDIF_MPI

   DO i=liminp1m2,limaxp2
     IF_MPI (ljmin==jmin) THEN
       IF (h0(i,jmin+1) /= -valmanq) THEN
         call ma7(tempis,hnsud(1,i),obc_tide_sshs(i),k)
         if(l_tide_M2harm) obc_tide_sshs(i)=obc_tide_sshs(i)*tide_coef
         IF (l_tide_surge) THEN
           obc_tide_sshs(i)=obc_tide_sshs(i)+mslshift
         ELSE
           obc_tide_sshs(i)=obc_tide_sshs(i)+mslshift-(pa(i,jmin+1)-paref)/10000.0_rsh
         END IF
       END IF
     ENDIF_MPI
     IF_MPI (ljmax==jmax) THEN
       IF (h0(i,jmax) /= -valmanq) THEN
         call ma7(tempis,hnnord(1,i),obc_tide_sshn(i),k)
         if(l_tide_M2harm) obc_tide_sshn(i)=obc_tide_sshn(i)*tide_coef
         IF (l_tide_surge) THEN
           obc_tide_sshn(i)=obc_tide_sshn(i)+mslshift
         ELSE
           obc_tide_sshn(i)=obc_tide_sshn(i)+mslshift-(pa(i,jmax)-paref)/10000.0_rsh
         END IF
       END IF
     ENDIF_MPI
   END DO

   IF_MPI (ljmin==jmin) THEN
      obc_tide_sshw(jmin)=0.0_rsh
      obc_tide_sshe(jmin)=0.0_rsh
   ENDIF_MPI

   DO j=ljminp1m2,ljmaxp2 !7,160         !ljminp1,ljmaxp1 !iman
if(j>=7 .and. j<=160) then
     IF_MPI (limin==imin) THEN
       IF (h0(imin+31,j) /= -valmanq) THEN          !iman imin+1
         call ma7(tempis,hnouest(1,j),obc_tide_sshw(j),k)
         if(l_tide_M2harm) obc_tide_sshw(j)=obc_tide_sshw(j)*tide_coef
         IF (l_tide_surge) THEN
           obc_tide_sshw(j)=obc_tide_sshw(j)+mslshift
         ELSE
           obc_tide_sshw(j)=obc_tide_sshw(j)+mslshift-(pa(imin+31,j)-paref)/10000.0_rsh       !iman imin+1!j'ai ajouter10cm a CL
         END IF
       END IF
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
       IF (h0(imax,j) /= -valmanq) THEN
         call ma7(tempis,hnest(1,j),obc_tide_sshe(j),k)
         if(l_tide_M2harm) obc_tide_sshe(j)=obc_tide_sshe(j)*tide_coef
         IF (l_tide_surge) THEN
           obc_tide_sshe(j)=obc_tide_sshe(j)+mslshift
         ELSE
           obc_tide_sshe(j)=obc_tide_sshe(j)+mslshift-(pa(imax,j)-paref)/10000.0_rsh
         END IF
       END IF
     ENDIF_MPI
else
   obc_tide_sshw(j)=0.0
endif
   END DO

   jourmem=ijour

  END SUBROUTINE tide_predic
!
!
!******************************************************
!
! Interpolations entres les 128 hauteurs precalculees
!
      subroutine ma7(seci,hn,haut,k)
!
!******************************************************
!
USE parameters
implicit none
!
 integer                      :: idk,l,k
 real(kind=rsh)               :: d1,haut,fk,d,hh,d2 
 real(kind=rsh),dimension(130):: hn
 real(kind=rsh)               :: seci 
!
      fk=seci/60.0_rsh/11.25_rsh+2.0_rsh
      k=fk+0.5_rsh
      d=fk-k
      if(k >= 130)then
       idk=k-129
       k=k-idk
       d=d+idk
      endif
!
      l=k
      if(l == 130)l=129
      hh=hn(l)+hn(l)
      d1=hn(l+1)-hn(l-1)
      d2=hn(l+1)+hn(l-1)-hh
!
      haut=(hh+d*(d*d2+d1))/800.0_rsh
!
      end subroutine ma7
!
!***********************************************************
!
!  calcul de hn :128 hauteurs de maree par jour
!
      subroutine ma1(t0,jour,mois,ia,hn,rr,gg,nno,fr1,fr2,nam)
!
!***********************************************************
!
USE parameters
implicit none
integer                        :: l,nj,i,julien2,nd0,n,m,jour,mois,ia,kg,ng,j,k
integer,dimension(8)              :: nam
integer,dimension(30)             :: noa
integer,dimension(30,8)           :: nno
real(kind=rsh)                    :: t0,fq,co,so
real(kind=rlg)                    :: fnjd
real(kind=rsh),dimension(30)      :: f,q,q0,r1,r2,v0
real(kind=rsh),dimension(31)      :: g,r
real(kind=rsh),dimension(130)     :: hn
real(kind=rsh),dimension(362)     :: c
real(kind=rsh),dimension(91)      :: cs=(/ &
       1.000000,0.999848,0.999391,0.998630,0.997564,0.996195,0.994522, &
       0.992546,0.990268,0.987688,0.984808,0.981627,0.978148,0.974370, &
       0.970296,0.965926,0.961262,0.956305,0.951057,0.945519,0.939693, &
       0.933580,0.927184,0.920505,0.913545,0.906308,0.898794,0.891007, &
       0.882948,0.874620,0.866025,0.857167,0.848048,0.838671,0.829038, &
       0.819152,0.809017,0.798636,0.788011,0.777146,0.766044,0.754710, &
       0.743145,0.731354,0.719340,0.707107,0.694658,0.681998,0.669131, &
       0.656059,0.642788,0.629320,0.615661,0.601815,0.587785,0.573576, &
       0.559193,0.544639,0.529919,0.515038,0.500000,0.484810,0.469472, &
       0.453991,0.438371,0.422618,0.406737,0.390731,0.374607,0.358368, &
       0.342020,0.325568,0.309017,0.292372,0.275637,0.258819,0.241922, &
       0.224951,0.207912,0.190809,0.173648,0.156434,0.139173,0.121869, &
       0.104528,0.087156,0.069756,0.052336,0.034899,0.017452,0.000000/)
real(kind=rsh),dimension(30,8)    :: rr,gg,fr1,fr2
real(kind=rsh),dimension(11,3)    :: x,y
!
      equivalence (c(1),cs(1))
!
! executable statements
!-------------------------
!
      l=1
!
      x(:,:)=0.0_rsh
      y(:,:)=0.0_rsh
!
      c(271)=0.0_rsh
      do i=1,90
       c(182-i)=-c(i)
       c(362-i)=c(i)
       c(180+i)=-c(i)
      end do
!
      nj=julien2(ia,jour,mois)-2415021
      fnjd=real(nj+l,rlg)-1.4921875_rlg
      nj=mod(nj,7)+1
      if(nj <= 0)nj=nj+7
      nd0=-1
      n=0
      m=1
!
      do kg=1,8
       n=nam(kg)
       if(n /= 0) then
        do i=1,n
         r(i)=rr(i,kg)
         g(i)=gg(i,kg)
         noa(i)=nno(i,kg)
         r1(i)=fr1(i,kg)
         r2(i)=fr2(i,kg)
        end do
!
        call ma3(f,v0,q,n,r,g,fnjd,ng,t0,r1,r2,noa)
!
        do i=1,n
         r(i)=r(i)*f(i)
         q0(i)=v0(i)-g(i)+360.0_rsh
         q(i)=mod(q(i),360.0_rsh)
         q0(i)=mod(q0(i),360.0_rsh)
         if(q0(i) < 0.0_rsh)q0(i)=q0(i)+360.0_rsh
        end do
!
        do i=l,3
         do j=1,n
          fq=q0(j)
          k=fq+1
          fq=fq-k+1
          co=((c(k+1)-c(k))*fq+c(k))*r(j)
          k=k-90
          if(k <= 0) k=k+360
          so=((c(k+1)-c(k))*fq+c(k))*r(j)
          x(ng,i)=x(ng,i)+co
          y(ng,i)=y(ng,i)+so
          q0(j)=q0(j)+q(j)
          if(q0(j) > 360.0_rsh)q0(j)=q0(j)-360.0_rsh
         end do
        end do
       end if
      end do
!
      call ma4(hn,l,x,y)
!
      end subroutine ma1
!
!**********************************************************************
!
!  Calcul des hauteurs reduites h pour les jours jour-1,jour et jour+1
!  et appel de la subroutine de calcul des 128 hauteurs par jour (ma5)
!
      subroutine ma4(hn,l,hm1,hm2)
!
!**********************************************************************
!
USE parameters
implicit none
integer                        :: k,l,i
real(kind=rsh)                 :: fnm
real(kind=rsh),dimension(130)  :: hn
real(kind=rsh),dimension(11,3) :: hm1,hm2
real(kind=rsh),dimension(128,3):: h
complex,dimension(128)         :: x
!
! executable statements
!-------------------------
!
      do  k=l,3
       h(:,k)=0.0_rsh
       x(:)=0.0_rsh
       x(1)=hm1(1,k)
       fnm=x(1)*2.0_rsh
       hm2(1,k)=0.0_rsh
       hm1(1,k)=0.0_rsh
       do i=2,11
        x(i)=cmplx(hm1(i,k),hm2(i,k))
        hm1(i,k)=0.0_rsh
        hm2(i,k)=0.0_rsh
        x(130-i)=conjg(x(i))
       end do
!
       call fft(x)
!
       h(:,k)=x(:)
      end do
!
      call ma5(h,hn)
!
      end subroutine ma4
!
!**********************************************************************
!
!  Calcul des 128 hauteurs par interpolation entre les hauteurs reduites
!
      subroutine ma5(h,x)
!
!**********************************************************************
!
USE parameters
implicit none
integer                        :: k
real(kind=rsh)                 :: b,hh
real(kind=rsh),dimension(128,3):: h
real(kind=rsh),dimension(130)  :: x
!
! executable statements
!-------------------------
!
      b=-0.5078125_rsh
      x(1)=x(129)
      x(2)=x(130)
      do k=1,128
       hh=h(k,2)+h(k,2)
       b=b+0.0078125_rsh
       x(k+2)=((h(k,1)+h(k,3)-hh)*b+h(k,3)-h(k,1))*b+hh
      end do
      if(x(1) == -valmanq)then
       x(2)=3.0_rsh*(x(3)-x(4))+x(5)
       x(1)=6.0_rsh*x(3)-8.0_rsh*x(4)+3.0_rsh*x(5)
      endif
!
      end subroutine ma5
!
!********************************************
!
!  Transformee de Fourier rapide
!
      subroutine fft(a)
!
!********************************************
! 
USE parameters
implicit none
integer                :: m,nv2,nm1,n,i,j,ip,l,k,le,le1
real(kind=rsh)         :: fij
complex,dimension(128) :: a
complex                :: u,w,t
!
! executable statements
!-------------------------
!
      m=7
      nv2=64
      nm1=127
      n=128
      j=1
!
      do  i=1,nm1
      if(i < j) then
       t=a(j)
       a(j)=a(i)
       a(i)=t
      end if
      k=nv2
    6 continue
      if(k < j) then
       j=j-k
       k=k/2
       goto 6
      end if
      j=j+k
      end do
!
      do l=1,m
      le=2**l
      le1=le/2
      u=(1.0_rsh,0.0_rsh)
      fij=pi/le1
      w=cexp(cmplx(0.0_rsh,fij))
      do j=1,le1
       do i=j,n,le
        ip=i+le1
        t=a(ip)*u
        a(ip)=a(i)-t
        a(i)=a(i)+t
       end do
       u=u*w
      end do
      end do
!
      end subroutine fft
!
!**************************************************************
!
!  Calcul des amplitudes et phases des composantes
!
      subroutine ma3(f,v0,q,n,r,g,fnjd,ng,t0,r1,r2,noa)
!
!**************************************************************
!
USE parameters
implicit none
integer                     :: i,ng,n
integer,dimension(6)        :: nd
integer,dimension(30)       :: noa
real(kind=rsh)              :: fa1,fb1,fa2,fb2,a1,a2,asfov,t0
real(kind=rlg)              :: fnjd
real(kind=rsh),dimension(30):: q,v0,f,r,g,r1,r2
real(kind=rlg)              :: dq
complex                     :: fc
!
!----------------------------------------------
!
      do i=1,n
        nd(6)=noa(i)
        nd(5)=noa(i)/10
        nd(4)=noa(i)/100
        nd(3)=noa(i)/1000
        nd(2)=noa(i)/10000
        nd(1)=noa(i)/100000
        nd(6)=nd(6)-nd(5)*10
        nd(5)=nd(5)-nd(4)*10
        nd(4)=nd(4)-nd(3)*10
        nd(3)=nd(3)-nd(2)*10
        nd(2)=nd(2)-nd(1)*10
        ng=nd(1)+1
        if(nd(1) == 0) r(i)=r(i)*2.0_rsh
        call masfo(nd,dq)
        q(i)=dq
        g(i)=g(i)+q(i)*t0/24.0_rsh
        g(i)=mod(g(i),360.0_rsh)
        if(g(i) < 0.0) g(i)=g(i)+360.0_rsh
        v0(i)=asfov(nd)+mod(REAL(dq*fnjd,rlg),360._rlg)
        fc=1.0_rsh
        if( r1(i) /= 0.0_rsh .or. r2(i) /= 0.0_rsh ) then
          fb2=9.242202e-4
          fb1=-fb2
          fa2=1.760045_rsh
          fa1=4.523139_rsh
          if(noa(i) == 275555)then
           fb1=2*fb2
           fa1=2*fa2
          endif
          a1=fa1+fb1*fnjd
          a2=fa2+fb2*fnjd
          fc=1.0_rsh+r1(i)*cmplx(cos(a1),sin(a1))+r2(i)*cmplx(cos(a2),sin(a2))
        end if
!
        v0(i)=v0(i)+atan2(aimag(fc),real(fc))*57.29578_rsh
        f(i)=cabs(fc)
        if(noa(i) == 355555.or.noa(i) == 382555 .or.noa(i) == 164555)  v0(i)=v0(i)-90.0_rsh
      end do
!
      end subroutine ma3
!
!***************************************************
!
! calcul des vitesses angulaires des composantes
!
      subroutine masfo(nd,asfo)
!
!***************************************************
!
USE parameters
implicit none
integer,dimension(6) :: nd
real(kind=rlg)       :: asfo
!
      asfo =(360._rlg-12.19074939_rlg)*(nd(1))+13.17639673_rlg*(nd(2)-5)+0.98564734_rlg*(nd(3)-5) &
             +0.11140408_rlg*(nd(4)-5)+0.05295392_rlg*(nd(5)-5)+0.00004707_rlg*(nd(6)-5)
!
      end subroutine masfo
!
!***************************************************
!
! Calcul des arguments astronomiques
!
      function asfov(nd)
!
!***************************************************
!
USE parameters
implicit none
integer,dimension(6) :: nd
real(kind=rsh)       :: fov,asfov
!
      fov=280.1895_rsh*(nd(1)+nd(3)-5)+277.0248_rsh*(nd(2)-nd(1)-5)+334.3853_rsh*(nd(4)-5) &
         +100.8432_rsh*(nd(5)-5)+281.2209_rsh*(nd(6)-5)+(mod(nd(1),2))*90.0_rsh
      asfov=(mod(fov,360.0_rsh))
!
      end function asfov
!
!**************************************************************
!
! Calcul du nombre de jours ecoules depuis le 1er Janvier 1900
!
      function julien2(ia,jou,moi)
!
!*********************************************************
!
USE parameters
implicit none 
integer                 :: i,julien2,ia,jou,moi,ibs,jour,mois,iy,m
integer,dimension(12,2),PARAMETER :: jo=reshape((/0,31,59,90,120,151,181,212,243,273,304,334,&
                                       0,31,60,91,121,152,182,213,244,274,305,335/),(/12,2/))
real(kind=rsh)          :: a,b
!
      jour=jou
      mois=moi
!
      if(mois==1) then
      ibs=1
      if( ia >= 1582 .and.(ia /= 1582 .or. jou >= 277) ) then
       ibs=mod(ia,4)+2
       if(ibs /= 2)ibs=1
      end if
       do i=1,12
        mois=12-i+1
        jour=jou-jo(mois,ibs)
        if(jour > 0) goto 20
       end do
   20 continue
      end if
!
      a=ia+mois/100.0_rsh+jour/10000.0_rsh
      b=0.0_rsh
      iy=ia
      m=mois
!
      if(mois <= 2)then
        iy=iy-1
        m=m+12
      endif
!
      if(a >= 1582.1015_rsh)then
       b=2-int(iy/100)+int(iy/100)/4
      endif
!
      julien2=int(365.25_rsh*iy)+int(30.6001_rsh*(m+1))+jour+1720995+b
!
      end function julien2
