Context Navigation

inigeom.F @ 2218

Last change on this file since 2218 was 2218, checked in by lguez, 9 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 lguez $
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