SUBROUTINE fxyhyp( clat,  rlatu,yprimu ,  rlatv,yprimv   ,
     *                         rlatu1 ,yprimu1, rlatu2,yprimu2   )
c
       IMPLICIT NONE
c   -----------------------
c     Auteur  :  P. Le Van 
c   -----------------------
c
c .. Objet : Calcul d'une fonction f(y)  ayant pour derivee une tangente
c        hyperbolique , donc des ordonnees y de la grille horizontale ,
c       ainsi que ses derivees yprim , utlisees dans  Inigeom  .
c
c
c ....  rlatu   =  fxyhyp (   j      )     avec 1  <= j  <= jjm + 1 
c ....  rlatv   =  fxyhyp ( j + 0.5  )     avec 1  <= j  <= jjm
c ....  rlatu1  =  fxyhyp ( j + 0.25 )     avec 1  <= j  <= jjm 
c ....  rlatu2  =  fxyhyp ( j + 0.75 )     avec 1  <= j  <= jjm
c

       INTEGER nmax
       PARAMETER ( nmax = 100 000 )
c
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
c
c    .................   Arguments     ................................
c
c   clat est (en degres) la latitude du centre du zoom : arg. d'entree
c    Les autres arguments     sont                des argum. de sortie
c
       REAL clat
       REAL rlatu(jjp1),rlatv(jjm),yprimu(jjp1),yprimv(jjm)
       REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm)
c
c    .................................................................
c
c             ...... .  Variables  locales   .........   
c
       REAL yprim( 0:nmax ),yi( 0:nmax ),yf( 0:nmax )
       REAL yyy(0:nmax)
       REAL delta,den,eps,phij,pijjm,pis2,pis2pyo,pism,psi,psi0
       REAL rappis2,unpdel,y,yint,yj,yo,yo1,yppp,ypr 
       REAL yprimin,pisnmax,ydeg
       INTEGER i,j,it,iter
c
c
c -------------------------------------------------------------------------
   
c   On veut calculer  y(Y) et y'(Y) , avec  Y ,latitude de travail =f(j),
c       et   y , latitude vraie .
c
c  Pour cela ,on se donne pour  Y'(y) une fonction hyperbolique calculee
c  sur  ( 1+ nmax ) points ( plus nmax est grand, plus la precision pour 
c  y(Y)  sera grande ) , on integre Y'(y) sur ces ( 1+nmax) pts  pour
c  obtenir  Y(y) apres normalisation  , puis on  aura y(Y)  en resolvant
c   l'equation  Y(yj) = Yj  par Newton_Raphson ,  avec Yj = f(j) =
c   pi/2 + pi/jjm ( 1.- j)  pour les scalaires  ( j =1, jjm+1 )
c   On en tire donc  yj  qui correspond  a  la ligne  j   .
c   On a ensuite pour les derivees :  y'(Y) = 1./ Y'(y) * ( pi/ jjm )
c
c   On calculera   y(Y) et y'(Y)  successivement   aux  latitudes  de  U 
c     ,  V  ,  U + 0.25 , U + 0.75   .
c
c   Notations ; 
c  -------------
c     yi    represente les nmax+1 latitudes de depart ( -pi/2 a pi/2 )
c  yyy(nmax)  repres.  l'integrale de Y'(y) de - pi/2 a pi/2
c    yprim  represente  Y'(y)
c     yf    represente  Y (y)
c
c     yo    represente la latitude du centre du zoom
c    delta  represente la demi largeur du zoom
c
c
c  N.B :  La fonction hyperbolique Y'(y) est  definie comme suit :
c ------
c
c    Y' = TANH ( ( psio - psi )/psi*(1- psi) )  avec  :
c
c          psi =  ( y - yo ) / ( pi/2 - yo )  si     yo < y < pi/2
c          psi =  (yo -  y )/  ( pi/2 + yo )  si - pi/2 < y < yo
c
c
c  ....  yo    represente la latitude du centre du zoom         .....
c
c  ....  psi0  et delta sont des parametres variant entre 0 et 1  ....
c----------------------------------------------------------------------
c
cc     psi0 = 0.3 et  delta = 0.5  sont les valeurs qui conviennent 
c      pour  *** iim = 96  et jjm = 72 ***  , d'apres des tests sur
c      les valeurs du Laplacien itere ( 2 fois) , analytiques et cal-
c      culees.   Pour d'autres valeurs de  iim et jjm , on peut tester
c      les nouvelles valeurs de psi0 et delta en lancant  testfxyhyp .
c
c ----------------------------------------------------------------------
c    Fxyhp  est appele par inigeom  a la place de  fxynew (  fonction
c    a derivee sinusoidale )  si la variable logique fxyhypb , lue par
c    defrun_new  est  . TRUE.    .
c ----------------------------------------------------------------------
c

       pi     =  2.*ASIN( 1. )
       pis2   =  pi/2.
       yo     =  clat  * pi/180.
       psi0   =  0.3
       delta  =  0.5
       unpdel =  1. + delta
       eps    =  .1e-7
c
       PRINT *,' ----------------------------------------------------'
       PRINT 1,    psi0,delta,clat
       PRINT *,' ----------------------------------------------------'

       pisnmax =  pi/FLOAT(nmax)
       yi(0)   = - pis2

       DO i = 1, nmax
        yi(i)= - pis2 + FLOAT(i)*pisnmax
       ENDDO

       yyy(0)=  0.

c
c     ***************************************************************
c     ************      Cas   ou   yo  <  0.     ********************
c     ***************************************************************
c
       IF( yo.LT. 0. )   THEN
c
       pis2pyo = - pis2 + yo
       rappis2 = ( -pis2-yo)/pis2pyo
       rappis2 = rappis2 * rappis2
c
        DO 5 i = 1, nmax

        y = 0.5 * ( yi(i)+ yi(i-1) )
        IF( y.LT.- pis2.OR.y.GT.pis2 )  STOP 0

        IF(y.GT.yo)  THEN
        psi = (y-yo)/(pis2-yo)
        ypr =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
         
        ELSE
        psi = (yo-y)/(pis2+yo)
        ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) )
        ENDIF

        yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1))

 5    CONTINUE

       den = unpdel + yyy(nmax)/pi

       yyy(0)=   0.
       yf(0) = - pis2

       yprim(nmax) = -1.
       yprim(0)    = 1.- rappis2 - rappis2

       DO 12 i =  0 ,nmax
        y = yi(i)
       IF(y.LT.-pis2.OR.y.GT.pis2)  STOP 1

       IF(y.LT.yo)  THEN
           psi = (yo -y)/(pis2+yo)
         IF(i.NE.0.AND.i.NE.nmax)  THEN 
         yprim(i) = 1.- rappis2 + rappis2
     *                 * TANH( (psi0 -psi)/(psi*(1.-psi)) )
         ENDIF
       ELSE
        psi = (y-yo)/(pis2-yo)
          IF(i.NE.0.AND.i.NE.nmax) THEN
           yprim(i) =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
          ENDIF
       ENDIF

c
c   ..................................................... 
c   .......   Normalisation de  Y'(y) et de  Y(y)   .....
c   .....................................................
c
        yprim(i)= ( yprim(i)+ unpdel )/ den
        yf(i)= (yyy(i)+  unpdel*( yi(i)+ pis2))/ den  - pis2
c
 12    CONTINUE
c
      ELSE
c
c     **************************************************************
c     **********          Cas  ou  yo  >=  0 .        **************
c     **************************************************************
c
       pis2pyo = pis2 + yo
       rappis2 = ( pis2 - yo )/pis2pyo
       rappis2 = rappis2 * rappis2
c

       DO 13 i = 1, nmax
        y= 0.5 * ( yi(i)+ yi(i-1) )
        IF(y.LT.- pis2.OR.y.GT.pis2)   STOP 2

       IF(y.LT.yo)  THEN
        psi = (yo -y)/(pis2+yo)
        ypr =  TANH( (psi0 - psi)/(psi*(1.-psi)) )

       ELSE
        psi = (y- yo)/(pis2-yo)
        ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) )

       ENDIF

        yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1))

 13    CONTINUE

       den = unpdel + yyy(nmax)/pi

       yyy(0)= 0.
       yf(0) = -pis2

       yprim(0)   = -1.
       yprim(nmax)= 1.- rappis2 - rappis2

       DO 15 i = 0,nmax
        y = yi(i)
       IF(y.LT.-pis2.OR.y.GT.pis2)   STOP 3

       IF(y.LT.yo)  THEN
        psi = (yo -y)/(pis2+yo)
        IF(i.NE.0.AND.i.NE.nmax) 
     *   yprim(i) =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
       ELSE
        psi = (y- yo)/(pis2-yo)
        IF(i.NE.0.AND.i.NE.nmax)
     *   yprim(i) = 1.- rappis2 + rappis2 
     *               * TANH( (psi0 -psi)/(psi*(1.-psi)) )                              
       ENDIF
c   ..................................................... 
c
c   .......   Normalisation de  Y'(y) et de  Y(y)   .....
c   .....................................................
c
        yprim(i)= ( yprim(i)+ unpdel )/ den
        yf(i)= (yyy(i)+  unpdel*( yi(i)+ pis2))/ den  - pis2

 15    CONTINUE

       ENDIF
c
c     ***************************************************************

      yf(0)   = - pis2
      yf(nmax)=   pis2
       
      pijjm = pi/FLOAT(jjm)

c   ...............................................................
c
c     Calculs  de y(Y) et y'(Y)   pour les latitudes  de  U   .....
c
      PRINT *,' ***************************** '
      PRINT *,' ***  Calcul de rlatu    ***** '
      PRINT *,' ***************************** '
c
c   ...............................................................
c
      PRINT *,' *** j   rlatu  yprimu **** '
c
      DO 1500 j=1, jjm +1

      phij = pis2 + pijjm* ( 1.-Float(j) )
      yo1  = 0.
      yj   = yo1
c
      DO 500 iter = 1,300

      DO 250 it = nmax,0,-1
      IF( yj.GE.yi(it))  GO TO 350
250   CONTINUE

      it= 0
      yj = yi(it)

350   CONTINUE
       IF(it.EQ.nmax)  THEN
         yint    = yf(nmax)
         yprimin = yprim(nmax)
      ELSE
c  .................................................................
c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yj)  
c      .....           et   Y'(yj)                             .....
c  .................................................................
      yint = (yf(it+1)-yf(it))/(yi(it+1)-yi(it))* ( yj-yi(it) )     +
     +                           yf(it)
      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
     +                           yprim(it)
      ENDIF
      yj = yj - (yint-phij)/yprimin

      IF( ABS(yj-yo1).LE.eps) GO TO 550
      yo1 = yj
c
500   CONTINUE
      PRINT *,' ***   PAS DE SOLUTION  ****  ',j,phij
        STOP 4
550   CONTINUE
       ydeg     = yj * 180./pi
       rlatu(j) = yj
       IF(j.EQ.1) rlatu(j) = pis2
c
c    .....  calcul de  yprimu .....
c

       IF(j.EQ.1.OR.j.EQ.jjm +1) GO TO 1500
       IF(yj.LT.yo)  THEN
        psi = (yo -yj)/(pis2+yo)
       ELSE
        psi = (yj-yo)/(pis2-yo)
       ENDIF

      DO  it = nmax,0,-1
      IF( yj.GE.yi(it))  GO TO 555
      ENDDO
555    CONTINUE
       IF(it.EQ.nmax) THEN
       PRINT *,' IT = nmax  Pbs ..',it,j
       STOP 5
       ENDIF
       yppp=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
     *               yprim(it)
       yprimu(j) = pijjm/yppp
c
       PRINT 3000,j,ydeg,yprimu(j)
1500   CONTINUE
       rlatu(1)     =   pis2
       rlatu(jjm+1) = - pis2
c
c     .............................................
c
c     .....   Calculs  pour  rlatv et yprimv  .....
c     .............................................
c

      PRINT *,' ***************************** '
      PRINT *,' ****  Calcul  de rlatv  ***** '
      PRINT *,' ***************************** '
      PRINT *,' *** j   rlatv  yprimv ****  '

      DO 1600 j=1, jjm 

      phij  = pis2+ pijjm*(1.-(FLOAT(j)+0.5))
      yo1   = 0.
      yj    = yo1

      DO 1550 iter =1,300

      DO 1520 it = nmax,0,-1
      IF( yj.GE.yi(it))  GO TO 1530
1520   CONTINUE
c
      it  = 0
      yj  = yi(it)
1530   CONTINUE
       IF(it.EQ.nmax)  THEN
       yint    = yf(nmax)
       yprimin = yprim(nmax)
      ELSE
      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj-yi(it) )  +
     +                          yf(it)
      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
     +                          yprim(it)
      ENDIF
      yj = yj - (yint-phij)/yprimin

      IF( ABS(yj-yo1).LE.eps ) GO TO 1560
      yo1 = yj
c
1550   CONTINUE
      PRINT *,'  ***  PAS DE SOLUTION  ***  ',j,phij
1560   CONTINUE
      ydeg    = yj * 180./pi
      rlatv(j)= yj
c
c    .....  calcul de  yprimv .....
c
       IF(yj.LT.-pis2.OR.yj.GT.pis2)   STOP 6

       IF(yj.LT.yo)  THEN
        psi = (yo -yj)/(pis2+yo)
       ELSE
        psi = (yj-yo)/(pis2-yo)
       ENDIF
      DO  it = nmax,0,-1
      IF( yj.GE.yi(it))   GO TO 1570
      ENDDO
1570    CONTINUE
       IF(it.EQ.nmax) THEN
       PRINT *,' IT = nmax  Pbs ..',it,j
       STOP 7
       ENDIF
       yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
     *                   yprim(it)
       yprimv(j)= pijjm/ yppp
c
       PRINT 3000,j,ydeg,yprimv(j)

1600  CONTINUE

c
c    ..................................................
c
      PRINT *,' ************************************ '
      PRINT *,' ****  Calcul pour f( j + 0.25 ) **** '
      PRINT *,' ************************************ '
      PRINT *,' *** j  rlatu1  yprimu1 **** '
c
c    ..................................................
c
      DO 1800 j=1, jjm 

      phij = pis2+ pijjm*( 1.- (FLOAT(j) + 0.25) )
      yo1  = 0.
      yj   = yo1
C

      DO 1620 iter = 1, 500

      DO 1610 it = nmax,0,-1
      IF( yj.GE.yi(it))  GO TO 1615
1610   CONTINUE
c
      it = 0
      yj = yi(it)
1615   CONTINUE
      IF(it.EQ.nmax)  THEN
       yint    = yf(nmax)
       yprimin = yprim(nmax)
      ELSE
      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it))*(yj-yi(it))+
     * yf(it)
      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it))+
     * yprim(it)
      ENDIF
      yj = yj - (yint-phij)/yprimin
c      IF( (- pis2-yj)*(yj- pis2 ).LT.0.)PRINT *,' Newton hors limites'

      IF( ABS(yj-yo1).LE.eps ) GO TO 1625
      yo1 = yj
c
1620   CONTINUE
      PRINT *,' ***  PAS DE SOLUTION  ****  ',j,phij
1625   CONTINUE
      ydeg     = yj * 180./pi
      rlatu1(j)= yj
c
c    .....  calcul de  yprimu1 .....
c
       IF( yj.LT.-pis2.OR.yj.GT.pis2 )  STOP 8
c
       IF(yj.LT.yo)  THEN
        psi = (yo -yj)/(pis2+yo)
       ELSE
        psi = (yj-yo)/(pis2-yo)
       ENDIF
      DO  it = nmax,0,-1
      IF( yj.GE.yi(it))   GO TO 1630
      ENDDO
1630    CONTINUE
       IF(it.EQ.nmax) THEN
       PRINT *,' IT = nmax  Pbs ..',it,j
       STOP 9
       ENDIF
       yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
     *          yprim(it)
       yprimu1(j) = pijjm/yppp
c
       PRINT 3000,j,ydeg,yprimu1(j)

1800  CONTINUE
    
c    ..............................................
c
      PRINT *,' ********************************** '
      PRINT *,' **** Calcul pour f( j + 3/4  ) *** '
      PRINT *,' ********************************** '
      PRINT *,' *** j  rlatu2  yprimu2 ***  '
c    ..............................................
c
      DO 1900 j=1, jjm 

      phij = pis2+ pijjm*(1.-(FLOAT(j)+0.75))
      yo1  = 0.
      yj   = yo1

      DO 1700 iter =1,500

      DO 1680 it = nmax,0,-1
      IF( yj.GE.yi(it))   GO TO 1690
1680   CONTINUE
      it  = 0
      yj  = yi(it)
1690   CONTINUE
      IF(it.EQ.nmax)  THEN
       yint     = yf(nmax)
       yprimin = yprim(nmax)
      ELSE
      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj -yi(it) ) +
     +                             yf(it)
      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
     +                             yprim(it)
      ENDIF
      yj = yj - (yint-phij)/yprimin

      IF( ABS(yj-yo1).LE.eps ) GO TO 1705
      yo1 = yj
c
1700   CONTINUE
      PRINT *,' PAS DE SOLUTION ',j,phij
        STOP 10
1705   CONTINUE
      ydeg      = yj * 180./pi
      rlatu2(j) = yj
c
c    .....  calcul de  yprimu2 .....
c
       IF(yj.LT.- pis2.OR.yj.GT.pis2)  STOP 11

       IF(yj.LT.yo)  THEN
        psi = (yo -yj)/(pis2+yo)
       ELSE
        psi = (yj-yo)/(pis2-yo)
       ENDIF
c
       DO  it = nmax,0,-1
       IF( yj.GE.yi(it))  GO TO 1710
       ENDDO
1710    CONTINUE
       IF(it.EQ.nmax) THEN
       PRINT *,' IT = nmax  Pbs ..',it,j
        STOP 12
       ENDIF
      yppp =(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj - yi(it) ) +
     +                       yprim(it)
       yprimu2(j) = pijjm/yppp
c
       PRINT 3000,j,ydeg,yprimu2(j)

1900  CONTINUE
c
       PRINT *,'  *** TEST DE COHERENCE DANS  FXYP **** '
c
       DO j = 1, jjm
c
       IF(rlatu1(j).LE.rlatu2(j))   THEN
         PRINT *,' Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
         STOP 13
       ENDIF
c
       IF(rlatu2(j).LE.rlatu(j+1))  THEN
        PRINT *,' Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
        STOP 14
       ENDIF
c
       IF(rlatu(j).LE.rlatu1(j))    THEN
        PRINT *,' rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
        STOP 15
       ENDIF
c
       IF(rlatv(j).LE.rlatu2(j))    THEN
        PRINT *,' rlatv > rlatu2 ',rlatv(j),rlatu2(j),j
        STOP 16
       ENDIF
c
       IF(rlatv(j).ge.rlatu1(j))    THEN
        PRINT *,' rlatv < rlatu1 ',rlatv(j),rlatu1(j),j
        STOP 17
       ENDIF
c
       IF(rlatv(j).ge.rlatu(j))     THEN
        PRINT *,' rlatv > rlatu ',rlatv(j),rlatu(j),j
        STOP 18
       ENDIF
c
       ENDDO
c
       PRINT *,' ***  Les latitudes yprimu,yprimv,yprimu1,yprimu2 ',
     * ' sont coherentes entre elles !  *** ' 

c
 1     FORMAT( 3x,'Fxyhyp.  psi0  delta clat ',3f7.2)
3000   FORMAT(2x,i4,2x,f8.2,f8.3)
       RETURN
       END