Context Navigation

fxhyp.F @ 285

Last change on this file since 285 was 251, checked in by lmdz, 23 years ago
Initialisations pour eviter des divisions par 0 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.5 KB

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

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

Download in other formats:

Original Format