source: trunk/LMDZ.COMMON/libf/dyn3d_common/fxhyp_m.F90 @ 1549

Last change on this file since 1549 was 1443, checked in by emillour, 10 years ago

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