Context Navigation

fxhyp.F @ 4100

Last change on this file since 4100 was 1674, checked in by Ehouarn Millour, 12 years ago

Modification pour activation du 2D latitude-pression
(pour pouvoir compiler en -d 1xjmxlm)
dyn3d/fxhyp.F : calcul des longitudes a la main pour iim=1
dyn3d/groupe.F : desactive si iim=1
dyn3d/paramet.h : iip1=iim+1 au lieu de iim+1-1/iim precedemment
phylmd/iophy.F90 : on enleve les -1/iim
phylmd/phyetat0.F90 : on enleve les -1/iim

Modification for activation of the 3D latitude-pressure version
(to be compiled with -d 1xjmxlm)
dyn3d/fxhyp.F : longitudes imposed for iim=1
dyn3d/groupe.F : desactived when iim=1
dyn3d/paramet.h : iip1=iim+1 instead of iim+1-1/iim previously
phylmd/iophy.F90 : -1/iim removed
phylmd/phyetat0.F90 : -1/iim removed

FH et EM

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 12.2 KB

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

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

Download in other formats:

Original Format