/[lmdze]/trunk/Sources/dyn3d/fxhyp.f
ViewVC logotype

Diff of /trunk/Sources/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 120 by guez, Tue Jan 13 14:56: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, 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.)  
   
     ! 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.  
14    
15      ! ...... arguments d'entree .......      ! On doit avoir grossismx \times dzoomx < pi (radians)
16    
17      REAL xzoomdeg, dzooma, tau, grossism      ! Le premier point scalaire pour une grille regulière (grossismx =
18      ! grossism etant le grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)      ! 1., taux=0., clon=0.) est à - 180 degrés.
     ! dzooma etant la distance totale de la zone du zoom  
     ! tau la raideur de la transition de l'interieur a l'exterieur du zoom  
19    
20      ! ...... arguments de sortie ......      use coefpoly_m, only: coefpoly
21        USE dimens_m, ONLY: iim
22        use nr_util, only: pi_d, twopi_d, arth
23        use serre, only: clon, grossismx, dzoomx, taux
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        DOUBLE PRECISION champmin, champmax
31        real rlonm025(iim + 1), rlonp025(iim + 1)
32        INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax
33      REAL dzoom      REAL dzoom
34      DOUBLE PRECISION xlon(iip1), xprimm(iip1), xuv      DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv
35      DOUBLE PRECISION xtild(0:nmax2)      DOUBLE PRECISION xtild(0:nmax2)
36      DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)      DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)
37      DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)      DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)
38      DOUBLE PRECISION xvrai(iip1), xxprim(iip1)      DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)
39      DOUBLE PRECISION pi, depi, epsilon, xzoom, fa, fb      DOUBLE PRECISION my_eps, xzoom, fa, fb
40      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2      DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
41      INTEGER i, it, ik, iter, ii, idif, ii1, ii2      INTEGER i, it, ik, iter, ii, idif, ii1, ii2
42      DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin      DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin
43      DOUBLE PRECISION champmin, champmax, decalx      DOUBLE PRECISION decalx
44      INTEGER is2      INTEGER is2
     SAVE is2  
   
     DOUBLE PRECISION heavyside  
45    
46      pi = 2. * ASIN(1.)      !----------------------------------------------------------------------
     depi = 2. * pi  
     epsilon = 1.e-3  
     xzoom = xzoomdeg * pi/180.  
47    
48      decalx = .75      print *, "Call sequence information: fxhyp"
     IF(grossism.EQ.1..AND.scal180) THEN  
        decalx = 1.  
     ENDIF  
49    
50      WRITE(6, *) 'FXHYP scal180, decalx', scal180, decalx      my_eps = 1e-3
51        xzoom = clon * pi_d / 180.
52    
53      IF(dzooma.LT.1.) THEN      IF (grossismx == 1.) THEN
54         dzoom = dzooma * depi         decalx = 1.
55      ELSEIF(dzooma.LT. 25.) THEN      else
56         WRITE(6, *) ' Le param. dzoomx pour fxhyp est trop petit ! L augmenter et relancer ! '         decalx = 0.75
57        END IF
58    
59        IF (dzoomx < 1.) THEN
60           dzoom = dzoomx * twopi_d
61        ELSE IF (dzoomx < 25.) THEN
62           print *, "dzoomx pour fxhyp est trop petit."
63         STOP 1         STOP 1
64      ELSE      ELSE
65         dzoom = dzooma * pi/180.         dzoom = dzoomx * pi_d / 180.
66      ENDIF      END IF
67    
68      WRITE(6, *) ' xzoom(rad.), grossism, tau, dzoom (radians)'      print *, 'dzoom (rad):', dzoom
     WRITE(6, 24) xzoom, grossism, tau, dzoom  
69    
70      DO i = 0, nmax2      xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)
        xtild(i) = - pi + FLOAT(i) * depi /nmax2  
     ENDDO  
71    
72      DO i = nmax, nmax2      DO i = nmax, nmax2
73           fa = taux * (dzoom / 2. - xtild(i))
74           fb = xtild(i) * (pi_d - xtild(i))
75    
76         fa = tau* (dzoom/2. - xtild(i))         IF (200. * fb < - fa) THEN
77         fb = xtild(i) * (pi - xtild(i))            fhyp(i) = - 1.
78           ELSE IF (200. * fb < fa) THEN
79         IF(200.* fb .LT. - fa) THEN            fhyp(i) = 1.
           fhyp (i) = - 1.  
        ELSEIF(200. * fb .LT. fa) THEN  
           fhyp (i) = 1.  
80         ELSE         ELSE
81            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
82               IF(200.*fb + fa.LT.1.e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
83                  fhyp (i) = - 1.                  fhyp(i) = - 1.
84               ELSEIF(200.*fb - fa.LT.1.e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
85                  fhyp (i) = 1.                  fhyp(i) = 1.
86               ENDIF               END IF
87            ELSE            ELSE
88               fhyp (i) = TANH (fa/fb)               fhyp(i) = TANH(fa / fb)
89            ENDIF            END IF
90         ENDIF         END IF
91         IF (xtild(i).EQ. 0.) fhyp(i) = 1.  
92         IF (xtild(i).EQ. pi) fhyp(i) = -1.         IF (xtild(i) == 0.) fhyp(i) = 1.
93           IF (xtild(i) == pi_d) fhyp(i) = -1.
94      ENDDO      END DO
95    
96      !c .... Calcul de beta ....      ! Calcul de beta
97    
98      ffdx = 0.      ffdx = 0.
99    
100      DO i = nmax +1, nmax2      DO i = nmax + 1, nmax2
   
101         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
102         fa = tau* (dzoom/2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
103         fb = xmoy * (pi - xmoy)         fb = xmoy * (pi_d - xmoy)
104    
105         IF(200.* fb .LT. - fa) THEN         IF (200. * fb < - fa) THEN
106            fxm = - 1.            fxm = - 1.
107         ELSEIF(200. * fb .LT. fa) THEN         ELSE IF (200. * fb < fa) THEN
108            fxm = 1.            fxm = 1.
109         ELSE         ELSE
110            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
111               IF(200.*fb + fa.LT.1.e-10) THEN               IF (200. * fb + fa < 1e-10) THEN
112                  fxm = - 1.                  fxm = - 1.
113               ELSEIF(200.*fb - fa.LT.1.e-10) THEN               ELSE IF (200. * fb - fa < 1e-10) THEN
114                  fxm = 1.                  fxm = 1.
115               ENDIF               END IF
116            ELSE            ELSE
117               fxm = TANH (fa/fb)               fxm = TANH(fa / fb)
118            ENDIF            END IF
119         ENDIF         END IF
120    
121         IF (xmoy.EQ. 0.) fxm = 1.         IF (xmoy == 0.) fxm = 1.
122         IF (xmoy.EQ. pi) fxm = -1.         IF (xmoy == pi_d) fxm = -1.
123    
124         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))         ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
125        END DO
126    
127      ENDDO      beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
   
     beta = (grossism * ffdx - pi) / (ffdx - pi)  
128    
129      IF(2.*beta - grossism.LE. 0.) THEN      IF (2. * beta - grossismx <= 0.) THEN
130         WRITE(6, *) ' ** Attention ! La valeur beta calculee dans la routine fxhyp est mauvaise ! '         print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
131         WRITE(6, *)'Modifier les valeurs de grossismx, tau ou dzoomx ', &         print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
             ' et relancer ! *** '  
132         STOP 1         STOP 1
133      ENDIF      END IF
   
     ! ..... calcul de Xprimt .....  
134    
135        ! calcul de Xprimt
136    
137      DO i = nmax, nmax2      DO i = nmax, nmax2
138         Xprimt(i) = beta + (grossism - beta) * fhyp(i)         Xprimt(i) = beta + (grossismx - beta) * fhyp(i)
139      ENDDO      END DO
140    
141      DO i = nmax+1, nmax2      DO i = nmax + 1, nmax2
142         Xprimt(nmax2 - i) = Xprimt(i)         Xprimt(nmax2 - i) = Xprimt(i)
143      ENDDO      END DO
144    
145        ! Calcul de Xf
146    
147      ! ..... Calcul de Xf ........      Xf(0) = - pi_d
   
     Xf(0) = - pi  
   
     DO i = nmax +1, nmax2  
148    
149        DO i = nmax + 1, nmax2
150         xmoy = 0.5 * (xtild(i-1) + xtild(i))         xmoy = 0.5 * (xtild(i-1) + xtild(i))
151         fa = tau* (dzoom/2. - xmoy)         fa = taux * (dzoom / 2. - xmoy)
152         fb = xmoy * (pi - xmoy)         fb = xmoy * (pi_d - xmoy)
153    
154         IF(200.* fb .LT. - fa) THEN         IF (200. * fb < - fa) THEN
155            fxm = - 1.            fxm = - 1.
156         ELSEIF(200. * fb .LT. fa) THEN         ELSE IF (200. * fb < fa) THEN
157            fxm = 1.            fxm = 1.
158         ELSE         ELSE
159            fxm = TANH (fa/fb)            fxm = TANH(fa / fb)
160         ENDIF         END IF
161    
162           IF (xmoy == 0.) fxm = 1.
163           IF (xmoy == pi_d) fxm = -1.
164           xxpr(i) = beta + (grossismx - beta) * fxm
165        END DO
166    
167         IF (xmoy.EQ. 0.) fxm = 1.      xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)
        IF (xmoy.EQ. pi) fxm = -1.  
        xxpr(i) = beta + (grossism - beta) * fxm  
   
     ENDDO  
   
     DO i = nmax+1, nmax2  
        xxpr(nmax2-i+1) = xxpr(i)  
     ENDDO  
168    
169      DO i=1, nmax2      DO i=1, nmax2
170         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))         Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
171      ENDDO      END DO
172    
173      ! *****************************************************************      is2 = 0
174    
175        loop_ik: DO ik = 1, 4
176           ! xuv = 0. si calcul aux points scalaires
177           ! xuv = 0.5 si calcul aux points U
178    
179      ! ..... xuv = 0. si calcul aux pts scalaires ........         IF (ik == 1) THEN
     ! ..... xuv = 0.5 si calcul aux pts U ........  
   
     WRITE(6, 18)  
   
     DO ik = 1, 4  
   
        IF(ik.EQ.1) THEN  
180            xuv = -0.25            xuv = -0.25
181         ELSE IF (ik.EQ.2) THEN         ELSE IF (ik == 2) THEN
182            xuv = 0.            xuv = 0.
183         ELSE IF (ik.EQ.3) THEN         ELSE IF (ik == 3) THEN
184            xuv = 0.50            xuv = 0.50
185         ELSE IF (ik.EQ.4) THEN         ELSE IF (ik == 4) THEN
186            xuv = 0.25            xuv = 0.25
187         ENDIF         END IF
188    
189         xo1 = 0.         xo1 = 0.
190    
191         ii1=1         IF (ik == 1 .and. grossismx == 1.) THEN
        ii2=iim  
        IF(ik.EQ.1.and.grossism.EQ.1.) THEN  
192            ii1 = 2            ii1 = 2
193            ii2 = iim+1            ii2 = iim + 1
194         ENDIF         else
195         DO i = ii1, ii2            ii1=1
196              ii2=iim
197            xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)         END IF
   
           Xfi = xlon2  
198    
199            DO it = nmax2, 0, -1         DO i = ii1, ii2
200               IF(Xfi.GE.Xf(it)) GO TO 350            Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)
           end DO  
   
           it = 0  
201    
202  350       CONTINUE            it = nmax2
203              do while (xfi < xf(it) .and. it >= 1)
204                 it = it - 1
205              end do
206    
207            ! ...... Calcul de Xf(xi) ......            ! Calcul de Xf(xi)
208    
209            xi = xtild(it)            xi = xtild(it)
210    
211            IF(it.EQ.nmax2) THEN            IF (it == nmax2) THEN
212               it = nmax2 -1               it = nmax2 -1
213               Xf(it+1) = pi               Xf(it + 1) = pi_d
214            ENDIF            END IF
           ! .....................................................................  
   
           ! 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))  
215    
216            CALL coefpoly (Xf(it), Xf(it+1), Xprimt(it), Xprimt(it+1), &            ! Appel de la routine qui calcule les coefficients a0, a1,
217                 xtild(it), xtild(it+1), a0, a1, a2, a3)            ! a2, a3 d'un polynome de degre 3 qui passe par les points
218              ! (Xf(it), xtild(it)) et (Xf(it + 1), xtild(it + 1))
219    
220              CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
221                   xtild(it), xtild(it + 1), a0, a1, a2, a3)
222    
223            Xf1 = Xf(it)            Xf1 = Xf(it)
224            Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi            Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi * xi
225    
226            DO iter = 1, 300            iter = 1
              xi = xi - (Xf1 - Xfi)/ Xprimin  
227    
228               IF(ABS(xi-xo1).LE.epsilon) GO TO 550            do
229                 xi = xi - (Xf1 - Xfi) / Xprimin
230                 IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
231               xo1 = xi               xo1 = xi
232               xi2 = xi * xi               xi2 = xi * xi
233               Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi               Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
234               Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2               Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2
235            end DO            end DO
           WRITE(6, *) ' Pas de solution ***** ', i, xlon2, iter  
           STOP 6  
 550       CONTINUE  
236    
237            xxprim(i) = depi/ (FLOAT(iim) * Xprimin)            if (ABS(xi - xo1) > my_eps) then
238            xvrai(i) = xi + xzoom               ! iter == 300
239                 print *, 'Pas de solution.'
240                 print *, i, xfi
241                 STOP 1
242              end if
243    
244              xxprim(i) = twopi_d / (REAL(iim) * Xprimin)
245              xvrai(i) = xi + xzoom
246         end DO         end DO
247    
248         IF(ik.EQ.1.and.grossism.EQ.1.) THEN         IF (ik == 1 .and. grossismx == 1.) THEN
249            xvrai(1) = xvrai(iip1)-depi            xvrai(1) = xvrai(iim + 1)-twopi_d
250            xxprim(1) = xxprim(iip1)            xxprim(1) = xxprim(iim + 1)
251         ENDIF         END IF
252    
253         DO i = 1, iim         DO i = 1, iim
254            xlon(i) = xvrai(i)            xlon(i) = xvrai(i)
255            xprimm(i) = xxprim(i)            xprimm(i) = xxprim(i)
256         ENDDO         END DO
        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  
257    
258         ! ... Reorganisation des longitudes pour les avoir entre - pi et pi ..         DO i = 1, iim -1
259         ! ........................................................................            IF (xvrai(i + 1) < xvrai(i)) THEN
260                 print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'
261                 STOP 1
262              END IF
263           END DO
264    
265           IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 0.1 &
266                .and. MAXval(xvrai(:iim)) <= pi_d + 0.1)) THEN
267              print *, &
268                   'Réorganisation des longitudes pour les avoir entre - pi et pi'
269    
270              IF (xzoom <= 0.) THEN
271                 IF (ik == 1) THEN
272                    i = 1
273    
274                    do while (xvrai(i) < - pi_d .and. i < iim)
275                       i = i + 1
276                    end do
277    
278                    if (xvrai(i) < - pi_d) then
279                       print *, 'Xvrai plus petit que - pi !'
280                       STOP 1
281                    end if
282    
        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  
283                  is2 = i                  is2 = i
284               ENDIF               END IF
285    
286               IF(is2.NE. 1) THEN               IF (is2 /= 1) THEN
287                  DO ii = is2, iim                  DO ii = is2, iim
288                     xlon (ii-is2+1) = xvrai(ii)                     xlon(ii-is2 + 1) = xvrai(ii)
289                     xprimm(ii-is2+1) = xxprim(ii)                     xprimm(ii-is2 + 1) = xxprim(ii)
290                  ENDDO                  END DO
291                  DO ii = 1, is2 -1                  DO ii = 1, is2 -1
292                     xlon (ii+iim-is2+1) = xvrai(ii) + depi                     xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d
293                     xprimm(ii+iim-is2+1) = xxprim(ii)                     xprimm(ii + iim-is2 + 1) = xxprim(ii)
294                  ENDDO                  END DO
295               ENDIF               END IF
296            ELSE            ELSE
297               IF(ik.EQ.1) THEN               IF (ik == 1) THEN
298                  DO i = iim, 1, -1                  i = iim
299                     IF(xvrai(i).LE. pi) GO TO 90  
300                  ENDDO                  do while (xvrai(i) > pi_d .and. i > 1)
301                  WRITE(6, *) ' PBS. 2 ! Xvrai plus grand que pi ! '                     i = i - 1
302                  STOP 9                  end do
303  90              CONTINUE  
304                    if (xvrai(i) > pi_d) then
305                       print *, 'Xvrai plus grand que pi !'
306                       STOP 1
307                    end if
308    
309                  is2 = i                  is2 = i
310               ENDIF               END IF
311    
312               idif = iim -is2               idif = iim -is2
313    
314               DO ii = 1, is2               DO ii = 1, is2
315                  xlon (ii+idif) = xvrai(ii)                  xlon(ii + idif) = xvrai(ii)
316                  xprimm(ii+idif) = xxprim(ii)                  xprimm(ii + idif) = xxprim(ii)
317               ENDDO               END DO
318    
319               DO ii = 1, idif               DO ii = 1, idif
320                  xlon (ii) = xvrai (ii+is2) - depi                  xlon(ii) = xvrai(ii + is2) - twopi_d
321                  xprimm(ii) = xxprim(ii+is2)                  xprimm(ii) = xxprim(ii + is2)
322               ENDDO               END DO
323            ENDIF            END IF
324         ENDIF         END IF
325    
326         ! ......... Fin de la reorganisation ............................         xlon(iim + 1) = xlon(1) + twopi_d
327           xprimm(iim + 1) = xprimm(1)
328         xlon (iip1) = xlon(1) + depi  
329         xprimm(iip1) = xprimm (1)         DO i = 1, iim + 1
330              xvrai(i) = xlon(i) * 180. / pi_d
331         DO i = 1, iim+1         END DO
           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  
332    
333            DO i = 1, iim +1         IF (ik == 1) THEN
334              DO i = 1, iim + 1
335               rlonm025(i) = xlon(i)               rlonm025(i) = xlon(i)
336               xprimm025(i) = xprimm(i)               xprimm025(i) = xprimm(i)
337            ENDDO            END DO
338         ELSE IF(ik.EQ.2) THEN         ELSE IF (ik == 2) THEN
339            ! WRITE(6, 18)            rlonv = xlon
340            ! WRITE(6, *) ' XLON aux pts. V apres (en deg.) '            xprimv = xprimm
341            ! WRITE(6, 68) xvrai         ELSE IF (ik == 3) THEN
           ! 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  
   
342            DO i = 1, iim + 1            DO i = 1, iim + 1
343               rlonu(i) = xlon(i)               rlonu(i) = xlon(i)
344               xprimu(i) = xprimm(i)               xprimu(i) = xprimm(i)
345            ENDDO            END DO
346           ELSE IF (ik == 4) THEN
347         ELSE IF(ik.EQ.4) THEN            rlonp025 = xlon
348            ! WRITE(6, 18)            xprimp025 = xprimm
349            ! WRITE(6, *) ' XLON aux pts. V+0.25 apres (en deg.) '         END IF
350            ! WRITE(6, 68) xvrai      end DO loop_ik
           ! 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  
351    
352      WRITE(6, 18)      print *
353    
354      DO i = 1, iim      DO i = 1, iim
355         xlon(i) = rlonv(i+1) - rlonv(i)         xlon(i) = rlonv(i + 1) - rlonv(i)
356      ENDDO      END DO
357      champmin = 1.e12      champmin = 1e12
358      champmax = -1.e12      champmax = -1e12
359      DO i = 1, iim      DO i = 1, iim
360         champmin = MIN(champmin, xlon(i))         champmin = MIN(champmin, xlon(i))
361         champmax = MAX(champmax, xlon(i))         champmax = MAX(champmax, xlon(i))
362      ENDDO      END DO
363      champmin = champmin * 180./pi      champmin = champmin * 180. / pi_d
364      champmax = champmax * 180./pi      champmax = champmax * 180. / pi_d
365    
366  18  FORMAT(/)      DO i = 1, iim + 1
367  24  FORMAT(2x, 'Parametres xzoom, gross, tau, dzoom pour fxhyp ', 4f8.3)         IF (rlonp025(i) < rlonv(i)) THEN
368  68  FORMAT(1x, 7f9.2)            print *, ' Attention ! rlonp025 < rlonv', i
369  566 FORMAT(1x, 7f9.4)            STOP 1
370           END IF
371    
372           IF (rlonv(i) < rlonm025(i)) THEN
373              print *, ' Attention ! rlonm025 > rlonv', i
374              STOP 1
375           END IF
376    
377           IF (rlonp025(i) > rlonu(i)) THEN
378              print *, 'rlonp025(', i, ') = ', rlonp025(i)
379              print *, "> rlonu(", i, ") = ", rlonu(i)
380              STOP 1
381           END IF
382        END DO
383    
384        print *, ' Longitudes '
385        print 3, champmin, champmax
386    
387    3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
388             ' d environ ', f0.2, ' degres ', /, &
389             ' alors que la maille en dehors de la zone du zoom est ', &
390             "d'environ ", f0.2, ' degres ')
391    
392    END SUBROUTINE fxhyp    END SUBROUTINE fxhyp
393    
Legend:
Removed from v.82  
changed lines
  Added in v.120
  ViewVC Help
Powered by ViewVC 1.1.21