/[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 91 by guez, Wed Mar 26 17:18:58 2014 UTC trunk/Sources/dyn3d/fxhyp.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC
# Line 4  module fxhyp_m Line 4  module fxhyp_m
4    
5  contains  contains
6    
7    SUBROUTINE fxhyp(xzoomdeg, grossism, dzooma, tau, rlonm025, xprimm025, &    SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
        rlonv, xprimv, rlonu, xprimu, rlonp025, xprimp025, champmin, champmax)  
8    
9      ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32      ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10      ! Author: P. Le Van      ! 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) à tangente hyperbolique.      ! une fonction f(x) à dérivée tangente hyperbolique.
14    
15      ! On doit avoir grossism \times dzoom < pi (radians), en longitude.      ! Il vaut mieux avoir : grossismx \times dzoom < pi
16    
17      USE dimens_m, ONLY: iim      ! Le premier point scalaire pour une grille regulière (grossismx =
18      USE paramet_m, ONLY: iip1      ! 1., taux=0., clon=0.) est à - 180 degrés.
   
     REAL, intent(in):: xzoomdeg  
   
     REAL, intent(in):: grossism  
     ! grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)  
   
     REAL, intent(in):: dzooma ! distance totale de la zone du zoom  
19    
20      REAL, intent(in):: tau      USE dimens_m, ONLY: iim
21      ! raideur de la transition de l'intérieur à l'extérieur du zoom      use invert_zoom_x_m, only: invert_zoom_x, nmax
22        use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23      ! arguments de sortie      use principal_cshift_m, only: principal_cshift
24        use serre, only: clon, grossismx, dzoomx, taux
     REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv  
     real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025  
25    
26      DOUBLE PRECISION, intent(out):: champmin, champmax      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
27        real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
28    
29      ! Local:      ! Local:
30        real rlonm025(iim + 1), rlonp025(iim + 1)
31        REAL dzoom, step
32        real d_rlonv(iim)
33        DOUBLE PRECISION xtild(0:2 * nmax)
34        DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35        DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36        DOUBLE PRECISION fa, fb
37        INTEGER i, is2
38        DOUBLE PRECISION xmoy, fxm
39    
40      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 (grossism = 1., tau=0., clon=0.) a  
     ! -180. degres. sinon scal180 = .FALSE.  
   
     REAL dzoom  
     DOUBLE PRECISION xlon(iip1), xprimm(iip1), 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(iip1), xxprim(iip1)  
     DOUBLE PRECISION pi, depi, epsilon, 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 is2  
     SAVE is2  
41    
42      DOUBLE PRECISION heavyside      print *, "Call sequence information: fxhyp"
43    
44      !----------------------------------------------------------------------      test_grossismx: if (grossismx == 1.) then
45           step = twopi / iim
46    
47      pi = 2. * ASIN(1.)         xprimm025(:iim) = step
48      depi = 2. * pi         xprimp025(:iim) = step
49      epsilon = 1.e-3         xprimv(:iim) = step
50      xzoom = xzoomdeg * pi/180.         xprimu(:iim) = step
51    
52      decalx = .75         rlonv(:iim) = arth(- pi + clon, step, iim)
53      IF (grossism == 1. .AND. scal180) THEN         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54         decalx = 1.         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55      ENDIF         rlonu(:iim) = rlonv(:iim) + 0.5 * step
56        else test_grossismx
57      print *, 'FXHYP scal180, decalx', scal180, decalx         dzoom = dzoomx * twopi_d
58           xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59      IF (dzooma.LT.1.) THEN  
60         dzoom = dzooma * depi         ! Compute fhyp:
61      ELSEIF (dzooma.LT. 25.) THEN         DO i = nmax, 2 * nmax
62         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &            fa = taux * (dzoom / 2. - xtild(i))
63              // "L'augmenter et relancer."            fb = xtild(i) * (pi_d - xtild(i))
64         STOP 1  
65      ELSE            IF (200. * fb < - fa) THEN
66         dzoom = dzooma * pi/180.               fhyp(i) = - 1.
67      ENDIF            ELSE IF (200. * fb < fa) THEN
68                 fhyp(i) = 1.
     print *, ' xzoom(rad), grossism, tau, dzoom (rad):'  
     print *, xzoom, grossism, tau, dzoom  
   
     DO i = 0, nmax2  
        xtild(i) = - pi + FLOAT(i) * depi /nmax2  
     ENDDO  
   
     DO i = nmax, nmax2  
        fa = tau* (dzoom/2. - xtild(i))  
        fb = xtild(i) * (pi - xtild(i))  
   
        IF (200.* fb .LT. - fa) THEN  
           fhyp (i) = - 1.  
        ELSEIF (200. * fb .LT. fa) THEN  
           fhyp (i) = 1.  
        ELSE  
           IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN  
              IF (200.*fb + fa.LT.1.e-10) THEN  
                 fhyp (i) = - 1.  
              ELSEIF (200.*fb - fa.LT.1.e-10) THEN  
                 fhyp (i) = 1.  
              ENDIF  
69            ELSE            ELSE
70               fhyp (i) = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71            ENDIF                  IF (200. * fb + fa < 1e-10) THEN
72         ENDIF                     fhyp(i) = - 1.
73                    ELSE IF (200. * fb - fa < 1e-10) THEN
74         IF (xtild(i) == 0.) fhyp(i) = 1.                     fhyp(i) = 1.
75         IF (xtild(i) == pi) fhyp(i) = -1.                  END IF
76      ENDDO               ELSE
77                    fhyp(i) = TANH(fa / fb)
78      ! Calcul de beta               END IF
79              END IF
80      ffdx = 0.  
81              IF (xtild(i) == 0.) fhyp(i) = 1.
82      DO i = nmax + 1, nmax2            IF (xtild(i) == pi_d) fhyp(i) = -1.
83         xmoy = 0.5 * (xtild(i-1) + xtild(i))         END DO
84         fa = tau* (dzoom/2. - xmoy)  
85         fb = xmoy * (pi - xmoy)         ! Calcul de beta
86    
87         IF (200.* fb .LT. - fa) THEN         ffdx = 0.
88            fxm = - 1.  
89         ELSEIF (200. * fb .LT. fa) THEN         DO i = nmax + 1, 2 * nmax
90            fxm = 1.            xmoy = 0.5 * (xtild(i-1) + xtild(i))
91         ELSE            fa = taux * (dzoom / 2. - xmoy)
92            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN            fb = xmoy * (pi_d - xmoy)
93               IF (200.*fb + fa.LT.1.e-10) THEN  
94                  fxm = - 1.            IF (200. * fb < - fa) THEN
95               ELSEIF (200.*fb - fa.LT.1.e-10) THEN               fxm = - 1.
96                  fxm = 1.            ELSE IF (200. * fb < fa) THEN
97               ENDIF               fxm = 1.
98            ELSE            ELSE
99               fxm = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100            ENDIF                  IF (200. * fb + fa < 1e-10) THEN
101         ENDIF                     fxm = - 1.
102                    ELSE IF (200. * fb - fa < 1e-10) THEN
103         IF (xmoy == 0.) fxm = 1.                     fxm = 1.
104         IF (xmoy == pi) fxm = -1.                  END IF
105                 ELSE
106         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))                  fxm = TANH(fa / fb)
107      ENDDO               END IF
108              END IF
109      beta = (grossism * ffdx - pi) / (ffdx - pi)  
110              IF (xmoy == 0.) fxm = 1.
111      IF (2.*beta - grossism <= 0.) THEN            IF (xmoy == pi_d) fxm = -1.
112         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'  
113         print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'            ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114         STOP 1         END DO
115      ENDIF  
116           print *, "ffdx = ", ffdx
117      ! calcul de Xprimt         beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118           print *, "beta = ", beta
119      DO i = nmax, nmax2  
120         Xprimt(i) = beta + (grossism - beta) * fhyp(i)         IF (2. * beta - grossismx <= 0.) THEN
121      ENDDO            print *, 'Bad choice of grossismx, taux, dzoomx.'
122              print *, 'Decrease dzoomx or grossismx.'
123      DO i = nmax + 1, nmax2            STOP 1
124         Xprimt(nmax2 - i) = Xprimt(i)         END IF
125      ENDDO  
126           ! calcul de Xprimt
127      ! Calcul de Xf         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128           xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129      Xf(0) = - pi  
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         xmoy = 0.5 * (xtild(i-1) + xtild(i))            IF (xmoy == pi_d) fxm = -1.
147         fa = tau* (dzoom/2. - xmoy)            xxpr(i) = beta + (grossismx - beta) * fxm
148         fb = xmoy * (pi - xmoy)         END DO
149    
150           xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151    
152           Xf(0) = - pi_d
153    
154           DO i=1, 2 * nmax - 1
155              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156           END DO
157    
158           Xf(2 * nmax) = pi_d
159    
160           call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161                xuv = - 0.25d0)
162           call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163                xuv = 0d0)
164           call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165                xuv = 0.5d0)
166           call invert_zoom_x(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         IF (200.* fb .LT. - fa) THEN            do while (rlonm025(is2) < - pi .and. is2 < iim)
178            fxm = - 1.               is2 = is2 + 1
        ELSEIF (200. * fb .LT. fa) THEN  
           fxm = 1.  
        ELSE  
           fxm = TANH (fa/fb)  
        ENDIF  
   
        IF (xmoy == 0.) fxm = 1.  
        IF (xmoy == pi) fxm = -1.  
        xxpr(i) = beta + (grossism - beta) * fxm  
     ENDDO  
   
     DO i = nmax + 1, nmax2  
        xxpr(nmax2-i + 1) = xxpr(i)  
     ENDDO  
   
     DO i=1, nmax2  
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     ENDDO  
   
     ! xuv = 0. si calcul aux pts scalaires  
     ! xuv = 0.5 si calcul aux pts U  
   
     print *  
   
     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  
        ENDIF  
   
        xo1 = 0.  
   
        ii1=1  
        ii2=iim  
        IF (ik == 1.and.grossism == 1.) THEN  
           ii1 = 2  
           ii2 = iim + 1  
        ENDIF  
   
        DO i = ii1, ii2  
           xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)  
           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  
           ENDIF  
   
           ! 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) <= epsilon .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) > epsilon) 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              do while (rlonm025(is2) > pi .and. is2 > 1)
189                 is2 = is2 - 1
190              end do
191    
192            xxprim(i) = depi/ (FLOAT(iim) * Xprimin)            if (rlonm025(is2) > pi) then
193            xvrai(i) = xi + xzoom               print *, 'Rlonm025 plus grand que pi !'
        end DO  
   
        IF (ik == 1.and.grossism == 1.) THEN  
           xvrai(1) = xvrai(iip1)-depi  
           xxprim(1) = xxprim(iip1)  
        ENDIF  
        DO i = 1, iim  
           xlon(i) = xvrai(i)  
           xprimm(i) = xxprim(i)  
        ENDDO  
        DO i = 1, iim -1  
           IF (xvrai(i + 1).LT. xvrai(i)) THEN  
              print *, 'Problème avec rlonu(', i + 1, &  
                   ') plus petit que rlonu(', i, ')'  
194               STOP 1               STOP 1
195            ENDIF            end if
196         ENDDO         END IF
197        END IF
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
198    
199         champmin = 1.e12      call principal_cshift(is2, rlonm025, xprimm025)
200         champmax = -1.e12      call principal_cshift(is2, rlonv, xprimv)
201         DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
202            champmin = MIN(champmin, xvrai(i))      call principal_cshift(is2, rlonp025, xprimp025)
203            champmax = MAX(champmax, xvrai(i))  
204         ENDDO      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
206         IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) THEN      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
207            print *, 'Reorganisation des longitudes pour avoir entre - pi', &  
208                 ' et pi '      ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
209        DO i = 1, iim + 1
210            IF (xzoom <= 0.) THEN         IF (rlonp025(i) < rlonv(i)) THEN
211               IF (ik == 1) THEN            print *, 'rlonp025(', i, ') = ', rlonp025(i)
212                  i = 1            print *, "< rlonv(", i, ") = ", rlonv(i)
213              STOP 1
214                  do while (xvrai(i) < - pi .and. i < iim)         END IF
215                     i = i + 1  
216                  end do         IF (rlonv(i) < rlonm025(i)) THEN
217              print *, 'rlonv(', i, ') = ', rlonv(i)
218                  if (xvrai(i) < - pi) then            print *, "< rlonm025(", i, ") = ", rlonm025(i)
219                     print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '            STOP 1
220                     STOP 1         END IF
221                  end if  
222           IF (rlonp025(i) > rlonu(i)) THEN
223                  is2 = i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
224               ENDIF            print *, "> rlonu(", i, ") = ", rlonu(i)
225              STOP 1
226               IF (is2.NE. 1) THEN         END IF
227                  DO ii = is2, iim      END DO
                    xlon (ii-is2 + 1) = xvrai(ii)  
                    xprimm(ii-is2 + 1) = xxprim(ii)  
                 ENDDO  
                 DO ii = 1, is2 -1  
                    xlon (ii + iim-is2 + 1) = xvrai(ii) + depi  
                    xprimm(ii + iim-is2 + 1) = xxprim(ii)  
                 ENDDO  
              ENDIF  
           ELSE  
              IF (ik == 1) THEN  
                 i = iim  
   
                 do while (xvrai(i) > pi .and. i > 1)  
                    i = i - 1  
                 end do  
   
                 if (xvrai(i) > pi) then  
                    print *, ' PBS. 2 ! Xvrai plus grand que pi ! '  
                    STOP 1  
                 end if  
   
                 is2 = i  
              ENDIF  
              idif = iim -is2  
              DO ii = 1, is2  
                 xlon (ii + idif) = xvrai(ii)  
                 xprimm(ii + idif) = xxprim(ii)  
              ENDDO  
              DO ii = 1, idif  
                 xlon (ii) = xvrai (ii + is2) - depi  
                 xprimm(ii) = xxprim(ii + is2)  
              ENDDO  
           ENDIF  
        ENDIF  
   
        ! Fin de la reorganisation  
   
        xlon (iip1) = xlon(1) + depi  
        xprimm(iip1) = xprimm (1)  
   
        DO i = 1, iim + 1  
           xvrai(i) = xlon(i)*180./pi  
        ENDDO  
   
        IF (ik == 1) THEN  
           DO i = 1, iim + 1  
              rlonm025(i) = xlon(i)  
              xprimm025(i) = xprimm(i)  
           ENDDO  
        ELSE IF (ik == 2) THEN  
           DO i = 1, iim + 1  
              rlonv(i) = xlon(i)  
              xprimv(i) = xprimm(i)  
           ENDDO  
        ELSE IF (ik == 3) THEN  
           DO i = 1, iim + 1  
              rlonu(i) = xlon(i)  
              xprimu(i) = xprimm(i)  
           ENDDO  
        ELSE IF (ik == 4) THEN  
           DO i = 1, iim + 1  
              rlonp025(i) = xlon(i)  
              xprimp025(i) = xprimm(i)  
           ENDDO  
        ENDIF  
     end DO  
   
     print *  
   
     DO i = 1, iim  
        xlon(i) = rlonv(i + 1) - rlonv(i)  
     ENDDO  
     champmin = 1.e12  
     champmax = -1.e12  
     DO i = 1, iim  
        champmin = MIN(champmin, xlon(i))  
        champmax = MAX(champmax, xlon(i))  
     ENDDO  
     champmin = champmin * 180./pi  
     champmax = champmax * 180./pi  
228    
229    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
230    
Legend:
Removed from v.91  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21