/[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 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      USE dimens_m, ONLY: iim      ! Le premier point scalaire pour une grille regulière (grossismx =
18      use nr_util, only: pi_d, twopi_d      ! 1) avec clon = 0 est à - 180 degrés.
     use serre, only: clon, grossismx, dzoomx, taux  
   
     REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)  
     real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)  
   
     ! Local:  
   
     DOUBLE PRECISION champmin, champmax  
     real rlonm025(iim + 1), rlonp025(iim + 1)  
     INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax  
   
     LOGICAL, PARAMETER:: scal180 = .TRUE.  
     ! 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  
   
     !----------------------------------------------------------------------  
   
     my_eps = 1e-3  
     xzoom = clon * pi_d / 180.  
   
     IF (grossismx == 1. .AND. scal180) THEN  
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
   
     IF (dzoomx < 1.) THEN  
        dzoom = dzoomx * twopi_d  
     ELSE IF (dzoomx < 25.) THEN  
        print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &  
             // "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  
19    
20      DO i = nmax, nmax2      USE dimens_m, ONLY: iim
21         fa = taux* (dzoom / 2. - xtild(i))      use dynetat0_m, only: clon, grossismx, dzoomx, taux
22         fb = xtild(i) * (pi_d - xtild(i))      use invert_zoom_x_m, only: invert_zoom_x, nmax
23        use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24         IF (200.* fb < - fa) THEN      use principal_cshift_m, only: principal_cshift
25            fhyp(i) = - 1.      use tanh_cautious_m, only: tanh_cautious
        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  
26    
27      ffdx = 0.      REAL, intent(out):: xprimm025(:) ! (iim + 1)
28    
29      DO i = nmax + 1, nmax2      REAL, intent(out):: rlonv(:) ! (iim + 1)
30         xmoy = 0.5 * (xtild(i-1) + xtild(i))      ! longitudes of points of the "scalar" and "v" grid, in rad
        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  
31    
32         IF (xmoy == 0.) fxm = 1.      REAL, intent(out):: xprimv(:) ! (iim + 1)
33         IF (xmoy == pi_d) fxm = -1.      ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonv)
34    
35         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))      real, intent(out):: rlonu(:) ! (iim + 1)
36      END DO      ! longitudes of points of the "u" grid, in rad
37    
38      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)      real, intent(out):: xprimu(:) ! (iim + 1)
39        ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonu)
40    
41      IF (2. * beta - grossismx <= 0.) THEN      real, intent(out):: xprimp025(:) ! (iim + 1)
        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  
42    
43      ! calcul de Xprimt      ! Local:
44        real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
45        REAL delta, h
46        DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
47        DOUBLE PRECISION beta
48        INTEGER i, is2
49        DOUBLE PRECISION xmoy(nmax), fxm(nmax)
50    
51      DO i = nmax, nmax2      !----------------------------------------------------------------------
        Xprimt(i) = beta + (grossismx - beta) * fhyp(i)  
     END DO  
52    
53      DO i = nmax + 1, nmax2      print *, "Call sequence information: fxhyp"
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
54    
55      ! Calcul de Xf      if (grossismx == 1.) then
56           h = twopi / iim
57    
58      Xf(0) = - pi_d         xprimm025(:iim) = h
59           xprimp025(:iim) = h
60           xprimv(:iim) = h
61           xprimu(:iim) = h
62    
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      DO i = nmax + 1, nmax2         ! Compute fhyp:
73         xmoy = 0.5 * (xtild(i-1) + xtild(i))         fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
74         fa = taux* (dzoom / 2. - xmoy)              - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
75         fb = xmoy * (pi_d - xmoy)              * (pi_d - xtild(1:nmax - 1)))
76           fhyp(0) = 1d0
77         IF (200.* fb < - fa) THEN         fhyp(nmax) = -1d0
           fxm = - 1.  
        ELSE IF (200. * fb < fa) THEN  
           fxm = 1.  
        ELSE  
           fxm = TANH(fa / fb)  
        END IF  
78    
79         IF (xmoy == 0.) fxm = 1.         fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))
        IF (xmoy == pi_d) fxm = -1.  
        xxpr(i) = beta + (grossismx - beta) * fxm  
     END DO  
80    
81      DO i = nmax + 1, nmax2         ! Compute \int_0 ^{\tilde x} F:
        xxpr(nmax2-i + 1) = xxpr(i)  
     END DO  
82    
83      DO i=1, nmax2         ffdx(0) = 0d0
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     END DO  
84    
85      ! xuv = 0. si calcul aux points scalaires         DO i = 1, nmax
86      ! xuv = 0.5 si calcul aux points U            ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
87           END DO
88    
89      loop_ik: DO ik = 1, 4         print *, "ffdx(nmax) = ", ffdx(nmax)
90         IF (ik == 1) THEN         beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
91            xuv = -0.25         print *, "beta = ", beta
92         ELSE IF (ik == 2) THEN  
93            xuv = 0.         IF (2d0 * beta - grossismx <= 0d0) THEN
94         ELSE IF (ik == 3) THEN            print *, 'Bad choice of grossismx, taux, dzoomx.'
95            xuv = 0.50            print *, 'Decrease dzoomx or grossismx.'
96         ELSE IF (ik == 4) THEN            STOP 1
           xuv = 0.25  
97         END IF         END IF
98    
99         xo1 = 0.         G = beta + (grossismx - beta) * fhyp
100    
101         ii1=1         Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
102         ii2=iim              * ffdx(:nmax - 1)
103         IF (ik == 1.and.grossismx == 1.) THEN         Xf(nmax) = pi_d
104            ii1 = 2  
105            ii2 = iim + 1         call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
106         END IF              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
114    
115        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
116             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
117           IF (clon <= 0.) THEN
118              is2 = 1
119    
120         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
121            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  
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, xlon2  
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 *, 'Problème avec rlonu(', i + 1, &  
                   ') plus petit que rlonu(', i, ')'  
137               STOP 1               STOP 1
138            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  
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    
165         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
166            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
167              print *, "> rlonu(", i, ") = ", rlonu(i)
168            STOP 1            STOP 1
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.119  
changed lines
  Added in v.156
  ViewVC Help
Powered by ViewVC 1.1.21