/[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 112 by guez, Thu Sep 18 13:36:51 2014 UTC revision 127 by guez, Tue Feb 10 17:58:56 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 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):'      test_grossismx: if (grossismx == 1.) then
45      print *, xzoom, grossism, tau, dzoom         step = twopi / iim
46    
47      DO i = 0, nmax2         xprimm025(:iim) = step
48         xtild(i) = - pi + REAL(i) * depi /nmax2         xprimp025(:iim) = step
49      ENDDO         xprimv(:iim) = step
50           xprimu(:iim) = step
51      DO i = nmax, nmax2  
52         fa = tau* (dzoom/2. - xtild(i))         rlonv(:iim) = arth(- pi + clon, step, iim)
53         fb = xtild(i) * (pi - xtild(i))         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
54           rlonp025(:iim) = rlonv(:iim) + 0.25 * step
55         IF (200.* fb .LT. - fa) THEN         rlonu(:iim) = rlonv(:iim) + 0.5 * step
56            fhyp (i) = - 1.      else
57         ELSEIF (200. * fb .LT. fa) THEN         dzoom = dzoomx * twopi_d
58            fhyp (i) = 1.         xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
59         ELSE  
60            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN         ! Compute fhyp:
61               IF (200.*fb + fa.LT.1.e-10) THEN         DO i = nmax, 2 * nmax
62                  fhyp (i) = - 1.            fa = taux * (dzoom / 2. - xtild(i))
63               ELSEIF (200.*fb - fa.LT.1.e-10) THEN            fb = xtild(i) * (pi_d - xtild(i))
64                  fhyp (i) = 1.  
65               ENDIF            IF (200. * fb < - fa) THEN
66                 fhyp(i) = - 1.
67              ELSE IF (200. * fb < fa) THEN
68                 fhyp(i) = 1.
69              ELSE
70                 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
71                    IF (200. * fb + fa < 1e-10) THEN
72                       fhyp(i) = - 1.
73                    ELSE IF (200. * fb - fa < 1e-10) THEN
74                       fhyp(i) = 1.
75                    END IF
76                 ELSE
77                    fhyp(i) = TANH(fa / fb)
78                 END IF
79              END IF
80    
81              IF (xtild(i) == 0.) fhyp(i) = 1.
82              IF (xtild(i) == pi_d) fhyp(i) = -1.
83           END DO
84    
85           ! Calcul de beta
86    
87           ffdx = 0.
88    
89           DO i = nmax + 1, 2 * nmax
90              xmoy = 0.5 * (xtild(i-1) + xtild(i))
91              fa = taux * (dzoom / 2. - xmoy)
92              fb = xmoy * (pi_d - xmoy)
93    
94              IF (200. * fb < - fa) THEN
95                 fxm = - 1.
96              ELSE IF (200. * fb < fa) THEN
97                 fxm = 1.
98            ELSE            ELSE
99               fhyp (i) = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
100            ENDIF                  IF (200. * fb + fa < 1e-10) THEN
101                       fxm = - 1.
102                    ELSE IF (200. * fb - fa < 1e-10) THEN
103                       fxm = 1.
104                    END IF
105                 ELSE
106                    fxm = TANH(fa / fb)
107                 END IF
108              END IF
109    
110              IF (xmoy == 0.) fxm = 1.
111              IF (xmoy == pi_d) fxm = -1.
112    
113              ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
114           END DO
115    
116           print *, "ffdx = ", ffdx
117           beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
118           print *, "beta = ", beta
119    
120           IF (2. * beta - grossismx <= 0.) THEN
121              print *, 'Bad choice of grossismx, taux, dzoomx.'
122              print *, 'Decrease dzoomx or grossismx.'
123              STOP 1
124         END IF         END IF
125    
126         IF (xtild(i) == 0.) fhyp(i) = 1.         ! calcul de Xprimt
127         IF (xtild(i) == pi) fhyp(i) = -1.         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
128      END DO         xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
129    
130      ! Calcul de beta         ! Calcul de Xf
131    
132      ffdx = 0.         DO i = nmax + 1, 2 * nmax
133              xmoy = 0.5 * (xtild(i-1) + xtild(i))
134      DO i = nmax + 1, nmax2            fa = taux * (dzoom / 2. - xmoy)
135         xmoy = 0.5 * (xtild(i-1) + xtild(i))            fb = xmoy * (pi_d - xmoy)
136         fa = tau* (dzoom/2. - xmoy)  
137         fb = xmoy * (pi - xmoy)            IF (200. * fb < - fa) THEN
138                 fxm = - 1.
139         IF (200.* fb .LT. - fa) THEN            ELSE IF (200. * fb < fa) THEN
140            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  
141            ELSE            ELSE
142               fxm = TANH (fa/fb)               fxm = TANH(fa / fb)
143            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  
144    
145      ! calcul de Xprimt            IF (xmoy == 0.) fxm = 1.
146              IF (xmoy == pi_d) fxm = -1.
147              xxpr(i) = beta + (grossismx - beta) * fxm
148           END DO
149    
150           xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
151    
152           Xf(0) = - pi_d
153    
154           DO i=1, 2 * nmax - 1
155              Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
156           END DO
157    
158           Xf(2 * nmax) = pi_d
159    
160           call fxhyp_loop_ik(xf, xtild, Xprimt, rlonm025(:iim), xprimm025(:iim), &
161                xuv = - 0.25d0)
162           call fxhyp_loop_ik(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
163                xuv = 0d0)
164           call fxhyp_loop_ik(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
165                xuv = 0.5d0)
166           call fxhyp_loop_ik(xf, xtild, Xprimt, rlonp025(:iim), xprimp025(:iim), &
167                xuv = 0.25d0)
168        end if test_grossismx
169    
170        is2 = 0
171    
172        IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
173             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
174           IF (clon <= 0.) THEN
175              is2 = 1
176    
177      DO i = nmax, nmax2            do while (rlonm025(is2) < - pi .and. is2 < iim)
178         Xprimt(i) = beta + (grossism - beta) * fhyp(i)               is2 = is2 + 1
179      END DO            end do
   
     DO i = nmax + 1, nmax2  
        Xprimt(nmax2 - i) = Xprimt(i)  
     END DO  
   
     ! Calcul de Xf  
   
     Xf(0) = - pi  
180    
181      DO i = nmax + 1, nmax2            if (rlonm025(is2) < - pi) then
182         xmoy = 0.5 * (xtild(i-1) + xtild(i))               print *, 'Rlonm025 plus petit que - pi !'
183         fa = tau* (dzoom/2. - xmoy)               STOP 1
184         fb = xmoy * (pi - xmoy)            end if
   
        IF (200.* fb .LT. - fa) THEN  
           fxm = - 1.  
        ELSEIF (200. * fb .LT. fa) THEN  
           fxm = 1.  
185         ELSE         ELSE
186            fxm = TANH (fa/fb)            is2 = iim
        ENDIF  
187    
188         IF (xmoy == 0.) fxm = 1.            do while (rlonm025(is2) > pi .and. is2 > 1)
189         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  
190            end do            end do
191    
192            ! Calcul de Xf(xi)            if (rlonm025(is2) > pi) then
193                 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  
194               STOP 1               STOP 1
195            end if            end if
196           END IF
197        END IF
198    
199        call principal_cshift(is2, rlonm025, xprimm025)
200        call principal_cshift(is2, rlonv, xprimv)
201        call principal_cshift(is2, rlonu, xprimu)
202        call principal_cshift(is2, rlonp025, xprimp025)
203    
204        forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
205        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
206        print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
207    
208        DO i = 1, iim + 1
209           IF (rlonp025(i) < rlonv(i)) THEN
210              print *, 'rlonp025(', i, ') = ', rlonp025(i)
211              print *, "< rlonv(", i, ") = ", rlonv(i)
212              STOP 1
213           END IF
214    
215            xxprim(i) = depi/ (REAL(iim) * Xprimin)         IF (rlonv(i) < rlonm025(i)) THEN
216            xvrai(i) = xi + xzoom            print *, 'rlonv(', i, ') = ', rlonv(i)
217         end DO            print *, "< rlonm025(", i, ") = ", rlonm025(i)
218              STOP 1
219         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  
220    
221         champmin = 1.e12         IF (rlonp025(i) > rlonu(i)) THEN
222         champmax = -1.e12            print *, 'rlonp025(', i, ') = ', rlonp025(i)
223         DO i = 1, iim            print *, "> rlonu(", i, ") = ", rlonu(i)
224            champmin = MIN(champmin, xvrai(i))            STOP 1
225            champmax = MAX(champmax, xvrai(i))         END IF
226         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  
227    
228    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
229    
Legend:
Removed from v.112  
changed lines
  Added in v.127
  ViewVC Help
Powered by ViewVC 1.1.21