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