Context Navigation

calfis.F @ 97

Last change on this file since 97 was 97, checked in by slebonnois, 14 years ago
Serie de modifs SL pour homogeneisation des phytitan et phyvenus Ca touche aussi aux liens phy/dyn (surtout a propos de clesphy0), a verifier avec les autres, donc...
File size: 20.0 KB

Line
1	!
2	! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
3	!
4	C
5	C
6	SUBROUTINE calfis(lafin,
7	$ jD_cur, jH_cur,
8	$ pucov,
9	$ pvcov,
10	$ pteta,
11	$ pq,
12	$ pmasse,
13	$ pps,
14	$ pp,
15	$ ppk,
16	$ pphis,
17	$ pphi,
18	$ pducov,
19	$ pdvcov,
20	$ pdteta,
21	$ pdq,
22	$ flxw,
23	$ pdufi,
24	$ pdvfi,
25	$ pdhfi,
26	$ pdqfi,
27	$ pdpsfi)
28	c
29	c Auteur : P. Le Van, F. Hourdin
30	c .........
31	USE infotrac
32	USE control_mod
33
34
35	IMPLICIT NONE
36	c=======================================================================
37	c
38	c 1. rearrangement des tableaux et transformation
39	c variables dynamiques > variables physiques
40	c 2. calcul des termes physiques
41	c 3. retransformation des tendances physiques en tendances dynamiques
42	c
43	c remarques:
44	c ----------
45	c
46	c - les vents sont donnes dans la physique par leurs composantes
47	c naturelles.
48	c - la variable thermodynamique de la physique est une variable
49	c intensive : T
50	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
51	c - les deux seules variables dependant de la geometrie necessaires
52	c pour la physique sont la latitude pour le rayonnement et
53	c l'aire de la maille quand on veut integrer une grandeur
54	c horizontalement.
55	c - les points de la physique sont les points scalaires de la
56	c la dynamique; numerotation:
57	c 1 pour le pole nord
58	c (jjm-1)*iim pour l'interieur du domaine
59	c ngridmx pour le pole sud
60	c ---> ngridmx=2+(jjm-1)*iim
61	c
62	c Input :
63	c -------
64	c pucov covariant zonal velocity
65	c pvcov covariant meridional velocity
66	c pteta potential temperature
67	c pps surface pressure
68	c pmasse masse d'air dans chaque maille
69	c pts surface temperature (K)
70	c callrad clef d'appel au rayonnement
71	c
72	c Output :
73	c --------
74	c pdufi tendency for the natural zonal velocity (ms-1)
75	c pdvfi tendency for the natural meridional velocity
76	c pdhfi tendency for the potential temperature (K/s)
77	c pdtsfi tendency for the surface temperature
78	c
79	c pdtrad radiative tendencies \ both input
80	c pfluxrad radiative fluxes / and output
81	c
82	c=======================================================================
83	c
84	c-----------------------------------------------------------------------
85	c
86	c 0. Declarations :
87	c ------------------
88
89	#include "dimensions.h"
90	#include "paramet.h"
91	#include "temps.h"
92
93	INTEGER ngridmx
94	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
95
96	#include "comconst.h"
97	#include "comvert.h"
98	#include "comgeom2.h"
99	#include "iniprint.h"
100
101	c Arguments :
102	c -----------
103	LOGICAL lafin
104
105
106	REAL pvcov(iip1,jjm,llm)
107	REAL pucov(iip1,jjp1,llm)
108	REAL pteta(iip1,jjp1,llm)
109	REAL pmasse(iip1,jjp1,llm)
110	REAL pq(iip1,jjp1,llm,nqtot)
111	REAL pphis(iip1,jjp1)
112	REAL pphi(iip1,jjp1,llm)
113	c
114	REAL pdvcov(iip1,jjm,llm)
115	REAL pducov(iip1,jjp1,llm)
116	REAL pdteta(iip1,jjp1,llm)
117	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
118	! qui lui meme ne sert a rien dans la routine telle qu'elle est
119	! ecrite, et que j'ai donc commente....
120	REAL pdq(iip1,jjp1,llm,nqtot)
121	c
122	REAL pps(iip1,jjp1)
123	REAL pp(iip1,jjp1,llmp1)
124	REAL ppk(iip1,jjp1,llm)
125	c
126	c TENDENCIES in */s
127	REAL pdvfi(iip1,jjm,llm)
128	REAL pdufi(iip1,jjp1,llm)
129	REAL pdhfi(iip1,jjp1,llm)
130	REAL pdqfi(iip1,jjp1,llm,nqtot)
131	REAL pdpsfi(iip1,jjp1)
132
133
134	c Local variables :
135	c -----------------
136
137	INTEGER i,j,l,ig0,ig,iq,iiq
138	REAL zpsrf(ngridmx)
139	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
140	REAL zphi(ngridmx,llm),zphis(ngridmx)
141	c
142	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
143	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
144	! ADAPTATION GCM POUR CP(T)
145	REAL zteta(ngridmx,llm)
146	REAL zpk(ngridmx,llm)
147	c
148	! RQ SL 13/10/10:
149	! Ces calculs ne servent pas.
150	! Si necessaire, decommenter ces variables et les calculs...
151	! REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
152	! REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
153	c
154	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
155	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
156	REAL zdpsrf(ngridmx)
157	c
158	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
159	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
160	REAL jH_cur_split,zdt_split
161	LOGICAL debut_split,lafin_split
162	INTEGER isplit
163
164	REAL zsin(iim),zcos(iim),z1(iim)
165	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
166	REAL unskap, pksurcp
167	c
168	cIM diagnostique PVteta, Amip2
169	INTEGER ntetaSTD
170	PARAMETER(ntetaSTD=3)
171	REAL rtetaSTD(ntetaSTD)
172	DATA rtetaSTD/350., 380., 405./
173	REAL PVteta(ngridmx,ntetaSTD)
174	c
175	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
176	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
177	c
178
179	REAL SSUM
180
181	LOGICAL firstcal, debut
182	DATA firstcal/.true./
183	SAVE firstcal,debut
184	! REAL rdayvrai
185	REAL, intent(in):: jD_cur, jH_cur
186	c
187	c-----------------------------------------------------------------------
188	c
189	c 1. Initialisations :
190	c --------------------
191	c
192	c
193	IF ( firstcal ) THEN
194	debut = .TRUE.
195	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
196	write(lunout,*) 'STOP dans calfis'
197	write(lunout,*)
198	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
199	write(lunout,*) ' ngridmx jjm iim '
200	write(lunout,*) ngridmx,jjm,iim
201	STOP
202	ENDIF
203	ELSE
204	debut = .FALSE.
205	ENDIF ! of IF (firstcal)
206
207	c
208	c
209	c-----------------------------------------------------------------------
210	c 40. transformation des variables dynamiques en variables physiques:
211	c ---------------------------------------------------------------
212
213	c 41. pressions au sol (en Pascals)
214	c ----------------------------------
215
216
217	zpsrf(1) = pps(1,1)
218
219	ig0 = 2
220	DO j = 2,jjm
221	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
222	ig0 = ig0+iim
223	ENDDO
224
225	zpsrf(ngridmx) = pps(1,jjp1)
226
227
228	c 42. pression intercouches et fonction d'Exner:
229	c
230	c -----------------------------------------------------------------
231	c .... zplev definis aux (llm +1) interfaces des couches ....
232	c .... zplay definis aux ( llm ) milieux des couches ....
233	c -----------------------------------------------------------------
234
235	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
236	c
237	unskap = 1./ kappa
238	c
239	! ADAPTATION GCM POUR CP(T)
240	DO l = 1, llm
241	zpk( 1,l ) = ppk(1,1,l)
242	zteta( 1,l ) = pteta(1,1,l)
243	zplev( 1,l ) = pp(1,1,l)
244	ig0 = 2
245	DO j = 2, jjm
246	DO i =1, iim
247	zpk( ig0,l ) = ppk(i,j,l)
248	zteta( ig0,l ) = pteta(i,j,l)
249	zplev( ig0,l ) = pp(i,j,l)
250	ig0 = ig0 +1
251	ENDDO
252	ENDDO
253	zpk( ngridmx,l ) = ppk(1,jjp1,l)
254	zteta( ngridmx,l ) = pteta(1,jjp1,l)
255	zplev( ngridmx,l ) = pp(1,jjp1,l)
256	ENDDO
257	zplev( 1,llmp1 ) = pp(1,1,llmp1)
258	ig0 = 2
259	DO j = 2, jjm
260	DO i =1, iim
261	zplev( ig0,llmp1 ) = pp(i,j,llmp1)
262	ig0 = ig0 +1
263	ENDDO
264	ENDDO
265	zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
266	c
267	c
268
269	c 43. temperature naturelle (en K) et pressions milieux couches .
270	c ---------------------------------------------------------------
271
272	! ADAPTATION GCM POUR CP(T)
273	call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
274
275	DO l=1,llm
276
277	pksurcp = ppk(1,1,l) / cpp
278	zplay(1,l) = preff * pksurcp ** unskap
279	! pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
280	ig0 = 2
281
282	DO j = 2, jjm
283	DO i = 1, iim
284	pksurcp = ppk(i,j,l) / cpp
285	zplay(ig0,l) = preff * pksurcp ** unskap
286	! pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
287	ig0 = ig0 + 1
288	ENDDO
289	ENDDO
290
291	pksurcp = ppk(1,jjp1,l) / cpp
292	zplay(ig0,l) = preff * pksurcp ** unskap
293	! pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
294
295	ENDDO
296
297	c 43.bis traceurs (tous intensifs)
298	c ---------------
299	c
300	DO iq=1,nqtot
301	iiq=niadv(iq)
302	DO l=1,llm
303	zqfi(1,l,iq) = pq(1,1,l,iiq)
304	ig0 = 2
305	DO j=2,jjm
306	DO i = 1, iim
307	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
308	ig0 = ig0 + 1
309	ENDDO
310	ENDDO
311	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
312	ENDDO
313	ENDDO ! boucle sur traceurs
314
315	!-----------------
316	! RQ SL 13/10/10:
317	! Ces calculs ne servent pas.
318	! Si necessaire, decommenter ces variables et les calculs...
319	!
320	! convergence dynamique pour les traceurs "EAU"
321	! Earth-specific treatment of first 2 tracers (water)
322	! if (planet_type=="earth") then
323	! DO iq=1,2
324	! DO l=1,llm
325	! pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
326	! ig0 = 2
327	! DO j=2,jjm
328	! DO i = 1, iim
329	! pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
330	! ig0 = ig0 + 1
331	! ENDDO
332	! ENDDO
333	! pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
334	! ENDDO
335	! ENDDO
336	! endif ! of if (planet_type=="earth")
337	!----------------
338
339	c Geopotentiel calcule par rapport a la surface locale:
340	c -----------------------------------------------------
341
342	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
343	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
344	DO l=1,llm
345	DO ig=1,ngridmx
346	zphi(ig,l)=zphi(ig,l)-zphis(ig)
347	ENDDO
348	ENDDO
349
350	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
351	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
352	c de masse est calclue dans advtrac.F
353	c DO l=1,llm
354	c pvervel(1,l)=pw(1,1,l) * g /apoln
355	c ig0=2
356	c DO j=2,jjm
357	c DO i = 1, iim
358	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
359	c ig0 = ig0 + 1
360	c ENDDO
361	c ENDDO
362	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
363	c ENDDO
364
365	c
366	c 45. champ u:
367	c ------------
368
369	DO 50 l=1,llm
370
371	DO 25 j=2,jjm
372	ig0 = 1+(j-2)*iim
373	zufi(ig0+1,l)= 0.5 *
374	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
375	! pcvgu(ig0+1,l)= 0.5 *
376	! $ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
377	DO 10 i=2,iim
378	zufi(ig0+i,l)= 0.5 *
379	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
380	! pcvgu(ig0+i,l)= 0.5 *
381	! $ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
382	10 CONTINUE
383	25 CONTINUE
384
385	50 CONTINUE
386
387
388	c 46.champ v:
389	c -----------
390
391	DO l=1,llm
392	DO j=2,jjm
393	ig0=1+(j-2)*iim
394	DO i=1,iim
395	zvfi(ig0+i,l)= 0.5 *
396	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
397	c pcvgv(ig0+i,l)= 0.5 *
398	c $ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
399	ENDDO
400	ENDDO
401	ENDDO
402
403
404	c 47. champs de vents aux pole nord
405	c ------------------------------
406	c U = 1 / pi * integrale [ v * cos(long) * d long ]
407	c V = 1 / pi * integrale [ v * sin(long) * d long ]
408
409	DO l=1,llm
410
411	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
412	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
413	DO i=2,iim
414	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
415	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
416	ENDDO
417
418	DO i=1,iim
419	zcos(i) = COS(rlonv(i))*z1(i)
420	zcosbis(i)= COS(rlonv(i))*z1bis(i)
421	zsin(i) = SIN(rlonv(i))*z1(i)
422	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
423	ENDDO
424
425	zufi(1,l) = SSUM(iim,zcos,1)/pi
426	! pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
427	zvfi(1,l) = SSUM(iim,zsin,1)/pi
428	! pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
429
430	ENDDO
431
432
433	c 48. champs de vents aux pole sud:
434	c ---------------------------------
435	c U = 1 / pi * integrale [ v * cos(long) * d long ]
436	c V = 1 / pi * integrale [ v * sin(long) * d long ]
437
438	DO l=1,llm
439
440	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
441	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
442	DO i=2,iim
443	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
444	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
445	ENDDO
446
447	DO i=1,iim
448	zcos(i) = COS(rlonv(i))*z1(i)
449	zcosbis(i) = COS(rlonv(i))*z1bis(i)
450	zsin(i) = SIN(rlonv(i))*z1(i)
451	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
452	ENDDO
453
454	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
455	! pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
456	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
457	! pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
458
459	ENDDO
460	c
461	if (planet_type=="earth") then
462	#ifdef CPP_EARTH
463	cIM calcul PV a teta=350, 380, 405K
464	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
465	$ ztfi,zplay,zplev,
466	$ ntetaSTD,rtetaSTD,PVteta)
467	#endif
468	endif
469	c
470	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
471	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
472
473	c-----------------------------------------------------------------------
474	c Appel de la physique:
475	c ---------------------
476
477	! Appel de la physique: pose probleme quand on tourne
478	! SANS physique, car physiq.F est dans le repertoire phy[]...
479	! Il faut une cle CPP_PHYS
480
481	! Le fait que les arguments de physiq soient differents selon les planetes
482	! ne pose pas de probleme a priori.
483
484	#ifdef CPP_PHYS
485
486	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
487	zdt_split=dtphys/nsplit_phys
488	zdufic(:,:)=0.
489	zdvfic(:,:)=0.
490	zdtfic(:,:)=0.
491	zdqfic(:,:,:)=0.
492
493	do isplit=1,nsplit_phys
494
495	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
496	debut_split=debut.and.isplit==1
497	lafin_split=lafin.and.isplit==nsplit_phys
498
499	if (planet_type.eq."earth") then
500	CALL physiq (ngridmx,
501	. llm,
502	. debut_split,
503	. lafin_split,
504	. jD_cur,
505	. jH_cur_split,
506	. zdt_split,
507	. zplev,
508	. zplay,
509	. zphi,
510	. zphis,
511	. presnivs,
512	. zufi,
513	. zvfi,
514	. ztfi,
515	. zqfi,
516	. flxwfi,
517	. zdufi,
518	. zdvfi,
519	. zdtfi,
520	. zdqfi,
521	. zdpsrf,
522	cIM diagnostique PVteta, Amip2
523	. pducov,
524	. PVteta)
525
526	else ! a moduler pour Mars !!
527
528	CALL physiq (ngridmx,
529	. llm,
530	. nqtot,
531	. debut_split,
532	. lafin_split,
533	. jD_cur,
534	. jH_cur_split,
535	. zdt_split,
536	. zplev,
537	. zplay,
538	. zpk,
539	. zphi,
540	. zphis,
541	. presnivs,
542	. zufi,
543	. zvfi,
544	. ztfi,
545	. zqfi,
546	. flxwfi,
547	. zdufi,
548	. zdvfi,
549	. zdtfi,
550	. zdqfi,
551	. zdpsrf)
552
553	endif ! planet_type
554
555	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
556	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
557	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
558	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
559
560	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
561	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
562	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
563	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
564
565	enddo
566	zdufi(:,:)=zdufic(:,:)/nsplit_phys
567	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
568	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
569	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
570
571	#endif
572	! #endif of #ifdef CPP_PHYS
573
574	500 CONTINUE
575
576	c-----------------------------------------------------------------------
577	c transformation des tendances physiques en tendances dynamiques:
578	c ---------------------------------------------------------------
579
580	c tendance sur la pression :
581	c -----------------------------------
582
583	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
584	c
585	c 62. enthalpie potentielle
586	c ---------------------
587
588	! ADAPTATION GCM POUR CP(T)
589	ztfi=ztfi+zdtfi*dtphys
590	call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
591
592	DO i=1,iip1
593	pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
594	pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
595	ENDDO
596
597	DO j=2,jjm
598	ig0=1+(j-2)*iim
599	DO i=1,iim
600	pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
601	ENDDO
602	pdhfi(iip1,j,:) = pdhfi(1,j,:)
603	ENDDO
604
605
606	c 62. humidite specifique
607	c ---------------------
608	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
609	! DO iq=1,nqtot
610	! DO l=1,llm
611	! DO i=1,iip1
612	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
613	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
614	! ENDDO
615	! DO j=2,jjm
616	! ig0=1+(j-2)*iim
617	! DO i=1,iim
618	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
619	! ENDDO
620	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
621	! ENDDO
622	! ENDDO
623	! ENDDO
624
625	c 63. traceurs (tous en intensifs)
626	c ------------
627	C initialisation des tendances
628	pdqfi(:,:,:,:)=0.
629	C
630	DO iq=1,nqtot
631	iiq=niadv(iq)
632	DO l=1,llm
633	DO i=1,iip1
634	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
635	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
636	ENDDO
637	DO j=2,jjm
638	ig0=1+(j-2)*iim
639	DO i=1,iim
640	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
641	ENDDO
642	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
643	ENDDO
644	ENDDO
645	ENDDO
646
647	c 65. champ u:
648	c ------------
649
650	DO l=1,llm
651
652	DO i=1,iip1
653	pdufi(i,1,l) = 0.
654	pdufi(i,jjp1,l) = 0.
655	ENDDO
656
657	DO j=2,jjm
658	ig0=1+(j-2)*iim
659	DO i=1,iim-1
660	pdufi(i,j,l)=
661	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
662	ENDDO
663	pdufi(iim,j,l)=
664	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
665	pdufi(iip1,j,l)=pdufi(1,j,l)
666	ENDDO
667
668	ENDDO
669
670
671	c 67. champ v:
672	c ------------
673
674	DO l=1,llm
675
676	DO j=2,jjm-1
677	ig0=1+(j-2)*iim
678	DO i=1,iim
679	pdvfi(i,j,l)=
680	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
681	ENDDO
682	pdvfi(iip1,j,l) = pdvfi(1,j,l)
683	ENDDO
684	ENDDO
685
686
687	c 68. champ v pres des poles:
688	c ---------------------------
689	c v = U * cos(long) + V * SIN(long)
690
691	DO l=1,llm
692
693	DO i=1,iim
694	pdvfi(i,1,l)=
695	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
696	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
697	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
698	pdvfi(i,1,l)=
699	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
700	pdvfi(i,jjm,l)=
701	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
702	ENDDO
703
704	pdvfi(iip1,1,l) = pdvfi(1,1,l)
705	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
706
707	ENDDO
708
709	c-----------------------------------------------------------------------
710
711	700 CONTINUE
712
713	firstcal = .FALSE.
714
715	RETURN
716	END

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

Download in other formats:

Original Format