Context Navigation

inigeom.F @ 3717

Last change on this file since 3717 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

Line
1	!
2	! $Id: inigeom.F 2218 2015-03-02 16:11:07Z 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	use fxhyp_m, only: fxhyp
19	use fyhyp_m, only: fyhyp
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
172	gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
173	ELSE
174	gamdi_gdiv = 0.
175	ENDIF
176	IF( nitergrot.NE.2 ) THEN
177	gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
178	ELSE
179	gamdi_grot = 0.
180	ENDIF
181	IF( niterh.NE.2 ) THEN
182	gamdi_h = coefdis/ ( REAL(niterh) -2. )
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
268	CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
269	CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
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)
382	ccc yprp = fyprim( REAL(j) - 0.25 )
383	ccc rlatp = fy ( REAL(j) - 0.25 )
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 )
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 )
421
422	coslatm = COS( rlatm )
423	coslatp = COS( rlatp )
424	radclatp = 0.5* rad * coslatp
425	radclatm = 0.5* rad * coslatm
426	c
427	ai14 = un4rad2 * coslatp * yprp
428	ai23 = un4rad2 * coslatm * yprm
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