/[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 112 by guez, Thu Sep 18 13:36:51 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  
   
     REAL, intent(in):: tau  
     ! raideur de la transition de l'intérieur à l'extérieur du zoom  
19    
20      ! arguments de sortie      USE dimens_m, ONLY: iim
21        use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22      REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv      use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23      real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025      use principal_cshift_m, only: principal_cshift
24        use serre, only: clon, grossismx, dzoomx, taux
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.  
     END IF  
43    
44      print *, ' xzoom(rad), grossism, tau, dzoom (rad):'      xzoom = clon * pi_d / 180d0
     print *, xzoom, grossism, tau, dzoom  
45    
46      DO i = 0, nmax2      test_grossismx: if (grossismx == 1.) then
47         xtild(i) = - pi + REAL(i) * depi /nmax2         step = twopi / iim
48      ENDDO  
49           xprimm025(:iim) = step
50      DO i = nmax, nmax2         xprimp025(:iim) = step
51         fa = tau* (dzoom/2. - xtild(i))         xprimv(:iim) = step
52         fb = xtild(i) * (pi - xtild(i))         xprimu(:iim) = step
53    
54         IF (200.* fb .LT. - fa) THEN         rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55            fhyp (i) = - 1.         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56         ELSEIF (200. * fb .LT. fa) THEN         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57            fhyp (i) = 1.         rlonu(:iim) = rlonv(:iim) + 0.5 * step
58         ELSE      else
59            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN         dzoom = dzoomx * twopi_d
60               IF (200.*fb + fa.LT.1.e-10) THEN         xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61                  fhyp (i) = - 1.  
62               ELSEIF (200.*fb - fa.LT.1.e-10) THEN         ! Compute fhyp:
63                  fhyp (i) = 1.         DO i = nmax, 2 * nmax
64               ENDIF            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            ELSE
101               fhyp (i) = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
102            ENDIF                  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         END IF
127    
128         IF (xtild(i) == 0.) fhyp(i) = 1.         ! calcul de Xprimt
129         IF (xtild(i) == pi) fhyp(i) = -1.         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130      END DO         xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131    
132      ! Calcul de beta         ! Calcul de Xf
133    
134      ffdx = 0.         DO i = nmax + 1, 2 * nmax
135              xmoy = 0.5 * (xtild(i-1) + xtild(i))
136      DO i = nmax + 1, nmax2            fa = taux * (dzoom / 2. - xmoy)
137         xmoy = 0.5 * (xtild(i-1) + xtild(i))            fb = xmoy * (pi_d - xmoy)
138         fa = tau* (dzoom/2. - xmoy)  
139         fb = xmoy * (pi - xmoy)            IF (200. * fb < - fa) THEN
140                 fxm = - 1.
141         IF (200.* fb .LT. - fa) THEN            ELSE IF (200. * fb < fa) THEN
142            fxm = - 1.               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  
143            ELSE            ELSE
144               fxm = TANH (fa/fb)               fxm = TANH(fa / fb)
145            ENDIF            END IF
        ENDIF  
   
        IF (xmoy == 0.) fxm = 1.  
        IF (xmoy == pi) fxm = -1.  
   
        ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))  
     ENDDO  
   
     beta = (grossism * ffdx - pi) / (ffdx - pi)  
   
     IF (2.*beta - grossism <= 0.) THEN  
        print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'  
        print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'  
        STOP 1  
     END IF  
146    
147      ! calcul de Xprimt            IF (xmoy == 0.) fxm = 1.
148              IF (xmoy == pi_d) fxm = -1.
149              xxpr(i) = beta + (grossismx - beta) * fxm
150           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      DO i = nmax, nmax2            do while (rlonm025(is2) < - pi .and. is2 < iim)
180         Xprimt(i) = beta + (grossism - beta) * fhyp(i)               is2 = is2 + 1
181      END DO            end do
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
   
     ! Calcul de Xf  
   
     Xf(0) = - pi  
182    
183      DO i = nmax + 1, nmax2            if (rlonm025(is2) < - pi) then
184         xmoy = 0.5 * (xtild(i-1) + xtild(i))               print *, 'Rlonm025 plus petit que - pi !'
185         fa = tau* (dzoom/2. - xmoy)               STOP 1
186         fb = xmoy * (pi - xmoy)            end if
   
        IF (200.* fb .LT. - fa) THEN  
           fxm = - 1.  
        ELSEIF (200. * fb .LT. fa) THEN  
           fxm = 1.  
187         ELSE         ELSE
188            fxm = TANH (fa/fb)            is2 = iim
        ENDIF  
189    
190         IF (xmoy == 0.) fxm = 1.            do while (rlonm025(is2) > pi .and. is2 > 1)
191         IF (xmoy == pi) fxm = -1.               is2 = is2 - 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 + (REAL(i) + xuv - decalx) * depi / REAL(iim)  
           Xfi = xlon2  
   
           it = nmax2  
           do while (xfi < xf(it) .and. it >= 1)  
              it = it - 1  
192            end do            end do
193    
194            ! Calcul de Xf(xi)            if (rlonm025(is2) > pi) then
195                 print *, 'Rlonm025 plus grand 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  
196               STOP 1               STOP 1
197            end if            end if
198           END IF
199        END IF
200    
201        call principal_cshift(is2, rlonm025, xprimm025)
202        call principal_cshift(is2, rlonv, xprimv)
203        call principal_cshift(is2, rlonu, xprimu)
204        call principal_cshift(is2, rlonp025, xprimp025)
205    
206        forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208        print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
209    
210        DO i = 1, iim + 1
211           IF (rlonp025(i) < rlonv(i)) THEN
212              print *, 'rlonp025(', i, ') = ', rlonp025(i)
213              print *, "< rlonv(", i, ") = ", rlonv(i)
214              STOP 1
215           END IF
216    
217            xxprim(i) = depi/ (REAL(iim) * Xprimin)         IF (rlonv(i) < rlonm025(i)) THEN
218            xvrai(i) = xi + xzoom            print *, 'rlonv(', i, ') = ', rlonv(i)
219         end DO            print *, "< rlonm025(", i, ") = ", rlonm025(i)
220              STOP 1
221         IF (ik == 1.and.grossism == 1.) THEN         END IF
           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, ')'  
              STOP 1  
           ENDIF  
        ENDDO  
   
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
222    
223         champmin = 1.e12         IF (rlonp025(i) > rlonu(i)) THEN
224         champmax = -1.e12            print *, 'rlonp025(', i, ') = ', rlonp025(i)
225         DO i = 1, iim            print *, "> rlonu(", i, ") = ", rlonu(i)
226            champmin = MIN(champmin, xvrai(i))            STOP 1
227            champmax = MAX(champmax, xvrai(i))         END IF
228         ENDDO      END DO
   
        IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) 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 .and. i < iim)  
                    i = i + 1  
                 end do  
   
                 if (xvrai(i) < - pi) then  
                    print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '  
                    STOP 1  
                 end if  
   
                 is2 = i  
              ENDIF  
   
              IF (is2.NE. 1) THEN  
                 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.112  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21