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

Diff of /trunk/dyn3d/fxhyp_loop_ik.f

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

revision 121 by guez, Wed Jan 28 16:10:02 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
# Line 109  contains Line 101  contains
101         END IF         END IF
102      END DO      END DO
103    
104      IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 1d-5 &      IF (ik == 1 .and. (MINval(xvrai(:iim)) < - pi_d - 0.1d0 &
105           .and. MAXval(xvrai(:iim)) <= pi_d + 1d-5)) THEN           .or. MAXval(xvrai(:iim)) > pi_d + 0.1d0)) THEN
106         IF (xzoom <= 0.) THEN         IF (xzoom <= 0.) THEN
107            IF (ik == 1) THEN            i = 1
              i = 1  
108    
109               do while (xvrai(i) < - pi_d .and. i < iim)            do while (xvrai(i) < - pi_d .and. i < iim)
110                  i = i + 1               i = i + 1
111               end do            end do
112    
113               if (xvrai(i) < - pi_d) then            if (xvrai(i) < - pi_d) then
114                  print *, 'Xvrai plus petit que - pi !'               print *, 'Xvrai plus petit que - pi !'
115                  STOP 1               STOP 1
116               end if            end if
117    
118               is2 = i            is2 = i
119            END IF         ELSE
120              i = iim
121    
122              do while (xvrai(i) > pi_d .and. i > 1)
123                 i = i - 1
124              end do
125    
126              if (xvrai(i) > pi_d) then
127                 print *, 'Xvrai plus grand que pi !'
128                 STOP 1
129              end if
130    
131              is2 = i
132           END IF
133        END IF
134    
135        if (is2 /= 0) then
136           IF (xzoom <= 0.) THEN
137            IF (is2 /= 1) THEN            IF (is2 /= 1) THEN
138               DO ii = is2, iim               DO ii = is2, iim
139                  xlon(ii-is2 + 1) = xvrai(ii)                  xlon(ii-is2 + 1) = xvrai(ii)
# Line 137  contains Line 144  contains
144                  xprimm(ii + iim-is2 + 1) = xxprim(ii)                  xprimm(ii + iim-is2 + 1) = xxprim(ii)
145               END DO               END DO
146            END IF            END IF
147         ELSE         else
           IF (ik == 1) THEN  
              i = iim  
   
              do while (xvrai(i) > pi_d .and. i > 1)  
                 i = i - 1  
              end do  
   
              if (xvrai(i) > pi_d) then  
                 print *, 'Xvrai plus grand que pi !'  
                 STOP 1  
              end if  
   
              is2 = i  
           END IF  
   
148            idif = iim -is2            idif = iim -is2
149    
150            DO ii = 1, is2            DO ii = 1, is2
# Line 164  contains Line 156  contains
156               xlon(ii) = xvrai(ii + is2) - twopi_d               xlon(ii) = xvrai(ii + is2) - twopi_d
157               xprimm(ii) = xxprim(ii + is2)               xprimm(ii) = xxprim(ii + is2)
158            END DO            END DO
159         END IF         end IF
160      END IF      end if
161    
162      xlon(iim + 1) = xlon(1) + twopi      xlon(iim + 1) = xlon(1) + twopi
163      xprimm(iim + 1) = xprimm(1)      xprimm(iim + 1) = xprimm(1)
Legend:
Removed from v.121  
changed lines
  Added in v.123
  ViewVC Help
Powered by ViewVC 1.1.21