My Project
 All Classes Files Functions Variables Macros
fxhyp.F
Go to the documentation of this file.
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,*)
81 ' Le param. dzoomx pour fxhyp est trop petit ! L aug ,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,*)
156 ' ** Attention ! La valeur beta calculee dans la rou ,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