/[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 122 by guez, Tue Feb 3 19:30:48 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 (ik == 1 .and. (MINval(xvrai(:iim)) < - pi_d - 0.1d0 &      xlon = xvrai + clon
66           .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  
67    
68            is2 = i    end subroutine invert_zoom_x
69         END IF  
70      END IF    !**********************************************************************
71    
72      SUBROUTINE funcd(x, fval, fderiv)
73    
74        use coefpoly_m, only: a0, a1, a2, a3
75    
76      if (is2 /= 0) then      DOUBLE PRECISION, INTENT(IN):: x
77         IF (xzoom <= 0.) THEN      DOUBLE PRECISION, INTENT(OUT):: fval, fderiv
           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  
           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.122  
changed lines
  Added in v.154
  ViewVC Help
Powered by ViewVC 1.1.21