/[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 124 by guez, Thu Feb 5 15:19:37 2015 UTC
# Line 11  contains Line 11  contains
11    
12      use coefpoly_m, only: coefpoly      use coefpoly_m, only: coefpoly
13      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
14      use nr_util, only: pi_d, twopi_d, twopi      use nr_util, only: pi_d, twopi_d
15      use serre, only: grossismx      use serre, only: grossismx
16    
17      INTEGER, intent(in):: ik      INTEGER, intent(in):: ik
18      DOUBLE PRECISION, intent(in):: decalx      DOUBLE PRECISION, intent(in):: decalx
19      DOUBLE PRECISION, intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax)      DOUBLE PRECISION, intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax)
20      DOUBLE PRECISION, intent(in):: xzoom      DOUBLE PRECISION, intent(in):: xzoom
21      real, intent(out):: xlon(:), xprimm(:) ! (iim + 1)      real, intent(out):: xlon(:), xprimm(:) ! (iim)
22    
23      DOUBLE PRECISION, intent(in):: xuv      DOUBLE PRECISION, intent(in):: xuv
24      ! 0. si calcul aux points scalaires      ! 0. si calcul aux points scalaires
25      ! 0.5 si calcul aux points U      ! 0.5 si calcul aux points U
26    
27      ! Local:      ! Local:
28        DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin
29      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin, xi2      integer ii1, ii2, i, it, iter
30      integer ii1, ii2, i, it, iter, idif, ii      DOUBLE PRECISION, parameter:: my_eps = 1e-6
     DOUBLE PRECISION, parameter:: my_eps = 1e-3  
31      DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)      DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)
     INTEGER:: is2 = 0  
32    
33      !------------------------------------------------------------------      !------------------------------------------------------------------
34    
     xo1 = 0.  
   
35      IF (ik == 1 .and. grossismx == 1.) THEN      IF (ik == 1 .and. grossismx == 1.) THEN
36         ii1 = 2         ii1 = 2
37         ii2 = iim + 1         ii2 = iim + 1
# Line 60  contains Line 56  contains
56            it = 2 * nmax -1            it = 2 * nmax -1
57         END IF         END IF
58    
        ! 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))  
   
59         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
60              xtild(it), xtild(it + 1), a0, a1, a2, a3)              xtild(it), xtild(it + 1), a0, a1, a2, a3)
   
61         Xf1 = Xf(it)         Xf1 = Xf(it)
62         Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi**2         Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)
63           xo1 = xi
64         iter = 1         iter = 1
65    
66         do         do
67            xi = xi - (Xf1 - Xfi) / Xprimin            xi = xi - (Xf1 - Xfi) / Xprimin
68            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
69            xo1 = xi            xo1 = xi
70            xi2 = xi * xi            Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))
71            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  
72         end DO         end DO
73    
74         if (ABS(xi - xo1) > my_eps) then         if (ABS(xi - xo1) > my_eps) then
# Line 97  contains Line 87  contains
87         xxprim(1) = xxprim(iim + 1)         xxprim(1) = xxprim(iim + 1)
88      END IF      END IF
89    
     DO i = 1, iim  
        xlon(i) = xvrai(i)  
        xprimm(i) = xxprim(i)  
     END DO  
   
90      DO i = 1, iim -1      DO i = 1, iim -1
91         IF (xvrai(i + 1) < xvrai(i)) THEN         IF (xvrai(i + 1) < xvrai(i)) THEN
92            print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
93            STOP 1            STOP 1
94         END IF         END IF
95      END DO      END DO
96    
97      IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 1d-5 &      DO i = 1, iim
98           .and. MAXval(xvrai(:iim)) <= pi_d + 1d-5)) THEN         xlon(i) = xvrai(i)
99         IF (xzoom <= 0.) THEN         xprimm(i) = xxprim(i)
100            IF (ik == 1) THEN      END DO
              i = 1  
   
              do while (xvrai(i) < - pi_d .and. i < iim)  
                 i = i + 1  
              end do  
   
              if (xvrai(i) < - pi_d) then  
                 print *, 'Xvrai plus petit que - pi !'  
                 STOP 1  
              end if  
   
              is2 = i  
           END IF  
   
           IF (is2 /= 1) THEN  
              DO ii = is2, iim  
                 xlon(ii-is2 + 1) = xvrai(ii)  
                 xprimm(ii-is2 + 1) = xxprim(ii)  
              END DO  
              DO ii = 1, is2 -1  
                 xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d  
                 xprimm(ii + iim-is2 + 1) = xxprim(ii)  
              END DO  
           END IF  
        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  
   
           idif = iim -is2  
   
           DO ii = 1, is2  
              xlon(ii + idif) = xvrai(ii)  
              xprimm(ii + idif) = xxprim(ii)  
           END DO  
   
           DO ii = 1, idif  
              xlon(ii) = xvrai(ii + is2) - twopi_d  
              xprimm(ii) = xxprim(ii + is2)  
           END DO  
        END IF  
     END IF  
   
     xlon(iim + 1) = xlon(1) + twopi  
     xprimm(iim + 1) = xprimm(1)  
101    
102    end subroutine fxhyp_loop_ik    end subroutine fxhyp_loop_ik
103    
Legend:
Removed from v.121  
changed lines
  Added in v.124
  ViewVC Help
Powered by ViewVC 1.1.21