source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/dynphy_lonlat/calfis.F @ 5134

Last change on this file since 5134 was 4464, checked in by lguez, 22 months ago

Replace stop by call to abort_(physiq|gcm)

  • 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 4464 2023-03-09 17:03:51Z abarral $
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, tracers
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,itr
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         call abort_gcm("calfis", "", 1)
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      itr=0
284      DO iq=1,nqtot
285         IF(.NOT.tracers(iq)%isAdvected) CYCLE
286         itr = itr + 1
287         DO l=1,llm
288            zqfi(1,l,itr) = pq(1,1,l,iq)
289            ig0           = 2
290            DO j=2,jjm
291               DO i = 1, iim
292                  zqfi(ig0,l,itr) = pq(i,j,l,iq)
293                  ig0             = ig0 + 1
294               ENDDO
295            ENDDO
296            zqfi(ig0,l,itr) = pq(1,jjp1,l,iq)
297         ENDDO
298      ENDDO
299
300c   convergence dynamique pour les traceurs "EAU"
301! Earth-specific treatment of first 2 tracers (water)
302       if (planet_type=="earth") then
303        DO iq=1,2
304         DO l=1,llm
305            pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
306            ig0          = 2
307            DO j=2,jjm
308               DO i = 1, iim
309                  pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
310                  ig0             = ig0 + 1
311               ENDDO
312            ENDDO
313            pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
314         ENDDO
315        ENDDO
316       endif ! of if (planet_type=="earth")
317
318
319c   Geopotentiel calcule par rapport a la surface locale:
320c   -----------------------------------------------------
321
322      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
323      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
324      DO l=1,llm
325         DO ig=1,ngridmx
326           zphi(ig,l)=zphi(ig,l)-zphis(ig)
327         ENDDO
328      ENDDO
329
330c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
331c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
332c    de masse est calclue dans advtrac.F 
333c      DO l=1,llm
334c        pvervel(1,l)=pw(1,1,l) * g /apoln
335c        ig0=2
336c       DO j=2,jjm
337c           DO i = 1, iim
338c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
339c              ig0 = ig0 + 1
340c           ENDDO
341c       ENDDO
342c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
343c      ENDDO
344
345c
346c   45. champ u:
347c   ------------
348
349      DO 50 l=1,llm
350
351         DO 25 j=2,jjm
352            ig0 = 1+(j-2)*iim
353            zufi(ig0+1,l)= 0.5 *
354     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
355            pcvgu(ig0+1,l)= 0.5 *
356     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
357            DO 10 i=2,iim
358               zufi(ig0+i,l)= 0.5 *
359     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
360               pcvgu(ig0+i,l)= 0.5 *
361     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
36210         CONTINUE
36325      CONTINUE
364
36550    CONTINUE
366
367
368C  Alvaro de la Camara (May 2014)
369C  46.1 Calcul de la vorticite et passage sur la grille physique
370C  --------------------------------------------------------------
371      DO l=1,llm
372        do i=1,iim
373          do j=1,jjm
374            zrot(i,j,l) = (pvcov(i+1,j,l) - pvcov(i,j,l)
375     $                   + pucov(i,j+1,l) - pucov(i,j,l))
376     $                   / (cu(i,j)+cu(i,j+1))
377     $                   / (cv(i+1,j)+cv(i,j)) *4
378          enddo
379        enddo
380      ENDDO
381
382c   46.champ v:
383c   -----------
384
385      DO l=1,llm
386         DO j=2,jjm
387            ig0=1+(j-2)*iim
388            DO i=1,iim
389               zvfi(ig0+i,l)= 0.5 *
390     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
391               pcvgv(ig0+i,l)= 0.5 *
392     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
393            ENDDO
394               zrfi(ig0 + 1,l)= 0.25 *(zrot(iim,j-1,l)+zrot(iim,j,l)
395     &                                +zrot(1,j-1,l)+zrot(1,j,l))
396            DO i=2,iim
397               zrfi(ig0 + i,l)= 0.25 *(zrot(i-1,j-1,l)+zrot(i-1,j,l)
398     $                   +zrot(i,j-1,l)+zrot(i,j,l))   !  AdlC MAY 2014
399            ENDDO
400         ENDDO
401      ENDDO
402
403
404c   47. champs de vents aux pole nord   
405c   ------------------------------
406c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
407c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
408
409      DO l=1,llm
410
411         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
412         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
413         DO i=2,iim
414            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
415            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
416         ENDDO
417
418         DO i=1,iim
419            zcos(i)   = COS(rlonv(i))*z1(i)
420            zcosbis(i)= COS(rlonv(i))*z1bis(i)
421            zsin(i)   = SIN(rlonv(i))*z1(i)
422            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
423         ENDDO
424
425         zufi(1,l)  = SSUM(iim,zcos,1)/pi
426         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
427         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
428         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
429         zrfi(1, l) = 0.
430      ENDDO
431
432
433c   48. champs de vents aux pole sud:
434c   ---------------------------------
435c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
436c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
437
438      DO l=1,llm
439
440         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
441         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
442         DO i=2,iim
443            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
444            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
445         ENDDO
446
447         DO i=1,iim
448            zcos(i)    = COS(rlonv(i))*z1(i)
449            zcosbis(i) = COS(rlonv(i))*z1bis(i)
450            zsin(i)    = SIN(rlonv(i))*z1(i)
451            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
452         ENDDO
453
454         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
455         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
456         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
457         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
458         zrfi(ngridmx, l) = 0.
459      ENDDO
460c
461c On change de grille, dynamique vers physiq, pour le flux de masse verticale
462      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
463
464c-----------------------------------------------------------------------
465c   Appel de la physique:
466c   ---------------------
467
468
469
470!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
471      zdt_split=dtphys/nsplit_phys
472      zdufic(:,:)=0.
473      zdvfic(:,:)=0.
474      zdtfic(:,:)=0.
475      zdqfic(:,:,:)=0.
476
477#ifdef CPP_PHYS
478
479       do isplit=1,nsplit_phys
480
481         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
482         debut_split=debut.and.isplit==1
483         lafin_split=lafin.and.isplit==nsplit_phys
484
485!      if (planet_type=="earth") then
486        CALL call_physiq(ngridmx,llm,nqtot,tracers(:)%name,
487     &                   debut_split,lafin_split,
488     &                   jD_cur,jH_cur_split,zdt_split,
489     &                   zplev,zplay,
490     &                   zpk,zphi,zphis,
491     &                   presnivs,
492     &                   zufi,zvfi,zrfi,ztfi,zqfi,
493     &                   flxwfi,pducov,
494     &                   zdufi,zdvfi,zdtfi,zdqfi,zdpsrf)
495!
496!      else if ( planet_type=="generic" ) then
497!
498!         CALL physiq (ngridmx,     !! ngrid
499!     .             llm,            !! nlayer
500!     .             nqtot,          !! nq
501!     .             tracers(:)%name,!! tracer names from dynamical core (given in infotrac)
502!     .             debut_split,    !! firstcall
503!     .             lafin_split,    !! lastcall
504!     .             jD_cur,         !! pday. see leapfrog
505!     .             jH_cur_split,   !! ptime "fraction of day"
506!     .             zdt_split,      !! ptimestep
507!     .             zplev,          !! pplev
508!     .             zplay,          !! pplay
509!     .             zphi,           !! pphi
510!     .             zufi,           !! pu
511!     .             zvfi,           !! pv
512!     .             ztfi,           !! pt
513!     .             zqfi,           !! pq
514!     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
515!     .             zdufi,          !! pdu
516!     .             zdvfi,          !! pdv
517!     .             zdtfi,          !! pdt
518!     .             zdqfi,          !! pdq
519!     .             zdpsrf,         !! pdpsrf
520!     .             tracerdyn)      !! tracerdyn <-- utilite ???
521!
522!      endif ! of if (planet_type=="earth")
523
524         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
525         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
526         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
527         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
528
529         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
530         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
531         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
532         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
533
534       enddo ! of do isplit=1,nsplit_phys
535
536#endif
537! of #ifdef CPP_PHYS
538
539      zdufi(:,:)=zdufic(:,:)/nsplit_phys
540      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
541      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
542      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
543
544
545500   CONTINUE
546
547c-----------------------------------------------------------------------
548c   transformation des tendances physiques en tendances dynamiques:
549c   ---------------------------------------------------------------
550
551c  tendance sur la pression :
552c  -----------------------------------
553
554      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
555c
556c   62. enthalpie potentielle
557c   ---------------------
558
559      DO l=1,llm
560
561         DO i=1,iip1
562          pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
563          pdhfi(i,jjp1,l) = cpp *  zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
564         ENDDO
565
566         DO j=2,jjm
567            ig0=1+(j-2)*iim
568            DO i=1,iim
569               pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
570            ENDDO
571               pdhfi(iip1,j,l) =  pdhfi(1,j,l)
572         ENDDO
573
574      ENDDO
575
576
577c   62. humidite specifique
578c   ---------------------
579! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
580!      DO iq=1,nqtot
581!         DO l=1,llm
582!            DO i=1,iip1
583!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
584!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
585!            ENDDO
586!            DO j=2,jjm
587!               ig0=1+(j-2)*iim
588!               DO i=1,iim
589!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
590!               ENDDO
591!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
592!            ENDDO
593!         ENDDO
594!      ENDDO
595
596c   63. traceurs
597c   ------------
598C     initialisation des tendances
599      pdqfi(:,:,:,:)=0.
600C
601      itr = 0
602      DO iq=1,nqtot
603         IF(.NOT.tracers(iq)%isAdvected) CYCLE
604         itr = itr + 1
605         DO l=1,llm
606            DO i=1,iip1
607               pdqfi(i,1,l,iq)    = zdqfi(1,l,itr)
608               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,itr)
609            ENDDO
610            DO j=2,jjm
611               ig0=1+(j-2)*iim
612               DO i=1,iim
613                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,itr)
614               ENDDO
615               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,itr)
616            ENDDO
617         ENDDO
618      ENDDO
619
620c   65. champ u:
621c   ------------
622
623      DO l=1,llm
624
625         DO i=1,iip1
626            pdufi(i,1,l)    = 0.
627            pdufi(i,jjp1,l) = 0.
628         ENDDO
629
630         DO j=2,jjm
631            ig0=1+(j-2)*iim
632            DO i=1,iim-1
633               pdufi(i,j,l)=
634     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
635            ENDDO
636            pdufi(iim,j,l)=
637     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
638            pdufi(iip1,j,l)=pdufi(1,j,l)
639         ENDDO
640
641      ENDDO
642
643
644c   67. champ v:
645c   ------------
646
647      DO l=1,llm
648
649         DO j=2,jjm-1
650            ig0=1+(j-2)*iim
651            DO i=1,iim
652               pdvfi(i,j,l)=
653     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
654            ENDDO
655            pdvfi(iip1,j,l) = pdvfi(1,j,l)
656         ENDDO
657      ENDDO
658
659
660c   68. champ v pres des poles:
661c   ---------------------------
662c      v = U * cos(long) + V * SIN(long)
663
664      DO l=1,llm
665
666         DO i=1,iim
667            pdvfi(i,1,l)=
668     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
669            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
670     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
671            pdvfi(i,1,l)=
672     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
673            pdvfi(i,jjm,l)=
674     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
675          ENDDO
676
677         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
678         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
679
680      ENDDO
681
682c-----------------------------------------------------------------------
683
684700   CONTINUE
685 
686      firstcal = .FALSE.
687
688      RETURN
689      END
Note: See TracBrowser for help on using the repository browser.