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

Annotation of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 144 - (hide annotations)
Wed Jun 10 16:46:46 2015 UTC (8 years, 11 months ago) by guez
Original Path: trunk/Sources/dyn3d/fxhyp.f
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 guez 78 module fxhyp_m
2 guez 3
3 guez 78 IMPLICIT NONE
4 guez 3
5 guez 78 contains
6 guez 3
7 guez 119 SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8 guez 3
9 guez 91 ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10 guez 119 ! Author: P. Le Van, from formulas by R. Sadourny
11 guez 3
12 guez 78 ! Calcule les longitudes et dérivées dans la grille du GCM pour
13 guez 119 ! une fonction f(x) à dérivée tangente hyperbolique.
14 guez 3
15 guez 121 ! Il vaut mieux avoir : grossismx \times dzoom < pi
16 guez 3
17 guez 120 ! Le premier point scalaire pour une grille regulière (grossismx =
18     ! 1., taux=0., clon=0.) est à - 180 degrés.
19    
20 guez 78 USE dimens_m, ONLY: iim
21 guez 139 use dynetat0_m, only: clon, grossismx, dzoomx, taux
22 guez 131 use invert_zoom_x_m, only: invert_zoom_x, nmax
23 guez 126 use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24 guez 124 use principal_cshift_m, only: principal_cshift
25 guez 144 use tanh_cautious_m, only: tanh_cautious
26 guez 3
27 guez 119 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
28     real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
29 guez 3
30 guez 91 ! Local:
31 guez 119 real rlonm025(iim + 1), rlonp025(iim + 1)
32 guez 126 REAL dzoom, step
33 guez 124 real d_rlonv(iim)
34 guez 121 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 guez 124 INTEGER i, is2
38 guez 144 DOUBLE PRECISION, dimension(nmax + 1:2 * nmax):: xmoy, fxm
39 guez 3
40 guez 91 !----------------------------------------------------------------------
41    
42 guez 120 print *, "Call sequence information: fxhyp"
43    
44 guez 126 test_grossismx: if (grossismx == 1.) then
45     step = twopi / iim
46 guez 78
47 guez 126 xprimm025(:iim) = step
48     xprimp025(:iim) = step
49     xprimv(:iim) = step
50     xprimu(:iim) = step
51    
52 guez 127 rlonv(:iim) = arth(- pi + clon, step, iim)
53 guez 126 rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54     rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55     rlonu(:iim) = rlonv(:iim) + 0.5 * step
56 guez 132 else test_grossismx
57 guez 126 dzoom = dzoomx * twopi_d
58     xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59 guez 144 forall (i = nmax + 1:2 * nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
60 guez 126
61     ! Compute fhyp:
62 guez 144 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 guez 126
68 guez 144 fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
69 guez 3
70 guez 126 ! Calcul de beta
71 guez 3
72 guez 126 ffdx = 0.
73 guez 3
74 guez 126 DO i = nmax + 1, 2 * nmax
75 guez 144 ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
76 guez 126 END DO
77 guez 3
78 guez 126 print *, "ffdx = ", ffdx
79 guez 144 beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
80 guez 126 print *, "beta = ", beta
81 guez 3
82 guez 126 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 guez 78
88 guez 126 ! calcul de Xprimt
89     Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
90     xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
91 guez 78
92 guez 126 ! Calcul de Xf
93 guez 78
94 guez 144 xxpr(nmax + 1:2 * nmax) = beta + (grossismx - beta) * fxm
95 guez 126 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
96 guez 3
97 guez 126 Xf(0) = - pi_d
98 guez 121
99 guez 126 DO i=1, 2 * nmax - 1
100     Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
101     END DO
102 guez 3
103 guez 126 Xf(2 * nmax) = pi_d
104 guez 3
105 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
106 guez 127 xuv = - 0.25d0)
107 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
108 guez 127 xuv = 0d0)
109 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
110 guez 127 xuv = 0.5d0)
111 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
112 guez 127 xuv = 0.25d0)
113 guez 126 end if test_grossismx
114 guez 123
115 guez 124 is2 = 0
116 guez 3
117 guez 124 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 guez 3
122 guez 124 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 guez 128 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
154 guez 119 DO i = 1, iim + 1
155     IF (rlonp025(i) < rlonv(i)) THEN
156 guez 121 print *, 'rlonp025(', i, ') = ', rlonp025(i)
157     print *, "< rlonv(", i, ") = ", rlonv(i)
158 guez 119 STOP 1
159     END IF
160    
161     IF (rlonv(i) < rlonm025(i)) THEN
162 guez 121 print *, 'rlonv(', i, ') = ', rlonv(i)
163     print *, "< rlonm025(", i, ") = ", rlonm025(i)
164 guez 119 STOP 1
165     END IF
166    
167     IF (rlonp025(i) > rlonu(i)) THEN
168 guez 120 print *, 'rlonp025(', i, ') = ', rlonp025(i)
169     print *, "> rlonu(", i, ") = ", rlonu(i)
170 guez 119 STOP 1
171     END IF
172     END DO
173    
174 guez 78 END SUBROUTINE fxhyp
175    
176     end module fxhyp_m
  ViewVC Help
Powered by ViewVC 1.1.21