/[lmdze]/trunk/Sources/dyn3d/fxhyp.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 124 by guez, Thu Feb 5 15:19:37 2015 UTC revision 126 by guez, Fri Feb 6 18:33:15 2015 UTC
# Line 19  contains Line 19  contains
19    
20      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
21      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22      use nr_util, only: pi, pi_d, twopi_d, arth      use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23      use principal_cshift_m, only: principal_cshift      use principal_cshift_m, only: principal_cshift
24      use serre, only: clon, grossismx, dzoomx, taux      use serre, only: clon, grossismx, dzoomx, taux
25    
# Line 28  contains Line 28  contains
28    
29      ! Local:      ! Local:
30      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
31      REAL dzoom      REAL dzoom, step
32      real d_rlonv(iim)      real d_rlonv(iim)
33      DOUBLE PRECISION xtild(0:2 * nmax)      DOUBLE PRECISION xtild(0:2 * nmax)
34      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
# Line 36  contains Line 36  contains
36      DOUBLE PRECISION xzoom, fa, fb      DOUBLE PRECISION xzoom, fa, fb
37      INTEGER i, is2      INTEGER i, is2
38      DOUBLE PRECISION xmoy, fxm      DOUBLE PRECISION xmoy, fxm
     DOUBLE PRECISION decalx  
39    
40      !----------------------------------------------------------------------      !----------------------------------------------------------------------
41    
42      print *, "Call sequence information: fxhyp"      print *, "Call sequence information: fxhyp"
43    
44      dzoom = dzoomx * twopi_d      xzoom = clon * pi_d / 180d0
     xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)  
45    
46      ! Compute fhyp:      test_grossismx: if (grossismx == 1.) then
47      DO i = nmax, 2 * nmax         step = twopi / iim
48         fa = taux * (dzoom / 2. - xtild(i))  
49         fb = xtild(i) * (pi_d - xtild(i))         xprimm025(:iim) = step
50           xprimp025(:iim) = step
51         IF (200. * fb < - fa) THEN         xprimv(:iim) = step
52            fhyp(i) = - 1.         xprimu(:iim) = step
53         ELSE IF (200. * fb < fa) THEN  
54            fhyp(i) = 1.         rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55         ELSE         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57               IF (200. * fb + fa < 1e-10) THEN         rlonu(:iim) = rlonv(:iim) + 0.5 * step
58                  fhyp(i) = - 1.      else
59               ELSE IF (200. * fb - fa < 1e-10) THEN         dzoom = dzoomx * twopi_d
60                  fhyp(i) = 1.         xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61               END IF  
62           ! Compute fhyp:
63           DO i = nmax, 2 * nmax
64              fa = taux * (dzoom / 2. - xtild(i))
65              fb = xtild(i) * (pi_d - xtild(i))
66    
67              IF (200. * fb < - fa) THEN
68                 fhyp(i) = - 1.
69              ELSE IF (200. * fb < fa) THEN
70                 fhyp(i) = 1.
71            ELSE            ELSE
72               fhyp(i) = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
73                    IF (200. * fb + fa < 1e-10) THEN
74                       fhyp(i) = - 1.
75                    ELSE IF (200. * fb - fa < 1e-10) THEN
76                       fhyp(i) = 1.
77                    END IF
78                 ELSE
79                    fhyp(i) = TANH(fa / fb)
80                 END IF
81            END IF            END IF
        END IF  
82    
83         IF (xtild(i) == 0.) fhyp(i) = 1.            IF (xtild(i) == 0.) fhyp(i) = 1.
84         IF (xtild(i) == pi_d) fhyp(i) = -1.            IF (xtild(i) == pi_d) fhyp(i) = -1.
85      END DO         END DO
86    
87      ! Calcul de beta         ! Calcul de beta
88    
89      ffdx = 0.         ffdx = 0.
90    
91      DO i = nmax + 1, 2 * nmax         DO i = nmax + 1, 2 * nmax
92         xmoy = 0.5 * (xtild(i-1) + xtild(i))            xmoy = 0.5 * (xtild(i-1) + xtild(i))
93         fa = taux * (dzoom / 2. - xmoy)            fa = taux * (dzoom / 2. - xmoy)
94         fb = xmoy * (pi_d - xmoy)            fb = xmoy * (pi_d - xmoy)
95    
96         IF (200. * fb < - fa) THEN            IF (200. * fb < - fa) THEN
97            fxm = - 1.               fxm = - 1.
98         ELSE IF (200. * fb < fa) THEN            ELSE IF (200. * fb < fa) THEN
99            fxm = 1.               fxm = 1.
        ELSE  
           IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN  
              IF (200. * fb + fa < 1e-10) THEN  
                 fxm = - 1.  
              ELSE IF (200. * fb - fa < 1e-10) THEN  
                 fxm = 1.  
              END IF  
100            ELSE            ELSE
101               fxm = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
102                    IF (200. * fb + fa < 1e-10) THEN
103                       fxm = - 1.
104                    ELSE IF (200. * fb - fa < 1e-10) THEN
105                       fxm = 1.
106                    END IF
107                 ELSE
108                    fxm = TANH(fa / fb)
109                 END IF
110            END IF            END IF
        END IF  
111    
112         IF (xmoy == 0.) fxm = 1.            IF (xmoy == 0.) fxm = 1.
113         IF (xmoy == pi_d) fxm = -1.            IF (xmoy == pi_d) fxm = -1.
114    
115         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))            ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
116      END DO         END DO
117    
118      print *, "ffdx = ", ffdx         print *, "ffdx = ", ffdx
119      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)         beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
120      print *, "beta = ", beta         print *, "beta = ", beta
121    
122      IF (2. * beta - grossismx <= 0.) THEN         IF (2. * beta - grossismx <= 0.) THEN
123         print *, 'Bad choice of grossismx, taux, dzoomx.'            print *, 'Bad choice of grossismx, taux, dzoomx.'
124         print *, 'Decrease dzoomx or grossismx.'            print *, 'Decrease dzoomx or grossismx.'
125         STOP 1            STOP 1
     END IF  
   
     ! calcul de Xprimt  
     Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp  
     xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)  
   
     ! Calcul de Xf  
   
     DO i = nmax + 1, 2 * nmax  
        xmoy = 0.5 * (xtild(i-1) + xtild(i))  
        fa = taux * (dzoom / 2. - xmoy)  
        fb = xmoy * (pi_d - xmoy)  
   
        IF (200. * fb < - fa) THEN  
           fxm = - 1.  
        ELSE IF (200. * fb < fa) THEN  
           fxm = 1.  
        ELSE  
           fxm = TANH(fa / fb)  
126         END IF         END IF
127    
128         IF (xmoy == 0.) fxm = 1.         ! calcul de Xprimt
129         IF (xmoy == pi_d) fxm = -1.         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130         xxpr(i) = beta + (grossismx - beta) * fxm         xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131      END DO  
132           ! Calcul de Xf
133      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)  
134           DO i = nmax + 1, 2 * nmax
135      Xf(0) = - pi_d            xmoy = 0.5 * (xtild(i-1) + xtild(i))
136              fa = taux * (dzoom / 2. - xmoy)
137      DO i=1, 2 * nmax - 1            fb = xmoy * (pi_d - xmoy)
138         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
139      END DO            IF (200. * fb < - fa) THEN
140                 fxm = - 1.
141      Xf(2 * nmax) = pi_d            ELSE IF (200. * fb < fa) THEN
142                 fxm = 1.
143      IF (grossismx == 1.) THEN            ELSE
144         decalx = 1d0               fxm = TANH(fa / fb)
145      else            END IF
        decalx = 0.75d0  
     END IF  
146    
147      xzoom = clon * pi_d / 180d0            IF (xmoy == 0.) fxm = 1.
148      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025(:iim), &            IF (xmoy == pi_d) fxm = -1.
149           xprimm025(:iim), xuv = - 0.25d0)            xxpr(i) = beta + (grossismx - beta) * fxm
150      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv(:iim), &         END DO
151           xprimv(:iim), xuv = 0d0)  
152      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu(:iim), &         xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153           xprimu(:iim), xuv = 0.5d0)  
154      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025(:iim), &         Xf(0) = - pi_d
155           xprimp025(:iim), xuv = 0.25d0)  
156           DO i=1, 2 * nmax - 1
157              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
158           END DO
159    
160           Xf(2 * nmax) = pi_d
161    
162           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonm025(:iim), &
163                xprimm025(:iim), xuv = - 0.25d0)
164           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonv(:iim), &
165                xprimv(:iim), xuv = 0d0)
166           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonu(:iim), &
167                xprimu(:iim), xuv = 0.5d0)
168           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonp025(:iim), &
169                xprimp025(:iim), xuv = 0.25d0)
170        end if test_grossismx
171    
172      is2 = 0      is2 = 0
173    
Legend:
Removed from v.124  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21