/[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 121 by guez, Wed Jan 28 16:10:02 2015 UTC trunk/dyn3d/invert_zoom_x.f revision 255 by guez, Tue Mar 6 13:39:57 2018 UTC
# Line 1  Line 1 
1  module fxhyp_loop_ik_m  module invert_zoom_x_m
2    
3    implicit none    implicit none
4    
5    INTEGER, PARAMETER:: nmax = 30000    INTEGER, PARAMETER:: nmax = 30000
6      DOUBLE PRECISION abs_y
7    
8      private abs_y, funcd
9    
10  contains  contains
11    
12    subroutine fxhyp_loop_ik(ik, decalx, xf, xtild, Xprimt, xzoom, xlon, &    subroutine invert_zoom_x(beta, xf, xtild, G, xlon, xprim, xuv)
        xprimm, xuv)  
13    
14      use coefpoly_m, only: coefpoly      use coefpoly_m, only: coefpoly, a1, a2, a3
15      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
16      use nr_util, only: pi_d, twopi_d, twopi      use dynetat0_m, only: clon, grossismx
17      use serre, only: grossismx      use nr_util, only: pi_d, twopi_d
18        use numer_rec_95, only: hunt, rtsafe
19    
20        DOUBLE PRECISION, intent(in):: beta, Xf(0:), xtild(0:), G(0:) ! (0:nmax)
21    
22      INTEGER, intent(in):: ik      real, intent(out):: xlon(:), xprim(:) ! (iim)
     DOUBLE PRECISION, intent(in):: decalx  
     DOUBLE PRECISION, intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax)  
     DOUBLE PRECISION, intent(in):: xzoom  
     real, intent(out):: xlon(:), xprimm(:) ! (iim + 1)  
23    
24      DOUBLE PRECISION, intent(in):: xuv      DOUBLE PRECISION, intent(in):: xuv
25        ! between - 0.25 and 0.5
26      ! 0. si calcul aux points scalaires      ! 0. si calcul aux points scalaires
27      ! 0.5 si calcul aux points U      ! 0.5 si calcul aux points U
28    
29      ! Local:      ! Local:
30        DOUBLE PRECISION Y
31        DOUBLE PRECISION h ! step of the uniform grid
32        integer i, it
33    
34      DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin, xi2      DOUBLE PRECISION xvrai(iim), Gvrai(iim)
35      integer ii1, ii2, i, it, iter, idif, ii      ! intermediary variables because xlon and xprim are single precision
     DOUBLE PRECISION, parameter:: my_eps = 1e-3  
     DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)  
     INTEGER:: is2 = 0  
36    
37      !------------------------------------------------------------------      !------------------------------------------------------------------
38    
39      xo1 = 0.      print *, "Call sequence information: invert_zoom_x"
40        it = 0 ! initial guess
41      IF (ik == 1 .and. grossismx == 1.) THEN      h = twopi_d / iim
        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  
   
        ! Calcul de Xf(xi)  
   
        xi = xtild(it)  
42    
43         IF (it == 2 * nmax) THEN      DO i = 1, iim
44            it = 2 * nmax -1         Y = - pi_d + (i + xuv - 0.75d0) * h
45         END IF         ! - pi <= y < pi
46           abs_y = abs(y)
47         ! Appel de la routine qui calcule les coefficients a0, a1,  
48         ! a2, a3 d'un polynome de degre 3 qui passe par les points         ! Distinguish boundaries in order to avoid roundoff error.
49         ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))         ! funcd should be exactly equal to 0 at xtild(it) or xtild(it +
50           ! 1) and could be very small with the wrong sign so rtsafe
51         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &         ! would fail.
52              xtild(it), xtild(it + 1), a0, a1, a2, a3)         if (abs_y == 0d0) then
53              xvrai(i) = 0d0
54         Xf1 = Xf(it)            gvrai(i) = grossismx
55         Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi**2         else if (abs_y == pi_d) then
56              xvrai(i) = pi_d
57         iter = 1            gvrai(i) = 2d0 * beta - grossismx
58           else
59         do            call hunt(xf, abs_y, it, my_lbound = 0)
60            xi = xi - (Xf1 - Xfi) / Xprimin            ! {0 <= it <= nmax - 1}
61            IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit  
62            xo1 = xi            ! Calcul de xvrai(i) et Gvrai(i)
63            xi2 = xi * xi            CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
64            Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi                 xtild(it + 1))
65            Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2            xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
66         end DO            Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
   
        if (ABS(xi - xo1) > my_eps) then  
           ! iter == 300  
           print *, 'Pas de solution.'  
           print *, i, xfi  
           STOP 1  
67         end if         end if
68    
69         xxprim(i) = twopi_d / (REAL(iim) * Xprimin)         if (y < 0d0) xvrai(i) = - xvrai(i)
        xvrai(i) = xi + xzoom  
70      end DO      end DO
71    
     IF (ik == 1 .and. grossismx == 1.) THEN  
        xvrai(1) = xvrai(iim + 1)-twopi_d  
        xxprim(1) = xxprim(iim + 1)  
     END IF  
   
     DO i = 1, iim  
        xlon(i) = xvrai(i)  
        xprimm(i) = xxprim(i)  
     END DO  
   
72      DO i = 1, iim -1      DO i = 1, iim -1
73         IF (xvrai(i + 1) < xvrai(i)) THEN         IF (xvrai(i + 1) < xvrai(i)) THEN
74            print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
75            STOP 1            STOP 1
76         END IF         END IF
77      END DO      END DO
78    
79      IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 1d-5 &      xlon = xvrai + clon
80           .and. MAXval(xvrai(:iim)) <= pi_d + 1d-5)) THEN      xprim = h / Gvrai
81         IF (xzoom <= 0.) THEN  
82            IF (ik == 1) THEN    end subroutine invert_zoom_x
83               i = 1  
84      !**********************************************************************
85               do while (xvrai(i) < - pi_d .and. i < iim)  
86                  i = i + 1    SUBROUTINE funcd(x, fval, fderiv)
87               end do  
88        use coefpoly_m, only: a0, a1, a2, a3
89               if (xvrai(i) < - pi_d) then  
90                  print *, 'Xvrai plus petit que - pi !'      DOUBLE PRECISION, INTENT(IN):: x
91                  STOP 1      DOUBLE PRECISION, INTENT(OUT):: fval, fderiv
              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  
92    
93      xlon(iim + 1) = xlon(1) + twopi      fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y
94      xprimm(iim + 1) = xprimm(1)      fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3)
95    
96    end subroutine fxhyp_loop_ik    END SUBROUTINE funcd
97    
98  end module fxhyp_loop_ik_m  end module invert_zoom_x_m
Legend:
Removed from v.121  
changed lines
  Added in v.255
  ViewVC Help
Powered by ViewVC 1.1.21