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

Diff of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 119 by guez, Wed Jan 7 14:34:57 2015 UTC revision 120 by guez, Tue Jan 13 14:56:15 2015 UTC
# Line 14  contains Line 14  contains
14    
15      ! On doit avoir grossismx \times dzoomx < pi (radians)      ! On doit avoir grossismx \times dzoomx < pi (radians)
16    
17        ! Le premier point scalaire pour une grille regulière (grossismx =
18        ! 1., taux=0., clon=0.) est à - 180 degrés.
19    
20        use coefpoly_m, only: coefpoly
21      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
22      use nr_util, only: pi_d, twopi_d      use nr_util, only: pi_d, twopi_d, arth
23      use serre, only: clon, grossismx, dzoomx, taux      use serre, only: clon, grossismx, dzoomx, taux
24    
25      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
# Line 25  contains Line 29  contains
29    
30      DOUBLE PRECISION champmin, champmax      DOUBLE PRECISION champmin, champmax
31      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
32      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax
   
     LOGICAL, PARAMETER:: scal180 = .TRUE.  
     ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour  
     ! une grille reguliere (grossismx = 1., taux=0., clon=0.) a  
     ! -180. degres. sinon scal180 = .FALSE.  
   
33      REAL dzoom      REAL dzoom
34      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv
35      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION xtild(0:nmax2)
36      DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)
37      DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)      DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)
38      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)
39      DOUBLE PRECISION my_eps, xzoom, fa, fb      DOUBLE PRECISION my_eps, xzoom, fa, fb
40      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
41      INTEGER i, it, ik, iter, ii, idif, ii1, ii2      INTEGER i, it, ik, iter, ii, idif, ii1, ii2
42      DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin      DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin
43      DOUBLE PRECISION decalx      DOUBLE PRECISION decalx
44      INTEGER, save:: is2      INTEGER is2
45    
46      !----------------------------------------------------------------------      !----------------------------------------------------------------------
47    
48        print *, "Call sequence information: fxhyp"
49    
50      my_eps = 1e-3      my_eps = 1e-3
51      xzoom = clon * pi_d / 180.      xzoom = clon * pi_d / 180.
52    
53      IF (grossismx == 1. .AND. scal180) THEN      IF (grossismx == 1.) THEN
54         decalx = 1.         decalx = 1.
55      else      else
56         decalx = 0.75         decalx = 0.75
# Line 59  contains Line 59  contains
59      IF (dzoomx < 1.) THEN      IF (dzoomx < 1.) THEN
60         dzoom = dzoomx * twopi_d         dzoom = dzoomx * twopi_d
61      ELSE IF (dzoomx < 25.) THEN      ELSE IF (dzoomx < 25.) THEN
62         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &         print *, "dzoomx pour fxhyp est trop petit."
             // "L'augmenter et relancer."  
63         STOP 1         STOP 1
64      ELSE      ELSE
65         dzoom = dzoomx * pi_d / 180.         dzoom = dzoomx * pi_d / 180.
# Line 68  contains Line 67  contains
67    
68      print *, 'dzoom (rad):', dzoom      print *, 'dzoom (rad):', dzoom
69    
70      DO i = 0, nmax2      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)
        xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2  
     END DO  
71    
72      DO i = nmax, nmax2      DO i = nmax, nmax2
73         fa = taux* (dzoom / 2. - xtild(i))         fa = taux * (dzoom / 2. - xtild(i))
74         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
75    
76         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
77            fhyp(i) = - 1.            fhyp(i) = - 1.
78         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
79            fhyp(i) = 1.            fhyp(i) = 1.
80         ELSE         ELSE
81            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
82               IF (200.*fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
83                  fhyp(i) = - 1.                  fhyp(i) = - 1.
84               ELSE IF (200.*fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
85                  fhyp(i) = 1.                  fhyp(i) = 1.
86               END IF               END IF
87            ELSE            ELSE
# Line 102  contains Line 99  contains
99    
100      DO i = nmax + 1, nmax2      DO i = nmax + 1, nmax2
101         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
102         fa = taux* (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
103         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
104    
105         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
106            fxm = - 1.            fxm = - 1.
107         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
108            fxm = 1.            fxm = 1.
109         ELSE         ELSE
110            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
111               IF (200.*fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
112                  fxm = - 1.                  fxm = - 1.
113               ELSE IF (200.*fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
114                  fxm = 1.                  fxm = 1.
115               END IF               END IF
116            ELSE            ELSE
# Line 151  contains Line 148  contains
148    
149      DO i = nmax + 1, nmax2      DO i = nmax + 1, nmax2
150         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
151         fa = taux* (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
152         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
153    
154         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
155            fxm = - 1.            fxm = - 1.
156         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
157            fxm = 1.            fxm = 1.
# Line 167  contains Line 164  contains
164         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
165      END DO      END DO
166    
167      DO i = nmax + 1, nmax2      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)
        xxpr(nmax2-i + 1) = xxpr(i)  
     END DO  
168    
169      DO i=1, nmax2      DO i=1, nmax2
170         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
171      END DO      END DO
172    
173      ! xuv = 0. si calcul aux points scalaires      is2 = 0
     ! xuv = 0.5 si calcul aux points U  
174    
175      loop_ik: DO ik = 1, 4      loop_ik: DO ik = 1, 4
176           ! xuv = 0. si calcul aux points scalaires
177           ! xuv = 0.5 si calcul aux points U
178    
179         IF (ik == 1) THEN         IF (ik == 1) THEN
180            xuv = -0.25            xuv = -0.25
181         ELSE IF (ik == 2) THEN         ELSE IF (ik == 2) THEN
# Line 191  contains Line 188  contains
188    
189         xo1 = 0.         xo1 = 0.
190    
191         ii1=1         IF (ik == 1 .and. grossismx == 1.) THEN
        ii2=iim  
        IF (ik == 1.and.grossismx == 1.) THEN  
192            ii1 = 2            ii1 = 2
193            ii2 = iim + 1            ii2 = iim + 1
194           else
195              ii1=1
196              ii2=iim
197         END IF         END IF
198    
199         DO i = ii1, ii2         DO i = ii1, ii2
200            xlon2 = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)            Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)
           Xfi = xlon2  
201    
202            it = nmax2            it = nmax2
203            do while (xfi < xf(it) .and. it >= 1)            do while (xfi < xf(it) .and. it >= 1)
# Line 224  contains Line 221  contains
221                 xtild(it), xtild(it + 1), a0, a1, a2, a3)                 xtild(it), xtild(it + 1), a0, a1, a2, a3)
222    
223            Xf1 = Xf(it)            Xf1 = Xf(it)
224            Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi            Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi * xi
225    
226            iter = 1            iter = 1
227    
# Line 234  contains Line 231  contains
231               xo1 = xi               xo1 = xi
232               xi2 = xi * xi               xi2 = xi * xi
233               Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi               Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
234               Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2               Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2
235            end DO            end DO
236    
237            if (ABS(xi - xo1) > my_eps) then            if (ABS(xi - xo1) > my_eps) then
238               ! iter == 300               ! iter == 300
239               print *, 'Pas de solution.'               print *, 'Pas de solution.'
240               print *, i, xlon2               print *, i, xfi
241               STOP 1               STOP 1
242            end if            end if
243    
# Line 260  contains Line 257  contains
257    
258         DO i = 1, iim -1         DO i = 1, iim -1
259            IF (xvrai(i + 1) < xvrai(i)) THEN            IF (xvrai(i + 1) < xvrai(i)) THEN
260               print *, 'Problème avec rlonu(', i + 1, &               print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'
                   ') plus petit que rlonu(', i, ')'  
261               STOP 1               STOP 1
262            END IF            END IF
263         END DO         END DO
264    
265         ! Réorganisation des longitudes pour les avoir entre - pi et pi         IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 0.1 &
266                .and. MAXval(xvrai(:iim)) <= pi_d + 0.1)) THEN
267         champmin = 1e12            print *, &
268         champmax = -1e12                 'Réorganisation des longitudes pour les avoir entre - pi et pi'
        DO i = 1, iim  
           champmin = MIN(champmin, xvrai(i))  
           champmax = MAX(champmax, xvrai(i))  
        END DO  
   
        IF (.not. (champmin >= -pi_d - 0.1 .and. champmax <= pi_d + 0.1)) THEN  
           print *, 'Reorganisation des longitudes pour avoir entre - pi', &  
                ' et pi '  
269    
270            IF (xzoom <= 0.) THEN            IF (xzoom <= 0.) THEN
271               IF (ik == 1) THEN               IF (ik == 1) THEN
# Line 295  contains Line 283  contains
283                  is2 = i                  is2 = i
284               END IF               END IF
285    
286               IF (is2.NE. 1) THEN               IF (is2 /= 1) THEN
287                  DO ii = is2, iim                  DO ii = is2, iim
288                     xlon(ii-is2 + 1) = xvrai(ii)                     xlon(ii-is2 + 1) = xvrai(ii)
289                     xprimm(ii-is2 + 1) = xxprim(ii)                     xprimm(ii-is2 + 1) = xxprim(ii)
# Line 335  contains Line 323  contains
323            END IF            END IF
324         END IF         END IF
325    
        ! Fin de la reorganisation  
   
326         xlon(iim + 1) = xlon(1) + twopi_d         xlon(iim + 1) = xlon(1) + twopi_d
327         xprimm(iim + 1) = xprimm(1)         xprimm(iim + 1) = xprimm(1)
328    
329         DO i = 1, iim + 1         DO i = 1, iim + 1
330            xvrai(i) = xlon(i)*180. / pi_d            xvrai(i) = xlon(i) * 180. / pi_d
331         END DO         END DO
332    
333         IF (ik == 1) THEN         IF (ik == 1) THEN
# Line 358  contains Line 344  contains
344               xprimu(i) = xprimm(i)               xprimu(i) = xprimm(i)
345            END DO            END DO
346         ELSE IF (ik == 4) THEN         ELSE IF (ik == 4) THEN
347            DO i = 1, iim + 1            rlonp025 = xlon
348               rlonp025(i) = xlon(i)            xprimp025 = xprimm
              xprimp025(i) = xprimm(i)  
           END DO  
349         END IF         END IF
350      end DO loop_ik      end DO loop_ik
351    
# Line 391  contains Line 375  contains
375         END IF         END IF
376    
377         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
378            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
379              print *, "> rlonu(", i, ") = ", rlonu(i)
380            STOP 1            STOP 1
381         END IF         END IF
382      END DO      END DO
# Line 402  contains Line 387  contains
387  3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &  3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
388           ' d environ ', f0.2, ' degres ', /, &           ' d environ ', f0.2, ' degres ', /, &
389           ' alors que la maille en dehors de la zone du zoom est ', &           ' alors que la maille en dehors de la zone du zoom est ', &
390           "d'environ", f0.2, ' degres ')           "d'environ ", f0.2, ' degres ')
391    
392    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
393    
Legend:
Removed from v.119  
changed lines
  Added in v.120
  ViewVC Help
Powered by ViewVC 1.1.21