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, 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
       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 fxhyp_loop_ik(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
            xuv = - 0.25d0)
       call fxhyp_loop_ik(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
            xuv = 0d0)
       call fxhyp_loop_ik(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
            xuv = 0.5d0)
       call fxhyp_loop_ik(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	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	126	use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23	guez	124	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	126	REAL dzoom, step
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	guez	127	DOUBLE PRECISION fa, fb
37	guez	124	INTEGER i, is2
38	guez	121	DOUBLE PRECISION xmoy, fxm
39	guez	3
40	guez	91	!----------------------------------------------------------------------
41
42	guez	120	print *, "Call sequence information: fxhyp"
43
44	guez	126	test_grossismx: if (grossismx == 1.) then
45			step = twopi / iim
46	guez	78
47	guez	126	xprimm025(:iim) = step
48			xprimp025(:iim) = step
49			xprimv(:iim) = step
50			xprimu(:iim) = step
51
52	guez	127	rlonv(:iim) = arth(- pi + clon, step, iim)
53	guez	126	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
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	guez	119	END IF
79			END IF
80	guez	3
81	guez	126	IF (xtild(i) == 0.) fhyp(i) = 1.
82			IF (xtild(i) == pi_d) fhyp(i) = -1.
83			END DO
84	guez	3
85	guez	126	! Calcul de beta
86	guez	3
87	guez	126	ffdx = 0.
88	guez	3
89	guez	126	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	guez	78
94	guez	126	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	guez	119	END IF
108			END IF
109	guez	3
110	guez	126	IF (xmoy == 0.) fxm = 1.
111			IF (xmoy == pi_d) fxm = -1.
112	guez	3
113	guez	126	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114			END DO
115	guez	3
116	guez	126	print *, "ffdx = ", ffdx
117			beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118			print *, "beta = ", beta
119	guez	3
120	guez	126	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	guez	78
126	guez	126	! calcul de Xprimt
127			Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128			xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129	guez	78
130	guez	126	! Calcul de Xf
131	guez	78
132	guez	126	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	guez	78
137	guez	126	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	guez	3
145	guez	126	IF (xmoy == 0.) fxm = 1.
146			IF (xmoy == pi_d) fxm = -1.
147			xxpr(i) = beta + (grossismx - beta) * fxm
148			END DO
149	guez	3
150	guez	126	xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151	guez	3
152	guez	126	Xf(0) = - pi_d
153	guez	121
154	guez	126	DO i=1, 2 * nmax - 1
155			Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156			END DO
157	guez	3
158	guez	126	Xf(2 * nmax) = pi_d
159	guez	3
160	guez	127	call fxhyp_loop_ik(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161			xuv = - 0.25d0)
162			call fxhyp_loop_ik(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163			xuv = 0d0)
164			call fxhyp_loop_ik(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165			xuv = 0.5d0)
166			call fxhyp_loop_ik(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167			xuv = 0.25d0)
168	guez	126	end if test_grossismx
169	guez	123
170	guez	124	is2 = 0
171	guez	3
172	guez	124	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	guez	3
177	guez	124	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	guez	128	! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209	guez	119	DO i = 1, iim + 1
210			IF (rlonp025(i) < rlonv(i)) THEN
211	guez	121	print *, 'rlonp025(', i, ') = ', rlonp025(i)
212			print *, "< rlonv(", i, ") = ", rlonv(i)
213	guez	119	STOP 1
214			END IF
215
216			IF (rlonv(i) < rlonm025(i)) THEN
217	guez	121	print *, 'rlonv(', i, ') = ', rlonv(i)
218			print *, "< rlonm025(", i, ") = ", rlonm025(i)
219	guez	119	STOP 1
220			END IF
221
222			IF (rlonp025(i) > rlonu(i)) THEN
223	guez	120	print *, 'rlonp025(', i, ') = ', rlonp025(i)
224			print *, "> rlonu(", i, ") = ", rlonu(i)
225	guez	119	STOP 1
226			END IF
227			END DO
228
229	guez	78	END SUBROUTINE fxhyp
230
231			end module fxhyp_m