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

Diff of /trunk/dyn3d/fxhyp.f

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

revision 120 by guez, Tue Jan 13 14:56:15 2015 UTC revision 121 by guez, Wed Jan 28 16:10:02 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 fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22      use nr_util, only: pi_d, twopi_d, arth      use nr_util, only: pi_d, twopi_d, arth
23      use serre, only: clon, grossismx, dzoomx, taux      use serre, only: clon, grossismx, dzoomx, taux
24    
# Line 27  contains Line 27  contains
27    
28      ! Local:      ! Local:
29    
     DOUBLE PRECISION champmin, champmax  
30      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
     INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax  
31      REAL dzoom      REAL dzoom
32      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      DOUBLE PRECISION xlon(iim)
33      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION xtild(0:2 * nmax)
34      DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35      DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)      DOUBLE PRECISION xzoom, fa, fb
37      DOUBLE PRECISION my_eps, xzoom, fa, fb      INTEGER i
38      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2      DOUBLE PRECISION xmoy, fxm
     INTEGER i, it, ik, iter, ii, idif, ii1, ii2  
     DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin  
39      DOUBLE PRECISION decalx      DOUBLE PRECISION decalx
     INTEGER is2  
40    
41      !----------------------------------------------------------------------      !----------------------------------------------------------------------
42    
43      print *, "Call sequence information: fxhyp"      print *, "Call sequence information: fxhyp"
44    
45      my_eps = 1e-3      xzoom = clon * pi_d / 180d0
     xzoom = clon * pi_d / 180.  
46    
47      IF (grossismx == 1.) THEN      IF (grossismx == 1.) THEN
48         decalx = 1.         decalx = 1d0
49      else      else
50         decalx = 0.75         decalx = 0.75d0
51      END IF      END IF
52    
53      IF (dzoomx < 1.) THEN      dzoom = dzoomx * twopi_d
54         dzoom = dzoomx * twopi_d      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
     ELSE IF (dzoomx < 25.) THEN  
        print *, "dzoomx pour fxhyp est trop petit."  
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
     END IF  
   
     print *, 'dzoom (rad):', dzoom  
55    
56      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)      ! Compute fhyp:
57        DO i = nmax, 2 * nmax
     DO i = nmax, nmax2  
58         fa = taux * (dzoom / 2. - xtild(i))         fa = taux * (dzoom / 2. - xtild(i))
59         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
60    
# Line 78  contains Line 63  contains
63         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
64            fhyp(i) = 1.            fhyp(i) = 1.
65         ELSE         ELSE
66            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
67               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
68                  fhyp(i) = - 1.                  fhyp(i) = - 1.
69               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 97  contains Line 82  contains
82    
83      ffdx = 0.      ffdx = 0.
84    
85      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
86         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
87         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
88         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 107  contains Line 92  contains
92         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
93            fxm = 1.            fxm = 1.
94         ELSE         ELSE
95            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
96               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
97                  fxm = - 1.                  fxm = - 1.
98               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 124  contains Line 109  contains
109         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
110      END DO      END DO
111    
112        print *, "ffdx = ", ffdx
113      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
114        print *, "beta = ", beta
115    
116      IF (2. * beta - grossismx <= 0.) THEN      IF (2. * beta - grossismx <= 0.) THEN
117         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
# Line 133  contains Line 120  contains
120      END IF      END IF
121    
122      ! calcul de Xprimt      ! calcul de Xprimt
123        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
124        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
125    
126      DO i = nmax, nmax2      ! Calcul de Xf
        Xprimt(i) = beta + (grossismx - beta) * fhyp(i)  
     END DO  
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
   
     ! Calcul de Xf  
   
     Xf(0) = - pi_d  
127    
128      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
129         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
130         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
131         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 164  contains Line 143  contains
143         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
144      END DO      END DO
145    
146      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
147    
148        Xf(0) = - pi_d
149    
150      DO i=1, nmax2      DO i=1, 2 * nmax - 1
151         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
152      END DO      END DO
153    
154      is2 = 0      Xf(2 * nmax) = pi_d
155    
156      loop_ik: DO ik = 1, 4      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
157         ! xuv = 0. si calcul aux points scalaires           xprimm025, xuv = - 0.25d0)
158         ! xuv = 0.5 si calcul aux points U      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
159             xuv = 0d0)
160         IF (ik == 1) THEN      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
161            xuv = -0.25           xuv = 0.5d0)
162         ELSE IF (ik == 2) THEN      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
163            xuv = 0.           xprimp025, xuv = 0.25d0)
        ELSE IF (ik == 3) THEN  
           xuv = 0.50  
        ELSE IF (ik == 4) THEN  
           xuv = 0.25  
        END IF  
   
        xo1 = 0.  
   
        IF (ik == 1 .and. grossismx == 1.) THEN  
           ii1 = 2  
           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  
           end do  
   
           ! Calcul de Xf(xi)  
   
           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  
              STOP 1  
           end if  
   
           xxprim(i) = twopi_d / (REAL(iim) * Xprimin)  
           xvrai(i) = xi + xzoom  
        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  
   
        DO i = 1, iim -1  
           IF (xvrai(i + 1) < xvrai(i)) THEN  
              print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'  
              STOP 1  
           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  
        END IF  
     end DO loop_ik  
164    
165      print *      print *
166    
167      DO i = 1, iim      forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
168         xlon(i) = rlonv(i + 1) - rlonv(i)      print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "°"
169      END DO      print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "°"
     champmin = 1e12  
     champmax = -1e12  
     DO i = 1, iim  
        champmin = MIN(champmin, xlon(i))  
        champmax = MAX(champmax, xlon(i))  
     END DO  
     champmin = champmin * 180. / pi_d  
     champmax = champmax * 180. / pi_d  
170    
171      DO i = 1, iim + 1      DO i = 1, iim + 1
172         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
173            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
174              print *, "< rlonv(", i, ") = ", rlonv(i)
175            STOP 1            STOP 1
176         END IF         END IF
177    
178         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
179            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
180              print *, "< rlonm025(", i, ") = ", rlonm025(i)
181            STOP 1            STOP 1
182         END IF         END IF
183    
# Line 381  contains Line 188  contains
188         END IF         END IF
189      END DO      END DO
190    
     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 ')  
   
191    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
192    
193  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.121
  ViewVC Help
Powered by ViewVC 1.1.21