Context Navigation

inigeom.F @ 1632

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

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 21.8 KB

Line
1	!
2	! $Id: inigeom.F 1299 2010-01-20 14:27:21Z 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	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