[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 |
---|
[2598] | 13 | use serre_mod, only: clon |
---|
[2218] | 14 | |
---|
| 15 | include "dimensions.h" |
---|
| 16 | ! for iim |
---|
| 17 | |
---|
[2228] | 18 | REAL(K8), intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax) |
---|
[2218] | 19 | real, intent(out):: xlon(:), xprimm(:) ! (iim) |
---|
| 20 | |
---|
[2228] | 21 | REAL(K8), intent(in):: xuv |
---|
[2218] | 22 | ! 0. si calcul aux points scalaires |
---|
| 23 | ! 0.5 si calcul aux points U |
---|
| 24 | |
---|
| 25 | ! Local: |
---|
[2228] | 26 | REAL(K8) xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin |
---|
[2218] | 27 | integer i, it, iter |
---|
[2228] | 28 | REAL(K8), parameter:: my_eps = 1e-6_k8 |
---|
[2218] | 29 | |
---|
[2228] | 30 | REAL(K8) xxprim(iim), xvrai(iim) |
---|
[2218] | 31 | ! intermediary variables because xlon and xprimm are simple precision |
---|
| 32 | |
---|
| 33 | !------------------------------------------------------------------ |
---|
| 34 | |
---|
| 35 | DO i = 1, iim |
---|
[2228] | 36 | Xfi = - pi_d + (i + xuv - 0.75_k8) * twopi_d / iim |
---|
[2218] | 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) |
---|
[2228] | 54 | Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) |
---|
[2218] | 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)) |
---|
[2228] | 63 | Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) |
---|
[2218] | 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 |
---|