/[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 119 by guez, Wed Jan 7 14:34:57 2015 UTC revision 122 by guez, Tue Feb 3 19:30:48 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 =
18        ! 1., taux=0., clon=0.) est à - 180 degrés.
19    
20      USE dimens_m, ONLY: iim      USE dimens_m, ONLY: iim
21      use nr_util, only: pi_d, twopi_d      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22        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    
25      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
# Line 23  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  
   
     LOGICAL, PARAMETER:: scal180 = .TRUE.  
     ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour  
     ! une grille reguliere (grossismx = 1., taux=0., clon=0.) a  
     ! -180. degres. sinon scal180 = .FALSE.  
   
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(0: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(0: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, xlon2, fxm, Xprimin  
39      DOUBLE PRECISION decalx      DOUBLE PRECISION decalx
     INTEGER, save:: is2  
40    
41      !----------------------------------------------------------------------      !----------------------------------------------------------------------
42    
43      my_eps = 1e-3      print *, "Call sequence information: fxhyp"
     xzoom = clon * pi_d / 180.  
44    
45      IF (grossismx == 1. .AND. scal180) THEN      xzoom = clon * pi_d / 180d0
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
46    
47      IF (dzoomx < 1.) THEN      IF (grossismx == 1.) THEN
48         dzoom = dzoomx * twopi_d         decalx = 1d0
49      ELSE IF (dzoomx < 25.) THEN      else
50         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &         decalx = 0.75d0
             // "L'augmenter et relancer."  
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
51      END IF      END IF
52    
53      print *, 'dzoom (rad):', dzoom      dzoom = dzoomx * twopi_d
54        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
     DO i = 0, nmax2  
        xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2  
     END DO  
55    
56      DO i = nmax, nmax2      ! Compute fhyp:
57         fa = taux* (dzoom / 2. - xtild(i))      DO i = nmax, 2 * nmax
58           fa = taux * (dzoom / 2. - xtild(i))
59         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
60    
61         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
62            fhyp(i) = - 1.            fhyp(i) = - 1.
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
70                  fhyp(i) = 1.                  fhyp(i) = 1.
71               END IF               END IF
72            ELSE            ELSE
# Line 100  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)
89    
90         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
91            fxm = - 1.            fxm = - 1.
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
99                  fxm = 1.                  fxm = 1.
100               END IF               END IF
101            ELSE            ELSE
# Line 127  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 *, 'Bad choice of grossismx, taux, dzoomx.'
118         print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'         print *, 'Decrease dzoomx or grossismx.'
119         STOP 1         STOP 1
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)
132    
133         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
134            fxm = - 1.            fxm = - 1.
135         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
136            fxm = 1.            fxm = 1.
# Line 167  contains Line 143  contains
143         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
144      END DO      END DO
145    
146      DO i = nmax + 1, nmax2      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
147         xxpr(nmax2-i + 1) = xxpr(i)  
148      END DO      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      ! xuv = 0. si calcul aux points scalaires      Xf(2 * nmax) = pi_d
     ! xuv = 0.5 si calcul aux points U  
   
     loop_ik: DO ik = 1, 4  
        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.  
   
        ii1=1  
        ii2=iim  
        IF (ik == 1.and.grossismx == 1.) THEN  
           ii1 = 2  
           ii2 = iim + 1  
        END IF  
   
        DO i = ii1, ii2  
           xlon2 = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)  
           Xfi = xlon2  
   
           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, xlon2  
              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 *, 'Problème avec rlonu(', i + 1, &  
                   ') plus petit que rlonu(', i, ')'  
              STOP 1  
           END IF  
        END DO  
   
        ! Réorganisation des longitudes pour les avoir entre - pi et pi  
   
        champmin = 1e12  
        champmax = -1e12  
        DO i = 1, iim  
           champmin = MIN(champmin, xvrai(i))  
           champmax = MAX(champmax, xvrai(i))  
        END DO  
   
        IF (.not. (champmin >= -pi_d - 0.1 .and. champmax <= pi_d + 0.1)) THEN  
           print *, 'Reorganisation des longitudes pour 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.NE. 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  
   
        ! Fin de la reorganisation  
155    
156         xlon(iim + 1) = xlon(1) + twopi_d      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
157         xprimm(iim + 1) = xprimm(1)           xprimm025, xuv = - 0.25d0)
158        call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
159         DO i = 1, iim + 1           xuv = 0d0)
160            xvrai(i) = xlon(i)*180. / pi_d      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
161         END DO           xuv = 0.5d0)
162        call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
163         IF (ik == 1) THEN           xprimp025, xuv = 0.25d0)
           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  
           DO i = 1, iim + 1  
              rlonp025(i) = xlon(i)  
              xprimp025(i) = xprimm(i)  
           END DO  
        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, "degrees"
169      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  
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    
184         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
185            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
186              print *, "> rlonu(", i, ") = ", rlonu(i)
187            STOP 1            STOP 1
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.119  
changed lines
  Added in v.122
  ViewVC Help
Powered by ViewVC 1.1.21