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