/[lmdze]/trunk/dyn3d/fxhyp.f
ViewVC logotype

Annotation of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 76 - (hide annotations)
Fri Nov 15 18:45:49 2013 UTC (10 years, 6 months ago) by guez
File size: 12115 byte(s)
Moved everything out of libf.
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/fxhyp.F,v 1.2 2005/06/03 09:11:32 fairhead Exp $
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     use dimens_m
13     use paramet_m
14     IMPLICIT NONE
15    
16     c Calcule les longitudes et derivees dans la grille du GCM pour une
17     c fonction f(x) a tangente hyperbolique .
18     c
19     c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
20     c dzoom etant la distance totale de la zone du zoom
21     c tau la raideur de la transition de l'interieur a l'exterieur du zoom
22     c
23     c On doit avoir grossism x dzoom < pi ( radians ) , en longitude.
24     c ********************************************************************
25    
26    
27     INTEGER nmax, nmax2
28     PARAMETER ( nmax = 30000, nmax2 = 2*nmax )
29     c
30     LOGICAL scal180
31     PARAMETER ( scal180 = .TRUE. )
32    
33     c scal180 = .TRUE. si on veut avoir le premier point scalaire pour
34     c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
35     c sinon scal180 = .FALSE.
36    
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 guez 71 DOUBLE PRECISION xlon(iip1),xprimm(iip1),xuv
51     DOUBLE PRECISION xtild(0:nmax2)
52     DOUBLE PRECISION fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
53     DOUBLE PRECISION Xf(0:nmax2),xxpr(0:nmax2)
54     DOUBLE PRECISION xvrai(iip1),xxprim(iip1)
55     DOUBLE PRECISION pi,depi,epsilon,xzoom,fa,fb
56     DOUBLE PRECISION Xf1, Xfi , a0,a1,a2,a3,xi2
57 guez 3 INTEGER i,it,ik,iter,ii,idif,ii1,ii2
58 guez 71 DOUBLE PRECISION xi,xo1,xmoy,xlon2,fxm,Xprimin
59     DOUBLE PRECISION champmin,champmax,decalx
60 guez 3 INTEGER is2
61     SAVE is2
62    
63 guez 71 DOUBLE PRECISION heavyside
64 guez 3
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,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
81     ,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 + FLOAT(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,*) ' ** Attention ! La valeur beta calculee dans la rou
156     ,tine fxhyp est mauvaise ! '
157     WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ',
158     , ' et relancer ! *** '
159     STOP 1
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 + (FLOAT(i) + xuv - decalx) * depi / FLOAT(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/ ( FLOAT(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

  ViewVC Help
Powered by ViewVC 1.1.21