/[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 119 by guez, Wed Jan 7 14:34:57 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 =
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
34      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      DOUBLE PRECISION ffdx, beta
35        INTEGER i, is2
36      LOGICAL, PARAMETER:: scal180 = .TRUE.      DOUBLE PRECISION xxpr(nmax - 1), 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  
58    
59      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)         ! 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      IF (2. * beta - grossismx <= 0.) THEN         fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
        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  
67    
68      ! calcul de Xprimt         ! Calcul de beta
69    
70      DO i = nmax, nmax2         ffdx = 0.
        Xprimt(i) = beta + (grossismx - beta) * fhyp(i)  
     END DO  
71    
72      DO i = nmax + 1, nmax2         DO i = 1, nmax
73         Xprimt(nmax2 - i) = Xprimt(i)            ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
74      END DO         END DO
   
     ! Calcul de Xf  
   
     Xf(0) = - pi_d  
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  
           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  
   
     DO i = nmax + 1, nmax2  
        xxpr(nmax2-i + 1) = xxpr(i)  
     END DO  
87    
88      DO i=1, nmax2         ! Calcul de Xf
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     END DO  
89    
90      ! xuv = 0. si calcul aux points scalaires         xxpr = beta + (grossismx - beta) * fxm(:nmax - 1)
91      ! xuv = 0.5 si calcul aux points U         Xf(0) = 0d0
92    
93      loop_ik: DO ik = 1, 4         DO i = 1, nmax - 1
94         IF (ik == 1) THEN            Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
95            xuv = -0.25         END DO
        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  
96    
97         xo1 = 0.         Xf(nmax) = pi_d
98    
99         ii1=1         call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
100         ii2=iim              xuv = - 0.25d0)
101         IF (ik == 1.and.grossismx == 1.) THEN         call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
102            ii1 = 2         call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
103            ii2 = iim + 1         call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
104         END IF              xuv = 0.25d0)
105        end if test_grossismx
106    
107        is2 = 0
108    
109        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
110             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
111           IF (clon <= 0.) THEN
112              is2 = 1
113    
114         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
115            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  
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, xlon2  
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 *, 'Problème avec rlonu(', i + 1, &  
                   ') plus petit que rlonu(', i, ')'  
131               STOP 1               STOP 1
132            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  
133         END IF         END IF
134      end DO loop_ik      END IF
135    
136      print *      call principal_cshift(is2, rlonm025, xprimm025)
137        call principal_cshift(is2, rlonv, xprimv)
138      DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
139         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonp025, xprimp025)
140      END DO  
141      champmin = 1e12      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
142      champmax = -1e12      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
143      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  
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    
159         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
160            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
161              print *, "> rlonu(", i, ") = ", rlonu(i)
162            STOP 1            STOP 1
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.119  
changed lines
  Added in v.146
  ViewVC Help
Powered by ViewVC 1.1.21