/[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 145 by guez, Tue Jun 16 15:23:29 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:
   
     DOUBLE PRECISION champmin, champmax  
31      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
32      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      REAL dzoom, step
33        real d_rlonv(iim)
34      LOGICAL, PARAMETER:: scal180 = .TRUE.      DOUBLE PRECISION xtild(0:2 * nmax)
35      ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
36      ! une grille reguliere (grossismx = 1., taux=0., clon=0.) a      DOUBLE PRECISION Xf(0:2 * nmax)
37      ! -180. degres. sinon scal180 = .FALSE.      INTEGER i, is2
38        DOUBLE PRECISION, dimension(nmax + 1:2 * nmax):: xxpr, xmoy, fxm
     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  
39    
40      !----------------------------------------------------------------------      !----------------------------------------------------------------------
41    
42      my_eps = 1e-3      print *, "Call sequence information: fxhyp"
     xzoom = clon * pi_d / 180.  
43    
44      IF (grossismx == 1. .AND. scal180) THEN      test_grossismx: if (grossismx == 1.) then
45         decalx = 1.         step = twopi / iim
     else  
        decalx = 0.75  
     END IF  
46    
47      IF (dzoomx < 1.) THEN         xprimm025(:iim) = step
48           xprimp025(:iim) = step
49           xprimv(:iim) = step
50           xprimu(:iim) = step
51    
52           rlonv(:iim) = arth(- pi + clon, step, iim)
53           rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55           rlonu(:iim) = rlonv(:iim) + 0.5 * step
56        else test_grossismx
57         dzoom = dzoomx * twopi_d         dzoom = dzoomx * twopi_d
58      ELSE IF (dzoomx < 25.) THEN         xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &         forall (i = nmax + 1:2 * 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  
60    
61      DO i = 0, nmax2         ! Compute fhyp:
62         xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2         fhyp(nmax + 1:2 * nmax - 1) = tanh_cautious(taux * (dzoom / 2. &
63      END DO              - xtild(nmax + 1:2 * nmax - 1)), xtild(nmax + 1:2 * nmax - 1) &
64                * (pi_d - xtild(nmax + 1:2 * nmax - 1)))
65           fhyp(nmax) = 1d0
66           fhyp(2 * nmax) = -1d0
67    
68      DO i = nmax, nmax2         fxm = tanh_cautious(taux * (dzoom / 2. - xmoy), xmoy * (pi_d - xmoy))
        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  
69    
70         IF (xtild(i) == 0.) fhyp(i) = 1.         ! Calcul de beta
        IF (xtild(i) == pi_d) fhyp(i) = -1.  
     END DO  
71    
72      ! Calcul de beta         ffdx = 0.
73    
74      ffdx = 0.         DO i = nmax + 1, 2 * nmax
75              ffdx = ffdx + fxm(i) * (xtild(i) - xtild(i-1))
76           END DO
77    
78      DO i = nmax + 1, nmax2         print *, "ffdx = ", ffdx
79         xmoy = 0.5 * (xtild(i-1) + xtild(i))         beta = (pi_d - grossismx * ffdx) / (pi_d - ffdx)
80         fa = taux* (dzoom / 2. - xmoy)         print *, "beta = ", beta
81         fb = xmoy * (pi_d - xmoy)  
82           IF (2. * beta - grossismx <= 0.) THEN
83         IF (200.* fb < - fa) THEN            print *, 'Bad choice of grossismx, taux, dzoomx.'
84            fxm = - 1.            print *, 'Decrease dzoomx or grossismx.'
85         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  
86         END IF         END IF
87    
88         IF (xmoy == 0.) fxm = 1.         ! calcul de Xprimt
89         IF (xmoy == pi_d) fxm = -1.         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
90           xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
        ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))  
     END DO  
91    
92      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)         ! Calcul de Xf
93    
94      IF (2. * beta - grossismx <= 0.) THEN         xxpr = beta + (grossismx - beta) * fxm
95         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         Xf(nmax) = 0d0
        print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'  
        STOP 1  
     END IF  
96    
97      ! calcul de Xprimt         DO i = nmax + 1, 2 * nmax - 1
98              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
99      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  
   
     DO i = nmax + 1, nmax2  
        xxpr(nmax2-i + 1) = xxpr(i)  
     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  
   
     loop_ik: DO ik = 1, 4  
        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  
100    
101         xo1 = 0.         Xf(2 * nmax) = pi_d
102           xf(:nmax - 1) = - xf(2 * nmax:nmax + 1:- 1)
103    
104         ii1=1         call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
105         ii2=iim              xuv = - 0.25d0)
106         IF (ik == 1.and.grossismx == 1.) THEN         call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
107            ii1 = 2              xuv = 0d0)
108            ii2 = iim + 1         call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
109         END IF              xuv = 0.5d0)
110           call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
111                xuv = 0.25d0)
112        end if test_grossismx
113    
114        is2 = 0
115    
116        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
117             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
118           IF (clon <= 0.) THEN
119              is2 = 1
120    
121         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
122            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  
123            end do            end do
124    
125            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
126                 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  
127               STOP 1               STOP 1
128            end if            end if
129           ELSE
130              is2 = iim
131    
132            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
133            xvrai(i) = xi + xzoom               is2 = is2 - 1
134         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  
135    
136         DO i = 1, iim -1            if (rlonm025(is2) > pi) then
137            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, ')'  
138               STOP 1               STOP 1
139            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  
140         END IF         END IF
141      end DO loop_ik      END IF
142    
143      print *      call principal_cshift(is2, rlonm025, xprimm025)
144        call principal_cshift(is2, rlonv, xprimv)
145      DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
146         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonp025, xprimp025)
147      END DO  
148      champmin = 1e12      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
149      champmax = -1e12      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
150      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  
151    
152        ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
153      DO i = 1, iim + 1      DO i = 1, iim + 1
154         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
155            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
156              print *, "< rlonv(", i, ") = ", rlonv(i)
157            STOP 1            STOP 1
158         END IF         END IF
159    
160         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
161            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
162              print *, "< rlonm025(", i, ") = ", rlonm025(i)
163            STOP 1            STOP 1
164         END IF         END IF
165    
166         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
167            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
168              print *, "> rlonu(", i, ") = ", rlonu(i)
169            STOP 1            STOP 1
170         END IF         END IF
171      END DO      END DO
172    
     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 ')  
   
173    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
174    
175  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.119  
changed lines
  Added in v.145
  ViewVC Help
Powered by ViewVC 1.1.21