source: LMDZ5/branches/AI-cosp/libf/dynphy_lonlat/calfis.F @ 5214

Last change on this file since 5214 was 2418, checked in by Ehouarn Millour, 9 years ago

Improving the physics/dynamics interface:

  • added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
  • turned "physiq.F90" into module "physiq_mod.F90" for better control of "physiq" arguments. Extracted embeded "gr_fi_ecrit" as self-standing routine (but note that this routine actually only works in serial mode).

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