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

Diff of /trunk/Sources/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 146 by guez, Tue Jun 16 17:27:33 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
34      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax      DOUBLE PRECISION ffdx, beta
35      REAL dzoom      INTEGER i, is2
36      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      DOUBLE PRECISION xxpr(nmax - 1), 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 / 2. &
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 / 2. - xmoy), xmoy * (pi_d - xmoy))
   
     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  
67    
68         IF (xtild(i) == 0.) fhyp(i) = 1.         ! Calcul de beta
        IF (xtild(i) == pi_d) fhyp(i) = -1.  
     END DO  
69    
70      ! Calcul de beta         ffdx = 0.
71    
72      ffdx = 0.         DO i = 1, nmax
73              ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
74           END DO
75    
76      DO i = nmax + 1, nmax2         print *, "ffdx = ", ffdx
77         xmoy = 0.5 * (xtild(i-1) + xtild(i))         beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
78         fa = taux * (dzoom / 2. - xmoy)         print *, "beta = ", beta
79         fb = xmoy * (pi_d - xmoy)  
80           IF (2. * beta - grossismx <= 0.) 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  
           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  
84         END IF         END IF
85    
86         IF (xmoy == 0.) fxm = 1.         G = beta + (grossismx - beta) * fhyp
        IF (xmoy == pi_d) fxm = -1.  
87    
88         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ! Calcul de Xf
     END DO  
   
     beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)  
89    
90      IF (2. * beta - grossismx <= 0.) THEN         xxpr = beta + (grossismx - beta) * fxm(:nmax - 1)
91         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         Xf(0) = 0d0
        print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'  
        STOP 1  
     END IF  
92    
93      ! calcul de Xprimt         DO i = 1, nmax - 1
94              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
95      DO i = nmax, nmax2         END DO
        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  
   
     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  
           fxm = TANH(fa / fb)  
        END IF  
   
        IF (xmoy == 0.) fxm = 1.  
        IF (xmoy == pi_d) fxm = -1.  
        xxpr(i) = beta + (grossismx - beta) * fxm  
     END DO  
96    
97      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)         Xf(nmax) = pi_d
98    
99      DO i=1, nmax2         call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
100         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))              xuv = - 0.25d0)
101      END DO         call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
102           call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
103           call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
104                xuv = 0.25d0)
105        end if test_grossismx
106    
107      is2 = 0      is2 = 0
108    
109      loop_ik: DO ik = 1, 4      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
110         ! xuv = 0. si calcul aux points scalaires           .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
111         ! xuv = 0.5 si calcul aux points U         IF (clon <= 0.) THEN
112              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  
113    
114         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
115            Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)               is2 = is2 + 1
   
           it = nmax2  
           do while (xfi < xf(it) .and. it >= 1)  
              it = it - 1  
116            end do            end do
117    
118            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
119                 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  
120               STOP 1               STOP 1
121            end if            end if
122           ELSE
123              is2 = iim
124    
125            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
126            xvrai(i) = xi + xzoom               is2 = is2 - 1
127         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  
128    
129         DO i = 1, iim -1            if (rlonm025(is2) > pi) then
130            IF (xvrai(i + 1) < xvrai(i)) THEN               print *, 'Rlonm025 plus grand que pi !'
              print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'  
131               STOP 1               STOP 1
132            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  
133         END IF         END IF
134      end DO loop_ik      END IF
   
     print *  
135    
136      DO i = 1, iim      call principal_cshift(is2, rlonm025, xprimm025)
137         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonv, xprimv)
138      END DO      call principal_cshift(is2, rlonu, xprimu)
139      champmin = 1e12      call principal_cshift(is2, rlonp025, xprimp025)
140      champmax = -1e12  
141      DO i = 1, iim      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
142         champmin = MIN(champmin, xlon(i))      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
143         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  
144    
145        ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
146      DO i = 1, iim + 1      DO i = 1, iim + 1
147         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
148            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
149              print *, "< rlonv(", i, ") = ", rlonv(i)
150            STOP 1            STOP 1
151         END IF         END IF
152    
153         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
154            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
155              print *, "< rlonm025(", i, ") = ", rlonm025(i)
156            STOP 1            STOP 1
157         END IF         END IF
158    
# Line 381  contains Line 163  contains
163         END IF         END IF
164      END DO      END DO
165    
     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 ')  
   
166    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
167    
168  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.146
  ViewVC Help
Powered by ViewVC 1.1.21