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	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, 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 xzoom, fa, fb
37	INTEGER i, is2
38	DOUBLE PRECISION xmoy, fxm
39
40	!----------------------------------------------------------------------
41
42	print *, "Call sequence information: fxhyp"
43
44	xzoom = clon * pi_d / 180d0
45
46	test_grossismx: if (grossismx == 1.) then
47	step = twopi / iim
48
49	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	END IF
81	END IF
82
83	IF (xtild(i) == 0.) fhyp(i) = 1.
84	IF (xtild(i) == pi_d) fhyp(i) = -1.
85	END DO
86
87	! Calcul de beta
88
89	ffdx = 0.
90
91	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
96	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	END IF
110	END IF
111
112	IF (xmoy == 0.) fxm = 1.
113	IF (xmoy == pi_d) fxm = -1.
114
115	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
116	END DO
117
118	print *, "ffdx = ", ffdx
119	beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
120	print *, "beta = ", beta
121
122	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
128	! calcul de Xprimt
129	Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130	xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131
132	! Calcul de Xf
133
134	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
139	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
147	IF (xmoy == 0.) fxm = 1.
148	IF (xmoy == pi_d) fxm = -1.
149	xxpr(i) = beta + (grossismx - beta) * fxm
150	END DO
151
152	xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153
154	Xf(0) = - pi_d
155
156	DO i=1, 2 * nmax - 1
157	Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
158	END DO
159
160	Xf(2 * nmax) = pi_d
161
162	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
172	is2 = 0
173
174	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
179	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	DO i = 1, iim + 1
211	IF (rlonp025(i) < rlonv(i)) THEN
212	print *, 'rlonp025(', i, ') = ', rlonp025(i)
213	print *, "< rlonv(", i, ") = ", rlonv(i)
214	STOP 1
215	END IF
216
217	IF (rlonv(i) < rlonm025(i)) THEN
218	print *, 'rlonv(', i, ') = ', rlonv(i)
219	print *, "< rlonm025(", i, ") = ", rlonm025(i)
220	STOP 1
221	END IF
222
223	IF (rlonp025(i) > rlonu(i)) THEN
224	print *, 'rlonp025(', i, ') = ', rlonp025(i)
225	print *, "> rlonu(", i, ") = ", rlonu(i)
226	STOP 1
227	END IF
228	END DO
229
230	END SUBROUTINE fxhyp
231
232	end module fxhyp_m