source: LMDZ5/branches/testing/libf/dyn3d/calfis.F @ 1999

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

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