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  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,*)
104 ' Le param. dzoomx pour fxhyp est trop petit ! L aug ,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,*)
179 ' ** Attention ! La valeur beta calculee dans la rou ,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