trunk/dyn3d/invert_zoom_x.f

module invert_zoom_x_m

  implicit none

  DOUBLE PRECISION, private, save:: abs_y
  private funcd

contains

  subroutine invert_zoom_x(beta, xf, xtild, G, xlon, xprim, xuv)

    ! Libraries:
    use nr_util, only: pi_d, twopi_d
    use numer_rec_95, only: hunt, rtsafe    
    
    use coefpoly_m, only: coefpoly, a1, a2, a3
    use dynetat0_chosen_m, only: clon, grossismx
    USE dimensions, ONLY: iim

    DOUBLE PRECISION, intent(in):: beta, Xf(0:), xtild(0:), G(0:) ! (0:nmax)

    real, intent(out):: xlon(:), xprim(:) ! (iim)

    DOUBLE PRECISION, intent(in):: xuv
    ! between - 0.25 and 0.5
    ! 0. si calcul aux points scalaires 
    ! 0.5 si calcul aux points U

    ! Local:
    
    DOUBLE PRECISION Y
    DOUBLE PRECISION h ! step of the uniform grid
    integer i, it

    DOUBLE PRECISION xvrai(iim), Gvrai(iim) 
    ! intermediary variables because xlon and xprim are single precision

    !------------------------------------------------------------------

    print *, "Call sequence information: invert_zoom_x"
    it = 0 ! initial guess
    h = twopi_d / iim

    DO i = 1, iim
       Y = - pi_d + (i + xuv - 0.75d0) * h
       ! - pi <= y < pi
       abs_y = abs(y)

       ! Distinguish boundaries in order to avoid roundoff error.
       ! funcd should be exactly equal to 0 at xtild(it) or xtild(it +
       ! 1) and could be very small with the wrong sign so rtsafe
       ! would fail.
       if (abs_y == 0d0) then
          xvrai(i) = 0d0
          gvrai(i) = grossismx
       else if (abs_y == pi_d) then
          xvrai(i) = pi_d
          gvrai(i) = 2d0 * beta - grossismx
       else
          call hunt(xf, abs_y, it, my_lbound = 0)
          ! {0 <= it <= nmax - 1}

          ! Calcul de xvrai(i) et Gvrai(i)
          CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
               xtild(it + 1))
          xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
          Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
       end if

       if (y < 0d0) xvrai(i) = - xvrai(i)
    end DO

    DO i = 1, iim - 1
       IF (xvrai(i + 1) < xvrai(i)) THEN
          print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
          STOP 1
       END IF
    END DO

    xlon = xvrai + clon
    xprim = h / Gvrai

  end subroutine invert_zoom_x

  !**********************************************************************

  SUBROUTINE funcd(x, fval, fderiv)

    use coefpoly_m, only: a0, a1, a2, a3

    DOUBLE PRECISION, INTENT(IN):: x
    DOUBLE PRECISION, INTENT(OUT):: fval, fderiv

    fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y
    fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3)

  END SUBROUTINE funcd

end module invert_zoom_x_m
1	guez	314	module invert_zoom_x_m
2
3			implicit none
4
5			DOUBLE PRECISION, private, save:: abs_y
6			private funcd
7
8			contains
9
10			subroutine invert_zoom_x(beta, xf, xtild, G, xlon, xprim, xuv)
11
12			! Libraries:
13			use nr_util, only: pi_d, twopi_d
14			use numer_rec_95, only: hunt, rtsafe
15
16			use coefpoly_m, only: coefpoly, a1, a2, a3
17			use dynetat0_chosen_m, only: clon, grossismx
18			USE dimensions, ONLY: iim
19
20			DOUBLE PRECISION, intent(in):: beta, Xf(0:), xtild(0:), G(0:) ! (0:nmax)
21
22			real, intent(out):: xlon(:), xprim(:) ! (iim)
23
24			DOUBLE PRECISION, intent(in):: xuv
25			! between - 0.25 and 0.5
26			! 0. si calcul aux points scalaires
27			! 0.5 si calcul aux points U
28
29			! Local:
30
31			DOUBLE PRECISION Y
32			DOUBLE PRECISION h ! step of the uniform grid
33			integer i, it
34
35			DOUBLE PRECISION xvrai(iim), Gvrai(iim)
36			! intermediary variables because xlon and xprim are single precision
37
38			!------------------------------------------------------------------
39
40			print *, "Call sequence information: invert_zoom_x"
41			it = 0 ! initial guess
42			h = twopi_d / iim
43
44			DO i = 1, iim
45			Y = - pi_d + (i + xuv - 0.75d0) * h
46			! - pi <= y < pi
47			abs_y = abs(y)
48
49			! Distinguish boundaries in order to avoid roundoff error.
50			! funcd should be exactly equal to 0 at xtild(it) or xtild(it +
51			! 1) and could be very small with the wrong sign so rtsafe
52			! would fail.
53			if (abs_y == 0d0) then
54			xvrai(i) = 0d0
55			gvrai(i) = grossismx
56			else if (abs_y == pi_d) then
57			xvrai(i) = pi_d
58			gvrai(i) = 2d0 * beta - grossismx
59			else
60			call hunt(xf, abs_y, it, my_lbound = 0)
61			! {0 <= it <= nmax - 1}
62
63			! Calcul de xvrai(i) et Gvrai(i)
64			CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
65			xtild(it + 1))
66			xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
67			Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
68			end if
69
70			if (y < 0d0) xvrai(i) = - xvrai(i)
71			end DO
72
73			DO i = 1, iim - 1
74			IF (xvrai(i + 1) < xvrai(i)) THEN
75			print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
76			STOP 1
77			END IF
78			END DO
79
80			xlon = xvrai + clon
81			xprim = h / Gvrai
82
83			end subroutine invert_zoom_x
84
85			!**********************************************************************
86
87			SUBROUTINE funcd(x, fval, fderiv)
88
89			use coefpoly_m, only: a0, a1, a2, a3
90
91			DOUBLE PRECISION, INTENT(IN):: x
92			DOUBLE PRECISION, INTENT(OUT):: fval, fderiv
93
94			fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y
95			fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3)
96
97			END SUBROUTINE funcd
98
99			end module invert_zoom_x_m