/[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 122 by guez, Tue Feb 3 19:30:48 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    
# Line 44  contains Line 43  contains
43    
44      xzoom = clon * pi_d / 180d0      xzoom = clon * pi_d / 180d0
45    
46      IF (grossismx == 1.) THEN      test_grossismx: if (grossismx == 1.) then
47         decalx = 1d0         step = twopi / iim
     else  
        decalx = 0.75d0  
     END IF  
48    
49      dzoom = dzoomx * twopi_d         xprimm025(:iim) = step
50      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)         xprimp025(:iim) = step
51           xprimv(:iim) = step
52           xprimu(:iim) = step
53    
54           rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55           rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57           rlonu(:iim) = rlonv(:iim) + 0.5 * step
58        else
59           dzoom = dzoomx * twopi_d
60           xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61    
62      ! Compute fhyp:         ! Compute fhyp:
63      DO i = nmax, 2 * nmax         DO i = nmax, 2 * nmax
64         fa = taux * (dzoom / 2. - xtild(i))            fa = taux * (dzoom / 2. - xtild(i))
65         fb = xtild(i) * (pi_d - xtild(i))            fb = xtild(i) * (pi_d - xtild(i))
66    
67         IF (200. * fb < - fa) THEN            IF (200. * fb < - fa) THEN
68            fhyp(i) = - 1.               fhyp(i) = - 1.
69         ELSE IF (200. * fb < fa) THEN            ELSE IF (200. * fb < fa) THEN
70            fhyp(i) = 1.               fhyp(i) = 1.
        ELSE  
           IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN  
              IF (200. * fb + fa < 1e-10) THEN  
                 fhyp(i) = - 1.  
              ELSE IF (200. * fb - fa < 1e-10) THEN  
                 fhyp(i) = 1.  
              END IF  
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
126      END IF         END IF
127    
128           ! calcul de Xprimt
129           Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130           xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131    
132           ! Calcul de Xf
133    
134           DO i = nmax + 1, 2 * nmax
135              xmoy = 0.5 * (xtild(i-1) + xtild(i))
136              fa = taux * (dzoom / 2. - xmoy)
137              fb = xmoy * (pi_d - xmoy)
138    
139              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      ! calcul de Xprimt            IF (xmoy == 0.) fxm = 1.
148      Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp            IF (xmoy == pi_d) fxm = -1.
149      xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)            xxpr(i) = beta + (grossismx - beta) * fxm
150           END DO
151      ! Calcul de Xf  
152           xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153      DO i = nmax + 1, 2 * nmax  
154         xmoy = 0.5 * (xtild(i-1) + xtild(i))         Xf(0) = - pi_d
155         fa = taux * (dzoom / 2. - xmoy)  
156         fb = xmoy * (pi_d - xmoy)         DO i=1, 2 * nmax - 1
157              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
158         IF (200. * fb < - fa) THEN         END DO
159            fxm = - 1.  
160         ELSE IF (200. * fb < fa) THEN         Xf(2 * nmax) = pi_d
161            fxm = 1.  
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         ELSE
188            fxm = TANH(fa / fb)            is2 = iim
189    
190              do while (rlonm025(is2) > pi .and. is2 > 1)
191                 is2 = is2 - 1
192              end do
193    
194              if (rlonm025(is2) > pi) then
195                 print *, 'Rlonm025 plus grand que pi !'
196                 STOP 1
197              end if
198         END IF         END IF
199        END IF
200    
201         IF (xmoy == 0.) fxm = 1.      call principal_cshift(is2, rlonm025, xprimm025)
202         IF (xmoy == pi_d) fxm = -1.      call principal_cshift(is2, rlonv, xprimv)
203         xxpr(i) = beta + (grossismx - beta) * fxm      call principal_cshift(is2, rlonu, xprimu)
204      END DO      call principal_cshift(is2, rlonp025, xprimp025)
205    
206      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208      Xf(0) = - pi_d      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
   
     DO i=1, 2 * nmax - 1  
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     END DO  
   
     Xf(2 * nmax) = pi_d  
   
     call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &  
          xprimm025, xuv = - 0.25d0)  
     call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &  
          xuv = 0d0)  
     call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &  
          xuv = 0.5d0)  
     call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &  
          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.122  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21