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

Diff of /trunk/dyn3d/invert_zoom_x.f

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

trunk/dyn3d/fxhyp_loop_ik.f revision 124 by guez, Thu Feb 5 15:19:37 2015 UTC trunk/Sources/dyn3d/invert_zoom_x.f revision 144 by guez, Wed Jun 10 16:46:46 2015 UTC
# Line 1  Line 1 
1  module fxhyp_loop_ik_m  module invert_zoom_x_m
2    
3    implicit none    implicit none
4    
# 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 invert_zoom_x(xf, xtild, Xprimt, 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 dynetat0_m, only: clon
14      use nr_util, only: pi_d, twopi_d      use nr_util, only: pi_d, twopi_d
     use serre, only: grossismx  
15    
     INTEGER, intent(in):: ik  
     DOUBLE PRECISION, intent(in):: decalx  
16      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)
     DOUBLE PRECISION, intent(in):: xzoom  
17      real, intent(out):: xlon(:), xprimm(:) ! (iim)      real, intent(out):: xlon(:), xprimm(:) ! (iim)
18    
19      DOUBLE PRECISION, intent(in):: xuv      DOUBLE PRECISION, intent(in):: xuv
# 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))
   
        xi = xtild(it)  
   
        IF (it == 2 * nmax) THEN  
           it = 2 * nmax -1  
        END IF  
42    
43           xvrai(i) = xtild(it)
44           IF (it == 2 * nmax) it = 2 * nmax -1
45         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
46              xtild(it), xtild(it + 1), a0, a1, a2, a3)              xtild(it), xtild(it + 1), a0, a1, a2, a3)
47         Xf1 = Xf(it)         Xf1 = Xf(it)
48         Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)         Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
49         xo1 = xi         xo1 = xvrai(i)
50         iter = 1         iter = 1
51    
52         do         do
53            xi = xi - (Xf1 - Xfi) / Xprimin            xvrai(i) = xvrai(i) - (Xf1 - Xfi) / Xprimin
54            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit            IF (ABS(xvrai(i) - xo1) <= my_eps .or. iter == 300) exit
55            xo1 = xi            xo1 = xvrai(i)
56            Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))            Xf1 = a0 + xvrai(i) * (a1 + xvrai(i) * (a2 + xvrai(i) * a3))
57            Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)            Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
58         end DO         end DO
59    
60         if (ABS(xi - xo1) > my_eps) then         if (ABS(xvrai(i) - xo1) > my_eps) then
61            ! iter == 300            ! iter == 300
62            print *, 'Pas de solution.'            print *, 'Pas de solution.'
63            print *, i, xfi            print *, i, xfi
64            STOP 1            STOP 1
65         end if         end if
66    
67         xxprim(i) = twopi_d / (REAL(iim) * Xprimin)         xxprim(i) = twopi_d / (iim * Xprimin)
        xvrai(i) = xi + xzoom  
68      end DO      end DO
69    
     IF (ik == 1 .and. grossismx == 1.) THEN  
        xvrai(1) = xvrai(iim + 1)-twopi_d  
        xxprim(1) = xxprim(iim + 1)  
     END IF  
   
70      DO i = 1, iim -1      DO i = 1, iim -1
71         IF (xvrai(i + 1) < xvrai(i)) THEN         IF (xvrai(i + 1) < xvrai(i)) THEN
72            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
# Line 94  contains Line 74  contains
74         END IF         END IF
75      END DO      END DO
76    
77      DO i = 1, iim      xlon = xvrai + clon
78         xlon(i) = xvrai(i)      xprimm = xxprim
        xprimm(i) = xxprim(i)  
     END DO  
79    
80    end subroutine fxhyp_loop_ik    end subroutine invert_zoom_x
81    
82  end module fxhyp_loop_ik_m  end module invert_zoom_x_m
Legend:
Removed from v.124  
changed lines
  Added in v.144
  ViewVC Help
Powered by ViewVC 1.1.21