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

Diff of /trunk/Sources/dyn3d/invert_zoom_x.f

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

revision 122 by guez, Tue Feb 3 19:30:48 2015 UTC revision 123 by guez, Thu Feb 5 12:41:08 2015 UTC
# Line 26  contains Line 26  contains
26    
27      ! Local:      ! Local:
28    
29      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin, xi2      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin
30      integer ii1, ii2, i, it, iter, idif, ii      integer ii1, ii2, i, it, iter, idif, ii
31      DOUBLE PRECISION, parameter:: my_eps = 1e-3      DOUBLE PRECISION, parameter:: my_eps = 1e-6
32      DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)      DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)
33      INTEGER:: is2 = 0      INTEGER:: is2 = 0
34    
35      !------------------------------------------------------------------      !------------------------------------------------------------------
36    
     xo1 = 0.  
   
37      IF (ik == 1 .and. grossismx == 1.) THEN      IF (ik == 1 .and. grossismx == 1.) THEN
38         ii1 = 2         ii1 = 2
39         ii2 = iim + 1         ii2 = iim + 1
# Line 60  contains Line 58  contains
58            it = 2 * nmax -1            it = 2 * nmax -1
59         END IF         END IF
60    
        ! Appel de la routine qui calcule les coefficients a0, a1,  
        ! a2, a3 d'un polynome de degre 3 qui passe par les points  
        ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))  
   
61         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
62              xtild(it), xtild(it + 1), a0, a1, a2, a3)              xtild(it), xtild(it + 1), a0, a1, a2, a3)
   
63         Xf1 = Xf(it)         Xf1 = Xf(it)
64         Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi**2         Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)
65           xo1 = xi
66         iter = 1         iter = 1
67    
68         do         do
69            xi = xi - (Xf1 - Xfi) / Xprimin            xi = xi - (Xf1 - Xfi) / Xprimin
70            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
71            xo1 = xi            xo1 = xi
72            xi2 = xi * xi            Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))
73            Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi            Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)
           Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2  
74         end DO         end DO
75    
76         if (ABS(xi - xo1) > my_eps) then         if (ABS(xi - xo1) > my_eps) then
Legend:
Removed from v.122  
changed lines
  Added in v.123
  ViewVC Help
Powered by ViewVC 1.1.21