/[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

trunk/dyn3d/fxhyp.f revision 120 by guez, Tue Jan 13 14:56:15 2015 UTC trunk/Sources/dyn3d/fxhyp.f revision 139 by guez, Tue May 26 17:46:03 2015 UTC
# Line 12  contains Line 12  contains
12      ! Calcule les longitudes et dérivées dans la grille du GCM pour      ! Calcule les longitudes et dérivées dans la grille du GCM pour
13      ! une fonction f(x) à dérivée tangente hyperbolique.      ! une fonction f(x) à dérivée tangente hyperbolique.
14    
15      ! On doit avoir grossismx \times dzoomx < pi (radians)      ! Il vaut mieux avoir : grossismx \times dzoom < pi
16    
17      ! Le premier point scalaire pour une grille regulière (grossismx =      ! Le premier point scalaire pour une grille regulière (grossismx =
18      ! 1., taux=0., clon=0.) est à - 180 degrés.      ! 1., taux=0., clon=0.) est à - 180 degrés.
19    
     use coefpoly_m, only: coefpoly  
20      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
21      use nr_util, only: pi_d, twopi_d, arth      use dynetat0_m, only: clon, grossismx, dzoomx, taux
22      use serre, only: clon, grossismx, dzoomx, taux      use invert_zoom_x_m, only: invert_zoom_x, nmax
23        use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24        use principal_cshift_m, only: principal_cshift
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:
   
     DOUBLE PRECISION champmin, champmax  
30      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
31      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax      REAL dzoom, step
32      REAL dzoom      real d_rlonv(iim)
33      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      DOUBLE PRECISION xtild(0:2 * nmax)
34      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35      DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)      DOUBLE PRECISION fa, fb
37      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)      INTEGER i, is2
38      DOUBLE PRECISION my_eps, xzoom, fa, fb      DOUBLE PRECISION xmoy, fxm
     DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2  
     INTEGER i, it, ik, iter, ii, idif, ii1, ii2  
     DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin  
     DOUBLE PRECISION decalx  
     INTEGER is2  
39    
40      !----------------------------------------------------------------------      !----------------------------------------------------------------------
41    
42      print *, "Call sequence information: fxhyp"      print *, "Call sequence information: fxhyp"
43    
44      my_eps = 1e-3      test_grossismx: if (grossismx == 1.) then
45      xzoom = clon * pi_d / 180.         step = twopi / iim
   
     IF (grossismx == 1.) THEN  
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
46    
47      IF (dzoomx < 1.) THEN         xprimm025(:iim) = step
48           xprimp025(:iim) = step
49           xprimv(:iim) = step
50           xprimu(:iim) = step
51    
52           rlonv(:iim) = arth(- pi + clon, step, iim)
53           rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55           rlonu(:iim) = rlonv(:iim) + 0.5 * step
56        else test_grossismx
57         dzoom = dzoomx * twopi_d         dzoom = dzoomx * twopi_d
58      ELSE IF (dzoomx < 25.) THEN         xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
        print *, "dzoomx pour fxhyp est trop petit."  
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
     END IF  
   
     print *, 'dzoom (rad):', dzoom  
59    
60      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)         ! Compute fhyp:
61           DO i = nmax, 2 * nmax
62      DO i = nmax, nmax2            fa = taux * (dzoom / 2. - xtild(i))
63         fa = taux * (dzoom / 2. - xtild(i))            fb = xtild(i) * (pi_d - xtild(i))
64         fb = xtild(i) * (pi_d - xtild(i))  
65              IF (200. * fb < - fa) THEN
66         IF (200. * fb < - fa) THEN               fhyp(i) = - 1.
67            fhyp(i) = - 1.            ELSE IF (200. * fb < fa) THEN
68         ELSE IF (200. * fb < fa) THEN               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  
69            ELSE            ELSE
70               fhyp(i) = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71                    IF (200. * fb + fa < 1e-10) THEN
72                       fhyp(i) = - 1.
73                    ELSE IF (200. * fb - fa < 1e-10) THEN
74                       fhyp(i) = 1.
75                    END IF
76                 ELSE
77                    fhyp(i) = TANH(fa / fb)
78                 END IF
79            END IF            END IF
        END IF  
80    
81         IF (xtild(i) == 0.) fhyp(i) = 1.            IF (xtild(i) == 0.) fhyp(i) = 1.
82         IF (xtild(i) == pi_d) fhyp(i) = -1.            IF (xtild(i) == pi_d) fhyp(i) = -1.
83      END DO         END DO
84    
85      ! Calcul de beta         ! Calcul de beta
86    
87      ffdx = 0.         ffdx = 0.
88    
89      DO i = nmax + 1, nmax2         DO i = nmax + 1, 2 * nmax
90         xmoy = 0.5 * (xtild(i-1) + xtild(i))            xmoy = 0.5 * (xtild(i-1) + xtild(i))
91         fa = taux * (dzoom / 2. - xmoy)            fa = taux * (dzoom / 2. - xmoy)
92         fb = xmoy * (pi_d - xmoy)            fb = xmoy * (pi_d - xmoy)
93    
94         IF (200. * fb < - fa) THEN            IF (200. * fb < - fa) THEN
95            fxm = - 1.               fxm = - 1.
96         ELSE IF (200. * fb < fa) THEN            ELSE IF (200. * fb < fa) THEN
97            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  
98            ELSE            ELSE
99               fxm = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100                    IF (200. * fb + fa < 1e-10) THEN
101                       fxm = - 1.
102                    ELSE IF (200. * fb - fa < 1e-10) THEN
103                       fxm = 1.
104                    END IF
105                 ELSE
106                    fxm = TANH(fa / fb)
107                 END IF
108            END IF            END IF
        END IF  
   
        IF (xmoy == 0.) fxm = 1.  
        IF (xmoy == pi_d) fxm = -1.  
   
        ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))  
     END DO  
   
     beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)  
109    
110      IF (2. * beta - grossismx <= 0.) THEN            IF (xmoy == 0.) fxm = 1.
111         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'            IF (xmoy == pi_d) fxm = -1.
        print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'  
        STOP 1  
     END IF  
112    
113      ! calcul de Xprimt            ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114           END DO
115    
116      DO i = nmax, nmax2         print *, "ffdx = ", ffdx
117         Xprimt(i) = beta + (grossismx - beta) * fhyp(i)         beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118      END DO         print *, "beta = ", beta
119    
120           IF (2. * beta - grossismx <= 0.) THEN
121              print *, 'Bad choice of grossismx, taux, dzoomx.'
122              print *, 'Decrease dzoomx or grossismx.'
123              STOP 1
124           END IF
125    
126      DO i = nmax + 1, nmax2         ! calcul de Xprimt
127         Xprimt(nmax2 - i) = Xprimt(i)         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128      END DO         xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129    
130           ! Calcul de Xf
131    
132           DO i = nmax + 1, 2 * nmax
133              xmoy = 0.5 * (xtild(i-1) + xtild(i))
134              fa = taux * (dzoom / 2. - xmoy)
135              fb = xmoy * (pi_d - xmoy)
136    
137              IF (200. * fb < - fa) THEN
138                 fxm = - 1.
139              ELSE IF (200. * fb < fa) THEN
140                 fxm = 1.
141              ELSE
142                 fxm = TANH(fa / fb)
143              END IF
144    
145      ! Calcul de Xf            IF (xmoy == 0.) fxm = 1.
146              IF (xmoy == pi_d) fxm = -1.
147              xxpr(i) = beta + (grossismx - beta) * fxm
148           END DO
149    
150      Xf(0) = - pi_d         xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151    
152      DO i = nmax + 1, nmax2         Xf(0) = - pi_d
        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)  
        END IF  
153    
154         IF (xmoy == 0.) fxm = 1.         DO i=1, 2 * nmax - 1
155         IF (xmoy == pi_d) fxm = -1.            Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156         xxpr(i) = beta + (grossismx - beta) * fxm         END DO
     END DO  
157    
158      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)         Xf(2 * nmax) = pi_d
159    
160      DO i=1, nmax2         call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))              xuv = - 0.25d0)
162      END DO         call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163                xuv = 0d0)
164           call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165                xuv = 0.5d0)
166           call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167                xuv = 0.25d0)
168        end if test_grossismx
169    
170      is2 = 0      is2 = 0
171    
172      loop_ik: DO ik = 1, 4      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
173         ! xuv = 0. si calcul aux points scalaires           .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
174         ! xuv = 0.5 si calcul aux points U         IF (clon <= 0.) THEN
175              is2 = 1
        IF (ik == 1) THEN  
           xuv = -0.25  
        ELSE IF (ik == 2) THEN  
           xuv = 0.  
        ELSE IF (ik == 3) THEN  
           xuv = 0.50  
        ELSE IF (ik == 4) THEN  
           xuv = 0.25  
        END IF  
   
        xo1 = 0.  
176    
177         IF (ik == 1 .and. grossismx == 1.) THEN            do while (rlonm025(is2) < - pi .and. is2 < iim)
178            ii1 = 2               is2 = is2 + 1
           ii2 = iim + 1  
        else  
           ii1=1  
           ii2=iim  
        END IF  
   
        DO i = ii1, ii2  
           Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)  
   
           it = nmax2  
           do while (xfi < xf(it) .and. it >= 1)  
              it = it - 1  
179            end do            end do
180    
181            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
182                 print *, 'Rlonm025 plus petit que - pi !'
           xi = xtild(it)  
   
           IF (it == nmax2) THEN  
              it = nmax2 -1  
              Xf(it + 1) = pi_d  
           END IF  
   
           ! Appel de la routine qui calcule les coefficients a0, a1,  
           ! a2, a3 d'un polynome de degre 3 qui passe par les points  
           ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))  
   
           CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &  
                xtild(it), xtild(it + 1), a0, a1, a2, a3)  
   
           Xf1 = Xf(it)  
           Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi * xi  
   
           iter = 1  
   
           do  
              xi = xi - (Xf1 - Xfi) / Xprimin  
              IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit  
              xo1 = xi  
              xi2 = xi * xi  
              Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi  
              Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2  
           end DO  
   
           if (ABS(xi - xo1) > my_eps) then  
              ! iter == 300  
              print *, 'Pas de solution.'  
              print *, i, xfi  
183               STOP 1               STOP 1
184            end if            end if
185           ELSE
186              is2 = iim
187    
188            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
189            xvrai(i) = xi + xzoom               is2 = is2 - 1
190         end DO            end do
   
        IF (ik == 1 .and. grossismx == 1.) THEN  
           xvrai(1) = xvrai(iim + 1)-twopi_d  
           xxprim(1) = xxprim(iim + 1)  
        END IF  
   
        DO i = 1, iim  
           xlon(i) = xvrai(i)  
           xprimm(i) = xxprim(i)  
        END DO  
191    
192         DO i = 1, iim -1            if (rlonm025(is2) > pi) then
193            IF (xvrai(i + 1) < xvrai(i)) THEN               print *, 'Rlonm025 plus grand que pi !'
              print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'  
194               STOP 1               STOP 1
195            END IF            end if
        END DO  
   
        IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 0.1 &  
             .and. MAXval(xvrai(:iim)) <= pi_d + 0.1)) THEN  
           print *, &  
                'Réorganisation des longitudes pour les avoir entre - pi et pi'  
   
           IF (xzoom <= 0.) THEN  
              IF (ik == 1) THEN  
                 i = 1  
   
                 do while (xvrai(i) < - pi_d .and. i < iim)  
                    i = i + 1  
                 end do  
   
                 if (xvrai(i) < - pi_d) then  
                    print *, 'Xvrai plus petit que - pi !'  
                    STOP 1  
                 end if  
   
                 is2 = i  
              END IF  
   
              IF (is2 /= 1) THEN  
                 DO ii = is2, iim  
                    xlon(ii-is2 + 1) = xvrai(ii)  
                    xprimm(ii-is2 + 1) = xxprim(ii)  
                 END DO  
                 DO ii = 1, is2 -1  
                    xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d  
                    xprimm(ii + iim-is2 + 1) = xxprim(ii)  
                 END DO  
              END IF  
           ELSE  
              IF (ik == 1) THEN  
                 i = iim  
   
                 do while (xvrai(i) > pi_d .and. i > 1)  
                    i = i - 1  
                 end do  
   
                 if (xvrai(i) > pi_d) then  
                    print *, 'Xvrai plus grand que pi !'  
                    STOP 1  
                 end if  
   
                 is2 = i  
              END IF  
   
              idif = iim -is2  
   
              DO ii = 1, is2  
                 xlon(ii + idif) = xvrai(ii)  
                 xprimm(ii + idif) = xxprim(ii)  
              END DO  
   
              DO ii = 1, idif  
                 xlon(ii) = xvrai(ii + is2) - twopi_d  
                 xprimm(ii) = xxprim(ii + is2)  
              END DO  
           END IF  
        END IF  
   
        xlon(iim + 1) = xlon(1) + twopi_d  
        xprimm(iim + 1) = xprimm(1)  
   
        DO i = 1, iim + 1  
           xvrai(i) = xlon(i) * 180. / pi_d  
        END DO  
   
        IF (ik == 1) THEN  
           DO i = 1, iim + 1  
              rlonm025(i) = xlon(i)  
              xprimm025(i) = xprimm(i)  
           END DO  
        ELSE IF (ik == 2) THEN  
           rlonv = xlon  
           xprimv = xprimm  
        ELSE IF (ik == 3) THEN  
           DO i = 1, iim + 1  
              rlonu(i) = xlon(i)  
              xprimu(i) = xprimm(i)  
           END DO  
        ELSE IF (ik == 4) THEN  
           rlonp025 = xlon  
           xprimp025 = xprimm  
196         END IF         END IF
197      end DO loop_ik      END IF
198    
199      print *      call principal_cshift(is2, rlonm025, xprimm025)
200        call principal_cshift(is2, rlonv, xprimv)
201      DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
202         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonp025, xprimp025)
203      END DO  
204      champmin = 1e12      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205      champmax = -1e12      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
206      DO i = 1, iim      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
        champmin = MIN(champmin, xlon(i))  
        champmax = MAX(champmax, xlon(i))  
     END DO  
     champmin = champmin * 180. / pi_d  
     champmax = champmax * 180. / pi_d  
207    
208        ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209      DO i = 1, iim + 1      DO i = 1, iim + 1
210         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
211            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
212              print *, "< rlonv(", i, ") = ", rlonv(i)
213            STOP 1            STOP 1
214         END IF         END IF
215    
216         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
217            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
218              print *, "< rlonm025(", i, ") = ", rlonm025(i)
219            STOP 1            STOP 1
220         END IF         END IF
221    
# Line 381  contains Line 226  contains
226         END IF         END IF
227      END DO      END DO
228    
     print *, ' Longitudes '  
     print 3, champmin, champmax  
   
 3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &  
          ' d environ ', f0.2, ' degres ', /, &  
          ' alors que la maille en dehors de la zone du zoom est ', &  
          "d'environ ", f0.2, ' degres ')  
   
229    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
230    
231  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.139
  ViewVC Help
Powered by ViewVC 1.1.21