Context Navigation

inigeom.F @ 1721

Last change on this file since 1721 was 1658, checked in by Laurent Fairhead, 12 years ago

Phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ4 (HEAD)
Validation effectuee par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modele donnee par Y. Meurdesoif et la version phasee avec la r1428 (fin du tronc LMDZ4)

Phasing of the localised (low memory) parallel dynamics package with the LMDZ4 trunk version of LMDZ
Validation consisted in comparing output debug files (u, v, t, q, masse, etc... ) of a no physics simulation
run with the version of the code given by Y. Meurdesoif and this version phased with r1428 (HEAD of the LMDZ4 trunk)

File size: 21.8 KB

Rev	Line
[1632]	1	!
[1658]	2	! $Id: inigeom.F 1403 2010-07-01 09:02:53Z fairhead $
[1632]	3	!
	4	c
	5	c
	6	SUBROUTINE inigeom
	7	c
	8	c Auteur : P. Le Van
	9	c
	10	c ............ Version du 01/04/2001 ........................
	11	c
	12	c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en-
	13	c endroits que les aires aireij1,..aireij4 .
	14
	15	c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
	16	c
	17	c
	18	IMPLICIT NONE
	19	c
	20	#include "dimensions.h"
	21	#include "paramet.h"
	22	#include "comconst.h"
	23	#include "comgeom2.h"
	24	#include "serre.h"
	25	#include "logic.h"
	26	#include "comdissnew.h"
	27
	28	c-----------------------------------------------------------------------
	29	c .... Variables locales ....
	30	c
	31	INTEGER i,j,itmax,itmay,iter
	32	REAL cvu(iip1,jjp1),cuv(iip1,jjm)
	33	REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
	34	REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
	35	REAL coslatm,coslatp,radclatm,radclatp
	36	REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
	37	* cuij4(iip1,jjp1)
	38	REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
	39	* cvij4(iip1,jjp1)
	40	REAL rlonvv(iip1),rlatuu(jjp1)
	41	REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
	42	* yprimv(jjm),yprimu(jjp1)
	43	REAL gamdi_gdiv, gamdi_grot, gamdi_h
	44
	45	REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
	46	, xprimp025(iip1)
	47	SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
	48	SAVE rlonm025,xprimm025,rlonp025,xprimp025
	49
	50	REAL SSUM
	51	c
	52	c
	53	c ------------------------------------------------------------------
	54	c - -
	55	c - calcul des coeff. ( cu, cv , 1./cu2, 1./cv2 ) -
	56	c - -
	57	c ------------------------------------------------------------------
	58	c
	59	c les coef. ( cu, cv ) permettent de passer des vitesses naturelles
	60	c aux vitesses covariantes et contravariantes , ou vice-versa ...
	61	c
	62	c
	63	c on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu
	64	c v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv
	65	c
	66	c on en tire : u(covariant) = cu * cu * u(contravariant)
	67	c v(covariant) = cv * cv * v(contravariant)
	68	c
	69	c
	70	c on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
	71	c = =
	72	c et - jm/2 < Y < jm/2
	73	c = =
	74	c
	75	c ...................................................
	76	c ...................................................
	77	c . x est la longitude du point en radians .
	78	c . y est la latitude du point en radians .
	79	c . .
	80	c . on a : cu(i,j) = rad * COS(y) * dx/dX .
	81	c . cv( j ) = rad * dy/dY .
	82	c . aire(i,j) = cu(i,j) * cv(j) .
	83	c . .
	84	c . y, dx/dX, dy/dY calcules aux points concernes .
	85	c . .
	86	c ...................................................
	87	c ...................................................
	88	c
	89	c
	90	c
	91	c ,
	92	c cv , bien que dependant de j uniquement,sera ici indice aussi en i
	93	c pour un adressage plus facile en ij .
	94	c
	95	c
	96	c
	97	c ************ aux points u et v , ***************
	98	c xprimu et xprimv sont respectivement les valeurs de dx/dX
	99	c yprimu et yprimv . . . . . . . . . . . dy/dY
	100	c rlatu et rlatv . . . . . . . . . . .la latitude
	101	c cvu et cv . . . . . . . . . . . cv
	102	c
	103	c ************ aux points u, v, scalaires, et z **************
	104	c cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu
	105	c
	106	c
	107	c
	108	c Exemple de distribution de variables sur la grille dans le
	109	c domaine de travail ( X,Y ) .
	110	c ................................................................
	111	c DX=DY= 1
	112	c
	113	c
	114	c + represente un point scalaire ( p.exp la pression )
	115	c > represente la composante zonale du vent
	116	c V represente la composante meridienne du vent
	117	c o represente la vorticite
	118	c
	119	c ---- , car aux poles , les comp.zonales covariantes sont nulles
	120	c
	121	c
	122	c
	123	c i ->
	124	c
	125	c 1 2 3 4 5 6 7 8
	126	c j
	127	c v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
	128	c
	129	c V o V o V o V o V o V o V o V o
	130	c
	131	c 2 + > + > + > + > + > + > + > + >
	132	c
	133	c V o V o V o V o V o V o V o V o
	134	c
	135	c 3 + > + > + > + > + > + > + > + >
	136	c
	137	c V o V o V o V o V o V o V o V o
	138	c
	139	c 4 + > + > + > + > + > + > + > + >
	140	c
	141	c V o V o V o V o V o V o V o V o
	142	c
	143	c 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
	144	c
	145	c
	146	c Ci-dessus, on voit que le nombre de pts.en longitude est egal
	147	c a IM = 8
	148	c De meme , le nombre d'intervalles entre les 2 poles est egal
	149	c a JM = 4
	150	c
	151	c Les points scalaires ( + ) correspondent donc a des valeurs
	152	c entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
	153	c
	154	c Les vents U ( > ) correspondent a des valeurs semi-
	155	c entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
	156	c
	157	c Les vents V ( V ) correspondent a des valeurs entieres
	158	c de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
	159	c
	160	c
	161	c
	162	WRITE(6,3)
	163	3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ',
	164	* //5x,' Calcul des elongations cu et cv comme sommes des 4 ' /
	165	* 5x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux
	166	* '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / )
	167	c
	168	c
	169	IF( nitergdiv.NE.2 ) THEN
	170	gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
	171	ELSE
	172	gamdi_gdiv = 0.
	173	ENDIF
	174	IF( nitergrot.NE.2 ) THEN
	175	gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
	176	ELSE
	177	gamdi_grot = 0.
	178	ENDIF
	179	IF( niterh.NE.2 ) THEN
	180	gamdi_h = coefdis/ ( REAL(niterh) -2. )
	181	ELSE
	182	gamdi_h = 0.
	183	ENDIF
	184
	185	WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
	186	* nitergdiv,nitergrot,niterh
	187	c
	188	pi = 2.* ASIN(1.)
	189	c
	190	WRITE(6,990)
	191
	192	c ----------------------------------------------------------------
	193	c
	194	IF( .NOT.fxyhypb ) THEN
	195	c
	196	c
	197	IF( ysinus ) THEN
	198	c
	199	WRITE(6,) ' Inigeom , Y = Sinus ( Latitude ) * '
	200	c
	201	c .... utilisation de f(x,y ) avec y = sinus de la latitude .....
	202
	203	CALL fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
	204	, rlatu2,yprimu2,
	205	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
	206
	207	ELSE
	208	c
	209	WRITE(6,) ' Inigeom , Y = Latitude , der. sinusoid . *'
	210
	211	c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
	212	c
	213
	214	pxo = clon *pi /180.
	215	pyo = 2.* clat* pi /180.
	216	c
	217	c .... determination de transx ( pour le zoom ) par Newton-Raphson ...
	218	c
	219	itmax = 10
	220	eps = .1e-7
	221	c
	222	xo1 = 0.
	223	DO 10 iter = 1, itmax
	224	x1 = xo1
	225	f = x1+ alphax *SIN(x1-pxo)
	226	df = 1.+ alphax *COS(x1-pxo)
	227	x1 = x1 - f/df
	228	xdm = ABS( x1- xo1 )
	229	IF( xdm.LE.eps )GO TO 11
	230	xo1 = x1
	231	10 CONTINUE
	232	11 CONTINUE
	233	c
	234	transx = xo1
	235
	236	itmay = 10
	237	eps = .1e-7
	238	C
	239	yo1 = 0.
	240	DO 15 iter = 1,itmay
	241	y1 = yo1
	242	f = y1 + alphay* SIN(y1-pyo)
	243	df = 1. + alphay* COS(y1-pyo)
	244	y1 = y1 -f/df
	245	ydm = ABS(y1-yo1)
	246	IF(ydm.LE.eps) GO TO 17
	247	yo1 = y1
	248	15 CONTINUE
	249	c
	250	17 CONTINUE
	251	transy = yo1
	252
	253	CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
	254	, rlatu2,yprimu2,
	255	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
	256
	257	ENDIF
	258	c
	259	ELSE
	260	c
	261	c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
	262	c .....................................................................
	263
	264	WRITE(6,)' Inigeom , Y = Latitude , der.tg. hyperbolique *'
	265
	266	CALL fxyhyper( clat, grossismy, dzoomy, tauy ,
	267	, clon, grossismx, dzoomx, taux ,
	268	, rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 ,
	269	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
	270
	271
	272	ENDIF
	273	c
	274	c -------------------------------------------------------------------
	275
	276	c
	277	rlatu(1) = ASIN(1.)
	278	rlatu(jjp1) = - rlatu(1)
	279	c
	280	c
	281	c .... calcul aux poles ....
	282	c
	283	yprimu(1) = 0.
	284	yprimu(jjp1) = 0.
	285	c
	286	c
	287	un4rad2 = 0.25 * rad * rad
	288	c
	289	c --------------------------------------------------------------------
	290	c --------------------------------------------------------------------
	291	c - -
	292	c - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) -
	293	c - et de fext , force de coriolis extensive . -
	294	c - -
	295	c --------------------------------------------------------------------
	296	c --------------------------------------------------------------------
	297	c
	298	c
	299	c
	300	c A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
	301	c affectees 4 aires entourant P , calculees respectivement aux points
	302	c ( i + 1/4, j - 1/4 ) : aireij1 (i,j)
	303	c ( i + 1/4, j + 1/4 ) : aireij2 (i,j)
	304	c ( i - 1/4, j + 1/4 ) : aireij3 (i,j)
	305	c ( i - 1/4, j - 1/4 ) : aireij4 (i,j)
	306	c
	307	c ,
	308	c Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
	309	c Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
	310	c des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au
	311	c point (i,j) .
	312	c On definit en outre les coefficients alpha comme etant egaux a
	313	c (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j)
	314	c
	315	c De meme, toute aire centree en 1 point U est egale a la somme des
	316	c 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
	317	c Idem pour airev, airez .
	318	c
	319	c On a ,pour chaque maille : dX = dY = 1
	320	c
	321	c
	322	c . V
	323	c
	324	c aireij4 . . aireij1
	325	c
	326	c U . . P . U
	327	c
	328	c aireij3 . . aireij2
	329	c
	330	c . V
	331	c
	332	c
	333	c
	334	c
	335	c
	336	c ....................................................................
	337	c
	338	c Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
	339	c qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
	340	c taires cuij et les 4 elongat. cvij qui sont calculees aux memes
	341	c endroits que les aireij .
	342	c
	343	c ....................................................................
	344	c
	345	c ....... do 35 : boucle sur les jjm + 1 latitudes .....
	346	c
	347	c
	348	DO 35 j = 1, jjp1
	349	c
	350	IF ( j. eq. 1 ) THEN
	351	c
	352	yprm = yprimu1(j)
	353	rlatm = rlatu1(j)
	354	c
	355	coslatm = COS( rlatm )
	356	radclatm = 0.5* rad * coslatm
	357	c
	358	DO 30 i = 1, iim
	359	xprp = xprimp025( i )
	360	xprm = xprimm025( i )
	361	aireij2( i,1 ) = un4rad2 * coslatm * xprp * yprm
	362	aireij3( i,1 ) = un4rad2 * coslatm * xprm * yprm
	363	cuij2 ( i,1 ) = radclatm * xprp
	364	cuij3 ( i,1 ) = radclatm * xprm
	365	cvij2 ( i,1 ) = 0.5* rad * yprm
	366	cvij3 ( i,1 ) = cvij2(i,1)
	367	30 CONTINUE
	368	c
	369	DO i = 1, iim
	370	aireij1( i,1 ) = 0.
	371	aireij4( i,1 ) = 0.
	372	cuij1 ( i,1 ) = 0.
	373	cuij4 ( i,1 ) = 0.
	374	cvij1 ( i,1 ) = 0.
	375	cvij4 ( i,1 ) = 0.
	376	ENDDO
	377	c
	378	END IF
	379	c
	380	IF ( j. eq. jjp1 ) THEN
	381	yprp = yprimu2(j-1)
	382	rlatp = rlatu2 (j-1)
	383	ccc yprp = fyprim( REAL(j) - 0.25 )
	384	ccc rlatp = fy ( REAL(j) - 0.25 )
	385	c
	386	coslatp = COS( rlatp )
	387	radclatp = 0.5* rad * coslatp
	388	c
	389	DO 31 i = 1,iim
	390	xprp = xprimp025( i )
	391	xprm = xprimm025( i )
	392	aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp
	393	aireij4( i,jjp1 ) = un4rad2 * coslatp * xprm * yprp
	394	cuij1(i,jjp1) = radclatp * xprp
	395	cuij4(i,jjp1) = radclatp * xprm
	396	cvij1(i,jjp1) = 0.5 * rad* yprp
	397	cvij4(i,jjp1) = cvij1(i,jjp1)
	398	31 CONTINUE
	399	c
	400	DO i = 1, iim
	401	aireij2( i,jjp1 ) = 0.
	402	aireij3( i,jjp1 ) = 0.
	403	cvij2 ( i,jjp1 ) = 0.
	404	cvij3 ( i,jjp1 ) = 0.
	405	cuij2 ( i,jjp1 ) = 0.
	406	cuij3 ( i,jjp1 ) = 0.
	407	ENDDO
	408	c
	409	END IF
	410	c
	411
	412	IF ( j .gt. 1 .AND. j .lt. jjp1 ) THEN
	413	c
	414	rlatp = rlatu2 ( j-1 )
	415	yprp = yprimu2( j-1 )
	416	rlatm = rlatu1 ( j )
	417	yprm = yprimu1( j )
	418	cc rlatp = fy ( REAL(j) - 0.25 )
	419	cc yprp = fyprim( REAL(j) - 0.25 )
	420	cc rlatm = fy ( REAL(j) + 0.25 )
	421	cc yprm = fyprim( REAL(j) + 0.25 )
	422
	423	coslatm = COS( rlatm )
	424	coslatp = COS( rlatp )
	425	radclatp = 0.5* rad * coslatp
	426	radclatm = 0.5* rad * coslatm
	427	c
	428	DO 32 i = 1,iim
	429	xprp = xprimp025( i )
	430	xprm = xprimm025( i )
	431
	432	ai14 = un4rad2 * coslatp * yprp
	433	ai23 = un4rad2 * coslatm * yprm
	434	aireij1 ( i,j ) = ai14 * xprp
	435	aireij2 ( i,j ) = ai23 * xprp
	436	aireij3 ( i,j ) = ai23 * xprm
	437	aireij4 ( i,j ) = ai14 * xprm
	438	cuij1 ( i,j ) = radclatp * xprp
	439	cuij2 ( i,j ) = radclatm * xprp
	440	cuij3 ( i,j ) = radclatm * xprm
	441	cuij4 ( i,j ) = radclatp * xprm
	442	cvij1 ( i,j ) = 0.5* rad * yprp
	443	cvij2 ( i,j ) = 0.5* rad * yprm
	444	cvij3 ( i,j ) = cvij2(i,j)
	445	cvij4 ( i,j ) = cvij1(i,j)
	446	32 CONTINUE
	447	c
	448	END IF
	449	c
	450	c ........ periodicite ............
	451	c
	452	cvij1 (iip1,j) = cvij1 (1,j)
	453	cvij2 (iip1,j) = cvij2 (1,j)
	454	cvij3 (iip1,j) = cvij3 (1,j)
	455	cvij4 (iip1,j) = cvij4 (1,j)
	456	cuij1 (iip1,j) = cuij1 (1,j)
	457	cuij2 (iip1,j) = cuij2 (1,j)
	458	cuij3 (iip1,j) = cuij3 (1,j)
	459	cuij4 (iip1,j) = cuij4 (1,j)
	460	aireij1 (iip1,j) = aireij1 (1,j )
	461	aireij2 (iip1,j) = aireij2 (1,j )
	462	aireij3 (iip1,j) = aireij3 (1,j )
	463	aireij4 (iip1,j) = aireij4 (1,j )
	464
	465	35 CONTINUE
	466	c
	467	c ..............................................................
	468	c
	469	DO 37 j = 1, jjp1
	470	DO 36 i = 1, iim
	471	aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
	472	* aireij4(i,j)
	473	alpha1 ( i,j ) = aireij1(i,j) / aire(i,j)
	474	alpha2 ( i,j ) = aireij2(i,j) / aire(i,j)
	475	alpha3 ( i,j ) = aireij3(i,j) / aire(i,j)
	476	alpha4 ( i,j ) = aireij4(i,j) / aire(i,j)
	477	alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j)
	478	alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j)
	479	alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j)
	480	alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j)
	481	36 CONTINUE
	482	c
	483	c
	484	aire (iip1,j) = aire (1,j)
	485	alpha1 (iip1,j) = alpha1 (1,j)
	486	alpha2 (iip1,j) = alpha2 (1,j)
	487	alpha3 (iip1,j) = alpha3 (1,j)
	488	alpha4 (iip1,j) = alpha4 (1,j)
	489	alpha1p2(iip1,j) = alpha1p2(1,j)
	490	alpha1p4(iip1,j) = alpha1p4(1,j)
	491	alpha2p3(iip1,j) = alpha2p3(1,j)
	492	alpha3p4(iip1,j) = alpha3p4(1,j)
	493	37 CONTINUE
	494	c
	495
	496	DO 42 j = 1,jjp1
	497	DO 41 i = 1,iim
	498	aireu (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
	499	* aireij3(i+1,j)
	500	unsaire ( i,j)= 1./ aire(i,j)
	501	unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv )
	502	unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h )
	503	airesurg ( i,j)= aire(i,j)/ g
	504	41 CONTINUE
	505	aireu (iip1,j) = aireu (1,j)
	506	unsaire (iip1,j) = unsaire(1,j)
	507	unsair_gam1(iip1,j) = unsair_gam1(1,j)
	508	unsair_gam2(iip1,j) = unsair_gam2(1,j)
	509	airesurg (iip1,j) = airesurg(1,j)
	510	42 CONTINUE
	511	c
	512	c
	513	DO 48 j = 1,jjm
	514	c
	515	DO i=1,iim
	516	airev (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
	517	* aireij4(i,j+1)
	518	ENDDO
	519	DO i=1,iim
	520	airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) +
	521	* aireij4(i+1,j+1)
	522	unsairez(i,j) = 1./ airez
	523	unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot )
	524	fext (i,j) = airez * SIN(rlatv(j))* 2.* omeg
	525	ENDDO
	526	airev (iip1,j) = airev(1,j)
	527	unsairez (iip1,j) = unsairez(1,j)
	528	fext (iip1,j) = fext(1,j)
	529	unsairz_gam(iip1,j) = unsairz_gam(1,j)
	530	c
	531	48 CONTINUE
	532	c
	533	c
	534	c ..... Calcul des elongations cu,cv, cvu .........
	535	c
	536	DO j = 1, jjm
	537	DO i = 1, iim
	538	cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1))
	539	cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) )
	540	cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1))
	541	unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) )
	542	ENDDO
	543	DO i = 1, iim
	544	cuvsurcv (i,j) = airev(i,j) * unscv2(i,j)
	545	cvsurcuv (i,j) = 1./cuvsurcv(i,j)
	546	cuvscvgam1(i,j) = cuvsurcv (i,j) ** ( - gamdi_gdiv )
	547	cuvscvgam2(i,j) = cuvsurcv (i,j) ** ( - gamdi_h )
	548	cvscuvgam(i,j) = cvsurcuv (i,j) ** ( - gamdi_grot )
	549	ENDDO
	550	cv (iip1,j) = cv (1,j)
	551	cvu (iip1,j) = cvu (1,j)
	552	unscv2 (iip1,j) = unscv2 (1,j)
	553	cuv (iip1,j) = cuv (1,j)
	554	cuvsurcv (iip1,j) = cuvsurcv (1,j)
	555	cvsurcuv (iip1,j) = cvsurcuv (1,j)
	556	cuvscvgam1(iip1,j) = cuvscvgam1(1,j)
	557	cuvscvgam2(iip1,j) = cuvscvgam2(1,j)
	558	cvscuvgam(iip1,j) = cvscuvgam(1,j)
	559	ENDDO
	560
	561	DO j = 2, jjm
	562	DO i = 1, iim
	563	cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j))
	564	unscu2 (i,j) = 1./ ( cu(i,j) * cu(i,j) )
	565	cvusurcu (i,j) = aireu(i,j) * unscu2(i,j)
	566	cusurcvu (i,j) = 1./ cvusurcu(i,j)
	567	cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv )
	568	cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h )
	569	cuscvugam (i,j) = cusurcvu(i,j) ** ( - gamdi_grot )
	570	ENDDO
	571	cu (iip1,j) = cu(1,j)
	572	unscu2 (iip1,j) = unscu2(1,j)
	573	cvusurcu (iip1,j) = cvusurcu(1,j)
	574	cusurcvu (iip1,j) = cusurcvu(1,j)
	575	cvuscugam1(iip1,j) = cvuscugam1(1,j)
	576	cvuscugam2(iip1,j) = cvuscugam2(1,j)
	577	cuscvugam (iip1,j) = cuscvugam(1,j)
	578	ENDDO
	579
	580	c
	581	c .... calcul aux poles ....
	582	c
	583	DO i = 1, iip1
	584	cu ( i, 1 ) = 0.
	585	unscu2( i, 1 ) = 0.
	586	cvu ( i, 1 ) = 0.
	587	c
	588	cu (i, jjp1) = 0.
	589	unscu2(i, jjp1) = 0.
	590	cvu (i, jjp1) = 0.
	591	ENDDO
	592	c
	593	c ..............................................................
	594	c
	595	DO j = 1, jjm
	596	DO i= 1, iim
	597	airvscu2 (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) )
	598	aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot )
	599	ENDDO
	600	airvscu2 (iip1,j) = airvscu2(1,j)
	601	aivscu2gam(iip1,j) = aivscu2gam(1,j)
	602	ENDDO
	603
	604	DO j=2,jjm
	605	DO i=1,iim
	606	airuscv2 (i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) )
	607	aiuscv2gam (i,j) = airuscv2(i,j)** ( - gamdi_grot )
	608	ENDDO
	609	airuscv2 (iip1,j) = airuscv2 (1,j)
	610	aiuscv2gam(iip1,j) = aiuscv2gam(1,j)
	611	ENDDO
	612
	613	c
	614	c calcul des aires aux poles :
	615	c -----------------------------
	616	c
	617	apoln = SSUM(iim,aire(1,1),1)
	618	apols = SSUM(iim,aire(1,jjp1),1)
	619	unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
	620	unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
	621	unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) )
	622	unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) )
	623	c
	624	c-----------------------------------------------------------------------
	625	c gtitre='Coriolis version ancienne'
	626	c gfichier='fext1'
	627	c CALL writestd(fext,iip1*jjm)
	628	c
	629	c changement F. Hourdin calcul conservatif pour fext
	630	c constang contient le produit a * cos ( latitude ) * omega
	631	c
	632	DO i=1,iim
	633	constang(i,1) = 0.
	634	ENDDO
	635	DO j=1,jjm-1
	636	DO i=1,iim
	637	constang(i,j+1) = radomegcu(i,j+1)*COS(rlatu(j+1))
	638	ENDDO
	639	ENDDO
	640	DO i=1,iim
	641	constang(i,jjp1) = 0.
	642	ENDDO
	643	c
	644	c periodicite en longitude
	645	c
	646	DO j=1,jjm
	647	fext(iip1,j) = fext(1,j)
	648	ENDDO
	649	DO j=1,jjp1
	650	constang(iip1,j) = constang(1,j)
	651	ENDDO
	652
	653	c fin du changement
	654
	655	c
	656	c-----------------------------------------------------------------------
	657	c
	658	WRITE(6,) ' Coordonnees de la grille * '
	659	WRITE(6,995)
	660	c
	661	WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) '
	662	WRITE(6,995)
	663	DO i=1,iip1
	664	rlonvv(i) = rlonv(i)*180./pi
	665	ENDDO
	666	WRITE(6,400) rlonvv
	667	c
	668	WRITE(6,995)
	669	WRITE(6,*) ' LATITUDES aux pts. V ( degres ) '
	670	WRITE(6,995)
	671	DO i=1,jjm
	672	rlatuu(i)=rlatv(i)*180./pi
	673	ENDDO
	674	WRITE(6,400) (rlatuu(i),i=1,jjm)
	675	c
	676	DO i=1,iip1
	677	rlonvv(i)=rlonu(i)*180./pi
	678	ENDDO
	679	WRITE(6,995)
	680	WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) '
	681	WRITE(6,995)
	682	WRITE(6,400) rlonvv
	683	WRITE(6,995)
	684
	685	WRITE(6,*) ' LATITUDES aux pts. U ( degres ) '
	686	WRITE(6,995)
	687	DO i=1,jjp1
	688	rlatuu(i)=rlatu(i)*180./pi
	689	ENDDO
	690	WRITE(6,400) (rlatuu(i),i=1,jjp1)
	691	WRITE(6,995)
	692	c
	693	444 format(f10.3,f6.0)
	694	400 FORMAT(1x,8f8.2)
	695	990 FORMAT(//)
	696	995 FORMAT(/)
	697	c
	698	RETURN
	699	END

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

Download in other formats:

Original Format