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

Diff of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 120 by guez, Tue Jan 13 14:56:15 2015 UTC revision 124 by guez, Thu Feb 5 15:19:37 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 fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22        use nr_util, only: pi, pi_d, twopi_d, arth
23        use principal_cshift_m, only: principal_cshift
24      use serre, only: clon, grossismx, dzoomx, taux      use serre, only: clon, grossismx, dzoomx, taux
25    
26      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
28    
29      ! Local:      ! Local:
   
     DOUBLE PRECISION champmin, champmax  
30      real rlonm025(iim + 1), rlonp025(iim + 1)      real rlonm025(iim + 1), rlonp025(iim + 1)
     INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax  
31      REAL dzoom      REAL dzoom
32      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv      real d_rlonv(iim)
33      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION xtild(0:2 * nmax)
34      DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35      DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)      DOUBLE PRECISION xzoom, fa, fb
37      DOUBLE PRECISION my_eps, xzoom, fa, fb      INTEGER i, is2
38      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2      DOUBLE PRECISION xmoy, fxm
     INTEGER i, it, ik, iter, ii, idif, ii1, ii2  
     DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin  
39      DOUBLE PRECISION decalx      DOUBLE PRECISION decalx
     INTEGER is2  
40    
41      !----------------------------------------------------------------------      !----------------------------------------------------------------------
42    
43      print *, "Call sequence information: fxhyp"      print *, "Call sequence information: fxhyp"
44    
45      my_eps = 1e-3      dzoom = dzoomx * twopi_d
46      xzoom = clon * pi_d / 180.      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
   
     IF (grossismx == 1.) THEN  
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
   
     IF (dzoomx < 1.) THEN  
        dzoom = dzoomx * twopi_d  
     ELSE IF (dzoomx < 25.) THEN  
        print *, "dzoomx pour fxhyp est trop petit."  
        STOP 1  
     ELSE  
        dzoom = dzoomx * pi_d / 180.  
     END IF  
   
     print *, 'dzoom (rad):', dzoom  
47    
48      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)      ! Compute fhyp:
49        DO i = nmax, 2 * nmax
     DO i = nmax, nmax2  
50         fa = taux * (dzoom / 2. - xtild(i))         fa = taux * (dzoom / 2. - xtild(i))
51         fb = xtild(i) * (pi_d - xtild(i))         fb = xtild(i) * (pi_d - xtild(i))
52    
# Line 78  contains Line 55  contains
55         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
56            fhyp(i) = 1.            fhyp(i) = 1.
57         ELSE         ELSE
58            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
59               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
60                  fhyp(i) = - 1.                  fhyp(i) = - 1.
61               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 97  contains Line 74  contains
74    
75      ffdx = 0.      ffdx = 0.
76    
77      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
78         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
79         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
80         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 107  contains Line 84  contains
84         ELSE IF (200. * fb < fa) THEN         ELSE IF (200. * fb < fa) THEN
85            fxm = 1.            fxm = 1.
86         ELSE         ELSE
87            IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
88               IF (200. * fb + fa < 1e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
89                  fxm = - 1.                  fxm = - 1.
90               ELSE IF (200. * fb - fa < 1e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
# Line 124  contains Line 101  contains
101         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
102      END DO      END DO
103    
104        print *, "ffdx = ", ffdx
105      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
106        print *, "beta = ", beta
107    
108      IF (2. * beta - grossismx <= 0.) THEN      IF (2. * beta - grossismx <= 0.) THEN
109         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         print *, 'Bad choice of grossismx, taux, dzoomx.'
110         print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'         print *, 'Decrease dzoomx or grossismx.'
111         STOP 1         STOP 1
112      END IF      END IF
113    
114      ! calcul de Xprimt      ! calcul de Xprimt
115        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
116        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
117    
118      DO i = nmax, nmax2      ! Calcul de Xf
        Xprimt(i) = beta + (grossismx - beta) * fhyp(i)  
     END DO  
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
   
     ! Calcul de Xf  
119    
120      Xf(0) = - pi_d      DO i = nmax + 1, 2 * nmax
   
     DO i = nmax + 1, nmax2  
121         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
122         fa = taux * (dzoom / 2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
123         fb = xmoy * (pi_d - xmoy)         fb = xmoy * (pi_d - xmoy)
# Line 164  contains Line 135  contains
135         xxpr(i) = beta + (grossismx - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
136      END DO      END DO
137    
138      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
139    
140        Xf(0) = - pi_d
141    
142      DO i=1, nmax2      DO i=1, 2 * nmax - 1
143         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
144      END DO      END DO
145    
146      is2 = 0      Xf(2 * nmax) = pi_d
147    
148      loop_ik: DO ik = 1, 4      IF (grossismx == 1.) THEN
149         ! xuv = 0. si calcul aux points scalaires         decalx = 1d0
150         ! xuv = 0.5 si calcul aux points U      else
151           decalx = 0.75d0
152         IF (ik == 1) THEN      END IF
           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  
153    
154         xo1 = 0.      xzoom = clon * pi_d / 180d0
155        call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025(:iim), &
156             xprimm025(:iim), xuv = - 0.25d0)
157        call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv(:iim), &
158             xprimv(:iim), xuv = 0d0)
159        call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu(:iim), &
160             xprimu(:iim), xuv = 0.5d0)
161        call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025(:iim), &
162             xprimp025(:iim), xuv = 0.25d0)
163    
164         IF (ik == 1 .and. grossismx == 1.) THEN      is2 = 0
           ii1 = 2  
           ii2 = iim + 1  
        else  
           ii1=1  
           ii2=iim  
        END IF  
165    
166         DO i = ii1, ii2      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
167            Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)           .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
168           IF (clon <= 0.) THEN
169              is2 = 1
170    
171            it = nmax2            do while (rlonm025(is2) < - pi .and. is2 < iim)
172            do while (xfi < xf(it) .and. it >= 1)               is2 = is2 + 1
              it = it - 1  
173            end do            end do
174    
175            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
176                 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  
177               STOP 1               STOP 1
178            end if            end if
179           ELSE
180              is2 = iim
181    
182            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
183            xvrai(i) = xi + xzoom               is2 = is2 - 1
184         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  
185    
186         DO i = 1, iim            if (rlonm025(is2) > pi) then
187            xlon(i) = xvrai(i)               print *, 'Rlonm025 plus grand que pi !'
           xprimm(i) = xxprim(i)  
        END DO  
   
        DO i = 1, iim -1  
           IF (xvrai(i + 1) < xvrai(i)) THEN  
              print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'  
188               STOP 1               STOP 1
189            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  
190         END IF         END IF
191      end DO loop_ik      END IF
   
     print *  
192    
193      DO i = 1, iim      call principal_cshift(is2, rlonm025, xprimm025)
194         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonv, xprimv)
195      END DO      call principal_cshift(is2, rlonu, xprimu)
196      champmin = 1e12      call principal_cshift(is2, rlonp025, xprimp025)
197      champmax = -1e12  
198      DO i = 1, iim      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
199         champmin = MIN(champmin, xlon(i))      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
200         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  
201    
202      DO i = 1, iim + 1      DO i = 1, iim + 1
203         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
204            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
205              print *, "< rlonv(", i, ") = ", rlonv(i)
206            STOP 1            STOP 1
207         END IF         END IF
208    
209         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
210            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
211              print *, "< rlonm025(", i, ") = ", rlonm025(i)
212            STOP 1            STOP 1
213         END IF         END IF
214    
# Line 381  contains Line 219  contains
219         END IF         END IF
220      END DO      END DO
221    
     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 ')  
   
222    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
223    
224  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.120  
changed lines
  Added in v.124
  ViewVC Help
Powered by ViewVC 1.1.21