/[lmdze]/trunk/Sources/dyn3d/fxhyp.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 144 - (show annotations)
Wed Jun 10 16:46:46 2015 UTC (8 years, 11 months ago) by guez
File size: 5491 byte(s) 
In procedure fxhyp, the convoluted computation of tanh(fa/fb) occurred
three times. Extracted it into a function. Also, the computation of
xmoy and fxm was repeated. So stored the values in arrays instead.

In procedure fxhyp, in the computation of fhyp, there were tests
xtild(i) == 0. and xtild(i) == pi_d. No use to do these tests at each
iteration. We now they are true for i == nmax and i == 2 * nmax,
respectively, and we know they are false for other values of
"i". Similarly, in the computations of ffdx and xxpr, there were the
tests xmoy == 0. and xmoy == pi_d, these could not be true.

Moved files from bibio to dyn3d, following LMDZ.
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 dynetat0_m, only: clon, grossismx, dzoomx, taux
22 use invert_zoom_x_m, only: invert_zoom_x, nmax
23 use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24 use principal_cshift_m, only: principal_cshift
25 use tanh_cautious_m, only: tanh_cautious
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 DOUBLE PRECISION xtild(0:2 * nmax)
35 DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
36 DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
37 INTEGER i, is2
38 DOUBLE PRECISION, dimension(nmax + 1:2 * nmax):: xmoy, fxm
39
40 !----------------------------------------------------------------------
41
42 print *, "Call sequence information: fxhyp"
43
44 test_grossismx: if (grossismx == 1.) then
45 step = twopi / iim
46
47 xprimm025(:iim) = step
48 xprimp025(:iim) = step
49 xprimv(:iim) = step
50 xprimu(:iim) = step
51
52 rlonv(:iim) = arth(- pi + clon, step, iim)
53 rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54 rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55 rlonu(:iim) = rlonv(:iim) + 0.5 * step
56 else test_grossismx
57 dzoom = dzoomx * twopi_d
58 xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59 forall (i = nmax + 1:2 * nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
60
61 ! Compute fhyp:
62 fhyp(nmax + 1:2 * nmax - 1) = tanh_cautious(taux * (dzoom / 2. &
63 - xtild(nmax + 1:2 * nmax - 1)), xtild(nmax + 1:2 * nmax - 1) &
64 * (pi_d - xtild(nmax + 1:2 * nmax - 1)))
65 fhyp(nmax) = 1d0
66 fhyp(2 * nmax) = -1d0
67
68 fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
69
70 ! Calcul de beta
71
72 ffdx = 0.
73
74 DO i = nmax + 1, 2 * nmax
75 ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
76 END DO
77
78 print *, "ffdx = ", ffdx
79 beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
80 print *, "beta = ", beta
81
82 IF (2. * beta - grossismx <= 0.) THEN
83 print *, 'Bad choice of grossismx, taux, dzoomx.'
84 print *, 'Decrease dzoomx or grossismx.'
85 STOP 1
86 END IF
87
88 ! calcul de Xprimt
89 Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
90 xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
91
92 ! Calcul de Xf
93
94 xxpr(nmax + 1:2 * nmax) = beta + (grossismx - beta) * fxm
95 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
96
97 Xf(0) = - pi_d
98
99 DO i=1, 2 * nmax - 1
100 Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
101 END DO
102
103 Xf(2 * nmax) = pi_d
104
105 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
106 xuv = - 0.25d0)
107 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
108 xuv = 0d0)
109 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
110 xuv = 0.5d0)
111 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
112 xuv = 0.25d0)
113 end if test_grossismx
114
115 is2 = 0
116
117 IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
118 .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
119 IF (clon <= 0.) THEN
120 is2 = 1
121
122 do while (rlonm025(is2) < - pi .and. is2 < iim)
123 is2 = is2 + 1
124 end do
125
126 if (rlonm025(is2) < - pi) then
127 print *, 'Rlonm025 plus petit que - pi !'
128 STOP 1
129 end if
130 ELSE
131 is2 = iim
132
133 do while (rlonm025(is2) > pi .and. is2 > 1)
134 is2 = is2 - 1
135 end do
136
137 if (rlonm025(is2) > pi) then
138 print *, 'Rlonm025 plus grand que pi !'
139 STOP 1
140 end if
141 END IF
142 END IF
143
144 call principal_cshift(is2, rlonm025, xprimm025)
145 call principal_cshift(is2, rlonv, xprimv)
146 call principal_cshift(is2, rlonu, xprimu)
147 call principal_cshift(is2, rlonp025, xprimp025)
148
149 forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
150 print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
151 print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
152
153 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
154 DO i = 1, iim + 1
155 IF (rlonp025(i) < rlonv(i)) THEN
156 print *, 'rlonp025(', i, ') = ', rlonp025(i)
157 print *, "< rlonv(", i, ") = ", rlonv(i)
158 STOP 1
159 END IF
160
161 IF (rlonv(i) < rlonm025(i)) THEN
162 print *, 'rlonv(', i, ') = ', rlonv(i)
163 print *, "< rlonm025(", i, ") = ", rlonm025(i)
164 STOP 1
165 END IF
166
167 IF (rlonp025(i) > rlonu(i)) THEN
168 print *, 'rlonp025(', i, ') = ', rlonp025(i)
169 print *, "> rlonu(", i, ") = ", rlonu(i)
170 STOP 1
171 END IF
172 END DO
173
174 END SUBROUTINE fxhyp
175
176 end module fxhyp_m
  ViewVC Help
Powered by ViewVC 1.1.21