Context Navigation

fxhyp_m.F90 @ 5114

Last change on this file since 5114 was 5113, checked in by abarral, 4 months ago
Rename modules in misc from _mod > lmdz_ Put cbrt.f90, ch.f90, pch.f90 in new lmdz_libmath_pch.f90
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 lmdz_arth, 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
248	end module fxhyp_m

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

Download in other formats:

Original Format