/[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 91 by guez, Wed Mar 26 17:18:58 2014 UTC revision 122 by guez, Tue Feb 3 19:30:48 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  
   
     ! arguments de sortie  
19    
20      REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv      USE dimens_m, ONLY: iim
21      real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22        use nr_util, only: pi_d, twopi_d, arth
23        use serre, only: clon, grossismx, dzoomx, taux
24    
25      DOUBLE PRECISION, intent(out):: champmin, champmax      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
26        real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
27    
28      ! Local:      ! Local:
29    
30      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      real rlonm025(iim + 1), rlonp025(iim + 1)
   
     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.  
   
31      REAL dzoom      REAL dzoom
32      DOUBLE PRECISION xlon(iip1), xprimm(iip1), xuv      DOUBLE PRECISION xlon(iim)
33      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION xtild(0:2 * nmax)
34      DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35      DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)      DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36      DOUBLE PRECISION xvrai(iip1), xxprim(iip1)      DOUBLE PRECISION xzoom, fa, fb
37      DOUBLE PRECISION pi, depi, epsilon, xzoom, fa, fb      INTEGER i
38      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2      DOUBLE PRECISION xmoy, fxm
     INTEGER i, it, ik, iter, ii, idif, ii1, ii2  
     DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin  
39      DOUBLE PRECISION decalx      DOUBLE PRECISION decalx
     INTEGER is2  
     SAVE is2  
   
     DOUBLE PRECISION heavyside  
40    
41      !----------------------------------------------------------------------      !----------------------------------------------------------------------
42    
43      pi = 2. * ASIN(1.)      print *, "Call sequence information: fxhyp"
44      depi = 2. * pi  
45      epsilon = 1.e-3      xzoom = clon * pi_d / 180d0
46      xzoom = xzoomdeg * pi/180.  
47        IF (grossismx == 1.) THEN
48      decalx = .75         decalx = 1d0
49      IF (grossism == 1. .AND. scal180) THEN      else
50         decalx = 1.         decalx = 0.75d0
51      ENDIF      END IF
52    
53      print *, 'FXHYP scal180, decalx', scal180, decalx      dzoom = dzoomx * twopi_d
54        xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
55      IF (dzooma.LT.1.) THEN  
56         dzoom = dzooma * depi      ! Compute fhyp:
57      ELSEIF (dzooma.LT. 25.) THEN      DO i = nmax, 2 * nmax
58         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &         fa = taux * (dzoom / 2. - xtild(i))
59              // "L'augmenter et relancer."         fb = xtild(i) * (pi_d - xtild(i))
60         STOP 1  
61      ELSE         IF (200. * fb < - fa) THEN
62         dzoom = dzooma * pi/180.            fhyp(i) = - 1.
63      ENDIF         ELSE IF (200. * fb < fa) THEN
64              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.  
65         ELSE         ELSE
66            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
67               IF (200.*fb + fa.LT.1.e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
68                  fhyp (i) = - 1.                  fhyp(i) = - 1.
69               ELSEIF (200.*fb - fa.LT.1.e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
70                  fhyp (i) = 1.                  fhyp(i) = 1.
71               ENDIF               END IF
72            ELSE            ELSE
73               fhyp (i) = TANH (fa/fb)               fhyp(i) = TANH(fa / fb)
74            ENDIF            END IF
75         ENDIF         END IF
76    
77         IF (xtild(i) == 0.) fhyp(i) = 1.         IF (xtild(i) == 0.) fhyp(i) = 1.
78         IF (xtild(i) == pi) fhyp(i) = -1.         IF (xtild(i) == pi_d) fhyp(i) = -1.
79      ENDDO      END DO
80    
81      ! Calcul de beta      ! Calcul de beta
82    
83      ffdx = 0.      ffdx = 0.
84    
85      DO i = nmax + 1, nmax2      DO i = nmax + 1, 2 * nmax
86         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
87         fa = tau* (dzoom/2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
88         fb = xmoy * (pi - xmoy)         fb = xmoy * (pi_d - xmoy)
89    
90         IF (200.* fb .LT. - fa) THEN         IF (200. * fb < - fa) THEN
91            fxm = - 1.            fxm = - 1.
92         ELSEIF (200. * fb .LT. fa) THEN         ELSE IF (200. * fb < fa) THEN
93            fxm = 1.            fxm = 1.
94         ELSE         ELSE
95            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN            IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
96               IF (200.*fb + fa.LT.1.e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
97                  fxm = - 1.                  fxm = - 1.
98               ELSEIF (200.*fb - fa.LT.1.e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
99                  fxm = 1.                  fxm = 1.
100               ENDIF               END IF
101            ELSE            ELSE
102               fxm = TANH (fa/fb)               fxm = TANH(fa / fb)
103            ENDIF            END IF
104         ENDIF         END IF
105    
106         IF (xmoy == 0.) fxm = 1.         IF (xmoy == 0.) fxm = 1.
107         IF (xmoy == pi) fxm = -1.         IF (xmoy == pi_d) fxm = -1.
108    
109         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
110      ENDDO      END DO
   
     beta = (grossism * ffdx - pi) / (ffdx - pi)  
111    
112      IF (2.*beta - grossism <= 0.) THEN      print *, "ffdx = ", ffdx
113         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
114         print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'      print *, "beta = ", beta
115    
116        IF (2. * beta - grossismx <= 0.) THEN
117           print *, 'Bad choice of grossismx, taux, dzoomx.'
118           print *, 'Decrease dzoomx or grossismx.'
119         STOP 1         STOP 1
120      ENDIF      END IF
121    
122      ! calcul de Xprimt      ! calcul de Xprimt
123        Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
124        xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
125    
126      DO i = nmax, nmax2      ! Calcul de Xf
        Xprimt(i) = beta + (grossism - beta) * fhyp(i)  
     ENDDO  
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     ENDDO  
   
     ! Calcul de Xf  
127    
128      Xf(0) = - pi      DO i = nmax + 1, 2 * nmax
   
     DO i = nmax + 1, nmax2  
129         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
130         fa = tau* (dzoom/2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
131         fb = xmoy * (pi - xmoy)         fb = xmoy * (pi_d - xmoy)
132    
133         IF (200.* fb .LT. - fa) THEN         IF (200. * fb < - fa) THEN
134            fxm = - 1.            fxm = - 1.
135         ELSEIF (200. * fb .LT. fa) THEN         ELSE IF (200. * fb < fa) THEN
136            fxm = 1.            fxm = 1.
137         ELSE         ELSE
138            fxm = TANH (fa/fb)            fxm = TANH(fa / fb)
139         ENDIF         END IF
140    
141         IF (xmoy == 0.) fxm = 1.         IF (xmoy == 0.) fxm = 1.
142         IF (xmoy == pi) fxm = -1.         IF (xmoy == pi_d) fxm = -1.
143         xxpr(i) = beta + (grossism - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
144      ENDDO      END DO
   
     DO i = nmax + 1, nmax2  
        xxpr(nmax2-i + 1) = xxpr(i)  
     ENDDO  
145    
146      DO i=1, nmax2      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
        Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))  
     ENDDO  
147    
148      ! xuv = 0. si calcul aux pts scalaires      Xf(0) = - pi_d
     ! xuv = 0.5 si calcul aux pts U  
149    
150      print *      DO i=1, 2 * nmax - 1
151           Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
152        END DO
153    
154        Xf(2 * nmax) = pi_d
155    
156      DO ik = 1, 4      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
157         IF (ik == 1) THEN           xprimm025, xuv = - 0.25d0)
158            xuv = -0.25      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
159         ELSE IF (ik == 2) THEN           xuv = 0d0)
160            xuv = 0.      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
161         ELSE IF (ik == 3) THEN           xuv = 0.5d0)
162            xuv = 0.50      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
163         ELSE IF (ik == 4) THEN           xprimp025, xuv = 0.25d0)
           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  
           end do  
   
           ! Calcul de Xf(xi)  
   
           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  
              STOP 1  
           end if  
   
   
           xxprim(i) = depi/ (FLOAT(iim) * Xprimin)  
           xvrai(i) = xi + xzoom  
        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, ')'  
              STOP 1  
           ENDIF  
        ENDDO  
   
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
   
        champmin = 1.e12  
        champmax = -1.e12  
        DO i = 1, iim  
           champmin = MIN(champmin, xvrai(i))  
           champmax = MAX(champmax, xvrai(i))  
        ENDDO  
   
        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  
164    
165      print *      print *
166    
167      DO i = 1, iim      forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
168         xlon(i) = rlonv(i + 1) - rlonv(i)      print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "degrees"
169      ENDDO      print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "degrees"
170      champmin = 1.e12  
171      champmax = -1.e12      DO i = 1, iim + 1
172      DO i = 1, iim         IF (rlonp025(i) < rlonv(i)) THEN
173         champmin = MIN(champmin, xlon(i))            print *, 'rlonp025(', i, ') = ', rlonp025(i)
174         champmax = MAX(champmax, xlon(i))            print *, "< rlonv(", i, ") = ", rlonv(i)
175      ENDDO            STOP 1
176      champmin = champmin * 180./pi         END IF
177      champmax = champmax * 180./pi  
178           IF (rlonv(i) < rlonm025(i)) THEN
179              print *, 'rlonv(', i, ') = ', rlonv(i)
180              print *, "< rlonm025(", i, ") = ", rlonm025(i)
181              STOP 1
182           END IF
183    
184           IF (rlonp025(i) > rlonu(i)) THEN
185              print *, 'rlonp025(', i, ') = ', rlonp025(i)
186              print *, "> rlonu(", i, ") = ", rlonu(i)
187              STOP 1
188           END IF
189        END DO
190    
191    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
192    
Legend:
Removed from v.91  
changed lines
  Added in v.122
  ViewVC Help
Powered by ViewVC 1.1.21