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

Annotation of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
File size: 11516 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 78 module fxhyp_m
2 guez 3
3 guez 78 IMPLICIT NONE
4 guez 3
5 guez 78 contains
6 guez 3
7 guez 78 SUBROUTINE fxhyp(xzoomdeg, grossism, dzooma, tau, rlonm025, xprimm025, &
8     rlonv, xprimv, rlonu, xprimu, rlonp025, xprimp025, champmin, champmax)
9 guez 3
10 guez 78 ! From LMDZ4/libf/dyn3d/fxhyp.F, v 1.2 2005/06/03 09:11:32 fairhead
11 guez 3
12 guez 78 ! Auteur : P. Le Van
13 guez 3
14 guez 78 ! Calcule les longitudes et dérivées dans la grille du GCM pour
15     ! une fonction f(x) à tangente hyperbolique.
16 guez 3
17 guez 78 ! On doit avoir grossism \times dzoom < pi (radians), en longitude.
18 guez 3
19 guez 78 USE dimens_m, ONLY: iim
20     USE paramet_m, ONLY: iip1
21 guez 3
22 guez 78 INTEGER nmax, nmax2
23     PARAMETER (nmax = 30000, nmax2 = 2*nmax)
24 guez 3
25 guez 78 LOGICAL scal180
26     PARAMETER (scal180 = .TRUE.)
27 guez 3
28 guez 78 ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour
29     ! une grille reguliere (grossism = 1., tau=0., clon=0.) a -180. degres.
30     ! sinon scal180 = .FALSE.
31 guez 3
32 guez 78 ! ...... arguments d'entree .......
33 guez 3
34 guez 78 REAL xzoomdeg, dzooma, tau, grossism
35     ! grossism etant le grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)
36     ! dzooma etant la distance totale de la zone du zoom
37     ! tau la raideur de la transition de l'interieur a l'exterieur du zoom
38 guez 3
39 guez 78 ! ...... arguments de sortie ......
40 guez 3
41 guez 78 REAL rlonm025(iip1), xprimm025(iip1), rlonv(iip1), xprimv(iip1), &
42     rlonu(iip1), xprimu(iip1), rlonp025(iip1), xprimp025(iip1)
43 guez 3
44 guez 78 ! .... variables locales ....
45 guez 3
46 guez 78 REAL dzoom
47     DOUBLE PRECISION xlon(iip1), xprimm(iip1), xuv
48     DOUBLE PRECISION xtild(0:nmax2)
49     DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)
50     DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)
51     DOUBLE PRECISION xvrai(iip1), xxprim(iip1)
52     DOUBLE PRECISION pi, depi, epsilon, xzoom, fa, fb
53     DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
54     INTEGER i, it, ik, iter, ii, idif, ii1, ii2
55     DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin
56     DOUBLE PRECISION champmin, champmax, decalx
57     INTEGER is2
58     SAVE is2
59 guez 3
60 guez 78 DOUBLE PRECISION heavyside
61 guez 3
62 guez 78 pi = 2. * ASIN(1.)
63     depi = 2. * pi
64     epsilon = 1.e-3
65     xzoom = xzoomdeg * pi/180.
66 guez 3
67 guez 78 decalx = .75
68     IF(grossism.EQ.1..AND.scal180) THEN
69     decalx = 1.
70     ENDIF
71 guez 3
72 guez 78 WRITE(6, *) 'FXHYP scal180, decalx', scal180, decalx
73 guez 3
74 guez 78 IF(dzooma.LT.1.) THEN
75     dzoom = dzooma * depi
76     ELSEIF(dzooma.LT. 25.) THEN
77     WRITE(6, *) ' Le param. dzoomx pour fxhyp est trop petit ! L augmenter et relancer ! '
78     STOP 1
79     ELSE
80     dzoom = dzooma * pi/180.
81     ENDIF
82 guez 3
83 guez 78 WRITE(6, *) ' xzoom(rad.), grossism, tau, dzoom (radians)'
84     WRITE(6, 24) xzoom, grossism, tau, dzoom
85 guez 3
86 guez 78 DO i = 0, nmax2
87     xtild(i) = - pi + FLOAT(i) * depi /nmax2
88     ENDDO
89    
90     DO i = nmax, nmax2
91    
92     fa = tau* (dzoom/2. - xtild(i))
93     fb = xtild(i) * (pi - xtild(i))
94    
95     IF(200.* fb .LT. - fa) THEN
96     fhyp (i) = - 1.
97     ELSEIF(200. * fb .LT. fa) THEN
98     fhyp (i) = 1.
99 guez 3 ELSE
100 guez 78 IF(ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
101     IF(200.*fb + fa.LT.1.e-10) THEN
102     fhyp (i) = - 1.
103     ELSEIF(200.*fb - fa.LT.1.e-10) THEN
104     fhyp (i) = 1.
105     ENDIF
106     ELSE
107     fhyp (i) = TANH (fa/fb)
108     ENDIF
109 guez 3 ENDIF
110 guez 78 IF (xtild(i).EQ. 0.) fhyp(i) = 1.
111     IF (xtild(i).EQ. pi) fhyp(i) = -1.
112 guez 3
113 guez 78 ENDDO
114 guez 3
115 guez 78 !c .... Calcul de beta ....
116 guez 3
117 guez 78 ffdx = 0.
118 guez 3
119 guez 78 DO i = nmax +1, nmax2
120 guez 3
121 guez 78 xmoy = 0.5 * (xtild(i-1) + xtild(i))
122     fa = tau* (dzoom/2. - xmoy)
123     fb = xmoy * (pi - xmoy)
124    
125     IF(200.* fb .LT. - fa) THEN
126     fxm = - 1.
127     ELSEIF(200. * fb .LT. fa) THEN
128     fxm = 1.
129     ELSE
130     IF(ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN
131     IF(200.*fb + fa.LT.1.e-10) THEN
132     fxm = - 1.
133     ELSEIF(200.*fb - fa.LT.1.e-10) THEN
134     fxm = 1.
135     ENDIF
136     ELSE
137     fxm = TANH (fa/fb)
138     ENDIF
139 guez 3 ENDIF
140    
141 guez 78 IF (xmoy.EQ. 0.) fxm = 1.
142     IF (xmoy.EQ. pi) fxm = -1.
143 guez 3
144 guez 78 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
145 guez 3
146 guez 78 ENDDO
147 guez 3
148 guez 78 beta = (grossism * ffdx - pi) / (ffdx - pi)
149 guez 3
150 guez 78 IF(2.*beta - grossism.LE. 0.) THEN
151     WRITE(6, *) ' ** Attention ! La valeur beta calculee dans la routine fxhyp est mauvaise ! '
152     WRITE(6, *)'Modifier les valeurs de grossismx, tau ou dzoomx ', &
153     ' et relancer ! *** '
154     STOP 1
155     ENDIF
156    
157     ! ..... calcul de Xprimt .....
158    
159    
160     DO i = nmax, nmax2
161     Xprimt(i) = beta + (grossism - beta) * fhyp(i)
162     ENDDO
163    
164     DO i = nmax+1, nmax2
165     Xprimt(nmax2 - i) = Xprimt(i)
166     ENDDO
167    
168    
169     ! ..... Calcul de Xf ........
170    
171     Xf(0) = - pi
172    
173     DO i = nmax +1, nmax2
174    
175     xmoy = 0.5 * (xtild(i-1) + xtild(i))
176     fa = tau* (dzoom/2. - xmoy)
177     fb = xmoy * (pi - xmoy)
178    
179     IF(200.* fb .LT. - fa) THEN
180     fxm = - 1.
181     ELSEIF(200. * fb .LT. fa) THEN
182     fxm = 1.
183 guez 3 ELSE
184 guez 78 fxm = TANH (fa/fb)
185 guez 3 ENDIF
186    
187 guez 78 IF (xmoy.EQ. 0.) fxm = 1.
188     IF (xmoy.EQ. pi) fxm = -1.
189     xxpr(i) = beta + (grossism - beta) * fxm
190 guez 3
191 guez 78 ENDDO
192 guez 3
193 guez 78 DO i = nmax+1, nmax2
194     xxpr(nmax2-i+1) = xxpr(i)
195     ENDDO
196 guez 3
197 guez 78 DO i=1, nmax2
198     Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
199     ENDDO
200 guez 3
201 guez 78 ! *****************************************************************
202 guez 3
203    
204 guez 78 ! ..... xuv = 0. si calcul aux pts scalaires ........
205     ! ..... xuv = 0.5 si calcul aux pts U ........
206 guez 3
207 guez 78 WRITE(6, 18)
208    
209     DO ik = 1, 4
210    
211     IF(ik.EQ.1) THEN
212     xuv = -0.25
213     ELSE IF (ik.EQ.2) THEN
214     xuv = 0.
215     ELSE IF (ik.EQ.3) THEN
216     xuv = 0.50
217     ELSE IF (ik.EQ.4) THEN
218     xuv = 0.25
219 guez 3 ENDIF
220    
221 guez 78 xo1 = 0.
222 guez 3
223 guez 78 ii1=1
224     ii2=iim
225     IF(ik.EQ.1.and.grossism.EQ.1.) THEN
226     ii1 = 2
227     ii2 = iim+1
228     ENDIF
229     DO i = ii1, ii2
230 guez 3
231 guez 78 xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)
232 guez 3
233 guez 78 Xfi = xlon2
234 guez 3
235 guez 78 DO it = nmax2, 0, -1
236     IF(Xfi.GE.Xf(it)) GO TO 350
237     end DO
238 guez 3
239 guez 78 it = 0
240 guez 3
241 guez 78 350 CONTINUE
242 guez 3
243 guez 78 ! ...... Calcul de Xf(xi) ......
244 guez 3
245 guez 78 xi = xtild(it)
246 guez 3
247 guez 78 IF(it.EQ.nmax2) THEN
248     it = nmax2 -1
249     Xf(it+1) = pi
250     ENDIF
251     ! .....................................................................
252 guez 3
253 guez 78 ! Appel de la routine qui calcule les coefficients a0, a1, a2, a3 d'un
254     ! polynome de degre 3 qui passe par les points (Xf(it), xtild(it))
255     ! et (Xf(it+1), xtild(it+1))
256 guez 3
257 guez 78 CALL coefpoly (Xf(it), Xf(it+1), Xprimt(it), Xprimt(it+1), &
258     xtild(it), xtild(it+1), a0, a1, a2, a3)
259    
260     Xf1 = Xf(it)
261     Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
262    
263     DO iter = 1, 300
264     xi = xi - (Xf1 - Xfi)/ Xprimin
265    
266     IF(ABS(xi-xo1).LE.epsilon) GO TO 550
267     xo1 = xi
268     xi2 = xi * xi
269     Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
270     Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
271     end DO
272     WRITE(6, *) ' Pas de solution ***** ', i, xlon2, iter
273 guez 3 STOP 6
274 guez 78 550 CONTINUE
275 guez 3
276 guez 78 xxprim(i) = depi/ (FLOAT(iim) * Xprimin)
277     xvrai(i) = xi + xzoom
278 guez 3
279 guez 78 end DO
280 guez 3
281 guez 78 IF(ik.EQ.1.and.grossism.EQ.1.) THEN
282     xvrai(1) = xvrai(iip1)-depi
283     xxprim(1) = xxprim(iip1)
284 guez 3 ENDIF
285 guez 78 DO i = 1, iim
286     xlon(i) = xvrai(i)
287     xprimm(i) = xxprim(i)
288 guez 3 ENDDO
289     DO i = 1, iim -1
290 guez 78 IF(xvrai(i+1).LT. xvrai(i)) THEN
291     WRITE(6, *) ' PBS. avec rlonu(', i+1, ') plus petit que rlonu(', i, &
292     ')'
293     STOP 7
294     ENDIF
295 guez 3 ENDDO
296    
297 guez 78 ! ... Reorganisation des longitudes pour les avoir entre - pi et pi ..
298     ! ........................................................................
299    
300     champmin = 1.e12
301 guez 3 champmax = -1.e12
302     DO i = 1, iim
303 guez 78 champmin = MIN(champmin, xvrai(i))
304     champmax = MAX(champmax, xvrai(i))
305 guez 3 ENDDO
306    
307 guez 78 IF(.not. (champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10)) THEN
308     WRITE(6, *) 'Reorganisation des longitudes pour avoir entre - pi', &
309     ' et pi '
310 guez 3
311 guez 78 IF(xzoom.LE.0.) THEN
312     IF(ik.EQ. 1) THEN
313     DO i = 1, iim
314     IF(xvrai(i).GE. - pi) GO TO 80
315     ENDDO
316     WRITE(6, *) ' PBS. 1 ! Xvrai plus petit que - pi ! '
317     STOP 8
318     80 CONTINUE
319     is2 = i
320     ENDIF
321    
322     IF(is2.NE. 1) THEN
323     DO ii = is2, iim
324     xlon (ii-is2+1) = xvrai(ii)
325     xprimm(ii-is2+1) = xxprim(ii)
326     ENDDO
327     DO ii = 1, is2 -1
328     xlon (ii+iim-is2+1) = xvrai(ii) + depi
329     xprimm(ii+iim-is2+1) = xxprim(ii)
330     ENDDO
331     ENDIF
332     ELSE
333     IF(ik.EQ.1) THEN
334     DO i = iim, 1, -1
335     IF(xvrai(i).LE. pi) GO TO 90
336     ENDDO
337     WRITE(6, *) ' PBS. 2 ! Xvrai plus grand que pi ! '
338     STOP 9
339     90 CONTINUE
340     is2 = i
341     ENDIF
342     idif = iim -is2
343     DO ii = 1, is2
344     xlon (ii+idif) = xvrai(ii)
345     xprimm(ii+idif) = xxprim(ii)
346     ENDDO
347     DO ii = 1, idif
348     xlon (ii) = xvrai (ii+is2) - depi
349     xprimm(ii) = xxprim(ii+is2)
350     ENDDO
351 guez 3 ENDIF
352 guez 78 ENDIF
353 guez 3
354 guez 78 ! ......... Fin de la reorganisation ............................
355 guez 3
356 guez 78 xlon (iip1) = xlon(1) + depi
357     xprimm(iip1) = xprimm (1)
358 guez 3
359 guez 78 DO i = 1, iim+1
360     xvrai(i) = xlon(i)*180./pi
361     ENDDO
362 guez 3
363 guez 78 IF(ik.EQ.1) THEN
364     ! WRITE(6, *) ' XLON aux pts. V-0.25 apres (en deg.) '
365     ! WRITE(6, 18)
366     ! WRITE(6, 68) xvrai
367     ! WRITE(6, *) ' XPRIM k ', ik
368     ! WRITE(6, 566) xprimm
369 guez 3
370 guez 78 DO i = 1, iim +1
371     rlonm025(i) = xlon(i)
372     xprimm025(i) = xprimm(i)
373     ENDDO
374     ELSE IF(ik.EQ.2) THEN
375     ! WRITE(6, 18)
376     ! WRITE(6, *) ' XLON aux pts. V apres (en deg.) '
377     ! WRITE(6, 68) xvrai
378     ! WRITE(6, *) ' XPRIM k ', ik
379     ! WRITE(6, 566) xprimm
380 guez 3
381 guez 78 DO i = 1, iim + 1
382     rlonv(i) = xlon(i)
383     xprimv(i) = xprimm(i)
384     ENDDO
385 guez 3
386 guez 78 ELSE IF(ik.EQ.3) THEN
387     ! WRITE(6, 18)
388     ! WRITE(6, *) ' XLON aux pts. U apres (en deg.) '
389     ! WRITE(6, 68) xvrai
390     ! WRITE(6, *) ' XPRIM ik ', ik
391     ! WRITE(6, 566) xprimm
392 guez 3
393 guez 78 DO i = 1, iim + 1
394     rlonu(i) = xlon(i)
395     xprimu(i) = xprimm(i)
396     ENDDO
397 guez 3
398 guez 78 ELSE IF(ik.EQ.4) THEN
399     ! WRITE(6, 18)
400     ! WRITE(6, *) ' XLON aux pts. V+0.25 apres (en deg.) '
401     ! WRITE(6, 68) xvrai
402     ! WRITE(6, *) ' XPRIM ik ', ik
403     ! WRITE(6, 566) xprimm
404 guez 3
405 guez 78 DO i = 1, iim + 1
406     rlonp025(i) = xlon(i)
407     xprimp025(i) = xprimm(i)
408     ENDDO
409 guez 3
410 guez 78 ENDIF
411 guez 3
412 guez 78 end DO
413 guez 3
414 guez 78 WRITE(6, 18)
415 guez 3
416 guez 78 DO i = 1, iim
417     xlon(i) = rlonv(i+1) - rlonv(i)
418     ENDDO
419     champmin = 1.e12
420     champmax = -1.e12
421     DO i = 1, iim
422     champmin = MIN(champmin, xlon(i))
423     champmax = MAX(champmax, xlon(i))
424     ENDDO
425     champmin = champmin * 180./pi
426     champmax = champmax * 180./pi
427 guez 3
428 guez 78 18 FORMAT(/)
429     24 FORMAT(2x, 'Parametres xzoom, gross, tau, dzoom pour fxhyp ', 4f8.3)
430     68 FORMAT(1x, 7f9.2)
431     566 FORMAT(1x, 7f9.4)
432    
433     END SUBROUTINE fxhyp
434    
435     end module fxhyp_m

  ViewVC Help
Powered by ViewVC 1.1.21