/[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

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 145 by guez, Tue Jun 16 15:23:29 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
15      use serre, only: grossismx      use numer_rec_95, only: hunt
16    
     INTEGER, intent(in):: ik  
     DOUBLE PRECISION, intent(in):: decalx  
17      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  
18      real, intent(out):: xlon(:), xprimm(:) ! (iim)      real, intent(out):: xlon(:), xprimm(:) ! (iim)
19    
20      DOUBLE PRECISION, intent(in):: xuv      DOUBLE PRECISION, intent(in):: xuv
21        ! between - 0.25 and 0.5
22      ! 0. si calcul aux points scalaires      ! 0. si calcul aux points scalaires
23      ! 0.5 si calcul aux points U      ! 0.5 si calcul aux points U
24    
25      ! Local:      ! Local:
26      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin      DOUBLE PRECISION xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin
27      integer ii1, ii2, i, it, iter      integer i, it, iter
28      DOUBLE PRECISION, parameter:: my_eps = 1e-6      DOUBLE PRECISION, parameter:: my_eps = 1d-6
29      DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)  
30        DOUBLE PRECISION xxprim(iim), xvrai(iim)
31        ! intermediary variables because xlon and xprimm are simple precision
32    
33      !------------------------------------------------------------------      !------------------------------------------------------------------
34    
35      IF (ik == 1 .and. grossismx == 1.) THEN      it = 0 ! initial guess
        ii1 = 2  
        ii2 = iim + 1  
     else  
        ii1=1  
        ii2=iim  
     END IF  
   
     DO i = ii1, ii2  
        Xfi = - pi_d + (i + xuv - decalx) * twopi_d / iim  
   
        it = 2 * nmax  
        do while (xfi < xf(it) .and. it >= 1)  
           it = it - 1  
        end do  
36    
37         ! Calcul de Xf(xi)      DO i = 1, iim
38           Xfi = - pi_d + (i + xuv - 0.75d0) * twopi_d / iim
39           ! - pi <= xfi < pi
40    
41         xi = xtild(it)         call hunt(xf, xfi, it, my_lbound = 0)
42           it = max(0, it)
43           ! In principle, xfi >= xf(0), max with 0 just in case of
44           ! roundoff error. {0 <= it <= 2 * nmax - 1}
45    
46         IF (it == 2 * nmax) THEN         ! Calcul de Xf(xvrai(i))
           it = 2 * nmax -1  
        END IF  
47    
48         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
49              xtild(it), xtild(it + 1), a0, a1, a2, a3)              xtild(it), xtild(it + 1), a0, a1, a2, a3)
50           xvrai(i) = xtild(it)
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 + clon
82         xlon(i) = xvrai(i)      xprimm = xxprim
        xprimm(i) = xxprim(i)  
     END DO  
83    
84    end subroutine fxhyp_loop_ik    end subroutine invert_zoom_x
85    
86  end module fxhyp_loop_ik_m  end module invert_zoom_x_m
Legend:
Removed from v.124  
changed lines
  Added in v.145
  ViewVC Help
Powered by ViewVC 1.1.21