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

Diff of /trunk/Sources/dyn3d/fxhyp.f

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

revision 105 by guez, Thu Sep 4 10:40:24 2014 UTC revision 126 by guez, Fri Feb 6 18:33:15 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 fxhyp_loop_ik_m, only: fxhyp_loop_ik, 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      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 (grossism = 1., tau=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(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  
39    
40      !----------------------------------------------------------------------      !----------------------------------------------------------------------
41    
42      pi = 2. * ASIN(1.)      print *, "Call sequence information: fxhyp"
     depi = 2. * pi  
     epsilon = 1.e-3  
     xzoom = xzoomdeg * pi/180.  
   
     decalx = .75  
     IF (grossism == 1. .AND. scal180) THEN  
        decalx = 1.  
     ENDIF  
   
     print *, 'FXHYP scal180, decalx', scal180, decalx  
   
     IF (dzooma.LT.1.) THEN  
        dzoom = dzooma * depi  
     ELSEIF (dzooma.LT. 25.) THEN  
        print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &  
             // "L'augmenter et relancer."  
        STOP 1  
     ELSE  
        dzoom = dzooma * pi/180.  
     ENDIF  
   
     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  
           ELSE  
              fhyp (i) = TANH (fa/fb)  
           ENDIF  
        ENDIF  
   
        IF (xtild(i) == 0.) fhyp(i) = 1.  
        IF (xtild(i) == pi) fhyp(i) = -1.  
     ENDDO  
   
     ! Calcul de beta  
   
     ffdx = 0.  
   
     DO i = nmax + 1, nmax2  
        xmoy = 0.5 * (xtild(i-1) + xtild(i))  
        fa = tau* (dzoom/2. - xmoy)  
        fb = xmoy * (pi - xmoy)  
   
        IF (200.* fb .LT. - fa) THEN  
           fxm = - 1.  
        ELSEIF (200. * fb .LT. fa) THEN  
           fxm = 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  
                 fxm = - 1.  
              ELSEIF (200.*fb - fa.LT.1.e-10) THEN  
                 fxm = 1.  
              ENDIF  
           ELSE  
              fxm = TANH (fa/fb)  
           ENDIF  
        ENDIF  
   
        IF (xmoy == 0.) fxm = 1.  
        IF (xmoy == pi) fxm = -1.  
   
        ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))  
     ENDDO  
43    
44      beta = (grossism * ffdx - pi) / (ffdx - pi)      xzoom = clon * pi_d / 180d0
45    
46      IF (2.*beta - grossism <= 0.) THEN      test_grossismx: if (grossismx == 1.) then
47         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'         step = twopi / iim
        print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'  
        STOP 1  
     ENDIF  
48    
49      ! calcul de Xprimt         xprimm025(:iim) = step
50           xprimp025(:iim) = step
51      DO i = nmax, nmax2         xprimv(:iim) = step
52         Xprimt(i) = beta + (grossism - beta) * fhyp(i)         xprimu(:iim) = step
53      ENDDO  
54           rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55      DO i = nmax + 1, nmax2         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56         Xprimt(nmax2 - i) = Xprimt(i)         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57      ENDDO         rlonu(:iim) = rlonv(:iim) + 0.5 * step
58        else
59      ! Calcul de Xf         dzoom = dzoomx * twopi_d
60           xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61      Xf(0) = - pi  
62           ! Compute fhyp:
63           DO i = nmax, 2 * nmax
64              fa = taux * (dzoom / 2. - xtild(i))
65              fb = xtild(i) * (pi_d - xtild(i))
66    
67              IF (200. * fb < - fa) THEN
68                 fhyp(i) = - 1.
69              ELSE IF (200. * fb < fa) THEN
70                 fhyp(i) = 1.
71              ELSE
72                 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
82    
83              IF (xtild(i) == 0.) fhyp(i) = 1.
84              IF (xtild(i) == pi_d) fhyp(i) = -1.
85           END DO
86    
87           ! Calcul de beta
88    
89           ffdx = 0.
90    
91           DO i = nmax + 1, 2 * nmax
92              xmoy = 0.5 * (xtild(i-1) + xtild(i))
93              fa = taux * (dzoom / 2. - xmoy)
94              fb = xmoy * (pi_d - xmoy)
95    
96              IF (200. * fb < - fa) THEN
97                 fxm = - 1.
98              ELSE IF (200. * fb < fa) THEN
99                 fxm = 1.
100              ELSE
101                 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
111    
112              IF (xmoy == 0.) fxm = 1.
113              IF (xmoy == pi_d) fxm = -1.
114    
115              ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
116           END DO
117    
118           print *, "ffdx = ", ffdx
119           beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
120           print *, "beta = ", beta
121    
122           IF (2. * beta - grossismx <= 0.) THEN
123              print *, 'Bad choice of grossismx, taux, dzoomx.'
124              print *, 'Decrease dzoomx or grossismx.'
125              STOP 1
126           END IF
127    
128           ! calcul de Xprimt
129           Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130           xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131    
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         xmoy = 0.5 * (xtild(i-1) + xtild(i))            IF (xmoy == pi_d) fxm = -1.
149         fa = tau* (dzoom/2. - xmoy)            xxpr(i) = beta + (grossismx - beta) * fxm
150         fb = xmoy * (pi - xmoy)         END DO
151    
152           xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153    
154           Xf(0) = - pi_d
155    
156           DO i=1, 2 * nmax - 1
157              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
158           END DO
159    
160           Xf(2 * nmax) = pi_d
161    
162           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonm025(:iim), &
163                xprimm025(:iim), xuv = - 0.25d0)
164           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonv(:iim), &
165                xprimv(:iim), xuv = 0d0)
166           call fxhyp_loop_ik(xf, xtild, Xprimt, xzoom, rlonu(:iim), &
167                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         IF (200.* fb .LT. - fa) THEN            do while (rlonm025(is2) < - pi .and. is2 < iim)
180            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  
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  
           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  
185               STOP 1               STOP 1
186            end if            end if
187           ELSE
188              is2 = iim
189    
190              do while (rlonm025(is2) > pi .and. is2 > 1)
191                 is2 = is2 - 1
192              end do
193    
194            xxprim(i) = depi/ (FLOAT(iim) * Xprimin)            if (rlonm025(is2) > pi) then
195            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, ')'  
196               STOP 1               STOP 1
197            ENDIF            end if
198         ENDDO         END IF
199        END IF
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
200    
201         champmin = 1.e12      call principal_cshift(is2, rlonm025, xprimm025)
202         champmax = -1.e12      call principal_cshift(is2, rlonv, xprimv)
203         DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
204            champmin = MIN(champmin, xvrai(i))      call principal_cshift(is2, rlonp025, xprimp025)
205            champmax = MAX(champmax, xvrai(i))  
206         ENDDO      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208         IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) THEN      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
209            print *, 'Reorganisation des longitudes pour avoir entre - pi', &  
210                 ' et pi '      DO i = 1, iim + 1
211           IF (rlonp025(i) < rlonv(i)) THEN
212            IF (xzoom <= 0.) THEN            print *, 'rlonp025(', i, ') = ', rlonp025(i)
213               IF (ik == 1) THEN            print *, "< rlonv(", i, ") = ", rlonv(i)
214                  i = 1            STOP 1
215           END IF
216                  do while (xvrai(i) < - pi .and. i < iim)  
217                     i = i + 1         IF (rlonv(i) < rlonm025(i)) THEN
218                  end do            print *, 'rlonv(', i, ') = ', rlonv(i)
219              print *, "< rlonm025(", i, ") = ", rlonm025(i)
220                  if (xvrai(i) < - pi) then            STOP 1
221                     print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '         END IF
222                     STOP 1  
223                  end if         IF (rlonp025(i) > rlonu(i)) THEN
224              print *, 'rlonp025(', i, ') = ', rlonp025(i)
225                  is2 = i            print *, "> rlonu(", i, ") = ", rlonu(i)
226               ENDIF            STOP 1
227           END IF
228               IF (is2.NE. 1) THEN      END DO
                 DO ii = is2, iim  
                    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  
229    
230    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
231    
Legend:
Removed from v.105  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21