/[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 82 by guez, Wed Mar 5 14:57:53 2014 UTC revision 121 by guez, Wed Jan 28 16:10:02 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, v 1.2 2005/06/03 09:11:32 fairhead      ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10        ! Author: P. Le Van, from formulas by R. Sadourny
     ! Auteur : P. Le Van  
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.
   
     ! On doit avoir grossism \times dzoom < pi (radians), en longitude.  
   
     USE dimens_m, ONLY: iim  
     USE paramet_m, ONLY: iip1  
   
     INTEGER nmax, nmax2  
     PARAMETER (nmax = 30000, nmax2 = 2*nmax)  
   
     LOGICAL scal180  
     PARAMETER (scal180 = .TRUE.)  
14    
15      ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour      ! Il vaut mieux avoir : grossismx \times dzoom < pi
     ! une grille reguliere (grossism = 1., tau=0., clon=0.) a -180. degres.  
     ! sinon scal180 = .FALSE.  
16    
17      ! ...... arguments d'entree .......      ! Le premier point scalaire pour une grille regulière (grossismx =
18        ! 1., taux=0., clon=0.) est à - 180 degrés.
19    
20      REAL xzoomdeg, dzooma, tau, grossism      USE dimens_m, ONLY: iim
21      ! grossism etant le grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22      ! dzooma etant la distance totale de la zone du zoom      use nr_util, only: pi_d, twopi_d, arth
23      ! tau la raideur de la transition de l'interieur a l'exterieur du zoom      use serre, only: clon, grossismx, dzoomx, taux
   
     ! ...... arguments de sortie ......  
24    
25      REAL rlonm025(iip1), xprimm025(iip1), rlonv(iip1), xprimv(iip1), &      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
26           rlonu(iip1), xprimu(iip1), rlonp025(iip1), xprimp025(iip1)      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
27    
28      ! .... variables locales ....      ! Local:
29    
30        real rlonm025(iim + 1), rlonp025(iim + 1)
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
39      INTEGER i, it, ik, iter, ii, idif, ii1, ii2      DOUBLE PRECISION decalx
40      DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin  
41      DOUBLE PRECISION champmin, champmax, decalx      !----------------------------------------------------------------------
42      INTEGER is2  
43      SAVE is2      print *, "Call sequence information: fxhyp"
44    
45      DOUBLE PRECISION heavyside      xzoom = clon * pi_d / 180d0
46    
47      pi = 2. * ASIN(1.)      IF (grossismx == 1.) THEN
48      depi = 2. * pi         decalx = 1d0
49      epsilon = 1.e-3      else
50      xzoom = xzoomdeg * pi/180.         decalx = 0.75d0
51        END IF
52      decalx = .75  
53      IF(grossism.EQ.1..AND.scal180) THEN      dzoom = dzoomx * twopi_d
54         decalx = 1.      xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
55      ENDIF  
56        ! Compute fhyp:
57      WRITE(6, *) 'FXHYP scal180, decalx', scal180, decalx      DO i = nmax, 2 * nmax
58           fa = taux * (dzoom / 2. - xtild(i))
59      IF(dzooma.LT.1.) THEN         fb = xtild(i) * (pi_d - xtild(i))
60         dzoom = dzooma * depi  
61      ELSEIF(dzooma.LT. 25.) THEN         IF (200. * fb < - fa) THEN
62         WRITE(6, *) ' Le param. dzoomx pour fxhyp est trop petit ! L augmenter et relancer ! '            fhyp(i) = - 1.
63         STOP 1         ELSE IF (200. * fb < fa) THEN
64      ELSE            fhyp(i) = 1.
        dzoom = dzooma * pi/180.  
     ENDIF  
   
     WRITE(6, *) ' xzoom(rad.), grossism, tau, dzoom (radians)'  
     WRITE(6, 24) 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         IF (xtild(i).EQ. 0.) fhyp(i) = 1.  
77         IF (xtild(i).EQ. pi) fhyp(i) = -1.         IF (xtild(i) == 0.) fhyp(i) = 1.
78           IF (xtild(i) == pi_d) fhyp(i) = -1.
79        END DO
80    
81      ENDDO      ! Calcul de beta
   
     !c .... 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.EQ. 0.) fxm = 1.         IF (xmoy == 0.) fxm = 1.
107         IF (xmoy.EQ. 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        END DO
111    
112      ENDDO      print *, "ffdx = ", ffdx
113        beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
114      beta = (grossism * ffdx - pi) / (ffdx - pi)      print *, "beta = ", beta
115    
116      IF(2.*beta - grossism.LE. 0.) THEN      IF (2. * beta - grossismx <= 0.) THEN
117         WRITE(6, *) ' ** Attention ! La valeur beta calculee dans la routine fxhyp est mauvaise ! '         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
118         WRITE(6, *)'Modifier les valeurs de grossismx, tau ou dzoomx ', &         print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
             ' et relancer ! *** '  
119         STOP 1         STOP 1
120      ENDIF      END IF
   
     ! ..... calcul de Xprimt .....  
   
   
     DO i = nmax, nmax2  
        Xprimt(i) = beta + (grossism - beta) * fhyp(i)  
     ENDDO  
121    
122      DO i = nmax+1, nmax2      ! calcul de Xprimt
123         Xprimt(nmax2 - i) = Xprimt(i)      Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
124      ENDDO      xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
125    
126        ! Calcul de Xf
127    
128      ! ..... Calcul de Xf ........      DO i = nmax + 1, 2 * nmax
   
     Xf(0) = - pi  
   
     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.EQ. 0.) fxm = 1.         IF (xmoy == 0.) fxm = 1.
142         IF (xmoy.EQ. pi) fxm = -1.         IF (xmoy == pi_d) fxm = -1.
143         xxpr(i) = beta + (grossism - beta) * fxm         xxpr(i) = beta + (grossismx - beta) * fxm
144        END DO
145    
146      ENDDO      xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
147    
148      DO i = nmax+1, nmax2      Xf(0) = - pi_d
        xxpr(nmax2-i+1) = xxpr(i)  
     ENDDO  
149    
150      DO i=1, nmax2      DO i=1, 2 * nmax - 1
151         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
152      ENDDO      END DO
   
     ! *****************************************************************  
   
153    
154      ! ..... xuv = 0. si calcul aux pts scalaires ........      Xf(2 * nmax) = pi_d
     ! ..... xuv = 0.5 si calcul aux pts U ........  
155    
156      WRITE(6, 18)      call fxhyp_loop_ik(1, decalx, xf, xtild, Xprimt, xzoom, rlonm025, &
157             xprimm025, xuv = - 0.25d0)
158      DO ik = 1, 4      call fxhyp_loop_ik(2, decalx, xf, xtild, Xprimt, xzoom, rlonv, xprimv, &
159             xuv = 0d0)
160         IF(ik.EQ.1) THEN      call fxhyp_loop_ik(3, decalx, xf, xtild, Xprimt, xzoom, rlonu, xprimu, &
161            xuv = -0.25           xuv = 0.5d0)
162         ELSE IF (ik.EQ.2) THEN      call fxhyp_loop_ik(4, decalx, xf, xtild, Xprimt, xzoom, rlonp025, &
163            xuv = 0.           xprimp025, xuv = 0.25d0)
164         ELSE IF (ik.EQ.3) THEN  
165            xuv = 0.50      print *
166         ELSE IF (ik.EQ.4) THEN  
167            xuv = 0.25      forall (i = 1: iim) xlon(i) = rlonv(i + 1) - rlonv(i)
168         ENDIF      print *, "Minimum longitude step:", MINval(xlon) * 180. / pi_d, "°"
169        print *, "Maximum longitude step:", MAXval(xlon) * 180. / pi_d, "°"
170         xo1 = 0.  
171        DO i = 1, iim + 1
172         ii1=1         IF (rlonp025(i) < rlonv(i)) THEN
173         ii2=iim            print *, 'rlonp025(', i, ') = ', rlonp025(i)
174         IF(ik.EQ.1.and.grossism.EQ.1.) THEN            print *, "< rlonv(", i, ") = ", rlonv(i)
175            ii1 = 2            STOP 1
176            ii2 = iim+1         END IF
177         ENDIF  
178         DO i = ii1, ii2         IF (rlonv(i) < rlonm025(i)) THEN
179              print *, 'rlonv(', i, ') = ', rlonv(i)
180            xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)            print *, "< rlonm025(", i, ") = ", rlonm025(i)
181              STOP 1
182            Xfi = xlon2         END IF
183    
184            DO it = nmax2, 0, -1         IF (rlonp025(i) > rlonu(i)) THEN
185               IF(Xfi.GE.Xf(it)) GO TO 350            print *, 'rlonp025(', i, ') = ', rlonp025(i)
186            end DO            print *, "> rlonu(", i, ") = ", rlonu(i)
187              STOP 1
188            it = 0         END IF
189        END DO
 350       CONTINUE  
   
           ! ...... Calcul de Xf(xi) ......  
   
           xi = xtild(it)  
   
           IF(it.EQ.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  
   
           DO iter = 1, 300  
              xi = xi - (Xf1 - Xfi)/ Xprimin  
   
              IF(ABS(xi-xo1).LE.epsilon) GO TO 550  
              xo1 = xi  
              xi2 = xi * xi  
              Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi  
              Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2  
           end DO  
           WRITE(6, *) ' Pas de solution ***** ', i, xlon2, iter  
           STOP 6  
 550       CONTINUE  
   
           xxprim(i) = depi/ (FLOAT(iim) * Xprimin)  
           xvrai(i) = xi + xzoom  
   
        end DO  
   
        IF(ik.EQ.1.and.grossism.EQ.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  
              WRITE(6, *) ' PBS. avec rlonu(', i+1, ') plus petit que rlonu(', i, &  
                   ')'  
              STOP 7  
           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 .GE.-pi-0.10.and.champmax.LE.pi+0.10)) THEN  
           WRITE(6, *) 'Reorganisation des longitudes pour avoir entre - pi', &  
                ' et pi '  
   
           IF(xzoom.LE.0.) THEN  
              IF(ik.EQ. 1) THEN  
                 DO i = 1, iim  
                    IF(xvrai(i).GE. - pi) GO TO 80  
                 ENDDO  
                 WRITE(6, *) ' PBS. 1 ! Xvrai plus petit que - pi ! '  
                 STOP 8  
 80              CONTINUE  
                 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.EQ.1) THEN  
                 DO i = iim, 1, -1  
                    IF(xvrai(i).LE. pi) GO TO 90  
                 ENDDO  
                 WRITE(6, *) ' PBS. 2 ! Xvrai plus grand que pi ! '  
                 STOP 9  
 90              CONTINUE  
                 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.EQ.1) THEN  
           ! WRITE(6, *) ' XLON aux pts. V-0.25 apres (en deg.) '  
           ! WRITE(6, 18)  
           ! WRITE(6, 68) xvrai  
           ! WRITE(6, *) ' XPRIM k ', ik  
           ! WRITE(6, 566) xprimm  
   
           DO i = 1, iim +1  
              rlonm025(i) = xlon(i)  
              xprimm025(i) = xprimm(i)  
           ENDDO  
        ELSE IF(ik.EQ.2) THEN  
           ! WRITE(6, 18)  
           ! WRITE(6, *) ' XLON aux pts. V apres (en deg.) '  
           ! WRITE(6, 68) xvrai  
           ! WRITE(6, *) ' XPRIM k ', ik  
           ! WRITE(6, 566) xprimm  
   
           DO i = 1, iim + 1  
              rlonv(i) = xlon(i)  
              xprimv(i) = xprimm(i)  
           ENDDO  
   
        ELSE IF(ik.EQ.3) THEN  
           ! WRITE(6, 18)  
           ! WRITE(6, *) ' XLON aux pts. U apres (en deg.) '  
           ! WRITE(6, 68) xvrai  
           ! WRITE(6, *) ' XPRIM ik ', ik  
           ! WRITE(6, 566) xprimm  
   
           DO i = 1, iim + 1  
              rlonu(i) = xlon(i)  
              xprimu(i) = xprimm(i)  
           ENDDO  
   
        ELSE IF(ik.EQ.4) THEN  
           ! WRITE(6, 18)  
           ! WRITE(6, *) ' XLON aux pts. V+0.25 apres (en deg.) '  
           ! WRITE(6, 68) xvrai  
           ! WRITE(6, *) ' XPRIM ik ', ik  
           ! WRITE(6, 566) xprimm  
   
           DO i = 1, iim + 1  
              rlonp025(i) = xlon(i)  
              xprimp025(i) = xprimm(i)  
           ENDDO  
   
        ENDIF  
   
     end DO  
   
     WRITE(6, 18)  
   
     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  
   
 18  FORMAT(/)  
 24  FORMAT(2x, 'Parametres xzoom, gross, tau, dzoom pour fxhyp ', 4f8.3)  
 68  FORMAT(1x, 7f9.2)  
 566 FORMAT(1x, 7f9.4)  
190    
191    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
192    
Legend:
Removed from v.82  
changed lines
  Added in v.121
  ViewVC Help
Powered by ViewVC 1.1.21