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.

    ! Nota bene : il vaut mieux avoir grossismy * dzoomy < pi / 2 (radians).

    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.

    IF (dzoomy<1.) THEN
       dzoom = dzoomy*pi
    ELSE IF (dzoomy<12.) THEN
       print *, "Le paramètre dzoomy pour fyhyp est trop petit. L'augmenter " &
            // "et relancer."
       STOP 1
    ELSE
       dzoom = dzoomy * pi/180.
    END IF

    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	119	! Nota bene : il vaut mieux avoir grossismy * dzoomy < pi / 2 (radians).
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	3
65	guez	119	IF (dzoomy<1.) THEN
66			dzoom = dzoomy*pi
67			ELSE IF (dzoomy<12.) THEN
68	guez	97	print *, "Le paramètre dzoomy pour fyhyp est trop petit. L'augmenter " &
69			// "et relancer."
70			STOP 1
71	guez	81	ELSE
72	guez	119	dzoom = dzoomy * pi/180.
73	guez	81	END IF
74
75	guez	119	print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
76			print *, y0, grossismy, tauy, dzoom
77	guez	81
78	guez	97	DO i = 0, nmax2
79			yt(i) = -pis2 + real(i)*pi/nmax2
80			END DO
81	guez	81
82	guez	97	heavyy0m = heavyside(-y0)
83			heavyy0 = heavyside(y0)
84			y0min = 2.y0heavyy0m - pis2
85			y0max = 2.y0heavyy0 + pis2
86	guez	81
87	guez	97	fa = 999.999
88			fb = 999.999
89	guez	81
90	guez	97	DO i = 0, nmax2
91			IF (yt(i)<y0) THEN
92	guez	119	fa(i) = tauy*(yt(i)-y0 + dzoom/2.)
93	guez	97	fb(i) = (yt(i)-2.y0heavyy0m + pis2)*(y0-yt(i))
94			ELSE IF (yt(i)>y0) THEN
95	guez	119	fa(i) = tauy*(y0-yt(i) + dzoom/2.)
96	guez	97	fb(i) = (2.y0heavyy0-yt(i) + pis2)*(yt(i)-y0)
97			END IF
98	guez	81
99	guez	97	IF (200.*fb(i)<-fa(i)) THEN
100			fhyp(i) = -1.
101			ELSE IF (200.*fb(i)<fa(i)) THEN
102			fhyp(i) = 1.
103			ELSE
104			fhyp(i) = tanh(fa(i)/fb(i))
105			END IF
106	guez	81
107	guez	97	IF (yt(i)==y0) fhyp(i) = 1.
108			IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
109			END DO
110	guez	81
111	guez	97	! Calcul de beta
112	guez	81
113	guez	97	ffdy = 0.
114	guez	81
115	guez	97	DO i = 1, nmax2
116			ymoy = 0.5*(yt(i-1) + yt(i))
117			IF (ymoy<y0) THEN
118	guez	119	fa(i) = tauy*(ymoy-y0 + dzoom/2.)
119	guez	97	fb(i) = (ymoy-2.y0heavyy0m + pis2)*(y0-ymoy)
120			ELSE IF (ymoy>y0) THEN
121	guez	119	fa(i) = tauy*(y0-ymoy + dzoom/2.)
122	guez	97	fb(i) = (2.y0heavyy0-ymoy + pis2)*(ymoy-y0)
123			END IF
124	guez	81
125	guez	97	IF (200.*fb(i)<-fa(i)) THEN
126			fxm(i) = -1.
127			ELSE IF (200.*fb(i)<fa(i)) THEN
128			fxm(i) = 1.
129			ELSE
130			fxm(i) = tanh(fa(i)/fb(i))
131			END IF
132			IF (ymoy==y0) fxm(i) = 1.
133			IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
134			ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
135			END DO
136	guez	81
137	guez	119	beta = (grossismy*ffdy-pi)/(ffdy-pi)
138	guez	81
139	guez	119	IF (2. * beta - grossismy <= 0.) THEN
140	guez	97	print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
141			// 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
142			// 'dzoomy et relancer.'
143			STOP 1
144	guez	81	END IF
145
146	guez	97	! calcul de Ytprim
147	guez	81
148	guez	97	DO i = 0, nmax2
149	guez	119	ytprim(i) = beta + (grossismy-beta)*fhyp(i)
150	guez	97	END DO
151	guez	81
152	guez	97	! Calcul de Yf
153	guez	81
154	guez	97	yf(0) = -pis2
155			DO i = 1, nmax2
156	guez	119	yypr(i) = beta + (grossismy-beta)*fxm(i)
157	guez	81	END DO
158	guez	3
159	guez	97	DO i = 1, nmax2
160			yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
161	guez	81	END DO
162	guez	3
163	guez	97	! yuv = 0. si calcul des latitudes aux pts. U
164			! yuv = 0.5 si calcul des latitudes aux pts. V
165	guez	3
166	guez	97	loop_ik: DO ik = 1, 4
167			IF (ik==1) THEN
168			yuv = 0.
169			jlat = jjm + 1
170			ELSE IF (ik==2) THEN
171			yuv = 0.5
172			jlat = jjm
173			ELSE IF (ik==3) THEN
174			yuv = 0.25
175			jlat = jjm
176			ELSE IF (ik==4) THEN
177			yuv = 0.75
178			jlat = jjm
179			END IF
180	guez	3
181	guez	97	yo1 = 0.
182			DO j = 1, jlat
183			yo1 = 0.
184			ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
185			yfi = ylon2
186	guez	3
187	guez	97	it = nmax2
188			DO while (it >= 1 .and. yfi < yf(it))
189			it = it - 1
190			END DO
191	guez	3
192	guez	97	yi = yt(it)
193			IF (it==nmax2) THEN
194			it = nmax2 - 1
195			yf(it + 1) = pis2
196			END IF
197	guez	3
198	guez	97	! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
199			! et Y'(yi)
200	guez	3
201	guez	97	CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
202			yt(it), yt(it + 1), a0, a1, a2, a3)
203	guez	3
204	guez	97	yf1 = yf(it)
205			yprimin = a1 + 2.a2yi + 3.a3yi*yi
206	guez	3
207	guez	97	iter = 1
208			DO
209			yi = yi - (yf1-yfi)/yprimin
210			IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
211			yo1 = yi
212			yi2 = yi*yi
213			yf1 = a0 + a1yi + a2yi2 + a3yi2yi
214			yprimin = a1 + 2.a2yi + 3.a3yi2
215			END DO
216			if (abs(yi-yo1) > epsilon) then
217			print *, 'Pas de solution.', j, ylon2
218			STOP 1
219			end if
220	guez	3
221	guez	97	yprimin = a1 + 2.a2yi + 3.a3yi*yi
222			yprim(j) = pi/(jjm*yprimin)
223			yvrai(j) = yi
224			END DO
225	guez	3
226	guez	97	DO j = 1, jlat - 1
227			IF (yvrai(j + 1)<yvrai(j)) THEN
228			print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
229			j, ')'
230			STOP 1
231			END IF
232			END DO
233	guez	3
234	guez	97	print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'
235	guez	3
236	guez	97	IF (ik==1) THEN
237			ypn = pis2
238	guez	118	DO j = jjm + 1, 1, -1
239	guez	97	IF (yvrai(j)<=ypn) exit
240			END DO
241	guez	3
242	guez	97	jpn = j
243			y00 = yvrai(jpn)
244			deply = pis2 - y00
245			END IF
246	guez	3
247	guez	97	DO j = 1, jjm + 1 - jpn
248			ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
249			yprimm(j) = yprim(jpn + j-1)
250			END DO
251	guez	3
252	guez	97	jjpn = jpn
253			IF (jlat==jjm) jjpn = jpn - 1
254	guez	3
255	guez	97	DO j = 1, jjpn
256			ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
257			yprimm(j + jjm + 1-jpn) = yprim(j)
258			END DO
259	guez	3
260	guez	97	! Fin de la reorganisation
261	guez	3
262	guez	97	DO j = 1, jlat
263			ylat(j) = ylatt(jlat + 1-j)
264			yprim(j) = yprimm(jlat + 1-j)
265			END DO
266	guez	81
267	guez	97	DO j = 1, jlat
268			yvrai(j) = ylat(j)*180./pi
269			END DO
270	guez	81
271	guez	97	IF (ik==1) THEN
272	guez	118	DO j = 1, jjm + 1
273	guez	119	rlatu(j) = ylat(j)
274	guez	97	yyprimu(j) = yprim(j)
275			END DO
276			ELSE IF (ik==2) THEN
277	guez	118	DO j = 1, jjm
278	guez	119	rlatv(j) = ylat(j)
279	guez	97	END DO
280			ELSE IF (ik==3) THEN
281	guez	118	DO j = 1, jjm
282	guez	97	rlatu2(j) = ylat(j)
283			yprimu2(j) = yprim(j)
284			END DO
285			ELSE IF (ik==4) THEN
286	guez	118	DO j = 1, jjm
287	guez	97	rlatu1(j) = ylat(j)
288			yprimu1(j) = yprim(j)
289			END DO
290			END IF
291			END DO loop_ik
292	guez	81
293	guez	97	DO j = 1, jjm
294	guez	119	ylat(j) = rlatu(j) - rlatu(j + 1)
295	guez	97	END DO
296			champmin = 1e12
297			champmax = -1e12
298			DO j = 1, jjm
299			champmin = min(champmin, ylat(j))
300			champmax = max(champmax, ylat(j))
301			END DO
302			champmin = champmin*180./pi
303			champmax = champmax*180./pi
304	guez	81
305	guez	119	DO j = 1, jjm
306			IF (rlatu1(j) <= rlatu2(j)) THEN
307			print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
308			STOP 13
309			ENDIF
310
311			IF (rlatu2(j) <= rlatu(j+1)) THEN
312			print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
313			STOP 14
314			ENDIF
315
316			IF (rlatu(j) <= rlatu1(j)) THEN
317			print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
318			STOP 15
319			ENDIF
320
321			IF (rlatv(j) <= rlatu2(j)) THEN
322			print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
323			STOP 16
324			ENDIF
325
326			IF (rlatv(j) >= rlatu1(j)) THEN
327			print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
328			STOP 17
329			ENDIF
330
331			IF (rlatv(j) >= rlatu(j)) THEN
332			print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
333			STOP 18
334			ENDIF
335			ENDDO
336
337			print *, 'Latitudes'
338			print 3, champmin, champmax
339
340			3 Format(1x, ' Au centre du zoom, la longueur de la maille est', &
341			' d environ ', f0.2, ' degres ', /, &
342			' alors que la maille en dehors de la zone du zoom est ', &
343	guez	120	"d'environ ", f0.2, ' degres ')
344	guez	119
345	guez	97	END SUBROUTINE fyhyp
346	guez	81
347	guez	97	end module fyhyp_m