LMDZ
fxhyp.F
Go to the documentation of this file.
1 !
2 ! $Id: fxhyp.F 1939 2014-01-21 14:23:17Z lguez $
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,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
102  ,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,*) ' ** Attention ! La valeur beta calculee dans la rou
177  ,tine fxhyp est mauvaise ! '
178  WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ',
179  , ' et relancer ! *** '
180  CALL abort_gcm("FXHYP", "", 1)
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  IF(ik.EQ.1.and.grossism.EQ.1.) THEN
310  xvrai(1) = xvrai(iip1)-depi
311  xxprim(1) = xxprim(iip1)
312  ENDIF
313 
314  DO i = 1 , iim
315  xlon(i) = xvrai(i)
316  xprimm(i) = xxprim(i)
317  ENDDO
318  DO i = 1, iim -1
319  IF( xvrai(i+1). lt. xvrai(i) ) THEN
320  WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
321  , ')'
322  stop 7
323  ENDIF
324  ENDDO
325 c
326 c ... Reorganisation des longitudes pour les avoir entre - pi et pi ..
327 c ........................................................................
328 
329  champmin = 1.e12
330  champmax = -1.e12
331  DO i = 1, iim
332  champmin = min( champmin,xvrai(i) )
333  champmax = max( champmax,xvrai(i) )
334  ENDDO
335 
336  IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN
337  GO TO 1600
338  ELSE
339  WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
340  , ' et pi '
341 c
342  IF( xzoom.LE.0.) THEN
343  IF( ik.EQ. 1 ) THEN
344  DO i = 1, iim
345  IF( xvrai(i).GE. - pi ) GO TO 80
346  ENDDO
347  WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! '
348  stop 8
349  80 CONTINUE
350  is2 = i
351  ENDIF
352 
353  IF( is2.NE. 1 ) THEN
354  DO ii = is2 , iim
355  xlon(ii-is2+1) = xvrai(ii)
356  xprimm(ii-is2+1) = xxprim(ii)
357  ENDDO
358  DO ii = 1 , is2 -1
359  xlon(ii+iim-is2+1) = xvrai(ii) + depi
360  xprimm(ii+iim-is2+1) = xxprim(ii)
361  ENDDO
362  ENDIF
363  ELSE
364  IF( ik.EQ.1 ) THEN
365  DO i = iim,1,-1
366  IF( xvrai(i).LE. pi ) GO TO 90
367  ENDDO
368  WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! '
369  stop 9
370  90 CONTINUE
371  is2 = i
372  ENDIF
373  idif = iim -is2
374  DO ii = 1, is2
375  xlon(ii+idif) = xvrai(ii)
376  xprimm(ii+idif) = xxprim(ii)
377  ENDDO
378  DO ii = 1, idif
379  xlon(ii) = xvrai(ii+is2) - depi
380  xprimm(ii) = xxprim(ii+is2)
381  ENDDO
382  ENDIF
383  ENDIF
384 c
385 c ......... Fin de la reorganisation ............................
386 
387  1600 CONTINUE
388 
389 
390  xlon( iip1) = xlon(1) + depi
391  xprimm( iip1 ) = xprimm(1 )
392 
393  DO i = 1, iim+1
394  xvrai(i) = xlon(i)*180./pi
395  ENDDO
396 
397  IF( ik.EQ.1 ) THEN
398 c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) '
399 c WRITE(6,18)
400 c WRITE(6,68) xvrai
401 c WRITE(6,*) ' XPRIM k ',ik
402 c WRITE(6,566) xprimm
403 
404  DO i = 1,iim +1
405  rlonm025(i) = xlon( i )
406  xprimm025(i) = xprimm(i)
407  ENDDO
408  ELSE IF( ik.EQ.2 ) THEN
409 c WRITE(6,18)
410 c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) '
411 c WRITE(6,68) xvrai
412 c WRITE(6,*) ' XPRIM k ',ik
413 c WRITE(6,566) xprimm
414 
415  DO i = 1,iim + 1
416  rlonv(i) = xlon( i )
417  xprimv(i) = xprimm(i)
418  ENDDO
419 
420  ELSE IF( ik.EQ.3) THEN
421 c WRITE(6,18)
422 c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) '
423 c WRITE(6,68) xvrai
424 c WRITE(6,*) ' XPRIM ik ',ik
425 c WRITE(6,566) xprimm
426 
427  DO i = 1,iim + 1
428  rlonu(i) = xlon( i )
429  xprimu(i) = xprimm(i)
430  ENDDO
431 
432  ELSE IF( ik.EQ.4 ) THEN
433 c WRITE(6,18)
434 c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) '
435 c WRITE(6,68) xvrai
436 c WRITE(6,*) ' XPRIM ik ',ik
437 c WRITE(6,566) xprimm
438 
439  DO i = 1,iim + 1
440  rlonp025(i) = xlon( i )
441  xprimp025(i) = xprimm(i)
442  ENDDO
443 
444  ENDIF
445 
446 5000 CONTINUE
447 c
448  WRITE(6,18)
449 c
450 c ........... fin de la boucle do 5000 ............
451 
452  DO i = 1, iim
453  xlon(i) = rlonv(i+1) - rlonv(i)
454  ENDDO
455  champmin = 1.e12
456  champmax = -1.e12
457  DO i = 1, iim
458  champmin = min( champmin, xlon(i) )
459  champmax = max( champmax, xlon(i) )
460  ENDDO
461  champmin = champmin * 180./pi
462  champmax = champmax * 180./pi
463 
464 18 FORMAT(/)
465 24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
466 68 FORMAT(1x,7f9.2)
467 566 FORMAT(1x,7f9.4)
468 
469  RETURN
470  END
xprimv
!$Header!CDK comgeom COMMON comgeom xprimv(iip1)!REAL &&cu
abort_gcm
subroutine abort_gcm(modname, message, ierr)
Definition: abort_gcm.F:7
xprimu
!$Header!CDK comgeom2 COMMON comgeom unsaire xprimu
Definition: comgeom2.h:26
parameter
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
rlonu
!$Header!CDK comgeom COMMON comgeom rlonu
Definition: comgeom.h:25
x
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
rlonv
!$Header!CDK comgeom COMMON comgeom rlonv
Definition: comgeom.h:25
fxhyp
subroutine fxhyp(xzoomdeg, grossism, dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu, rlonp025, xprimp025, champmin, champmax)
Definition: fxhyp.F:9