/[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 123 by guez, Thu Feb 5 12:41:08 2015 UTC trunk/dyn3d/invert_zoom_x.f revision 254 by guez, Mon Feb 5 10:39:38 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      DOUBLE PRECISION xvrai(iim), Gvrai(iim)
35      integer ii1, ii2, i, it, iter, idif, ii      ! intermediary variables because xlon and xprim are simple precision
     DOUBLE PRECISION, parameter:: my_eps = 1e-6  
     DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)  
     INTEGER:: is2 = 0  
36    
37      !------------------------------------------------------------------      !------------------------------------------------------------------
38    
39      IF (ik == 1 .and. grossismx == 1.) THEN      print *, "Call sequence information: invert_zoom_x"
40         ii1 = 2      it = 0 ! initial guess
41         ii2 = iim + 1      h = twopi_d / iim
     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    
48           ! Distinguish boundaries in order to avoid roundoff error.
49           ! 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           ! would fail.
52           if (abs_y == 0d0) then
53              xvrai(i) = 0d0
54              gvrai(i) = grossismx
55           else if (abs_y == pi_d) then
56              xvrai(i) = pi_d
57              gvrai(i) = 2d0 * beta - grossismx
58           else
59              call hunt(xf, abs_y, it, my_lbound = 0)
60              ! {0 <= it <= nmax - 1}
61    
62         CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &            ! Calcul de xvrai(i) et Gvrai(i)
63              xtild(it), xtild(it + 1), a0, a1, a2, a3)            CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
64         Xf1 = Xf(it)                 xtild(it + 1))
65         Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)            xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
66         xo1 = xi            Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
        iter = 1  
   
        do  
           xi = xi - (Xf1 - Xfi) / Xprimin  
           IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit  
           xo1 = xi  
           Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))  
           Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)  
        end DO  
   
        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 (ik == 1 .and. (MINval(xvrai(:iim)) < - pi_d - 0.1d0 &      xlon = xvrai + clon
80           .or. MAXval(xvrai(:iim)) > pi_d + 0.1d0)) THEN      xprim = h / Gvrai
        IF (xzoom <= 0.) THEN  
           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  
        ELSE  
           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  
81    
82            is2 = i    end subroutine invert_zoom_x
        END IF  
     END IF  
83    
84      if (is2 /= 0) then    !**********************************************************************
85         IF (xzoom <= 0.) THEN  
86            IF (is2 /= 1) THEN    SUBROUTINE funcd(x, fval, fderiv)
87               DO ii = is2, iim  
88                  xlon(ii-is2 + 1) = xvrai(ii)      use coefpoly_m, only: a0, a1, a2, a3
                 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  
           idif = iim -is2  
89    
90            DO ii = 1, is2      DOUBLE PRECISION, INTENT(IN):: x
91               xlon(ii + idif) = xvrai(ii)      DOUBLE PRECISION, INTENT(OUT):: fval, fderiv
              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.123  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21