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