Context Navigation

fyhyp.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 10.1 KB

Line
1	!
2	! $Id: fyhyp.F 1907 2013-11-26 13:10:46Z lguez $
3	!
4	c
5	c
6	SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau ,
7	, rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
8	, champmin,champmax )
9
10	cc ... Version du 01/04/2001 ....
11
12	IMPLICIT NONE
13	c
14	c ... Auteur : P. Le Van ...
15	c
16	c ....... d'apres formulations de R. Sadourny .......
17	c
18	c Calcule les latitudes et derivees dans la grille du GCM pour une
19	c fonction f(y) a tangente hyperbolique .
20	c
21	c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
22	c dzoom etant la distance totale de la zone du zoom ( en radians )
23	c tau la raideur de la transition de l'interieur a l'exterieur du zoom
24	c
25	c
26	c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati.
27	c ********************************************************************
28	c
29	c
30	#include "dimensions.h"
31	#include "paramet.h"
32
33	INTEGER nmax , nmax2
34	PARAMETER ( nmax = 30000, nmax2 = 2*nmax )
35	c
36	c
37	c ....... arguments d'entree .......
38	c
39	REAL yzoomdeg, grossism,dzooma,tau
40	c ( rentres par run.def )
41
42	c ....... arguments de sortie .......
43	c
44	REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm),
45	, rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
46
47	c
48	c ..... champs locaux .....
49	c
50
51	REAL dzoom
52	REAL(KIND=8) ylat(jjp1), yprim(jjp1)
53	REAL(KIND=8) yuv
54	REAL(KIND=8) yt(0:nmax2)
55	REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
56	SAVE Ytprim, yt,Yf
57	REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2)
58	REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
59	REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm
60	REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
61	REAL(KIND=8) yfi,Yf1,ffdy
62	REAL(KIND=8) ypn,deply,y00
63	SAVE y00, deply
64
65	INTEGER i,j,it,ik,iter,jlat
66	INTEGER jpn,jjpn
67	SAVE jpn
68	REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m
69	REAL(KIND=8) fa(0:nmax2),fb(0:nmax2)
70	REAL y0min,y0max
71
72	REAL(KIND=8) heavyside
73
74	pi = 2. * ASIN(1.)
75	depi = 2. * pi
76	pis2 = pi/2.
77	pisjm = pi/ REAL(jjm)
78	epsilon = 1.e-3
79	y0 = yzoomdeg * pi/180.
80
81	IF( dzooma.LT.1.) THEN
82	dzoom = dzooma * pi
83	ELSEIF( dzooma.LT. 12. ) THEN
84	WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
85	,menter et relancer ! '
86	STOP 1
87	ELSE
88	dzoom = dzooma * pi/180.
89	ENDIF
90
91	WRITE(6,18)
92	WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
93	WRITE(6,24) y0,grossism,tau,dzoom
94
95	DO i = 0, nmax2
96	yt(i) = - pis2 + REAL(i)* pi /nmax2
97	ENDDO
98
99	heavyy0m = heavyside( -y0 )
100	heavyy0 = heavyside( y0 )
101	y0min = 2.y0heavyy0m - pis2
102	y0max = 2.y0heavyy0 + pis2
103
104	fa = 999.999
105	fb = 999.999
106
107	DO i = 0, nmax2
108	IF( yt(i).LT.y0 ) THEN
109	fa (i) = tau* (yt(i)-y0+dzoom/2. )
110	fb(i) = (yt(i)-2.y0heavyy0m +pis2) * ( y0 - yt(i) )
111	ELSEIF ( yt(i).GT.y0 ) THEN
112	fa(i) = tau *(y0-yt(i)+dzoom/2. )
113	fb(i) = (2.y0heavyy0 -yt(i)+pis2) * ( yt(i) - y0 )
114	ENDIF
115
116	IF( 200.* fb(i) .LT. - fa(i) ) THEN
117	fhyp ( i) = - 1.
118	ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
119	fhyp ( i) = 1.
120	ELSE
121	fhyp(i) = TANH ( fa(i)/fb(i) )
122	ENDIF
123
124	IF( yt(i).EQ.y0 ) fhyp(i) = 1.
125	IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
126
127	ENDDO
128
129	cc .... Calcul de beta ....
130	c
131	ffdy = 0.
132
133	DO i = 1, nmax2
134	ymoy = 0.5 * ( yt(i-1) + yt( i ) )
135	IF( ymoy.LT.y0 ) THEN
136	fa(i)= tau * ( ymoy-y0+dzoom/2.)
137	fb(i) = (ymoy-2.y0heavyy0m +pis2) * ( y0 - ymoy )
138	ELSEIF ( ymoy.GT.y0 ) THEN
139	fa(i)= tau * ( y0-ymoy+dzoom/2. )
140	fb(i) = (2.y0heavyy0 -ymoy+pis2) * ( ymoy - y0 )
141	ENDIF
142
143	IF( 200.* fb(i) .LT. - fa(i) ) THEN
144	fxm ( i) = - 1.
145	ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
146	fxm ( i) = 1.
147	ELSE
148	fxm(i) = TANH ( fa(i)/fb(i) )
149	ENDIF
150	IF( ymoy.EQ.y0 ) fxm(i) = 1.
151	IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
152	ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
153
154	ENDDO
155
156	beta = ( grossism * ffdy - pi ) / ( ffdy - pi )
157
158	IF( 2.*beta - grossism.LE. 0.) THEN
159
160	WRITE(6,) ' * Attention ! La valeur beta calculee dans la rou
161	,tine fyhyp est mauvaise ! '
162	WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy',
163	, ' et relancer ! *** '
164	CALL ABORT
165
166	ENDIF
167	c
168	c ..... calcul de Ytprim .....
169	c
170
171	DO i = 0, nmax2
172	Ytprim(i) = beta + ( grossism - beta ) * fhyp(i)
173	ENDDO
174
175	c ..... Calcul de Yf ........
176
177	Yf(0) = - pis2
178	DO i = 1, nmax2
179	yypr(i) = beta + ( grossism - beta ) * fxm(i)
180	ENDDO
181
182	DO i=1,nmax2
183	Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
184	ENDDO
185
186	c ****************************************************************
187	c
188	c ..... yuv = 0. si calcul des latitudes aux pts. U .....
189	c ..... yuv = 0.5 si calcul des latitudes aux pts. V .....
190	c
191	WRITE(6,18)
192	c
193	DO 5000 ik = 1,4
194
195	IF( ik.EQ.1 ) THEN
196	yuv = 0.
197	jlat = jjm + 1
198	ELSE IF ( ik.EQ.2 ) THEN
199	yuv = 0.5
200	jlat = jjm
201	ELSE IF ( ik.EQ.3 ) THEN
202	yuv = 0.25
203	jlat = jjm
204	ELSE IF ( ik.EQ.4 ) THEN
205	yuv = 0.75
206	jlat = jjm
207	ENDIF
208	c
209	yo1 = 0.
210	DO 1500 j = 1,jlat
211	yo1 = 0.
212	ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.)
213	yfi = ylon2
214	c
215	DO 250 it = nmax2,0,-1
216	IF( yfi.GE.Yf(it)) GO TO 350
217	250 CONTINUE
218	it = 0
219	350 CONTINUE
220
221	yi = yt(it)
222	IF(it.EQ.nmax2) THEN
223	it = nmax2 -1
224	Yf(it+1) = pis2
225	ENDIF
226	c .................................................................
227	c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yi)
228	c ..... et Y'(yi) .....
229	c .................................................................
230
231	CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),
232	, yt(it),yt(it+1) , a0,a1,a2,a3 )
233
234	Yf1 = Yf(it)
235	Yprimin = a1 + 2.* a2 * yi + 3.a3 yi *yi
236
237	DO 500 iter = 1,300
238	yi = yi - ( Yf1 - yfi )/ Yprimin
239
240	IF( ABS(yi-yo1).LE.epsilon) GO TO 550
241	yo1 = yi
242	yi2 = yi * yi
243	Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi
244	Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2
245	500 CONTINUE
246	WRITE(6,) ' Pas de solution **** ',j,ylon2,iter
247	STOP 2
248	550 CONTINUE
249	c
250	Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi
251	yprim(j) = pi / ( jjm * Yprimin )
252	yvrai(j) = yi
253
254	1500 CONTINUE
255
256	DO j = 1, jlat -1
257	IF( yvrai(j+1). LT. yvrai(j) ) THEN
258	WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j,
259	, ')'
260	STOP 3
261	ENDIF
262	ENDDO
263
264	WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
265	, ,' et pi/2 '
266	c
267	IF( ik.EQ.1 ) THEN
268	ypn = pis2
269	DO j = jlat,1,-1
270	IF( yvrai(j).LE. ypn ) GO TO 1502
271	ENDDO
272	1502 CONTINUE
273
274	jpn = j
275	y00 = yvrai(jpn)
276	deply = pis2 - y00
277	ENDIF
278
279	DO j = 1, jjm +1 - jpn
280	ylatt (j) = -pis2 - y00 + yvrai(jpn+j-1)
281	yprimm(j) = yprim(jpn+j-1)
282	ENDDO
283
284	jjpn = jpn
285	IF( jlat.EQ. jjm ) jjpn = jpn -1
286
287	DO j = 1,jjpn
288	ylatt (j + jjm+1 -jpn) = yvrai(j) + deply
289	yprimm(j + jjm+1 -jpn) = yprim(j)
290	ENDDO
291
292	c ********* Fin de la reorganisation ***********
293	c
294	1600 CONTINUE
295
296	DO j = 1, jlat
297	ylat(j) = ylatt( jlat +1 -j )
298	yprim(j) = yprimm( jlat +1 -j )
299	ENDDO
300
301	DO j = 1, jlat
302	yvrai(j) = ylat(j)*180./pi
303	ENDDO
304
305	IF( ik.EQ.1 ) THEN
306	c WRITE(6,18)
307	c WRITE(6,*) ' YLAT en U apres ( en deg. ) '
308	c WRITE(6,68) (yvrai(j),j=1,jlat)
309	cc WRITE(6,*) ' YPRIM '
310	cc WRITE(6,445) ( yprim(j),j=1,jlat)
311
312	DO j = 1, jlat
313	rrlatu(j) = ylat( j )
314	yyprimu(j) = yprim( j )
315	ENDDO
316
317	ELSE IF ( ik.EQ. 2 ) THEN
318	c WRITE(6,18)
319	c WRITE(6,*) ' YLAT en V apres ( en deg. ) '
320	c WRITE(6,68) (yvrai(j),j=1,jlat)
321	cc WRITE(6,*)' YPRIM '
322	cc WRITE(6,445) ( yprim(j),j=1,jlat)
323
324	DO j = 1, jlat
325	rrlatv(j) = ylat( j )
326	yyprimv(j) = yprim( j )
327	ENDDO
328
329	ELSE IF ( ik.EQ. 3 ) THEN
330	c WRITE(6,18)
331	c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) '
332	c WRITE(6,68) (yvrai(j),j=1,jlat)
333	cc WRITE(6,*) ' YPRIM '
334	cc WRITE(6,445) ( yprim(j),j=1,jlat)
335
336	DO j = 1, jlat
337	rlatu2(j) = ylat( j )
338	yprimu2(j) = yprim( j )
339	ENDDO
340
341	ELSE IF ( ik.EQ. 4 ) THEN
342	c WRITE(6,18)
343	c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) '
344	c WRITE(6,68)(yvrai(j),j=1,jlat)
345	cc WRITE(6,*) ' YPRIM '
346	cc WRITE(6,68) ( yprim(j),j=1,jlat)
347
348	DO j = 1, jlat
349	rlatu1(j) = ylat( j )
350	yprimu1(j) = yprim( j )
351	ENDDO
352
353	ENDIF
354
355	5000 CONTINUE
356	c
357	WRITE(6,18)
358	c
359	c ..... fin de la boucle do 5000 .....
360
361	DO j = 1, jjm
362	ylat(j) = rrlatu(j) - rrlatu(j+1)
363	ENDDO
364	champmin = 1.e12
365	champmax = -1.e12
366	DO j = 1, jjm
367	champmin = MIN( champmin, ylat(j) )
368	champmax = MAX( champmax, ylat(j) )
369	ENDDO
370	champmin = champmin * 180./pi
371	champmax = champmax * 180./pi
372
373	24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
374	18 FORMAT(/)
375	68 FORMAT(1x,7f9.2)
376
377	RETURN
378	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format