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