/[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 124 by guez, Thu Feb 5 15:19:37 2015 UTC revision 126 by guez, Fri Feb 6 18:33:15 2015 UTC
# Line 6  module fxhyp_loop_ik_m Line 6  module fxhyp_loop_ik_m
6    
7  contains  contains
8    
9    subroutine fxhyp_loop_ik(ik, decalx, xf, xtild, Xprimt, xzoom, xlon, &    subroutine fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, xlon, xprimm, xuv)
        xprimm, xuv)  
10    
11      use coefpoly_m, only: coefpoly      use coefpoly_m, only: coefpoly
12      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
13      use nr_util, only: pi_d, twopi_d      use nr_util, only: pi_d, twopi_d
     use serre, only: grossismx  
14    
     INTEGER, intent(in):: ik  
     DOUBLE PRECISION, intent(in):: decalx  
15      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)
16      DOUBLE PRECISION, intent(in):: xzoom      DOUBLE PRECISION, intent(in):: xzoom
17      real, intent(out):: xlon(:), xprimm(:) ! (iim)      real, intent(out):: xlon(:), xprimm(:) ! (iim)
# Line 25  contains Line 21  contains
21      ! 0.5 si calcul aux points U      ! 0.5 si calcul aux points U
22    
23      ! Local:      ! Local:
24      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin      DOUBLE PRECISION xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin
25      integer ii1, ii2, i, it, iter      integer i, it, iter
26      DOUBLE PRECISION, parameter:: my_eps = 1e-6      DOUBLE PRECISION, parameter:: my_eps = 1d-6
     DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)  
27    
28      !------------------------------------------------------------------      DOUBLE PRECISION xxprim(iim), xvrai(iim)
29        ! intermediary variables because xlon and xprimm are simple precision
30    
31      IF (ik == 1 .and. grossismx == 1.) THEN      !------------------------------------------------------------------
        ii1 = 2  
        ii2 = iim + 1  
     else  
        ii1=1  
        ii2=iim  
     END IF  
32    
33      DO i = ii1, ii2      DO i = 1, iim
34         Xfi = - pi_d + (i + xuv - decalx) * twopi_d / iim         Xfi = - pi_d + (i + xuv - 0.75d0) * twopi_d / iim
35    
36         it = 2 * nmax         it = 2 * nmax
37         do while (xfi < xf(it) .and. it >= 1)         do while (xfi < xf(it) .and. it >= 1)
38            it = it - 1            it = it - 1
39         end do         end do
40    
41         ! Calcul de Xf(xi)         ! Calcul de Xf(xvrai(i))
42    
43         xi = xtild(it)         xvrai(i) = xtild(it)
44    
45         IF (it == 2 * nmax) THEN         IF (it == 2 * nmax) THEN
46            it = 2 * nmax -1            it = 2 * nmax -1
# Line 59  contains Line 49  contains
49         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
50              xtild(it), xtild(it + 1), a0, a1, a2, a3)              xtild(it), xtild(it + 1), a0, a1, a2, a3)
51         Xf1 = Xf(it)         Xf1 = Xf(it)
52         Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)         Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
53         xo1 = xi         xo1 = xvrai(i)
54         iter = 1         iter = 1
55    
56         do         do
57            xi = xi - (Xf1 - Xfi) / Xprimin            xvrai(i) = xvrai(i) - (Xf1 - Xfi) / Xprimin
58            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit            IF (ABS(xvrai(i) - xo1) <= my_eps .or. iter == 300) exit
59            xo1 = xi            xo1 = xvrai(i)
60            Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))            Xf1 = a0 + xvrai(i) * (a1 + xvrai(i) * (a2 + xvrai(i) * a3))
61            Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)            Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
62         end DO         end DO
63    
64         if (ABS(xi - xo1) > my_eps) then         if (ABS(xvrai(i) - xo1) > my_eps) then
65            ! iter == 300            ! iter == 300
66            print *, 'Pas de solution.'            print *, 'Pas de solution.'
67            print *, i, xfi            print *, i, xfi
68            STOP 1            STOP 1
69         end if         end if
70    
71         xxprim(i) = twopi_d / (REAL(iim) * Xprimin)         xxprim(i) = twopi_d / (iim * Xprimin)
        xvrai(i) = xi + xzoom  
72      end DO      end DO
73    
     IF (ik == 1 .and. grossismx == 1.) THEN  
        xvrai(1) = xvrai(iim + 1)-twopi_d  
        xxprim(1) = xxprim(iim + 1)  
     END IF  
   
74      DO i = 1, iim -1      DO i = 1, iim -1
75         IF (xvrai(i + 1) < xvrai(i)) THEN         IF (xvrai(i + 1) < xvrai(i)) THEN
76            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
# Line 94  contains Line 78  contains
78         END IF         END IF
79      END DO      END DO
80    
81      DO i = 1, iim      xlon = xvrai + xzoom
82         xlon(i) = xvrai(i)      xprimm = xxprim
        xprimm(i) = xxprim(i)  
     END DO  
83    
84    end subroutine fxhyp_loop_ik    end subroutine fxhyp_loop_ik
85    
Legend:
Removed from v.124  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21