/[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 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 =
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      dzoom = dzoomx * twopi_d
46         decalx = 1.      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
     else  
        decalx = 0.75  
     END IF  
47    
48      IF (dzoomx < 1.) THEN      ! Compute fhyp:
49         dzoom = dzoomx * twopi_d      DO i = nmax, 2 * nmax
50      ELSE IF (dzoomx < 25.) THEN         fa = taux * (dzoom / 2. - xtild(i))
        print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &  
             // "L'augmenter et relancer."  
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
     END IF  
   
     print *, 'dzoom (rad):', dzoom  
   
     DO i = 0, nmax2  
        xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2  
     END DO  
   
     DO i = nmax, nmax2  
        fa = taux* (dzoom / 2. - xtild(i))  
51         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
52    
53         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
54            fhyp(i) = - 1.            fhyp(i) = - 1.
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
62                  fhyp(i) = 1.                  fhyp(i) = 1.
63               END IF               END IF
64            ELSE            ELSE
# Line 100  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)
81    
82         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
83            fxm = - 1.            fxm = - 1.
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
91                  fxm = 1.                  fxm = 1.
92               END IF               END IF
93            ELSE            ELSE
# Line 127  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  
   
     Xf(0) = - pi_d  
119    
120      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
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)
124    
125         IF (200.* fb < - fa) THEN         IF (200. * fb < - fa) THEN
126            fxm = - 1.            fxm = - 1.
127         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
128            fxm = 1.            fxm = 1.
# Line 167  contains Line 135  contains
135         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
136      END DO      END DO
137    
138      DO i = nmax + 1, nmax2      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
139         xxpr(nmax2-i + 1) = xxpr(i)  
140      END DO      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      ! xuv = 0. si calcul aux points scalaires      Xf(2 * nmax) = pi_d
     ! xuv = 0.5 si calcul aux points U  
147    
148      loop_ik: DO ik = 1, 4      IF (grossismx == 1.) THEN
149         IF (ik == 1) THEN         decalx = 1d0
150            xuv = -0.25      else
151         ELSE IF (ik == 2) THEN         decalx = 0.75d0
152            xuv = 0.      END IF
        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  
   
        xlon(iim + 1) = xlon(1) + twopi_d  
        xprimm(iim + 1) = xprimm(1)  
153    
154         DO i = 1, iim + 1      xzoom = clon * pi_d / 180d0
155            xvrai(i) = xlon(i)*180. / pi_d      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
156         END DO           xprimm025, xuv = - 0.25d0)
157        call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
158         IF (ik == 1) THEN           xuv = 0d0)
159            DO i = 1, iim + 1      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
160               rlonm025(i) = xlon(i)           xuv = 0.5d0)
161               xprimm025(i) = xprimm(i)      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
162            END DO           xprimp025, xuv = 0.25d0)
        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  
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    
183         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
184            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
185              print *, "> rlonu(", i, ") = ", rlonu(i)
186            STOP 1            STOP 1
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.119  
changed lines
  Added in v.123
  ViewVC Help
Powered by ViewVC 1.1.21