source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90 @ 5195

Last change on this file since 5195 was 5160, checked in by abarral, 4 months ago

Put .h into modules
  • 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.1 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 lmdz_dimensions, ONLY: iim
21    USE lmdz_arth, ONLY: arth
22    USE invert_zoom_x_m, ONLY: invert_zoom_x, nmax
23    USE lmdz_physical_constants, ONLY: pi, pi_d, twopi, twopi_d, k8
24    USE principal_cshift_m, ONLY: principal_cshift
25    USE serre_mod, ONLY: clon, grossismx, dzoomx, taux
26
27    REAL, INTENT(OUT) :: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
28    REAL, INTENT(OUT) :: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
29
30    ! Local:
31    REAL rlonm025(iim + 1), rlonp025(iim + 1)
32    REAL dzoom, step
33    REAL d_rlonv(iim)
34    REAL(K8) xtild(0:2 * nmax)
35    REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
36    REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax)
37    REAL(K8) fa, fb
38    INTEGER i, is2
39    REAL(K8) xmoy, fxm
40
41    !----------------------------------------------------------------------
42
43    PRINT *, "Call sequence information: fxhyp"
44
45    test_iim: if (iim==1) THEN
46      rlonv(1) = 0.
47      rlonu(1) = pi
48      rlonv(2) = rlonv(1) + twopi
49      rlonu(2) = rlonu(1) + twopi
50
51      xprimm025(:) = 1.
52      xprimv(:) = 1.
53      xprimu(:) = 1.
54      xprimp025(:) = 1.
55    else test_iim
56      test_grossismx: if (grossismx == 1.) THEN
57        step = twopi / iim
58
59        xprimm025(:iim) = step
60        xprimp025(:iim) = step
61        xprimv(:iim) = step
62        xprimu(:iim) = step
63
64        rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim)
65        rlonm025(:iim) = rlonv(:iim) - 0.25 * step
66        rlonp025(:iim) = rlonv(:iim) + 0.25 * step
67        rlonu(:iim) = rlonv(:iim) + 0.5 * step
68      else test_grossismx
69        dzoom = dzoomx * twopi_d
70        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
71
72        ! Compute fhyp:
73        DO i = nmax, 2 * nmax
74          fa = taux * (dzoom / 2. - xtild(i))
75          fb = xtild(i) * (pi_d - xtild(i))
76
77          IF (200. * fb < - fa) THEN
78            fhyp(i) = - 1.
79          ELSE IF (200. * fb < fa) THEN
80            fhyp(i) = 1.
81          ELSE
82            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
83              IF (200. * fb + fa < 1e-10) THEN
84                fhyp(i) = - 1.
85              ELSE IF (200. * fb - fa < 1e-10) THEN
86                fhyp(i) = 1.
87              END IF
88            ELSE
89              fhyp(i) = TANH(fa / fb)
90            END IF
91          END IF
92
93          IF (xtild(i) == 0.) fhyp(i) = 1.
94          IF (xtild(i) == pi_d) fhyp(i) = -1.
95        END DO
96
97        ! Calcul de beta
98
99        ffdx = 0.
100
101        DO i = nmax + 1, 2 * nmax
102          xmoy = 0.5 * (xtild(i - 1) + xtild(i))
103          fa = taux * (dzoom / 2. - xmoy)
104          fb = xmoy * (pi_d - xmoy)
105
106          IF (200. * fb < - fa) THEN
107            fxm = - 1.
108          ELSE IF (200. * fb < fa) THEN
109            fxm = 1.
110          ELSE
111            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
112              IF (200. * fb + fa < 1e-10) THEN
113                fxm = - 1.
114              ELSE IF (200. * fb - fa < 1e-10) THEN
115                fxm = 1.
116              END IF
117            ELSE
118              fxm = TANH(fa / fb)
119            END IF
120          END IF
121
122          IF (xmoy == 0.) fxm = 1.
123          IF (xmoy == pi_d) fxm = -1.
124
125          ffdx = ffdx + fxm * (xtild(i) - xtild(i - 1))
126        END DO
127
128        PRINT *, "ffdx = ", ffdx
129        beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
130        PRINT *, "beta = ", beta
131
132        IF (2. * beta - grossismx <= 0.) THEN
133          PRINT *, 'Bad choice of grossismx, taux, dzoomx.'
134          PRINT *, 'Decrease dzoomx or grossismx.'
135          STOP 1
136        END IF
137
138        ! calcul de Xprimt
139        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
140        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
141
142        ! Calcul de Xf
143
144        DO i = nmax + 1, 2 * nmax
145          xmoy = 0.5 * (xtild(i - 1) + xtild(i))
146          fa = taux * (dzoom / 2. - xmoy)
147          fb = xmoy * (pi_d - xmoy)
148
149          IF (200. * fb < - fa) THEN
150            fxm = - 1.
151          ELSE IF (200. * fb < fa) THEN
152            fxm = 1.
153          ELSE
154            fxm = TANH(fa / fb)
155          END IF
156
157          IF (xmoy == 0.) fxm = 1.
158          IF (xmoy == pi_d) fxm = -1.
159          xxpr(i) = beta + (grossismx - beta) * fxm
160        END DO
161
162        xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
163
164        Xf(0) = - pi_d
165
166        DO i = 1, 2 * nmax - 1
167          Xf(i) = Xf(i - 1) + xxpr(i) * (xtild(i) - xtild(i - 1))
168        END DO
169
170        Xf(2 * nmax) = pi_d
171
172        CALL invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), &
173                xprimm025(:iim), xuv = - 0.25_k8)
174        CALL invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
175                xuv = 0._k8)
176        CALL invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
177                xuv = 0.5_k8)
178        CALL invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
179                xprimp025(:iim), xuv = 0.25_k8)
180      end if test_grossismx
181
182      is2 = 0
183
184      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
185              .OR. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
186        IF (clon <= 0.) THEN
187          is2 = 1
188
189          DO while (rlonm025(is2) < - pi .AND. is2 < iim)
190            is2 = is2 + 1
191          END DO
192
193          IF (rlonm025(is2) < - pi) THEN
194            PRINT *, 'Rlonm025 plus petit que - pi !'
195            STOP 1
196          end if
197        ELSE
198          is2 = iim
199
200          DO while (rlonm025(is2) > pi .AND. is2 > 1)
201            is2 = is2 - 1
202          END DO
203
204          IF (rlonm025(is2) > pi) THEN
205            PRINT *, 'Rlonm025 plus grand que pi !'
206            STOP 1
207          end if
208        END IF
209      END IF
210
211      CALL principal_cshift(is2, rlonm025, xprimm025)
212      CALL principal_cshift(is2, rlonv, xprimv)
213      CALL principal_cshift(is2, rlonu, xprimu)
214      CALL principal_cshift(is2, rlonp025, xprimp025)
215
216      forall (i = 1:iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
217      PRINT *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, &
218              "degrees"
219      PRINT *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, &
220              "degrees"
221
222      ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
223      DO i = 1, iim + 1
224        IF (rlonp025(i) < rlonv(i)) THEN
225          PRINT *, 'rlonp025(', i, ') = ', rlonp025(i)
226          PRINT *, "< rlonv(", i, ") = ", rlonv(i)
227          STOP 1
228        END IF
229
230        IF (rlonv(i) < rlonm025(i)) THEN
231          PRINT *, 'rlonv(', i, ') = ', rlonv(i)
232          PRINT *, "< rlonm025(", i, ") = ", rlonm025(i)
233          STOP 1
234        END IF
235
236        IF (rlonp025(i) > rlonu(i)) THEN
237          PRINT *, 'rlonp025(', i, ') = ', rlonp025(i)
238          PRINT *, "> rlonu(", i, ") = ", rlonu(i)
239          STOP 1
240        END IF
241      END DO
242    end if test_iim
243
244  END SUBROUTINE fxhyp
245
246END MODULE fxhyp_m
Note: See TracBrowser for help on using the repository browser.