--- trunk/dyn3d/fxhyp.f	2015/03/20 16:31:06	132
+++ trunk/Sources/dyn3d/fxhyp.f	2015/06/16 15:23:29	145
@@ -15,13 +15,14 @@
     ! Il vaut mieux avoir : grossismx \times dzoom < pi
 
     ! Le premier point scalaire pour une grille regulière (grossismx =
-    ! 1., taux=0., clon=0.) est à - 180 degrés.
+    ! 1., taux = 0., clon = 0.) est à - 180 degrés.
 
     USE dimens_m, ONLY: iim
+    use dynetat0_m, only: clon, grossismx, dzoomx, taux
     use invert_zoom_x_m, only: invert_zoom_x, nmax
     use nr_util, only: pi, pi_d, twopi, twopi_d, arth
     use principal_cshift_m, only: principal_cshift
-    use serre, only: clon, grossismx, dzoomx, taux
+    use tanh_cautious_m, only: tanh_cautious
 
     REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
     real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
@@ -32,10 +33,9 @@
     real d_rlonv(iim)
     DOUBLE PRECISION xtild(0:2 * nmax)
     DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
-    DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
-    DOUBLE PRECISION fa, fb
+    DOUBLE PRECISION Xf(0:2 * nmax)
     INTEGER i, is2
-    DOUBLE PRECISION xmoy, fxm
+    DOUBLE PRECISION, dimension(nmax + 1:2 * nmax):: xxpr, xmoy, fxm
 
     !----------------------------------------------------------------------
 
@@ -56,65 +56,27 @@
     else test_grossismx
        dzoom = dzoomx * twopi_d
        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
+       forall (i = nmax + 1:2 * nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
 
        ! Compute fhyp:
-       DO i = nmax, 2 * nmax
-          fa = taux * (dzoom / 2. - xtild(i))
-          fb = xtild(i) * (pi_d - xtild(i))
-
-          IF (200. * fb < - fa) THEN
-             fhyp(i) = - 1.
-          ELSE IF (200. * fb < fa) THEN
-             fhyp(i) = 1.
-          ELSE
-             IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
-                IF (200. * fb + fa < 1e-10) THEN
-                   fhyp(i) = - 1.
-                ELSE IF (200. * fb - fa < 1e-10) THEN
-                   fhyp(i) = 1.
-                END IF
-             ELSE
-                fhyp(i) = TANH(fa / fb)
-             END IF
-          END IF
+       fhyp(nmax + 1:2 * nmax - 1) = tanh_cautious(taux * (dzoom / 2. &
+            - xtild(nmax + 1:2 * nmax - 1)), xtild(nmax + 1:2 * nmax - 1) &
+            * (pi_d - xtild(nmax + 1:2 * nmax - 1)))
+       fhyp(nmax) = 1d0
+       fhyp(2 * nmax) = -1d0
 
-          IF (xtild(i) == 0.) fhyp(i) = 1.
-          IF (xtild(i) == pi_d) fhyp(i) = -1.
-       END DO
+       fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
 
        ! Calcul de beta 
 
        ffdx = 0.
 
        DO i = nmax + 1, 2 * nmax
-          xmoy = 0.5 * (xtild(i-1) + xtild(i))
-          fa = taux * (dzoom / 2. - xmoy)
-          fb = xmoy * (pi_d - xmoy)
-
-          IF (200. * fb < - fa) THEN
-             fxm = - 1.
-          ELSE IF (200. * fb < fa) THEN
-             fxm = 1.
-          ELSE
-             IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
-                IF (200. * fb + fa < 1e-10) THEN
-                   fxm = - 1.
-                ELSE IF (200. * fb - fa < 1e-10) THEN
-                   fxm = 1.
-                END IF
-             ELSE
-                fxm = TANH(fa / fb)
-             END IF
-          END IF
-
-          IF (xmoy == 0.) fxm = 1.
-          IF (xmoy == pi_d) fxm = -1.
-
-          ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
+          ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
        END DO
 
        print *, "ffdx = ", ffdx
-       beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
+       beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
        print *, "beta = ", beta
 
        IF (2. * beta - grossismx <= 0.) THEN
@@ -129,33 +91,15 @@
 
        ! Calcul de Xf
 
-       DO i = nmax + 1, 2 * nmax
-          xmoy = 0.5 * (xtild(i-1) + xtild(i))
-          fa = taux * (dzoom / 2. - xmoy)
-          fb = xmoy * (pi_d - xmoy)
-
-          IF (200. * fb < - fa) THEN
-             fxm = - 1.
-          ELSE IF (200. * fb < fa) THEN
-             fxm = 1.
-          ELSE
-             fxm = TANH(fa / fb)
-          END IF
-
-          IF (xmoy == 0.) fxm = 1.
-          IF (xmoy == pi_d) fxm = -1.
-          xxpr(i) = beta + (grossismx - beta) * fxm
-       END DO
-
-       xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
-
-       Xf(0) = - pi_d
+       xxpr = beta + (grossismx - beta) * fxm
+       Xf(nmax) = 0d0
 
-       DO i=1, 2 * nmax - 1
+       DO i = nmax + 1, 2 * nmax - 1
           Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
        END DO
 
        Xf(2 * nmax) = pi_d
+       xf(:nmax - 1) = - xf(2 * nmax:nmax + 1:- 1)
 
        call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
             xuv = - 0.25d0)