/[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/Sources/dyn3d/invert_zoom_x.f revision 154 by guez, Tue Jul 7 17:49:23 2015 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(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
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):: 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 simple 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.      it = 0 ! initial guess
40        h = twopi_d / iim
     IF (ik == 1 .and. grossismx == 1.) THEN  
        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)  
   
        IF (it == 2 * nmax) THEN  
           it = 2 * nmax -1  
        END IF  
   
        ! Appel de la routine qui calcule les coefficients a0, a1,  
        ! a2, a3 d'un polynome de degre 3 qui passe par les points  
        ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))  
   
        CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &  
             xtild(it), xtild(it + 1), a0, a1, a2, a3)  
   
        Xf1 = Xf(it)  
        Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi**2  
   
        iter = 1  
   
        do  
           xi = xi - (Xf1 - Xfi) / Xprimin  
           IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit  
           xo1 = xi  
           xi2 = xi * xi  
           Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi  
           Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2  
        end DO  
   
        if (ABS(xi - xo1) > my_eps) then  
           ! iter == 300  
           print *, 'Pas de solution.'  
           print *, i, xfi  
           STOP 1  
        end if  
   
        xxprim(i) = twopi_d / (REAL(iim) * Xprimin)  
        xvrai(i) = xi + xzoom  
     end DO  
   
     IF (ik == 1 .and. grossismx == 1.) THEN  
        xvrai(1) = xvrai(iim + 1)-twopi_d  
        xxprim(1) = xxprim(iim + 1)  
     END IF  
41    
42      DO i = 1, iim      DO i = 1, iim
43         xlon(i) = xvrai(i)         Y = - pi_d + (i + xuv - 0.75d0) * h
44         xprimm(i) = xxprim(i)         ! - pi <= y < pi
45      END DO         abs_y = abs(y)
46    
47           call hunt(xf, abs_y, it, my_lbound = 0)
48           ! {0 <= it <= nmax - 1}
49    
50           ! Calcul de xvrai(i) et Gvrai(i)
51           CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
52                xtild(it + 1))
53           xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
54           Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
55           if (y < 0d0) xvrai(i) = - xvrai(i)
56        end DO
57    
58      DO i = 1, iim -1      DO i = 1, iim -1
59         IF (xvrai(i + 1) < xvrai(i)) THEN         IF (xvrai(i + 1) < xvrai(i)) THEN
60            print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'            print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
61            STOP 1            STOP 1
62         END IF         END IF
63      END DO      END DO
64    
65      IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 1d-5 &      xlon = xvrai + clon
66           .and. MAXval(xvrai(:iim)) <= pi_d + 1d-5)) THEN      xprim = h / Gvrai
67         IF (xzoom <= 0.) THEN  
68            IF (ik == 1) THEN    end subroutine invert_zoom_x
69               i = 1  
70      !**********************************************************************
71               do while (xvrai(i) < - pi_d .and. i < iim)  
72                  i = i + 1    SUBROUTINE funcd(x, fval, fderiv)
73               end do  
74        use coefpoly_m, only: a0, a1, a2, a3
75               if (xvrai(i) < - pi_d) then  
76                  print *, 'Xvrai plus petit que - pi !'      DOUBLE PRECISION, INTENT(IN):: x
77                  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  
78    
79      xlon(iim + 1) = xlon(1) + twopi      fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y
80      xprimm(iim + 1) = xprimm(1)      fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3)
81    
82    end subroutine fxhyp_loop_ik    END SUBROUTINE funcd
83    
84  end module fxhyp_loop_ik_m  end module invert_zoom_x_m
Legend:
Removed from v.121  
changed lines
  Added in v.154
  ViewVC Help
Powered by ViewVC 1.1.21