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