/[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 119 by guez, Wed Jan 7 14:34:57 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 =
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 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        INTEGER i, is2
36      LOGICAL, PARAMETER:: scal180 = .TRUE.      DOUBLE PRECISION xmoy(nmax), fxm(nmax)
     ! 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.  
   
     REAL dzoom  
     DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv  
     DOUBLE PRECISION xtild(0:nmax2)  
     DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)  
     DOUBLE PRECISION Xf(0:nmax2), xxpr(0: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, xlon2, fxm, Xprimin  
     DOUBLE PRECISION decalx  
     INTEGER, save:: is2  
37    
38      !----------------------------------------------------------------------      !----------------------------------------------------------------------
39    
40      my_eps = 1e-3      print *, "Call sequence information: fxhyp"
     xzoom = clon * pi_d / 180.  
41    
42      IF (grossismx == 1. .AND. scal180) THEN      test_grossismx: if (grossismx == 1.) then
43         decalx = 1.         step = twopi / iim
     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 *, "Le paramètre dzoomx pour fxhyp est trop petit. " &         forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
             // "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))  
        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  
   
     ! Calcul de beta  
   
     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  
   
     beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)  
   
     IF (2. * beta - grossismx <= 0.) THEN  
        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  
58    
59      DO i = nmax, nmax2         ! Compute fhyp:
60         Xprimt(i) = beta + (grossismx - beta) * fhyp(i)         fhyp(1:nmax - 1) = tanh_cautious(taux * (dzoom / 2d0 &
61      END DO              - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
62                * (pi_d - xtild(1:nmax - 1)))
63      DO i = nmax + 1, nmax2         fhyp(0) = 1d0
64         Xprimt(nmax2 - i) = Xprimt(i)         fhyp(nmax) = -1d0
     END DO  
   
     ! Calcul de Xf  
65    
66      Xf(0) = - pi_d         fxm = tanh_cautious(taux * (dzoom / 2d0 - xmoy), xmoy * (pi_d - xmoy))
67    
68      DO i = nmax + 1, nmax2         ! Compute \int_0 ^{\tilde x} F:
        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  
69    
70         IF (xmoy == 0.) fxm = 1.         ffdx(0) = 0d0
        IF (xmoy == pi_d) fxm = -1.  
        xxpr(i) = beta + (grossismx - beta) * fxm  
     END DO  
71    
72      DO i = nmax + 1, nmax2         DO i = 1, nmax
73         xxpr(nmax2-i + 1) = xxpr(i)            ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
74      END DO         END DO
   
     DO i=1, nmax2  
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     END DO  
   
     ! xuv = 0. si calcul aux points scalaires  
     ! xuv = 0.5 si calcul aux points U  
75    
76      loop_ik: DO ik = 1, 4         print *, "ffdx(nmax) = ", ffdx(nmax)
77         IF (ik == 1) THEN         beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
78            xuv = -0.25         print *, "beta = ", beta
79         ELSE IF (ik == 2) THEN  
80            xuv = 0.         IF (2d0 * beta - grossismx <= 0d0) THEN
81         ELSE IF (ik == 3) THEN            print *, 'Bad choice of grossismx, taux, dzoomx.'
82            xuv = 0.50            print *, 'Decrease dzoomx or grossismx.'
83         ELSE IF (ik == 4) THEN            STOP 1
           xuv = 0.25  
84         END IF         END IF
85    
86         xo1 = 0.         G = beta + (grossismx - beta) * fhyp
87    
88         ii1=1         Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
89         ii2=iim              * ffdx(:nmax - 1)
90         IF (ik == 1.and.grossismx == 1.) THEN         Xf(nmax) = pi_d
91            ii1 = 2  
92            ii2 = iim + 1         call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
93         END IF              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
101    
102        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
103             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
104           IF (clon <= 0.) THEN
105              is2 = 1
106    
107         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
108            xlon2 = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)               is2 = is2 + 1
           Xfi = xlon2  
   
           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, xlon2  
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 *, 'Problème avec rlonu(', i + 1, &  
                   ') plus petit que rlonu(', i, ')'  
124               STOP 1               STOP 1
125            END IF            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)  
   
        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  
           DO i = 1, iim + 1  
              rlonp025(i) = xlon(i)  
              xprimp025(i) = xprimm(i)  
           END DO  
126         END IF         END IF
127      end DO loop_ik      END IF
128    
129      print *      call principal_cshift(is2, rlonm025, xprimm025)
130        call principal_cshift(is2, rlonv, xprimv)
131      DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
132         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonp025, xprimp025)
133      END DO  
134      champmin = 1e12      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
135      champmax = -1e12      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
136      DO i = 1, iim      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
        champmin = MIN(champmin, xlon(i))  
        champmax = MAX(champmax, xlon(i))  
     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    
152         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
153            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
154              print *, "> rlonu(", i, ") = ", rlonu(i)
155            STOP 1            STOP 1
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.119  
changed lines
  Added in v.147
  ViewVC Help
Powered by ViewVC 1.1.21