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

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

    print *, "Call sequence information: fxhyp"

    xzoom = clon * pi_d / 180d0

    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 * pi / 180., 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, xzoom, rlonm025(:iim), &
            xprimm025(:iim), xuv = - 0.25d0)
       call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonv(:iim), &
            xprimv(:iim), xuv = 0d0)
       call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonu(:iim), &
            xprimu(:iim), xuv = 0.5d0)
       call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, 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"

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