/[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 123 by guez, Thu Feb 5 12:41:08 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      dzoom = dzoomx * twopi_d
46      xzoom = clon * pi_d / 180.      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
47    
48      IF (grossismx == 1.) THEN      ! Compute fhyp:
49         decalx = 1.      DO i = nmax, 2 * nmax
     else  
        decalx = 0.75  
     END IF  
   
     IF (dzoomx < 1.) THEN  
        dzoom = dzoomx * twopi_d  
     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  
   
     xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)  
   
     DO i = nmax, nmax2  
50         fa = taux * (dzoom / 2. - xtild(i))         fa = taux * (dzoom / 2. - xtild(i))
51         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
52    
# Line 78  contains Line 55  contains
55         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
56            fhyp(i) = 1.            fhyp(i) = 1.
57         ELSE         ELSE
58            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
59               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
60                  fhyp(i) = - 1.                  fhyp(i) = - 1.
61               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 97  contains Line 74  contains
74    
75      ffdx = 0.      ffdx = 0.
76    
77      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
78         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
79         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
80         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 107  contains Line 84  contains
84         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
85            fxm = 1.            fxm = 1.
86         ELSE         ELSE
87            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
88               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
89                  fxm = - 1.                  fxm = - 1.
90               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 124  contains Line 101  contains
101         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
102      END DO      END DO
103    
104        print *, "ffdx = ", ffdx
105      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
106        print *, "beta = ", beta
107    
108      IF (2. * beta - grossismx <= 0.) THEN      IF (2. * beta - grossismx <= 0.) THEN
109         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         print *, 'Bad choice of grossismx, taux, dzoomx.'
110         print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'         print *, 'Decrease dzoomx or grossismx.'
111         STOP 1         STOP 1
112      END IF      END IF
113    
114      ! calcul de Xprimt      ! calcul de Xprimt
115        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
116        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
117    
118      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  
119    
120      Xf(0) = - pi_d      DO i = nmax + 1, 2 * nmax
   
     DO i = nmax + 1, nmax2  
121         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
122         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
123         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 164  contains Line 135  contains
135         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
136      END DO      END DO
137    
138      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
139    
140        Xf(0) = - pi_d
141    
142      DO i=1, nmax2      DO i=1, 2 * nmax - 1
143         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
144      END DO      END DO
145    
146      is2 = 0      Xf(2 * nmax) = pi_d
   
     loop_ik: DO ik = 1, 4  
        ! xuv = 0. si calcul aux points scalaires  
        ! xuv = 0.5 si calcul aux points U  
   
        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.  
   
        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)  
147    
148            IF (it == nmax2) THEN      IF (grossismx == 1.) THEN
149               it = nmax2 -1         decalx = 1d0
150               Xf(it + 1) = pi_d      else
151            END IF         decalx = 0.75d0
152        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  
153    
154         xlon(iim + 1) = xlon(1) + twopi_d      xzoom = clon * pi_d / 180d0
155         xprimm(iim + 1) = xprimm(1)      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
156             xprimm025, xuv = - 0.25d0)
157         DO i = 1, iim + 1      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
158            xvrai(i) = xlon(i) * 180. / pi_d           xuv = 0d0)
159         END DO      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
160             xuv = 0.5d0)
161         IF (ik == 1) THEN      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
162            DO i = 1, iim + 1           xprimp025, xuv = 0.25d0)
              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  
163    
164      print *      print *
165    
166      DO i = 1, iim      forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
167         xlon(i) = rlonv(i + 1) - rlonv(i)      print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "degrees"
168      END DO      print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "degrees"
     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  
169    
170      DO i = 1, iim + 1      DO i = 1, iim + 1
171         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
172            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
173              print *, "< rlonv(", i, ") = ", rlonv(i)
174            STOP 1            STOP 1
175         END IF         END IF
176    
177         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
178            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
179              print *, "< rlonm025(", i, ") = ", rlonm025(i)
180            STOP 1            STOP 1
181         END IF         END IF
182    
# Line 381  contains Line 187  contains
187         END IF         END IF
188      END DO      END DO
189    
     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 ')  
   
190    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
191    
192  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.123
  ViewVC Help
Powered by ViewVC 1.1.21