source: trunk/LMDZ.COMMON/libf/dyn3d/calfis.F @ 801

Last change on this file since 801 was 776, checked in by emillour, 12 years ago

Common dynamics: updates to keep up with LMDZ5 Earth (rev 1649)
See file "DOC/chantiers/commit_importants.log" for details.
EM

File size: 20.1 KB
Line 
1!
2! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
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     $                  pdufi,
24     $                  pdvfi,
25     $                  pdhfi,
26     $                  pdqfi,
27     $                  pdpsfi)
28c
29c    Auteur :  P. Le Van, F. Hourdin
30c   .........
31      USE infotrac
32      USE control_mod
33 
34
35      IMPLICIT NONE
36c=======================================================================
37c
38c   1. rearrangement des tableaux et transformation
39c      variables dynamiques  >  variables physiques
40c   2. calcul des termes physiques
41c   3. retransformation des tendances physiques en tendances dynamiques
42c
43c   remarques:
44c   ----------
45c
46c    - les vents sont donnes dans la physique par leurs composantes
47c      naturelles.
48c    - la variable thermodynamique de la physique est une variable
49c      intensive :   T
50c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
51c    - les deux seules variables dependant de la geometrie necessaires
52c      pour la physique sont la latitude pour le rayonnement et
53c      l'aire de la maille quand on veut integrer une grandeur
54c      horizontalement.
55c    - les points de la physique sont les points scalaires de la
56c      la dynamique; numerotation:
57c          1 pour le pole nord
58c          (jjm-1)*iim pour l'interieur du domaine
59c          ngridmx pour le pole sud
60c      ---> ngridmx=2+(jjm-1)*iim
61c
62c     Input :
63c     -------
64c       pucov           covariant zonal velocity
65c       pvcov           covariant meridional velocity
66c       pteta           potential temperature
67c       pps             surface pressure
68c       pmasse          masse d'air dans chaque maille
69c       pts             surface temperature  (K)
70c       callrad         clef d'appel au rayonnement
71c
72c    Output :
73c    --------
74c        pdufi          tendency for the natural zonal velocity (ms-1)
75c        pdvfi          tendency for the natural meridional velocity
76c        pdhfi          tendency for the potential temperature (K/s)
77c        pdtsfi         tendency for the surface temperature
78c
79c        pdtrad         radiative tendencies  \  both input
80c        pfluxrad       radiative fluxes      /  and output
81c
82c=======================================================================
83c
84c-----------------------------------------------------------------------
85c
86c    0.  Declarations :
87c    ------------------
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
101c    Arguments :
102c    -----------
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)
113c
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)
121c
122      REAL pps(iip1,jjp1)
123      REAL pp(iip1,jjp1,llmp1)
124      REAL ppk(iip1,jjp1,llm)
125c
126c 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
134c    Local variables :
135c    -----------------
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)
141c
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)
147c
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)
153c
154      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
155      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
156      REAL zdpsrf(ngridmx)
157c
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
167c
168cIM diagnostique PVteta, Amip2
169      INTEGER ntetaSTD
170      PARAMETER(ntetaSTD=3)
171      REAL rtetaSTD(ntetaSTD)
172      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
173      REAL PVteta(ngridmx,ntetaSTD)
174c
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
177c
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
186c
187c-----------------------------------------------------------------------
188c
189c    1. Initialisations :
190c    --------------------
191c
192c
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
207c
208c
209c-----------------------------------------------------------------------
210c   40. transformation des variables dynamiques en variables physiques:
211c   ---------------------------------------------------------------
212
213c   41. pressions au sol (en Pascals)
214c   ----------------------------------
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
228c   42. pression intercouches et fonction d'Exner:
229c
230c   -----------------------------------------------------------------
231c     .... zplev  definis aux (llm +1) interfaces des couches  ....
232c     .... zplay  definis aux (  llm )    milieux des couches  ....
233c   -----------------------------------------------------------------
234
235c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
236c
237       unskap   = 1./ kappa
238c
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)
266c
267c
268
269c   43. temperature naturelle (en K) et pressions milieux couches .
270c   ---------------------------------------------------------------
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
297c   43.bis traceurs (tous intensifs)
298c   ---------------
299c
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
339c   Geopotentiel calcule par rapport a la surface locale:
340c   -----------------------------------------------------
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
350c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
351c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
352c    de masse est calclue dans advtrac.F 
353c      DO l=1,llm
354c        pvervel(1,l)=pw(1,1,l) * g /apoln
355c        ig0=2
356c       DO j=2,jjm
357c           DO i = 1, iim
358c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
359c              ig0 = ig0 + 1
360c           ENDDO
361c       ENDDO
362c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
363c      ENDDO
364
365c
366c   45. champ u:
367c   ------------
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) )
38210         CONTINUE
38325      CONTINUE
384
38550    CONTINUE
386
387
388c   46.champ v:
389c   -----------
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) )
397c               pcvgv(ig0+i,l)= 0.5 *
398c     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
399            ENDDO
400         ENDDO
401      ENDDO
402
403
404c   47. champs de vents aux pole nord   
405c   ------------------------------
406c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
407c        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
433c   48. champs de vents aux pole sud:
434c   ---------------------------------
435c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
436c        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
460c
461      if (planet_type=="earth") then
462#ifdef CPP_EARTH
463! PVtheta calls tetalevel, which is in the (Earth) physics
464cIM calcul PV a teta=350, 380, 405K
465      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
466     $           ztfi,zplay,zplev,
467     $           ntetaSTD,rtetaSTD,PVteta)
468#endif
469      endif
470c
471c On change de grille, dynamique vers physiq, pour le flux de masse verticale
472      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
473
474c-----------------------------------------------------------------------
475c   Appel de la physique:
476c   ---------------------
477
478! Appel de la physique: pose probleme quand on tourne
479! SANS physique, car physiq.F est dans le repertoire phy[]...
480! Il faut une cle CPP_PHYS
481
482! Le fait que les arguments de physiq soient differents selon les planetes
483! ne pose pas de probleme a priori.
484
485!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
486      zdt_split=dtphys/nsplit_phys
487      zdufic(:,:)=0.
488      zdvfic(:,:)=0.
489      zdtfic(:,:)=0.
490      zdqfic(:,:,:)=0.
491
492#ifdef CPP_PHYS
493
494      do isplit=1,nsplit_phys
495
496         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
497         debut_split=debut.and.isplit==1
498         lafin_split=lafin.and.isplit==nsplit_phys
499
500      if (planet_type.eq."earth") then
501         CALL physiq (ngridmx,
502     .             llm,
503     .             debut_split,
504     .             lafin_split,
505     .             jD_cur,
506     .             jH_cur_split,
507     .             zdt_split,
508     .             zplev,
509     .             zplay,
510     .             zphi,
511     .             zphis,
512     .             presnivs,
513     .             zufi,
514     .             zvfi,
515     .             ztfi,
516     .             zqfi,
517     .             flxwfi,
518     .             zdufi,
519     .             zdvfi,
520     .             zdtfi,
521     .             zdqfi,
522     .             zdpsrf,
523cIM diagnostique PVteta, Amip2         
524     .             pducov,
525     .             PVteta)
526
527      else  ! a moduler pour Mars !!
528
529         CALL physiq (ngridmx,
530     .             llm,
531     .             nqtot,
532     .             debut_split,
533     .             lafin_split,
534     .             jD_cur,
535     .             jH_cur_split,
536     .             zdt_split,
537     .             zplev,
538     .             zplay,
539     .             zpk,
540     .             zphi,
541     .             zphis,
542     .             presnivs,
543     .             zufi,
544     .             zvfi,
545     .             ztfi,
546     .             zqfi,
547     .             flxwfi,
548     .             zdufi,
549     .             zdvfi,
550     .             zdtfi,
551     .             zdqfi,
552     .             zdpsrf)
553
554      endif ! planet_type
555
556         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
557         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
558         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
559         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
560
561         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
562         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
563         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
564         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
565
566      enddo ! of do isplit=1,nsplit_phys
567
568#endif
569! #endif of #ifdef CPP_PHYS
570
571      zdufi(:,:)=zdufic(:,:)/nsplit_phys
572      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
573      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
574      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
575
576
577500   CONTINUE
578
579c-----------------------------------------------------------------------
580c   transformation des tendances physiques en tendances dynamiques:
581c   ---------------------------------------------------------------
582
583c  tendance sur la pression :
584c  -----------------------------------
585
586      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
587c
588c   62. enthalpie potentielle
589c   ---------------------
590
591! ADAPTATION GCM POUR CP(T)
592      call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
593
594         DO i=1,iip1
595      pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
596      pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
597         ENDDO
598
599         DO j=2,jjm
600            ig0=1+(j-2)*iim
601            DO i=1,iim
602      pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
603            ENDDO
604               pdhfi(iip1,j,:) =  pdhfi(1,j,:)
605         ENDDO
606
607
608c   62. humidite specifique
609c   ---------------------
610! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
611!      DO iq=1,nqtot
612!         DO l=1,llm
613!            DO i=1,iip1
614!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
615!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
616!            ENDDO
617!            DO j=2,jjm
618!               ig0=1+(j-2)*iim
619!               DO i=1,iim
620!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
621!               ENDDO
622!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
623!            ENDDO
624!         ENDDO
625!      ENDDO
626
627c   63. traceurs (tous en intensifs)
628c   ------------
629C     initialisation des tendances
630      pdqfi(:,:,:,:)=0.
631C
632       DO iq=1,nqtot
633         iiq=niadv(iq)
634         DO l=1,llm
635            DO i=1,iip1
636               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
637               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
638            ENDDO
639            DO j=2,jjm
640               ig0=1+(j-2)*iim
641               DO i=1,iim
642                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
643               ENDDO
644               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
645            ENDDO
646         ENDDO
647       ENDDO
648
649c   65. champ u:
650c   ------------
651
652      DO l=1,llm
653
654         DO i=1,iip1
655            pdufi(i,1,l)    = 0.
656            pdufi(i,jjp1,l) = 0.
657         ENDDO
658
659         DO j=2,jjm
660            ig0=1+(j-2)*iim
661            DO i=1,iim-1
662               pdufi(i,j,l)=
663     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
664            ENDDO
665            pdufi(iim,j,l)=
666     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
667            pdufi(iip1,j,l)=pdufi(1,j,l)
668         ENDDO
669
670      ENDDO
671
672
673c   67. champ v:
674c   ------------
675
676      DO l=1,llm
677
678         DO j=2,jjm-1
679            ig0=1+(j-2)*iim
680            DO i=1,iim
681               pdvfi(i,j,l)=
682     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
683            ENDDO
684            pdvfi(iip1,j,l) = pdvfi(1,j,l)
685         ENDDO
686      ENDDO
687
688
689c   68. champ v pres des poles:
690c   ---------------------------
691c      v = U * cos(long) + V * SIN(long)
692
693      DO l=1,llm
694
695         DO i=1,iim
696            pdvfi(i,1,l)=
697     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
698            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
699     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
700            pdvfi(i,1,l)=
701     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
702            pdvfi(i,jjm,l)=
703     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
704          ENDDO
705
706         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
707         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
708
709      ENDDO
710
711c-----------------------------------------------------------------------
712
713700   CONTINUE
714 
715      firstcal = .FALSE.
716
717      RETURN
718      END
Note: See TracBrowser for help on using the repository browser.