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