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

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

    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
    print *, 'Si cette derniere est trop lache, modifiez les parametres'
    print *, 'grossismy, tauy, dzoom pour Y et repasser ! '

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