trunk/dyn3d/fyhyp.f

module fyhyp_m

 IMPLICIT NONE

contains

  SUBROUTINE fyhyp(yzoomdeg, grossism, dzooma, tau, rrlatu, yyprimu, rrlatv, &
       yyprimv, rlatu2, yprimu2, rlatu1, yprimu1, champmin, champmax)

    ! 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) à tangente hyperbolique.

    ! Nota bene : il vaut mieux avoir grossism * dzoom < pi / 2 (rad),
    ! en latitude.

    USE dimens_m, only: jjm
    USE paramet_m, only: JJP1

    REAL, intent(in):: yzoomdeg

    REAL, intent(in):: grossism
    ! grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)

    REAL, intent(in):: dzooma

    REAL, intent(in):: tau
    ! raideur de la transition de l'intérieur à l'extérieur du zoom

    ! arguments de sortie 

    REAL rrlatu(jjp1), yyprimu(jjp1), rrlatv(jjm), yyprimv(jjm)
    real rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)
    DOUBLE PRECISION champmin, champmax

    ! Local: 

    INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
    REAL dzoom ! distance totale de la zone du zoom (en radians)
    DOUBLE PRECISION ylat(jjp1), yprim(jjp1)
    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(jjp1), yprimm(jjp1), ylatt(jjp1)
    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 = yzoomdeg*pi/180.

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

    print *, 'yzoom(rad), grossism, tau, dzoom (rad):'
    print *, y0, grossism, tau, 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) = tau*(yt(i)-y0 + dzoom/2.)
          fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))
       ELSE IF (yt(i)>y0) THEN
          fa(i) = tau*(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) = tau*(ymoy-y0 + dzoom/2.)
          fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)
       ELSE IF (ymoy>y0) THEN
          fa(i) = tau*(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 = (grossism*ffdy-pi)/(ffdy-pi)

    IF (2. * beta - grossism <= 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 + (grossism-beta)*fhyp(i)
    END DO

    ! Calcul de Yf 

    yf(0) = -pis2
    DO i = 1, nmax2
       yypr(i) = beta + (grossism-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 = jlat, 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, jlat
             rrlatu(j) = ylat(j)
             yyprimu(j) = yprim(j)
          END DO
       ELSE IF (ik==2) THEN
          DO j = 1, jlat
             rrlatv(j) = ylat(j)
             yyprimv(j) = yprim(j)
          END DO
       ELSE IF (ik==3) THEN
          DO j = 1, jlat
             rlatu2(j) = ylat(j)
             yprimu2(j) = yprim(j)
          END DO
       ELSE IF (ik==4) THEN
          DO j = 1, jlat
             rlatu1(j) = ylat(j)
             yprimu1(j) = yprim(j)
          END DO
       END IF
    END DO loop_ik

    DO j = 1, jjm
       ylat(j) = rrlatu(j) - rrlatu(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

  END SUBROUTINE fyhyp

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