Context Navigation

fxhyp.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

File size: 11.7 KB

Line
1	!
2	! $Id: fxhyp.F 1403 2010-07-01 09:02:53Z fairhead $
3	!
4	c
5	c
6	SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau ,
7	, rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
8	, champmin,champmax )
9
10	c Auteur : P. Le Van
11
12	IMPLICIT NONE
13
14	c Calcule les longitudes et derivees dans la grille du GCM pour une
15	c fonction f(x) a tangente hyperbolique .
16	c
17	c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
18	c dzoom etant la distance totale de la zone du zoom
19	c tau la raideur de la transition de l'interieur a l'exterieur du zoom
20	c
21	c On doit avoir grossism x dzoom < pi ( radians ) , en longitude.
22	c ********************************************************************
23
24
25	INTEGER nmax, nmax2
26	PARAMETER ( nmax = 30000, nmax2 = 2*nmax )
27	c
28	LOGICAL scal180
29	PARAMETER ( scal180 = .TRUE. )
30
31	c scal180 = .TRUE. si on veut avoir le premier point scalaire pour
32	c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
33	c sinon scal180 = .FALSE.
34
35	#include "dimensions.h"
36	#include "paramet.h"
37
38	c ...... arguments d'entree .......
39	c
40	REAL xzoomdeg,dzooma,tau,grossism
41
42	c ...... arguments de sortie ......
43
44	REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
45	, rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
46
47	c .... variables locales ....
48	c
49	REAL dzoom
50	REAL*8 xlon(iip1),xprimm(iip1),xuv
51	REAL*8 xtild(0:nmax2)
52	REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
53	REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
54	REAL*8 xvrai(iip1),xxprim(iip1)
55	REAL*8 pi,depi,epsilon,xzoom,fa,fb
56	REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
57	INTEGER i,it,ik,iter,ii,idif,ii1,ii2
58	REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
59	REAL*8 champmin,champmax,decalx
60	INTEGER is2
61	SAVE is2
62
63	REAL*8 heavyside
64
65	pi = 2. * ASIN(1.)
66	depi = 2. * pi
67	epsilon = 1.e-3
68	xzoom = xzoomdeg * pi/180.
69	c
70	decalx = .75
71	IF( grossism.EQ.1..AND.scal180 ) THEN
72	decalx = 1.
73	ENDIF
74
75	WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx
76	c
77	IF( dzooma.LT.1.) THEN
78	dzoom = dzooma * depi
79	ELSEIF( dzooma.LT. 25. ) THEN
80	WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
81	,menter et relancer ! '
82	STOP 1
83	ELSE
84	dzoom = dzooma * pi/180.
85	ENDIF
86
87	WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
88	WRITE(6,24) xzoom,grossism,tau,dzoom
89
90	DO i = 0, nmax2
91	xtild(i) = - pi + REAL(i) * depi /nmax2
92	ENDDO
93
94	DO i = nmax, nmax2
95
96	fa = tau* ( dzoom/2. - xtild(i) )
97	fb = xtild(i) * ( pi - xtild(i) )
98
99	IF( 200.* fb .LT. - fa ) THEN
100	fhyp ( i) = - 1.
101	ELSEIF( 200. * fb .LT. fa ) THEN
102	fhyp ( i) = 1.
103	ELSE
104	IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
105	IF( 200.*fb + fa.LT.1.e-10 ) THEN
106	fhyp ( i ) = - 1.
107	ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN
108	fhyp ( i ) = 1.
109	ENDIF
110	ELSE
111	fhyp ( i ) = TANH ( fa/fb )
112	ENDIF
113	ENDIF
114	IF ( xtild(i).EQ. 0. ) fhyp(i) = 1.
115	IF ( xtild(i).EQ. pi ) fhyp(i) = -1.
116
117	ENDDO
118
119	cc .... Calcul de beta ....
120
121	ffdx = 0.
122
123	DO i = nmax +1,nmax2
124
125	xmoy = 0.5 * ( xtild(i-1) + xtild( i ) )
126	fa = tau* ( dzoom/2. - xmoy )
127	fb = xmoy * ( pi - xmoy )
128
129	IF( 200.* fb .LT. - fa ) THEN
130	fxm = - 1.
131	ELSEIF( 200. * fb .LT. fa ) THEN
132	fxm = 1.
133	ELSE
134	IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
135	IF( 200.*fb + fa.LT.1.e-10 ) THEN
136	fxm = - 1.
137	ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN
138	fxm = 1.
139	ENDIF
140	ELSE
141	fxm = TANH ( fa/fb )
142	ENDIF
143	ENDIF
144
145	IF ( xmoy.EQ. 0. ) fxm = 1.
146	IF ( xmoy.EQ. pi ) fxm = -1.
147
148	ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
149
150	ENDDO
151
152	beta = ( grossism * ffdx - pi ) / ( ffdx - pi )
153
154	IF( 2.*beta - grossism.LE. 0.) THEN
155	WRITE(6,) ' * Attention ! La valeur beta calculee dans la rou
156	,tine fxhyp est mauvaise ! '
157	WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ',
158	, ' et relancer ! *** '
159	CALL ABORT
160	ENDIF
161	c
162	c ..... calcul de Xprimt .....
163	c
164
165	DO i = nmax, nmax2
166	Xprimt(i) = beta + ( grossism - beta ) * fhyp(i)
167	ENDDO
168	c
169	DO i = nmax+1, nmax2
170	Xprimt( nmax2 - i ) = Xprimt( i )
171	ENDDO
172	c
173
174	c ..... Calcul de Xf ........
175
176	Xf(0) = - pi
177
178	DO i = nmax +1, nmax2
179
180	xmoy = 0.5 * ( xtild(i-1) + xtild( i ) )
181	fa = tau* ( dzoom/2. - xmoy )
182	fb = xmoy * ( pi - xmoy )
183
184	IF( 200.* fb .LT. - fa ) THEN
185	fxm = - 1.
186	ELSEIF( 200. * fb .LT. fa ) THEN
187	fxm = 1.
188	ELSE
189	fxm = TANH ( fa/fb )
190	ENDIF
191
192	IF ( xmoy.EQ. 0. ) fxm = 1.
193	IF ( xmoy.EQ. pi ) fxm = -1.
194	xxpr(i) = beta + ( grossism - beta ) * fxm
195
196	ENDDO
197
198	DO i = nmax+1, nmax2
199	xxpr(nmax2-i+1) = xxpr(i)
200	ENDDO
201
202	DO i=1,nmax2
203	Xf(i) = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
204	ENDDO
205
206
207	c *****************************************************************
208	c
209
210	c ..... xuv = 0. si calcul aux pts scalaires ........
211	c ..... xuv = 0.5 si calcul aux pts U ........
212	c
213	WRITE(6,18)
214	c
215	DO 5000 ik = 1, 4
216
217	IF( ik.EQ.1 ) THEN
218	xuv = -0.25
219	ELSE IF ( ik.EQ.2 ) THEN
220	xuv = 0.
221	ELSE IF ( ik.EQ.3 ) THEN
222	xuv = 0.50
223	ELSE IF ( ik.EQ.4 ) THEN
224	xuv = 0.25
225	ENDIF
226
227	xo1 = 0.
228
229	ii1=1
230	ii2=iim
231	IF(ik.EQ.1.and.grossism.EQ.1.) THEN
232	ii1 = 2
233	ii2 = iim+1
234	ENDIF
235	DO 1500 i = ii1, ii2
236
237	xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim)
238
239	Xfi = xlon2
240	c
241	DO 250 it = nmax2,0,-1
242	IF( Xfi.GE.Xf(it)) GO TO 350
243	250 CONTINUE
244
245	it = 0
246
247	350 CONTINUE
248
249	c ...... Calcul de Xf(xi) ......
250	c
251	xi = xtild(it)
252
253	IF(it.EQ.nmax2) THEN
254	it = nmax2 -1
255	Xf(it+1) = pi
256	ENDIF
257	c .....................................................................
258	c
259	c Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
260	c polynome de degre 3 qui passe par les points (Xf(it),xtild(it) )
261	c et (Xf(it+1),xtild(it+1) )
262
263	CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
264	, xtild(it),xtild(it+1), a0, a1, a2, a3 )
265
266	Xf1 = Xf(it)
267	Xprimin = a1 + 2.* a2 * xi + 3.a3 xi *xi
268
269	DO 500 iter = 1,300
270	xi = xi - ( Xf1 - Xfi )/ Xprimin
271
272	IF( ABS(xi-xo1).LE.epsilon) GO TO 550
273	xo1 = xi
274	xi2 = xi * xi
275	Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
276	Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
277	500 CONTINUE
278	WRITE(6,) ' Pas de solution **** ',i,xlon2,iter
279	STOP 6
280	550 CONTINUE
281
282	xxprim(i) = depi/ ( REAL(iim) * Xprimin )
283	xvrai(i) = xi + xzoom
284
285	1500 CONTINUE
286
287
288	IF(ik.EQ.1.and.grossism.EQ.1.) THEN
289	xvrai(1) = xvrai(iip1)-depi
290	xxprim(1) = xxprim(iip1)
291	ENDIF
292	DO i = 1 , iim
293	xlon(i) = xvrai(i)
294	xprimm(i) = xxprim(i)
295	ENDDO
296	DO i = 1, iim -1
297	IF( xvrai(i+1). LT. xvrai(i) ) THEN
298	WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
299	, ')'
300	STOP 7
301	ENDIF
302	ENDDO
303	c
304	c ... Reorganisation des longitudes pour les avoir entre - pi et pi ..
305	c ........................................................................
306
307	champmin = 1.e12
308	champmax = -1.e12
309	DO i = 1, iim
310	champmin = MIN( champmin,xvrai(i) )
311	champmax = MAX( champmax,xvrai(i) )
312	ENDDO
313
314	IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN
315	GO TO 1600
316	ELSE
317	WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
318	, ' et pi '
319	c
320	IF( xzoom.LE.0.) THEN
321	IF( ik.EQ. 1 ) THEN
322	DO i = 1, iim
323	IF( xvrai(i).GE. - pi ) GO TO 80
324	ENDDO
325	WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! '
326	STOP 8
327	80 CONTINUE
328	is2 = i
329	ENDIF
330
331	IF( is2.NE. 1 ) THEN
332	DO ii = is2 , iim
333	xlon (ii-is2+1) = xvrai(ii)
334	xprimm(ii-is2+1) = xxprim(ii)
335	ENDDO
336	DO ii = 1 , is2 -1
337	xlon (ii+iim-is2+1) = xvrai(ii) + depi
338	xprimm(ii+iim-is2+1) = xxprim(ii)
339	ENDDO
340	ENDIF
341	ELSE
342	IF( ik.EQ.1 ) THEN
343	DO i = iim,1,-1
344	IF( xvrai(i).LE. pi ) GO TO 90
345	ENDDO
346	WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! '
347	STOP 9
348	90 CONTINUE
349	is2 = i
350	ENDIF
351	idif = iim -is2
352	DO ii = 1, is2
353	xlon (ii+idif) = xvrai(ii)
354	xprimm(ii+idif) = xxprim(ii)
355	ENDDO
356	DO ii = 1, idif
357	xlon (ii) = xvrai (ii+is2) - depi
358	xprimm(ii) = xxprim(ii+is2)
359	ENDDO
360	ENDIF
361	ENDIF
362	c
363	c ......... Fin de la reorganisation ............................
364
365	1600 CONTINUE
366
367
368	xlon ( iip1) = xlon(1) + depi
369	xprimm( iip1 ) = xprimm (1 )
370
371	DO i = 1, iim+1
372	xvrai(i) = xlon(i)*180./pi
373	ENDDO
374
375	IF( ik.EQ.1 ) THEN
376	c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) '
377	c WRITE(6,18)
378	c WRITE(6,68) xvrai
379	c WRITE(6,*) ' XPRIM k ',ik
380	c WRITE(6,566) xprimm
381
382	DO i = 1,iim +1
383	rlonm025(i) = xlon( i )
384	xprimm025(i) = xprimm(i)
385	ENDDO
386	ELSE IF( ik.EQ.2 ) THEN
387	c WRITE(6,18)
388	c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) '
389	c WRITE(6,68) xvrai
390	c WRITE(6,*) ' XPRIM k ',ik
391	c WRITE(6,566) xprimm
392
393	DO i = 1,iim + 1
394	rlonv(i) = xlon( i )
395	xprimv(i) = xprimm(i)
396	ENDDO
397
398	ELSE IF( ik.EQ.3) THEN
399	c WRITE(6,18)
400	c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) '
401	c WRITE(6,68) xvrai
402	c WRITE(6,*) ' XPRIM ik ',ik
403	c WRITE(6,566) xprimm
404
405	DO i = 1,iim + 1
406	rlonu(i) = xlon( i )
407	xprimu(i) = xprimm(i)
408	ENDDO
409
410	ELSE IF( ik.EQ.4 ) THEN
411	c WRITE(6,18)
412	c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) '
413	c WRITE(6,68) xvrai
414	c WRITE(6,*) ' XPRIM ik ',ik
415	c WRITE(6,566) xprimm
416
417	DO i = 1,iim + 1
418	rlonp025(i) = xlon( i )
419	xprimp025(i) = xprimm(i)
420	ENDDO
421
422	ENDIF
423
424	5000 CONTINUE
425	c
426	WRITE(6,18)
427	c
428	c ........... fin de la boucle do 5000 ............
429
430	DO i = 1, iim
431	xlon(i) = rlonv(i+1) - rlonv(i)
432	ENDDO
433	champmin = 1.e12
434	champmax = -1.e12
435	DO i = 1, iim
436	champmin = MIN( champmin, xlon(i) )
437	champmax = MAX( champmax, xlon(i) )
438	ENDDO
439	champmin = champmin * 180./pi
440	champmax = champmax * 180./pi
441
442	18 FORMAT(/)
443	24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
444	68 FORMAT(1x,7f9.2)
445	566 FORMAT(1x,7f9.4)
446
447	RETURN
448	END

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

Download in other formats:

Original Format