source: LMDZ5/trunk/libf/dynlonlat_phylonlat/calfis.F @ 2246

Last change on this file since 2246 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

  • what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
  • 1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
  • "bibio" directory is now "misc" and should only contain autonomous utilities.
  • "cosp" is now a subdirectory of phylmd.

EM

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