Context Navigation

fxhyp.F @ 2329

Last change on this file since 2329 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 11.7 KB

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

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

Download in other formats:

Original Format