Context Navigation

fyhyp.F @ 212

Last change on this file since 212 was 212, checked in by lmdz, 24 years ago
Modifications des sorties de controle Levan LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 9.9 KB

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

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

Download in other formats:

Original Format