source: LMDZ5/trunk/libf/dyn3d_common/invert_zoom_x_m.F90 @ 2218

Last change on this file since 2218 was 2218, checked in by lguez, 10 years ago

Bug fix in fxhyp. There was some bad logic for the shifting of
longitude grids rlonuv, rlonv, rlonm025 and rlonp025 toward [- pi,
pi]. In some cases, one of the four grids was not shifted and the
others were. For example, you could see the bug with: iim = 48, clon =
20, grossismx = 3, dzoomx = 0.15. The bad logic involved the variable
is2 in the loop on ik = 1, 4. The loop included some tests depending
on ik. The whole thing was quite convoluted. Rewrote fxhyp. In
particular, replaced the loop on ik by a call to a new procedure,
invert_zoom_x. fxhyp.F was in dyn3d, it is replaced by fxhyp_m.F90
which is in dyn3d_common. Removed several arguments of fxhyp: zoom
parameters are accessed through serre.h; rlonm025 and rlonp025 were
not used in inigeom; min and max of longitude steps are written inside
fxhyp.

Some simplifications and modernizations in fyhyp. In particular,
several arguments are removed: zoom parameters are accessed through
serre.h; yprimv was not used in inigeom; min and max of latitude steps
are written inside fyhyp.

Removed now useless intermediary procedure fxyhyper.

Changed default value of dzoomx and dzoomy from 0 to 0.2. dzoom[xy] is
only needed when grossism[xy] > 1 and then it should be > 0.

dzoom[xy] can no longer be the extent of the zoom in degrees: it must
be expressed as the fraction of the total domain.
File size: 2.2 KB
Line 
1module invert_zoom_x_m
2
3  implicit none
4
5  INTEGER, PARAMETER:: nmax = 30000
6
7contains
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
13
14    include "dimensions.h"
15    ! for iim
16
17    include "serre.h"
18    ! for clon
19
20    DOUBLE PRECISION, intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax)
21    real, intent(out):: xlon(:), xprimm(:) ! (iim)
22
23    DOUBLE PRECISION, intent(in):: xuv
24    ! 0. si calcul aux points scalaires
25    ! 0.5 si calcul aux points U
26
27    ! Local:
28    DOUBLE PRECISION xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin
29    integer i, it, iter
30    DOUBLE PRECISION, parameter:: my_eps = 1d-6
31
32    DOUBLE PRECISION xxprim(iim), xvrai(iim)
33    ! intermediary variables because xlon and xprimm are simple precision
34
35    !------------------------------------------------------------------
36
37    DO i = 1, iim
38       Xfi = - pi_d + (i + xuv - 0.75d0) * twopi_d / iim
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)
56       Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
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))
65          Xprimin = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3)
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
90end module invert_zoom_x_m
Note: See TracBrowser for help on using the repository browser.