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, pi_d, twopi_d, arth
    use principal_cshift_m, only: principal_cshift
    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
    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 xzoom, fa, fb
    INTEGER i, is2
    DOUBLE PRECISION xmoy, fxm
    DOUBLE PRECISION decalx

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

    print *, "Call sequence information: fxhyp"

    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

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

    xzoom = clon * pi_d / 180d0
    call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025(:iim), &
         xprimm025(:iim), xuv = - 0.25d0)
    call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv(:iim), &
         xprimv(:iim), xuv = 0d0)
    call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu(:iim), &
         xprimu(:iim), xuv = 0.5d0)
    call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025(:iim), &
         xprimp025(:iim), xuv = 0.25d0)

    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"

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