source: trunk/LMDZ.COMMON/libf/dyn3d/calfis.F @ 1242

Last change on this file since 1242 was 1189, checked in by emillour, 11 years ago

Common dynamics: a couple of bug fixes:

  • Correctly account for the change in pressure, mass, etc. after modifying surface pressure following a call to the physics.
  • Corrected tracer advection, which is computed using values at the beginning of the time step, so it is done at Matsuno forward step.

EM

File size: 30.0 KB
Line 
1!
2! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
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      USE write_field
34      USE cpdet_mod, only: t2tpot,tpot2t
35 
36! used only for zonal averages
37      USE moyzon_mod
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 (K/s)
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#include "temps.h"
96#include "logic.h"
97
98      INTEGER ngridmx
99      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
100
101#include "comconst.h"
102#include "comvert.h"
103#include "comgeom2.h"
104#include "iniprint.h"
105
106c    Arguments :
107c    -----------
108      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
109      REAL,INTENT(IN) :: jD_cur, jH_cur
110      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
111      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
112      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
113      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
114      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
115      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
116      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
117
118      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
119      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
120      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
121! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
122! qui lui meme ne sert a rien dans la routine telle qu'elle est
123! ecrite, et que j'ai donc commente....
124      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
125      ! NB: pdq is only used to compute pcvgq which is in fact not used...
126
127      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
128      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
129      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
130      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm)  ! Vertical mass flux on dynamics grid
131
132      ! tendencies (in */s) from the physics
133      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
134      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
135      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
136      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
137      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
138
139
140c    Local variables :
141c    -----------------
142
143      INTEGER i,j,l,ig0,ig,iq,iiq
144      REAL zpsrf(ngridmx)
145      REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
146      REAL zphi(ngridmx,llm),zphis(ngridmx)
147
148      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
149      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
150! ADAPTATION GCM POUR CP(T)
151      REAL zteta(ngridmx,llm)
152      REAL zpk(ngridmx,llm)
153
154! RQ SL 13/10/10:
155! Ces calculs ne servent pas.
156! Si necessaire, decommenter ces variables et les calculs...
157!      REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
158!      REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
159
160      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
161      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
162      REAL zdpsrf(ngridmx)
163
164      REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
165      REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
166      REAL jH_cur_split,zdt_split
167      LOGICAL debut_split,lafin_split
168      INTEGER isplit
169
170      REAL zsin(iim),zcos(iim),z1(iim)
171      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
172      REAL unskap, pksurcp
173      save unskap
174
175cIM diagnostique PVteta, Amip2
176      INTEGER,PARAMETER :: ntetaSTD=3
177      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
178      REAL PVteta(ngridmx,ntetaSTD)
179
180      REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
181     
182      REAL SSUM
183
184      LOGICAL,SAVE :: firstcal=.true., debut=.true.
185!      REAL rdayvrai
186
187      LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
188
189! For Titan only right now:
190! to allow for 2D computation of microphys and chemistry
191      LOGICAL,save :: flag_moyzon
192      REAL,allocatable,save :: tmpvar(:,:)
193      REAL,allocatable,save :: tmpvarp1(:,:)
194      REAL,allocatable,save :: tmpvarbar(:)
195      REAL,allocatable,save :: tmpvarbarp1(:)
196      real :: zz1,zz2
197
198c-----------------------------------------------------------------------
199
200c    1. Initialisations :
201c    --------------------
202
203
204      IF ( firstcal )  THEN
205        debut = .TRUE.
206        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
207         write(lunout,*) 'STOP dans calfis'
208         write(lunout,*)
209     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
210         write(lunout,*) '  ngridmx  jjm   iim   '
211         write(lunout,*) ngridmx,jjm,iim
212         STOP
213        ENDIF
214
215        unskap   = 1./ kappa
216
217        flag_moyzon = .false.
218        if(moyzon_ch.or.moyzon_mu) then
219         flag_moyzon = .true.
220         allocate(tmpvar(iip1,llm))
221         allocate(tmpvarp1(iip1,llmp1))
222         allocate(tmpvarbar(llm))
223         allocate(tmpvarbarp1(llmp1))
224        endif
225
226        if (flag_moyzon) call moyzon_init
227
228c----------------------------------------------
229c moyennes globales pour le profil de pression
230       if(planet_type.eq."titan".or.planet_type.eq."venus") then
231        ALLOCATE(plevmoy(llm+1))
232        ALLOCATE(playmoy(llm))
233        ALLOCATE(tmoy(llm))
234        ALLOCATE(tetamoy(llm))
235        ALLOCATE(pkmoy(llm))
236        ALLOCATE(phimoy(0:llm))
237        ALLOCATE(zlevmoy(llm+1))
238        ALLOCATE(zlaymoy(llm))
239        plevmoy=0.
240        do l=1,llmp1
241         do i=1,iip1
242          do j=1,jjp1
243            plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
244          enddo
245         enddo
246        enddo
247        tetamoy=0.
248        pkmoy=0.
249        phimoy=0.
250        do i=1,iip1
251         do j=1,jjp1
252            phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
253         enddo
254        enddo
255        do l=1,llm
256         do i=1,iip1
257          do j=1,jjp1
258            tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
259            pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
260            phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
261          enddo
262         enddo
263        enddo
264        playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
265        call tpot2t(llm,tetamoy,tmoy,pkmoy)
266c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
267        zlaymoy(:) = g*rad*rad/(g*rad-phimoy(:))-rad
268        zlevmoy(1) = phimoy(0)/g
269        DO l=2,llm
270            zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
271            zz2=(plevmoy(l)  +playmoy(l))/(plevmoy(l)  -playmoy(l))
272            zlevmoy(l)=(zz1*zlaymoy(l-1)+zz2*zlaymoy(l))/(zz1+zz2)
273        ENDDO
274        zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
275c-------------------
276c + lat index
277        allocate(klat(ngridmx))
278        klat=0
279        klat(1)  = 1
280        ig0  = 2
281        DO j = 2,jjm
282         do i=0,iim-1
283          klat(ig0+i) = j
284         enddo
285         ig0 = ig0+iim
286        ENDDO
287        klat(ngridmx)  = jjp1
288       endif   ! planet_type=titan
289c----------------------------------------------
290      ELSE
291        debut = .FALSE.
292      ENDIF ! of IF (firstcal)
293
294
295c-----------------------------------------------------------------------
296c   40. transformation des variables dynamiques en variables physiques:
297c   ---------------------------------------------------------------
298
299c   41. pressions au sol (en Pascals)
300c   ----------------------------------
301       
302      zpsrf(1) = pps(1,1)
303
304      ig0  = 2
305      DO j = 2,jjm
306         CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
307         ig0 = ig0+iim
308      ENDDO
309
310      zpsrf(ngridmx) = pps(1,jjp1)
311
312c   42. pression intercouches et fonction d'Exner:
313
314c   -----------------------------------------------------------------
315c     .... zplev  definis aux (llm +1) interfaces des couches  ....
316c     .... zplay  definis aux (  llm )    milieux des couches  ....
317c   -----------------------------------------------------------------
318
319c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
320
321! ADAPTATION GCM POUR CP(T)
322      DO l = 1, llm
323        zpk(   1,l ) = ppk(1,1,l)
324        zteta( 1,l ) = pteta(1,1,l)
325        zplev( 1,l ) = pp(1,1,l)
326        ig0 = 2
327          DO j = 2, jjm
328             DO i =1, iim
329              zpk(   ig0,l ) = ppk(i,j,l)
330              zteta( ig0,l ) = pteta(i,j,l)
331              zplev( ig0,l ) = pp(i,j,l)
332              ig0 = ig0 +1
333             ENDDO
334          ENDDO
335        zpk(   ngridmx,l ) = ppk(1,jjp1,l)
336        zteta( ngridmx,l ) = pteta(1,jjp1,l)
337        zplev( ngridmx,l ) = pp(1,jjp1,l)
338      ENDDO
339        zplev( 1,llmp1 ) = pp(1,1,llmp1)
340        ig0 = 2
341          DO j = 2, jjm
342             DO i =1, iim
343              zplev( ig0,llmp1 ) = pp(i,j,llmp1)
344              ig0 = ig0 +1
345             ENDDO
346          ENDDO
347        zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
348
349      if (flag_moyzon) then
350        tmpvarp1(:,:) = pp(:,1,:)
351        call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
352        zplevbar(1,:) = tmpvarbarp1
353        tmpvar(:,:) = ppk(:,1,:)
354        call moyzon(llm,tmpvar,tmpvarbar)
355        zpkbar(1,:) = tmpvarbar
356        tmpvar(:,:) = pteta(:,1,:)
357        call moyzon(llm,tmpvar,tmpvarbar)
358        ztetabar(1,:) = tmpvarbar
359        call tpot2t(llm,ztetabar(1,:),ztfibar(1,:),zpkbar(1,:))
360        ig0 = 2
361         do j = 2, jjm
362          tmpvarp1(:,:) = pp(:,j,:)
363          call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
364          zplevbar(ig0,:) = tmpvarbarp1
365          tmpvar(:,:) = ppk(:,j,:)
366          call moyzon(llm,tmpvar,tmpvarbar)
367          zpkbar(ig0,:) = tmpvarbar
368          tmpvar(:,:) = pteta(:,j,:)
369          call moyzon(llm,tmpvar,tmpvarbar)
370          ztetabar(ig0,:) = tmpvarbar
371          call tpot2t(llm,ztetabar(ig0,:),ztfibar(ig0,:),zpkbar(ig0,:))
372          ig0 = ig0+1
373          do i=2,iim
374            zplevbar(ig0,:) = zplevbar(ig0-1,:)
375            zpkbar(ig0,:)   = zpkbar(ig0-1,:)
376            ztetabar(ig0,:) = ztetabar(ig0-1,:)
377            ztfibar(ig0,:)  = ztfibar(ig0-1,:)
378            ig0 = ig0+1
379          enddo
380         enddo
381        tmpvarp1(:,:) = pp(:,jjp1,:)
382        call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
383        zplevbar(ngridmx,:) = tmpvarbarp1
384        tmpvar(:,:) = ppk(:,jjp1,:)
385        call moyzon(llm,tmpvar,tmpvarbar)
386        zpkbar(ngridmx,:) = tmpvarbar
387        tmpvar(:,:) = pteta(:,jjp1,:)
388        call moyzon(llm,tmpvar,tmpvarbar)
389        ztetabar(ngridmx,:) = tmpvarbar
390        call tpot2t(llm,ztetabar(ngridmx,:),
391     .                  ztfibar(ngridmx,:),zpkbar(ngridmx,:))
392      endif
393
394c   43. temperature naturelle (en K) et pressions milieux couches .
395c   ---------------------------------------------------------------
396
397! ADAPTATION GCM POUR CP(T)
398         call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
399
400      DO l=1,llm
401
402         pksurcp     =  ppk(1,1,l) / cpp
403         zplay(1,l)  =  preff * pksurcp ** unskap
404!         pcvgt(1,l)  =  pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
405         ig0         = 2
406
407         DO j = 2, jjm
408            DO i = 1, iim
409              pksurcp        = ppk(i,j,l) / cpp
410              zplay(ig0,l)   = preff * pksurcp ** unskap
411!              pcvgt(ig0,l)   = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
412              ig0            = ig0 + 1
413            ENDDO
414         ENDDO
415
416         pksurcp       = ppk(1,jjp1,l) / cpp
417         zplay(ig0,l)  = preff * pksurcp ** unskap
418!         pcvgt(ig0,l)  = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
419
420      ENDDO
421
422      if (flag_moyzon) then
423        zplaybar(:,:) = preff * (zpkbar(:,:)/cpp)**unskap
424      endif
425
426c   43.bis traceurs (tous intensifs)
427c   ---------------
428
429      DO iq=1,nqtot
430         DO l=1,llm
431            zqfi(1,l,iq) = pq(1,1,l,iq)
432            ig0          = 2
433            DO j=2,jjm
434               DO i = 1, iim
435                  zqfi(ig0,l,iq)  = pq(i,j,l,iq)
436                  ig0             = ig0 + 1
437               ENDDO
438            ENDDO
439            zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
440         ENDDO
441      ENDDO  ! boucle sur traceurs
442
443      if (flag_moyzon) then
444       DO iq=1,nqtot
445! RQ: REVOIR A QUOI CA SERT... ET VERIFIER...
446!       iiq=niadv(iq)
447! en fait, iiq=iq...
448! FIN RQ
449        tmpvar(:,:) = pq(:,1,:,iq)
450        call moyzon(llm,tmpvar,tmpvarbar)
451        zqfibar(1,:,iq) = tmpvarbar
452        ig0 = 2
453         do j = 2, jjm
454          tmpvar(:,:) = pq(:,j,:,iq)
455          call moyzon(llm,tmpvar,tmpvarbar)
456          zqfibar(ig0,:,iq) = tmpvarbar
457          ig0 = ig0+1
458          do i=2,iim
459            zqfibar(ig0,:,iq) = zqfibar(ig0-1,:,iq)
460            ig0 = ig0+1
461          enddo
462         enddo
463        tmpvar(:,:) = pq(:,jjp1,:,iq)
464        call moyzon(llm,tmpvar,tmpvarbar)
465        zqfibar(ngridmx,:,iq) = tmpvarbar
466       ENDDO ! of DO iq=1,nqtot
467      endif
468
469! DEBUG
470!     do ig0=1,ngridmx
471!       write(*,'(6(e13.5,1x))') zqfibar(ig0,1,10),zqfi(ig0,1,10),
472!    .                         zqfibar(ig0,llm/2,10),zqfi(ig0,llm/2,10),
473!    .                           zqfibar(ig0,llm,10),zqfi(ig0,llm,10)
474!     enddo
475!     stop
476
477!-----------------
478! RQ SL 13/10/10:
479! Ces calculs ne servent pas.
480! Si necessaire, decommenter ces variables et les calculs...
481!
482!   convergence dynamique pour les traceurs "EAU"
483! Earth-specific treatment of first 2 tracers (water)
484!      if (planet_type=="earth") then
485!       DO iq=1,2
486!        DO l=1,llm
487!           pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
488!           ig0          = 2
489!           DO j=2,jjm
490!              DO i = 1, iim
491!                 pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
492!                 ig0             = ig0 + 1
493!              ENDDO
494!           ENDDO
495!           pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
496!        ENDDO
497!       ENDDO
498!      endif ! of if (planet_type=="earth")
499!----------------
500
501c   Geopotentiel calcule par rapport a la surface locale:
502c   -----------------------------------------------------
503
504      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
505      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
506      DO l=1,llm
507         DO ig=1,ngridmx
508           zphi(ig,l)=zphi(ig,l)-zphis(ig)
509         ENDDO
510      ENDDO
511
512      if (flag_moyzon) then
513        tmpvar(:,1) = pphis(:,1)
514        call moyzon(1,tmpvar(:,1),tmpvarbar(1))
515        zphisbar(1) = tmpvarbar(1)
516        tmpvar(:,:) = pphi(:,1,:)
517        call moyzon(llm,tmpvar,tmpvarbar)
518        zphibar(1,:) = tmpvarbar
519        ig0 = 2
520         do j = 2, jjm
521          tmpvar(:,1) = pphis(:,j)
522          call moyzon(1,tmpvar(:,1),tmpvarbar(1))
523          zphisbar(ig0) = tmpvarbar(1)
524          tmpvar(:,:) = pphi(:,j,:)
525          call moyzon(llm,tmpvar,tmpvarbar)
526          zphibar(ig0,:) = tmpvarbar
527          ig0 = ig0+1
528          do i=2,iim
529            zphisbar(ig0)  = zphisbar(ig0-1)
530            zphibar(ig0,:) = zphibar(ig0-1,:)
531            ig0 = ig0+1
532          enddo
533         enddo
534        tmpvar(:,1) = pphis(:,jjp1)
535        call moyzon(1,tmpvar(:,1),tmpvarbar(1))
536        zphisbar(ngridmx) = tmpvarbar(1)
537        tmpvar(:,:) = pphi(:,jjp1,:)
538        call moyzon(llm,tmpvar,tmpvarbar)
539        zphibar(ngridmx,:) = tmpvarbar
540      endif
541
542c   ....  Calcul de la vitesse  verticale  ( en Pa*m*s  ou Kg/s )  ....
543c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
544c    de masse est calclue dans advtrac.F 
545c      DO l=1,llm
546c        pvervel(1,l)=pw(1,1,l) * g /apoln
547c        ig0=2
548c       DO j=2,jjm
549c           DO i = 1, iim
550c              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
551c              ig0 = ig0 + 1
552c           ENDDO
553c       ENDDO
554c        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
555c      ENDDO
556
557c
558c   45. champ u:
559c   ------------
560
561      DO 50 l=1,llm
562
563         DO 25 j=2,jjm
564            ig0 = 1+(j-2)*iim
565            zufi(ig0+1,l)= 0.5 *
566     $      ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
567!            pcvgu(ig0+1,l)= 0.5 *
568!     $      ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
569            DO 10 i=2,iim
570               zufi(ig0+i,l)= 0.5 *
571     $         ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
572!               pcvgu(ig0+i,l)= 0.5 *
573!     $         ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
57410         CONTINUE
57525      CONTINUE
576
57750    CONTINUE
578
579
580c   46.champ v:
581c   -----------
582
583      DO l=1,llm
584         DO j=2,jjm
585            ig0=1+(j-2)*iim
586            DO i=1,iim
587               zvfi(ig0+i,l)= 0.5 *
588     $         ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
589c               pcvgv(ig0+i,l)= 0.5 *
590c     $         ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
591            ENDDO
592         ENDDO
593      ENDDO
594
595
596c   47. champs de vents aux pole nord   
597c   ------------------------------
598c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
599c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
600
601      DO l=1,llm
602
603         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
604         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1)
605         DO i=2,iim
606            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
607            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
608         ENDDO
609
610         DO i=1,iim
611            zcos(i)   = COS(rlonv(i))*z1(i)
612            zcosbis(i)= COS(rlonv(i))*z1bis(i)
613            zsin(i)   = SIN(rlonv(i))*z1(i)
614            zsinbis(i)= SIN(rlonv(i))*z1bis(i)
615         ENDDO
616
617         zufi(1,l)  = SSUM(iim,zcos,1)/pi
618!         pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
619         zvfi(1,l)  = SSUM(iim,zsin,1)/pi
620!         pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
621
622      ENDDO
623
624
625c   48. champs de vents aux pole sud:
626c   ---------------------------------
627c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
628c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
629
630      DO l=1,llm
631
632         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
633         z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm)
634         DO i=2,iim
635            z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
636            z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
637         ENDDO
638
639         DO i=1,iim
640            zcos(i)    = COS(rlonv(i))*z1(i)
641            zcosbis(i) = COS(rlonv(i))*z1bis(i)
642            zsin(i)    = SIN(rlonv(i))*z1(i)
643            zsinbis(i) = SIN(rlonv(i))*z1bis(i)
644         ENDDO
645
646         zufi(ngridmx,l)  = SSUM(iim,zcos,1)/pi
647!         pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
648         zvfi(ngridmx,l)  = SSUM(iim,zsin,1)/pi
649!         pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
650
651      ENDDO
652c
653      if (planet_type=="earth") then
654#ifdef CPP_EARTH
655! PVtheta calls tetalevel, which is in the (Earth) physics
656cIM calcul PV a teta=350, 380, 405K
657      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
658     $           ztfi,zplay,zplev,
659     $           ntetaSTD,rtetaSTD,PVteta)
660#endif
661      endif
662c
663c On change de grille, dynamique vers physiq, pour le flux de masse verticale
664      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
665
666c-----------------------------------------------------------------------
667c   Appel de la physique:
668c   ---------------------
669
670! Appel de la physique: pose probleme quand on tourne
671! SANS physique, car physiq.F est dans le repertoire phy[]...
672! Il faut une cle CPP_PHYS
673
674! Le fait que les arguments de physiq soient differents selon les planetes
675! ne pose pas de probleme a priori.
676
677!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
678      zdt_split=dtphys/nsplit_phys
679      zdufic(:,:)=0.
680      zdvfic(:,:)=0.
681      zdtfic(:,:)=0.
682      zdqfic(:,:,:)=0.
683
684#ifdef CPP_PHYS
685
686      do isplit=1,nsplit_phys
687
688         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
689         debut_split=debut.and.isplit==1
690         lafin_split=lafin.and.isplit==nsplit_phys
691
692      if (planet_type.eq."earth") then
693         CALL physiq (ngridmx,
694     .             llm,
695     .             debut_split,
696     .             lafin_split,
697     .             jD_cur,
698     .             jH_cur_split,
699     .             zdt_split,
700     .             zplev,
701     .             zplay,
702     .             zphi,
703     .             zphis,
704     .             presnivs,
705     .             zufi,
706     .             zvfi,
707     .             ztfi,
708     .             zqfi,
709     .             flxwfi,
710     .             zdufi,
711     .             zdvfi,
712     .             zdtfi,
713     .             zdqfi,
714     .             zdpsrf,
715cIM diagnostique PVteta, Amip2         
716     .             pducov,
717     .             PVteta)
718
719      else if ( planet_type=="generic" ) then
720
721         CALL physiq (ngridmx,     !! ngrid
722     .             llm,            !! nlayer
723     .             nqtot,          !! nq
724     .             tname,          !! tracer names from dynamical core (given in infotrac)
725     .             debut_split,    !! firstcall
726     .             lafin_split,    !! lastcall
727     .             jD_cur,         !! pday. see leapfrog
728     .             jH_cur_split,   !! ptime "fraction of day"
729     .             zdt_split,      !! ptimestep
730     .             zplev,          !! pplev
731     .             zplay,          !! pplay
732     .             zphi,           !! pphi
733     .             zufi,           !! pu
734     .             zvfi,           !! pv
735     .             ztfi,           !! pt
736     .             zqfi,           !! pq
737     .             flxwfi,         !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
738     .             zdufi,          !! pdu
739     .             zdvfi,          !! pdv
740     .             zdtfi,          !! pdt
741     .             zdqfi,          !! pdq
742     .             zdpsrf,         !! pdpsrf
743     .             tracerdyn)      !! tracerdyn <-- utilite ???
744
745      else if ( planet_type=="mars" ) then
746
747        CALL physiq (ngridmx,    ! ngrid
748     .             llm,          ! nlayer
749     .             nqtot,        ! nq
750     .             debut_split,  ! firstcall
751     .             lafin_split,  ! lastcall
752     .             jD_cur,       ! pday
753     .             jH_cur_split, ! ptime
754     .             zdt_split,    ! ptimestep
755     .             zplev,        ! pplev
756     .             zplay,        ! pplay
757     .             zphi,         ! pphi
758     .             zufi,         ! pu
759     .             zvfi,         ! pv
760     .             ztfi,         ! pt
761     .             zqfi,         ! pq
762     .             flxwfi,       ! pw
763     .             zdufi,        ! pdu
764     .             zdvfi,        ! pdv
765     .             zdtfi,        ! pdt
766     .             zdqfi,        ! pdq
767     .             zdpsrf,       ! pdpsrf
768     .             tracerdyn)    ! tracerdyn (somewhat obsolete)
769
770      else if ((planet_type=="titan").or.(planet_type=="venus")) then
771
772         CALL physiq (ngridmx,
773     .             llm,
774     .             nqtot,
775     .             debut_split,
776     .             lafin_split,
777     .             jD_cur,
778     .             jH_cur_split,
779     .             zdt_split,
780     .             zplev,
781     .             zplay,
782     .             zpk,
783     .             zphi,
784     .             zphis,
785     .             presnivs,
786     .             zufi,
787     .             zvfi,
788     .             ztfi,
789     .             zqfi,
790     .             flxwfi,
791     .             zdufi,
792     .             zdvfi,
793     .             zdtfi,
794     .             zdqfi,
795     .             zdpsrf)
796
797      else ! unknown "planet_type"
798
799        write(lunout,*) "calfis_p: error, unknown planet_type: ",
800     &                  trim(planet_type)
801        stop
802
803      endif ! planet_type
804
805         zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
806         zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
807         ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
808         zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
809
810         zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
811         zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
812         zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
813         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
814
815      enddo ! of do isplit=1,nsplit_phys
816
817! ATTENTION...
818      if (flag_moyzon.and.(nsplit_phys.ne.1)) then
819       print*,"WARNING ! flag_moyzon + nsplit_phys"
820       print*,"zqfibar n'est pas implemente au cours des iterations"
821       print*,"Donc a revoir..."
822       stop
823      endif
824
825#endif
826! #endif of #ifdef CPP_PHYS
827
828      zdufi(:,:)=zdufic(:,:)/nsplit_phys
829      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
830      zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
831      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
832
833
834500   CONTINUE
835
836c-----------------------------------------------------------------------
837c   transformation des tendances physiques en tendances dynamiques:
838c   ---------------------------------------------------------------
839
840c  tendance sur la pression :
841c  -----------------------------------
842
843      CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
844c
845c   62. enthalpie potentielle
846c   ---------------------
847
848! ADAPTATION GCM POUR CP(T)
849      call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
850
851         DO i=1,iip1
852      pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
853      pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
854         ENDDO
855
856         DO j=2,jjm
857            ig0=1+(j-2)*iim
858            DO i=1,iim
859      pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
860            ENDDO
861               pdhfi(iip1,j,:) =  pdhfi(1,j,:)
862         ENDDO
863
864
865c   62. humidite specifique
866c   ---------------------
867! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
868!      DO iq=1,nqtot
869!         DO l=1,llm
870!            DO i=1,iip1
871!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)
872!               pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
873!            ENDDO
874!            DO j=2,jjm
875!               ig0=1+(j-2)*iim
876!               DO i=1,iim
877!                  pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
878!               ENDDO
879!               pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
880!            ENDDO
881!         ENDDO
882!      ENDDO
883
884c   63. traceurs (tous en intensifs)
885c   ------------
886C     initialisation des tendances
887      pdqfi(:,:,:,:)=0.
888C
889       DO iq=1,nqtot
890         iiq=niadv(iq)
891         DO l=1,llm
892            DO i=1,iip1
893               pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
894               pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
895            ENDDO
896            DO j=2,jjm
897               ig0=1+(j-2)*iim
898               DO i=1,iim
899                  pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
900               ENDDO
901               pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
902            ENDDO
903         ENDDO
904       ENDDO
905
906c   65. champ u:
907c   ------------
908
909      DO l=1,llm
910
911         DO i=1,iip1
912            pdufi(i,1,l)    = 0.
913            pdufi(i,jjp1,l) = 0.
914         ENDDO
915
916         DO j=2,jjm
917            ig0=1+(j-2)*iim
918            DO i=1,iim-1
919               pdufi(i,j,l)=
920     $         0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j)
921            ENDDO
922            pdufi(iim,j,l)=
923     $      0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j)
924            pdufi(iip1,j,l)=pdufi(1,j,l)
925         ENDDO
926
927      ENDDO
928
929
930c   67. champ v:
931c   ------------
932
933      DO l=1,llm
934
935         DO j=2,jjm-1
936            ig0=1+(j-2)*iim
937            DO i=1,iim
938               pdvfi(i,j,l)=
939     $         0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j)
940            ENDDO
941            pdvfi(iip1,j,l) = pdvfi(1,j,l)
942         ENDDO
943      ENDDO
944
945
946c   68. champ v pres des poles:
947c   ---------------------------
948c      v = U * cos(long) + V * SIN(long)
949
950      DO l=1,llm
951
952         DO i=1,iim
953            pdvfi(i,1,l)=
954     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
955            pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
956     $      +zdvfi(ngridmx,l)*SIN(rlonv(i))
957            pdvfi(i,1,l)=
958     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
959            pdvfi(i,jjm,l)=
960     $      0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm)
961          ENDDO
962
963         pdvfi(iip1,1,l)  = pdvfi(1,1,l)
964         pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
965
966      ENDDO
967
968c-----------------------------------------------------------------------
969
970700   CONTINUE
971 
972      firstcal = .FALSE.
973
974      RETURN
975      END
Note: See TracBrowser for help on using the repository browser.