Context Navigation

inigeom.F @ 2337

Last change on this file since 2337 was 2218, checked in by lguez, 10 years ago

Bug fix in fxhyp. There was some bad logic for the shifting of
longitude grids rlonuv, rlonv, rlonm025 and rlonp025 toward [- pi,
pi]. In some cases, one of the four grids was not shifted and the
others were. For example, you could see the bug with: iim = 48, clon =
20, grossismx = 3, dzoomx = 0.15. The bad logic involved the variable
is2 in the loop on ik = 1, 4. The loop included some tests depending
on ik. The whole thing was quite convoluted. Rewrote fxhyp. In
particular, replaced the loop on ik by a call to a new procedure,
invert_zoom_x. fxhyp.F was in dyn3d, it is replaced by fxhyp_m.F90
which is in dyn3d_common. Removed several arguments of fxhyp: zoom
parameters are accessed through serre.h; rlonm025 and rlonp025 were
not used in inigeom; min and max of longitude steps are written inside
fxhyp.

Some simplifications and modernizations in fyhyp. In particular,
several arguments are removed: zoom parameters are accessed through
serre.h; yprimv was not used in inigeom; min and max of latitude steps
are written inside fyhyp.

Removed now useless intermediary procedure fxyhyper.

Changed default value of dzoomx and dzoomy from 0 to 0.2. dzoom[xy] is
only needed when grossism[xy] > 1 and then it should be > 0.

dzoom[xy] can no longer be the extent of the zoom in degrees: it must
be expressed as the fraction of the total domain.

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.7 KB

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

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

Download in other formats:

Original Format