Context Navigation

fxhyp.F @ 204

Last change on this file since 204 was 203, checked in by lmdz, 24 years ago
Mise en oeuvre du nouveau zoom. Modifs P. 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: 10.8 KB

Rev	Line
[203]	1	c
	2	c $Header
	3	c
[2]	4	SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
	5	, rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025)
	6
	7	c Auteur : P. Le Van
	8
	9	IMPLICIT NONE
	10
	11	c Calcule les longitudes et derivees dans la grille du GCM pour une
	12	c fonction f(x) a tangente hyperbolique .
	13	c
	14	c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
	15	c dzoom etant la distance totale de la zone du zoom
[203]	16	c tau la raideur de la transition de l'interieur a l'exterieur du zoom
[2]	17	c
[203]	18	c On doit avoir grossism x dzoom < pi ( radians ) , en longitude.
	19	c ********************************************************************
[2]	20
[203]	21
[2]	22	INTEGER nmax, nmax2
[203]	23	PARAMETER ( nmax = 30000, nmax2 = 2*nmax )
[2]	24
	25	#include "dimensions.h"
	26	#include "paramet.h"
	27
	28	c ...... arguments d'entree .......
	29	c
	30	REAL xzoomdeg,dzoom,tau,grossism
[203]	31
	32	c ...... arguments de sortie ......
	33
[2]	34	REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
	35	, rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
	36
[203]	37	c .... variables locales ....
[2]	38	c
[203]	39	REAL*8 xlon(iip1),xprimm(iip1),xuv
	40	REAL*8 xtild(0:nmax2)
	41	REAL*8 fhyp(0:nmax),ffdx,beta,Xprimt(0:nmax2)
	42	REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
	43	REAL*8 xvrai(iip1),xxprim(iip1)
	44	REAL*8 pi,depi,epsilon,xzoom,fa,fb
	45	REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
[2]	46	INTEGER i,it,ik,iter,ii,idif
[203]	47	REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
	48	REAL*8 champmin,champmax
[2]	49	INTEGER is2
	50	SAVE is2
[203]	51	REAL*8 heavyside
	52	EXTERNAL coefpoly,heavyside
[2]	53
	54	pi = 2. * ASIN(1.)
	55	depi = 2. * pi
[203]	56	epsilon = 1.e-3
[2]	57	xzoom = xzoomdeg * pi/180.
	58
[203]	59	WRITE(6,24) xzoomdeg,grossism,tau,dzoom
	60	IF( dzoom.LT.1.) THEN
	61	dzoom = dzoom * depi
	62	ELSEIF( dzoom.LT. 25. ) THEN
	63	WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
	64	,menter et relancer ! '
	65	STOP 1
	66	ELSE
	67	dzoom = dzoom * pi/180.
	68	ENDIF
[2]	69
	70	DO i = 0, nmax2
[203]	71	xtild(i) = - pi + FLOAT(i) * depi /nmax2
[2]	72	ENDDO
	73
[203]	74	DO i = 0, nmax2
	75
	76	fa = tau* ( dzoom/2. - xtild(i) )
	77	fb = xtild(i) * ( pi - xtild(i) )
	78
	79
	80	IF( 200.* fb .LT. - fa ) THEN
	81	fhyp ( i) = - 1.
	82	ELSEIF( 200. * fb .LT. fa ) THEN
	83	fhyp ( i) = 1.
	84	ELSE
	85	fhyp(i) = TANH ( fa/fb )
	86	ENDIF
	87
	88	IF ( xtild(i).EQ. 0. ) fhyp(i) = 1.
	89	IF ( xtild(i).EQ. pi ) fhyp(i) = -1.
	90
[2]	91	ENDDO
	92
	93	cc .... Calcul de beta ....
[203]	94	c ............................
[2]	95
[203]	96	ffdx = 0.
	97
	98	DO i = nmax +1,nmax2
	99
	100	xmoy = 0.5 * ( xtild(i-1) + xtild( i ) )
	101	fa = tau* ( dzoom/2. - xmoy )
	102	fb = xmoy * ( pi - xmoy )
	103
	104	IF( 200.* fb .LT. - fa ) THEN
	105	fxm = - 1.
	106	ELSEIF( 200. * fb .LT. fa ) THEN
	107	fxm = 1.
	108	ELSE
	109	fxm = TANH ( fa/fb )
	110	ENDIF
	111
	112	IF ( xmoy.EQ. 0. ) fhyp(i) = 1.
	113	IF ( xmoy.EQ. pi ) fhyp(i) = -1.
	114	ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
	115
[2]	116	ENDDO
	117
[203]	118	beta = ( grossism * ffdx - pi ) / ( ffdx - pi )
	119	ccc WRITE(6,) ' * X beta **',beta,grossism
	120
	121	IF( 2.*beta - grossism.LE. 0.) THEN
	122	WRITE(6,) ' * Attention ! La valeur beta calculee dans la rou
	123	,tine fxhyp est mauvaise ! '
	124	WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ',
	125	, ' et relancer ! *** '
	126	CALL ABORT
	127	ENDIF
[2]	128	c
	129	c ..... calcul de Xprimt .....
	130	c
	131
[203]	132	DO i = nmax, nmax2
[2]	133	Xprimt(i) = beta + ( grossism - beta ) * fhyp(i)
	134	ENDDO
	135	c
[203]	136	DO i = nmax+1, nmax2
[2]	137	Xprimt( nmax2 - i ) = Xprimt( i )
	138	ENDDO
	139	c
	140
	141	c ..... Calcul de Xf ........
	142
[203]	143	Xf(0) = - pi
	144
	145	DO i = nmax +1, nmax2
	146
	147	xmoy = 0.5 * ( xtild(i-1) + xtild( i ) )
	148	fa = tau* ( dzoom/2. - xmoy )
	149	fb = xmoy * ( pi - xmoy )
	150
	151	IF( 200.* fb .LT. - fa ) THEN
	152	fxm = - 1.
	153	ELSEIF( 200. * fb .LT. fa ) THEN
	154	fxm = 1.
	155	ELSE
	156	fxm = TANH ( fa/fb )
	157	ENDIF
	158
	159	IF ( xmoy.EQ. 0. ) fxm = 1.
	160	IF ( xmoy.EQ. pi ) fxm = -1.
	161	xxpr(i) = beta + ( grossism - beta ) * fxm
	162
[2]	163	ENDDO
	164
[203]	165	DO i = nmax+1, nmax2
	166	xxpr(nmax2-i+1) = xxpr(i)
[2]	167	ENDDO
	168
	169	DO i=1,nmax2
	170	Xf(i) = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
	171	ENDDO
	172
	173
	174	c *****************************************************************
	175	c
	176
	177	c ..... xuv = 0. si calcul aux pts scalaires ........
	178	c ..... xuv = 0.5 si calcul aux pts U ........
	179	c
	180	c
	181	DO 5000 ik = 1, 4
	182
	183	IF( ik.EQ.1 ) THEN
	184	xuv = - 0.25
	185	ELSE IF ( ik.EQ.2 ) THEN
	186	xuv = 0.
	187	ELSE IF ( ik.EQ.3 ) THEN
	188	xuv = 0.5
	189	ELSE IF ( ik.EQ.4 ) THEN
	190	xuv = 0.25
	191	ENDIF
	192
[203]	193	xo1 = 0.
[2]	194
	195	DO 1500 i = 1, iim
	196
[203]	197	xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim)
[2]	198
[203]	199	Xfi = xlon2
[2]	200	c
	201	DO 250 it = nmax2,0,-1
[203]	202	IF( Xfi.GE.Xf(it)) GO TO 350
[2]	203	250 CONTINUE
	204
	205	it = 0
	206
	207	350 CONTINUE
	208
[203]	209	c ...... Calcul de Xf(xi) ......
	210	c
	211	xi = xtild(it)
[2]	212
[203]	213	IF(it.EQ.nmax2) THEN
	214	it = nmax2 -1
	215	Xf(it+1) = pi
	216	ENDIF
	217	c .....................................................................
	218	c
	219	c Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
	220	c polynome de degre 3 qui passe par les points (Xf(it),xtild(it) )
	221	c et (Xf(it+1),xtild(it+1) )
[2]	222
[203]	223	CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
	224	, xtild(it),xtild(it+1), a0, a1, a2, a3 )
	225
	226	Xf1 = Xf(it)
	227	Xprimin = a1 + 2.* a2 * xi + 3.a3 xi *xi
	228
	229	DO 500 iter = 1,300
	230	xi = xi - ( Xf1 - Xfi )/ Xprimin
	231
	232	IF( ABS(xi-xo1).LE.epsilon) GO TO 550
	233	xo1 = xi
	234	xi2 = xi * xi
	235	Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
	236	Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
[2]	237	500 CONTINUE
[203]	238	WRITE(6,) ' Pas de solution **** ',i,xlon2,iter
	239	STOP 6
[2]	240	550 CONTINUE
	241
[203]	242	xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
	243	xvrai(i) = xi + xzoom
[2]	244
	245	1500 CONTINUE
	246
	247	DO i = 1 , iim
[203]	248	xlon(i) = xvrai(i)
[2]	249	xprimm(i) = xxprim(i)
	250	ENDDO
	251
	252	DO i = 1, iim -1
	253	IF( xvrai(i+1). LT. xvrai(i) ) THEN
[203]	254	WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
	255	, ')'
	256	STOP 7
[2]	257	ENDIF
	258	ENDDO
	259	c
	260	c ... Reorganisation des longitudes pour les avoir entre - pi et pi ..
	261	c ........................................................................
	262
	263	champmin = 1.e12
	264	champmax = -1.e12
	265	DO i = 1, iim
	266	champmin = MIN( champmin,xvrai(i) )
	267	champmax = MAX( champmax,xvrai(i) )
	268	ENDDO
	269
	270	IF(champmin .GE. - pi .AND. champmax.LE. pi ) THEN
	271	GO TO 1600
	272	ELSE
[203]	273	WRITE(6,18)
	274	WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
	275	, ' et pi '
[2]	276	c
	277	IF( xzoom.LE.0.) THEN
	278	IF( ik.EQ. 1 ) THEN
	279	DO i = 1, iim
	280	IF( xvrai(i).GE. - pi ) GO TO 80
	281	ENDDO
[203]	282	WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! '
	283	STOP 8
[2]	284	80 CONTINUE
	285	is2 = i
	286	ENDIF
	287
	288	IF( is2.NE. 1 ) THEN
	289	DO ii = is2 , iim
	290	xlon (ii-is2+1) = xvrai(ii)
	291	xprimm(ii-is2+1) = xxprim(ii)
	292	ENDDO
	293	DO ii = 1 , is2 -1
	294	xlon (ii+iim-is2+1) = xvrai(ii) + depi
	295	xprimm(ii+iim-is2+1) = xxprim(ii)
	296	ENDDO
	297	ENDIF
	298	ELSE
	299	IF( ik.EQ.1 ) THEN
	300	DO i = iim,1,-1
	301	IF( xvrai(i).LE. pi ) GO TO 90
	302	ENDDO
[203]	303	WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! '
	304	STOP 9
[2]	305	90 CONTINUE
	306	is2 = i
	307	ENDIF
	308	idif = iim -is2
	309	DO ii = 1, is2
	310	xlon (ii+idif) = xvrai(ii)
	311	xprimm(ii+idif) = xxprim(ii)
	312	ENDDO
	313	DO ii = 1, idif
	314	xlon (ii) = xvrai (ii+is2) - depi
	315	xprimm(ii) = xxprim(ii+is2)
	316	ENDDO
	317	ENDIF
	318	ENDIF
	319	c
	320	c ......... Fin de la reorganisation ............................
	321
	322	1600 CONTINUE
	323
	324
	325	xlon ( iip1) = xlon(1) + depi
	326	xprimm( iip1 ) = xprimm (1 )
	327
	328	DO i = 1, iim+1
	329	xvrai(i) = xlon(i)*180./pi
	330	ENDDO
	331
[203]	332	IF( ik.EQ.1 ) THEN
	333	WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) '
	334	WRITE(6,18)
	335	WRITE(6,68) xvrai
	336	ccc WRITE(6,*) ' XPRIM k ',ik
	337	ccc WRITE(6,566) xprimm
[2]	338
	339	DO i = 1,iim + 1
	340	rlonm025(i) = xlon( i )
	341	xprimm025(i) = xprimm(i)
	342	ENDDO
[203]	343
[2]	344	ELSE IF( ik.EQ.2 ) THEN
[203]	345	WRITE(6,18)
	346	WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) '
	347	WRITE(6,68) xvrai
	348	ccc WRITE(6,*) ' XPRIM k ',ik
	349	ccc WRITE(6,566) xprimm
	350
[2]	351	DO i = 1,iim + 1
	352	rlonv(i) = xlon( i )
	353	xprimv(i) = xprimm(i)
	354	ENDDO
[203]	355
	356	ELSE IF( ik.EQ.3) THEN
	357	WRITE(6,18)
	358	WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) '
	359	WRITE(6,68) xvrai
	360	ccc WRITE(6,*) ' XPRIM ik ',ik
	361	ccc WRITE(6,566) xprimm
	362
[2]	363	DO i = 1,iim + 1
	364	rlonu(i) = xlon( i )
	365	xprimu(i) = xprimm(i)
	366	ENDDO
[203]	367
[2]	368	ELSE IF( ik.EQ.4 ) THEN
[203]	369	WRITE(6,18)
	370	WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) '
	371	WRITE(6,68) xvrai
	372	ccc WRITE(6,*) ' XPRIM ik ',ik
	373	ccc WRITE(6,566) xprimm
	374
[2]	375	DO i = 1,iim + 1
	376	rlonp025(i) = xlon( i )
	377	xprimp025(i) = xprimm(i)
	378	ENDDO
[203]	379
[2]	380	ENDIF
	381
	382	5000 CONTINUE
	383	c
[203]	384	c ........... fin de la boucle do 5000 ............
[2]	385
[203]	386	DO i = 1, iim
	387	xlon(i) = rlonv(i+1) - rlonv(i)
[2]	388	ENDDO
[203]	389	champmin = 1.e12
	390	champmax = -1.e12
[2]	391	DO i = 1, iim
[203]	392	champmin = MIN( champmin, xlon(i) )
	393	champmax = MAX( champmax, xlon(i) )
[2]	394	ENDDO
[203]	395	champmin = champmin * 180./pi
	396	champmax = champmax * 180./pi
[2]	397
[203]	398	WRITE(6,18)
	399	WRITE(6,*) ' Longitudes '
	400	WRITE(6,18)
	401	WRITE(6,3) champmin, champmax
	402	WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
	403	,ametres grossism , tau , dzoom pour X et repasser ! '
	404	WRITE(6,18)
[2]	405
[203]	406	3 Format(1x, ' Au centre du zoom , la longueur de la maille est',
	407	, ' d environ ',f8.2 ,' degres ',
	408	, ' alors que la maille en dehors de la zone du zoom est d environ
	409	, ', f8.2,' degres ' )
	410	18 FORMAT(/)
	411	24 FORMAT(2x,' Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
	412	68 FORMAT(1x,7f9.2)
[2]	413
[203]	414	RETURN
	415	END

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

Download in other formats:

Original Format