/[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 123 by guez, Thu Feb 5 12:41:08 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_d, twopi_d, arth      use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23        use principal_cshift_m, only: principal_cshift
24      use serre, only: clon, grossismx, dzoomx, taux      use serre, only: clon, grossismx, dzoomx, taux
25    
26      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
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      DOUBLE PRECISION xlon(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)
35      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      DOUBLE PRECISION xzoom, fa, fb      DOUBLE PRECISION xzoom, fa, fb
37      INTEGER i      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              ELSE IF (200. * fb < fa) THEN
142                 fxm = 1.
143              ELSE
144                 fxm = TANH(fa / fb)
145              END IF
146    
147      Xf(2 * nmax) = pi_d            IF (xmoy == 0.) fxm = 1.
148              IF (xmoy == pi_d) fxm = -1.
149              xxpr(i) = beta + (grossismx - beta) * fxm
150           END DO
151    
152           xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153    
154           Xf(0) = - pi_d
155    
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
173    
174        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
175             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
176           IF (clon <= 0.) THEN
177              is2 = 1
178    
179              do while (rlonm025(is2) < - pi .and. is2 < iim)
180                 is2 = is2 + 1
181              end do
182    
183              if (rlonm025(is2) < - pi) then
184                 print *, 'Rlonm025 plus petit que - pi !'
185                 STOP 1
186              end if
187           ELSE
188              is2 = iim
189    
190      IF (grossismx == 1.) THEN            do while (rlonm025(is2) > pi .and. is2 > 1)
191         decalx = 1d0               is2 = is2 - 1
192      else            end do
193         decalx = 0.75d0  
194              if (rlonm025(is2) > pi) then
195                 print *, 'Rlonm025 plus grand que pi !'
196                 STOP 1
197              end if
198           END IF
199      END IF      END IF
200    
201      xzoom = clon * pi_d / 180d0      call principal_cshift(is2, rlonm025, xprimm025)
202      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &      call principal_cshift(is2, rlonv, xprimv)
203           xprimm025, xuv = - 0.25d0)      call principal_cshift(is2, rlonu, xprimu)
204      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &      call principal_cshift(is2, rlonp025, xprimp025)
205           xuv = 0d0)  
206      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207           xuv = 0.5d0)      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
          xprimp025, xuv = 0.25d0)  
   
     print *  
   
     forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)  
     print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "degrees"  
     print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "degrees"  
209    
210      DO i = 1, iim + 1      DO i = 1, iim + 1
211         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
Legend:
Removed from v.123  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21