Context Navigation

inigeom.F @ 3751

Last change on this file since 3751 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 21.7 KB

Line
1	!
2	! $Id: inigeom.F 1403 2010-07-01 09:02:53Z adurocher $
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