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

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

Bug fix in fxhyp. There was some bad logic for the shifting of
longitude grids rlonuv, rlonv, rlonm025 and rlonp025 toward [- pi,
pi]. In some cases, one of the four grids was not shifted and the
others were. For example, you could see the bug with: iim = 48, clon =
20, grossismx = 3, dzoomx = 0.15. The bad logic involved the variable
is2 in the loop on ik = 1, 4. The loop included some tests depending
on ik. The whole thing was quite convoluted. Rewrote fxhyp. In
particular, replaced the loop on ik by a call to a new procedure,
invert_zoom_x. fxhyp.F was in dyn3d, it is replaced by fxhyp_m.F90
which is in dyn3d_common. Removed several arguments of fxhyp: zoom
parameters are accessed through serre.h; rlonm025 and rlonp025 were
not used in inigeom; min and max of longitude steps are written inside
fxhyp.

Some simplifications and modernizations in fyhyp. In particular,
several arguments are removed: zoom parameters are accessed through
serre.h; yprimv was not used in inigeom; min and max of latitude steps
are written inside fyhyp.

Removed now useless intermediary procedure fxyhyper.

Changed default value of dzoomx and dzoomy from 0 to 0.2. dzoom[xy] is
only needed when grossism[xy] > 1 and then it should be > 0.

dzoom[xy] can no longer be the extent of the zoom in degrees: it must
be expressed as the fraction of the total domain.
  • 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.5 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
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    DOUBLE PRECISION xtild(0:2 * nmax)
39    DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
40    DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
41    DOUBLE PRECISION fa, fb
42    INTEGER i, is2
43    DOUBLE PRECISION 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.25d0)
178          call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
179               xuv = 0d0)
180          call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
181               xuv = 0.5d0)
182          call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
183               xprimp025(:iim), xuv = 0.25d0)
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.