source: LMDZ5/trunk/libf/dynphy_lonlat/calfis.F @ 2601

Last change on this file since 2601 was 2601, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn temps.h into module temps_mod.F90
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: 21.3 KB
Line 
1!
2! $Id: calfis.F 2601 2016-07-24 09:51:55Z emillour $
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#ifdef CPP_PHYS
34      USE callphysiq_mod, ONLY: call_physiq
35#endif
36      USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
37      USE comvert_mod, ONLY: preff, presnivs
38     
39      IMPLICIT NONE
40c=======================================================================
41c
42c   1. rearrangement des tableaux et transformation
43c      variables dynamiques  >  variables physiques
44c   2. calcul des termes physiques
45c   3. retransformation des tendances physiques en tendances dynamiques
46c
47c   remarques:
48c   ----------
49c
50c    - les vents sont donnes dans la physique par leurs composantes
51c      naturelles.
52c    - la variable thermodynamique de la physique est une variable
53c      intensive :   T
54c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
55c    - les deux seules variables dependant de la geometrie necessaires
56c      pour la physique sont la latitude pour le rayonnement et
57c      l'aire de la maille quand on veut integrer une grandeur
58c      horizontalement.
59c    - les points de la physique sont les points scalaires de la
60c      la dynamique; numerotation:
61c          1 pour le pole nord
62c          (jjm-1)*iim pour l'interieur du domaine
63c          ngridmx pour le pole sud
64c      ---> ngridmx=2+(jjm-1)*iim
65c
66c     Input :
67c     -------
68c       pucov           covariant zonal velocity
69c       pvcov           covariant meridional velocity
70c       pteta           potential temperature
71c       pps             surface pressure
72c       pmasse          masse d'air dans chaque maille
73c       pts             surface temperature  (K)
74c       callrad         clef d'appel au rayonnement
75c
76c    Output :
77c    --------
78c        pdufi          tendency for the natural zonal velocity (ms-1)
79c        pdvfi          tendency for the natural meridional velocity
80c        pdhfi          tendency for the potential temperature
81c        pdtsfi         tendency for the surface temperature
82c
83c        pdtrad         radiative tendencies  \  both input
84c        pfluxrad       radiative fluxes      /  and output
85c
86c=======================================================================
87c
88c-----------------------------------------------------------------------
89c
90c    0.  Declarations :
91c    ------------------
92
93      include "dimensions.h"
94      include "paramet.h"
95
96      INTEGER ngridmx
97      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
98
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 lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
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
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 zrot(iip1,jjm,llm) ! AdlC May 2014
143      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm), zrfi(ngridmx,llm)
144      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
145c
146      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
147      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
148c
149      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
150      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
151      REAL zdpsrf(ngridmx)
152c
153      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
154      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
155      REAL jH_cur_split,zdt_split
156      LOGICAL debut_split,lafin_split
157      INTEGER isplit
158
159      REAL zsin(iim),zcos(iim),z1(iim)
160      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
161      REAL unskap, pksurcp
162c
163      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
164c
165     
166      REAL SSUM
167
168      LOGICAL,SAVE :: firstcal=.true., debut=.true.
169!      REAL rdayvrai
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  Alvaro de la Camara (May 2014)
353C  46.1 Calcul de la vorticite et passage sur la grille physique
354C  --------------------------------------------------------------
355      DO l=1,llm
356        do i=1,iim
357          do j=1,jjm
358            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
359     $                   + pucov(i,j+1,l) - pucov(i,j,l))
360     $                   / (cu(i,j)+cu(i,j+1))
361     $                   / (cv(i+1,j)+cv(i,j)) *4
362          enddo
363        enddo
364      ENDDO
365
366c   46.champ v:
367c   -----------
368
369      DO l=1,llm
370         DO j=2,jjm
371            ig0=1+(j-2)*iim
372            DO i=1,iim
373               zvfi(ig0+i,l)= 0.5 *
374     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
375               pcvgv(ig0+i,l)= 0.5 *
376     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
377            ENDDO
378               zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
379     &                                +zrot(1,j-1,l)+zrot(1,j,l))
380            DO i=2,iim
381               zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
382     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
383            ENDDO
384         ENDDO
385      ENDDO
386
387
388c   47. champs de vents aux pole nord   
389c   ------------------------------
390c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
391c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
392
393      DO l=1,llm
394
395         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
396         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
397         DO i=2,iim
398            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
399            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
400         ENDDO
401
402         DO i=1,iim
403            zcos(i)   = COS(rlonv(i))*z1(i)
404            zcosbis(i)= COS(rlonv(i))*z1bis(i)
405            zsin(i)   = SIN(rlonv(i))*z1(i)
406            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
407         ENDDO
408
409         zufi(1,l)  = SSUM(iim,zcos,1)/pi
410         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
411         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
412         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
413         zrfi(1, l) = 0.
414      ENDDO
415
416
417c   48. champs de vents aux pole sud:
418c   ---------------------------------
419c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
420c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
421
422      DO l=1,llm
423
424         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
425         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
426         DO i=2,iim
427            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
428            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
429         ENDDO
430
431         DO i=1,iim
432            zcos(i)    = COS(rlonv(i))*z1(i)
433            zcosbis(i) = COS(rlonv(i))*z1bis(i)
434            zsin(i)    = SIN(rlonv(i))*z1(i)
435            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
436         ENDDO
437
438         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
439         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
440         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
441         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
442         zrfi(ngridmx, l) = 0.
443      ENDDO
444c
445c On change de grille, dynamique vers physiq, pour le flux de masse verticale
446      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
447
448c-----------------------------------------------------------------------
449c   Appel de la physique:
450c   ---------------------
451
452
453
454!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
455      zdt_split=dtphys/nsplit_phys
456      zdufic(:,:)=0.
457      zdvfic(:,:)=0.
458      zdtfic(:,:)=0.
459      zdqfic(:,:,:)=0.
460
461#ifdef CPP_PHYS
462
463       do isplit=1,nsplit_phys
464
465         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
466         debut_split=debut.and.isplit==1
467         lafin_split=lafin.and.isplit==nsplit_phys
468
469        CALL call_physiq(ngridmx,llm,nqtot,tname,
470     &                   debut_split,lafin_split,
471     &                   jD_cur,jH_cur_split,zdt_split,
472     &                   zplev,zplay,
473     &                   zphi,zphis,
474     &                   presnivs,
475     &                   zufi,zvfi,zrfi,ztfi,zqfi,
476     &                   flxwfi,pducov,
477     &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
478                             
479!      if (planet_type=="earth") then
480!
481!         CALL physiq (ngridmx,
482!     .             llm,
483!     .             debut_split,
484!     .             lafin_split,
485!     .             jD_cur,
486!     .             jH_cur_split,
487!     .             zdt_split,
488!     .             zplev,
489!     .             zplay,
490!     .             zphi,
491!     .             zphis,
492!     .             presnivs,
493!     .             zufi,
494!     .             zvfi, zrfi,
495!     .             ztfi,
496!     .             zqfi,
497!     .             flxwfi,
498!     .             zdufi,
499!     .             zdvfi,
500!     .             zdtfi,
501!     .             zdqfi,
502!     .             zdpsrf,
503!     .             pducov)
504!
505!      else if ( planet_type=="generic" ) then
506!
507!         CALL physiq (ngridmx,     !! ngrid
508!     .             llm,            !! nlayer
509!     .             nqtot,          !! nq
510!     .             tname,          !! tracer names from dynamical core (given in infotrac)
511!     .             debut_split,    !! firstcall
512!     .             lafin_split,    !! lastcall
513!     .             jD_cur,         !! pday. see leapfrog
514!     .             jH_cur_split,   !! ptime "fraction of day"
515!     .             zdt_split,      !! ptimestep
516!     .             zplev,          !! pplev
517!     .             zplay,          !! pplay
518!     .             zphi,           !! pphi
519!     .             zufi,           !! pu
520!     .             zvfi,           !! pv
521!     .             ztfi,           !! pt
522!     .             zqfi,           !! pq
523!     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
524!     .             zdufi,          !! pdu
525!     .             zdvfi,          !! pdv
526!     .             zdtfi,          !! pdt
527!     .             zdqfi,          !! pdq
528!     .             zdpsrf,         !! pdpsrf
529!     .             tracerdyn)      !! tracerdyn <-- utilite ???
530!
531!      endif ! of if (planet_type=="earth")
532
533         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
534         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
535         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
536         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
537
538         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
539         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
540         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
541         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
542
543       enddo ! of do isplit=1,nsplit_phys
544
545#endif
546! of #ifdef CPP_PHYS
547
548      zdufi(:,:)=zdufic(:,:)/nsplit_phys
549      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
550      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
551      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
552
553
554500   CONTINUE
555
556c-----------------------------------------------------------------------
557c   transformation des tendances physiques en tendances dynamiques:
558c   ---------------------------------------------------------------
559
560c  tendance sur la pression :
561c  -----------------------------------
562
563      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
564c
565c   62. enthalpie potentielle
566c   ---------------------
567
568      DO l=1,llm
569
570         DO i=1,iip1
571          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
572          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
573         ENDDO
574
575         DO j=2,jjm
576            ig0=1+(j-2)*iim
577            DO i=1,iim
578               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
579            ENDDO
580               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
581         ENDDO
582
583      ENDDO
584
585
586c   62. humidite specifique
587c   ---------------------
588! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
589!      DO iq=1,nqtot
590!         DO l=1,llm
591!            DO i=1,iip1
592!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
593!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
594!            ENDDO
595!            DO j=2,jjm
596!               ig0=1+(j-2)*iim
597!               DO i=1,iim
598!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
599!               ENDDO
600!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
601!            ENDDO
602!         ENDDO
603!      ENDDO
604
605c   63. traceurs
606c   ------------
607C     initialisation des tendances
608      pdqfi(:,:,:,:)=0.
609C
610      DO iq=1,nqtot
611         iiq=niadv(iq)
612         DO l=1,llm
613            DO i=1,iip1
614               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
615               pdqfi(i,jjp1,l,iiq) = 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,iiq) = zdqfi(ig0+i,l,iq)
621               ENDDO
622               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
623            ENDDO
624         ENDDO
625      ENDDO
626
627c   65. champ u:
628c   ------------
629
630      DO l=1,llm
631
632         DO i=1,iip1
633            pdufi(i,1,l)    = 0.
634            pdufi(i,jjp1,l) = 0.
635         ENDDO
636
637         DO j=2,jjm
638            ig0=1+(j-2)*iim
639            DO i=1,iim-1
640               pdufi(i,j,l)=
641     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
642            ENDDO
643            pdufi(iim,j,l)=
644     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
645            pdufi(iip1,j,l)=pdufi(1,j,l)
646         ENDDO
647
648      ENDDO
649
650
651c   67. champ v:
652c   ------------
653
654      DO l=1,llm
655
656         DO j=2,jjm-1
657            ig0=1+(j-2)*iim
658            DO i=1,iim
659               pdvfi(i,j,l)=
660     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
661            ENDDO
662            pdvfi(iip1,j,l) = pdvfi(1,j,l)
663         ENDDO
664      ENDDO
665
666
667c   68. champ v pres des poles:
668c   ---------------------------
669c      v = U * cos(long) + V * SIN(long)
670
671      DO l=1,llm
672
673         DO i=1,iim
674            pdvfi(i,1,l)=
675     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
676            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
677     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
678            pdvfi(i,1,l)=
679     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
680            pdvfi(i,jjm,l)=
681     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
682          ENDDO
683
684         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
685         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
686
687      ENDDO
688
689c-----------------------------------------------------------------------
690
691700   CONTINUE
692 
693      firstcal = .FALSE.
694
695      RETURN
696      END
Note: See TracBrowser for help on using the repository browser.