trunk/dyn3d/fxhyp.f

module fxhyp_m

  IMPLICIT NONE

contains

  SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)

    ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
    ! Author: P. Le Van, from formulas by R. Sadourny

    ! Calcule les longitudes et dérivées dans la grille du GCM pour
    ! une fonction f(x) à dérivée tangente hyperbolique.

    ! On doit avoir grossismx \times dzoomx < pi (radians)

    ! Le premier point scalaire pour une grille regulière (grossismx =
    ! 1., taux=0., clon=0.) est à - 180 degrés.

    use coefpoly_m, only: coefpoly
    USE dimens_m, ONLY: iim
    use nr_util, only: pi_d, twopi_d, arth
    use serre, only: clon, grossismx, dzoomx, taux

    REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
    real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)

    ! Local:

    DOUBLE PRECISION champmin, champmax
    real rlonm025(iim + 1), rlonp025(iim + 1)
    INTEGER, PARAMETER:: nmax = 30000, nmax2 = 2 * nmax
    REAL dzoom
    DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv
    DOUBLE PRECISION xtild(0:nmax2)
    DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)
    DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)
    DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1) 
    DOUBLE PRECISION my_eps, 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, fxm, Xprimin
    DOUBLE PRECISION decalx
    INTEGER is2

    !----------------------------------------------------------------------

    print *, "Call sequence information: fxhyp"

    my_eps = 1e-3
    xzoom = clon * pi_d / 180. 

    IF (grossismx == 1.) THEN
       decalx = 1.
    else
       decalx = 0.75
    END IF

    IF (dzoomx < 1.) THEN
       dzoom = dzoomx * twopi_d
    ELSE IF (dzoomx < 25.) THEN
       print *, "dzoomx pour fxhyp est trop petit."
       STOP 1
    ELSE
       dzoom = dzoomx * pi_d / 180.
    END IF

    print *, 'dzoom (rad):', dzoom

    xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)

    DO i = nmax, nmax2
       fa = taux * (dzoom / 2. - xtild(i))
       fb = xtild(i) * (pi_d - xtild(i))

       IF (200. * fb < - fa) THEN
          fhyp(i) = - 1.
       ELSE IF (200. * fb < fa) THEN
          fhyp(i) = 1.
       ELSE
          IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
             IF (200. * fb + fa < 1e-10) THEN
                fhyp(i) = - 1.
             ELSE IF (200. * fb - fa < 1e-10) THEN
                fhyp(i) = 1.
             END IF
          ELSE
             fhyp(i) = TANH(fa / fb)
          END IF
       END IF

       IF (xtild(i) == 0.) fhyp(i) = 1.
       IF (xtild(i) == pi_d) fhyp(i) = -1.
    END DO

    ! Calcul de beta 

    ffdx = 0.

    DO i = nmax + 1, nmax2
       xmoy = 0.5 * (xtild(i-1) + xtild(i))
       fa = taux * (dzoom / 2. - xmoy)
       fb = xmoy * (pi_d - xmoy)

       IF (200. * fb < - fa) THEN
          fxm = - 1.
       ELSE IF (200. * fb < fa) THEN
          fxm = 1.
       ELSE
          IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
             IF (200. * fb + fa < 1e-10) THEN
                fxm = - 1.
             ELSE IF (200. * fb - fa < 1e-10) THEN
                fxm = 1.
             END IF
          ELSE
             fxm = TANH(fa / fb)
          END IF
       END IF

       IF (xmoy == 0.) fxm = 1.
       IF (xmoy == pi_d) fxm = -1.

       ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
    END DO

    beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)

    IF (2. * beta - grossismx <= 0.) THEN
       print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
       print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
       STOP 1
    END IF

    ! calcul de Xprimt 

    DO i = nmax, nmax2
       Xprimt(i) = beta + (grossismx - beta) * fhyp(i)
    END DO

    DO i = nmax + 1, nmax2
       Xprimt(nmax2 - i) = Xprimt(i)
    END DO

    ! Calcul de Xf 

    Xf(0) = - pi_d

    DO i = nmax + 1, nmax2
       xmoy = 0.5 * (xtild(i-1) + xtild(i))
       fa = taux * (dzoom / 2. - xmoy)
       fb = xmoy * (pi_d - xmoy)

       IF (200. * fb < - fa) THEN
          fxm = - 1.
       ELSE IF (200. * fb < fa) THEN
          fxm = 1.
       ELSE
          fxm = TANH(fa / fb)
       END IF

       IF (xmoy == 0.) fxm = 1.
       IF (xmoy == pi_d) fxm = -1.
       xxpr(i) = beta + (grossismx - beta) * fxm
    END DO

    xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)

    DO i=1, nmax2
       Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
    END DO

    is2 = 0

    loop_ik: DO ik = 1, 4
       ! xuv = 0. si calcul aux points scalaires 
       ! xuv = 0.5 si calcul aux points U 

       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
       END IF

       xo1 = 0.

       IF (ik == 1 .and. grossismx == 1.) THEN
          ii1 = 2 
          ii2 = iim + 1
       else
          ii1=1
          ii2=iim
       END IF

       DO i = ii1, ii2
          Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)

          it = nmax2
          do while (xfi < xf(it) .and. it >= 1)
             it = it - 1
          end do

          ! Calcul de Xf(xi) 

          xi = xtild(it)

          IF (it == nmax2) THEN
             it = nmax2 -1
             Xf(it + 1) = pi_d
          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))

          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) <= my_eps .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) > my_eps) then
             ! iter == 300
             print *, 'Pas de solution.'
             print *, i, xfi
             STOP 1
          end if

          xxprim(i) = twopi_d / (REAL(iim) * Xprimin)
          xvrai(i) = xi + xzoom
       end DO

       IF (ik == 1 .and. grossismx == 1.) THEN
          xvrai(1) = xvrai(iim + 1)-twopi_d
          xxprim(1) = xxprim(iim + 1)
       END IF

       DO i = 1, iim
          xlon(i) = xvrai(i)
          xprimm(i) = xxprim(i)
       END DO

       DO i = 1, iim -1
          IF (xvrai(i + 1) < xvrai(i)) THEN
             print *, 'rlonu(', i + 1, ') < rlonu(', i, ')'
             STOP 1
          END IF
       END DO

       IF (.not. (MINval(xvrai(:iim)) >= - pi_d - 0.1 &
            .and. MAXval(xvrai(:iim)) <= pi_d + 0.1)) THEN
          print *, &
               'Réorganisation des longitudes pour les avoir entre - pi et pi'

          IF (xzoom <= 0.) THEN
             IF (ik == 1) THEN
                i = 1

                do while (xvrai(i) < - pi_d .and. i < iim)
                   i = i + 1
                end do

                if (xvrai(i) < - pi_d) then
                   print *, 'Xvrai plus petit que - pi !'
                   STOP 1
                end if

                is2 = i
             END IF

             IF (is2 /= 1) THEN
                DO ii = is2, iim
                   xlon(ii-is2 + 1) = xvrai(ii)
                   xprimm(ii-is2 + 1) = xxprim(ii)
                END DO
                DO ii = 1, is2 -1
                   xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d
                   xprimm(ii + iim-is2 + 1) = xxprim(ii) 
                END DO
             END IF
          ELSE 
             IF (ik == 1) THEN
                i = iim

                do while (xvrai(i) > pi_d .and. i > 1)
                   i = i - 1
                end do

                if (xvrai(i) > pi_d) then
                   print *, 'Xvrai plus grand que pi !'
                   STOP 1
                end if

                is2 = i
             END IF

             idif = iim -is2

             DO ii = 1, is2
                xlon(ii + idif) = xvrai(ii)
                xprimm(ii + idif) = xxprim(ii)
             END DO

             DO ii = 1, idif
                xlon(ii) = xvrai(ii + is2) - twopi_d
                xprimm(ii) = xxprim(ii + is2) 
             END DO
          END IF
       END IF

       xlon(iim + 1) = xlon(1) + twopi_d
       xprimm(iim + 1) = xprimm(1)

       DO i = 1, iim + 1
          xvrai(i) = xlon(i) * 180. / pi_d
       END DO

       IF (ik == 1) THEN
          DO i = 1, iim + 1
             rlonm025(i) = xlon(i)
             xprimm025(i) = xprimm(i)
          END DO
       ELSE IF (ik == 2) THEN
          rlonv = xlon
          xprimv = xprimm
       ELSE IF (ik == 3) THEN
          DO i = 1, iim + 1
             rlonu(i) = xlon(i)
             xprimu(i) = xprimm(i)
          END DO
       ELSE IF (ik == 4) THEN
          rlonp025 = xlon
          xprimp025 = xprimm
       END IF
    end DO loop_ik

    print *

    DO i = 1, iim
       xlon(i) = rlonv(i + 1) - rlonv(i)
    END DO
    champmin = 1e12
    champmax = -1e12
    DO i = 1, iim
       champmin = MIN(champmin, xlon(i))
       champmax = MAX(champmax, xlon(i))
    END DO
    champmin = champmin * 180. / pi_d
    champmax = champmax * 180. / pi_d

    DO i = 1, iim + 1
       IF (rlonp025(i) < rlonv(i)) THEN
          print *, ' Attention ! rlonp025 < rlonv', i
          STOP 1
       END IF

       IF (rlonv(i) < rlonm025(i)) THEN 
          print *, ' Attention ! rlonm025 > rlonv', i
          STOP 1
       END IF

       IF (rlonp025(i) > rlonu(i)) THEN
          print *, 'rlonp025(', i, ') = ', rlonp025(i)
          print *, "> rlonu(", i, ") = ", rlonu(i)
          STOP 1
       END IF
    END DO

    print *, ' Longitudes '
    print 3, champmin, champmax

3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
         ' d environ ', f0.2, ' degres ', /, &
         ' alors que la maille en dehors de la zone du zoom est ', &
         "d'environ ", f0.2, ' degres ')

  END SUBROUTINE fxhyp

end module fxhyp_m
1	module fxhyp_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8
9	! 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
11
12	! Calcule les longitudes et dérivées dans la grille du GCM pour
13	! une fonction f(x) à dérivée tangente hyperbolique.
14
15	! On doit avoir grossismx \times dzoomx < pi (radians)
16
17	! Le premier point scalaire pour une grille regulière (grossismx =
18	! 1., taux=0., clon=0.) est à - 180 degrés.
19
20	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, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
26	real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
27
28	! 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
34	DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv
35	DOUBLE PRECISION xtild(0:nmax2)
36	DOUBLE PRECISION fhyp(nmax:nmax2), ffdx, beta, Xprimt(0:nmax2)
37	DOUBLE PRECISION Xf(0:nmax2), xxpr(nmax2)
38	DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)
39	DOUBLE PRECISION my_eps, xzoom, fa, fb
40	DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
41	INTEGER i, it, ik, iter, ii, idif, ii1, ii2
42	DOUBLE PRECISION xi, xo1, xmoy, fxm, Xprimin
43	DOUBLE PRECISION decalx
44	INTEGER is2
45
46	!----------------------------------------------------------------------
47
48	print *, "Call sequence information: fxhyp"
49
50	my_eps = 1e-3
51	xzoom = clon * pi_d / 180.
52
53	IF (grossismx == 1.) THEN
54	decalx = 1.
55	else
56	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
64	ELSE
65	dzoom = dzoomx * pi_d / 180.
66	END IF
67
68	print *, 'dzoom (rad):', dzoom
69
70	xtild = arth(- pi_d, twopi_d / nmax2, nmax2 + 1)
71
72	DO i = nmax, nmax2
73	fa = taux * (dzoom / 2. - xtild(i))
74	fb = xtild(i) * (pi_d - xtild(i))
75
76	IF (200. * fb < - fa) THEN
77	fhyp(i) = - 1.
78	ELSE IF (200. * fb < fa) THEN
79	fhyp(i) = 1.
80	ELSE
81	IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
82	IF (200. * fb + fa < 1e-10) THEN
83	fhyp(i) = - 1.
84	ELSE IF (200. * fb - fa < 1e-10) THEN
85	fhyp(i) = 1.
86	END IF
87	ELSE
88	fhyp(i) = TANH(fa / fb)
89	END IF
90	END IF
91
92	IF (xtild(i) == 0.) fhyp(i) = 1.
93	IF (xtild(i) == pi_d) fhyp(i) = -1.
94	END DO
95
96	! Calcul de beta
97
98	ffdx = 0.
99
100	DO i = nmax + 1, nmax2
101	xmoy = 0.5 * (xtild(i-1) + xtild(i))
102	fa = taux * (dzoom / 2. - xmoy)
103	fb = xmoy * (pi_d - xmoy)
104
105	IF (200. * fb < - fa) THEN
106	fxm = - 1.
107	ELSE IF (200. * fb < fa) THEN
108	fxm = 1.
109	ELSE
110	IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
111	IF (200. * fb + fa < 1e-10) THEN
112	fxm = - 1.
113	ELSE IF (200. * fb - fa < 1e-10) THEN
114	fxm = 1.
115	END IF
116	ELSE
117	fxm = TANH(fa / fb)
118	END IF
119	END IF
120
121	IF (xmoy == 0.) fxm = 1.
122	IF (xmoy == pi_d) fxm = -1.
123
124	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
125	END DO
126
127	beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
128
129	IF (2. * beta - grossismx <= 0.) THEN
130	print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
131	print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
132	STOP 1
133	END IF
134
135	! calcul de Xprimt
136
137	DO i = nmax, nmax2
138	Xprimt(i) = beta + (grossismx - beta) * fhyp(i)
139	END DO
140
141	DO i = nmax + 1, nmax2
142	Xprimt(nmax2 - i) = Xprimt(i)
143	END DO
144
145	! Calcul de Xf
146
147	Xf(0) = - pi_d
148
149	DO i = nmax + 1, nmax2
150	xmoy = 0.5 * (xtild(i-1) + xtild(i))
151	fa = taux * (dzoom / 2. - xmoy)
152	fb = xmoy * (pi_d - xmoy)
153
154	IF (200. * fb < - fa) THEN
155	fxm = - 1.
156	ELSE IF (200. * fb < fa) THEN
157	fxm = 1.
158	ELSE
159	fxm = TANH(fa / fb)
160	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	xxpr(:nmax) = xxpr(nmax2:nmax + 1:- 1)
168
169	DO i=1, nmax2
170	Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
171	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	IF (ik == 1) THEN
180	xuv = -0.25
181	ELSE IF (ik == 2) THEN
182	xuv = 0.
183	ELSE IF (ik == 3) THEN
184	xuv = 0.50
185	ELSE IF (ik == 4) THEN
186	xuv = 0.25
187	END IF
188
189	xo1 = 0.
190
191	IF (ik == 1 .and. grossismx == 1.) THEN
192	ii1 = 2
193	ii2 = iim + 1
194	else
195	ii1=1
196	ii2=iim
197	END IF
198
199	DO i = ii1, ii2
200	Xfi = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)
201
202	it = nmax2
203	do while (xfi < xf(it) .and. it >= 1)
204	it = it - 1
205	end do
206
207	! Calcul de Xf(xi)
208
209	xi = xtild(it)
210
211	IF (it == nmax2) THEN
212	it = nmax2 -1
213	Xf(it + 1) = pi_d
214	END IF
215
216	! Appel de la routine qui calcule les coefficients a0, a1,
217	! 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)
224	Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi * xi
225
226	iter = 1
227
228	do
229	xi = xi - (Xf1 - Xfi) / Xprimin
230	IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
231	xo1 = xi
232	xi2 = xi * xi
233	Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
234	Xprimin = a1 + 2. * a2 * xi + 3. * a3 * xi2
235	end DO
236
237	if (ABS(xi - xo1) > my_eps) then
238	! 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
247
248	IF (ik == 1 .and. grossismx == 1.) THEN
249	xvrai(1) = xvrai(iim + 1)-twopi_d
250	xxprim(1) = xxprim(iim + 1)
251	END IF
252
253	DO i = 1, iim
254	xlon(i) = xvrai(i)
255	xprimm(i) = xxprim(i)
256	END DO
257
258	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
283	is2 = i
284	END IF
285
286	IF (is2 /= 1) THEN
287	DO ii = is2, iim
288	xlon(ii-is2 + 1) = xvrai(ii)
289	xprimm(ii-is2 + 1) = xxprim(ii)
290	END DO
291	DO ii = 1, is2 -1
292	xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d
293	xprimm(ii + iim-is2 + 1) = xxprim(ii)
294	END DO
295	END IF
296	ELSE
297	IF (ik == 1) THEN
298	i = iim
299
300	do while (xvrai(i) > pi_d .and. i > 1)
301	i = i - 1
302	end do
303
304	if (xvrai(i) > pi_d) then
305	print *, 'Xvrai plus grand que pi !'
306	STOP 1
307	end if
308
309	is2 = i
310	END IF
311
312	idif = iim -is2
313
314	DO ii = 1, is2
315	xlon(ii + idif) = xvrai(ii)
316	xprimm(ii + idif) = xxprim(ii)
317	END DO
318
319	DO ii = 1, idif
320	xlon(ii) = xvrai(ii + is2) - twopi_d
321	xprimm(ii) = xxprim(ii + is2)
322	END DO
323	END IF
324	END IF
325
326	xlon(iim + 1) = xlon(1) + twopi_d
327	xprimm(iim + 1) = xprimm(1)
328
329	DO i = 1, iim + 1
330	xvrai(i) = xlon(i) * 180. / pi_d
331	END DO
332
333	IF (ik == 1) THEN
334	DO i = 1, iim + 1
335	rlonm025(i) = xlon(i)
336	xprimm025(i) = xprimm(i)
337	END DO
338	ELSE IF (ik == 2) THEN
339	rlonv = xlon
340	xprimv = xprimm
341	ELSE IF (ik == 3) THEN
342	DO i = 1, iim + 1
343	rlonu(i) = xlon(i)
344	xprimu(i) = xprimm(i)
345	END DO
346	ELSE IF (ik == 4) THEN
347	rlonp025 = xlon
348	xprimp025 = xprimm
349	END IF
350	end DO loop_ik
351
352	print *
353
354	DO i = 1, iim
355	xlon(i) = rlonv(i + 1) - rlonv(i)
356	END DO
357	champmin = 1e12
358	champmax = -1e12
359	DO i = 1, iim
360	champmin = MIN(champmin, xlon(i))
361	champmax = MAX(champmax, xlon(i))
362	END DO
363	champmin = champmin * 180. / pi_d
364	champmax = champmax * 180. / pi_d
365
366	DO i = 1, iim + 1
367	IF (rlonp025(i) < rlonv(i)) THEN
368	print *, ' Attention ! rlonp025 < rlonv', i
369	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
393
394	end module fxhyp_m