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