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

    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)
    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
    INTEGER i, is2
    DOUBLE PRECISION xmoy, fxm

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

    print *, "Call sequence information: fxhyp"

    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, 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.25d0)
       call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
            xuv = 0d0)
       call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
            xuv = 0.5d0)
       call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
            xuv = 0.25d0)
    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 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 invert_zoom_x_m, only: invert_zoom_x, nmax
22	use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23	use principal_cshift_m, only: principal_cshift
24	use serre, only: clon, grossismx, dzoomx, taux
25
26	REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27	real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
28
29	! Local:
30	real rlonm025(iim + 1), rlonp025(iim + 1)
31	REAL dzoom, step
32	real d_rlonv(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 fa, fb
37	INTEGER i, is2
38	DOUBLE PRECISION xmoy, fxm
39
40	!----------------------------------------------------------------------
41
42	print *, "Call sequence information: fxhyp"
43
44	test_grossismx: if (grossismx == 1.) then
45	step = twopi / iim
46
47	xprimm025(:iim) = step
48	xprimp025(:iim) = step
49	xprimv(:iim) = step
50	xprimu(:iim) = step
51
52	rlonv(:iim) = arth(- pi + clon, step, iim)
53	rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54	rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55	rlonu(:iim) = rlonv(:iim) + 0.5 * step
56	else test_grossismx
57	dzoom = dzoomx * twopi_d
58	xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59
60	! Compute fhyp:
61	DO i = nmax, 2 * nmax
62	fa = taux * (dzoom / 2. - xtild(i))
63	fb = xtild(i) * (pi_d - xtild(i))
64
65	IF (200. * fb < - fa) THEN
66	fhyp(i) = - 1.
67	ELSE IF (200. * fb < fa) THEN
68	fhyp(i) = 1.
69	ELSE
70	IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71	IF (200. * fb + fa < 1e-10) THEN
72	fhyp(i) = - 1.
73	ELSE IF (200. * fb - fa < 1e-10) THEN
74	fhyp(i) = 1.
75	END IF
76	ELSE
77	fhyp(i) = TANH(fa / fb)
78	END IF
79	END IF
80
81	IF (xtild(i) == 0.) fhyp(i) = 1.
82	IF (xtild(i) == pi_d) fhyp(i) = -1.
83	END DO
84
85	! Calcul de beta
86
87	ffdx = 0.
88
89	DO i = nmax + 1, 2 * nmax
90	xmoy = 0.5 * (xtild(i-1) + xtild(i))
91	fa = taux * (dzoom / 2. - xmoy)
92	fb = xmoy * (pi_d - xmoy)
93
94	IF (200. * fb < - fa) THEN
95	fxm = - 1.
96	ELSE IF (200. * fb < fa) THEN
97	fxm = 1.
98	ELSE
99	IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100	IF (200. * fb + fa < 1e-10) THEN
101	fxm = - 1.
102	ELSE IF (200. * fb - fa < 1e-10) THEN
103	fxm = 1.
104	END IF
105	ELSE
106	fxm = TANH(fa / fb)
107	END IF
108	END IF
109
110	IF (xmoy == 0.) fxm = 1.
111	IF (xmoy == pi_d) fxm = -1.
112
113	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114	END DO
115
116	print *, "ffdx = ", ffdx
117	beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118	print *, "beta = ", beta
119
120	IF (2. * beta - grossismx <= 0.) THEN
121	print *, 'Bad choice of grossismx, taux, dzoomx.'
122	print *, 'Decrease dzoomx or grossismx.'
123	STOP 1
124	END IF
125
126	! calcul de Xprimt
127	Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128	xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129
130	! Calcul de Xf
131
132	DO i = nmax + 1, 2 * nmax
133	xmoy = 0.5 * (xtild(i-1) + xtild(i))
134	fa = taux * (dzoom / 2. - xmoy)
135	fb = xmoy * (pi_d - xmoy)
136
137	IF (200. * fb < - fa) THEN
138	fxm = - 1.
139	ELSE IF (200. * fb < fa) THEN
140	fxm = 1.
141	ELSE
142	fxm = TANH(fa / fb)
143	END IF
144
145	IF (xmoy == 0.) fxm = 1.
146	IF (xmoy == pi_d) fxm = -1.
147	xxpr(i) = beta + (grossismx - beta) * fxm
148	END DO
149
150	xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151
152	Xf(0) = - pi_d
153
154	DO i=1, 2 * nmax - 1
155	Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156	END DO
157
158	Xf(2 * nmax) = pi_d
159
160	call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161	xuv = - 0.25d0)
162	call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163	xuv = 0d0)
164	call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165	xuv = 0.5d0)
166	call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167	xuv = 0.25d0)
168	end if test_grossismx
169
170	is2 = 0
171
172	IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
173	.or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
174	IF (clon <= 0.) THEN
175	is2 = 1
176
177	do while (rlonm025(is2) < - pi .and. is2 < iim)
178	is2 = is2 + 1
179	end do
180
181	if (rlonm025(is2) < - pi) then
182	print *, 'Rlonm025 plus petit que - pi !'
183	STOP 1
184	end if
185	ELSE
186	is2 = iim
187
188	do while (rlonm025(is2) > pi .and. is2 > 1)
189	is2 = is2 - 1
190	end do
191
192	if (rlonm025(is2) > pi) then
193	print *, 'Rlonm025 plus grand que pi !'
194	STOP 1
195	end if
196	END IF
197	END IF
198
199	call principal_cshift(is2, rlonm025, xprimm025)
200	call principal_cshift(is2, rlonv, xprimv)
201	call principal_cshift(is2, rlonu, xprimu)
202	call principal_cshift(is2, rlonp025, xprimp025)
203
204	forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205	print , "Minimum longitude step:", MINval(d_rlonv) 180. / pi, "degrees"
206	print , "Maximum longitude step:", MAXval(d_rlonv) 180. / pi, "degrees"
207
208	! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209	DO i = 1, iim + 1
210	IF (rlonp025(i) < rlonv(i)) THEN
211	print *, 'rlonp025(', i, ') = ', rlonp025(i)
212	print *, "< rlonv(", i, ") = ", rlonv(i)
213	STOP 1
214	END IF
215
216	IF (rlonv(i) < rlonm025(i)) THEN
217	print *, 'rlonv(', i, ') = ', rlonv(i)
218	print *, "< rlonm025(", i, ") = ", rlonm025(i)
219	STOP 1
220	END IF
221
222	IF (rlonp025(i) > rlonu(i)) THEN
223	print *, 'rlonp025(', i, ') = ', rlonp025(i)
224	print *, "> rlonu(", i, ") = ", rlonu(i)
225	STOP 1
226	END IF
227	END DO
228
229	END SUBROUTINE fxhyp
230
231	end module fxhyp_m