source: LMDZ5/trunk/libf/dyn3d_common/fxhyp_m.F90 @ 2311

Last change on this file since 2311 was 2228, checked in by lguez, 10 years ago

Correcting a problem from revision 2218. The type double precision
with option "-fdefault-real-8" of gfortran is promoted to 16-byte
precision and there is no specific procedure in arth with this
precision. Could not add a specific procedure in arth with double
precision because, with ifort, the option "-real-size 64" does not
promote the double precision, so that would make two identical
specific procedures in arth.

In module nrtype, replaced double precision by a parameterized real
kind so that the effective precision does not depend on a compiler
option.

In coefpoly, fxhyp, fyhyp and invert_zoom_x, use the parameterized
real kind defined in nrtype, instead of double precision.

Also, in module nrtype, removed unused derived types sprs2_sp and
sprs2_dp.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.4 KB
Line 
1module fxhyp_m
2
3  IMPLICIT NONE
4
5contains
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 arth_m, only: arth
21    use invert_zoom_x_m, only: invert_zoom_x, nmax
22    use nrtype, only: pi, pi_d, twopi, twopi_d, k8
23    use principal_cshift_m, only: principal_cshift
24
25    include "dimensions.h"
26    ! for iim
27
28    include "serre.h"
29    ! for clon, grossismx, dzoomx, taux
30
31    REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
32    real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
33
34    ! Local:
35    real rlonm025(iim + 1), rlonp025(iim + 1)
36    REAL dzoom, step
37    real d_rlonv(iim)
38    REAL(K8) xtild(0:2 * nmax)
39    REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
40    REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax)
41    REAL(K8) fa, fb
42    INTEGER i, is2
43    REAL(K8) xmoy, fxm
44
45    !----------------------------------------------------------------------
46
47    print *, "Call sequence information: fxhyp"
48
49    test_iim: if (iim==1) then
50       rlonv(1)=0.
51       rlonu(1)=pi
52       rlonv(2)=rlonv(1)+twopi
53       rlonu(2)=rlonu(1)+twopi
54
55       xprimm025(:)=1.
56       xprimv(:)=1.
57       xprimu(:)=1.
58       xprimp025(:)=1.
59    else test_iim
60       test_grossismx: if (grossismx == 1.) then
61          step = twopi / iim
62
63          xprimm025(:iim) = step
64          xprimp025(:iim) = step
65          xprimv(:iim) = step
66          xprimu(:iim) = step
67
68          rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim)
69          rlonm025(:iim) = rlonv(:iim) - 0.25 * step
70          rlonp025(:iim) = rlonv(:iim) + 0.25 * step
71          rlonu(:iim) = rlonv(:iim) + 0.5 * step
72       else test_grossismx
73          dzoom = dzoomx * twopi_d
74          xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
75
76          ! Compute fhyp:
77          DO i = nmax, 2 * nmax
78             fa = taux * (dzoom / 2. - xtild(i))
79             fb = xtild(i) * (pi_d - xtild(i))
80
81             IF (200. * fb < - fa) THEN
82                fhyp(i) = - 1.
83             ELSE IF (200. * fb < fa) THEN
84                fhyp(i) = 1.
85             ELSE
86                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
87                   IF (200. * fb + fa < 1e-10) THEN
88                      fhyp(i) = - 1.
89                   ELSE IF (200. * fb - fa < 1e-10) THEN
90                      fhyp(i) = 1.
91                   END IF
92                ELSE
93                   fhyp(i) = TANH(fa / fb)
94                END IF
95             END IF
96
97             IF (xtild(i) == 0.) fhyp(i) = 1.
98             IF (xtild(i) == pi_d) fhyp(i) = -1.
99          END DO
100
101          ! Calcul de beta
102
103          ffdx = 0.
104
105          DO i = nmax + 1, 2 * nmax
106             xmoy = 0.5 * (xtild(i-1) + xtild(i))
107             fa = taux * (dzoom / 2. - xmoy)
108             fb = xmoy * (pi_d - xmoy)
109
110             IF (200. * fb < - fa) THEN
111                fxm = - 1.
112             ELSE IF (200. * fb < fa) THEN
113                fxm = 1.
114             ELSE
115                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
116                   IF (200. * fb + fa < 1e-10) THEN
117                      fxm = - 1.
118                   ELSE IF (200. * fb - fa < 1e-10) THEN
119                      fxm = 1.
120                   END IF
121                ELSE
122                   fxm = TANH(fa / fb)
123                END IF
124             END IF
125
126             IF (xmoy == 0.) fxm = 1.
127             IF (xmoy == pi_d) fxm = -1.
128
129             ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
130          END DO
131
132          print *, "ffdx = ", ffdx
133          beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
134          print *, "beta = ", beta
135
136          IF (2. * beta - grossismx <= 0.) THEN
137             print *, 'Bad choice of grossismx, taux, dzoomx.'
138             print *, 'Decrease dzoomx or grossismx.'
139             STOP 1
140          END IF
141
142          ! calcul de Xprimt
143          Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
144          xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
145
146          ! Calcul de Xf
147
148          DO i = nmax + 1, 2 * nmax
149             xmoy = 0.5 * (xtild(i-1) + xtild(i))
150             fa = taux * (dzoom / 2. - xmoy)
151             fb = xmoy * (pi_d - xmoy)
152
153             IF (200. * fb < - fa) THEN
154                fxm = - 1.
155             ELSE IF (200. * fb < fa) THEN
156                fxm = 1.
157             ELSE
158                fxm = TANH(fa / fb)
159             END IF
160
161             IF (xmoy == 0.) fxm = 1.
162             IF (xmoy == pi_d) fxm = -1.
163             xxpr(i) = beta + (grossismx - beta) * fxm
164          END DO
165
166          xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
167
168          Xf(0) = - pi_d
169
170          DO i=1, 2 * nmax - 1
171             Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
172          END DO
173
174          Xf(2 * nmax) = pi_d
175
176          call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), &
177               xprimm025(:iim), xuv = - 0.25_k8)
178          call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
179               xuv = 0._k8)
180          call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
181               xuv = 0.5_k8)
182          call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
183               xprimp025(:iim), xuv = 0.25_k8)
184       end if test_grossismx
185
186       is2 = 0
187
188       IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
189            .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
190          IF (clon <= 0.) THEN
191             is2 = 1
192
193             do while (rlonm025(is2) < - pi .and. is2 < iim)
194                is2 = is2 + 1
195             end do
196
197             if (rlonm025(is2) < - pi) then
198                print *, 'Rlonm025 plus petit que - pi !'
199                STOP 1
200             end if
201          ELSE
202             is2 = iim
203
204             do while (rlonm025(is2) > pi .and. is2 > 1)
205                is2 = is2 - 1
206             end do
207
208             if (rlonm025(is2) > pi) then
209                print *, 'Rlonm025 plus grand que pi !'
210                STOP 1
211             end if
212          END IF
213       END IF
214
215       call principal_cshift(is2, rlonm025, xprimm025)
216       call principal_cshift(is2, rlonv, xprimv)
217       call principal_cshift(is2, rlonu, xprimu)
218       call principal_cshift(is2, rlonp025, xprimp025)
219
220       forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
221       print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, &
222            "degrees"
223       print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, &
224            "degrees"
225
226       ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
227       DO i = 1, iim + 1
228          IF (rlonp025(i) < rlonv(i)) THEN
229             print *, 'rlonp025(', i, ') = ', rlonp025(i)
230             print *, "< rlonv(", i, ") = ", rlonv(i)
231             STOP 1
232          END IF
233
234          IF (rlonv(i) < rlonm025(i)) THEN
235             print *, 'rlonv(', i, ') = ', rlonv(i)
236             print *, "< rlonm025(", i, ") = ", rlonm025(i)
237             STOP 1
238          END IF
239
240          IF (rlonp025(i) > rlonu(i)) THEN
241             print *, 'rlonp025(', i, ') = ', rlonp025(i)
242             print *, "> rlonu(", i, ") = ", rlonu(i)
243             STOP 1
244          END IF
245       END DO
246    end if test_iim
247
248  END SUBROUTINE fxhyp
249
250end module fxhyp_m
Note: See TracBrowser for help on using the repository browser.