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

Diff of /trunk/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 147 by guez, Wed Jun 17 14:20:14 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        use tanh_cautious_m, only: tanh_cautious
26    
27      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
28      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
29    
30      ! Local:      ! Local:
31        real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
32      DOUBLE PRECISION champmin, champmax      REAL dzoom, step
33      real rlonm025(iim + 1), rlonp025(iim + 1)      DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
34      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax      DOUBLE PRECISION beta
35      REAL dzoom      INTEGER i, is2
36      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      DOUBLE PRECISION xmoy(nmax), fxm(nmax)
     DOUBLE PRECISION xtild(0:nmax2)  
     DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)  
     DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)  
     DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)  
     DOUBLE PRECISION my_eps, xzoom, fa, fb  
     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  
37    
38      !----------------------------------------------------------------------      !----------------------------------------------------------------------
39    
40      print *, "Call sequence information: fxhyp"      print *, "Call sequence information: fxhyp"
41    
42      my_eps = 1e-3      test_grossismx: if (grossismx == 1.) then
43      xzoom = clon * pi_d / 180.         step = twopi / iim
   
     IF (grossismx == 1.) THEN  
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
44    
45      IF (dzoomx < 1.) THEN         xprimm025(:iim) = step
46           xprimp025(:iim) = step
47           xprimv(:iim) = step
48           xprimu(:iim) = step
49    
50           rlonv(:iim) = arth(- pi + clon, step, iim)
51           rlonm025(:iim) = rlonv(:iim) - 0.25 * step
52           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
53           rlonu(:iim) = rlonv(:iim) + 0.5 * step
54        else test_grossismx
55         dzoom = dzoomx * twopi_d         dzoom = dzoomx * twopi_d
56      ELSE IF (dzoomx < 25.) THEN         xtild = arth(0d0, pi_d / nmax, nmax + 1)
57         print *, "dzoomx pour fxhyp est trop petit."         forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
     END IF  
58    
59      print *, 'dzoom (rad):', dzoom         ! Compute fhyp:
60           fhyp(1:nmax - 1) = tanh_cautious(taux * (dzoom / 2d0 &
61                - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
62                * (pi_d - xtild(1:nmax - 1)))
63           fhyp(0) = 1d0
64           fhyp(nmax) = -1d0
65    
66      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)         fxm = tanh_cautious(taux * (dzoom / 2d0 - xmoy), xmoy * (pi_d - xmoy))
67    
68      DO i = nmax, nmax2         ! Compute \int_0 ^{\tilde x} F:
        fa = taux * (dzoom / 2. - xtild(i))  
        fb = xtild(i) * (pi_d - xtild(i))  
   
        IF (200. * fb < - fa) THEN  
           fhyp(i) = - 1.  
        ELSE IF (200. * fb < fa) THEN  
           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  
           ELSE  
              fhyp(i) = TANH(fa / fb)  
           END IF  
        END IF  
   
        IF (xtild(i) == 0.) fhyp(i) = 1.  
        IF (xtild(i) == pi_d) fhyp(i) = -1.  
     END DO  
69    
70      ! Calcul de beta         ffdx(0) = 0d0
   
     ffdx = 0.  
   
     DO i = nmax + 1, nmax2  
        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  
           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  
           ELSE  
              fxm = TANH(fa / fb)  
           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  
71    
72      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)         DO i = 1, nmax
73              ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
74      IF (2. * beta - grossismx <= 0.) THEN         END DO
        print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'  
        print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'  
        STOP 1  
     END IF  
   
     ! calcul de Xprimt  
   
     DO i = nmax, nmax2  
        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  
75    
76      DO i = nmax + 1, nmax2         print *, "ffdx(nmax) = ", ffdx(nmax)
77         xmoy = 0.5 * (xtild(i-1) + xtild(i))         beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
78         fa = taux * (dzoom / 2. - xmoy)         print *, "beta = ", beta
79         fb = xmoy * (pi_d - xmoy)  
80           IF (2d0 * beta - grossismx <= 0d0) THEN
81         IF (200. * fb < - fa) THEN            print *, 'Bad choice of grossismx, taux, dzoomx.'
82            fxm = - 1.            print *, 'Decrease dzoomx or grossismx.'
83         ELSE IF (200. * fb < fa) THEN            STOP 1
           fxm = 1.  
        ELSE  
           fxm = TANH(fa / fb)  
84         END IF         END IF
85    
86         IF (xmoy == 0.) fxm = 1.         G = beta + (grossismx - beta) * fhyp
        IF (xmoy == pi_d) fxm = -1.  
        xxpr(i) = beta + (grossismx - beta) * fxm  
     END DO  
87    
88      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)         Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
89                * ffdx(:nmax - 1)
90      DO i=1, nmax2         Xf(nmax) = pi_d
91         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
92      END DO         call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
93                xuv = - 0.25d0)
94           call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
95           call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
96           call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
97                xuv = 0.25d0)
98        end if test_grossismx
99    
100      is2 = 0      is2 = 0
101    
102      loop_ik: DO ik = 1, 4      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
103         ! xuv = 0. si calcul aux points scalaires           .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
104         ! xuv = 0.5 si calcul aux points U         IF (clon <= 0.) THEN
105              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.  
106    
107         IF (ik == 1 .and. grossismx == 1.) THEN            do while (rlonm025(is2) < - pi .and. is2 < iim)
108            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  
109            end do            end do
110    
111            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
112                 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  
113               STOP 1               STOP 1
114            end if            end if
115           ELSE
116              is2 = iim
117    
118            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
119            xvrai(i) = xi + xzoom               is2 = is2 - 1
120         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  
121    
122         DO i = 1, iim -1            if (rlonm025(is2) > pi) then
123            IF (xvrai(i + 1) < xvrai(i)) THEN               print *, 'Rlonm025 plus grand que pi !'
              print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'  
124               STOP 1               STOP 1
125            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  
126         END IF         END IF
127      end DO loop_ik      END IF
   
     print *  
128    
129      DO i = 1, iim      call principal_cshift(is2, rlonm025, xprimm025)
130         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonv, xprimv)
131      END DO      call principal_cshift(is2, rlonu, xprimu)
132      champmin = 1e12      call principal_cshift(is2, rlonp025, xprimp025)
133      champmax = -1e12  
134      DO i = 1, iim      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
135         champmin = MIN(champmin, xlon(i))      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
136         champmax = MAX(champmax, xlon(i))      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
     END DO  
     champmin = champmin * 180. / pi_d  
     champmax = champmax * 180. / pi_d  
137    
138        ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
139      DO i = 1, iim + 1      DO i = 1, iim + 1
140         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
141            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
142              print *, "< rlonv(", i, ") = ", rlonv(i)
143            STOP 1            STOP 1
144         END IF         END IF
145    
146         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
147            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
148              print *, "< rlonm025(", i, ") = ", rlonm025(i)
149            STOP 1            STOP 1
150         END IF         END IF
151    
# Line 381  contains Line 156  contains
156         END IF         END IF
157      END DO      END DO
158    
     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 ')  
   
159    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
160    
161  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.147
  ViewVC Help
Powered by ViewVC 1.1.21