/[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 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  
19    
20      REAL, intent(in):: tau      USE dimens_m, ONLY: iim
21      ! raideur de la transition de l'intérieur à l'extérieur du zoom      use fxhyp_loop_ik_m, only: fxhyp_loop_ik, nmax
22        use nr_util, only: pi, pi_d, twopi, twopi_d, arth
23      ! arguments de sortie      use principal_cshift_m, only: principal_cshift
24        use serre, only: clon, grossismx, dzoomx, taux
     REAL, dimension(iip1):: rlonm025, xprimm025, rlonv, xprimv  
     real, dimension(iip1):: rlonu, xprimu, rlonp025, xprimp025  
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        REAL dzoom, step
32        real d_rlonv(iim)
33        DOUBLE PRECISION xtild(0:2 * nmax)
34        DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
35        DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
36        DOUBLE PRECISION xzoom, fa, fb
37        INTEGER i, is2
38        DOUBLE PRECISION xmoy, fxm
39    
40      INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2*nmax      !----------------------------------------------------------------------
41    
42      LOGICAL, PARAMETER:: scal180 = .TRUE.      print *, "Call sequence information: fxhyp"
     ! 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.  
   
     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  
43    
44      DOUBLE PRECISION heavyside      xzoom = clon * pi_d / 180d0
45    
46      !----------------------------------------------------------------------      test_grossismx: if (grossismx == 1.) then
47           step = twopi / iim
48    
49      pi = 2. * ASIN(1.)         xprimm025(:iim) = step
50      depi = 2. * pi         xprimp025(:iim) = step
51      epsilon = 1.e-3         xprimv(:iim) = step
52      xzoom = xzoomdeg * pi/180.         xprimu(:iim) = step
53    
54      decalx = .75         rlonv(:iim) = arth(- pi + clon * pi / 180., step, iim)
55      IF (grossism == 1. .AND. scal180) THEN         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
56         decalx = 1.         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
57      ENDIF         rlonu(:iim) = rlonv(:iim) + 0.5 * step
58        else
59      print *, 'FXHYP scal180, decalx', scal180, decalx         dzoom = dzoomx * twopi_d
60           xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
61      IF (dzooma.LT.1.) THEN  
62         dzoom = dzooma * depi         ! Compute fhyp:
63      ELSEIF (dzooma.LT. 25.) THEN         DO i = nmax, 2 * nmax
64         print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &            fa = taux * (dzoom / 2. - xtild(i))
65              // "L'augmenter et relancer."            fb = xtild(i) * (pi_d - xtild(i))
66         STOP 1  
67      ELSE            IF (200. * fb < - fa) THEN
68         dzoom = dzooma * pi/180.               fhyp(i) = - 1.
69      ENDIF            ELSE IF (200. * fb < fa) THEN
70                 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.  
        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  
71            ELSE            ELSE
72               fhyp (i) = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
73            ENDIF                  IF (200. * fb + fa < 1e-10) THEN
74         ENDIF                     fhyp(i) = - 1.
75                    ELSE IF (200. * fb - fa < 1e-10) THEN
76         IF (xtild(i) == 0.) fhyp(i) = 1.                     fhyp(i) = 1.
77         IF (xtild(i) == pi) fhyp(i) = -1.                  END IF
78      ENDDO               ELSE
79                    fhyp(i) = TANH(fa / fb)
80      ! Calcul de beta               END IF
81              END IF
82      ffdx = 0.  
83              IF (xtild(i) == 0.) fhyp(i) = 1.
84      DO i = nmax + 1, nmax2            IF (xtild(i) == pi_d) fhyp(i) = -1.
85         xmoy = 0.5 * (xtild(i-1) + xtild(i))         END DO
86         fa = tau* (dzoom/2. - xmoy)  
87         fb = xmoy * (pi - xmoy)         ! Calcul de beta
88    
89         IF (200.* fb .LT. - fa) THEN         ffdx = 0.
90            fxm = - 1.  
91         ELSEIF (200. * fb .LT. fa) THEN         DO i = nmax + 1, 2 * nmax
92            fxm = 1.            xmoy = 0.5 * (xtild(i-1) + xtild(i))
93         ELSE            fa = taux * (dzoom / 2. - xmoy)
94            IF (ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN            fb = xmoy * (pi_d - xmoy)
95               IF (200.*fb + fa.LT.1.e-10) THEN  
96                  fxm = - 1.            IF (200. * fb < - fa) THEN
97               ELSEIF (200.*fb - fa.LT.1.e-10) THEN               fxm = - 1.
98                  fxm = 1.            ELSE IF (200. * fb < fa) THEN
99               ENDIF               fxm = 1.
100            ELSE            ELSE
101               fxm = TANH (fa/fb)               IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
102            ENDIF                  IF (200. * fb + fa < 1e-10) THEN
103         ENDIF                     fxm = - 1.
104                    ELSE IF (200. * fb - fa < 1e-10) THEN
105         IF (xmoy == 0.) fxm = 1.                     fxm = 1.
106         IF (xmoy == pi) fxm = -1.                  END IF
107                 ELSE
108         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))                  fxm = TANH(fa / fb)
109      ENDDO               END IF
110              END IF
111      beta = (grossism * ffdx - pi) / (ffdx - pi)  
112              IF (xmoy == 0.) fxm = 1.
113      IF (2.*beta - grossism <= 0.) THEN            IF (xmoy == pi_d) fxm = -1.
114         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'  
115         print *, 'Modifier les valeurs de grossismx, tau ou dzoomx et relancer.'            ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
116         STOP 1         END DO
117      ENDIF  
118           print *, "ffdx = ", ffdx
119      ! calcul de Xprimt         beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
120           print *, "beta = ", beta
121      DO i = nmax, nmax2  
122         Xprimt(i) = beta + (grossism - beta) * fhyp(i)         IF (2. * beta - grossismx <= 0.) THEN
123      ENDDO            print *, 'Bad choice of grossismx, taux, dzoomx.'
124              print *, 'Decrease dzoomx or grossismx.'
125      DO i = nmax + 1, nmax2            STOP 1
126         Xprimt(nmax2 - i) = Xprimt(i)         END IF
127      ENDDO  
128           ! calcul de Xprimt
129      ! Calcul de Xf         Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
130           xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
131      Xf(0) = - pi  
132           ! Calcul de Xf
133      DO i = nmax + 1, nmax2  
134         xmoy = 0.5 * (xtild(i-1) + xtild(i))         DO i = nmax + 1, 2 * nmax
135         fa = tau* (dzoom/2. - xmoy)            xmoy = 0.5 * (xtild(i-1) + xtild(i))
136         fb = xmoy * (pi - xmoy)            fa = taux * (dzoom / 2. - xmoy)
137              fb = xmoy * (pi_d - xmoy)
138    
139              IF (200. * fb < - fa) THEN
140                 fxm = - 1.
141              ELSE IF (200. * fb < fa) THEN
142                 fxm = 1.
143              ELSE
144                 fxm = TANH(fa / fb)
145              END IF
146    
147         IF (200.* fb .LT. - fa) THEN            IF (xmoy == 0.) fxm = 1.
148            fxm = - 1.            IF (xmoy == pi_d) fxm = -1.
149         ELSEIF (200. * fb .LT. fa) THEN            xxpr(i) = beta + (grossismx - beta) * fxm
150            fxm = 1.         END DO
151         ELSE  
152            fxm = TANH (fa/fb)         xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
153         ENDIF  
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         IF (xmoy == 0.) fxm = 1.            do while (rlonm025(is2) < - pi .and. is2 < iim)
180         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 + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)  
           Xfi = xlon2  
   
           it = nmax2  
           do while (xfi < xf(it) .and. it >= 1)  
              it = it - 1  
181            end do            end do
182    
183            ! Calcul de Xf(xi)            if (rlonm025(is2) < - pi) then
184                 print *, 'Rlonm025 plus petit 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  
185               STOP 1               STOP 1
186            end if            end if
187           ELSE
188              is2 = iim
189    
190              do while (rlonm025(is2) > pi .and. is2 > 1)
191                 is2 = is2 - 1
192              end do
193    
194            xxprim(i) = depi/ (FLOAT(iim) * Xprimin)            if (rlonm025(is2) > pi) then
195            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, ')'  
196               STOP 1               STOP 1
197            ENDIF            end if
198         ENDDO         END IF
199        END IF
        ! Reorganisation des longitudes pour les avoir entre - pi et pi  
200    
201         champmin = 1.e12      call principal_cshift(is2, rlonm025, xprimm025)
202         champmax = -1.e12      call principal_cshift(is2, rlonv, xprimv)
203         DO i = 1, iim      call principal_cshift(is2, rlonu, xprimu)
204            champmin = MIN(champmin, xvrai(i))      call principal_cshift(is2, rlonp025, xprimp025)
205            champmax = MAX(champmax, xvrai(i))  
206         ENDDO      forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
207        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
208         IF (.not. (champmin >= -pi-0.10.and.champmax <= pi + 0.10)) THEN      print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
209            print *, 'Reorganisation des longitudes pour avoir entre - pi', &  
210                 ' et pi '      DO i = 1, iim + 1
211           IF (rlonp025(i) < rlonv(i)) THEN
212            IF (xzoom <= 0.) THEN            print *, 'rlonp025(', i, ') = ', rlonp025(i)
213               IF (ik == 1) THEN            print *, "< rlonv(", i, ") = ", rlonv(i)
214                  i = 1            STOP 1
215           END IF
216                  do while (xvrai(i) < - pi .and. i < iim)  
217                     i = i + 1         IF (rlonv(i) < rlonm025(i)) THEN
218                  end do            print *, 'rlonv(', i, ') = ', rlonv(i)
219              print *, "< rlonm025(", i, ") = ", rlonm025(i)
220                  if (xvrai(i) < - pi) then            STOP 1
221                     print *, ' PBS. 1 ! Xvrai plus petit que - pi ! '         END IF
222                     STOP 1  
223                  end if         IF (rlonp025(i) > rlonu(i)) THEN
224              print *, 'rlonp025(', i, ') = ', rlonp025(i)
225                  is2 = i            print *, "> rlonu(", i, ") = ", rlonu(i)
226               ENDIF            STOP 1
227           END IF
228               IF (is2.NE. 1) THEN      END DO
                 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.91  
changed lines
  Added in v.126
  ViewVC Help
Powered by ViewVC 1.1.21