Context Navigation

fxhyp_m.f90 @ 5319

Last change on this file since 5319 was 5271, checked in by abarral, 2 weeks ago
Move dimensions.h into a module Nb: doesn't compile yet
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.4 KB

Line
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 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	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
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
246	end module fxhyp_m

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format