Context Navigation

inigeom.F @ 5425

Last change on this file since 5425 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.9 KB

Rev	Line
[524]	1	!
[1403]	2	! $Id: inigeom.F 2641 2016-09-29 21:26:46Z jyg $
[524]	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
[2258]	18	use fxhyp_m, only: fxhyp
	19	use fyhyp_m, only: fyhyp
[2641]	20	USE comconst_mod, ONLY: pi, g, omeg, rad
	21	USE logic_mod, ONLY: fxyhypb, ysinus
	22	USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy,
	23	& alphax,alphay,taux,tauy,transx,transy,pxo,pyo
[524]	24	IMPLICIT NONE
	25	c
[2641]	26	include "dimensions.h"
	27	include "paramet.h"
	28	include "comgeom2.h"
	29	include "comdissnew.h"
[524]	30
	31	c-----------------------------------------------------------------------
	32	c .... Variables locales ....
	33	c
	34	INTEGER i,j,itmax,itmay,iter
	35	REAL cvu(iip1,jjp1),cuv(iip1,jjm)
	36	REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
	37	REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
	38	REAL coslatm,coslatp,radclatm,radclatp
	39	REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
	40	* cuij4(iip1,jjp1)
	41	REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
	42	* cvij4(iip1,jjp1)
	43	REAL rlonvv(iip1),rlatuu(jjp1)
	44	REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
	45	* yprimv(jjm),yprimu(jjp1)
	46	REAL gamdi_gdiv, gamdi_grot, gamdi_h
	47
	48	REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
	49	, xprimp025(iip1)
	50	SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
	51	SAVE rlonm025,xprimm025,rlonp025,xprimp025
	52
	53	REAL SSUM
	54	c
	55	c
	56	c ------------------------------------------------------------------
	57	c - -
	58	c - calcul des coeff. ( cu, cv , 1./cu2, 1./cv2 ) -
	59	c - -
	60	c ------------------------------------------------------------------
	61	c
	62	c les coef. ( cu, cv ) permettent de passer des vitesses naturelles
	63	c aux vitesses covariantes et contravariantes , ou vice-versa ...
	64	c
	65	c
	66	c on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu
	67	c v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv
	68	c
	69	c on en tire : u(covariant) = cu * cu * u(contravariant)
	70	c v(covariant) = cv * cv * v(contravariant)
	71	c
	72	c
	73	c on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
	74	c = =
	75	c et - jm/2 < Y < jm/2
	76	c = =
	77	c
	78	c ...................................................
	79	c ...................................................
	80	c . x est la longitude du point en radians .
	81	c . y est la latitude du point en radians .
	82	c . .
	83	c . on a : cu(i,j) = rad * COS(y) * dx/dX .
	84	c . cv( j ) = rad * dy/dY .
	85	c . aire(i,j) = cu(i,j) * cv(j) .
	86	c . .
	87	c . y, dx/dX, dy/dY calcules aux points concernes .
	88	c . .
	89	c ...................................................
	90	c ...................................................
	91	c
	92	c
	93	c
	94	c ,
	95	c cv , bien que dependant de j uniquement,sera ici indice aussi en i
	96	c pour un adressage plus facile en ij .
	97	c
	98	c
	99	c
	100	c ************ aux points u et v , ***************
	101	c xprimu et xprimv sont respectivement les valeurs de dx/dX
	102	c yprimu et yprimv . . . . . . . . . . . dy/dY
	103	c rlatu et rlatv . . . . . . . . . . .la latitude
	104	c cvu et cv . . . . . . . . . . . cv
	105	c
	106	c ************ aux points u, v, scalaires, et z **************
	107	c cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu
	108	c
	109	c
	110	c
	111	c Exemple de distribution de variables sur la grille dans le
	112	c domaine de travail ( X,Y ) .
	113	c ................................................................
	114	c DX=DY= 1
	115	c
	116	c
	117	c + represente un point scalaire ( p.exp la pression )
	118	c > represente la composante zonale du vent
	119	c V represente la composante meridienne du vent
	120	c o represente la vorticite
	121	c
	122	c ---- , car aux poles , les comp.zonales covariantes sont nulles
	123	c
	124	c
	125	c
	126	c i ->
	127	c
	128	c 1 2 3 4 5 6 7 8
	129	c j
	130	c v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
	131	c
	132	c V o V o V o V o V o V o V o V o
	133	c
	134	c 2 + > + > + > + > + > + > + > + >
	135	c
	136	c V o V o V o V o V o V o V o V o
	137	c
	138	c 3 + > + > + > + > + > + > + > + >
	139	c
	140	c V o V o V o V o V o V o V o V o
	141	c
	142	c 4 + > + > + > + > + > + > + > + >
	143	c
	144	c V o V o V o V o V o V o V o V o
	145	c
	146	c 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
	147	c
	148	c
	149	c Ci-dessus, on voit que le nombre de pts.en longitude est egal
	150	c a IM = 8
	151	c De meme , le nombre d'intervalles entre les 2 poles est egal
	152	c a JM = 4
	153	c
	154	c Les points scalaires ( + ) correspondent donc a des valeurs
	155	c entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
	156	c
	157	c Les vents U ( > ) correspondent a des valeurs semi-
	158	c entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
	159	c
	160	c Les vents V ( V ) correspondent a des valeurs entieres
	161	c de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
	162	c
	163	c
	164	c
	165	WRITE(6,3)
	166	3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ',
	167	* //5x,' Calcul des elongations cu et cv comme sommes des 4 ' /
	168	* 5x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux
	169	* '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / )
	170	c
	171	c
	172	IF( nitergdiv.NE.2 ) THEN
[1403]	173	gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
[524]	174	ELSE
	175	gamdi_gdiv = 0.
	176	ENDIF
	177	IF( nitergrot.NE.2 ) THEN
[1403]	178	gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
[524]	179	ELSE
	180	gamdi_grot = 0.
	181	ENDIF
	182	IF( niterh.NE.2 ) THEN
[1403]	183	gamdi_h = coefdis/ ( REAL(niterh) -2. )
[524]	184	ELSE
	185	gamdi_h = 0.
	186	ENDIF
	187
	188	WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
	189	* nitergdiv,nitergrot,niterh
	190	c
	191	pi = 2.* ASIN(1.)
	192	c
	193	WRITE(6,990)
	194
	195	c ----------------------------------------------------------------
	196	c
	197	IF( .NOT.fxyhypb ) THEN
	198	c
	199	c
	200	IF( ysinus ) THEN
	201	c
	202	WRITE(6,) ' Inigeom , Y = Sinus ( Latitude ) * '
	203	c
	204	c .... utilisation de f(x,y ) avec y = sinus de la latitude .....
	205
	206	CALL fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
	207	, rlatu2,yprimu2,
	208	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
	209
	210	ELSE
	211	c
	212	WRITE(6,) ' Inigeom , Y = Latitude , der. sinusoid . *'
	213
	214	c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
	215	c
	216
	217	pxo = clon *pi /180.
	218	pyo = 2.* clat* pi /180.
	219	c
	220	c .... determination de transx ( pour le zoom ) par Newton-Raphson ...
	221	c
	222	itmax = 10
	223	eps = .1e-7
	224	c
	225	xo1 = 0.
	226	DO 10 iter = 1, itmax
	227	x1 = xo1
	228	f = x1+ alphax *SIN(x1-pxo)
	229	df = 1.+ alphax *COS(x1-pxo)
	230	x1 = x1 - f/df
	231	xdm = ABS( x1- xo1 )
	232	IF( xdm.LE.eps )GO TO 11
	233	xo1 = x1
	234	10 CONTINUE
	235	11 CONTINUE
	236	c
	237	transx = xo1
	238
	239	itmay = 10
	240	eps = .1e-7
	241	C
	242	yo1 = 0.
	243	DO 15 iter = 1,itmay
	244	y1 = yo1
	245	f = y1 + alphay* SIN(y1-pyo)
	246	df = 1. + alphay* COS(y1-pyo)
	247	y1 = y1 -f/df
	248	ydm = ABS(y1-yo1)
	249	IF(ydm.LE.eps) GO TO 17
	250	yo1 = y1
	251	15 CONTINUE
	252	c
	253	17 CONTINUE
	254	transy = yo1
	255
	256	CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
	257	, rlatu2,yprimu2,
	258	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
	259
	260	ENDIF
	261	c
	262	ELSE
	263	c
	264	c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
	265	c .....................................................................
	266
	267	WRITE(6,)' Inigeom , Y = Latitude , der.tg. hyperbolique *'
	268
[2258]	269	CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
	270	CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
[524]	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)
[1403]	383	ccc yprp = fyprim( REAL(j) - 0.25 )
	384	ccc rlatp = fy ( REAL(j) - 0.25 )
[524]	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 )
[1403]	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 )
[524]	422
	423	coslatm = COS( rlatm )
	424	coslatp = COS( rlatp )
	425	radclatp = 0.5* rad * coslatp
	426	radclatm = 0.5* rad * coslatm
	427	c
[1783]	428	ai14 = un4rad2 * coslatp * yprp
	429	ai23 = un4rad2 * coslatm * yprm
[524]	430	DO 32 i = 1,iim
	431	xprp = xprimp025( i )
	432	xprm = xprimm025( i )
	433
	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