!!======================================================================

  SUBROUTINE obc2d_rank_default_surge

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE obc2d_rank_default_surge ***
!&E
!&E ** Purpose : Read a sea surface height from an obc file issued from
!&E              a simulation without wind and atmospherical pressure
!&E              Remove this signal at the open boundaries
!&E
!&E ** Description : At open boundaries, the final signal is
!&E                  - the tide issued from a file of harmonic components (SHOM)
!&E                  - plus the influence of the wind and atmospherical pressure
!&E                  (difference of obc issued from two previous ranks,
!&E                  without and whithout wind and Pa)
!&E
!&E                  Straight boundaries case
!&E
!&E ** Called by : obc2drank.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : obc_tide_sshw,obc_tide_sshe,
!&E                         obc_tide_sshs,obc_tide_sshn
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2008-02  (F. Dumas)  Original code
!&E       !  2010-05 (V. Garnier, F. Vandermeirsch) Interpolation of velocities
!&E       !  2010-11 (F. Dumas, S. Theetten) management of child ranks with coarser resolution (imxx,jmxx)
!&E       !          (dirichlet and characteristics methodology
!&E       !  2011-10 (V. Garnier) Evolution required by the spatial extension of ssh(liminm1:limaxp1,ljminm1:ljmaxp1)
!&E       !  2011-10 (V. Garnier) Correction of the inconsistent test on the location on cell (child/mother)
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters, ONLY : valmanq,liminp1m2,limaxm1,limax,ljminp1m2,       &
                          ljmaxm1,ljmax,liminp1,ljminp1
   USE comvars2d
   USE comstatic,  ONLY : imxx,jmxx,mx,my,idh,jhh,xlim,ylim,               &
                          ulim_n,ulim_s,ulim_w,ulim_e, &
                          vlim_n,vlim_s,vlim_w,vlim_e, &
                          hlim_n,hlim_s,hlim_w,hlim_e, &
                          sv1_n,sv1_s,sv2_n,sv2_s,su1_n,su1_s,su2_n,su2_s, &
                          sv1_w,sv1_e,sv2_w,sv2_e,su1_w,su1_e,su2_w,su2_e, &
                          sxe1_n,sxe2_n,sxe1_s,sxe2_s, &
                          sxe1_e,sxe2_e,sxe1_w,sxe2_w, &
                          st1,st2
   USE obccombine, ONLY : obc_tide_sshs,obc_tide_sshn,obc_tide_sshe,obc_tide_sshw, &
                          obc_tide_us,obc_tide_vs,obc_tide_uw,obc_tide_vw,         &
                          obc_tide_un,obc_tide_vn,obc_tide_ue,obc_tide_ve
   USE toolmpi,ONLY : MASTER

   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   INTEGER              :: i,j,l
   REAL(kind=rsh)       :: alf,x,h1,h2,hautd,yy,k1,k2,k3,k4
   REAL(kind=rsh),DIMENSION(imxx+jmxx)           :: an,bn,cn,as,bs,cs
   CHARACTER (len=19) :: tool_sectodat,ccdate

!!----------------------------------------------------------------------
!! * Executable part

   IF(nbouc == 0) THEN

     OPEN(14,file = file_surge_nometeo,form='unformatted')
     READ(14) st1,(sxe1_n(i),sxe1_s(i),su1_n(i),sv1_n(i,1),sv1_n(i,2),su1_s(i),sv1_s(i,1),sv1_s(i,2),i=1,mx)
     READ(14)     (sxe1_w(j),sxe1_e(j),su1_w(1,j),su1_w(2,j),sv1_w(j),su1_e(1,j),su1_e(2,j),sv1_e(j),j=1,my)

     IF(t < st1) THEN
       WRITE(ierrorlog,*) 'ERROR: In obc2d_rank_default_surge.F90,'
       WRITE(ierrorlog,*) 'the simulation starts before the 1st available SSH condition'
       WRITE(ierrorlog,*) 'If this is not rank 0, modify the "date_startshift" parameter'
       WRITE(ierrorlog,*) 'in the "paraspec.txt" file.'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       STOP
     ENDIF

     ccdate=tool_sectodat(st1)
     IF (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'obc2d_rank_default_surge.F90:'
     WRITE(iscreenlog,*) 'The 1st open boundary condition is available at: '
     WRITE(iscreenlog,*) 'time = ',st1,' in sec., corresponding to the date: ',ccdate
     ENDIF
     ccdate=tool_sectodat(t)
     IF (MASTER) &
     WRITE(iscreenlog,*) 'The start date of the simulation is: ',ccdate

     READ(14,END=34,err=35)st2,(sxe2_n(i),sxe2_s(i),su2_n(i),sv2_n(i,1),sv2_n(i,2),su2_s(i),sv2_s(i,1),sv2_s(i,2),i=1,mx)
     READ(14)                 (sxe2_w(j),sxe2_e(j),su2_w(1,j),su2_w(2,j),sv2_w(j),su2_e(1,j),su2_e(2,j),sv2_e(j),j=1,my)

   END IF
!---------------------------
9 CONTINUE                 !     deroulement normal
   IF(t-st1 <= 0.0_rlg)THEN

     ulim_n(1:mx)=su1_n(1:mx)
     ulim_s(1:mx)=su1_s(1:mx)
     ulim_e(1:2,1:my)=su1_e(1:2,1:my)
     ulim_w(1:2,1:my)=su1_w(1:2,1:my)

     vlim_n(1:mx,1:2)=sv1_n(1:mx,1:2)
     vlim_s(1:mx,1:2)=sv1_s(1:mx,1:2)
     vlim_e(1:my)=sv1_e(1:my)
     vlim_w(1:my)=sv1_w(1:my)

     hlim_s(1:mx)=sxe1_s(1:mx)
     hlim_n(1:mx)=sxe1_n(1:mx)
     hlim_w(1:my)=sxe1_w(1:my)
     hlim_e(1:my)=sxe1_e(1:my)

   ELSE

     IF(t-st2 <= 0.0_rlg) THEN
       alf=(st2-t)/(st2-st1)

       ulim_s(1:mx)=alf*su1_s(1:mx)+(1.0_rlg-alf)*su2_s(1:mx)
       ulim_n(1:mx)=alf*su1_n(1:mx)+(1.0_rlg-alf)*su2_n(1:mx)
       ulim_w(1:2,1:my)=alf*su1_w(1:2,1:my)+(1.0_rlg-alf)*su2_w(1:2,1:my)
       ulim_e(1:2,1:my)=alf*su1_e(1:2,1:my)+(1.0_rlg-alf)*su2_e(1:2,1:my)

       vlim_s(1:mx,1:2)=alf*sv1_s(1:mx,1:2)+(1.0_rlg-alf)*sv2_s(1:mx,1:2)
       vlim_n(1:mx,1:2)=alf*sv1_n(1:mx,1:2)+(1.0_rlg-alf)*sv2_n(1:mx,1:2)
       vlim_w(1:my)=alf*sv1_w(1:my)+(1.0_rlg-alf)*sv2_w(1:my)
       vlim_e(1:my)=alf*sv1_e(1:my)+(1.0_rlg-alf)*sv2_e(1:my)

       hlim_s(1:mx)=alf*sxe1_s(1:mx)+(1.0_rlg-alf)*sxe2_s(1:mx)
       hlim_n(1:mx)=alf*sxe1_n(1:mx)+(1.0_rlg-alf)*sxe2_n(1:mx)
       hlim_w(1:my)=alf*sxe1_w(1:my)+(1.0_rlg-alf)*sxe2_w(1:my)
       hlim_e(1:my)=alf*sxe1_e(1:my)+(1.0_rlg-alf)*sxe2_e(1:my)

     ELSE
       st1=st2

       su1_s(1:mx)=su2_s(1:mx)
       su1_n(1:mx)=su2_n(1:mx)
       su1_w(1:2,1:my)=su2_w(1:2,1:my)
       su1_e(1:2,1:my)=su2_e(1:2,1:my)

       sv1_s(1:mx,1:2)=sv2_s(1:mx,1:2)
       sv1_n(1:mx,1:2)=sv2_n(1:mx,1:2)
       sv1_w(1:my)=sv2_w(1:my)
       sv1_e(1:my)=sv2_e(1:my)

       sxe1_s(1:mx)=sxe2_s(1:mx)
       sxe1_n(1:mx)=sxe2_n(1:mx)
       sxe1_w(1:my)=sxe2_w(1:my)
       sxe1_e(1:my)=sxe2_e(1:my)

       READ(14,END=34,err=35)st2,(sxe2_n(i),sxe2_s(i),su2_n(i),sv2_n(i,1),sv2_n(i,2),su2_s(i),sv2_s(i,1),sv2_s(i,2),i=1,mx)
       READ(14)                  (sxe2_w(j),sxe2_e(j),su2_w(1,j),su2_w(2,j),sv2_w(j),su2_e(1,j),su2_e(2,j),sv2_e(j),j=1,my)

       GOTO 9
     END IF

   END IF

   DO i=2,mx-1
     as(i)=(hlim_s(i+1)+hlim_s(i-1)-2.0_rsh*hlim_s(i))/(2.0_rlg*dgp*dgp)
     an(i)=(hlim_n(i+1)+hlim_n(i-1)-2.0_rsh*hlim_n(i))/(2.0_rlg*dgp*dgp)
     bs(i)=(hlim_s(i)-hlim_s(i-1))/dgp-as(i)*(xlim(i)+xlim(i-1))
     bn(i)=(hlim_n(i)-hlim_n(i-1))/dgp-an(i)*(xlim(i)+xlim(i-1))
     cs(i)=hlim_s(i)-(as(i)*xlim(i)+bs(i))*xlim(i)
     cn(i)=hlim_n(i)-(an(i)*xlim(i)+bn(i))*xlim(i)
   END DO

   IF (limin==imin .AND. ljmin==jmin) THEN
     obc_tide_sshw(1)=obc_tide_sshw(1)-hlim_s(1)
     obc_tide_sshs(1)=obc_tide_sshs(1)-hlim_s(1)
   ENDIF

   l=2
   DO i=liminp1m2,limaxp2
     x=(i-1)*dgac
     IF(x <= xlim(2))THEN
       obc_tide_sshn(i)=obc_tide_sshn(i)-((an(2)*x+bn(2))*x+cn(2))
       obc_tide_sshs(i)=obc_tide_sshs(i)-((as(2)*x+bs(2))*x+cs(2))
     ELSEIF(x >= xlim(mx-1))THEN
       obc_tide_sshs(i)=obc_tide_sshs(i)-((as(mx-1)*x+bs(mx-1))*x+cs(mx-1))
       obc_tide_sshn(i)=obc_tide_sshn(i)-((an(mx-1)*x+bn(mx-1))*x+cn(mx-1))
     ELSE
       DO WHILE(x > xlim(l+1))
         l=l+1
       END DO
       h1=(an(l  )*x+bn(l  ))*x+cn(l  )
       h2=(an(l+1)*x+bn(l+1))*x+cn(l+1)
       obc_tide_sshn(i)=obc_tide_sshn(i)-(h1*(xlim(l+1)-x)+(x-xlim(l))*h2)/dgp
       h1=(as(l  )*x+bs(l  ))*x+cs(l  )
       h2=(as(l+1)*x+bs(l+1))*x+cs(l+1)
       obc_tide_sshs(i)=obc_tide_sshs(i)-(h1*(xlim(l+1)-x)+(x-xlim(l))*h2)/dgp
     END IF
   END DO

   IF (limax==imax .AND. ljmax==jmax) THEN
     hautd=obc_tide_sshn(imax)
   ENDIF

   DO i=2,my-1
     as(i)=(hlim_w(i+1)+hlim_w(i-1)-2.0_rsh*hlim_w(i))/(2.0_rlg*dfip*dfip)
     an(i)=(hlim_e(i+1)+hlim_e(i-1)-2.0_rsh*hlim_e(i))/(2.0_rlg*dfip*dfip)
     bs(i)=(hlim_w(i)-hlim_w(i-1))/dfip -as(i)*(ylim(i)+ylim(i-1))
     bn(i)=(hlim_e(i)-hlim_e(i-1))/dfip -an(i)*(ylim(i)+ylim(i-1))
     cs(i)=hlim_w(i)-(as(i)*ylim(i)+bs(i))*ylim(i)
     cn(i)=hlim_e(i)-(an(i)*ylim(i)+bn(i))*ylim(i)
   END DO

   l=2
   DO j=ljminp1m2,ljmaxp2
     yy=(j-1)*dfiac
     IF(yy <= ylim(2))THEN
       obc_tide_sshe(j)=obc_tide_sshe(j)-((an(2)*yy+bn(2))*yy+cn(2))
       obc_tide_sshw(j)=obc_tide_sshw(j)-((as(2)*yy+bs(2))*yy+cs(2))
     ELSEIF(yy >= ylim(my-1))THEN
       obc_tide_sshw(j)=obc_tide_sshw(j)-((as(my-1)*yy+bs(my-1))*yy+cs(my-1))
       obc_tide_sshe(j)=obc_tide_sshe(j)-((an(my-1)*yy+bn(my-1))*yy+cn(my-1))
     ELSE
       DO WHILE(yy > ylim(l+1))
       l=l+1
       END DO
       h1=(an(l  )*yy+bn(l  ))*yy+cn(l  )
       h2=(an(l+1)*yy+bn(l+1))*yy+cn(l+1)
       obc_tide_sshe(j)=obc_tide_sshe(j)-(h1*(ylim(l+1)-yy)+(yy-ylim(l))*h2)/dfip

       h1=(as(l  )*yy+bs(l  ))*yy+cs(l  )
       h2=(as(l+1)*yy+bs(l+1))*yy+cs(l+1)
       obc_tide_sshw(j)=obc_tide_sshw(j)-(h1*(ylim(l+1)-yy)+(yy-ylim(l))*h2)/dfip
     END IF
   END DO

   IF (limax==imax .AND. ljmax==jmax) THEN
     obc_tide_sshe(jmax)=0.5*(obc_tide_sshe(jmax)+hautd)
     obc_tide_sshn(imax)=obc_tide_sshe(jmax)
   ENDIF

   IF (ljmax==jmax) THEN
     WHERE(h0(liminp1m2:limaxp2,jmax)==-valmanq) obc_tide_sshn(liminp1m2:limaxp2)=valmanq
   ENDIF
   IF (ljmin==jmin) THEN
     WHERE(h0(liminp1m2:limaxp2,jmin+1)==-valmanq) obc_tide_sshs(liminp1m2:limaxp2)=valmanq
   ENDIF
   IF (limax==imax) THEN
     WHERE(h0(imax,ljminp1m2:ljmaxp2)==-valmanq) obc_tide_sshe(ljminp1m2:ljmaxp2)=valmanq
   ENDIF
   IF (limin==imin) THEN
     WHERE(h0(imin+1,ljminp1m2:ljmaxp2)==-valmanq) obc_tide_sshw(ljminp1m2:ljmaxp2)=valmanq
   ENDIF

! linear interpolation of velocities

   idh=2
 
   DO i=liminp1,limaxm1
     x=(i-0.5)*dgac
     IF(x  > xlim(idh)+.5_rsh*dgp)THEN
!DO WHILE (x >= xlim(idh)+.5_rsh*dgp.AND. idh.LE.mx) ! to make sure seq and mpi choose the same jhh when x=ylim
       DO WHILE (x > xlim(idh)+.5_rsh*dgp.AND. idh.LE.mx)
        idh=idh+1
       END DO
     END IF
     obc_tide_un(i,jmax)=  ((x-xlim(idh-1)-.5_rsh*dgp)*ulim_n(idh)+(xlim(idh)+.5_rsh*dgp-x)*ulim_n(idh-1))/dgp
     obc_tide_us(i,jmin+1)=((x-xlim(idh-1)-.5_rsh*dgp)*ulim_s(idh)+(xlim(idh)+.5_rsh*dgp-x)*ulim_s(idh-1))/dgp
   END DO

   idh=2
   DO i=liminp1,limaxm1
     x=(i-1)*dgac
     IF(x  > xlim(idh))THEN
!DO WHILE (x >= xlim(idh).AND. idh.LE.mx) ! to make sure seq and mpi choose the same jhh when x=ylim
       DO WHILE (x > xlim(idh).AND. idh.LE.mx)
        idh=idh+1
       END DO
     END IF
     obc_tide_vn(i,jmax-1:jmax)=((x-xlim(idh-1))*vlim_n(idh,1:2)+(xlim(idh)-x)*vlim_n(idh-1,1:2))/dgp
     obc_tide_vs(i,jmin:jmin+1)=((x-xlim(idh-1))*vlim_s(idh,1:2)+(xlim(idh)-x)*vlim_s(idh-1,1:2))/dgp
   END DO

   jhh=2
   DO j=ljminp1,ljmaxm1
     x=(j-0.5)*dfiac
     IF(x  > ylim(jhh)+.5_rsh*dfip)THEN
!DO WHILE (x >= ylim(jhh)+.5_rsh*dfip.AND. jhh.LE.my) ! to make sure seq and mpi choose the same jhh when x=ylim
       DO WHILE (x > ylim(jhh)+.5_rsh*dfip.AND. jhh.LE.my)
        jhh=jhh+1
       END DO
     END IF
     obc_tide_ve(imax,j)  =((x-ylim(jhh-1)-.5_rsh*dfip)*vlim_e(jhh)+(ylim(jhh)+.5_rsh*dfip-x)*vlim_e(jhh-1))/dfip
     obc_tide_vw(imin+1,j)=((x-ylim(jhh-1)-.5_rsh*dfip)*vlim_w(jhh)+(ylim(jhh)+.5_rsh*dfip-x)*vlim_w(jhh-1))/dfip
   END DO

   jhh=2
   DO j=ljminp1,ljmaxm1
     x=(j-1)*dfiac
     IF(x  > ylim(jhh))THEN
!DO WHILE (x >= ylim(jhh).AND. jhh.LE.my )
       DO WHILE (x > ylim(jhh).AND. jhh.LE.my ) ! to make sure seq and mpi choose the same jhh when x=ylim
        jhh=jhh+1
       END DO
     END IF
     obc_tide_ue(imax-1:imax,j)=((x-ylim(jhh-1))*ulim_e(1:2,jhh)+(ylim(jhh)-x)*ulim_e(1:2,jhh-1))/dfip
     obc_tide_uw(imin:imin+1,j)=((x-ylim(jhh-1))*ulim_w(1:2,jhh)+(ylim(jhh)-x)*ulim_w(1:2,jhh-1))/dfip
   END DO


   RETURN
34 WRITE(ierrorlog,*) 'ERROR: In OBC2D_RANK.F90,'
   WRITE(ierrorlog,*) 'The open boundary file (from previous rank) is not long enough.'
   WRITE(ierrorlog,*) 'Simulation stopped.'
   STOP
35 WRITE(ierrorlog,*) 'ERROR: In OBC2D_RANK.F90,'
   WRITE(ierrorlog,*) 'Problem while reading the 2nd open boundary file (from previous rank).'
   WRITE(ierrorlog,*) 'Simulation stopped.'
   STOP

  END SUBROUTINE obc2d_rank_default_surge
