/[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 167 by guez, Mon Aug 24 16:30:33 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 x_f(\tilde x) à dérivée tangente hyperbolique.
14    
15      ! On doit avoir grossism \times dzoom < pi (radians), en longitude.      ! Il vaut mieux avoir : grossismx \times delta < 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) avec clon = 0 est à - 180 degrés.
19    
20      REAL, intent(in):: xzoomdeg      USE dimens_m, ONLY: iim
21        use dynetat0_m, only: clon, grossismx, dzoomx, taux
22        use invert_zoom_x_m, only: invert_zoom_x, nmax
23        use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24        use principal_cshift_m, only: principal_cshift
25        use tanh_cautious_m, only: tanh_cautious
26    
27      REAL, intent(in):: grossism      REAL, intent(out):: xprimm025(:) ! (iim + 1)
     ! grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)  
28    
29      REAL, intent(in):: dzooma ! distance totale de la zone du zoom      REAL, intent(out):: rlonv(:) ! (iim + 1)
30        ! longitudes of points of the "scalar" and "v" grid, in rad
31    
32      REAL, intent(in):: tau      REAL, intent(out):: xprimv(:) ! (iim + 1)
33      ! raideur de la transition de l'intérieur à l'extérieur du zoom      ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonv)
34    
35      ! arguments de sortie      real, intent(out):: rlonu(:) ! (iim + 1)
36        ! longitudes of points of the "u" grid, in rad
37    
38      REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv      real, intent(out):: xprimu(:) ! (iim + 1)
39      real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025      ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlonu)
40    
41      DOUBLE PRECISION, intent(out):: champmin, champmax      real, intent(out):: xprimp025(:) ! (iim + 1)
42    
43      ! Local:      ! Local:
44        real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
45      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      REAL delta, h
46        DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
47      LOGICAL, PARAMETER:: scal180 = .TRUE.      DOUBLE PRECISION beta
48      ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour      INTEGER i, is2
49      ! une grille reguliere (grossism = 1., tau=0., clon=0.) a      DOUBLE PRECISION xmoy(nmax), fxm(nmax)
     ! -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  
   
     DOUBLE PRECISION heavyside  
50    
51      !----------------------------------------------------------------------      !----------------------------------------------------------------------
52    
53      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  
   
     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  
     ENDIF  
   
     ! calcul de Xprimt  
   
     DO i = nmax, nmax2  
        Xprimt(i) = beta + (grossism - beta) * fhyp(i)  
     ENDDO  
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     ENDDO  
   
     ! Calcul de Xf  
   
     Xf(0) = - pi  
   
     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  
           fxm = TANH (fa/fb)  
        ENDIF  
54    
55         IF (xmoy == 0.) fxm = 1.      if (grossismx == 1.) then
56         IF (xmoy == pi) fxm = -1.         h = twopi / iim
        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  
           end do  
57    
58            ! Calcul de Xf(xi)         xprimm025(:iim) = h
59           xprimp025(:iim) = h
60           xprimv(:iim) = h
61           xprimu(:iim) = h
62    
63           rlonv(:iim) = arth(- pi + clon, h, iim)
64           rlonm025(:iim) = rlonv(:iim) - 0.25 * h
65           rlonp025(:iim) = rlonv(:iim) + 0.25 * h
66           rlonu(:iim) = rlonv(:iim) + 0.5 * h
67        else
68           delta = dzoomx * twopi_d
69           xtild = arth(0d0, pi_d / nmax, nmax + 1)
70           forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
71    
72           ! Compute fhyp:
73           fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
74                - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
75                * (pi_d - xtild(1:nmax - 1)))
76           fhyp(0) = 1d0
77           fhyp(nmax) = -1d0
78    
79           fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))
80    
81           ! Compute \int_0 ^{\tilde x} F:
82    
83           ffdx(0) = 0d0
84    
85           DO i = 1, nmax
86              ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
87           END DO
88    
89           print *, "ffdx(nmax) = ", ffdx(nmax)
90           beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
91           print *, "beta = ", beta
92    
93           IF (2d0 * beta - grossismx <= 0d0) THEN
94              print *, 'Bad choice of grossismx, taux, dzoomx.'
95              print *, 'Decrease dzoomx or grossismx.'
96              STOP 1
97           END IF
98    
99           G = beta + (grossismx - beta) * fhyp
100    
101           Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
102                * ffdx(:nmax - 1)
103           Xf(nmax) = pi_d
104    
105           call invert_zoom_x(beta, xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
106                xuv = - 0.25d0)
107           call invert_zoom_x(beta, xf, xtild, G, rlonv(:iim), xprimv(:iim), &
108                xuv = 0d0)
109           call invert_zoom_x(beta, xf, xtild, G, rlonu(:iim), xprimu(:iim), &
110                xuv = 0.5d0)
111           call invert_zoom_x(beta, xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
112                xuv = 0.25d0)
113        end if
114    
115        is2 = 0
116    
117        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
118             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
119           IF (clon <= 0.) THEN
120              is2 = 1
121    
122            xi = xtild(it)            do while (rlonm025(is2) < - pi .and. is2 < iim)
123                 is2 = is2 + 1
124              end do
125    
126            IF (it == nmax2) THEN            if (rlonm025(is2) < - pi) then
127               it = nmax2 -1               print *, 'Rlonm025 plus petit que - pi !'
              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  
128               STOP 1               STOP 1
129            end if            end if
130           ELSE
131              is2 = iim
132    
133              do while (rlonm025(is2) > pi .and. is2 > 1)
134                 is2 = is2 - 1
135              end do
136    
137            xxprim(i) = depi/ (FLOAT(iim) * Xprimin)            if (rlonm025(is2) > pi) then
138            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, ')'  
139               STOP 1               STOP 1
140            ENDIF            end if
141         ENDDO         END IF
142        END IF
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
143    
144         champmin = 1.e12      call principal_cshift(is2, rlonm025, xprimm025)
145         champmax = -1.e12      call principal_cshift(is2, rlonv, xprimv)
146         DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
147            champmin = MIN(champmin, xvrai(i))      call principal_cshift(is2, rlonp025, xprimp025)
148            champmax = MAX(champmax, xvrai(i))  
149         ENDDO      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
150        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
151         IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) THEN      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
152            print *, 'Reorganisation des longitudes pour avoir entre - pi', &  
153                 ' et pi '      ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
154        DO i = 1, iim + 1
155            IF (xzoom <= 0.) THEN         IF (rlonp025(i) < rlonv(i)) THEN
156               IF (ik == 1) THEN            print *, 'rlonp025(', i, ') = ', rlonp025(i)
157                  i = 1            print *, "< rlonv(", i, ") = ", rlonv(i)
158              STOP 1
159                  do while (xvrai(i) < - pi .and. i < iim)         END IF
160                     i = i + 1  
161                  end do         IF (rlonv(i) < rlonm025(i)) THEN
162              print *, 'rlonv(', i, ') = ', rlonv(i)
163                  if (xvrai(i) < - pi) then            print *, "< rlonm025(", i, ") = ", rlonm025(i)
164                     print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '            STOP 1
165                     STOP 1         END IF
166                  end if  
167           IF (rlonp025(i) > rlonu(i)) THEN
168                  is2 = i            print *, 'rlonp025(', i, ') = ', rlonp025(i)
169               ENDIF            print *, "> rlonu(", i, ") = ", rlonu(i)
170              STOP 1
171               IF (is2.NE. 1) THEN         END IF
172                  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  
173    
174    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
175    
Legend:
Removed from v.91  
changed lines
  Added in v.167
  ViewVC Help
Powered by ViewVC 1.1.21