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

Annotation of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 139 - (hide annotations)
Tue May 26 17:46:03 2015 UTC (8 years, 11 months ago) by guez
Original Path: trunk/Sources/dyn3d/fxhyp.f
File size: 6843 byte(s) 
dynetat0 read rlonu, rlatu, rlonv, rlatv, cu_2d, cv_2d, aire_2d from
"start.nc" and then these variables were overwritten by
inigeom. Corrected this. Now, inigeom does not compute rlonu, rlatu,
rlonv and rlatv. Moreover, cu_2d, cv_2d, aire_2d are not written to
"restart.nc". Since xprimu, xprimv, xprimm025, xprimp025, rlatu1,
rlatu2, yprimu1, yprimu2 are computed at the same time as rlonu,
rlatu, rlonv, rlatv, and since it would not be convenient to separate
those computations, we decide to write xprimu, xprimv, xprimm025,
xprimp025, rlatu1, rlatu2, yprimu1, yprimu2 into "restart.nc", read
them from "start.nc" and not compute them in inigeom. So, in summary,
"start.nc" contains all the coordinates and their derivatives, and
inigeom only computes the 2D-variables.

Technical details:

Moved variables rlatu, rlonv, rlonu, rlatv, xprimu, xprimv from module
comgeom to module dynetat0_m. Upgraded local variables rlatu1,
yprimu1, rlatu2, yprimu2, xprimm025, xprimp025 of procedure inigeom to
variables of module dynetat0_m.

Removed unused local variable yprimu of procedure inigeom and
corresponding argument yyprimu of fyhyp.

Moved variables clat, clon, grossismx, grossismy, dzoomx, dzoomy,
taux, tauy from module serre to module dynetat0_m (since they are read
from "start.nc"). The default values are now defined in read_serre
instead of in the declarations. Changed name of module serre to
read_serre_m, no more module variable here.

The calls to fxhyp and fyhyp are moved from inigeom to etat0.

Side effects in programs other than gcm: etat0 and read_serre write
variables of module dynetat0; the programs test_fxyp and
test_inter_barxy need more source files.

Removed unused arguments len and nd of cv3_tracer. Removed unused
argument PPSOL of LWU.

Bug fix in test_inter_barxy: forgotten call to read_serre.
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 3
26 guez 119 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27     real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
28 guez 3
29 guez 91 ! Local:
30 guez 119 real rlonm025(iim + 1), rlonp025(iim + 1)
31 guez 126 REAL dzoom, step
32 guez 124 real d_rlonv(iim)
33 guez 121 DOUBLE PRECISION xtild(0:2 * nmax)
34     DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35     DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36 guez 127 DOUBLE PRECISION fa, fb
37 guez 124 INTEGER i, is2
38 guez 121 DOUBLE PRECISION 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    
60     ! Compute fhyp:
61     DO i = nmax, 2 * nmax
62     fa = taux * (dzoom / 2. - xtild(i))
63     fb = xtild(i) * (pi_d - xtild(i))
64    
65     IF (200. * fb < - fa) THEN
66     fhyp(i) = - 1.
67     ELSE IF (200. * fb < fa) THEN
68     fhyp(i) = 1.
69     ELSE
70     IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71     IF (200. * fb + fa < 1e-10) THEN
72     fhyp(i) = - 1.
73     ELSE IF (200. * fb - fa < 1e-10) THEN
74     fhyp(i) = 1.
75     END IF
76     ELSE
77     fhyp(i) = TANH(fa / fb)
78 guez 119 END IF
79     END IF
80 guez 3
81 guez 126 IF (xtild(i) == 0.) fhyp(i) = 1.
82     IF (xtild(i) == pi_d) fhyp(i) = -1.
83     END DO
84 guez 3
85 guez 126 ! Calcul de beta
86 guez 3
87 guez 126 ffdx = 0.
88 guez 3
89 guez 126 DO i = nmax + 1, 2 * nmax
90     xmoy = 0.5 * (xtild(i-1) + xtild(i))
91     fa = taux * (dzoom / 2. - xmoy)
92     fb = xmoy * (pi_d - xmoy)
93 guez 78
94 guez 126 IF (200. * fb < - fa) THEN
95     fxm = - 1.
96     ELSE IF (200. * fb < fa) THEN
97     fxm = 1.
98     ELSE
99     IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100     IF (200. * fb + fa < 1e-10) THEN
101     fxm = - 1.
102     ELSE IF (200. * fb - fa < 1e-10) THEN
103     fxm = 1.
104     END IF
105     ELSE
106     fxm = TANH(fa / fb)
107 guez 119 END IF
108     END IF
109 guez 3
110 guez 126 IF (xmoy == 0.) fxm = 1.
111     IF (xmoy == pi_d) fxm = -1.
112 guez 3
113 guez 126 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114     END DO
115 guez 3
116 guez 126 print *, "ffdx = ", ffdx
117     beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118     print *, "beta = ", beta
119 guez 3
120 guez 126 IF (2. * beta - grossismx <= 0.) THEN
121     print *, 'Bad choice of grossismx, taux, dzoomx.'
122     print *, 'Decrease dzoomx or grossismx.'
123     STOP 1
124     END IF
125 guez 78
126 guez 126 ! calcul de Xprimt
127     Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128     xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129 guez 78
130 guez 126 ! Calcul de Xf
131 guez 78
132 guez 126 DO i = nmax + 1, 2 * nmax
133     xmoy = 0.5 * (xtild(i-1) + xtild(i))
134     fa = taux * (dzoom / 2. - xmoy)
135     fb = xmoy * (pi_d - xmoy)
136 guez 78
137 guez 126 IF (200. * fb < - fa) THEN
138     fxm = - 1.
139     ELSE IF (200. * fb < fa) THEN
140     fxm = 1.
141     ELSE
142     fxm = TANH(fa / fb)
143     END IF
144 guez 3
145 guez 126 IF (xmoy == 0.) fxm = 1.
146     IF (xmoy == pi_d) fxm = -1.
147     xxpr(i) = beta + (grossismx - beta) * fxm
148     END DO
149 guez 3
150 guez 126 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151 guez 3
152 guez 126 Xf(0) = - pi_d
153 guez 121
154 guez 126 DO i=1, 2 * nmax - 1
155     Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156     END DO
157 guez 3
158 guez 126 Xf(2 * nmax) = pi_d
159 guez 3
160 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161 guez 127 xuv = - 0.25d0)
162 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163 guez 127 xuv = 0d0)
164 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165 guez 127 xuv = 0.5d0)
166 guez 131 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167 guez 127 xuv = 0.25d0)
168 guez 126 end if test_grossismx
169 guez 123
170 guez 124 is2 = 0
171 guez 3
172 guez 124 IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
173     .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
174     IF (clon <= 0.) THEN
175     is2 = 1
176 guez 3
177 guez 124 do while (rlonm025(is2) < - pi .and. is2 < iim)
178     is2 = is2 + 1
179     end do
180    
181     if (rlonm025(is2) < - pi) then
182     print *, 'Rlonm025 plus petit que - pi !'
183     STOP 1
184     end if
185     ELSE
186     is2 = iim
187    
188     do while (rlonm025(is2) > pi .and. is2 > 1)
189     is2 = is2 - 1
190     end do
191    
192     if (rlonm025(is2) > pi) then
193     print *, 'Rlonm025 plus grand que pi !'
194     STOP 1
195     end if
196     END IF
197     END IF
198    
199     call principal_cshift(is2, rlonm025, xprimm025)
200     call principal_cshift(is2, rlonv, xprimv)
201     call principal_cshift(is2, rlonu, xprimu)
202     call principal_cshift(is2, rlonp025, xprimp025)
203    
204     forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205     print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
206     print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
207    
208 guez 128 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209 guez 119 DO i = 1, iim + 1
210     IF (rlonp025(i) < rlonv(i)) THEN
211 guez 121 print *, 'rlonp025(', i, ') = ', rlonp025(i)
212     print *, "< rlonv(", i, ") = ", rlonv(i)
213 guez 119 STOP 1
214     END IF
215    
216     IF (rlonv(i) < rlonm025(i)) THEN
217 guez 121 print *, 'rlonv(', i, ') = ', rlonv(i)
218     print *, "< rlonm025(", i, ") = ", rlonm025(i)
219 guez 119 STOP 1
220     END IF
221    
222     IF (rlonp025(i) > rlonu(i)) THEN
223 guez 120 print *, 'rlonp025(', i, ') = ', rlonp025(i)
224     print *, "> rlonu(", i, ") = ", rlonu(i)
225 guez 119 STOP 1
226     END IF
227     END DO
228    
229 guez 78 END SUBROUTINE fxhyp
230    
231     end module fxhyp_m
  ViewVC Help
Powered by ViewVC 1.1.21