Context Navigation

inigeom.F @ 2757

Last change on this file since 2757 was 1441, checked in by emillour, 9 years ago

Updates in common dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2250):

compilation:
added test in grid/dimension/makdim to check that # of longitudes is a multiple of 8

dyn3d_common:

Bug correction concerning zoom (cf LMDZ5 rev 2218)

coefpoly.F becomes coefpoly_m.F90 (in misc)
fxhyp.F => fxhyp_m.F90 , fyhyp.F => fyhyp_m.F90
new routines for zoom: invert_zoom_x_m.F90 and principal_cshift_m.F90
inigeom.F adapted to new zoom definition routines
fluxstokenc.F : got rid of calls to initial0()

dyn3d:
advtrac.F90 : got rid of calls to initial0()
conf_gcm.F90 : cosmetic changes and change in default dzoomx,dzoomy values
guide_mod.F90 : followed updates from Earth Model
gcm.F is now gcm.F90

dyn3dpar:
advtrac_p.F90, covcont_p.F90, mod_hallo.F90 : cosmetic changes
conf_gcm.F90 : cosmetic and changed in default dzoomx,dzoomy values
parallel_lmdz.F90 : updates to keep up with Earth model

misc:
arth.F90 becomes arth_m.F90
wxios.F90 updated wrt Earth model changes
nrtype.F90 and coefpoly_m.F90 added
ran1.F, sort.F, minmax.F, minmax2.F, juldate.F moved over from dyn3d_common

File size: 21.9 KB

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

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

Download in other formats:

Original Format