[2218] | 1 | module invert_zoom_x_m |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | INTEGER, PARAMETER:: nmax = 30000 |
---|
| 6 | |
---|
| 7 | contains |
---|
| 8 | |
---|
| 9 | subroutine invert_zoom_x(xf, xtild, Xprimt, xlon, xprimm, xuv) |
---|
| 10 | |
---|
| 11 | use coefpoly_m, only: coefpoly |
---|
[2228] | 12 | use nrtype, only: pi, pi_d, twopi_d, k8 |
---|
[2218] | 13 | |
---|
| 14 | include "dimensions.h" |
---|
| 15 | ! for iim |
---|
| 16 | |
---|
| 17 | include "serre.h" |
---|
| 18 | ! for clon |
---|
| 19 | |
---|
[2228] | 20 | REAL(K8), intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax) |
---|
[2218] | 21 | real, intent(out):: xlon(:), xprimm(:) ! (iim) |
---|
| 22 | |
---|
[2228] | 23 | REAL(K8), intent(in):: xuv |
---|
[2218] | 24 | ! 0. si calcul aux points scalaires |
---|
| 25 | ! 0.5 si calcul aux points U |
---|
| 26 | |
---|
| 27 | ! Local: |
---|
[2228] | 28 | REAL(K8) xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin |
---|
[2218] | 29 | integer i, it, iter |
---|
[2228] | 30 | REAL(K8), parameter:: my_eps = 1e-6_k8 |
---|
[2218] | 31 | |
---|
[2228] | 32 | REAL(K8) xxprim(iim), xvrai(iim) |
---|
[2218] | 33 | ! intermediary variables because xlon and xprimm are simple precision |
---|
| 34 | |
---|
| 35 | !------------------------------------------------------------------ |
---|
| 36 | |
---|
| 37 | DO i = 1, iim |
---|
[2228] | 38 | Xfi = - pi_d + (i + xuv - 0.75_k8) * twopi_d / iim |
---|
[2218] | 39 | |
---|
| 40 | it = 2 * nmax |
---|
| 41 | do while (xfi < xf(it) .and. it >= 1) |
---|
| 42 | it = it - 1 |
---|
| 43 | end do |
---|
| 44 | |
---|
| 45 | ! Calcul de Xf(xvrai(i)) |
---|
| 46 | |
---|
| 47 | xvrai(i) = xtild(it) |
---|
| 48 | |
---|
| 49 | IF (it == 2 * nmax) THEN |
---|
| 50 | it = 2 * nmax -1 |
---|
| 51 | END IF |
---|
| 52 | |
---|
| 53 | CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), & |
---|
| 54 | xtild(it), xtild(it + 1), a0, a1, a2, a3) |
---|
| 55 | Xf1 = Xf(it) |
---|
[2228] | 56 | Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) |
---|
[2218] | 57 | xo1 = xvrai(i) |
---|
| 58 | iter = 1 |
---|
| 59 | |
---|
| 60 | do |
---|
| 61 | xvrai(i) = xvrai(i) - (Xf1 - Xfi) / Xprimin |
---|
| 62 | IF (ABS(xvrai(i) - xo1) <= my_eps .or. iter == 300) exit |
---|
| 63 | xo1 = xvrai(i) |
---|
| 64 | Xf1 = a0 + xvrai(i) * (a1 + xvrai(i) * (a2 + xvrai(i) * a3)) |
---|
[2228] | 65 | Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) |
---|
[2218] | 66 | end DO |
---|
| 67 | |
---|
| 68 | if (ABS(xvrai(i) - xo1) > my_eps) then |
---|
| 69 | ! iter == 300 |
---|
| 70 | print *, 'Pas de solution.' |
---|
| 71 | print *, i, xfi |
---|
| 72 | STOP 1 |
---|
| 73 | end if |
---|
| 74 | |
---|
| 75 | xxprim(i) = twopi_d / (iim * Xprimin) |
---|
| 76 | end DO |
---|
| 77 | |
---|
| 78 | DO i = 1, iim -1 |
---|
| 79 | IF (xvrai(i + 1) < xvrai(i)) THEN |
---|
| 80 | print *, 'xvrai(', i + 1, ') < xvrai(', i, ')' |
---|
| 81 | STOP 1 |
---|
| 82 | END IF |
---|
| 83 | END DO |
---|
| 84 | |
---|
| 85 | xlon = xvrai + clon / 180. * pi |
---|
| 86 | xprimm = xxprim |
---|
| 87 | |
---|
| 88 | end subroutine invert_zoom_x |
---|
| 89 | |
---|
| 90 | end module invert_zoom_x_m |
---|