trunk/dyn3d/invert_zoom_x.f

module invert_zoom_x_m

  implicit none

  INTEGER, PARAMETER:: nmax = 30000
  DOUBLE PRECISION abs_y

  private abs_y, funcd

contains

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

    use coefpoly_m, only: coefpoly, a1, a2, a3
    USE dimens_m, ONLY: iim
    use dynetat0_m, only: clon
    use nr_util, only: pi_d, twopi_d
    use numer_rec_95, only: hunt, rtsafe

    DOUBLE PRECISION, intent(in):: 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
    integer i, it

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

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

    it = 0 ! initial guess

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

       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)
       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 = twopi_d / (iim * 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	module invert_zoom_x_m
2
3	implicit none
4
5	INTEGER, PARAMETER:: nmax = 30000
6	DOUBLE PRECISION abs_y
7
8	private abs_y, funcd
9
10	contains
11
12	subroutine invert_zoom_x(xf, xtild, G, xlon, xprim, xuv)
13
14	use coefpoly_m, only: coefpoly, a1, a2, a3
15	USE dimens_m, ONLY: iim
16	use dynetat0_m, only: clon
17	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	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	DOUBLE PRECISION Y
31	integer i, it
32
33	DOUBLE PRECISION xvrai(iim), Gvrai(iim)
34	! intermediary variables because xlon and xprim are simple precision
35
36	!------------------------------------------------------------------
37
38	it = 0 ! initial guess
39
40	DO i = 1, iim
41	Y = - pi_d + (i + xuv - 0.75d0) * twopi_d / iim
42	! - pi <= y < pi
43	abs_y = abs(y)
44
45	call hunt(xf, abs_y, it, my_lbound = 0)
46	! {0 <= it <= nmax - 1}
47
48	! Calcul de xvrai(i) et Gvrai(i)
49	CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), &
50	xtild(it + 1))
51	xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6)
52	Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
53	if (y < 0d0) xvrai(i) = - xvrai(i)
54	end DO
55
56	DO i = 1, iim -1
57	IF (xvrai(i + 1) < xvrai(i)) THEN
58	print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
59	STOP 1
60	END IF
61	END DO
62
63	xlon = xvrai + clon
64	xprim = twopi_d / (iim * Gvrai)
65
66	end subroutine invert_zoom_x
67
68	!**********************************************************************
69
70	SUBROUTINE funcd(x, fval, fderiv)
71
72	use coefpoly_m, only: a0, a1, a2, a3
73
74	DOUBLE PRECISION, INTENT(IN):: x
75	DOUBLE PRECISION, INTENT(OUT):: fval, fderiv
76
77	fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y
78	fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3)
79
80	END SUBROUTINE funcd
81
82	end module invert_zoom_x_m