/[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 119 by guez, Wed Jan 7 14:34:57 2015 UTC revision 126 by guez, Fri Feb 6 18:33:15 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 fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22        use nr_util, only: pi, pi_d, twopi, 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)
31      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      REAL dzoom, step
32        real d_rlonv(iim)
33      LOGICAL, PARAMETER:: scal180 = .TRUE.      DOUBLE PRECISION xtild(0:2 * nmax)
34      ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35      ! une grille reguliere (grossismx = 1., taux=0., clon=0.) a      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      ! -180. degres. sinon scal180 = .FALSE.      DOUBLE PRECISION xzoom, fa, fb
37        INTEGER i, is2
38      REAL dzoom      DOUBLE PRECISION xmoy, fxm
     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      xzoom = clon * pi_d / 180d0
        decalx = 1.  
     else  
        decalx = 0.75  
     END IF  
45    
46      IF (dzoomx < 1.) THEN      test_grossismx: if (grossismx == 1.) then
47         dzoom = dzoomx * twopi_d         step = twopi / iim
     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  
48    
49      DO i = 0, nmax2         xprimm025(:iim) = step
50         xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2         xprimp025(:iim) = step
51      END DO         xprimv(:iim) = step
52           xprimu(:iim) = step
53    
54           rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55           rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57           rlonu(:iim) = rlonv(:iim) + 0.5 * step
58        else
59           dzoom = dzoomx * twopi_d
60           xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61    
62      DO i = nmax, nmax2         ! Compute fhyp:
63         fa = taux* (dzoom / 2. - xtild(i))         DO i = nmax, 2 * nmax
64         fb = xtild(i) * (pi_d - xtild(i))            fa = taux * (dzoom / 2. - xtild(i))
65              fb = xtild(i) * (pi_d - xtild(i))
66         IF (200.* fb < - fa) THEN  
67            fhyp(i) = - 1.            IF (200. * fb < - fa) THEN
68         ELSE IF (200. * fb < fa) THEN               fhyp(i) = - 1.
69            fhyp(i) = 1.            ELSE IF (200. * fb < fa) THEN
70         ELSE               fhyp(i) = 1.
           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  
71            ELSE            ELSE
72               fhyp(i) = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
73                    IF (200. * fb + fa < 1e-10) THEN
74                       fhyp(i) = - 1.
75                    ELSE IF (200. * fb - fa < 1e-10) THEN
76                       fhyp(i) = 1.
77                    END IF
78                 ELSE
79                    fhyp(i) = TANH(fa / fb)
80                 END IF
81            END IF            END IF
        END IF  
82    
83         IF (xtild(i) == 0.) fhyp(i) = 1.            IF (xtild(i) == 0.) fhyp(i) = 1.
84         IF (xtild(i) == pi_d) fhyp(i) = -1.            IF (xtild(i) == pi_d) fhyp(i) = -1.
85      END DO         END DO
86    
87      ! Calcul de beta         ! Calcul de beta
88    
89      ffdx = 0.         ffdx = 0.
90    
91      DO i = nmax + 1, nmax2         DO i = nmax + 1, 2 * nmax
92         xmoy = 0.5 * (xtild(i-1) + xtild(i))            xmoy = 0.5 * (xtild(i-1) + xtild(i))
93         fa = taux* (dzoom / 2. - xmoy)            fa = taux * (dzoom / 2. - xmoy)
94         fb = xmoy * (pi_d - xmoy)            fb = xmoy * (pi_d - xmoy)
95    
96         IF (200.* fb < - fa) THEN            IF (200. * fb < - fa) THEN
97            fxm = - 1.               fxm = - 1.
98         ELSE IF (200. * fb < fa) THEN            ELSE IF (200. * fb < fa) THEN
99            fxm = 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  
100            ELSE            ELSE
101               fxm = TANH(fa / fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
102                    IF (200. * fb + fa < 1e-10) THEN
103                       fxm = - 1.
104                    ELSE IF (200. * fb - fa < 1e-10) THEN
105                       fxm = 1.
106                    END IF
107                 ELSE
108                    fxm = TANH(fa / fb)
109                 END IF
110            END IF            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)  
111    
112      IF (2. * beta - grossismx <= 0.) THEN            IF (xmoy == 0.) fxm = 1.
113         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'            IF (xmoy == pi_d) fxm = -1.
        print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'  
        STOP 1  
     END IF  
   
     ! calcul de Xprimt  
   
     DO i = nmax, nmax2  
        Xprimt(i) = beta + (grossismx - beta) * fhyp(i)  
     END DO  
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
114    
115      ! Calcul de Xf            ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
116           END DO
     Xf(0) = - pi_d  
117    
118      DO i = nmax + 1, nmax2         print *, "ffdx = ", ffdx
119         xmoy = 0.5 * (xtild(i-1) + xtild(i))         beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
120         fa = taux* (dzoom / 2. - xmoy)         print *, "beta = ", beta
121         fb = xmoy * (pi_d - xmoy)  
122           IF (2. * beta - grossismx <= 0.) THEN
123         IF (200.* fb < - fa) THEN            print *, 'Bad choice of grossismx, taux, dzoomx.'
124            fxm = - 1.            print *, 'Decrease dzoomx or grossismx.'
125         ELSE IF (200. * fb < fa) THEN            STOP 1
           fxm = 1.  
        ELSE  
           fxm = TANH(fa / fb)  
126         END IF         END IF
127    
128         IF (xmoy == 0.) fxm = 1.         ! calcul de Xprimt
129         IF (xmoy == pi_d) fxm = -1.         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130         xxpr(i) = beta + (grossismx - beta) * fxm         xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131      END DO  
132           ! Calcul de Xf
133    
134           DO i = nmax + 1, 2 * nmax
135              xmoy = 0.5 * (xtild(i-1) + xtild(i))
136              fa = taux * (dzoom / 2. - xmoy)
137              fb = xmoy * (pi_d - xmoy)
138    
139              IF (200. * fb < - fa) THEN
140                 fxm = - 1.
141              ELSE IF (200. * fb < fa) THEN
142                 fxm = 1.
143              ELSE
144                 fxm = TANH(fa / fb)
145              END IF
146    
147      DO i = nmax + 1, nmax2            IF (xmoy == 0.) fxm = 1.
148         xxpr(nmax2-i + 1) = xxpr(i)            IF (xmoy == pi_d) fxm = -1.
149      END DO            xxpr(i) = beta + (grossismx - beta) * fxm
150           END DO
151    
152      DO i=1, nmax2         xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     END DO  
153    
154      ! xuv = 0. si calcul aux points scalaires         Xf(0) = - pi_d
     ! xuv = 0.5 si calcul aux points U  
155    
156      loop_ik: DO ik = 1, 4         DO i=1, 2 * nmax - 1
157         IF (ik == 1) THEN            Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
158            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  
159    
160         xo1 = 0.         Xf(2 * nmax) = pi_d
161    
162         ii1=1         call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonm025(:iim), &
163         ii2=iim              xprimm025(:iim), xuv = - 0.25d0)
164         IF (ik == 1.and.grossismx == 1.) THEN         call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonv(:iim), &
165            ii1 = 2              xprimv(:iim), xuv = 0d0)
166            ii2 = iim + 1         call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonu(:iim), &
167         END IF              xprimu(:iim), xuv = 0.5d0)
168           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonp025(:iim), &
169                xprimp025(:iim), xuv = 0.25d0)
170        end if test_grossismx
171    
172        is2 = 0
173    
174        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
175             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
176           IF (clon <= 0.) THEN
177              is2 = 1
178    
179         DO i = ii1, ii2            do while (rlonm025(is2) < - pi .and. is2 < iim)
180            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  
181            end do            end do
182    
183            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
184                 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  
185               STOP 1               STOP 1
186            end if            end if
187           ELSE
188              is2 = iim
189    
190            xxprim(i) = twopi_d / (REAL(iim) * Xprimin)            do while (rlonm025(is2) > pi .and. is2 > 1)
191            xvrai(i) = xi + xzoom               is2 = is2 - 1
192         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  
193    
194         DO i = 1, iim -1            if (rlonm025(is2) > pi) then
195            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, ')'  
196               STOP 1               STOP 1
197            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  
198         END IF         END IF
199      end DO loop_ik      END IF
   
     print *  
200    
201      DO i = 1, iim      call principal_cshift(is2, rlonm025, xprimm025)
202         xlon(i) = rlonv(i + 1) - rlonv(i)      call principal_cshift(is2, rlonv, xprimv)
203      END DO      call principal_cshift(is2, rlonu, xprimu)
204      champmin = 1e12      call principal_cshift(is2, rlonp025, xprimp025)
205      champmax = -1e12  
206      DO i = 1, iim      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207         champmin = MIN(champmin, xlon(i))      print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208         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  
209    
210      DO i = 1, iim + 1      DO i = 1, iim + 1
211         IF (rlonp025(i) < rlonv(i)) THEN         IF (rlonp025(i) < rlonv(i)) THEN
212            print *, ' Attention ! rlonp025 < rlonv', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
213              print *, "< rlonv(", i, ") = ", rlonv(i)
214            STOP 1            STOP 1
215         END IF         END IF
216    
217         IF (rlonv(i) < rlonm025(i)) THEN         IF (rlonv(i) < rlonm025(i)) THEN
218            print *, ' Attention ! rlonm025 > rlonv', i            print *, 'rlonv(', i, ') = ', rlonv(i)
219              print *, "< rlonm025(", i, ") = ", rlonm025(i)
220            STOP 1            STOP 1
221         END IF         END IF
222    
223         IF (rlonp025(i) > rlonu(i)) THEN         IF (rlonp025(i) > rlonu(i)) THEN
224            print *, ' Attention ! rlonp025 > rlonu', i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
225              print *, "> rlonu(", i, ") = ", rlonu(i)
226            STOP 1            STOP 1
227         END IF         END IF
228      END DO      END DO
229    
     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 ')  
   
230    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
231    
232  end module fxhyp_m  end module fxhyp_m
Legend:
Removed from v.119  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21