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

Annotation of /trunk/Sources/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 91 - (hide annotations)
Wed Mar 26 17:18:58 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/dyn3d/fxhyp.f
File size: 10520 byte(s) 
Removed unused variables lock_startdate and time_stamp of module
calendar.

Noticed that physiq does not change the surface pressure. So removed
arguments ps and dpfi of subroutine addfi. dpfi was always 0. The
computation of ps in addfi included some averaging at the poles. In
principle, this does not change ps but in practice it does because of
finite numerical precision. So the results of the simulation are
changed. Removed arguments ps and dpfi of calfis. Removed argument
d_ps of physiq.

du at the poles is not computed by dudv1, so declare only the
corresponding latitudes in dudv1. caldyn passes only a section of the
array dudyn as argument.

Removed variable niadv of module iniadvtrac_m.

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