/[lmdze]/trunk/dyn3d/fxhyp_loop_ik.f
ViewVC logotype

Annotation of /trunk/dyn3d/fxhyp_loop_ik.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 124 - (hide annotations)
Thu Feb 5 15:19:37 2015 UTC (9 years, 3 months ago) by guez
File size: 2574 byte(s) 
Moved some processing from fxhyp_loop_ik to fxhyp. Now fxhyp_loop_ik
does not necessarily give longitudes near [-pi, pi]. In fxhyp, we look
in rlonm025 whether we need to move the array toward [-pi, pi]. If so,
we apply the same move to all grids: rlonm025, rlonv, rlonp025, rlonu
and the corresponding derivatives. The move itself is done by the new
procedure principal_cshift. This revision makes the logic
clearer. (For example, we do not have a saved variable is2 in
fxhyp_loop_ik any longer and we remove a test on ik in fxhyp_loop_ik.)

Fixed a bad error message in fxhyp_loop_ik: talked about rlonu when
xvrai is not always rlonu.
1 guez 121 module fxhyp_loop_ik_m
2    
3     implicit none
4    
5     INTEGER, PARAMETER:: nmax = 30000
6    
7     contains
8    
9     subroutine fxhyp_loop_ik(ik, decalx, xf, xtild, Xprimt, xzoom, xlon, &
10     xprimm, xuv)
11    
12     use coefpoly_m, only: coefpoly
13     USE dimens_m, ONLY: iim
14 guez 124 use nr_util, only: pi_d, twopi_d
15 guez 121 use serre, only: grossismx
16    
17     INTEGER, intent(in):: ik
18     DOUBLE PRECISION, intent(in):: decalx
19     DOUBLE PRECISION, intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax)
20     DOUBLE PRECISION, intent(in):: xzoom
21 guez 124 real, intent(out):: xlon(:), xprimm(:) ! (iim)
22 guez 121
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 guez 123 DOUBLE PRECISION xo1, Xfi, xi, a0, a1, a2, a3, Xf1, Xprimin
29 guez 124 integer ii1, ii2, i, it, iter
30 guez 123 DOUBLE PRECISION, parameter:: my_eps = 1e-6
31 guez 121 DOUBLE PRECISION xxprim(iim + 1), xvrai(iim + 1)
32    
33     !------------------------------------------------------------------
34    
35     IF (ik == 1 .and. grossismx == 1.) THEN
36     ii1 = 2
37     ii2 = iim + 1
38     else
39     ii1=1
40     ii2=iim
41     END IF
42    
43     DO i = ii1, ii2
44     Xfi = - pi_d + (i + xuv - decalx) * twopi_d / iim
45    
46     it = 2 * nmax
47     do while (xfi < xf(it) .and. it >= 1)
48     it = it - 1
49     end do
50    
51     ! Calcul de Xf(xi)
52    
53     xi = xtild(it)
54    
55     IF (it == 2 * nmax) THEN
56     it = 2 * nmax -1
57     END IF
58    
59     CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
60     xtild(it), xtild(it + 1), a0, a1, a2, a3)
61     Xf1 = Xf(it)
62 guez 123 Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)
63     xo1 = xi
64 guez 121 iter = 1
65    
66     do
67     xi = xi - (Xf1 - Xfi) / Xprimin
68     IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
69     xo1 = xi
70 guez 123 Xf1 = a0 + xi * (a1 + xi * (a2 + xi * a3))
71     Xprimin = a1 + xi * (2d0 * a2 + xi * 3d0 * a3)
72 guez 121 end DO
73    
74     if (ABS(xi - xo1) > my_eps) then
75     ! iter == 300
76     print *, 'Pas de solution.'
77     print *, i, xfi
78     STOP 1
79     end if
80    
81     xxprim(i) = twopi_d / (REAL(iim) * Xprimin)
82     xvrai(i) = xi + xzoom
83     end DO
84    
85     IF (ik == 1 .and. grossismx == 1.) THEN
86     xvrai(1) = xvrai(iim + 1)-twopi_d
87     xxprim(1) = xxprim(iim + 1)
88     END IF
89    
90     DO i = 1, iim -1
91     IF (xvrai(i + 1) < xvrai(i)) THEN
92 guez 124 print *, 'xvrai(', i + 1, ') < xvrai(', i, ')'
93 guez 121 STOP 1
94     END IF
95     END DO
96    
97 guez 124 DO i = 1, iim
98     xlon(i) = xvrai(i)
99     xprimm(i) = xxprim(i)
100     END DO
101 guez 121
102     end subroutine fxhyp_loop_ik
103    
104     end module fxhyp_loop_ik_m
  ViewVC Help
Powered by ViewVC 1.1.21