trunk/dyn3d/fyhyp.f

module fyhyp_m

  IMPLICIT NONE

contains

  SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)

    ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32

    ! Author: P. Le Van, from analysis by R. Sadourny 

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

    ! Il vaut mieux avoir : grossismy * dzoom < pi / 2

    use coefpoly_m, only: coefpoly
    USE dimens_m, only: jjm
    use serre, only: clat, grossismy, dzoomy, tauy

    REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1)
    REAL, intent(out):: rlatv(jjm)
    real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)

    ! Local: 

    DOUBLE PRECISION champmin, champmax
    INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
    REAL dzoom ! distance totale de la zone du zoom (en radians)
    DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
    DOUBLE PRECISION yuv
    DOUBLE PRECISION, save:: yt(0:nmax2)
    DOUBLE PRECISION fhyp(0:nmax2), beta
    DOUBLE PRECISION, save:: ytprim(0:nmax2)
    DOUBLE PRECISION fxm(0:nmax2)
    DOUBLE PRECISION, save:: yf(0:nmax2)
    DOUBLE PRECISION yypr(0:nmax2)
    DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
    DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm
    DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
    DOUBLE PRECISION yfi, yf1, ffdy
    DOUBLE PRECISION ypn, deply, y00
    SAVE y00, deply

    INTEGER i, j, it, ik, iter, jlat
    INTEGER jpn, jjpn
    SAVE jpn
    DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m
    DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
    REAL y0min, y0max

    DOUBLE PRECISION heavyside

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

    print *, "Call sequence information: fyhyp"

    pi = 2.*asin(1.)
    pis2 = pi/2.
    pisjm = pi/real(jjm)
    epsilon = 1e-3
    y0 = clat*pi/180.
    dzoom = dzoomy*pi
    print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
    print *, y0, grossismy, tauy, dzoom

    DO i = 0, nmax2
       yt(i) = -pis2 + real(i)*pi/nmax2
    END DO

    heavyy0m = heavyside(-y0)
    heavyy0 = heavyside(y0)
    y0min = 2.*y0*heavyy0m - pis2
    y0max = 2.*y0*heavyy0 + pis2

    fa = 999.999
    fb = 999.999

    DO i = 0, nmax2
       IF (yt(i)<y0) THEN
          fa(i) = tauy*(yt(i)-y0 + dzoom/2.)
          fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))
       ELSE IF (yt(i)>y0) THEN
          fa(i) = tauy*(y0-yt(i) + dzoom/2.)
          fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0)
       END IF

       IF (200.*fb(i)<-fa(i)) THEN
          fhyp(i) = -1.
       ELSE IF (200.*fb(i)<fa(i)) THEN
          fhyp(i) = 1.
       ELSE
          fhyp(i) = tanh(fa(i)/fb(i))
       END IF

       IF (yt(i)==y0) fhyp(i) = 1.
       IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
    END DO

    ! Calcul de beta 

    ffdy = 0.

    DO i = 1, nmax2
       ymoy = 0.5*(yt(i-1) + yt(i))
       IF (ymoy<y0) THEN
          fa(i) = tauy*(ymoy-y0 + dzoom/2.)
          fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)
       ELSE IF (ymoy>y0) THEN
          fa(i) = tauy*(y0-ymoy + dzoom/2.)
          fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0)
       END IF

       IF (200.*fb(i)<-fa(i)) THEN
          fxm(i) = -1.
       ELSE IF (200.*fb(i)<fa(i)) THEN
          fxm(i) = 1.
       ELSE
          fxm(i) = tanh(fa(i)/fb(i))
       END IF
       IF (ymoy==y0) fxm(i) = 1.
       IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
       ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
    END DO

    beta = (grossismy*ffdy-pi)/(ffdy-pi)

    IF (2. * beta - grossismy <= 0.) THEN
       print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
            // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
            // 'dzoomy et relancer.'
       STOP 1
    END IF

    ! calcul de Ytprim 

    DO i = 0, nmax2
       ytprim(i) = beta + (grossismy-beta)*fhyp(i)
    END DO

    ! Calcul de Yf 

    yf(0) = -pis2
    DO i = 1, nmax2
       yypr(i) = beta + (grossismy-beta)*fxm(i)
    END DO

    DO i = 1, nmax2
       yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
    END DO

    ! yuv = 0. si calcul des latitudes aux pts. U 
    ! yuv = 0.5 si calcul des latitudes aux pts. V 

    loop_ik: DO ik = 1, 4
       IF (ik==1) THEN
          yuv = 0.
          jlat = jjm + 1
       ELSE IF (ik==2) THEN
          yuv = 0.5
          jlat = jjm
       ELSE IF (ik==3) THEN
          yuv = 0.25
          jlat = jjm
       ELSE IF (ik==4) THEN
          yuv = 0.75
          jlat = jjm
       END IF

       yo1 = 0.
       DO j = 1, jlat
          yo1 = 0.
          ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
          yfi = ylon2

          it = nmax2
          DO while (it >= 1 .and. yfi < yf(it))
             it = it - 1
          END DO

          yi = yt(it)
          IF (it==nmax2) THEN
             it = nmax2 - 1
             yf(it + 1) = pis2
          END IF

          ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
          ! et Y'(yi) 

          CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
               yt(it), yt(it + 1), a0, a1, a2, a3)

          yf1 = yf(it)
          yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi

          iter = 1
          DO
             yi = yi - (yf1-yfi)/yprimin
             IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
             yo1 = yi
             yi2 = yi*yi
             yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi
             yprimin = a1 + 2.*a2*yi + 3.*a3*yi2
          END DO
          if (abs(yi-yo1) > epsilon) then
             print *, 'Pas de solution.', j, ylon2
             STOP 1
          end if

          yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
          yprim(j) = pi/(jjm*yprimin)
          yvrai(j) = yi
       END DO

       DO j = 1, jlat - 1
          IF (yvrai(j + 1)<yvrai(j)) THEN
             print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
                  j, ')'
             STOP 1
          END IF
       END DO

       print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'

       IF (ik==1) THEN
          ypn = pis2
          DO j = jjm + 1, 1, -1
             IF (yvrai(j)<=ypn) exit
          END DO

          jpn = j
          y00 = yvrai(jpn)
          deply = pis2 - y00
       END IF

       DO j = 1, jjm + 1 - jpn
          ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
          yprimm(j) = yprim(jpn + j-1)
       END DO

       jjpn = jpn
       IF (jlat==jjm) jjpn = jpn - 1

       DO j = 1, jjpn
          ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
          yprimm(j + jjm + 1-jpn) = yprim(j)
       END DO

       ! Fin de la reorganisation

       DO j = 1, jlat
          ylat(j) = ylatt(jlat + 1-j)
          yprim(j) = yprimm(jlat + 1-j)
       END DO

       DO j = 1, jlat
          yvrai(j) = ylat(j)*180./pi
       END DO

       IF (ik==1) THEN
          DO j = 1, jjm + 1
             rlatu(j) = ylat(j)
             yyprimu(j) = yprim(j)
          END DO
       ELSE IF (ik==2) THEN
          DO j = 1, jjm
             rlatv(j) = ylat(j)
          END DO
       ELSE IF (ik==3) THEN
          DO j = 1, jjm
             rlatu2(j) = ylat(j)
             yprimu2(j) = yprim(j)
          END DO
       ELSE IF (ik==4) THEN
          DO j = 1, jjm
             rlatu1(j) = ylat(j)
             yprimu1(j) = yprim(j)
          END DO
       END IF
    END DO loop_ik

    DO j = 1, jjm
       ylat(j) = rlatu(j) - rlatu(j + 1)
    END DO
    champmin = 1e12
    champmax = -1e12
    DO j = 1, jjm
       champmin = min(champmin, ylat(j))
       champmax = max(champmax, ylat(j))
    END DO
    champmin = champmin*180./pi
    champmax = champmax*180./pi

    DO j = 1, jjm
       IF (rlatu1(j) <= rlatu2(j)) THEN
          print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
          STOP 13
       ENDIF

       IF (rlatu2(j) <= rlatu(j+1)) THEN
          print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
          STOP 14
       ENDIF

       IF (rlatu(j) <= rlatu1(j)) THEN
          print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
          STOP 15
       ENDIF

       IF (rlatv(j) <= rlatu2(j)) THEN
          print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
          STOP 16
       ENDIF

       IF (rlatv(j) >= rlatu1(j)) THEN
          print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
          STOP 17
       ENDIF

       IF (rlatv(j) >= rlatu(j)) THEN
          print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
          STOP 18
       ENDIF
    ENDDO

    print *, 'Latitudes'
    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 fyhyp

end module fyhyp_m
1	guez	97	module fyhyp_m
2	guez	3
3	guez	112	IMPLICIT NONE
4	guez	3
5	guez	97	contains
6	guez	3
7	guez	119	SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
8	guez	3
9	guez	97	! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
10	guez	3
11	guez	97	! Author: P. Le Van, from analysis by R. Sadourny
12	guez	3
13	guez	97	! Calcule les latitudes et dérivées dans la grille du GCM pour une
14	guez	119	! fonction f(y) à dérivée tangente hyperbolique.
15	guez	3
16	guez	121	! Il vaut mieux avoir : grossismy * dzoom < pi / 2
17	guez	3
18	guez	120	use coefpoly_m, only: coefpoly
19	guez	97	USE dimens_m, only: jjm
20	guez	119	use serre, only: clat, grossismy, dzoomy, tauy
21	guez	3
22	guez	119	REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1)
23			REAL, intent(out):: rlatv(jjm)
24	guez	118	real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)
25	guez	3
26	guez	97	! Local:
27	guez	3
28	guez	119	DOUBLE PRECISION champmin, champmax
29	guez	97	INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
30			REAL dzoom ! distance totale de la zone du zoom (en radians)
31	guez	118	DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
32	guez	97	DOUBLE PRECISION yuv
33			DOUBLE PRECISION, save:: yt(0:nmax2)
34			DOUBLE PRECISION fhyp(0:nmax2), beta
35			DOUBLE PRECISION, save:: ytprim(0:nmax2)
36			DOUBLE PRECISION fxm(0:nmax2)
37			DOUBLE PRECISION, save:: yf(0:nmax2)
38			DOUBLE PRECISION yypr(0:nmax2)
39	guez	118	DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
40	guez	97	DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm
41			DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
42			DOUBLE PRECISION yfi, yf1, ffdy
43			DOUBLE PRECISION ypn, deply, y00
44			SAVE y00, deply
45	guez	3
46	guez	97	INTEGER i, j, it, ik, iter, jlat
47			INTEGER jpn, jjpn
48			SAVE jpn
49			DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m
50			DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
51			REAL y0min, y0max
52	guez	3
53	guez	97	DOUBLE PRECISION heavyside
54	guez	3
55	guez	97	!-------------------------------------------------------------------
56	guez	3
57	guez	120	print *, "Call sequence information: fyhyp"
58
59	guez	97	pi = 2.*asin(1.)
60			pis2 = pi/2.
61			pisjm = pi/real(jjm)
62			epsilon = 1e-3
63	guez	119	y0 = clat*pi/180.
64	guez	121	dzoom = dzoomy*pi
65	guez	119	print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
66			print *, y0, grossismy, tauy, dzoom
67	guez	81
68	guez	97	DO i = 0, nmax2
69			yt(i) = -pis2 + real(i)*pi/nmax2
70			END DO
71	guez	81
72	guez	97	heavyy0m = heavyside(-y0)
73			heavyy0 = heavyside(y0)
74			y0min = 2.y0heavyy0m - pis2
75			y0max = 2.y0heavyy0 + pis2
76	guez	81
77	guez	97	fa = 999.999
78			fb = 999.999
79	guez	81
80	guez	97	DO i = 0, nmax2
81			IF (yt(i)<y0) THEN
82	guez	119	fa(i) = tauy*(yt(i)-y0 + dzoom/2.)
83	guez	97	fb(i) = (yt(i)-2.y0heavyy0m + pis2)*(y0-yt(i))
84			ELSE IF (yt(i)>y0) THEN
85	guez	119	fa(i) = tauy*(y0-yt(i) + dzoom/2.)
86	guez	97	fb(i) = (2.y0heavyy0-yt(i) + pis2)*(yt(i)-y0)
87			END IF
88	guez	81
89	guez	97	IF (200.*fb(i)<-fa(i)) THEN
90			fhyp(i) = -1.
91			ELSE IF (200.*fb(i)<fa(i)) THEN
92			fhyp(i) = 1.
93			ELSE
94			fhyp(i) = tanh(fa(i)/fb(i))
95			END IF
96	guez	81
97	guez	97	IF (yt(i)==y0) fhyp(i) = 1.
98			IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
99			END DO
100	guez	81
101	guez	97	! Calcul de beta
102	guez	81
103	guez	97	ffdy = 0.
104	guez	81
105	guez	97	DO i = 1, nmax2
106			ymoy = 0.5*(yt(i-1) + yt(i))
107			IF (ymoy<y0) THEN
108	guez	119	fa(i) = tauy*(ymoy-y0 + dzoom/2.)
109	guez	97	fb(i) = (ymoy-2.y0heavyy0m + pis2)*(y0-ymoy)
110			ELSE IF (ymoy>y0) THEN
111	guez	119	fa(i) = tauy*(y0-ymoy + dzoom/2.)
112	guez	97	fb(i) = (2.y0heavyy0-ymoy + pis2)*(ymoy-y0)
113			END IF
114	guez	81
115	guez	97	IF (200.*fb(i)<-fa(i)) THEN
116			fxm(i) = -1.
117			ELSE IF (200.*fb(i)<fa(i)) THEN
118			fxm(i) = 1.
119			ELSE
120			fxm(i) = tanh(fa(i)/fb(i))
121			END IF
122			IF (ymoy==y0) fxm(i) = 1.
123			IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
124			ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
125			END DO
126	guez	81
127	guez	119	beta = (grossismy*ffdy-pi)/(ffdy-pi)
128	guez	81
129	guez	119	IF (2. * beta - grossismy <= 0.) THEN
130	guez	97	print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
131			// 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
132			// 'dzoomy et relancer.'
133			STOP 1
134	guez	81	END IF
135
136	guez	97	! calcul de Ytprim
137	guez	81
138	guez	97	DO i = 0, nmax2
139	guez	119	ytprim(i) = beta + (grossismy-beta)*fhyp(i)
140	guez	97	END DO
141	guez	81
142	guez	97	! Calcul de Yf
143	guez	81
144	guez	97	yf(0) = -pis2
145			DO i = 1, nmax2
146	guez	119	yypr(i) = beta + (grossismy-beta)*fxm(i)
147	guez	81	END DO
148	guez	3
149	guez	97	DO i = 1, nmax2
150			yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
151	guez	81	END DO
152	guez	3
153	guez	97	! yuv = 0. si calcul des latitudes aux pts. U
154			! yuv = 0.5 si calcul des latitudes aux pts. V
155	guez	3
156	guez	97	loop_ik: DO ik = 1, 4
157			IF (ik==1) THEN
158			yuv = 0.
159			jlat = jjm + 1
160			ELSE IF (ik==2) THEN
161			yuv = 0.5
162			jlat = jjm
163			ELSE IF (ik==3) THEN
164			yuv = 0.25
165			jlat = jjm
166			ELSE IF (ik==4) THEN
167			yuv = 0.75
168			jlat = jjm
169			END IF
170	guez	3
171	guez	97	yo1 = 0.
172			DO j = 1, jlat
173			yo1 = 0.
174			ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
175			yfi = ylon2
176	guez	3
177	guez	97	it = nmax2
178			DO while (it >= 1 .and. yfi < yf(it))
179			it = it - 1
180			END DO
181	guez	3
182	guez	97	yi = yt(it)
183			IF (it==nmax2) THEN
184			it = nmax2 - 1
185			yf(it + 1) = pis2
186			END IF
187	guez	3
188	guez	97	! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
189			! et Y'(yi)
190	guez	3
191	guez	97	CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
192			yt(it), yt(it + 1), a0, a1, a2, a3)
193	guez	3
194	guez	97	yf1 = yf(it)
195			yprimin = a1 + 2.a2yi + 3.a3yi*yi
196	guez	3
197	guez	97	iter = 1
198			DO
199			yi = yi - (yf1-yfi)/yprimin
200			IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
201			yo1 = yi
202			yi2 = yi*yi
203			yf1 = a0 + a1yi + a2yi2 + a3yi2yi
204			yprimin = a1 + 2.a2yi + 3.a3yi2
205			END DO
206			if (abs(yi-yo1) > epsilon) then
207			print *, 'Pas de solution.', j, ylon2
208			STOP 1
209			end if
210	guez	3
211	guez	97	yprimin = a1 + 2.a2yi + 3.a3yi*yi
212			yprim(j) = pi/(jjm*yprimin)
213			yvrai(j) = yi
214			END DO
215	guez	3
216	guez	97	DO j = 1, jlat - 1
217			IF (yvrai(j + 1)<yvrai(j)) THEN
218			print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
219			j, ')'
220			STOP 1
221			END IF
222			END DO
223	guez	3
224	guez	97	print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'
225	guez	3
226	guez	97	IF (ik==1) THEN
227			ypn = pis2
228	guez	118	DO j = jjm + 1, 1, -1
229	guez	97	IF (yvrai(j)<=ypn) exit
230			END DO
231	guez	3
232	guez	97	jpn = j
233			y00 = yvrai(jpn)
234			deply = pis2 - y00
235			END IF
236	guez	3
237	guez	97	DO j = 1, jjm + 1 - jpn
238			ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
239			yprimm(j) = yprim(jpn + j-1)
240			END DO
241	guez	3
242	guez	97	jjpn = jpn
243			IF (jlat==jjm) jjpn = jpn - 1
244	guez	3
245	guez	97	DO j = 1, jjpn
246			ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
247			yprimm(j + jjm + 1-jpn) = yprim(j)
248			END DO
249	guez	3
250	guez	97	! Fin de la reorganisation
251	guez	3
252	guez	97	DO j = 1, jlat
253			ylat(j) = ylatt(jlat + 1-j)
254			yprim(j) = yprimm(jlat + 1-j)
255			END DO
256	guez	81
257	guez	97	DO j = 1, jlat
258			yvrai(j) = ylat(j)*180./pi
259			END DO
260	guez	81
261	guez	97	IF (ik==1) THEN
262	guez	118	DO j = 1, jjm + 1
263	guez	119	rlatu(j) = ylat(j)
264	guez	97	yyprimu(j) = yprim(j)
265			END DO
266			ELSE IF (ik==2) THEN
267	guez	118	DO j = 1, jjm
268	guez	119	rlatv(j) = ylat(j)
269	guez	97	END DO
270			ELSE IF (ik==3) THEN
271	guez	118	DO j = 1, jjm
272	guez	97	rlatu2(j) = ylat(j)
273			yprimu2(j) = yprim(j)
274			END DO
275			ELSE IF (ik==4) THEN
276	guez	118	DO j = 1, jjm
277	guez	97	rlatu1(j) = ylat(j)
278			yprimu1(j) = yprim(j)
279			END DO
280			END IF
281			END DO loop_ik
282	guez	81
283	guez	97	DO j = 1, jjm
284	guez	119	ylat(j) = rlatu(j) - rlatu(j + 1)
285	guez	97	END DO
286			champmin = 1e12
287			champmax = -1e12
288			DO j = 1, jjm
289			champmin = min(champmin, ylat(j))
290			champmax = max(champmax, ylat(j))
291			END DO
292			champmin = champmin*180./pi
293			champmax = champmax*180./pi
294	guez	81
295	guez	119	DO j = 1, jjm
296			IF (rlatu1(j) <= rlatu2(j)) THEN
297			print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
298			STOP 13
299			ENDIF
300
301			IF (rlatu2(j) <= rlatu(j+1)) THEN
302			print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
303			STOP 14
304			ENDIF
305
306			IF (rlatu(j) <= rlatu1(j)) THEN
307			print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
308			STOP 15
309			ENDIF
310
311			IF (rlatv(j) <= rlatu2(j)) THEN
312			print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
313			STOP 16
314			ENDIF
315
316			IF (rlatv(j) >= rlatu1(j)) THEN
317			print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
318			STOP 17
319			ENDIF
320
321			IF (rlatv(j) >= rlatu(j)) THEN
322			print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
323			STOP 18
324			ENDIF
325			ENDDO
326
327			print *, 'Latitudes'
328			print 3, champmin, champmax
329
330			3 Format(1x, ' Au centre du zoom, la longueur de la maille est', &
331			' d environ ', f0.2, ' degres ', /, &
332			' alors que la maille en dehors de la zone du zoom est ', &
333	guez	120	"d'environ ", f0.2, ' degres ')
334	guez	119
335	guez	97	END SUBROUTINE fyhyp
336	guez	81
337	guez	97	end module fyhyp_m