source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90

module fxhyp_m

  IMPLICIT NONE

CONTAINS

  SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)

    ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
    ! Author: P. Le Van, from formulas by R. Sadourny

    ! Calcule les longitudes et dérivées dans la grille du GCM pour
    ! une fonction f(x) à dérivée tangente hyperbolique.

    ! 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.

    USE lmdz_dimensions, ONLY: iim
    USE lmdz_arth, ONLY: arth
    USE invert_zoom_x_m, ONLY: invert_zoom_x, nmax
    USE lmdz_physical_constants, ONLY: pi, pi_d, twopi, twopi_d, k8
    USE principal_cshift_m, ONLY: principal_cshift
    USE serre_mod, ONLY: clon, grossismx, dzoomx, taux

    REAL, INTENT(OUT) :: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
    REAL, INTENT(OUT) :: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)

    ! Local:
    REAL rlonm025(iim + 1), rlonp025(iim + 1)
    REAL dzoom, step
    REAL d_rlonv(iim)
    REAL(K8) xtild(0:2 * nmax)
    REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
    REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax)
    REAL(K8) fa, fb
    INTEGER i, is2
    REAL(K8) xmoy, fxm

    !----------------------------------------------------------------------

    PRINT *, "Call sequence information: fxhyp"

    test_iim: if (iim==1) THEN
      rlonv(1) = 0.
      rlonu(1) = pi
      rlonv(2) = rlonv(1) + twopi
      rlonu(2) = rlonu(1) + twopi

      xprimm025(:) = 1.
      xprimv(:) = 1.
      xprimu(:) = 1.
      xprimp025(:) = 1.
    else test_iim
      test_grossismx: if (grossismx == 1.) THEN
        step = twopi / iim

        xprimm025(:iim) = step
        xprimp025(:iim) = step
        xprimv(:iim) = step
        xprimu(:iim) = step

        rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim)
        rlonm025(:iim) = rlonv(:iim) - 0.25 * step
        rlonp025(:iim) = rlonv(:iim) + 0.25 * step
        rlonu(:iim) = rlonv(:iim) + 0.5 * step
      else test_grossismx
        dzoom = dzoomx * twopi_d
        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)

        ! 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

          IF (xtild(i) == 0.) fhyp(i) = 1.
          IF (xtild(i) == pi_d) fhyp(i) = -1.
        END DO

        ! 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))
        END DO

        PRINT *, "ffdx = ", ffdx
        beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
        PRINT *, "beta = ", beta

        IF (2. * beta - grossismx <= 0.) THEN
          PRINT *, 'Bad choice of grossismx, taux, dzoomx.'
          PRINT *, 'Decrease dzoomx or grossismx.'
          STOP 1
        END IF

        ! calcul de Xprimt
        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)

        ! 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

        DO i = 1, 2 * nmax - 1
          Xf(i) = Xf(i - 1) + xxpr(i) * (xtild(i) - xtild(i - 1))
        END DO

        Xf(2 * nmax) = pi_d

        CALL invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), &
                xprimm025(:iim), xuv = - 0.25_k8)
        CALL invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
                xuv = 0._k8)
        CALL invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
                xuv = 0.5_k8)
        CALL invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
                xprimp025(:iim), xuv = 0.25_k8)
      end if test_grossismx

      is2 = 0

      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
              .OR. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
        IF (clon <= 0.) THEN
          is2 = 1

          DO while (rlonm025(is2) < - pi .AND. is2 < iim)
            is2 = is2 + 1
          END DO

          IF (rlonm025(is2) < - pi) THEN
            PRINT *, 'Rlonm025 plus petit que - pi !'
            STOP 1
          end if
        ELSE
          is2 = iim

          DO while (rlonm025(is2) > pi .AND. is2 > 1)
            is2 = is2 - 1
          END DO

          IF (rlonm025(is2) > pi) THEN
            PRINT *, 'Rlonm025 plus grand que pi !'
            STOP 1
          end if
        END IF
      END IF

      CALL principal_cshift(is2, rlonm025, xprimm025)
      CALL principal_cshift(is2, rlonv, xprimv)
      CALL principal_cshift(is2, rlonu, xprimu)
      CALL principal_cshift(is2, rlonp025, xprimp025)

      forall (i = 1:iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
      PRINT *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, &
              "degrees"
      PRINT *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, &
              "degrees"

      ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
      DO i = 1, iim + 1
        IF (rlonp025(i) < rlonv(i)) THEN
          PRINT *, 'rlonp025(', i, ') = ', rlonp025(i)
          PRINT *, "< rlonv(", i, ") = ", rlonv(i)
          STOP 1
        END IF

        IF (rlonv(i) < rlonm025(i)) THEN
          PRINT *, 'rlonv(', i, ') = ', rlonv(i)
          PRINT *, "< rlonm025(", i, ") = ", rlonm025(i)
          STOP 1
        END IF

        IF (rlonp025(i) > rlonu(i)) THEN
          PRINT *, 'rlonp025(', i, ') = ', rlonp025(i)
          PRINT *, "> rlonu(", i, ") = ", rlonu(i)
          STOP 1
        END IF
      END DO
    end if test_iim

  END SUBROUTINE fxhyp

END MODULE fxhyp_m