Context Navigation

calfis.F @ 108

Last change on this file since 108 was 108, checked in by slebonnois, 14 years ago
Sebastien Lebonnois: sponge layer et dissip horizontale.
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	call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
590
591	DO i=1,iip1
592	pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
593	pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
594	ENDDO
595
596	DO j=2,jjm
597	ig0=1+(j-2)*iim
598	DO i=1,iim
599	pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
600	ENDDO
601	pdhfi(iip1,j,:) = pdhfi(1,j,:)
602	ENDDO
603
604
605	c 62. humidite specifique
606	c ---------------------
607	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
608	! DO iq=1,nqtot
609	! DO l=1,llm
610	! DO i=1,iip1
611	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
612	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
613	! ENDDO
614	! DO j=2,jjm
615	! ig0=1+(j-2)*iim
616	! DO i=1,iim
617	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
618	! ENDDO
619	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
620	! ENDDO
621	! ENDDO
622	! ENDDO
623
624	c 63. traceurs (tous en intensifs)
625	c ------------
626	C initialisation des tendances
627	pdqfi(:,:,:,:)=0.
628	C
629	DO iq=1,nqtot
630	iiq=niadv(iq)
631	DO l=1,llm
632	DO i=1,iip1
633	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
634	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
635	ENDDO
636	DO j=2,jjm
637	ig0=1+(j-2)*iim
638	DO i=1,iim
639	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
640	ENDDO
641	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
642	ENDDO
643	ENDDO
644	ENDDO
645
646	c 65. champ u:
647	c ------------
648
649	DO l=1,llm
650
651	DO i=1,iip1
652	pdufi(i,1,l) = 0.
653	pdufi(i,jjp1,l) = 0.
654	ENDDO
655
656	DO j=2,jjm
657	ig0=1+(j-2)*iim
658	DO i=1,iim-1
659	pdufi(i,j,l)=
660	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
661	ENDDO
662	pdufi(iim,j,l)=
663	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
664	pdufi(iip1,j,l)=pdufi(1,j,l)
665	ENDDO
666
667	ENDDO
668
669
670	c 67. champ v:
671	c ------------
672
673	DO l=1,llm
674
675	DO j=2,jjm-1
676	ig0=1+(j-2)*iim
677	DO i=1,iim
678	pdvfi(i,j,l)=
679	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
680	ENDDO
681	pdvfi(iip1,j,l) = pdvfi(1,j,l)
682	ENDDO
683	ENDDO
684
685
686	c 68. champ v pres des poles:
687	c ---------------------------
688	c v = U * cos(long) + V * SIN(long)
689
690	DO l=1,llm
691
692	DO i=1,iim
693	pdvfi(i,1,l)=
694	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
695	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
696	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
697	pdvfi(i,1,l)=
698	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
699	pdvfi(i,jjm,l)=
700	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
701	ENDDO
702
703	pdvfi(iip1,1,l) = pdvfi(1,1,l)
704	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
705
706	ENDDO
707
708	c-----------------------------------------------------------------------
709
710	700 CONTINUE
711
712	firstcal = .FALSE.
713
714	RETURN
715	END

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

Download in other formats:

Original Format