Sources/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)

    USE dimens_m, ONLY: iim
    use nr_util, only: pi_d, twopi_d
    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

    LOGICAL, PARAMETER:: scal180 = .TRUE.
    ! scal180 = .TRUE. si on veut avoir le premier point scalaire pour
    ! une grille reguliere (grossismx = 1., taux=0., clon=0.) a
    ! -180. degres. sinon scal180 = .FALSE.

    REAL dzoom
    DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), 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(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, xlon2, fxm, Xprimin
    DOUBLE PRECISION decalx
    INTEGER, save:: is2

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

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

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

    IF (dzoomx < 1.) THEN
       dzoom = dzoomx * twopi_d
    ELSE IF (dzoomx < 25.) THEN
       print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &
            // "L'augmenter et relancer."
       STOP 1
    ELSE
       dzoom = dzoomx * pi_d / 180.
    END IF

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

    DO i = 0, nmax2 
       xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2
    END DO

    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

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

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

    ! xuv = 0. si calcul aux points scalaires 
    ! xuv = 0.5 si calcul aux points U 

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

       xo1 = 0.

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

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

          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, xlon2
             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 *, 'Problème avec rlonu(', i + 1, &
                  ') plus petit que rlonu(', i, ')'
             STOP 1
          END IF
       END DO

       ! Réorganisation des longitudes pour les avoir entre - pi et pi 

       champmin = 1e12
       champmax = -1e12
       DO i = 1, iim
          champmin = MIN(champmin, xvrai(i))
          champmax = MAX(champmax, xvrai(i))
       END DO

       IF (.not. (champmin >= -pi_d - 0.1 .and. champmax <= pi_d + 0.1)) THEN
          print *, 'Reorganisation des longitudes pour 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.NE. 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

       ! Fin de la reorganisation 

       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
          DO i = 1, iim + 1
             rlonp025(i) = xlon(i)
             xprimp025(i) = xprimm(i)
          END DO
       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 *, ' Attention ! rlonp025 > 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	guez	78	module fxhyp_m
2	guez	3
3	guez	78	IMPLICIT NONE
4	guez	3
5	guez	78	contains
6	guez	3
7	guez	119	SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8	guez	3
9	guez	91	! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10	guez	119	! Author: P. Le Van, from formulas by R. Sadourny
11	guez	3
12	guez	78	! Calcule les longitudes et dérivées dans la grille du GCM pour
13	guez	119	! une fonction f(x) à dérivée tangente hyperbolique.
14	guez	3
15	guez	119	! On doit avoir grossismx \times dzoomx < pi (radians)
16	guez	3
17	guez	78	USE dimens_m, ONLY: iim
18	guez	119	use nr_util, only: pi_d, twopi_d
19			use serre, only: clon, grossismx, dzoomx, taux
20	guez	3
21	guez	119	REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
22			real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
23	guez	3
24	guez	91	! Local:
25	guez	3
26	guez	119	DOUBLE PRECISION champmin, champmax
27			real rlonm025(iim + 1), rlonp025(iim + 1)
28	guez	91	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	guez	119	! une grille reguliere (grossismx = 1., taux=0., clon=0.) a
33	guez	91	! -180. degres. sinon scal180 = .FALSE.
34
35	guez	78	REAL dzoom
36	guez	119	DOUBLE PRECISION xlon(iim + 1), xprimm(iim + 1), xuv
37	guez	78	DOUBLE PRECISION xtild(0:nmax2)
38			DOUBLE PRECISION fhyp(0:nmax2), ffdx, beta, Xprimt(0:nmax2)
39			DOUBLE PRECISION Xf(0:nmax2), xxpr(0:nmax2)
40	guez	119	DOUBLE PRECISION xvrai(iim + 1), xxprim(iim + 1)
41			DOUBLE PRECISION my_eps, xzoom, fa, fb
42	guez	78	DOUBLE PRECISION Xf1, Xfi, a0, a1, a2, a3, xi2
43			INTEGER i, it, ik, iter, ii, idif, ii1, ii2
44			DOUBLE PRECISION xi, xo1, xmoy, xlon2, fxm, Xprimin
45	guez	91	DOUBLE PRECISION decalx
46	guez	119	INTEGER, save:: is2
47	guez	3
48	guez	91	!----------------------------------------------------------------------
49
50	guez	119	my_eps = 1e-3
51			xzoom = clon * pi_d / 180.
52	guez	3
53	guez	119	IF (grossismx == 1. .AND. scal180) THEN
54	guez	78	decalx = 1.
55	guez	119	else
56			decalx = 0.75
57			END IF
58	guez	3
59	guez	119	IF (dzoomx < 1.) THEN
60			dzoom = dzoomx * twopi_d
61			ELSE IF (dzoomx < 25.) THEN
62	guez	91	print *, "Le paramètre dzoomx pour fxhyp est trop petit. " &
63			// "L'augmenter et relancer."
64	guez	78	STOP 1
65			ELSE
66	guez	119	dzoom = dzoomx * pi_d / 180.
67	guez	112	END IF
68	guez	3
69	guez	119	print *, 'dzoom (rad):', dzoom
70	guez	3
71	guez	78	DO i = 0, nmax2
72	guez	119	xtild(i) = - pi_d + REAL(i) * twopi_d / nmax2
73			END DO
74	guez	78
75			DO i = nmax, nmax2
76	guez	119	fa = taux* (dzoom / 2. - xtild(i))
77			fb = xtild(i) * (pi_d - xtild(i))
78	guez	78
79	guez	119	IF (200.* fb < - fa) THEN
80			fhyp(i) = - 1.
81			ELSE IF (200. * fb < fa) THEN
82			fhyp(i) = 1.
83	guez	3	ELSE
84	guez	119	IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
85			IF (200.*fb + fa < 1e-10) THEN
86			fhyp(i) = - 1.
87			ELSE IF (200.*fb - fa < 1e-10) THEN
88			fhyp(i) = 1.
89			END IF
90	guez	78	ELSE
91	guez	119	fhyp(i) = TANH(fa / fb)
92			END IF
93	guez	112	END IF
94	guez	3
95	guez	91	IF (xtild(i) == 0.) fhyp(i) = 1.
96	guez	119	IF (xtild(i) == pi_d) fhyp(i) = -1.
97	guez	112	END DO
98	guez	3
99	guez	91	! Calcul de beta
100	guez	3
101	guez	78	ffdx = 0.
102	guez	3
103	guez	91	DO i = nmax + 1, nmax2
104	guez	78	xmoy = 0.5 * (xtild(i-1) + xtild(i))
105	guez	119	fa = taux* (dzoom / 2. - xmoy)
106			fb = xmoy * (pi_d - xmoy)
107	guez	78
108	guez	119	IF (200.* fb < - fa) THEN
109	guez	78	fxm = - 1.
110	guez	119	ELSE IF (200. * fb < fa) THEN
111	guez	78	fxm = 1.
112			ELSE
113	guez	119	IF (ABS(fa) < 1e-13.AND.ABS(fb) < 1e-13) THEN
114			IF (200.*fb + fa < 1e-10) THEN
115	guez	78	fxm = - 1.
116	guez	119	ELSE IF (200.*fb - fa < 1e-10) THEN
117	guez	78	fxm = 1.
118	guez	119	END IF
119	guez	78	ELSE
120	guez	119	fxm = TANH(fa / fb)
121			END IF
122			END IF
123	guez	3
124	guez	91	IF (xmoy == 0.) fxm = 1.
125	guez	119	IF (xmoy == pi_d) fxm = -1.
126	guez	3
127	guez	78	ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
128	guez	119	END DO
129	guez	3
130	guez	119	beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
131	guez	3
132	guez	119	IF (2. * beta - grossismx <= 0.) THEN
133	guez	91	print *, 'Attention ! La valeur beta calculée dans fxhyp est mauvaise.'
134	guez	119	print *, 'Modifier les valeurs de grossismx, taux ou dzoomx et relancer.'
135	guez	78	STOP 1
136	guez	112	END IF
137	guez	78
138	guez	91	! calcul de Xprimt
139	guez	78
140			DO i = nmax, nmax2
141	guez	119	Xprimt(i) = beta + (grossismx - beta) * fhyp(i)
142	guez	112	END DO
143	guez	78
144	guez	91	DO i = nmax + 1, nmax2
145	guez	78	Xprimt(nmax2 - i) = Xprimt(i)
146	guez	112	END DO
147	guez	78
148	guez	91	! Calcul de Xf
149	guez	78
150	guez	119	Xf(0) = - pi_d
151	guez	78
152	guez	91	DO i = nmax + 1, nmax2
153	guez	78	xmoy = 0.5 * (xtild(i-1) + xtild(i))
154	guez	119	fa = taux* (dzoom / 2. - xmoy)
155			fb = xmoy * (pi_d - xmoy)
156	guez	78
157	guez	119	IF (200.* fb < - fa) THEN
158	guez	78	fxm = - 1.
159	guez	119	ELSE IF (200. * fb < fa) THEN
160	guez	78	fxm = 1.
161	guez	3	ELSE
162	guez	119	fxm = TANH(fa / fb)
163			END IF
164	guez	3
165	guez	91	IF (xmoy == 0.) fxm = 1.
166	guez	119	IF (xmoy == pi_d) fxm = -1.
167			xxpr(i) = beta + (grossismx - beta) * fxm
168			END DO
169	guez	3
170	guez	91	DO i = nmax + 1, nmax2
171			xxpr(nmax2-i + 1) = xxpr(i)
172	guez	119	END DO
173	guez	3
174	guez	78	DO i=1, nmax2
175			Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
176	guez	119	END DO
177	guez	3
178	guez	119	! xuv = 0. si calcul aux points scalaires
179			! xuv = 0.5 si calcul aux points U
180	guez	3
181	guez	119	loop_ik: DO ik = 1, 4
182	guez	91	IF (ik == 1) THEN
183	guez	78	xuv = -0.25
184	guez	91	ELSE IF (ik == 2) THEN
185	guez	78	xuv = 0.
186	guez	91	ELSE IF (ik == 3) THEN
187	guez	78	xuv = 0.50
188	guez	91	ELSE IF (ik == 4) THEN
189	guez	78	xuv = 0.25
190	guez	119	END IF
191	guez	3
192	guez	78	xo1 = 0.
193	guez	3
194	guez	78	ii1=1
195			ii2=iim
196	guez	119	IF (ik == 1.and.grossismx == 1.) THEN
197	guez	78	ii1 = 2
198	guez	91	ii2 = iim + 1
199	guez	119	END IF
200	guez	91
201	guez	78	DO i = ii1, ii2
202	guez	119	xlon2 = - pi_d + (REAL(i) + xuv - decalx) * twopi_d / REAL(iim)
203	guez	78	Xfi = xlon2
204	guez	3
205	guez	91	it = nmax2
206			do while (xfi < xf(it) .and. it >= 1)
207			it = it - 1
208			end do
209	guez	3
210	guez	91	! Calcul de Xf(xi)
211	guez	3
212	guez	78	xi = xtild(it)
213	guez	3
214	guez	91	IF (it == nmax2) THEN
215	guez	78	it = nmax2 -1
216	guez	119	Xf(it + 1) = pi_d
217			END IF
218	guez	3
219	guez	91	! Appel de la routine qui calcule les coefficients a0, a1,
220			! 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	guez	3
223	guez	91	CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), &
224			xtild(it), xtild(it + 1), a0, a1, a2, a3)
225	guez	78
226			Xf1 = Xf(it)
227			Xprimin = a1 + 2.* a2 * xi + 3.a3 xi *xi
228
229	guez	91	iter = 1
230
231			do
232	guez	119	xi = xi - (Xf1 - Xfi) / Xprimin
233			IF (ABS(xi - xo1) <= my_eps .or. iter == 300) exit
234	guez	78	xo1 = xi
235			xi2 = xi * xi
236			Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
237			Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
238			end DO
239	guez	3
240	guez	119	if (ABS(xi - xo1) > my_eps) then
241	guez	91	! iter == 300
242			print *, 'Pas de solution.'
243			print *, i, xlon2
244			STOP 1
245			end if
246
247	guez	119	xxprim(i) = twopi_d / (REAL(iim) * Xprimin)
248	guez	78	xvrai(i) = xi + xzoom
249			end DO
250	guez	3
251	guez	119	IF (ik == 1 .and. grossismx == 1.) THEN
252			xvrai(1) = xvrai(iim + 1)-twopi_d
253			xxprim(1) = xxprim(iim + 1)
254			END IF
255
256	guez	78	DO i = 1, iim
257			xlon(i) = xvrai(i)
258			xprimm(i) = xxprim(i)
259	guez	119	END DO
260
261	guez	3	DO i = 1, iim -1
262	guez	119	IF (xvrai(i + 1) < xvrai(i)) THEN
263	guez	91	print *, 'Problème avec rlonu(', i + 1, &
264			') plus petit que rlonu(', i, ')'
265			STOP 1
266	guez	119	END IF
267			END DO
268	guez	3
269	guez	119	! Réorganisation des longitudes pour les avoir entre - pi et pi
270	guez	78
271	guez	119	champmin = 1e12
272			champmax = -1e12
273	guez	3	DO i = 1, iim
274	guez	78	champmin = MIN(champmin, xvrai(i))
275			champmax = MAX(champmax, xvrai(i))
276	guez	119	END DO
277	guez	3
278	guez	119	IF (.not. (champmin >= -pi_d - 0.1 .and. champmax <= pi_d + 0.1)) THEN
279	guez	91	print *, 'Reorganisation des longitudes pour avoir entre - pi', &
280	guez	78	' et pi '
281	guez	3
282	guez	91	IF (xzoom <= 0.) THEN
283			IF (ik == 1) THEN
284			i = 1
285
286	guez	119	do while (xvrai(i) < - pi_d .and. i < iim)
287	guez	91	i = i + 1
288			end do
289
290	guez	119	if (xvrai(i) < - pi_d) then
291			print *, 'Xvrai plus petit que - pi !'
292	guez	91	STOP 1
293			end if
294
295	guez	78	is2 = i
296	guez	119	END IF
297	guez	78
298	guez	91	IF (is2.NE. 1) THEN
299	guez	78	DO ii = is2, iim
300	guez	119	xlon(ii-is2 + 1) = xvrai(ii)
301	guez	91	xprimm(ii-is2 + 1) = xxprim(ii)
302	guez	119	END DO
303	guez	78	DO ii = 1, is2 -1
304	guez	119	xlon(ii + iim-is2 + 1) = xvrai(ii) + twopi_d
305	guez	91	xprimm(ii + iim-is2 + 1) = xxprim(ii)
306	guez	119	END DO
307			END IF
308	guez	78	ELSE
309	guez	91	IF (ik == 1) THEN
310			i = iim
311
312	guez	119	do while (xvrai(i) > pi_d .and. i > 1)
313	guez	91	i = i - 1
314			end do
315
316	guez	119	if (xvrai(i) > pi_d) then
317			print *, 'Xvrai plus grand que pi !'
318	guez	91	STOP 1
319			end if
320
321	guez	78	is2 = i
322	guez	119	END IF
323
324	guez	78	idif = iim -is2
325	guez	119
326	guez	78	DO ii = 1, is2
327	guez	119	xlon(ii + idif) = xvrai(ii)
328	guez	91	xprimm(ii + idif) = xxprim(ii)
329	guez	119	END DO
330
331	guez	78	DO ii = 1, idif
332	guez	119	xlon(ii) = xvrai(ii + is2) - twopi_d
333	guez	91	xprimm(ii) = xxprim(ii + is2)
334	guez	119	END DO
335			END IF
336			END IF
337	guez	3
338	guez	91	! Fin de la reorganisation
339	guez	3
340	guez	119	xlon(iim + 1) = xlon(1) + twopi_d
341			xprimm(iim + 1) = xprimm(1)
342	guez	3
343	guez	91	DO i = 1, iim + 1
344	guez	119	xvrai(i) = xlon(i)*180. / pi_d
345			END DO
346	guez	3
347	guez	91	IF (ik == 1) THEN
348			DO i = 1, iim + 1
349	guez	78	rlonm025(i) = xlon(i)
350			xprimm025(i) = xprimm(i)
351	guez	119	END DO
352	guez	91	ELSE IF (ik == 2) THEN
353	guez	119	rlonv = xlon
354			xprimv = xprimm
355	guez	91	ELSE IF (ik == 3) THEN
356	guez	78	DO i = 1, iim + 1
357			rlonu(i) = xlon(i)
358			xprimu(i) = xprimm(i)
359	guez	119	END DO
360	guez	91	ELSE IF (ik == 4) THEN
361	guez	78	DO i = 1, iim + 1
362			rlonp025(i) = xlon(i)
363			xprimp025(i) = xprimm(i)
364	guez	119	END DO
365			END IF
366			end DO loop_ik
367	guez	3
368	guez	91	print *
369	guez	3
370	guez	78	DO i = 1, iim
371	guez	91	xlon(i) = rlonv(i + 1) - rlonv(i)
372	guez	119	END DO
373			champmin = 1e12
374			champmax = -1e12
375	guez	78	DO i = 1, iim
376			champmin = MIN(champmin, xlon(i))
377			champmax = MAX(champmax, xlon(i))
378	guez	119	END DO
379			champmin = champmin * 180. / pi_d
380			champmax = champmax * 180. / pi_d
381	guez	3
382	guez	119	DO i = 1, iim + 1
383			IF (rlonp025(i) < rlonv(i)) THEN
384			print *, ' Attention ! rlonp025 < rlonv', i
385			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	guez	78	END SUBROUTINE fxhyp
408
409			end module fxhyp_m