Context Navigation

fyhyp.F @ 4220

Last change on this file since 4220 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

File size: 10.1 KB

Rev	Line
[524]	1	!
[1403]	2	! $Id: fyhyp.F 1403 2010-07-01 09:02:53Z lguez $
[524]	3	!
	4	c
	5	c
[650]	6	SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau ,
[524]	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
[650]	39	REAL yzoomdeg, grossism,dzooma,tau
[524]	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
[650]	51	REAL dzoom
[1279]	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)
[524]	56	SAVE Ytprim, yt,Yf
[1279]	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
[524]	63	SAVE y00, deply
	64
	65	INTEGER i,j,it,ik,iter,jlat
	66	INTEGER jpn,jjpn
	67	SAVE jpn
[1279]	68	REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m
	69	REAL(KIND=8) fa(0:nmax2),fb(0:nmax2)
[524]	70	REAL y0min,y0max
	71
[1279]	72	REAL(KIND=8) heavyside
[524]	73
	74	pi = 2. * ASIN(1.)
	75	depi = 2. * pi
	76	pis2 = pi/2.
[1403]	77	pisjm = pi/ REAL(jjm)
[524]	78	epsilon = 1.e-3
	79	y0 = yzoomdeg * pi/180.
	80
[650]	81	IF( dzooma.LT.1.) THEN
	82	dzoom = dzooma * pi
	83	ELSEIF( dzooma.LT. 12. ) THEN
[524]	84	WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
	85	,menter et relancer ! '
	86	STOP 1
	87	ELSE
[650]	88	dzoom = dzooma * pi/180.
[524]	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
[1403]	96	yt(i) = - pis2 + REAL(i)* pi /nmax2
[524]	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.
[1403]	212	ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.)
[524]	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