source: LMDZ5/branches/LMDZ6_rc0/libf/dyn3d/calfis.F @ 2461

Last change on this file since 2461 was 2056, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1997:2055 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 KB
Line 
1!
2! $Id: calfis.F 2056 2014-06-11 13:46:46Z fhourdin $
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, ONLY: nqtot, niadv, tname
33      USE control_mod, ONLY: planet_type, nsplit_phys
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,INTENT(IN) ::  lafin ! .true. for the very last call to physics
105      REAL,INTENT(IN):: jD_cur, jH_cur
106      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
107      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
108      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
109      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
110      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
111      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
112      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
113
114      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
115      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
116      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
117      ! NB: pdteta is used only to compute pcvgt which is in fact not used...
118      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
119      ! NB: pdq is only used to compute pcvgq which is in fact not used...
120
121      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
122      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
123      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
124      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm)  ! Vertical mass flux on dynamics grid
125
126      ! tendencies (in */s) from the physics
127      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
128      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
129      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
130      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
131      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
132
133      INTEGER,PARAMETER :: longcles = 20
134      REAL,INTENT(IN) :: clesphy0( longcles ) ! unused
135
136
137c    Local variables :
138c    -----------------
139
140      INTEGER i,j,l,ig0,ig,iq,iiq
141      REAL zpsrf(ngridmx)
142      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
143      REAL zphi(ngridmx,llm),zphis(ngridmx)
144c
145      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
146      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
147c
148      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
149      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
150c
151      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
152      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
153      REAL zdpsrf(ngridmx)
154c
155      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
156      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
157      REAL jH_cur_split,zdt_split
158      LOGICAL debut_split,lafin_split
159      INTEGER isplit
160
161      REAL zsin(iim),zcos(iim),z1(iim)
162      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
163      REAL unskap, pksurcp
164c
165      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
166c
167     
168      REAL SSUM
169
170      LOGICAL,SAVE :: firstcal=.true., debut=.true.
171!      REAL rdayvrai
172
173      LOGICAL tracerdyn
174
175c
176c-----------------------------------------------------------------------
177c
178c    1. Initialisations :
179c    --------------------
180c
181c
182      IF ( firstcal )  THEN
183        debut = .TRUE.
184        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
185         write(lunout,*) 'STOP dans calfis'
186         write(lunout,*)
187     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
188         write(lunout,*) '  ngridmx  jjm   iim   '
189         write(lunout,*) ngridmx,jjm,iim
190         STOP
191        ENDIF
192      ELSE
193        debut = .FALSE.
194      ENDIF ! of IF (firstcal)
195
196c
197c
198c-----------------------------------------------------------------------
199c   40. transformation des variables dynamiques en variables physiques:
200c   ---------------------------------------------------------------
201
202c   41. pressions au sol (en Pascals)
203c   ----------------------------------
204
205       
206      zpsrf(1) = pps(1,1)
207
208      ig0  = 2
209      DO j = 2,jjm
210         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
211         ig0 = ig0+iim
212      ENDDO
213
214      zpsrf(ngridmx) = pps(1,jjp1)
215
216
217c   42. pression intercouches :
218c
219c   -----------------------------------------------------------------
220c     .... zplev  definis aux (llm +1) interfaces des couches  ....
221c     .... zplay  definis aux (  llm )    milieux des couches  ....
222c   -----------------------------------------------------------------
223
224c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
225c
226       unskap   = 1./ kappa
227c
228      DO l = 1, llmp1
229        zplev( 1,l ) = pp(1,1,l)
230        ig0 = 2
231          DO j = 2, jjm
232             DO i =1, iim
233              zplev( ig0,l ) = pp(i,j,l)
234              ig0 = ig0 +1
235             ENDDO
236          ENDDO
237        zplev( ngridmx,l ) = pp(1,jjp1,l)
238      ENDDO
239c
240c
241
242c   43. temperature naturelle (en K) et pressions milieux couches .
243c   ---------------------------------------------------------------
244
245      DO l=1,llm
246
247         pksurcp     =  ppk(1,1,l) / cpp
248         zplay(1,l)  =  preff * pksurcp ** unskap
249         ztfi(1,l)   =  pteta(1,1,l) *  pksurcp
250         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
251         ig0         = 2
252
253         DO j = 2, jjm
254            DO i = 1, iim
255              pksurcp        = ppk(i,j,l) / cpp
256              zplay(ig0,l)   = preff * pksurcp ** unskap
257              ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
258              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
259              ig0            = ig0 + 1
260            ENDDO
261         ENDDO
262
263         pksurcp       = ppk(1,jjp1,l) / cpp
264         zplay(ig0,l)  = preff * pksurcp ** unskap
265         ztfi (ig0,l)  = pteta(1,jjp1,l)  * pksurcp
266         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
267
268      ENDDO
269
270c   43.bis traceurs
271c   ---------------
272c
273      DO iq=1,nqtot
274          iiq=niadv(iq)
275         DO l=1,llm
276            zqfi(1,l,iq) = pq(1,1,l,iiq)
277            ig0          = 2
278            DO j=2,jjm
279               DO i = 1, iim
280                  zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
281                  ig0             = ig0 + 1
282               ENDDO
283            ENDDO
284            zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
285         ENDDO
286      ENDDO
287
288c   convergence dynamique pour les traceurs "EAU"
289! Earth-specific treatment of first 2 tracers (water)
290       if (planet_type=="earth") then
291        DO iq=1,2
292         DO l=1,llm
293            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
294            ig0          = 2
295            DO j=2,jjm
296               DO i = 1, iim
297                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
298                  ig0             = ig0 + 1
299               ENDDO
300            ENDDO
301            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
302         ENDDO
303        ENDDO
304       endif ! of if (planet_type=="earth")
305
306
307c   Geopotentiel calcule par rapport a la surface locale:
308c   -----------------------------------------------------
309
310      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
311      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
312      DO l=1,llm
313         DO ig=1,ngridmx
314           zphi(ig,l)=zphi(ig,l)-zphis(ig)
315         ENDDO
316      ENDDO
317
318c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
319c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
320c    de masse est calclue dans advtrac.F 
321c      DO l=1,llm
322c        pvervel(1,l)=pw(1,1,l) * g /apoln
323c        ig0=2
324c       DO j=2,jjm
325c           DO i = 1, iim
326c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
327c              ig0 = ig0 + 1
328c           ENDDO
329c       ENDDO
330c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
331c      ENDDO
332
333c
334c   45. champ u:
335c   ------------
336
337      DO 50 l=1,llm
338
339         DO 25 j=2,jjm
340            ig0 = 1+(j-2)*iim
341            zufi(ig0+1,l)= 0.5 *
342     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
343            pcvgu(ig0+1,l)= 0.5 *
344     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
345            DO 10 i=2,iim
346               zufi(ig0+i,l)= 0.5 *
347     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
348               pcvgu(ig0+i,l)= 0.5 *
349     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
35010         CONTINUE
35125      CONTINUE
352
35350    CONTINUE
354
355
356c   46.champ v:
357c   -----------
358
359      DO l=1,llm
360         DO j=2,jjm
361            ig0=1+(j-2)*iim
362            DO i=1,iim
363               zvfi(ig0+i,l)= 0.5 *
364     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
365               pcvgv(ig0+i,l)= 0.5 *
366     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
367            ENDDO
368         ENDDO
369      ENDDO
370
371
372c   47. champs de vents aux pole nord   
373c   ------------------------------
374c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
375c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
376
377      DO l=1,llm
378
379         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
380         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
381         DO i=2,iim
382            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
383            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
384         ENDDO
385
386         DO i=1,iim
387            zcos(i)   = COS(rlonv(i))*z1(i)
388            zcosbis(i)= COS(rlonv(i))*z1bis(i)
389            zsin(i)   = SIN(rlonv(i))*z1(i)
390            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
391         ENDDO
392
393         zufi(1,l)  = SSUM(iim,zcos,1)/pi
394         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
395         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
396         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
397
398      ENDDO
399
400
401c   48. champs de vents aux pole sud:
402c   ---------------------------------
403c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
404c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
405
406      DO l=1,llm
407
408         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
409         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
410         DO i=2,iim
411            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
412            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
413         ENDDO
414
415         DO i=1,iim
416            zcos(i)    = COS(rlonv(i))*z1(i)
417            zcosbis(i) = COS(rlonv(i))*z1bis(i)
418            zsin(i)    = SIN(rlonv(i))*z1(i)
419            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
420         ENDDO
421
422         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
423         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
424         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
425         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
426
427      ENDDO
428c
429c On change de grille, dynamique vers physiq, pour le flux de masse verticale
430      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
431
432c-----------------------------------------------------------------------
433c   Appel de la physique:
434c   ---------------------
435
436
437
438!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
439      zdt_split=dtphys/nsplit_phys
440      zdufic(:,:)=0.
441      zdvfic(:,:)=0.
442      zdtfic(:,:)=0.
443      zdqfic(:,:,:)=0.
444
445#ifdef CPP_PHYS
446
447       do isplit=1,nsplit_phys
448
449         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
450         debut_split=debut.and.isplit==1
451         lafin_split=lafin.and.isplit==nsplit_phys
452
453      if (planet_type=="earth") then
454
455         CALL physiq (ngridmx,
456     .             llm,
457     .             debut_split,
458     .             lafin_split,
459     .             jD_cur,
460     .             jH_cur_split,
461     .             zdt_split,
462     .             zplev,
463     .             zplay,
464     .             zphi,
465     .             zphis,
466     .             presnivs,
467     .             clesphy0,
468     .             zufi,
469     .             zvfi,
470     .             ztfi,
471     .             zqfi,
472     .             flxwfi,
473     .             zdufi,
474     .             zdvfi,
475     .             zdtfi,
476     .             zdqfi,
477     .             zdpsrf,
478     .             pducov)
479
480      else if ( planet_type=="generic" ) then
481
482         CALL physiq (ngridmx,     !! ngrid
483     .             llm,            !! nlayer
484     .             nqtot,          !! nq
485     .             tname,          !! tracer names from dynamical core (given in infotrac)
486     .             debut_split,    !! firstcall
487     .             lafin_split,    !! lastcall
488     .             jD_cur,         !! pday. see leapfrog
489     .             jH_cur_split,   !! ptime "fraction of day"
490     .             zdt_split,      !! ptimestep
491     .             zplev,          !! pplev
492     .             zplay,          !! pplay
493     .             zphi,           !! pphi
494     .             zufi,           !! pu
495     .             zvfi,           !! pv
496     .             ztfi,           !! pt
497     .             zqfi,           !! pq
498     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
499     .             zdufi,          !! pdu
500     .             zdvfi,          !! pdv
501     .             zdtfi,          !! pdt
502     .             zdqfi,          !! pdq
503     .             zdpsrf,         !! pdpsrf
504     .             tracerdyn)      !! tracerdyn <-- utilite ???
505
506      endif ! of if (planet_type=="earth")
507
508         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
509         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
510         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
511         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
512
513         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
514         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
515         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
516         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
517
518       enddo ! of do isplit=1,nsplit_phys
519
520#endif
521! of #ifdef CPP_PHYS
522
523      zdufi(:,:)=zdufic(:,:)/nsplit_phys
524      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
525      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
526      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
527
528
529500   CONTINUE
530
531c-----------------------------------------------------------------------
532c   transformation des tendances physiques en tendances dynamiques:
533c   ---------------------------------------------------------------
534
535c  tendance sur la pression :
536c  -----------------------------------
537
538      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
539c
540c   62. enthalpie potentielle
541c   ---------------------
542
543      DO l=1,llm
544
545         DO i=1,iip1
546          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
547          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
548         ENDDO
549
550         DO j=2,jjm
551            ig0=1+(j-2)*iim
552            DO i=1,iim
553               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
554            ENDDO
555               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
556         ENDDO
557
558      ENDDO
559
560
561c   62. humidite specifique
562c   ---------------------
563! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
564!      DO iq=1,nqtot
565!         DO l=1,llm
566!            DO i=1,iip1
567!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
568!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
569!            ENDDO
570!            DO j=2,jjm
571!               ig0=1+(j-2)*iim
572!               DO i=1,iim
573!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
574!               ENDDO
575!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
576!            ENDDO
577!         ENDDO
578!      ENDDO
579
580c   63. traceurs
581c   ------------
582C     initialisation des tendances
583      pdqfi(:,:,:,:)=0.
584C
585      DO iq=1,nqtot
586         iiq=niadv(iq)
587         DO l=1,llm
588            DO i=1,iip1
589               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
590               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
591            ENDDO
592            DO j=2,jjm
593               ig0=1+(j-2)*iim
594               DO i=1,iim
595                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
596               ENDDO
597               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
598            ENDDO
599         ENDDO
600      ENDDO
601
602c   65. champ u:
603c   ------------
604
605      DO l=1,llm
606
607         DO i=1,iip1
608            pdufi(i,1,l)    = 0.
609            pdufi(i,jjp1,l) = 0.
610         ENDDO
611
612         DO j=2,jjm
613            ig0=1+(j-2)*iim
614            DO i=1,iim-1
615               pdufi(i,j,l)=
616     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
617            ENDDO
618            pdufi(iim,j,l)=
619     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
620            pdufi(iip1,j,l)=pdufi(1,j,l)
621         ENDDO
622
623      ENDDO
624
625
626c   67. champ v:
627c   ------------
628
629      DO l=1,llm
630
631         DO j=2,jjm-1
632            ig0=1+(j-2)*iim
633            DO i=1,iim
634               pdvfi(i,j,l)=
635     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
636            ENDDO
637            pdvfi(iip1,j,l) = pdvfi(1,j,l)
638         ENDDO
639      ENDDO
640
641
642c   68. champ v pres des poles:
643c   ---------------------------
644c      v = U * cos(long) + V * SIN(long)
645
646      DO l=1,llm
647
648         DO i=1,iim
649            pdvfi(i,1,l)=
650     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
651            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
652     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
653            pdvfi(i,1,l)=
654     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
655            pdvfi(i,jjm,l)=
656     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
657          ENDDO
658
659         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
660         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
661
662      ENDDO
663
664c-----------------------------------------------------------------------
665
666700   CONTINUE
667 
668      firstcal = .FALSE.
669
670      RETURN
671      END
Note: See TracBrowser for help on using the repository browser.