source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3d/calfis.F @ 5394

Last change on this file since 5394 was 1676, checked in by aslmd, 12 years ago

  • Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation


  • Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon


Author : AS equipe planeto

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.4 KB
Line 
1!
2! $Id: calfis.F 1676 2012-11-08 20:03:30Z evignon $
3!
4C
5C
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)
29c
30c    Auteur :  P. Le Van, F. Hourdin
31c   .........
32      USE infotrac
33      USE control_mod
34 
35
36      IMPLICIT NONE
37c=======================================================================
38c
39c   1. rearrangement des tableaux et transformation
40c      variables dynamiques  >  variables physiques
41c   2. calcul des termes physiques
42c   3. retransformation des tendances physiques en tendances dynamiques
43c
44c   remarques:
45c   ----------
46c
47c    - les vents sont donnes dans la physique par leurs composantes
48c      naturelles.
49c    - la variable thermodynamique de la physique est une variable
50c      intensive :   T
51c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
52c    - les deux seules variables dependant de la geometrie necessaires
53c      pour la physique sont la latitude pour le rayonnement et
54c      l'aire de la maille quand on veut integrer une grandeur
55c      horizontalement.
56c    - les points de la physique sont les points scalaires de la
57c      la dynamique; numerotation:
58c          1 pour le pole nord
59c          (jjm-1)*iim pour l'interieur du domaine
60c          ngridmx pour le pole sud
61c      ---> ngridmx=2+(jjm-1)*iim
62c
63c     Input :
64c     -------
65c       pucov           covariant zonal velocity
66c       pvcov           covariant meridional velocity
67c       pteta           potential temperature
68c       pps             surface pressure
69c       pmasse          masse d'air dans chaque maille
70c       pts             surface temperature  (K)
71c       callrad         clef d'appel au rayonnement
72c
73c    Output :
74c    --------
75c        pdufi          tendency for the natural zonal velocity (ms-1)
76c        pdvfi          tendency for the natural meridional velocity
77c        pdhfi          tendency for the potential temperature
78c        pdtsfi         tendency for the surface temperature
79c
80c        pdtrad         radiative tendencies  \  both input
81c        pfluxrad       radiative fluxes      /  and output
82c
83c=======================================================================
84c
85c-----------------------------------------------------------------------
86c
87c    0.  Declarations :
88c    ------------------
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
102c    Arguments :
103c    -----------
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)
114c
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)
119c
120      REAL pps(iip1,jjp1)
121      REAL pp(iip1,jjp1,llmp1)
122      REAL ppk(iip1,jjp1,llm)
123c
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
135c    Local variables :
136c    -----------------
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)
142c
143      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
144      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
145c
146      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
147      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
148c
149      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
150      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
151      REAL zdpsrf(ngridmx)
152c
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
162c
163cIM 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)
169c
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
172c
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
184c
185c-----------------------------------------------------------------------
186c
187c    1. Initialisations :
188c    --------------------
189c
190c
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
205c
206c
207c-----------------------------------------------------------------------
208c   40. transformation des variables dynamiques en variables physiques:
209c   ---------------------------------------------------------------
210
211c   41. pressions au sol (en Pascals)
212c   ----------------------------------
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
226c   42. pression intercouches :
227c
228c   -----------------------------------------------------------------
229c     .... zplev  definis aux (llm +1) interfaces des couches  ....
230c     .... zplay  definis aux (  llm )    milieux des couches  ....
231c   -----------------------------------------------------------------
232
233c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
234c
235       unskap   = 1./ kappa
236c
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
248c
249c
250
251c   43. temperature naturelle (en K) et pressions milieux couches .
252c   ---------------------------------------------------------------
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
279c   43.bis traceurs
280c   ---------------
281c
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
297c   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
316c   Geopotentiel calcule par rapport a la surface locale:
317c   -----------------------------------------------------
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
327c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
328c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
329c    de masse est calclue dans advtrac.F 
330c      DO l=1,llm
331c        pvervel(1,l)=pw(1,1,l) * g /apoln
332c        ig0=2
333c       DO j=2,jjm
334c           DO i = 1, iim
335c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
336c              ig0 = ig0 + 1
337c           ENDDO
338c       ENDDO
339c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
340c      ENDDO
341
342c
343c   45. champ u:
344c   ------------
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) )
35910         CONTINUE
36025      CONTINUE
361
36250    CONTINUE
363
364
365c   46.champ v:
366c   -----------
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
381c   47. champs de vents aux pole nord   
382c   ------------------------------
383c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
384c        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
410c   48. champs de vents aux pole sud:
411c   ---------------------------------
412c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
413c        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
437c
438      if (planet_type=="earth") then
439#ifdef CPP_PHYS
440! PVtheta calls tetalevel, which is in the physics
441cIM 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
447c
448c 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
451c-----------------------------------------------------------------------
452c   Appel de la physique:
453c   ---------------------
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,
497cIM 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
550500   CONTINUE
551
552c-----------------------------------------------------------------------
553c   transformation des tendances physiques en tendances dynamiques:
554c   ---------------------------------------------------------------
555
556c  tendance sur la pression :
557c  -----------------------------------
558
559      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
560c
561c   62. enthalpie potentielle
562c   ---------------------
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
582c   62. humidite specifique
583c   ---------------------
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
601c   63. traceurs
602c   ------------
603C     initialisation des tendances
604      pdqfi(:,:,:,:)=0.
605C
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
623c   65. champ u:
624c   ------------
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
647c   67. champ v:
648c   ------------
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
663c   68. champ v pres des poles:
664c   ---------------------------
665c      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
685c-----------------------------------------------------------------------
686
687700   CONTINUE
688 
689      firstcal = .FALSE.
690
691      RETURN
692      END
Note: See TracBrowser for help on using the repository browser.