trunk/dyn3d/fxhyp.f

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 dimens_m, ONLY: iim
    use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
    use nr_util, only: pi_d, twopi_d, arth
    use serre, 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
    DOUBLE PRECISION xlon(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 xzoom, fa, fb
    INTEGER i
    DOUBLE PRECISION xmoy, fxm
    DOUBLE PRECISION decalx

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

    print *, "Call sequence information: fxhyp"

    xzoom = clon * pi_d / 180d0

    IF (grossismx == 1.) THEN
       decalx = 1d0
    else
       decalx = 0.75d0
    END IF

    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 *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
       print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
       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 fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
         xprimm025, xuv = - 0.25d0)
    call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
         xuv = 0d0)
    call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
         xuv = 0.5d0)
    call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
         xprimp025, xuv = 0.25d0)

    print *

    forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
    print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "°"
    print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "°"

    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 SUBROUTINE fxhyp

end module fxhyp_m
1	module fxhyp_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8
9	! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10	! Author: P. Le Van, from formulas by R. Sadourny
11
12	! Calcule les longitudes et dérivées dans la grille du GCM pour
13	! une fonction f(x) à dérivée tangente hyperbolique.
14
15	! Il vaut mieux avoir : grossismx \times dzoom < pi
16
17	! Le premier point scalaire pour une grille regulière (grossismx =
18	! 1., taux=0., clon=0.) est à - 180 degrés.
19
20	USE dimens_m, ONLY: iim
21	use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22	use nr_util, only: pi_d, twopi_d, arth
23	use serre, only: clon, grossismx, dzoomx, taux
24
25	REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
26	real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
27
28	! Local:
29
30	real rlonm025(iim + 1), rlonp025(iim + 1)
31	REAL dzoom
32	DOUBLE PRECISION xlon(iim)
33	DOUBLE PRECISION xtild(0:2 * nmax)
34	DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35	DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36	DOUBLE PRECISION xzoom, fa, fb
37	INTEGER i
38	DOUBLE PRECISION xmoy, fxm
39	DOUBLE PRECISION decalx
40
41	!----------------------------------------------------------------------
42
43	print *, "Call sequence information: fxhyp"
44
45	xzoom = clon * pi_d / 180d0
46
47	IF (grossismx == 1.) THEN
48	decalx = 1d0
49	else
50	decalx = 0.75d0
51	END IF
52
53	dzoom = dzoomx * twopi_d
54	xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
55
56	! Compute fhyp:
57	DO i = nmax, 2 * nmax
58	fa = taux * (dzoom / 2. - xtild(i))
59	fb = xtild(i) * (pi_d - xtild(i))
60
61	IF (200. * fb < - fa) THEN
62	fhyp(i) = - 1.
63	ELSE IF (200. * fb < fa) THEN
64	fhyp(i) = 1.
65	ELSE
66	IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
67	IF (200. * fb + fa < 1e-10) THEN
68	fhyp(i) = - 1.
69	ELSE IF (200. * fb - fa < 1e-10) THEN
70	fhyp(i) = 1.
71	END IF
72	ELSE
73	fhyp(i) = TANH(fa / fb)
74	END IF
75	END IF
76
77	IF (xtild(i) == 0.) fhyp(i) = 1.
78	IF (xtild(i) == pi_d) fhyp(i) = -1.
79	END DO
80
81	! Calcul de beta
82
83	ffdx = 0.
84
85	DO i = nmax + 1, 2 * nmax
86	xmoy = 0.5 * (xtild(i-1) + xtild(i))
87	fa = taux * (dzoom / 2. - xmoy)
88	fb = xmoy * (pi_d - xmoy)
89
90	IF (200. * fb < - fa) THEN
91	fxm = - 1.
92	ELSE IF (200. * fb < fa) THEN
93	fxm = 1.
94	ELSE
95	IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
96	IF (200. * fb + fa < 1e-10) THEN
97	fxm = - 1.
98	ELSE IF (200. * fb - fa < 1e-10) THEN
99	fxm = 1.
100	END IF
101	ELSE
102	fxm = TANH(fa / fb)
103	END IF
104	END IF
105
106	IF (xmoy == 0.) fxm = 1.
107	IF (xmoy == pi_d) fxm = -1.
108
109	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
110	END DO
111
112	print *, "ffdx = ", ffdx
113	beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
114	print *, "beta = ", beta
115
116	IF (2. * beta - grossismx <= 0.) THEN
117	print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
118	print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
119	STOP 1
120	END IF
121
122	! calcul de Xprimt
123	Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
124	xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
125
126	! Calcul de Xf
127
128	DO i = nmax + 1, 2 * nmax
129	xmoy = 0.5 * (xtild(i-1) + xtild(i))
130	fa = taux * (dzoom / 2. - xmoy)
131	fb = xmoy * (pi_d - xmoy)
132
133	IF (200. * fb < - fa) THEN
134	fxm = - 1.
135	ELSE IF (200. * fb < fa) THEN
136	fxm = 1.
137	ELSE
138	fxm = TANH(fa / fb)
139	END IF
140
141	IF (xmoy == 0.) fxm = 1.
142	IF (xmoy == pi_d) fxm = -1.
143	xxpr(i) = beta + (grossismx - beta) * fxm
144	END DO
145
146	xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
147
148	Xf(0) = - pi_d
149
150	DO i=1, 2 * nmax - 1
151	Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
152	END DO
153
154	Xf(2 * nmax) = pi_d
155
156	call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
157	xprimm025, xuv = - 0.25d0)
158	call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
159	xuv = 0d0)
160	call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
161	xuv = 0.5d0)
162	call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
163	xprimp025, xuv = 0.25d0)
164
165	print *
166
167	forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
168	print , "Minimum longitude step:", MINval(xlon) 180. / pi_d, "°"
169	print , "Maximum longitude step:", MAXval(xlon) 180. / pi_d, "°"
170
171	DO i = 1, iim + 1
172	IF (rlonp025(i) < rlonv(i)) THEN
173	print *, 'rlonp025(', i, ') = ', rlonp025(i)
174	print *, "< rlonv(", i, ") = ", rlonv(i)
175	STOP 1
176	END IF
177
178	IF (rlonv(i) < rlonm025(i)) THEN
179	print *, 'rlonv(', i, ') = ', rlonv(i)
180	print *, "< rlonm025(", i, ") = ", rlonm025(i)
181	STOP 1
182	END IF
183
184	IF (rlonp025(i) > rlonu(i)) THEN
185	print *, 'rlonp025(', i, ') = ', rlonp025(i)
186	print *, "> rlonu(", i, ") = ", rlonu(i)
187	STOP 1
188	END IF
189	END DO
190
191	END SUBROUTINE fxhyp
192
193	end module fxhyp_m