source: LMDZ5/branches/testing/libf/dynphy_lonlat/calfis.F @ 3363

Last change on this file since 3363 was 2641, checked in by Laurent Fairhead, 8 years ago

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