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