Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/inigeom.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 29.5 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	IMPLICIT NONE
19	c
20	!-----------------------------------------------------------------------
21	! INCLUDE 'dimensions.h'
22	!
23	! dimensions.h contient les dimensions du modele
24	! ndm est tel que iim=2**ndm
25	!-----------------------------------------------------------------------
26
27	INTEGER iim,jjm,llm,ndm
28
29	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
30
31	!-----------------------------------------------------------------------
32	!
33	! $Header$
34	!
35	!
36	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
37	! veillez n'utiliser que des ! pour les commentaires
38	! et bien positionner les & des lignes de continuation
39	! (les placer en colonne 6 et en colonne 73)
40	!
41	!
42	!-----------------------------------------------------------------------
43	! INCLUDE 'paramet.h'
44
45	INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
46	INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
47	INTEGER ijmllm,mvar
48	INTEGER jcfil,jcfllm
49
50	PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
51	& ,jjp1=jjm+1-1/jjm)
52	PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
53	PARAMETER( kftd = iim/2 -ndm )
54	PARAMETER( ip1jm = iip1jjm, ip1jmp1= iip1jjp1 )
55	PARAMETER( ip1jmi1= ip1jm - iip1 )
56	PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
57	PARAMETER( mvar= ip1jmp1( 2llm+1) + ijmllm )
58	PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
59
60	!-----------------------------------------------------------------------
61	!
62	! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
63	!
64	!-----------------------------------------------------------------------
65	! INCLUDE comconst.h
66
67	COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, &
68	& iflag_top_bound,mode_top_bound
69	COMMON/comconstr/dtvr,daysec, &
70	& pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg &
71	& ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta &
72	& ,dissip_pupstart ,tau_top_bound, &
73	& daylen,molmass, ihf
74	COMMON/cpdetvenus/cpp,nu_venus,t0_venus
75
76	INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
77	REAL dtvr ! dynamical time step (in s)
78	REAL daysec !length (in s) of a standard day
79	REAL pi ! something like 3.14159....
80	REAL dtphys ! (s) time step for the physics
81	REAL dtdiss ! (s) time step for the dissipation
82	REAL rad ! (m) radius of the planet
83	REAL r ! Reduced Gas constant r=R/mu
84	! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
85	REAL cpp ! Cp
86	REAL kappa ! kappa=R/Cp
87	REAL cotot
88	REAL unsim ! = 1./iim
89	REAL g ! (m/s2) gravity
90	REAL omeg ! (rad/s) rotation rate of the planet
91	! Dissipation factors, for Earth model:
92	REAL dissip_factz,dissip_zref !dissip_deltaz
93	! Dissipation factors, for other planets:
94	REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
95	REAL dissip_pupstart
96	INTEGER iflag_top_bound,mode_top_bound
97	REAL tau_top_bound
98	REAL daylen ! length of solar day, in 'standard' day length
99	REAL molmass ! (g/mol) molar mass of the atmosphere
100
101	REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
102	REAL ihf ! (W/m2) intrinsic heat flux for giant planets
103
104
105	!-----------------------------------------------------------------------
106	!
107	! $Header$
108	!
109	!CDK comgeom2
110	COMMON/comgeom/ &
111	& cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm) , &
112	& aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1) , &
113	& airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols , &
114	& unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm) , &
115	& aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1) , &
116	& aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1) , &
117	& alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1) , &
118	& alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1) , &
119	& fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm), &
120	& rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm) , &
121	& cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1) , &
122	& cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
123	& cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
124	& unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2 , &
125	& unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1) , &
126	& unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm) &
127	& , xprimu(iip1),xprimv(iip1)
128
129
130	REAL &
131	& cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
132	& ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4 , &
133	& alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
134	& fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2 , &
135	& cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1 , &
136	& unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2 , &
137	& unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu , &
138	& cusurcvu,xprimu,xprimv
139	!
140	! $Header$
141	!
142	!c
143	!c
144	!c..include serre.h
145	!c
146	REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo, &
147	& grossismx, grossismy, dzoomx, dzoomy,taux,tauy
148	COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo , &
149	& grossismx, grossismy, dzoomx, dzoomy,taux,tauy
150	!
151	! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
152	!
153	!
154	! NB: keep items of different kinds in seperate common blocs to avoid
155	! "misaligned commons" issues
156	!-----------------------------------------------------------------------
157	! INCLUDE 'logic.h'
158
159	COMMON/logicl/ purmats,forward,leapf,apphys, &
160	& statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus &
161	& ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile &
162	& ,ok_limit,ok_etat0,hybrid &
163	& ,moyzon_mu,moyzon_ch
164
165	COMMON/logici/ iflag_phys,iflag_trac
166
167	LOGICAL purmats,forward,leapf,apphys,statcl,conser, &
168	& apdiss,apdelq,saison,ecripar,fxyhypb,ysinus &
169	& ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile &
170	& ,ok_limit,ok_etat0
171	logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
172	! (only used if disvert_type==2)
173	logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
174
175	integer iflag_phys,iflag_trac
176	!$OMP THREADPRIVATE(/logicl/)
177	!$OMP THREADPRIVATE(/logici/)
178	!-----------------------------------------------------------------------
179	!
180	! $Id: comdissnew.h 1319 2010-02-23 21:29:54Z fairhead $
181	!
182	!
183	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
184	! veillez à n'utiliser que des ! pour les commentaires
185	! et à bien positionner les & des lignes de continuation
186	! (les placer en colonne 6 et en colonne 73)
187	!
188	!-----------------------------------------------------------------------
189	! INCLUDE 'comdissnew.h'
190
191	COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv, &
192	& tetagrot,tetatemp,coefdis, vert_prof_dissip
193
194	LOGICAL lstardis
195	INTEGER nitergdiv, nitergrot, niterh
196
197	! For the Earth model:
198	integer vert_prof_dissip ! vertical profile of horizontal dissipation
199	! Allowed values:
200	! 0: rational fraction, function of pressure
201	! 1: tanh of altitude
202
203	REAL tetagdiv, tetagrot, tetatemp, coefdis
204
205	!
206	! ... Les parametres de ce common comdissnew sont lues par defrun_new
207	! sur le fichier run.def ....
208	!
209	!-----------------------------------------------------------------------
210
211	c-----------------------------------------------------------------------
212	c .... Variables locales ....
213	c
214	INTEGER i,j,itmax,itmay,iter
215	REAL cvu(iip1,jjp1),cuv(iip1,jjm)
216	REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm
217	REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm
218	REAL coslatm,coslatp,radclatm,radclatp
219	REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1),
220	* cuij4(iip1,jjp1)
221	REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1),
222	* cvij4(iip1,jjp1)
223	REAL rlonvv(iip1),rlatuu(jjp1)
224	REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) ,
225	* yprimv(jjm),yprimu(jjp1)
226	REAL gamdi_gdiv, gamdi_grot, gamdi_h
227
228	REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),
229	, xprimp025(iip1)
230	SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu
231	SAVE rlonm025,xprimm025,rlonp025,xprimp025
232
233	REAL SSUM
234	c
235	c
236	c ------------------------------------------------------------------
237	c - -
238	c - calcul des coeff. ( cu, cv , 1./cu2, 1./cv2 ) -
239	c - -
240	c ------------------------------------------------------------------
241	c
242	c les coef. ( cu, cv ) permettent de passer des vitesses naturelles
243	c aux vitesses covariantes et contravariantes , ou vice-versa ...
244	c
245	c
246	c on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu
247	c v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv
248	c
249	c on en tire : u(covariant) = cu * cu * u(contravariant)
250	c v(covariant) = cv * cv * v(contravariant)
251	c
252	c
253	c on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
254	c = =
255	c et - jm/2 < Y < jm/2
256	c = =
257	c
258	c ...................................................
259	c ...................................................
260	c . x est la longitude du point en radians .
261	c . y est la latitude du point en radians .
262	c . .
263	c . on a : cu(i,j) = rad * COS(y) * dx/dX .
264	c . cv( j ) = rad * dy/dY .
265	c . aire(i,j) = cu(i,j) * cv(j) .
266	c . .
267	c . y, dx/dX, dy/dY calcules aux points concernes .
268	c . .
269	c ...................................................
270	c ...................................................
271	c
272	c
273	c
274	c ,
275	c cv , bien que dependant de j uniquement,sera ici indice aussi en i
276	c pour un adressage plus facile en ij .
277	c
278	c
279	c
280	c ************ aux points u et v , ***************
281	c xprimu et xprimv sont respectivement les valeurs de dx/dX
282	c yprimu et yprimv . . . . . . . . . . . dy/dY
283	c rlatu et rlatv . . . . . . . . . . .la latitude
284	c cvu et cv . . . . . . . . . . . cv
285	c
286	c ************ aux points u, v, scalaires, et z **************
287	c cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu
288	c
289	c
290	c
291	c Exemple de distribution de variables sur la grille dans le
292	c domaine de travail ( X,Y ) .
293	c ................................................................
294	c DX=DY= 1
295	c
296	c
297	c + represente un point scalaire ( p.exp la pression )
298	c > represente la composante zonale du vent
299	c V represente la composante meridienne du vent
300	c o represente la vorticite
301	c
302	c ---- , car aux poles , les comp.zonales covariantes sont nulles
303	c
304	c
305	c
306	c i ->
307	c
308	c 1 2 3 4 5 6 7 8
309	c j
310	c v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
311	c
312	c V o V o V o V o V o V o V o V o
313	c
314	c 2 + > + > + > + > + > + > + > + >
315	c
316	c V o V o V o V o V o V o V o V o
317	c
318	c 3 + > + > + > + > + > + > + > + >
319	c
320	c V o V o V o V o V o V o V o V o
321	c
322	c 4 + > + > + > + > + > + > + > + >
323	c
324	c V o V o V o V o V o V o V o V o
325	c
326	c 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
327	c
328	c
329	c Ci-dessus, on voit que le nombre de pts.en longitude est egal
330	c a IM = 8
331	c De meme , le nombre d'intervalles entre les 2 poles est egal
332	c a JM = 4
333	c
334	c Les points scalaires ( + ) correspondent donc a des valeurs
335	c entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
336	c
337	c Les vents U ( > ) correspondent a des valeurs semi-
338	c entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
339	c
340	c Les vents V ( V ) correspondent a des valeurs entieres
341	c de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
342	c
343	c
344	c
345	WRITE(6,3)
346	3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ',
347	* //5x,' Calcul des elongations cu et cv comme sommes des 4 ' /
348	* 5x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux
349	* '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / )
350	c
351	c
352	IF( nitergdiv.NE.2 ) THEN
353	gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. )
354	ELSE
355	gamdi_gdiv = 0.
356	ENDIF
357	IF( nitergrot.NE.2 ) THEN
358	gamdi_grot = coefdis/ ( REAL(nitergrot) -2. )
359	ELSE
360	gamdi_grot = 0.
361	ENDIF
362	IF( niterh.NE.2 ) THEN
363	gamdi_h = coefdis/ ( REAL(niterh) -2. )
364	ELSE
365	gamdi_h = 0.
366	ENDIF
367
368	WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
369	* nitergdiv,nitergrot,niterh
370	c
371	pi = 2.* ASIN(1.)
372	c
373	WRITE(6,990)
374
375	c ----------------------------------------------------------------
376	c
377	IF( .NOT.fxyhypb ) THEN
378	c
379	c
380	IF( ysinus ) THEN
381	c
382	WRITE(6,) ' Inigeom , Y = Sinus ( Latitude ) * '
383	c
384	c .... utilisation de f(x,y ) avec y = sinus de la latitude .....
385
386	CALL fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
387	, rlatu2,yprimu2,
388	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
389
390	ELSE
391	c
392	WRITE(6,) ' Inigeom , Y = Latitude , der. sinusoid . *'
393
394	c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
395	c
396
397	pxo = clon *pi /180.
398	pyo = 2.* clat* pi /180.
399	c
400	c .... determination de transx ( pour le zoom ) par Newton-Raphson ...
401	c
402	itmax = 10
403	eps = .1e-7
404	c
405	xo1 = 0.
406	DO 10 iter = 1, itmax
407	x1 = xo1
408	f = x1+ alphax *SIN(x1-pxo)
409	df = 1.+ alphax *COS(x1-pxo)
410	x1 = x1 - f/df
411	xdm = ABS( x1- xo1 )
412	IF( xdm.LE.eps )GO TO 11
413	xo1 = x1
414	10 CONTINUE
415	11 CONTINUE
416	c
417	transx = xo1
418
419	itmay = 10
420	eps = .1e-7
421	C
422	yo1 = 0.
423	DO 15 iter = 1,itmay
424	y1 = yo1
425	f = y1 + alphay* SIN(y1-pyo)
426	df = 1. + alphay* COS(y1-pyo)
427	y1 = y1 -f/df
428	ydm = ABS(y1-yo1)
429	IF(ydm.LE.eps) GO TO 17
430	yo1 = y1
431	15 CONTINUE
432	c
433	17 CONTINUE
434	transy = yo1
435
436	CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1,
437	, rlatu2,yprimu2,
438	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
439
440	ENDIF
441	c
442	ELSE
443	c
444	c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
445	c .....................................................................
446
447	WRITE(6,)' Inigeom , Y = Latitude , der.tg. hyperbolique *'
448
449	CALL fxyhyper( clat, grossismy, dzoomy, tauy ,
450	, clon, grossismx, dzoomx, taux ,
451	, rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 ,
452	, rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
453
454
455	ENDIF
456	c
457	c -------------------------------------------------------------------
458
459	c
460	rlatu(1) = ASIN(1.)
461	rlatu(jjp1) = - rlatu(1)
462	c
463	c
464	c .... calcul aux poles ....
465	c
466	yprimu(1) = 0.
467	yprimu(jjp1) = 0.
468	c
469	c
470	un4rad2 = 0.25 * rad * rad
471	c
472	c --------------------------------------------------------------------
473	c --------------------------------------------------------------------
474	c - -
475	c - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) -
476	c - et de fext , force de coriolis extensive . -
477	c - -
478	c --------------------------------------------------------------------
479	c --------------------------------------------------------------------
480	c
481	c
482	c
483	c A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
484	c affectees 4 aires entourant P , calculees respectivement aux points
485	c ( i + 1/4, j - 1/4 ) : aireij1 (i,j)
486	c ( i + 1/4, j + 1/4 ) : aireij2 (i,j)
487	c ( i - 1/4, j + 1/4 ) : aireij3 (i,j)
488	c ( i - 1/4, j - 1/4 ) : aireij4 (i,j)
489	c
490	c ,
491	c Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
492	c Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
493	c des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au
494	c point (i,j) .
495	c On definit en outre les coefficients alpha comme etant egaux a
496	c (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j)
497	c
498	c De meme, toute aire centree en 1 point U est egale a la somme des
499	c 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
500	c Idem pour airev, airez .
501	c
502	c On a ,pour chaque maille : dX = dY = 1
503	c
504	c
505	c . V
506	c
507	c aireij4 . . aireij1
508	c
509	c U . . P . U
510	c
511	c aireij3 . . aireij2
512	c
513	c . V
514	c
515	c
516	c
517	c
518	c
519	c ....................................................................
520	c
521	c Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
522	c qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
523	c taires cuij et les 4 elongat. cvij qui sont calculees aux memes
524	c endroits que les aireij .
525	c
526	c ....................................................................
527	c
528	c ....... do 35 : boucle sur les jjm + 1 latitudes .....
529	c
530	c
531	DO 35 j = 1, jjp1
532	c
533	IF ( j. eq. 1 ) THEN
534	c
535	yprm = yprimu1(j)
536	rlatm = rlatu1(j)
537	c
538	coslatm = COS( rlatm )
539	radclatm = 0.5* rad * coslatm
540	c
541	DO 30 i = 1, iim
542	xprp = xprimp025( i )
543	xprm = xprimm025( i )
544	aireij2( i,1 ) = un4rad2 * coslatm * xprp * yprm
545	aireij3( i,1 ) = un4rad2 * coslatm * xprm * yprm
546	cuij2 ( i,1 ) = radclatm * xprp
547	cuij3 ( i,1 ) = radclatm * xprm
548	cvij2 ( i,1 ) = 0.5* rad * yprm
549	cvij3 ( i,1 ) = cvij2(i,1)
550	30 CONTINUE
551	c
552	DO i = 1, iim
553	aireij1( i,1 ) = 0.
554	aireij4( i,1 ) = 0.
555	cuij1 ( i,1 ) = 0.
556	cuij4 ( i,1 ) = 0.
557	cvij1 ( i,1 ) = 0.
558	cvij4 ( i,1 ) = 0.
559	ENDDO
560	c
561	END IF
562	c
563	IF ( j. eq. jjp1 ) THEN
564	yprp = yprimu2(j-1)
565	rlatp = rlatu2 (j-1)
566	ccc yprp = fyprim( REAL(j) - 0.25 )
567	ccc rlatp = fy ( REAL(j) - 0.25 )
568	c
569	coslatp = COS( rlatp )
570	radclatp = 0.5* rad * coslatp
571	c
572	DO 31 i = 1,iim
573	xprp = xprimp025( i )
574	xprm = xprimm025( i )
575	aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp
576	aireij4( i,jjp1 ) = un4rad2 * coslatp * xprm * yprp
577	cuij1(i,jjp1) = radclatp * xprp
578	cuij4(i,jjp1) = radclatp * xprm
579	cvij1(i,jjp1) = 0.5 * rad* yprp
580	cvij4(i,jjp1) = cvij1(i,jjp1)
581	31 CONTINUE
582	c
583	DO i = 1, iim
584	aireij2( i,jjp1 ) = 0.
585	aireij3( i,jjp1 ) = 0.
586	cvij2 ( i,jjp1 ) = 0.
587	cvij3 ( i,jjp1 ) = 0.
588	cuij2 ( i,jjp1 ) = 0.
589	cuij3 ( i,jjp1 ) = 0.
590	ENDDO
591	c
592	END IF
593	c
594
595	IF ( j .gt. 1 .AND. j .lt. jjp1 ) THEN
596	c
597	rlatp = rlatu2 ( j-1 )
598	yprp = yprimu2( j-1 )
599	rlatm = rlatu1 ( j )
600	yprm = yprimu1( j )
601	cc rlatp = fy ( REAL(j) - 0.25 )
602	cc yprp = fyprim( REAL(j) - 0.25 )
603	cc rlatm = fy ( REAL(j) + 0.25 )
604	cc yprm = fyprim( REAL(j) + 0.25 )
605
606	coslatm = COS( rlatm )
607	coslatp = COS( rlatp )
608	radclatp = 0.5* rad * coslatp
609	radclatm = 0.5* rad * coslatm
610	c
611	ai14 = un4rad2 * coslatp * yprp
612	ai23 = un4rad2 * coslatm * yprm
613	DO 32 i = 1,iim
614	xprp = xprimp025( i )
615	xprm = xprimm025( i )
616
617	aireij1 ( i,j ) = ai14 * xprp
618	aireij2 ( i,j ) = ai23 * xprp
619	aireij3 ( i,j ) = ai23 * xprm
620	aireij4 ( i,j ) = ai14 * xprm
621	cuij1 ( i,j ) = radclatp * xprp
622	cuij2 ( i,j ) = radclatm * xprp
623	cuij3 ( i,j ) = radclatm * xprm
624	cuij4 ( i,j ) = radclatp * xprm
625	cvij1 ( i,j ) = 0.5* rad * yprp
626	cvij2 ( i,j ) = 0.5* rad * yprm
627	cvij3 ( i,j ) = cvij2(i,j)
628	cvij4 ( i,j ) = cvij1(i,j)
629	32 CONTINUE
630	c
631	END IF
632	c
633	c ........ periodicite ............
634	c
635	cvij1 (iip1,j) = cvij1 (1,j)
636	cvij2 (iip1,j) = cvij2 (1,j)
637	cvij3 (iip1,j) = cvij3 (1,j)
638	cvij4 (iip1,j) = cvij4 (1,j)
639	cuij1 (iip1,j) = cuij1 (1,j)
640	cuij2 (iip1,j) = cuij2 (1,j)
641	cuij3 (iip1,j) = cuij3 (1,j)
642	cuij4 (iip1,j) = cuij4 (1,j)
643	aireij1 (iip1,j) = aireij1 (1,j )
644	aireij2 (iip1,j) = aireij2 (1,j )
645	aireij3 (iip1,j) = aireij3 (1,j )
646	aireij4 (iip1,j) = aireij4 (1,j )
647
648	35 CONTINUE
649	c
650	c ..............................................................
651	c
652	DO 37 j = 1, jjp1
653	DO 36 i = 1, iim
654	aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
655	* aireij4(i,j)
656	alpha1 ( i,j ) = aireij1(i,j) / aire(i,j)
657	alpha2 ( i,j ) = aireij2(i,j) / aire(i,j)
658	alpha3 ( i,j ) = aireij3(i,j) / aire(i,j)
659	alpha4 ( i,j ) = aireij4(i,j) / aire(i,j)
660	alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j)
661	alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j)
662	alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j)
663	alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j)
664	36 CONTINUE
665	c
666	c
667	aire (iip1,j) = aire (1,j)
668	alpha1 (iip1,j) = alpha1 (1,j)
669	alpha2 (iip1,j) = alpha2 (1,j)
670	alpha3 (iip1,j) = alpha3 (1,j)
671	alpha4 (iip1,j) = alpha4 (1,j)
672	alpha1p2(iip1,j) = alpha1p2(1,j)
673	alpha1p4(iip1,j) = alpha1p4(1,j)
674	alpha2p3(iip1,j) = alpha2p3(1,j)
675	alpha3p4(iip1,j) = alpha3p4(1,j)
676	37 CONTINUE
677	c
678
679	DO 42 j = 1,jjp1
680	DO 41 i = 1,iim
681	aireu (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
682	* aireij3(i+1,j)
683	unsaire ( i,j)= 1./ aire(i,j)
684	unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv )
685	unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h )
686	airesurg ( i,j)= aire(i,j)/ g
687	41 CONTINUE
688	aireu (iip1,j) = aireu (1,j)
689	unsaire (iip1,j) = unsaire(1,j)
690	unsair_gam1(iip1,j) = unsair_gam1(1,j)
691	unsair_gam2(iip1,j) = unsair_gam2(1,j)
692	airesurg (iip1,j) = airesurg(1,j)
693	42 CONTINUE
694	c
695	c
696	DO 48 j = 1,jjm
697	c
698	DO i=1,iim
699	airev (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
700	* aireij4(i,j+1)
701	ENDDO
702	DO i=1,iim
703	airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) +
704	* aireij4(i+1,j+1)
705	unsairez(i,j) = 1./ airez
706	unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot )
707	fext (i,j) = airez * SIN(rlatv(j))* 2.* omeg
708	ENDDO
709	airev (iip1,j) = airev(1,j)
710	unsairez (iip1,j) = unsairez(1,j)
711	fext (iip1,j) = fext(1,j)
712	unsairz_gam(iip1,j) = unsairz_gam(1,j)
713	c
714	48 CONTINUE
715	c
716	c
717	c ..... Calcul des elongations cu,cv, cvu .........
718	c
719	DO j = 1, jjm
720	DO i = 1, iim
721	cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1))
722	cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) )
723	cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1))
724	unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) )
725	ENDDO
726	DO i = 1, iim
727	cuvsurcv (i,j) = airev(i,j) * unscv2(i,j)
728	cvsurcuv (i,j) = 1./cuvsurcv(i,j)
729	cuvscvgam1(i,j) = cuvsurcv (i,j) ** ( - gamdi_gdiv )
730	cuvscvgam2(i,j) = cuvsurcv (i,j) ** ( - gamdi_h )
731	cvscuvgam(i,j) = cvsurcuv (i,j) ** ( - gamdi_grot )
732	ENDDO
733	cv (iip1,j) = cv (1,j)
734	cvu (iip1,j) = cvu (1,j)
735	unscv2 (iip1,j) = unscv2 (1,j)
736	cuv (iip1,j) = cuv (1,j)
737	cuvsurcv (iip1,j) = cuvsurcv (1,j)
738	cvsurcuv (iip1,j) = cvsurcuv (1,j)
739	cuvscvgam1(iip1,j) = cuvscvgam1(1,j)
740	cuvscvgam2(iip1,j) = cuvscvgam2(1,j)
741	cvscuvgam(iip1,j) = cvscuvgam(1,j)
742	ENDDO
743
744	DO j = 2, jjm
745	DO i = 1, iim
746	cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j))
747	unscu2 (i,j) = 1./ ( cu(i,j) * cu(i,j) )
748	cvusurcu (i,j) = aireu(i,j) * unscu2(i,j)
749	cusurcvu (i,j) = 1./ cvusurcu(i,j)
750	cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv )
751	cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h )
752	cuscvugam (i,j) = cusurcvu(i,j) ** ( - gamdi_grot )
753	ENDDO
754	cu (iip1,j) = cu(1,j)
755	unscu2 (iip1,j) = unscu2(1,j)
756	cvusurcu (iip1,j) = cvusurcu(1,j)
757	cusurcvu (iip1,j) = cusurcvu(1,j)
758	cvuscugam1(iip1,j) = cvuscugam1(1,j)
759	cvuscugam2(iip1,j) = cvuscugam2(1,j)
760	cuscvugam (iip1,j) = cuscvugam(1,j)
761	ENDDO
762
763	c
764	c .... calcul aux poles ....
765	c
766	DO i = 1, iip1
767	cu ( i, 1 ) = 0.
768	unscu2( i, 1 ) = 0.
769	cvu ( i, 1 ) = 0.
770	c
771	cu (i, jjp1) = 0.
772	unscu2(i, jjp1) = 0.
773	cvu (i, jjp1) = 0.
774	ENDDO
775	c
776	c ..............................................................
777	c
778	DO j = 1, jjm
779	DO i= 1, iim
780	airvscu2 (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) )
781	aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot )
782	ENDDO
783	airvscu2 (iip1,j) = airvscu2(1,j)
784	aivscu2gam(iip1,j) = aivscu2gam(1,j)
785	ENDDO
786
787	DO j=2,jjm
788	DO i=1,iim
789	airuscv2 (i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) )
790	aiuscv2gam (i,j) = airuscv2(i,j)** ( - gamdi_grot )
791	ENDDO
792	airuscv2 (iip1,j) = airuscv2 (1,j)
793	aiuscv2gam(iip1,j) = aiuscv2gam(1,j)
794	ENDDO
795
796	c
797	c calcul des aires aux poles :
798	c -----------------------------
799	c
800	apoln = SSUM(iim,aire(1,1),1)
801	apols = SSUM(iim,aire(1,jjp1),1)
802	unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
803	unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
804	unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) )
805	unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) )
806	c
807	c-----------------------------------------------------------------------
808	c gtitre='Coriolis version ancienne'
809	c gfichier='fext1'
810	c CALL writestd(fext,iip1*jjm)
811	c
812	c changement F. Hourdin calcul conservatif pour fext
813	c constang contient le produit a * cos ( latitude ) * omega
814	c
815	DO i=1,iim
816	constang(i,1) = 0.
817	ENDDO
818	DO j=1,jjm-1
819	DO i=1,iim
820	constang(i,j+1) = radomegcu(i,j+1)*COS(rlatu(j+1))
821	ENDDO
822	ENDDO
823	DO i=1,iim
824	constang(i,jjp1) = 0.
825	ENDDO
826	c
827	c periodicite en longitude
828	c
829	DO j=1,jjm
830	fext(iip1,j) = fext(1,j)
831	ENDDO
832	DO j=1,jjp1
833	constang(iip1,j) = constang(1,j)
834	ENDDO
835
836	c fin du changement
837
838	c
839	c-----------------------------------------------------------------------
840	c
841	WRITE(6,) ' Coordonnees de la grille * '
842	WRITE(6,995)
843	c
844	WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) '
845	WRITE(6,995)
846	DO i=1,iip1
847	rlonvv(i) = rlonv(i)*180./pi
848	ENDDO
849	WRITE(6,400) rlonvv
850	c
851	WRITE(6,995)
852	WRITE(6,*) ' LATITUDES aux pts. V ( degres ) '
853	WRITE(6,995)
854	DO i=1,jjm
855	rlatuu(i)=rlatv(i)*180./pi
856	ENDDO
857	WRITE(6,400) (rlatuu(i),i=1,jjm)
858	c
859	DO i=1,iip1
860	rlonvv(i)=rlonu(i)*180./pi
861	ENDDO
862	WRITE(6,995)
863	WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) '
864	WRITE(6,995)
865	WRITE(6,400) rlonvv
866	WRITE(6,995)
867
868	WRITE(6,*) ' LATITUDES aux pts. U ( degres ) '
869	WRITE(6,995)
870	DO i=1,jjp1
871	rlatuu(i)=rlatu(i)*180./pi
872	ENDDO
873	WRITE(6,400) (rlatuu(i),i=1,jjp1)
874	WRITE(6,995)
875	c
876	444 format(f10.3,f6.0)
877	400 FORMAT(1x,8f8.2)
878	990 FORMAT(//)
879	995 FORMAT(/)
880	c
881	RETURN
882	END

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

Download in other formats: